Changeset 1325 for framspy/dissimilarity


Ignore:
Timestamp:
08/14/24 02:48:39 (5 months ago)
Author:
Maciej Komosinski
Message:

Minor optimizations

Location:
framspy/dissimilarity
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • framspy/dissimilarity/alignmodel.py

    r1208 r1325  
    22import numpy as np
    33
     4
    45def wcentre(matrix, weights):
    5     sw = weights.sum()
    6     swx = (matrix*weights).sum(axis=1)
    7     swx /= sw
    8     return (matrix.transpose()-swx).transpose()*np.sqrt(weights)
    9    
     6        sw = weights.sum()
     7        swx = (matrix * weights).sum(axis=1)
     8        swx /= sw
     9        return (matrix.transpose() - swx).transpose() * np.sqrt(weights)
     10
     11
    1012def weightedMDS(distances, weights):
    11     n = len(weights)
    12     distances = distances**2
    13     for i in range(2):
    14         distances = wcentre(distances, weights)
    15         distances = distances.T
    16     distances *= -0.5
    17     _, eigenvalues, vh = np.linalg.svd(distances)
    18     W = (vh/np.sqrt(weights)).T
    19     S = np.zeros((n,n))
    20     np.fill_diagonal(S, eigenvalues)
    21     S = S**0.5
    22     dcoords = W.dot(S)
    23     coords = np.zeros((n, 3))
    24     coords[:,0]=dcoords[:,0]
    25     for i in range(1,3):
    26         if n>i:
    27             coords[:,i]=dcoords[:,i]
    28     return coords
    29    
     13        n = len(weights)
     14        distances = distances ** 2
     15        for i in range(2):
     16                distances = wcentre(distances, weights)
     17                distances = distances.T
     18        distances *= -0.5
     19        _, eigenvalues, vh = np.linalg.svd(distances)
     20        W = (vh / np.sqrt(weights)).T
     21        S = np.zeros((n, n))
     22        np.fill_diagonal(S, eigenvalues)
     23        S = S ** 0.5
     24        dcoords = W.dot(S)
     25        coords = np.zeros((n, 3))
     26        coords[:, 0] = dcoords[:, 0]
     27        for i in range(1, 3):
     28                if n > i:
     29                        coords[:, i] = dcoords[:, i]
     30        return coords
     31
    3032
    3133def align(model, fixedZaxis=False):
    32         numparts=model.numparts._value()
     34        numparts = model.numparts._value()
    3335        distmatrix = np.zeros((numparts, numparts), dtype=float)
    3436        for p1 in range(numparts):
    35                 for p2 in range(numparts): #TODO optimize, only calculate a triangle
    36                         P1=model.getPart(p1)
    37                         P2=model.getPart(p2)
     37                for p2 in range(p1 + 1, numparts):  # only calculate a triangle since Euclidean distance is symmetrical
     38                        P1 = model.getPart(p1)
     39                        P2 = model.getPart(p2)
    3840                        if fixedZaxis:
    39                                 #fixed vertical axis, so pretend all points are on the xy plane
     41                                # fixed vertical axis, so pretend all points are on the xy plane
    4042                                z_dist = 0
    4143                        else:
    42                                 z_dist = (P1.z._value()-P2.z._value())**2
    43                         distmatrix[p1,p2]=math.sqrt((P1.x._value()-P2.x._value())**2+(P1.y._value()-P2.y._value())**2+z_dist)
    44        
     44                                z_dist = (P1.z._value() - P2.z._value()) ** 2
     45                        distmatrix[p1, p2] = distmatrix[p2, p1] = math.sqrt((P1.x._value() - P2.x._value()) ** 2 + (P1.y._value() - P2.y._value()) ** 2 + z_dist)
     46
    4547        if model.numjoints._value() > 0:
    46                 weightvector=np.zeros((numparts), dtype=int)
     48                weightvector = np.zeros((numparts), dtype=int)
    4749        else:
    48                 weightvector=np.ones((numparts), dtype=int)
    49        
     50                weightvector = np.ones((numparts), dtype=int)
     51
    5052        for j in range(model.numjoints._value()):
    51                 J=model.getJoint(j)
    52                 weightvector[J.p1._value()]+=1
    53                 weightvector[J.p2._value()]+=1
    54         weightvector=weightvector.astype(float) # convert to float once, since later it would be promoted to float so many times anyway...
     53                J = model.getJoint(j)
     54                weightvector[J.p1._value()] += 1
     55                weightvector[J.p2._value()] += 1
     56        weightvector = weightvector.astype(float) # convert to float once, since later it would be promoted to float so many times anyway...
    5557        coords = weightedMDS(distmatrix, weightvector)
    5658
    … …  
    6668                                P.z = coords[p, 2]
    6769
    68 
    6970        if fixedZaxis:
    7071                if np.shape(coords)[1] > 2:
    71                 #restore original z coordinate
     72                        # restore original z coordinate
    7273                        for p in range(numparts):
    73                                 P=model.getPart(p)
    74                                 coords[p,2]=P.z._value()
    75 
     74                                P = model.getPart(p)
     75                                coords[p, 2] = P.z._value()

  • framspy/dissimilarity/density_distribution.py

    r1322 r1325  
    4141
    4242
    43     def calculateNeighberhood(self,array,mean_coords):
     43    def calculateNeighborhood(self,array,mean_coords):
    4444        """ Calculates number of elements for given sample and set ups the center of this sample
    4545        to the center of mass (calculated by mean of every coordinate)
    … …  
    5555        weight = len(array)
    5656        if weight > 0:
    57             point = [np.mean(array[:,0]),np.mean(array[:,1]),np.mean(array[:,2])]
     57            point = np.mean(array, axis=0)  # equivalent to [np.mean(array[:,0]),np.mean(array[:,1]),np.mean(array[:,2])]
    5858            return weight, point
    5959        else:
    … …  
    6262
    6363    def calculateDistPoints(self,point1, point2):
    64         """ Returns euclidean distance between two points
     64        """ Returns Euclidean distance between two points
    6565        Args (distribution):
    6666            point1 ([float,float,float]) - coordinates of first point
    … …  
    7171
    7272        Returns:
    73             [float]: euclidean distance
     73            [float]: Euclidean distance
    7474        """
    7575        if self.frequency:
    76             return abs(point1-point2)
     76            return abs(point1-point2) # TODO vector instead of scalar returned?
    7777        else:
    78             return np.sqrt(np.sum(np.square(point1-point2)))
     78            return np.linalg.norm(point1-point2, ord=2)
    7979
    8080
    … …  
    107107        """
    108108        lens = len(s1)
    109         indices = []
    110         for i in range(lens):
    111             if s1[i]==0 and s2[i]==0:
    112                     indices.append(i)
     109        indices = [i for i in range(lens) if s1[i]==0 and s2[i]==0]
    113110
    114111        return np.delete(s1, indices), np.delete(s2, indices)
    … …  
    126123        """
    127124        lens = len(s1[0])
    128         indices = []
    129         for i in range(lens):
    130             if s1[1][i]==0 and s2[1][i]==0:
    131                 indices.append(i)
     125        indices = [i for i in range(lens) if s1[1][i]==0 and s2[1][i]==0]
    132126
    133127        s1 = [np.delete(s1[0], indices, axis=0), np.delete(s1[1], indices, axis=0)]
    … …  
    168162                        feature_array.append(len(array[rows]))
    169163                    else:
    170                         weight, point = self.calculateNeighberhood(array[rows],[edges_x[x]+step_x_half,edges_y[y]+step_y_half,edges_z[z]+step_z_half])
     164                        weight, point = self.calculateNeighborhood(array[rows],[edges_x[x]+step_x_half,edges_y[y]+step_y_half,edges_z[z]+step_z_half])
    171165                        feature_array.append(point)
    172166                        weight_array.append(weight)
    … …  
    180174
    181175    def getSignaturesForPair(self,array1,array2):
    182         """Generates signatures for given pair of models represented by array of voxels.
     176        """Generates signatures for a given pair of models represented by array of voxels.
    183177        We calculate space for given models by taking the extremas for each axis and dividing the space by the resolution.
    184178        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.
    … …  
    248242
    249243    def calculateDissimforVoxels(self, voxels1, voxels2):
    250         """Calculates EMD for pair of voxels representing models.
     244        """Calculates EMD for a pair of voxels representing models.
    251245        Args:
    252246            voxels1 np.array([np.array(,dtype=float)]: list of voxels representing model1.
    … …  
    254248
    255249        Returns:
    256             float: dissim for pair of list of voxels
     250            float: dissim for a pair of list of voxels
    257251        """
    258252        numvox1 = len(voxels1)
    … …  
    278272        if self.metric == 'l1':
    279273            if self.frequency:
    280                 out = np.linalg.norm((s1-s2), ord=1)
     274                out = np.linalg.norm(s1-s2, ord=1)
    281275            else:
    282                 out = np.linalg.norm((s1[1]-s2[1]), ord=1)
     276                out = np.linalg.norm(s1[1]-s2[1], ord=1)
    283277
    284278        elif self.metric == 'l2':
    285279            if self.frequency:
    286                 out = np.linalg.norm((s1-s2))
     280                out = np.linalg.norm(s1-s2)
    287281            else:
    288                 out = np.linalg.norm((s1[1]-s2[1]))
     282                out = np.linalg.norm(s1[1]-s2[1])
    289283
    290284        elif self.metric == 'emd':
    … …  
    320314
    321315        Returns:
    322             float: dissim for pair of strings representing models.
     316            float: dissim for a pair of strings representing models.
    323317        """     
    324318
    … …  
    349343        listOfVoxels = [self.getVoxels(g) for g in listOfGeno]
    350344        for i in range(numOfGeno):
    351             for j in range(numOfGeno):
     345            for j in range(numOfGeno): # could only calculate a triangle if the definition of similarity and its calculation guarantees symmetry
    352346                dissimMatrix[i,j] = self.calculateDissimforVoxels(listOfVoxels[i], listOfVoxels[j])
    353347        return dissimMatrix
Note: See TracChangeset for help on using the changeset viewer.