Index: /framspy/dissimilarity/density-distribution.py
===================================================================
--- /framspy/dissimilarity/density-distribution.py	(revision 1223)
+++ /framspy/dissimilarity/density-distribution.py	(revision 1224)
@@ -6,15 +6,17 @@
 
 class DensityDistribution:
-    """ Dissimilarity measures based on the distribution. The structure's bounding box is divided into equal-sized cuboids, the number of which depends on the 'steps' parameter. Then the structure's surface is covered with points, the density of the surface's sampling depends on the 'density' parameter. There are two versions of the measure. In the default version ('frequency'=False) a signature is computed as centroids and a number of samples in each cuboid. In the 'frequency' version FFT is computed from the vector containing the number of samples in each cuboid. The distance between signatures can be computed using EMD, L1, or L2 norms.
+    """Two dissimilarity measures based on the spatial distribution of two Models. The Model bounding box is divided into a grid of equally-sized cuboids, the number of which is the 'resolution' parameter cubed. Then the Model surface is covered with points; the density of the surface sampling is determined by the 'density' parameter. There are two versions of the measure. In the default version ('frequency'=False), a signature of each cuboid is the centroid and the number of samples. In the 'frequency'=True version, FFT is computed from the vector containing the number of samples in each cuboid. The final result of the dissimilarity measure is the distance between the signatures and it can be computed using EMD, L1, or L2 norms (the 'metric' parameter).
     """
-    libm = cdll.LoadLibrary(find_library('m'))
+    
+    libm = cdll.LoadLibrary(find_library('m')) # for disabling/enabling floating point exceptions (division by zero occurs in the EMD library)
     EPSILON = 0.0001
-    def __init__(self, frams_module=None, density = 10, steps = 3, reduce=True, frequency=False, metric = 'emd', fixedZaxis=False, verbose=False):
+    
+    def __init__(self, frams_module=None, density = 10, resolution = 3, reduce_empty=True, frequency=False, metric = 'emd', fixedZaxis=False, verbose=False):
         """ __init__
         Args:
-            density (int, optional): density of samplings for frams.ModelGeometry . Defaults to 10.
-            steps (int, optional): How many steps are used for sampling the space of voxels, 
-                The higher the value, the more accurate the sampling and the longer the calculations. Defaults to 3.
-            reduce (bool, optional): If we should use reduction to remove blank samples. Defaults to True.
+            density (int, optional): density of samplings for frams.ModelGeometry. Defaults to 10.
+            resolution (int, optional): How many intervals are used in each dimension to partition surface samples of Models in the 3D space. 
+                The higher the value, the more detailed the comparison and the longer the calculations. Defaults to 3.
+            reduce_empty (bool, optional): If we should use reduction to remove blank samples. Defaults to True.
             frequency (bool, optional): If we should use frequency distribution. Defaults to False.
             metric (string, optional): The distance metric that should be used ('emd', 'l1', or 'l2'). Defaults to 'emd'.
@@ -27,7 +29,7 @@
 
         self.density = density
-        self.steps = steps
+        self.resolution = resolution
         self.verbose = verbose
-        self.reduce = reduce
+        self.reduce_empty = reduce_empty
         self.frequency = frequency
         self.metric = metric
@@ -76,9 +78,9 @@
         """
         Args:
-            array1 ([type]): array of size n with points representing firsts model 
-            array2 ([type]): array of size n with points representing second model
-
-        Returns:
-            np.array(np.array(,dtype=float)): distance matrix n x n 
+            array1 ([type]): array of size n with points representing the first Model 
+            array2 ([type]): array of size n with points representing the second Model
+
+        Returns:
+            np.array(np.array(,dtype=float)): distance matrix n*n 
         """
         n = len(array1)
@@ -90,5 +92,5 @@
 
 
-    def reduceSignaturesFreq(self,s1,s2):
+    def reduceEmptySignatures_Frequency(self,s1,s2):
         """Removes samples from signatures if corresponding samples for both models have weight 0.
         Args:
@@ -109,5 +111,5 @@
 
 
-    def reduceSignaturesDens(self,s1,s2):
+    def reduceEmptySignatures_Density(self,s1,s2):
         """Removes samples from signatures if corresponding samples for both models have weight 0. 
         Args:
@@ -130,11 +132,11 @@
 
 
-    def getSignatures(self,array,steps_all,step_all):
-        """Generates signature for array representing model. Signature is composed of list of points [x,y,z] (float) and list of weights (int).
-
-        Args:
-            array (np.array(np.array(,dtype=float))): array with voxels representing model
-            steps_all ([np.array(,dtype=float),np.array(,dtype=float),np.array(,dtype=float)]): lists with edges for each step for each axis in order x,y,z
-            step_all ([float,float,float]): [size of step for x axis, size of step for y axis, size of step for y axis] 
+    def getSignatures(self,array,edges3,steps3):
+        """Generates signature for array representing the Model. Signature is composed of list of points [x,y,z] (float) and list of weights (int).
+
+        Args:
+            array (np.array(np.array(,dtype=float))): array with voxels representing the Model
+            edges3 ([np.array(,dtype=float),np.array(,dtype=float),np.array(,dtype=float)]): lists with edges for each step for each axis in order x,y,z
+            steps3 ([float,float,float]): [size of interval for x axis, size of interval for y axis, size of interval for y axis] 
 
         Returns (distribution):
@@ -143,24 +145,24 @@
            signature np.array(,dtype=np.float64): returns signatuere np.array of coefficients
         """
-        x_steps,y_steps,z_steps = steps_all
-        x_step,y_step,z_step=step_all
+        edges_x,edges_y,edges_z = edges3
+        step_x,step_y,step_z=steps3
         feature_array = []
         weight_array = []
-        step_half_x = x_step/2
-        step_half_y = y_step/2
-        step_half_z = z_step/2
-        for x in range(len(x_steps[:-1])):
-            for y in range(len(y_steps[:-1])) :
-                for z in range(len(z_steps[:-1])):
-                    rows=np.where((array[:,0]> x_steps[x]) &
-                                  (array[:,0]<= x_steps[x+1]) &
-                                  (array[:,1]> y_steps[y]) &
-                                  (array[:,1]<= y_steps[y+1]) &
-                                  (array[:,2]> z_steps[z]) &
-                                  (array[:,2]<= z_steps[z+1]))
+        step_x_half = step_x/2
+        step_y_half = step_y/2
+        step_z_half = step_z/2
+        for x in range(len(edges_x[:-1])):
+            for y in range(len(edges_y[:-1])) :
+                for z in range(len(edges_z[:-1])):
+                    rows=np.where((array[:,0]> edges_x[x]) &
+                                  (array[:,0]<= edges_x[x+1]) &
+                                  (array[:,1]> edges_y[y]) &
+                                  (array[:,1]<= edges_y[y+1]) &
+                                  (array[:,2]> edges_z[z]) &
+                                  (array[:,2]<= edges_z[z+1]))
                     if self.frequency:
                         feature_array.append(len(array[rows]))
                     else:
-                        weight, point = self.calculateNeighberhood(array[rows],[x_steps[x]+step_half_x,y_steps[y]+step_half_y,z_steps[z]+step_half_z])
+                        weight, point = self.calculateNeighberhood(array[rows],[edges_x[x]+step_x_half,edges_y[y]+step_y_half,edges_z[z]+step_z_half])
                         feature_array.append(point)
                         weight_array.append(weight)
@@ -175,13 +177,11 @@
     def getSignaturesForPair(self,array1,array2):
         """Generates signatures for given pair of models represented by array of voxels.
-        We calculate space for given models by taking the extremas for each axis and dividing the space by the number of steps.
-        This divided space generate us samples which contains points. Each sample will have new coordinates which are mean of all points from it and weight
-        which equals to the number of points.
+        We calculate space for given models by taking the extremas for each axis and dividing the space by the resolution.
+        This divided space generate us samples which contains points. Each sample will have new coordinates which are mean of all points from it and weight which equals to the number of points.
        
         Args:
             array1 (np.array(np.array(,dtype=float))): array with voxels representing model1
             array2 (np.array(np.array(,dtype=float))): array with voxels representing model2
-            steps (int, optional): How many steps is used for sampling space of voxels. Defaults to self.steps (3).
-                
+
         Returns:
             s1 ([np.array(,dtype=np.float64),np.array(,dtype=np.float64)]): [coordinates of samples, weights] 
@@ -196,17 +196,17 @@
         max_z = np.max([np.max(array1[:,2]),np.max(array2[:,2])])
 
-        # We request self.steps+1 samples since we need self.steps intervals
-        x_steps,x_step = np.linspace(min_x,max_x,self.steps+1,retstep=True)
-        y_steps,y_step = np.linspace(min_y,max_y,self.steps+1,retstep=True)
-        z_steps,z_step = np.linspace(min_z,max_z,self.steps+1,retstep=True)
+        # We request self.resolution+1 samples since we need self.resolution intervals
+        edges_x,step_x = np.linspace(min_x,max_x,self.resolution+1,retstep=True)
+        edges_y,step_y = np.linspace(min_y,max_y,self.resolution+1,retstep=True)
+        edges_z,step_z = np.linspace(min_z,max_z,self.resolution+1,retstep=True)
         
-        for intervals in (x_steps, y_steps, z_steps):  # EPSILON subtracted to deal with boundary voxels (one-sided open intervals and comparisons in loops in function getSignatures())
-            intervals[0] -= self.EPSILON
-
-        steps_all = (x_steps,y_steps,z_steps)
-        step_all = (x_step,y_step,z_step)
+        for edges in (edges_x, edges_y, edges_z):  # EPSILON subtracted to deal with boundary voxels (one-sided open intervals and comparisons in loops in function getSignatures())
+            edges[0] -= self.EPSILON
+
+        edges3 = (edges_x,edges_y,edges_z)
+        steps3 = (step_x,step_y,step_z)
         
-        s1 = self.getSignatures(array1,steps_all,step_all)
-        s2 = self.getSignatures(array2,steps_all,step_all)    
+        s1 = self.getSignatures(array1,edges3,steps3)
+        s2 = self.getSignatures(array2,edges3,steps3)    
         
         return s1,s2
@@ -217,8 +217,8 @@
 
         Args:
-            geno (string): representation of model in one of the formats handled by frams http://www.framsticks.com/a/al_genotype.html
-
-        Returns:
-            np.array([np.array(,dtype=float)]: list of voxels representing model.
+            geno (string): representation of Model in one of the formats supported by Framsticks, http://www.framsticks.com/a/al_genotype.html
+
+        Returns:
+            np.array([np.array(,dtype=float)]: list of voxels representing the Model.
         """
         model = self.frams.Model.newFromString(geno)
@@ -236,5 +236,4 @@
             voxels1 np.array([np.array(,dtype=float)]: list of voxels representing model1.
             voxels2 np.array([np.array(,dtype=float)]: list of voxels representing model2.
-            steps (int, optional): How many steps is used for sampling space of voxels. Defaults to self.steps (3).
 
         Returns:
@@ -250,15 +249,16 @@
             print("Base voxels fig1: ", numvox1, " fig2: ", numvox2)
             print("After reduction voxels fig1: ", sum(s1[1]), " fig2: ", sum(s2[1]))
-            raise ValueError("Bad signature!")
-
-        reduce_fun = self.reduceSignaturesFreq if self.frequency else self.reduceSignaturesDens
-        if self.reduce:
+            raise RuntimeError("Bad signature!")
+
+        reduce_fun = self.reduceEmptySignatures_Frequency if self.frequency else self.reduceEmptySignatures_Density
+        if self.reduce_empty:
             s1, s2 = reduce_fun(s1,s2)
 
             if not self.frequency:
                 if numvox1 != sum(s1[1]) or numvox2 != sum(s2[1]):
-                    print("Voxel reduction didnt work properly")
+                    print("Voxel reduction didn't work properly")
                     print("Base voxels fig1: ", numvox1, " fig2: ", numvox2)
                     print("After reduction voxels fig1: ", sum(s1[1]), " fig2: ", sum(s2[1]))
+                    raise RuntimeError("Voxel reduction error!")
         
         if self.metric == 'l1':
@@ -281,5 +281,5 @@
                 dist_matrix = self.calculateDistanceMatrix(s1[0],s2[0])
 
-            self.libm.fedisableexcept(0x04)  # allowing for operation divide by 0 because pyemd requiers it.
+            self.libm.fedisableexcept(0x04)  # change default flag value - don't cause exceptions when dividing by 0 (pyemd does it)
 
             if self.frequency:
@@ -288,5 +288,5 @@
                 out = emd(s1[1],s2[1],dist_matrix)
 
-            self.libm.feclearexcept(0x04) # disabling operation divide by 0 because framsticks doesnt like it.
+            self.libm.feclearexcept(0x04) # restoring default flag values...
             self.libm.feenableexcept(0x04)
 
@@ -300,7 +300,6 @@
         """Calculates EMD for a pair of genotypes.
         Args:
-            geno1 (string): representation of model1 in one of the formats handled by frams http://www.framsticks.com/a/al_genotype.html
-            geno2 (string): representation of model2 in one of the formats handled by frams http://www.framsticks.com/a/al_genotype.html
-            steps (int, optional): How many steps is used for sampling space of voxels. Defaults to self.steps (3).
+            geno1 (string): representation of model1 in one of the formats supported by Framsticks, http://www.framsticks.com/a/al_genotype.html
+            geno2 (string): representation of model2 in one of the formats supported by Framsticks, http://www.framsticks.com/a/al_genotype.html
 
         Returns:
@@ -314,5 +313,5 @@
 
         if self.verbose == True:
-            print("Steps: ", self.steps)
+            print("Intervals: ", self.resolution)
             print("Geno1:\n",geno1)
             print("Geno2:\n",geno2)
@@ -325,5 +324,5 @@
         """
         Args:
-            listOfGeno ([string]): list of strings representing genotypes in one of the formats handled by frams http://www.framsticks.com/a/al_genotype.html
+            listOfGeno ([string]): list of strings representing genotypes in one of the formats supported by Framsticks, http://www.framsticks.com/a/al_genotype.html
 
         Returns:
