Package mvpa :: Package datasets :: Module metric
[hide private]
[frames] | no frames]

Source Code for Module mvpa.datasets.metric

  1  #emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  #ex: set sts=4 ts=4 sw=4 et: 
  3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  4  # 
  5  #   See COPYING file distributed along with the PyMVPA package for the 
  6  #   copyright and license terms. 
  7  # 
  8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  9  """Classes and functions to provide sense of distances between sample points""" 
 10   
 11  __docformat__ = 'restructuredtext' 
 12   
 13  import numpy as N 
 14   
 15  # 
 16  # Distance functions. 
 17  # 
18 -def cartesianDistance(a, b):
19 """Return Cartesian distance between a and b 20 """ 21 return N.linalg.norm(a-b)
22
23 -def absminDistance(a, b):
24 """Returns dinstance max(\|a-b\|) 25 XXX There must be better name! 26 27 Useful to select a whole cube of a given "radius" 28 """ 29 return max(abs(a-b))
30
31 -def manhattenDistance(a, b):
32 """Return Manhatten distance between a and b 33 """ 34 return sum(abs(a-b))
35
36 -def mahalanobisDistance(x, y=None, w=None):
37 """Caclulcate Mahalanobis distance of the pairs of points. 38 39 :Parameters: 40 `x` 41 first list of points. Rows are samples, columns are 42 features. 43 `y` 44 second list of points (optional) 45 `w` : N.ndarray 46 optional inverse covariance matrix between the points. It is 47 computed if not given 48 49 Inverse covariance matrix can be calculated with the following 50 51 w = N.linalg.solve(N.cov(x.T),N.identity(x.shape[1])) 52 53 or 54 55 w = N.linalg.inv(N.cov(x.T)) 56 """ 57 # see if pairwise between two matrices or just within a single matrix 58 if y is None: 59 # pairwise distances of single matrix 60 # calculate the inverse correlation matrix if necessary 61 if w is None: 62 w = N.linalg.inv(N.cov(x.T)) 63 64 # get some shapes of the data 65 mx, nx = x.shape 66 #mw, nw = w.shape 67 68 # allocate for the matrix to fill 69 d = N.zeros((mx, mx), dtype=N.float32) 70 for i in range(mx-1): 71 # get the current row to compare 72 xi = x[i, :] 73 # replicate the row 74 xi = xi[N.newaxis, :].repeat(mx-i-1, axis=0) 75 # take the distance between all the matrices 76 dc = x[i+1:mx, :] - xi 77 # scale the distance by the correlation 78 d[i+1:mx, i] = N.real(N.sum((N.inner(dc, w) * N.conj(dc)), 1)) 79 # fill the other direction of the matrix 80 d[i, i+1:mx] = d[i+1:mx, i].T 81 else: 82 # is between two matrixes 83 # calculate the inverse correlation matrix if necessary 84 if w is None: 85 # calculate over all points 86 w = N.linalg.inv(N.cov(N.concatenate((x, y)).T)) 87 88 # get some shapes of the data 89 mx, nx = x.shape 90 my, ny = y.shape 91 92 # allocate for the matrix to fill 93 d = N.zeros((mx, my), dtype=N.float32) 94 95 # loop over shorter of two dimensions 96 if mx <= my: 97 # loop over the x patterns 98 for i in range(mx): 99 # get the current row to compare 100 xi = x[i, :] 101 # replicate the row 102 xi = xi[N.newaxis, :].repeat(my, axis=0) 103 # take the distance between all the matrices 104 dc = xi - y 105 # scale the distance by the correlation 106 d[i, :] = N.real(N.sum((N.inner(dc, w) * N.conj(dc)), 1)) 107 else: 108 # loop over the y patterns 109 for j in range(my): 110 # get the current row to compare 111 yj = y[j, :] 112 # replicate the row 113 yj = yj[N.newaxis, :].repeat(mx, axis=0) 114 # take the distance between all the matrices 115 dc = x - yj 116 # scale the distance by the correlation 117 d[:, j] = N.real(N.sum((N.inner(dc, w) * N.conj(dc)), 1)) 118 119 # return the dist 120 return N.sqrt(d)
121 122
123 -class Metric(object):
124 """Abstract class for any finder. 125 126 Classes subclasses from this class show know about structure of 127 the data and thus be able to provide information about the 128 neighbors. 129 At least one of the methods (getNeighbors, getNeighbor) has to be 130 overriden in the derived class. 131 NOTE: derived #2 from derived class #1 has to override all methods 132 which were overrident in class #1 133 """ 134
135 - def getNeighbors(self, *args, **kwargs):
136 """Return the list of coordinates for the neighbors. 137 138 By default it simply constracts the list based on 139 the generator getNeighbor 140 """ 141 return [ x for x in self.getNeighbor(*args, **kwargs) ]
142 143
144 - def getNeighbor(self, *args, **kwargs):
145 """Generator to return coordinate of the neighbor. 146 147 Base class contains the simplest implementation, assuming that 148 getNeighbors returns iterative structure to spit out neighbors 149 1-by-1 150 """ 151 for neighbor in self.getNeighbors(*args, **kwargs): 152 yield neighbor
153 154 155
156 -class DescreteMetric(Metric):
157 """Find neighboring points in descretized space 158 159 If input space is descretized and all points fill in 160 N-dimensional cube, this finder returns list of neighboring 161 points for a given distance. 162 163 As input points it operates on discretized values, not absolute 164 coordinates (which are e.g. in mm) 165 """ 166
167 - def __init__(self, 168 elementsize=1, 169 distance_function=cartesianDistance):
170 """ 171 Initialize the class provided @elementsize and @distance_function 172 """ 173 Metric.__init__(self) 174 self.__filter_radius = None 175 self.__filter_coord = None 176 self.__distance_function = distance_function 177 178 # XXX might not need assume compatible spacemetric 179 self.__elementsize = N.array(elementsize, ndmin=1) 180 self.__Ndims = len(self.__elementsize)
181 182
183 - def _computeFilter(self, radius):
184 """ (Re)Computer filter_coord based on given radius 185 """ 186 # compute radius in units of elementsize per axis 187 elementradius_per_axis = float(radius) / self.__elementsize 188 189 # build prototype search space 190 filter_radiuses = N.array([int(N.ceil(abs(x))) 191 for x in elementradius_per_axis]) 192 filter_center = filter_radiuses 193 filter_mask = N.ones( ( filter_radiuses * 2 ) + 1 ) 194 195 # now zero out all too distant elements 196 f_coords = N.transpose( filter_mask.nonzero() ) 197 198 # check all filter element 199 for coord in f_coords: 200 dist = self.__distance_function(coord*self.__elementsize, 201 filter_center*self.__elementsize) 202 # compare with radius 203 if radius < dist: 204 # zero too distant 205 filter_mask[N.array(coord, ndmin=2).T.tolist()] = 0 206 207 self.__filter_coord = N.array( filter_mask.nonzero() ).T \ 208 - filter_center 209 self.__filter_radius = radius
210 211
212 - def getNeighbors(self, origin, radius=0):
213 """Returns coordinates of the neighbors which are within 214 distance from coord 215 216 XXX radius might need to be not a scalar but a vector of 217 scalars to specify search distance in different dimensions 218 differently... but then may be it has to be a tensor to 219 specify orientation etc? :-) so it might not be necessary for 220 now 221 """ 222 if len(origin) != self.__Ndims: 223 raise ValueError("Obtained coordinates [%s] which have different" 224 + " number of dimensions (%d) from known " 225 + " elementsize" % (`origin`, self.__Ndims)) 226 if radius != self.__filter_radius: 227 self._computeFilter(radius) 228 229 # for the ease of future references, it is better to transform 230 # coordinates into tuples 231 return origin + self.__filter_coord
232 233
234 - def _setFilter(self, filter_coord):
235 """Lets allow to specify some custom filter to use 236 """ 237 self.__filter_coord = filter_coord
238 239
240 - def _getFilter(self):
241 """Lets allow to specify some custom filter to use 242 """ 243 return self.__filter_coord
244 245 filter_coord = property(fget=_getFilter, fset=_setFilter) 246 elementsize = property(fget=lambda self: self.__elementsize)
247 248 # Template for future classes 249 # 250 # class MeshMetric(Metric): 251 # """Return list of neighboring points on a mesh 252 # """ 253 # def getNeighbors(self, origin, distance=0): 254 # """Return neighbors""" 255 # raise NotImplementedError 256