changes for Python 3.7+ compatibility; minor bug fixes

Readme.md

@@ -19,7 +19,7 @@ The algorithms are written in  Python and at the moment work only with SWC files
 With Conda (Linux or Windows):
 
 1. Create a new environment: 
->conda create --name regmaxsn python=2.7 numpy scipy pillow matplotlib scikit-learn pandas seaborn openpyxl xlrd statsmodels
+>conda create --name regmaxsn python=3.7 numpy scipy pillow matplotlib scikit-learn pandas seaborn openpyxl xlrd statsmodels
 2. Activate environment: 
 >(on Linux) source activate regmaxsn 
 >(on Windows) activate regmaxsn

regmaxsn/core/SWCTransforms.py

@@ -1,7 +1,7 @@
 import os
 import numpy as np
-from transforms import compose_matrix
-from swcFuncs import readSWC_numpy, writeSWC_numpy
+from regmaxsn.core.transforms import compose_matrix
+from regmaxsn.core.swcFuncs import readSWC_numpy, writeSWC_numpy
 
 
 def three32BitInt2complexList(arr):
@@ -39,7 +39,7 @@ def __init__(self, refSWC, SWC2Align, gridSize):
         elif type(refSWC) == np.ndarray:
             self.refSWCPts = refSWC[:, 2:5]
         else:
-            raise(ValueError('Unknown data in SWC2Align'))
+            raise ValueError('Unknown data in SWC2Align')
 
         self.refCenter = self.refSWCPts.mean(axis=0)
         refVox = np.array(np.round(self.refSWCPts / gridSize), np.int32)
@@ -53,7 +53,7 @@ def __init__(self, refSWC, SWC2Align, gridSize):
             self.SWC2AlignFull = SWC2Align
             self.headr = ''
         else:
-            raise(ValueError('Unknown data in SWC2Align'))
+            raise ValueError('Unknown data in SWC2Align')
         self.SWC2AlignPts = self.SWC2AlignFull[:, 2:5]
         self.center = self.SWC2AlignPts.mean(axis=0)
 
@@ -111,7 +111,7 @@ def __init__(self, refSWC, SWC2Align, gridSize):
         elif type(refSWC) == np.ndarray:
             self.refSWCPts = refSWC[:, 2:5]
         else:
-            raise(ValueError('Unknown data in SWC2Align'))
+            raise ValueError('Unknown data in SWC2Align')
 
         refCenter = self.refSWCPts.mean(axis=0)
         refSWCPtsCentered = self.refSWCPts - refCenter
@@ -125,7 +125,7 @@ def __init__(self, refSWC, SWC2Align, gridSize):
             self.SWC2AlignFull = SWC2Align
             self.headr = ''
         else:
-            raise(ValueError('Unknown data in SWC2Align'))
+            raise ValueError('Unknown data in SWC2Align')
 
         self.SWC2AlignPts = self.SWC2AlignFull[:, 2:5].copy()
         self.center = self.SWC2AlignPts.mean(axis=0)
@@ -165,7 +165,7 @@ def __iter__(self):
         self.pointsDone = 0
         return self
 
-    def next(self):
+    def __next__(self):
 
         if self.pointsDone < len(self.arg1):
             toReturn = (self.arg1[self.pointsDone], self.arg2)

regmaxsn/core/farthestPointStats.py

@@ -25,7 +25,7 @@ def maxDistStats(swcFiles):
 
     for swcInd1, swcData1 in enumerate(swcDatas):
 
-        swcInds = range(len(swcFiles))
+        swcInds = list(range(len(swcFiles)))
         swcInds.remove(swcInd1)
         swcDataSeries1 = swcDataSeries[swcInd1]
 

regmaxsn/core/iterativeRegistration.py

@@ -5,6 +5,9 @@
 import shutil
 import json
 import subprocess
+from functools import reduce
+import pathlib as pl
+
 
 def transPreference(x, y):
     """
@@ -71,8 +74,6 @@ def __init__(self, refSWC, gridSizes, rotBounds, transBounds,
         self.nCPU = nCPU
         self.allFuncs = {'trans': self.transOnce, 'rot': self.rotOnce, 'scale': self.scaleOnce}
 
-
-
     def rotOnce(self, SWC2Align, outFiles, ipParFile):
         """
         Runs exhaustive search to find the best rotation euler angles about XYZ axes that maximize the volume overlap
@@ -102,7 +103,7 @@ def rotOnce(self, SWC2Align, outFiles, ipParFile):
             bestSol = out['bestSol']
             done = out['done']
             bestVal = out['bestVal']
-            print(bestSol, bestVal, done)
+            print((bestSol, bestVal, done))
 
         return bestSol, bestVal, done
 
@@ -135,7 +136,7 @@ def transOnce(self, SWC2Align, outFiles, ipParFile):
             bestSol = out['bestSol']
             done = out['done']
             bestVal = out['bestVal']
-            print(bestSol, bestVal, done)
+            print((bestSol, bestVal, done))
 
         return bestSol, bestVal, done
 
@@ -168,7 +169,7 @@ def scaleOnce(self, SWC2Align, outFiles, ipParFile, scaleBounds):
             bestSol = out['bestSol']
             done = out['done']
             bestVal = out['bestVal']
-            print(bestSol, bestVal, done)
+            print((bestSol, bestVal, done))
 
         return bestSol, bestVal, done
 
@@ -204,8 +205,7 @@ def compare(self, srts, SWC2Align, tempOutFiles, ipParFile, scaleBounds):
             elif g == 'trans':
                 bestSol, bestVal, done = self.transOnce(SWC2Align, tempOutFiles[g], ipParFile)
             else:
-                raise('Invalid transformation type ' + g)
-
+                raise ValueError(f'Invalid transformation type {g}')
 
             tempDones[g] = done
 
@@ -216,13 +216,8 @@ def compare(self, srts, SWC2Align, tempOutFiles, ipParFile, scaleBounds):
                 presBestSol = bestSol
                 presBestDone = done
 
-
         return tempDones, presBestSol, presBestVal, presBestDone, presBestTrans
 
-
-
-
-
     def performReg(self, SWC2Align, resFile, scaleBounds,
                    inPartsDir=None, outPartsDir=None,
                    initGuessType='just_centroids',
@@ -246,6 +241,7 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
         """
 
         resDir, expName = os.path.split(resFile[:-4])
+        pl.Path(resDir).mkdir(exist_ok=True)
 
         ipParFile = os.path.join(resDir, 'tmp.json')
         vals = ['trans', 'rot', 'scale']
@@ -280,8 +276,7 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
             totalTranslation = SWC2AlignMean
 
         else:
-            raise(ValueError('Unknown value for argument \'initGuessType\''))
-
+            raise ValueError('Unknown value for argument \'initGuessType\'')
 
         SWC2AlignT = SWC2AlignLocal
 
@@ -293,7 +288,6 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
             done = False
             srts = ['rot', 'trans']
 
-
             while not done:
 
                 tempDones, bestSol, bestVal, lDone, g = self.compare(srts, SWC2AlignT, tempOutFiles, ipParFile, None)
@@ -313,7 +307,7 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
                     else:
                         totalTransform = np.dot(presTrans, totalTransform)
 
-                print(str(iterationNo) + g)
+                print((str(iterationNo) + g))
 
                 bestVals[bestVal] = {"outFile": outFile, "outFileSol": outFileSol,
                                      "totalTransform": totalTransform,
@@ -342,7 +336,7 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
                 presTrans = np.array(pars['transMat'])
                 totalTransform = np.dot(presTrans, totalTransform)
 
-            print(str(iterationNo) + 's')
+            print((str(iterationNo) + 's'))
 
             bestVals[bestVal] = {"outFile": outFile, "outFileSol": outFileSol,
                                  "totalTransform": totalTransform,
@@ -387,7 +381,7 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
                     else:
                         totalTransform = np.dot(presTrans, totalTransform)
 
-                print(str(iterationNo) + g)
+                print((str(iterationNo) + g))
 
                 bestVals[bestVal] = {"outFile": outFile, "outFileSol": outFileSol,
                                      "totalTransform": totalTransform,
@@ -404,7 +398,7 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
         totalTranslation = bestVals[championBestVal]["totalTranslation"]
         bestIterIndicator = bestVals[championBestVal]["iterationIndicator"]
 
-        print("bestIter: {}, bestVal: {}".format(bestIterIndicator, championBestVal))
+        print(("bestIter: {}, bestVal: {}".format(bestIterIndicator, championBestVal)))
 
         totalTransform[:3, 3] += totalTranslation
 
@@ -447,7 +441,7 @@ def performReg(self, SWC2Align, resFile, scaleBounds,
                                            )
 
             else:
-                print('Specified partsDir {} not found'.format(inPartsDir))
+                print(('Specified partsDir {} not found'.format(inPartsDir)))
 
 
         return finalFile, finalSolFile
@@ -525,7 +519,7 @@ def writeFakeSWC(data, fName, extraCol=None):
         toWrite = np.empty((data.shape[0], 7))
 
     toWrite[:, 2:5] = data
-    toWrite[:, 0] = range(1, data.shape[0] + 1)
+    toWrite[:, 0] = list(range(1, data.shape[0] + 1))
     toWrite[:, 1] = 3
     toWrite[:, 5] = 1
     toWrite[:, 6] = -np.arange(1, data.shape[0] + 1)

regmaxsn/core/maxDistanceBasedMetric.py

@@ -43,8 +43,9 @@ def calcMaxDistances(swcList):
 
     unionWithDuplicates = np.concatenate(swcPointSets, axis=0)
     if any(np.abs(unionWithDuplicates).max(axis=0) == 0):
-        raise(ValueError("The list of SWCs all lie on a plane or on  a line and hence do not "
-                         "for a 3D point cloud. Such SWCs are not supported."))
+        raise ValueError(
+            "The list of SWCs all lie on a plane or on a line and hence do not "
+            "for a 3D point cloud. Such SWCs are not supported.")
 
     hull = ConvexHull(unionWithDuplicates)
 

regmaxsn/core/occupancyBasedMeasure.py

@@ -2,6 +2,7 @@
 from collections import Counter
 from pyemd import emd
 
+
 def calcOccupancyDistribution(swcList, voxelSize):
     """
     Returns the distribution of the sum of voxel occupancies across swcs in swcList.
@@ -19,7 +20,7 @@ def calcOccupancyDistribution(swcList, voxelSize):
         voxels.extend(list(aVoxSet))
 
     voxelCounter = Counter(voxels)
-    counts = voxelCounter.values()
+    counts = list(voxelCounter.values())
 
     bins = np.arange(1, len(swcList) + 2) - 0.5
 
@@ -29,7 +30,7 @@ def calcOccupancyDistribution(swcList, voxelSize):
 
     histNormed = histWeighted / float(sum(histWeighted))
 
-    return dict(zip(np.arange(1, len(swcList) + 1), histNormed))
+    return {k + 1: v for k, v in enumerate(histNormed)}
 
 
 def occupancyEMD(swcList, voxelSize):

regmaxsn/core/plotDensities.py

@@ -4,6 +4,7 @@
 from scipy.ndimage import gaussian_filter
 import tifffile
 
+
 class DensityVizualizations(object):
 
     def __init__(self, swcSet, gridUnitSizes, resampleLen,
@@ -25,13 +26,13 @@ def __init__(self, swcSet, gridUnitSizes, resampleLen,
         datas = {}
 
         for swcInd, swcFile in enumerate(swcSet):
-            print('Resamping ' + swcFile)
+            print(('Resamping ' + swcFile))
             totalLen, data = resampleSWC(swcFile, resampleLen, mask=masks[swcInd])
             dataT = np.dot(initTrans, data[:, :3].T).T
             datas[swcFile] = dataT
         self.transMat[:3, :3] = initTrans
 
-        allData = np.concatenate(tuple(datas.itervalues()), axis=0)
+        allData = np.concatenate(tuple(datas.values()), axis=0)
         self.allDataMean = allData.mean(axis=0)
 
         if pcaView == 'closestPCMatch':
@@ -42,7 +43,7 @@ def __init__(self, swcSet, gridUnitSizes, resampleLen,
                 refEvecs, thrash = getPCADetails(refSWC, center=True)
                 fEvecs = np.empty_like(refEvecs)
                 coreff = np.dot(refEvecs.T, evecs)
-                possInds = range(refEvecs.shape[1])
+                possInds = list(range(refEvecs.shape[1]))
                 for rowInd in range(refEvecs.shape[1]):
                     bestCorrInd = np.argmax(np.abs(coreff[rowInd, possInds]))
                     fEvecs[:, rowInd] = np.sign(coreff[rowInd, possInds[bestCorrInd]]) * evecs[:, possInds[bestCorrInd]]
@@ -61,16 +62,16 @@ def __init__(self, swcSet, gridUnitSizes, resampleLen,
                 mean2Use = np.loadtxt(refSWC)[:, 2:5].mean(axis=0)
 
             else:
-                raise(ValueError('RefSWC must be specified when pcaView == \'assumeRegistered\''))
+                raise ValueError('RefSWC must be specified when pcaView == \'assumeRegistered\'')
 
         else:
             fEvecs = np.eye(3)
             mean2Use = self.allDataMean
 
 
         self.digDatas = {}
-        for swcFile, data in datas.iteritems():
-            print('Digitizing ' + swcFile)
+        for swcFile, data in datas.items():
+            print(('Digitizing ' + swcFile))
             data -= mean2Use
             data = np.dot(fEvecs.T, data.T).T
             digData = digitizeSWCXYZ(data + mean2Use, gridUnitSizes)
@@ -106,13 +107,13 @@ def calculateDensity(self, swcFiles, sigmas):
 
         for swcFile in swcFiles:
 
-            print('Calculating Density for ' + swcFile)
+            print(('Calculating Density for ' + swcFile))
             densityMat = np.zeros_like(densityMatSum)
-            print('Doing ' + os.path.split(swcFile)[1])
+            print(('Doing ' + os.path.split(swcFile)[1]))
             if swcFile in self.digDatas:
                 digDataTranslated = self.digDatas[swcFile][:, :3] - self.minXYZ
             else:
-                raise(ValueError(swcFile + ' not initialized in constructing DensityVizualizations object'))
+                raise ValueError(swcFile + ' not initialized in constructing DensityVizualizations object')
             densityMat[digDataTranslated[:, 0], digDataTranslated[:, 1], digDataTranslated[:, 2]] = 1
             densityMatSum += densityMat
             del densityMat

regmaxsn/core/rotOnce.py

@@ -1,4 +1,4 @@
-from SWCTransforms import SWCRotate, ArgGenIterator, objFun
+from regmaxsn.core.SWCTransforms import SWCRotate, ArgGenIterator, objFun
 import multiprocessing as mp
 import numpy as np
 import json
@@ -31,16 +31,16 @@
 for gridInd, gridSize in enumerate(gridSizes):
 
     if debugging:
-        print('Gridsize:' + str(gridSize))
+        print(('Gridsize:' + str(gridSize)))
     stepSize = stepSizes[gridInd]
     if debugging:
-        print('Stepsize: ' + str(np.rad2deg(stepSize)))
+        print(('Stepsize: ' + str(np.rad2deg(stepSize))))
     bounds = (np.array(bounds).T - np.array(bestSol)).T
     boundsRoundedUp = np.sign(bounds) * np.ceil(np.abs(bounds) / stepSize) * stepSize
     possiblePts1D = [np.round(bestSol[ind] + np.arange(x[0], x[1] + stepSize, stepSize), 3).tolist()
                      for ind, x in enumerate(boundsRoundedUp)]
     if debugging:
-        print(np.rad2deg([bestSol[ind] + x for ind, x in enumerate(boundsRoundedUp)]))
+        print((np.rad2deg([bestSol[ind] + x for ind, x in enumerate(boundsRoundedUp)])))
     possiblePts3D = np.round(list(product(*possiblePts1D)), 6).tolist()
     argGen = ArgGenIterator(possiblePts3D, SWCDatas[gridInd])
     funcVals = pool.map_async(objFun, argGen).get(1800)
@@ -54,21 +54,21 @@
 
         prevVals = [objFun((x, SWCDatas[gridInd - 1])) for x in minimzers]
         bestSol = minimzers[np.argmin(prevVals)]
-    bounds = map(lambda x: [x - np.sqrt(2) * stepSize, x + np.sqrt(2) * stepSize], bestSol)
+    bounds = [[x - np.sqrt(2) * stepSize, x + np.sqrt(2) * stepSize] for x in bestSol]
 
     if debugging:
         bestVal = objFun((bestSol, SWCDatas[gridInd]))
-        print(np.rad2deg(bestSol), bestVal)
+        print((np.rad2deg(bestSol), bestVal))
 
 
 if minRes < stepSizes[-1]:
 
     if debugging:
-        print('Stepsize: ' + str(np.rad2deg(minRes)))
+        print(('Stepsize: ' + str(np.rad2deg(minRes))))
     bounds = (np.array(bounds).T - np.array(bestSol)).T
     boundsRoundedUp = np.sign(bounds) * np.ceil(np.abs(bounds) / minRes) * minRes
     if debugging:
-        print(np.rad2deg([bestSol[ind] + x for ind, x in enumerate(boundsRoundedUp)]))
+        print((np.rad2deg([bestSol[ind] + x for ind, x in enumerate(boundsRoundedUp)])))
     possiblePts1D = [np.round(bestSol[ind] + np.arange(x[0], x[1] + minRes, minRes), 3).tolist()
                      for ind, x in enumerate(boundsRoundedUp)]
     possiblePts3D = np.round(list(product(*possiblePts1D)), 6).tolist()
@@ -81,12 +81,12 @@
     bestSol = minimzers[np.argmin(prevVals)]
     if debugging:
         bestVal = objFun((bestSol, SWCDatas[-1]))
-        print(np.rad2deg(bestSol), bestVal)
+        print((np.rad2deg(bestSol), bestVal))
 
 bestVal = objFun((bestSol, SWCDatas[-1]))
 nochange = objFun(([0, 0, 0], SWCDatas[-1]))
 if debugging:
-    print(np.rad2deg(bestSol), bestVal, nochange)
+    print((np.rad2deg(bestSol), bestVal, nochange))
 
 
 done = False