1
2
3
4
5
6
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
17
19 """Return Cartesian distance between a and b
20 """
21 return N.linalg.norm(a-b)
22
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
32 """Return Manhatten distance between a and b
33 """
34 return sum(abs(a-b))
35
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
58 if y is None:
59
60
61 if w is None:
62 w = N.linalg.inv(N.cov(x.T))
63
64
65 mx, nx = x.shape
66
67
68
69 d = N.zeros((mx, mx), dtype=N.float32)
70 for i in range(mx-1):
71
72 xi = x[i, :]
73
74 xi = xi[N.newaxis, :].repeat(mx-i-1, axis=0)
75
76 dc = x[i+1:mx, :] - xi
77
78 d[i+1:mx, i] = N.real(N.sum((N.inner(dc, w) * N.conj(dc)), 1))
79
80 d[i, i+1:mx] = d[i+1:mx, i].T
81 else:
82
83
84 if w is None:
85
86 w = N.linalg.inv(N.cov(N.concatenate((x, y)).T))
87
88
89 mx, nx = x.shape
90 my, ny = y.shape
91
92
93 d = N.zeros((mx, my), dtype=N.float32)
94
95
96 if mx <= my:
97
98 for i in range(mx):
99
100 xi = x[i, :]
101
102 xi = xi[N.newaxis, :].repeat(my, axis=0)
103
104 dc = xi - y
105
106 d[i, :] = N.real(N.sum((N.inner(dc, w) * N.conj(dc)), 1))
107 else:
108
109 for j in range(my):
110
111 yj = y[j, :]
112
113 yj = yj[N.newaxis, :].repeat(mx, axis=0)
114
115 dc = x - yj
116
117 d[:, j] = N.real(N.sum((N.inner(dc, w) * N.conj(dc)), 1))
118
119
120 return N.sqrt(d)
121
122
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
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
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
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
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
179 self.__elementsize = N.array(elementsize, ndmin=1)
180 self.__Ndims = len(self.__elementsize)
181
182
184 """ (Re)Computer filter_coord based on given radius
185 """
186
187 elementradius_per_axis = float(radius) / self.__elementsize
188
189
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
196 f_coords = N.transpose( filter_mask.nonzero() )
197
198
199 for coord in f_coords:
200 dist = self.__distance_function(coord*self.__elementsize,
201 filter_center*self.__elementsize)
202
203 if radius < dist:
204
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
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
230
231 return origin + self.__filter_coord
232
233
235 """Lets allow to specify some custom filter to use
236 """
237 self.__filter_coord = filter_coord
238
239
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
249
250
251
252
253
254
255
256