Merge pull request #144 from lucasimi/feature/improve-vptree

Feature/improve vptree

Author: lucasimi

LICENSE

@@ -186,7 +186,7 @@
       same "printed page" as the copyright notice for easier
       identification within third-party archives.
 
-   Copyright [yyyy] [name of copyright owner]
+   Copyright 2020 Luca Simi
 
    Licensed under the Apache License, Version 2.0 (the "License");
    you may not use this file except in compliance with the License.

src/tdamapper/utils/heap.py

@@ -51,7 +51,7 @@ def __len__(self):
 
     def top(self):
         if not self.__heap:
-            return None
+            return (None, None)
         return self.__heap[0].get()
 
     def pop(self):

src/tdamapper/utils/vptree_flat.py

@@ -1,7 +1,6 @@
 from random import randrange
 
 from tdamapper.utils.metrics import get_metric
-from tdamapper.utils.quickselect import quickselect
 from tdamapper.utils.heap import MaxHeap
 
 
@@ -13,6 +12,31 @@ def _mid(start, end):
     return (start + end) // 2
 
 
+def _partition(data, start, end, p_ord):
+    higher = start
+    for j in range(start, end):
+        j_ord, _, _ = data[j]
+        if j_ord < p_ord:
+            _swap(data, higher, j)
+            higher += 1
+    return higher
+
+
+def _quickselect(data, start, end, k):
+    if (k < start) or (k >= end):
+        return
+    start_, end_, higher = start, end, None
+    while higher != k + 1:
+        p, _, _ = data[k]
+        _swap(data, start_, k)
+        higher = _partition(data, start_ + 1, end_, p)
+        _swap(data, start_, higher - 1)
+        if k <= higher - 1:
+            end_ = higher
+        else:
+            start_ = higher
+
+
 class VPTree:
 
     def __init__(
@@ -57,7 +81,7 @@ class _Build:
 
         def __init__(self, vpt, X):
             self.__distance = vpt._get_distance()
-            self.__dataset = [(0.0, x) for x in X]
+            self.__dataset = [(0.0, x, False) for x in X]
             self.__leaf_capacity = vpt.get_leaf_capacity()
             self.__leaf_radius = vpt.get_leaf_radius()
             pivoting = vpt.get_pivoting()
@@ -71,7 +95,7 @@ def _pivoting_disabled(self, start, end):
             pass
 
         def _pivoting_random(self, start, end):
-            if end - start < 2:
+            if end <= start:
                 return
             pivot = randrange(start, end)
             if pivot > start:
@@ -80,17 +104,17 @@ def _pivoting_random(self, start, end):
         def _furthest(self, start, end, i):
             furthest_dist = 0.0
             furthest = start
-            _, i_point = self.__dataset[i]
+            _, i_point, _ = self.__dataset[i]
             for j in range(start, end):
-                _, j_point = self.__dataset[j]
+                _, j_point, _ = self.__dataset[j]
                 j_dist = self.__distance(i_point, j_point)
                 if j_dist > furthest_dist:
                     furthest = j
                     furthest_dist = j_dist
             return furthest
 
         def _pivoting_furthest(self, start, end):
-            if end - start < 2:
+            if end <= start:
                 return
             rnd = randrange(start, end)
             furthest_rnd = self._furthest(start, end, rnd)
@@ -100,10 +124,10 @@ def _pivoting_furthest(self, start, end):
 
         def _update(self, start, end):
             self.__pivoting(start, end)
-            _, v_point = self.__dataset[start]
-            for i in range(start, end):
-                _, point = self.__dataset[i]
-                self.__dataset[i] = self.__distance(v_point, point), point
+            _, v_point, is_terminal = self.__dataset[start]
+            for i in range(start + 1, end):
+                _, point, _ = self.__dataset[i]
+                self.__dataset[i] = self.__distance(v_point, point), point, is_terminal
 
         def build(self):
             self._build_iter()
@@ -115,159 +139,122 @@ def _build_iter(self):
                 start, end = stack.pop()
                 mid = _mid(start, end)
                 self._update(start, end)
-                _, v_point = self.__dataset[start]
-                quickselect(self.__dataset, start + 1, end, mid)
-                v_radius, _ = self.__dataset[mid]
-                self.__dataset[start] = (v_radius, v_point)
-                if (end - start > 2 * self.__leaf_capacity) and (v_radius > self.__leaf_radius):
+                _, v_point, _ = self.__dataset[start]
+                _quickselect(self.__dataset, start + 1, end, mid)
+                v_radius, _, _ = self.__dataset[mid]
+                if (
+                    (end - start > 2 * self.__leaf_capacity) and
+                    (v_radius > self.__leaf_radius)
+                ):
+                    self.__dataset[start] = (v_radius, v_point, False)
                     stack.append((mid, end))
                     stack.append((start + 1, mid))
+                else:
+                    self.__dataset[start] = (v_radius, v_point, True)
 
     def ball_search(self, point, eps, inclusive=True):
         return self._BallSearch(self, point, eps, inclusive).search()
 
     class _BallSearch:
 
-        def __init__(self, vpt, center, radius, inclusive):
-            self.__distance = vpt._get_distance()
+        def __init__(self, vpt, point, eps, inclusive=True):
             self.__dataset = vpt._get_dataset()
-            self.__leaf_capacity = vpt._get_leaf_capacity()
-            self.__leaf_radius = vpt._get_leaf_radius()
-            self.__center = center
-            self.__radius = radius
-            self.__items = []
+            self.__distance = vpt._get_distance()
+            self.__point = point
+            self.__eps = eps
             self.__inclusive = inclusive
 
-        class _BSVisit:
-
-            def __init__(self, start, end, m_radius):
-                self.__start = start
-                self.__end = end
-                self.__m_radius = m_radius
+        def search(self):
+            return self._search_iter()
 
-            def bounds(self):
-                return self.__start, self.__end, self.__m_radius
+        def _inside(self, dist):
+            if self.__inclusive:
+                return dist <= self.__eps
+            return dist < self.__eps
 
-        def search(self):
-            stack = [self._BSVisit(0, len(self.__dataset), float('inf'))]
+        def _search_iter(self):
+            stack = [(0, len(self.__dataset))]
+            result = []
             while stack:
-                visit = stack.pop()
-                start, end, m_radius = visit.bounds()
-                v_radius, v_point = self.__dataset[start]
-                if (end - start <= 2 * self.__leaf_capacity) or (m_radius <= self.__leaf_radius) or (v_radius <= self.__leaf_radius):
-                    for _, x in self.__dataset[start:end]:
-                        dist = self.__distance(self.__center, x)
+                start, end = stack.pop()
+                v_radius, v_point, is_terminal = self.__dataset[start]
+                if is_terminal:
+                    for _, x, _ in self.__dataset[start:end]:
+                        dist = self.__distance(self.__point, x)
                         if self._inside(dist):
-                            self.__items.append(x)
+                            result.append(x)
                 else:
-                    dist = self.__distance(self.__center, v_point)
-                    if self._inside(dist):
-                        self.__items.append(v_point)
+                    dist = self.__distance(self.__point, v_point)
                     mid = _mid(start, end)
+                    if self._inside(dist):
+                        result.append(v_point)
                     if dist <= v_radius:
-                        fst_start, fst_end, fst_radius = start + 1, mid, v_radius
-                        snd_start, snd_end, snd_radius = mid, end, float('inf')
+                        fst = (start + 1, mid)
+                        snd = (mid, end)
                     else:
-                        fst_start, fst_end, fst_radius = mid, end, float('inf')
-                        snd_start, snd_end, snd_radius = start + 1, mid, v_radius
-                    if abs(dist - v_radius) <= self.__radius:
-                        stack.append(self._BSVisit(snd_start, snd_end, snd_radius))
-                    stack.append(self._BSVisit(fst_start, fst_end, fst_radius))
-            return self.__items
+                        fst = (mid, end)
+                        snd = (start + 1, mid)
+                    if abs(dist - v_radius) <= self.__eps:
+                        stack.append(snd)
+                    stack.append(fst)
+            return result
 
-        def _inside(self, dist):
-            if self.__inclusive:
-                return dist <= self.__radius
-            return dist < self.__radius
+    def knn_search(self, point, k):
+        return self._KnnSearch(self, point, k).search()
 
-    def knn_search(self, point, neighbors):
-        return self._KNNSearch(self, point, neighbors).search()
+    class _KnnSearch:
 
-    class _KNNSearch:
-
-        def __init__(self, vpt, center, neighbors):
-            self.__distance = vpt._VPTree__distance
-            self.__dataset = vpt._VPTree__dataset
-            self.__leaf_capacity = vpt._VPTree__leaf_capacity
-            self.__leaf_radius = vpt._VPTree__leaf_radius
-            self.__center = center
+        def __init__(self, vpt, point, neighbors):
+            self.__dataset = vpt._get_dataset()
+            self.__distance = vpt._get_distance()
+            self.__point = point
             self.__neighbors = neighbors
-            self.__items = MaxHeap()
+            self.__radius = float('inf')
+            self.__result = MaxHeap()
 
         def _get_items(self):
-            while len(self.__items) > self.__neighbors:
-                self.__items.pop()
-            return [x for (_, x) in self.__items]
+            while len(self.__result) > self.__neighbors:
+                self.__result.pop()
+            return [x for (_, x) in self.__result]
 
-        def _get_radius(self):
-            if len(self.__items) < self.__neighbors:
-                return float('inf')
-            furthest_dist, _ = self.__items.top()
-            return furthest_dist
+        def search(self):
+            self._search_iter()
+            return self._get_items()
 
         def _process(self, x):
-            dist = self.__distance(self.__center, x)
-            if dist >= self._get_radius():
+            dist = self.__distance(self.__point, x)
+            if dist >= self.__radius:
                 return dist
-            self.__items.add(dist, x)
-            while len(self.__items) > self.__neighbors:
-                self.__items.pop()
+            self.__result.add(dist, x)
+            while len(self.__result) > self.__neighbors:
+                self.__result.pop()
+            if len(self.__result) == self.__neighbors:
+                self.__radius, _ = self.__result.top()
             return dist
 
-        def pre(self, pre, stack):
-            start, end, _ = pre.bounds()
-            v_radius, v_point = self.__dataset[start]
-            dist = self._process(v_point)
-            mid = _mid(start, end)
-            if dist <= v_radius:
-                fst_start, fst_end, fst_radius = start + 1, mid, v_radius
-                snd_start, snd_end, snd_radius = mid, end, float('inf')
-            else:
-                fst_start, fst_end, fst_radius = mid, end, float('inf')
-                snd_start, snd_end, snd_radius = start + 1, mid, v_radius
-            stack.append((self._KVPost(snd_start, snd_end, snd_radius, dist, v_radius), self.post))
-            stack.append((self._KVPre(fst_start, fst_end, fst_radius), self.pre))
-
-        def post(self, post, stack):
-            start, end, _ = post.bounds()
-            m_radius, dist, v_radius = post.rad()
-            if abs(dist - v_radius) <= self._get_radius():
-                stack.append((self._KVPre(start, end, m_radius), self.pre))
-
-        def search(self):
-            stack = [(self._KVPre(0, len(self.__dataset), float('inf')), self.pre)]
+        def _search_iter(self):
+            PRE, POST = 0, 1
+            self.__result = MaxHeap()
+            stack = [(0, len(self.__dataset), 0.0, PRE)]
             while stack:
-                visit, after = stack.pop()
-                start, end, m_radius = visit.bounds()
-                v_radius, _ = self.__dataset[start]
-                if (end - start <= 2 * self.__leaf_capacity) or (m_radius <= self.__leaf_radius) or (v_radius <= self.__leaf_radius):
-                    for _, x in self.__dataset[start:end]:
+                start, end, thr, action = stack.pop()
+                v_radius, v_point, is_terminal = self.__dataset[start]
+                if is_terminal:
+                    for _, x, _ in self.__dataset[start:end]:
                         self._process(x)
                 else:
-                    after(visit, stack)
+                    if action == PRE:
+                        mid = _mid(start, end)
+                        dist = self._process(v_point)
+                        if dist <= v_radius:
+                            fst_start, fst_end = start + 1, mid
+                            snd_start, snd_end = mid, end
+                        else:
+                            fst_start, fst_end = mid, end
+                            snd_start, snd_end = start + 1, mid
+                        stack.append((snd_start, snd_end, abs(v_radius - dist), POST))
+                        stack.append((fst_start, fst_end, 0.0, PRE))
+                    elif action == POST:
+                        if self.__radius > thr:
+                            stack.append((start, end, 0.0, PRE))
             return self._get_items()
-
-        class _KVPre:
-
-            def __init__(self, start, end, m_radius):
-                self.__start = start
-                self.__end = end
-                self.__m_radius = m_radius
-
-            def bounds(self):
-                return self.__start, self.__end, self.__m_radius
-
-        class _KVPost:
-
-            def __init__(self, start, end, m_radius, dist, v_radius):
-                self.__start = start
-                self.__end = end
-                self.__m_radius = m_radius
-                self.__dist = dist
-                self.__v_radius = v_radius
-
-            def bounds(self):
-                return self.__start, self.__end, self.__m_radius
-
-            def rad(self):
-                return self.__m_radius, self.__dist, self.__v_radius

src/tdamapper/utils/vptree_hier.py

@@ -120,7 +120,10 @@ def _build_rec(self, start, end):
             quickselect(self.__dataset, start + 1, end, mid)
             v_radius, _ = self.__dataset[mid]
             self.__dataset[start] = (v_radius, v_point)
-            if (end - start <= 2 * self.__leaf_capacity) or (v_radius <= self.__leaf_radius):
+            if (
+                (end - start <= 2 * self.__leaf_capacity) or
+                (v_radius <= self.__leaf_radius)
+            ):
                 left = _Leaf(start + 1, mid)
                 right = _Leaf(mid, end)
             else:
@@ -201,6 +204,10 @@ def _get_radius(self):
             furthest_dist, _ = self.__items.top()
             return furthest_dist
 
+        def search(self):
+            self._search_rec(self.__tree)
+            return self._get_items()
+
         def _search_rec(self, tree):
             if tree.is_terminal():
                 start, end = tree.get_bounds()
@@ -221,10 +228,6 @@ def _search_rec(self, tree):
                 if abs(dist - v_radius) <= self._get_radius():
                     self._search_rec(snd)
 
-        def search(self):
-            self._search_rec(self.__tree)
-            return self._get_items()
-
 
 class _Node: