import math from ctypes import cdll, c_uint, c_double, c_void_p import tensorflow as tf import numpy as np import numpy.ctypeslib as npct array_2d_uint = npct.ndpointer(dtype=np.uint, ndim=2, flags='CONTIGUOUS') array_1d_double = npct.ndpointer(dtype=np.double, ndim=1, flags='CONTIGUOUS') array_2d_double = npct.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS') from sklearn.neighbors import NearestNeighbors from library.timing import timing from library.MaxHeap import MaxHeap nbhLib = cdll.LoadLibrary("./library/c/libNeighborhood.so") nbhLib.Neighborhood.rettype = None nbhLib.Neighborhood.argtypes = [c_uint, c_uint, c_uint, array_2d_double, array_2d_uint] nbhLib.NeighborhoodHeuristic.rettype = None nbhLib.NeighborhoodHeuristic.argtypes = [c_uint, c_uint, c_uint, array_2d_double, array_2d_uint] def dist(x,y): return sum(map(lambda z: (z[0] - z[1])*(z[0] - z[1]), zip(x, y))) def maxby(data, fn, startValue=0.0): m = startValue for v in data: m = max(m, fn(v)) return m def distancesToPoint(p, points): w = np.array(np.repeat([p], len(points), axis=0)) d = tf.keras.layers.Subtract()([w, np.array(points)]) t = tf.keras.layers.Dot(axes=(1,1))([d,d]) # As the concrete distance is not needed and sqrt(x) is strict monotone # we avoid here unneccessary calculating of expensive roots. return t.numpy() def calculateCenter(points): if points.shape[0] == 1: return points[0] return tf.keras.layers.Average()(list(points)).numpy() class Ball: def __init__(self, points=None, indices=None, parent=None, center=None): if center is not None: self.center = center elif points is not None: self.center = calculateCenter(np.array(points)) else: raise ParameterError("Missing points or center") self.radius = 0 self.points = [] self.indices = set() self.childs = [] self.parent = parent if points is not None and indices is not None: for (i, r) in zip(indices, distancesToPoint(self.center, points)): self.add(i, r) elif points is not None: raise ParameterError("Missing indices") def findIt(self, r): if r == self.points[0][1]: return 0 upper = len(self.points) - 1 lower = 0 while upper > lower + 1: h = (upper + lower) // 2 if self.points[h][1] >= r: upper = h else: lower = h return upper def add(self, xi, r): if xi in self.indices: return False # Here we know, that x is not in points. newEntry = (xi, r) self.indices.add(xi) # Special case: empty list or new element will extend the radius: # Here we can avoid the search if self.points == [] or r >= self.radius: self.points.append(newEntry) self.radius = r return True # Special case: r <= min radius # Here we can avoid the search if self.points[0][1] >= r: self.points = [newEntry] + self.points return True # Here we know that min radius < r < max radius. # So len(points) >= 2. pos = self.findIt(r) # here shoul be: r(pos) >= r > r(pos - 1) self.points = self.points[:pos] + [newEntry] + self.points[pos: ] return True def remove(self, xi, r): if xi not in self.indices: return False # special case: remove the element with the heighest radius: if self.points[-1][0] == xi: self.indices.remove(xi) self.points = self.points[: -1] if self.points == []: self.radius = 0.0 else: self.radius = self.points[-1][1] return True pos = self.findIt(r) nPoints = len(self.points) # pos is the smalest position with: # r(pos - 1) < r <= r(pos) # Just in case: check that we remove the correct index. while pos < nPoints and r == self.points[pos]: if self.points[pos][0] == xi: self.indices.remove(xi) self.points = self.points[:pos] + self.points[(pos + 1): ] return True pos += 1 return False def divideBall(self, X): indices = [i for (i, _r) in self.points] points = np.array([X[i] for i in indices]) distances = tf.keras.layers.Dot(axes=(1,1))([points, points]).numpy() ball = Ball(points, indices, center=np.zeros(len(points[0]))) h = len(ball) // 2 indicesA = [i for (i, _) in ball.points[:h]] indicesB = [i for (i, _) in ball.points[h:]] pointsA = np.array([X[i] for i in indicesA]) pointsB = np.array([X[i] for i in indicesB]) self.childs = [] print(f"{len(points)} -> <{len(pointsA)}|{len(pointsB)}>") self.childs.append(Ball(pointsA, indicesA, self)) self.childs.append(Ball(pointsB, indicesB, self)) return self.childs def __len__(self): return len(self.points) def smalestBallFor(self, i): if i not in self.indices: return None for c in self.childs: b = c.smalestBallFor(i) if b is not None: return b return self class NNSearch: def __init__(self, nebSize=5, timingDict=None): self.nebSize = nebSize self.neighbourhoods = [] self.timingDict = timingDict def timerStart(self, name): if self.timingDict is not None: if name not in self.timingDict: self.timingDict[name] = timing(name) self.timingDict[name].start() def timerStop(self, name): if self.timingDict is not None: if name in self.timingDict: self.timingDict[name].stop() def neighbourhoodOfItem(self, i): return self.neighbourhoods[i] def fit(self, X, nebSize=None): self.timerStart("NN_fit_all") self.fit_chained(X, nebSize) #self.fit_bruteForce_np(X, nebSize) self.timerStop("NN_fit_all") def fit_optimized(self, X, nebSize=None): self.timerStart("NN_fit_all") if nebSize == None: nebSize = self.nebSize nPoints = len(X) nFeatures = len(X[0]) if nFeatures > 15 or nebSize >= (nPoints // 2): print("Using brute force") self.fit_bruteForce_np(X, nebSize) else: print("Using chained") self.fit_chained(X, nebSize) self.timerStop("NN_fit_all") def fit_bruteForce(self, X, nebSize=None): if nebSize == None: nebSize = self.nebSize self.timerStart("NN_fit_bf_init") isGreaterThan = lambda x, y: x[1] > y[1] self.neighbourhoods = [MaxHeap(nebSize, isGreaterThan, (i, 0.0)) for i in range(len(X))] self.timerStop("NN_fit_bf_init") self.timerStart("NN_fit_bf_loop") for (i, x) in enumerate(X): nbh = self.neighbourhoods[i] for (j, y) in enumerate(X[i+1:]): j += i + 1 self.timerStart("NN_fit_bf_dist") d = dist(x,y) self.timerStop("NN_fit_bf_dist") self.timerStart("NN_fit_bf_insert") nbh.insert((j,d)) self.neighbourhoods[j].insert((i,d)) self.timerStop("NN_fit_bf_insert") self.timerStart("NN_fit_bf_toList") self.neighbourhoods[i] = nbh.toArray(lambda v: v[0]) self.timerStop("NN_fit_bf_toList") self.timerStop("NN_fit_bf_loop") def fit_bruteForce_np(self, X, nebSize=None): self.timerStart("NN_fit_bfnp_init") numOfPoints = len(X) nFeatures = len(X[0]) tX = tf.convert_to_tensor(X) def distancesTo(x): w = np.repeat([x], numOfPoints, axis=0) d = tf.keras.layers.Subtract()([w,tX]) t = tf.keras.layers.Dot(axes=(1,1))([d,d]) return t.numpy() if nebSize == None: nebSize = self.nebSize isGreaterThan = lambda x, y: x[1] > y[1] self.neighbourhoods = [MaxHeap(nebSize, isGreaterThan, (i, 0.0)) for i in range(len(X))] self.timerStop("NN_fit_bfnp_init") self.timerStart("NN_fit_bfnp_loop") for (i, x) in enumerate(X): self.timerStart("NN_fit_bfnp_dist") distances = distancesTo(x) self.timerStop("NN_fit_bfnp_dist") nbh = self.neighbourhoods[i] for (j, y) in enumerate(X[i+1:]): j += i + 1 d = distances[j] self.timerStart("NN_fit_bfnp_insert") nbh.insert((j,d)) self.neighbourhoods[j].insert((i,d)) self.timerStop("NN_fit_bfnp_insert") self.timerStart("NN_fit_bfnp_toList") self.neighbourhoods[i] = nbh.toArray(lambda v: v[0]) self.timerStop("NN_fit_bfnp_toList") self.timerStop("NN_fit_bfnp_loop") def fit_chained(self, X, nebSize=None): self.timerStart("NN_fit_chained_init") if nebSize == None: nebSize = self.nebSize nPoints = len(X) nFeatures = len(X[0]) neigh = NearestNeighbors(n_neighbors=nebSize) neigh.fit(X) self.timerStop("NN_fit_chained_init") self.timerStart("NN_fit_chained_toList") self.neighbourhoods = [ (neigh.kneighbors([x], nebSize, return_distance=False))[0] for (i, x) in enumerate(X) ] self.timerStop("NN_fit_chained_toList") def fit_cLib(self, X, nebSize=None): self.timerStart("NN_fit_cLib_init") if nebSize == None: nebSize = self.nebSize nbh = np.array([np.zeros(nebSize, dtype=np.uint) for i in range(X.shape[0])]) self.timerStop("NN_fit_cLib_init") self.timerStart("NN_fit_cLib_call") nbhLib.Neighborhood(nebSize, X.shape[0], X.shape[1], X, nbh) self.timerStop("NN_fit_cLib_call") self.timerStart("NN_fit_cLib_list") self.neighbourhoods = list(nbh) self.timerStop("NN_fit_cLib_list") def fit_cLibHeuristic(self, X, nebSize=None): self.timerStart("NN_fit_cLib_init") if nebSize == None: nebSize = self.nebSize nbh = np.array([np.zeros(nebSize, dtype=np.uint) for i in range(X.shape[0])]) self.timerStop("NN_fit_cLib_init") self.timerStart("NN_fit_cLib_call") nbhLib.NeighborhoodHeuristic(nebSize, X.shape[0], X.shape[1], X, nbh) self.timerStop("NN_fit_cLib_call") self.timerStart("NN_fit_cLib_list") self.neighbourhoods = list(nbh) self.timerStop("NN_fit_cLib_list") # =============================================================== # Heuristic search # =============================================================== def fit_heuristic(self, X, nebSize=None): if nebSize == None: nebSize = self.nebSize nPoints = len(X) def walkUp(nbh, ball, x, i): while ball.parent is not None: print(f"{i}: up (r: {nbh.getMax()})") oldBall = ball ball = ball.parent for c in ball.childs: if c != oldBall: walkDown(nbh, c, x) def walkDown(nbh, ball, x): if ball is None: return print(f"{i}: down (r: {nbh.getMax()})") if dist(x, ball.center) - ball.radius < nbh.getMax()[1]: if ball.childs == []: for (j, _) in ball.points: nbh.insert((j, dist(x, X[j]))) else: for c in ball.childs: walkDown(nbh, c, x) def countBoles(b): if b is None: return 0 return 1 + sum(map(countBoles, b.childs)) root = Ball(X, range(len(X))) queue = [root] while queue != []: ball = queue[0] queue = queue[1:] if len(ball) <= nebSize: continue queue = ball.divideBall(X) + queue isGreaterThan = lambda x, y: x[1] > y[1] self.neighbourhoods = [MaxHeap(nPoints, isGreaterThan, (i, 0.0)) for i in range(len(X))] print("#B: " + str(countBoles(root))) exit() z = X[0] for (i, x) in enumerate(X): nbh = self.neighbourhoods[i] b = root.smalestBallFor(i) if b.parent is not None: b = b.parent for (j, _) in b.points: d = dist(x, X[j]) nbh.insert((j, d)) walkUp(nbh, b, x, i)