# scipy/scipy

### Subversion checkout URL

You can clone with
or
.
Fetching contributors…

Cannot retrieve contributors at this time

676 lines (582 sloc) 27.579 kB
 # Copyright Anne M. Archibald 2008 # Released under the scipy license import numpy as np from heapq import heappush, heappop import scipy.sparse def minkowski_distance_p(x,y,p=2): """Compute the pth power of the L**p distance between x and y For efficiency, this function computes the L**p distance but does not extract the pth root. If p is 1 or infinity, this is equal to the actual L**p distance. """ x = np.asarray(x) y = np.asarray(y) if p==np.inf: return np.amax(np.abs(y-x),axis=-1) elif p==1: return np.sum(np.abs(y-x),axis=-1) else: return np.sum(np.abs(y-x)**p,axis=-1) def minkowski_distance(x,y,p=2): """Compute the L**p distance between x and y""" x = np.asarray(x) y = np.asarray(y) if p==np.inf or p==1: return minkowski_distance_p(x,y,p) else: return minkowski_distance_p(x,y,p)**(1./p) class Rectangle(object): """Hyperrectangle class. Represents a Cartesian product of intervals. """ def __init__(self, maxes, mins): """Construct a hyperrectangle.""" self.maxes = np.maximum(maxes,mins).astype(np.float) self.mins = np.minimum(maxes,mins).astype(np.float) self.m, = self.maxes.shape def __repr__(self): return "" % zip(self.mins, self.maxes) def volume(self): """Total volume.""" return np.prod(self.maxes-self.mins) def split(self, d, split): """Produce two hyperrectangles by splitting along axis d. In general, if you need to compute maximum and minimum distances to the children, it can be done more efficiently by updating the maximum and minimum distances to the parent. """ # FIXME: do this mid = np.copy(self.maxes) mid[d] = split less = Rectangle(self.mins, mid) mid = np.copy(self.mins) mid[d] = split greater = Rectangle(mid, self.maxes) return less, greater def min_distance_point(self, x, p=2.): """Compute the minimum distance between x and a point in the hyperrectangle.""" return minkowski_distance(0, np.maximum(0,np.maximum(self.mins-x,x-self.maxes)),p) def max_distance_point(self, x, p=2.): """Compute the maximum distance between x and a point in the hyperrectangle.""" return minkowski_distance(0, np.maximum(self.maxes-x,x-self.mins),p) def min_distance_rectangle(self, other, p=2.): """Compute the minimum distance between points in the two hyperrectangles.""" return minkowski_distance(0, np.maximum(0,np.maximum(self.mins-other.maxes,other.mins-self.maxes)),p) def max_distance_rectangle(self, other, p=2.): """Compute the maximum distance between points in the two hyperrectangles.""" return minkowski_distance(0, np.maximum(self.maxes-other.mins,other.maxes-self.mins),p) class KDTree(object): """kd-tree for quick nearest-neighbor lookup This class provides an index into a set of k-dimensional points which can be used to rapidly look up the nearest neighbors of any point. The algorithm used is described in Maneewongvatana and Mount 1999. The general idea is that the kd-tree is a binary trie, each of whose nodes represents an axis-aligned hyperrectangle. Each node specifies an axis and splits the set of points based on whether their coordinate along that axis is greater than or less than a particular value. During construction, the axis and splitting point are chosen by the "sliding midpoint" rule, which ensures that the cells do not all become long and thin. The tree can be queried for the r closest neighbors of any given point (optionally returning only those within some maximum distance of the point). It can also be queried, with a substantial gain in efficiency, for the r approximate closest neighbors. For large dimensions (20 is already large) do not expect this to run significantly faster than brute force. High-dimensional nearest-neighbor queries are a substantial open problem in computer science. The tree also supports all-neighbors queries, both with arrays of points and with other kd-trees. These do use a reasonably efficient algorithm, but the kd-tree is not necessarily the best data structure for this sort of calculation. """ def __init__(self, data, leafsize=10): """Construct a kd-tree. Parameters: =========== data : array-like, shape (n,k) The data points to be indexed. This array is not copied, and so modifying this data will result in bogus results. leafsize : positive integer The number of points at which the algorithm switches over to brute-force. """ self.data = np.asarray(data) self.n, self.m = np.shape(self.data) self.leafsize = int(leafsize) if self.leafsize<1: raise ValueError("leafsize must be at least 1") self.maxes = np.amax(self.data,axis=0) self.mins = np.amin(self.data,axis=0) self.tree = self.__build(np.arange(self.n), self.maxes, self.mins) class node(object): pass class leafnode(node): def __init__(self, idx): self.idx = idx self.children = len(idx) class innernode(node): def __init__(self, split_dim, split, less, greater): self.split_dim = split_dim self.split = split self.less = less self.greater = greater self.children = less.children+greater.children def __build(self, idx, maxes, mins): if len(idx)<=self.leafsize: return KDTree.leafnode(idx) else: data = self.data[idx] #maxes = np.amax(data,axis=0) #mins = np.amin(data,axis=0) d = np.argmax(maxes-mins) maxval = maxes[d] minval = mins[d] if maxval==minval: # all points are identical; warn user? return KDTree.leafnode(idx) data = data[:,d] # sliding midpoint rule; see Maneewongvatana and Mount 1999 # for arguments that this is a good idea. split = (maxval+minval)/2 less_idx = np.nonzero(data<=split)[0] greater_idx = np.nonzero(data>split)[0] if len(less_idx)==0: split = np.amin(data) less_idx = np.nonzero(data<=split)[0] greater_idx = np.nonzero(data>split)[0] if len(greater_idx)==0: split = np.amax(data) less_idx = np.nonzero(data=split)[0] if len(less_idx)==0: # _still_ zero? all must have the same value assert np.all(data==data[0]), "Troublesome data array: %s" % data split = data[0] less_idx = np.arange(len(data)-1) greater_idx = np.array([len(data)-1]) lessmaxes = np.copy(maxes) lessmaxes[d] = split greatermins = np.copy(mins) greatermins[d] = split return KDTree.innernode(d, split, self.__build(idx[less_idx],lessmaxes,mins), self.__build(idx[greater_idx],maxes,greatermins)) def __query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf): side_distances = np.maximum(0,np.maximum(x-self.maxes,self.mins-x)) if p!=np.inf: side_distances**=p min_distance = np.sum(side_distances) else: min_distance = np.amax(side_distances) # priority queue for chasing nodes # entries are: # minimum distance between the cell and the target # distances between the nearest side of the cell and the target # the head node of the cell q = [(min_distance, tuple(side_distances), self.tree)] # priority queue for the nearest neighbors # furthest known neighbor first # entries are (-distance**p, i) neighbors = [] if eps==0: epsfac=1 elif p==np.inf: epsfac = 1/(1+eps) else: epsfac = 1/(1+eps)**p if p!=np.inf and distance_upper_bound!=np.inf: distance_upper_bound = distance_upper_bound**p while q: min_distance, side_distances, node = heappop(q) if isinstance(node, KDTree.leafnode): # brute-force data = self.data[node.idx] ds = minkowski_distance_p(data,x[np.newaxis,:],p) for i in range(len(ds)): if ds[i]distance_upper_bound*epsfac: # since this is the nearest cell, we're done, bail out break # compute minimum distances to the children and push them on if x[node.split_dim]1: dd = np.empty(retshape+(k,),dtype=np.float) dd.fill(np.inf) ii = np.empty(retshape+(k,),dtype=np.int) ii.fill(self.n) elif k==1: dd = np.empty(retshape,dtype=np.float) dd.fill(np.inf) ii = np.empty(retshape,dtype=np.int) ii.fill(self.n) elif k is None: dd = np.empty(retshape,dtype=np.object) ii = np.empty(retshape,dtype=np.object) else: raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None") for c in np.ndindex(retshape): hits = self.__query(x[c], k=k, p=p, distance_upper_bound=distance_upper_bound) if k>1: for j in range(len(hits)): dd[c+(j,)], ii[c+(j,)] = hits[j] elif k==1: if len(hits)>0: dd[c], ii[c] = hits[0] else: dd[c] = np.inf ii[c] = self.n elif k is None: dd[c] = [d for (d,i) in hits] ii[c] = [i for (d,i) in hits] return dd, ii else: hits = self.__query(x, k=k, p=p, distance_upper_bound=distance_upper_bound) if k==1: if len(hits)>0: return hits[0] else: return np.inf, self.n elif k>1: dd = np.empty(k,dtype=np.float) dd.fill(np.inf) ii = np.empty(k,dtype=np.int) ii.fill(self.n) for j in range(len(hits)): dd[j], ii[j] = hits[j] return dd, ii elif k is None: return [d for (d,i) in hits], [i for (d,i) in hits] else: raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None") def __query_ball_point(self, x, r, p=2., eps=0): R = Rectangle(self.maxes, self.mins) def traverse_checking(node, rect): if rect.min_distance_point(x,p)>=r/(1.+eps): return [] elif rect.max_distance_point(x,p)r/(1.+eps): return elif rect1.max_distance_rectangle(rect2, p)max_r result[idx[c_greater]] += node1.children*node2.children idx = idx[(min_r<=r[idx]) & (r[idx]<=max_r)] if len(idx)==0: return if isinstance(node1,KDTree.leafnode): if isinstance(node2,KDTree.leafnode): ds = minkowski_distance(self.data[node1.idx][:,np.newaxis,:], other.data[node2.idx][np.newaxis,:,:], p).ravel() ds.sort() result[idx] += np.searchsorted(ds,r[idx],side='right') else: less, greater = rect2.split(node2.split_dim, node2.split) traverse(node1, rect1, node2.less, less, idx) traverse(node1, rect1, node2.greater, greater, idx) else: if isinstance(node2,KDTree.leafnode): less, greater = rect1.split(node1.split_dim, node1.split) traverse(node1.less, less, node2, rect2, idx) traverse(node1.greater, greater, node2, rect2, idx) else: less1, greater1 = rect1.split(node1.split_dim, node1.split) less2, greater2 = rect2.split(node2.split_dim, node2.split) traverse(node1.less,less1,node2.less,less2,idx) traverse(node1.less,less1,node2.greater,greater2,idx) traverse(node1.greater,greater1,node2.less,less2,idx) traverse(node1.greater,greater1,node2.greater,greater2,idx) R1 = Rectangle(self.maxes, self.mins) R2 = Rectangle(other.maxes, other.mins) if np.shape(r) == (): r = np.array([r]) result = np.zeros(1,dtype=int) traverse(self.tree, R1, other.tree, R2, np.arange(1)) return result[0] elif len(np.shape(r))==1: r = np.asarray(r) n, = r.shape result = np.zeros(n,dtype=int) traverse(self.tree, R1, other.tree, R2, np.arange(n)) return result else: raise ValueError("r must be either a single value or a one-dimensional array of values") def sparse_distance_matrix(self, other, max_distance, p=2.): """Compute a sparse distance matrix Computes a distance matrix between two KDTrees, leaving as zero any distance greater than max_distance. Parameters ========== other : KDTree max_distance : positive float Returns ======= result : dok_matrix Sparse matrix representing the results in "dictionary of keys" format. """ result = scipy.sparse.dok_matrix((self.n,other.n)) def traverse(node1, rect1, node2, rect2): if rect1.min_distance_rectangle(rect2, p)>max_distance: return elif isinstance(node1, KDTree.leafnode): if isinstance(node2, KDTree.leafnode): for i in node1.idx: for j in node2.idx: d = minkowski_distance(self.data[i],other.data[j],p) if d<=max_distance: result[i,j] = d else: less, greater = rect2.split(node2.split_dim, node2.split) traverse(node1,rect1,node2.less,less) traverse(node1,rect1,node2.greater,greater) elif isinstance(node2, KDTree.leafnode): less, greater = rect1.split(node1.split_dim, node1.split) traverse(node1.less,less,node2,rect2) traverse(node1.greater,greater,node2,rect2) else: less1, greater1 = rect1.split(node1.split_dim, node1.split) less2, greater2 = rect2.split(node2.split_dim, node2.split) traverse(node1.less,less1,node2.less,less2) traverse(node1.less,less1,node2.greater,greater2) traverse(node1.greater,greater1,node2.less,less2) traverse(node1.greater,greater1,node2.greater,greater2) traverse(self.tree, Rectangle(self.maxes, self.mins), other.tree, Rectangle(other.maxes, other.mins)) return result def distance_matrix(x,y,p=2,threshold=1000000): """Compute the distance matrix. Computes the matrix of all pairwise distances. Parameters ========== x : array-like, m by k y : array-like, n by k p : float 1<=p<=infinity Which Minkowski p-norm to use. threshold : positive integer If m*n*k>threshold use a python loop instead of creating a very large temporary. Returns ======= result : array-like, m by n """ x = np.asarray(x) m, k = x.shape y = np.asarray(y) n, kk = y.shape if k != kk: raise ValueError("x contains %d-dimensional vectors but y contains %d-dimensional vectors" % (k, kk)) if m*n*k <= threshold: return minkowski_distance(x[:,np.newaxis,:],y[np.newaxis,:,:],p) else: result = np.empty((m,n),dtype=np.float) #FIXME: figure out the best dtype if m
Something went wrong with that request. Please try again.