Fix root-node precision problems

CHANGELOG.md

@@ -15,9 +15,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
   - Fixed tree consistency (single-entry leaf after remove)
   - Fixed tree consistency (nValues) -> verify
   - Fixed bug in qt2.contains()
-- Fixed QT2 inconsistency after root resizing after insert(). [#42](https://github.com/tzaeschke/tinspin-indexes/issues/42)
+- Fixed QT2 inconsistency after root resizing after insert(). 
   Essentially, we force all radii and the center of the root to be powers of two.
-  This should immensely reduce precision problems. 
+  This should immensely reduce precision problems.
+  [#42](https://github.com/tzaeschke/tinspin-indexes/issues/42)
 - Removed unnecessary JUnit test console output. [#45](https://github.com/tzaeschke/tinspin-indexes/pull/45)
 
 ## [2.1.3] - 2023-11-19

src/main/java/org/tinspin/index/BoxMap.java

@@ -147,6 +147,7 @@ static <T> BoxMap<T> createQuadtree(int dims) {
         }
 
         /**
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
          * Create a plain Quadtree.
          * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
          * span an area that is somewhat larger rather than smaller than the actual data.
@@ -157,11 +158,32 @@ static <T> BoxMap<T> createQuadtree(int dims) {
          * @param max             Estimated maximum of all coordinates.
          * @param <T>             Value type
          * @return New Quadtree
+         * @deprecated Please use {@link #createQuadtree(double[], double[], boolean, int)}
          */
+        @Deprecated
         static <T> BoxMap<T> createQuadtree(int dims, int maxNodeCapacity, double[] min, double[] max) {
             return QuadTreeRKD0.create(dims, maxNodeCapacity, min, max);
         }
 
+        /**
+         * Create a plain Quadtree.
+         * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param min             Estimated minimum of all coordinates.
+         * @param max             Estimated maximum of all coordinates.
+         * @param align           Whether center and radius should be aligned to powers of two. Aligning considerably
+         *                        reduces risk of precision problems. Recommended: "true".
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param <T>             Value type
+         * @return New Quadtree
+         */
+        static <T> BoxMap<T> createQuadtree(double[] min, double[] max, boolean align, int maxNodeCapacity) {
+            return QuadTreeRKD0.create(min, max, align, maxNodeCapacity);
+        }
+
         /**
          * Create a Quadtree with hypercube navigation.
          *
@@ -174,7 +196,8 @@ static <T> BoxMap<T> createQuadtreeHC(int dims) {
         }
 
         /**
-         * Create a plain Quadtree.
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
+         * Create a Quadtree with hypercube navigation.
          * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
          * span an area that is somewhat larger rather than smaller than the actual data.
          *
@@ -184,11 +207,32 @@ static <T> BoxMap<T> createQuadtreeHC(int dims) {
          * @param max             Estimated maximum of all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC
+         * @deprecated Please use {@link #createQuadtree(double[], double[], boolean, int)}
          */
+        @Deprecated
         static <T> BoxMap<T> createQuadtreeHC(int dims, int maxNodeCapacity, double[] min, double[] max) {
             return QuadTreeRKD.create(dims, maxNodeCapacity, min, max);
         }
 
+        /**
+         * Create a Quadtree with hypercube navigation.
+         * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param min             Estimated minimum of all coordinates.
+         * @param max             Estimated maximum of all coordinates.
+         * @param align           Whether center and radius should be aligned to powers of two. Aligning considerably
+         *                        reduces risk of precision problems. Recommended: "true".
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param <T>             Value type
+         * @return New QuadtreeHC
+         */
+        static <T> BoxMap<T> createQuadtreeHC(double[] min, double[] max, boolean align, int maxNodeCapacity) {
+            return QuadTreeRKD.create(min, max, align, maxNodeCapacity);
+        }
+
         /**
          * Create an R*Tree.
          *

src/main/java/org/tinspin/index/BoxMultimap.java

@@ -145,7 +145,7 @@ default BoxEntryKnn<T> query1nn(double[] center) {
 
     interface Factory {
         /**
-         * Create an array backed BoxMap. This is only for testing and rather inefficient for large data sets.
+         * Create an array backed BoxMultiMap. This is only for testing and rather inefficient for large data sets.
          *
          * @param dims Number of dimensions.
          * @param size Number of entries.
@@ -168,6 +168,7 @@ static <T> BoxMultimap<T> createQuadtree(int dims) {
         }
 
         /**
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
          * Create a plain Quadtree.
          * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
          * span an area that is somewhat larger rather than smaller than the actual data.
@@ -178,11 +179,32 @@ static <T> BoxMultimap<T> createQuadtree(int dims) {
          * @param max             Estimated maximum of all coordinates.
          * @param <T>             Value type
          * @return New Quadtree
+         * @deprecated Please use {@link #createQuadtree(double[], double[], boolean, int)}
          */
+        @Deprecated
         static <T> BoxMultimap<T> createQuadtree(int dims, int maxNodeCapacity, double[] min, double[] max) {
             return QuadTreeRKD0.create(dims, maxNodeCapacity, min, max);
         }
 
+        /**
+         * Create a plain Quadtree.
+         * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param min             Estimated minimum of all coordinates.
+         * @param max             Estimated maximum of all coordinates.
+         * @param align           Whether center and radius should be aligned to powers of two. Aligning considerably
+         *                        reduces risk of precision problems. Recommended: "true".
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param <T>             Value type
+         * @return New Quadtree
+         */
+        static <T> BoxMultimap<T> createQuadtree(double[] min, double[] max, boolean align, int maxNodeCapacity) {
+            return QuadTreeRKD0.create(min, max, align, maxNodeCapacity);
+        }
+
         /**
          * Create a Quadtree with hypercube navigation.
          *
@@ -195,7 +217,8 @@ static <T> BoxMultimap<T> createQuadtreeHC(int dims) {
         }
 
         /**
-         * Create a plain Quadtree.
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
+         * Create a Quadtree with hypercube navigation.
          * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
          * span an area that is somewhat larger rather than smaller than the actual data.
          *
@@ -205,11 +228,32 @@ static <T> BoxMultimap<T> createQuadtreeHC(int dims) {
          * @param max             Estimated maximum of all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC
+         * @deprecated Please use {@link #createQuadtree(double[], double[], boolean, int)}
          */
+        @Deprecated
         static <T> BoxMultimap<T> createQuadtreeHC(int dims, int maxNodeCapacity, double[] min, double[] max) {
             return QuadTreeRKD.create(dims, maxNodeCapacity, min, max);
         }
 
+        /**
+         * Create a Quadtree with hypercube navigation.
+         * Min/max are used to find a good initial root. They do not need to be exact. If possible, min/max should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param min             Estimated minimum of all coordinates.
+         * @param max             Estimated maximum of all coordinates.
+         * @param align           Whether center and radius should be aligned to powers of two. Aligning considerably
+         *                        reduces risk of precision problems. Recommended: "true".
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param <T>             Value type
+         * @return New QuadtreeHC
+         */
+        static <T> BoxMultimap<T> createQuadtreeHC(double[] min, double[] max, boolean align, int maxNodeCapacity) {
+            return QuadTreeRKD.create(min, max, align, maxNodeCapacity);
+        }
+
         /**
          * Create an R*Tree.
          *

src/main/java/org/tinspin/index/PointMap.java

@@ -283,6 +283,7 @@ static <T> PointMap<T> createQuadtreeHC2(int dims) {
         }
 
         /**
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
          * Create a Quadtree with extended hypercube navigation.
          * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
          * span an area that is somewhat larger rather than smaller than the actual data.
@@ -293,7 +294,7 @@ static <T> PointMap<T> createQuadtreeHC2(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC2
-         * @deprecated PLease use {@link #createQuadtreeHC2(double[], double, boolean, int)}
+         * @deprecated Please use {@link #createQuadtreeHC2(double[], double, boolean, int)}
          */
         @Deprecated
         static <T> PointMap<T> createQuadtreeHC2(int dims, int maxNodeCapacity, double[] center, double radius) {

src/main/java/org/tinspin/index/PointMultimap.java

@@ -143,7 +143,7 @@ default PointIteratorKnn<T> queryKnn(double[] center, int k) {
 
     interface Factory {
         /**
-         * Create an array backed PointMap. This is only for testing and rather inefficient for large data sets.
+         * Create an array backed PointMultiMap. This is only for testing and rather inefficient for large data sets.
          *
          * @param dims Number of dimensions.
          * @param size Number of entries.
@@ -307,7 +307,7 @@ static <T> PointMultimap<T> createQuadtreeHC2(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC2
-         * @deprecated PLease use {@link #createQuadtreeHC2(double[], double, boolean, int)}
+         * @deprecated Please use {@link #createQuadtreeHC2(double[], double, boolean, int)}
          */
         @Deprecated
         static <T> PointMultimap<T> createQuadtreeHC2(int dims, int maxNodeCapacity, double[] center, double radius) {

src/main/java/org/tinspin/index/qthypercube/QuadTreeRKD.java

@@ -23,6 +23,7 @@
 import org.tinspin.index.*;
 import org.tinspin.index.qthypercube.QuadTreeKD.QStats;
 import org.tinspin.index.util.BoxIteratorWrapper;
+import org.tinspin.index.util.MathTools;
 import org.tinspin.index.util.StringBuilderLn;
 
 /**
@@ -61,33 +62,74 @@ public static <T> QuadTreeRKD<T> create(int dims, int maxNodeSize) {
 	}
 
 	/**
-	 *
 	 * @param dims Number of dimensions per coordinate, usually 2 or 3
-	 * @param maxNodeSize Maximum node capacity before a split occurs
+	 * @param maxNodeSize Maximum node capacity before a split occurs. Default is 10.
 	 * @param min Estimated global minimum
 	 * @param max Estimated global minimum
 	 * @return New quadtree
 	 * @param <T> Value type
+	 * @deprecated Please use {@link #create(double[], double[], boolean, int)}
 	 */
-	public static <T> QuadTreeRKD<T> create(int dims, int maxNodeSize,
-			double[] min, double[] max) {
+	@Deprecated
+	public static <T> QuadTreeRKD<T> create(int dims, int maxNodeSize, double[] min, double[] max) {
+		return create(min, max, false, maxNodeSize);
+	}
+
+	/**
+	 * @param min Estimated global minimum
+	 * @param max Estimated global minimum
+	 * @param align Whether min and max should be aligned to powers of two. Aligning considerably
+	 *              reduces risk of precision problems. Recommended: "true".
+	 * @param maxNodeSize Maximum node capacity before a split occurs. Default is 10.
+	 * @return New quadtree
+	 * @param <T> Value type
+	 */
+	public static <T> QuadTreeRKD<T> create(double[] min, double[] max, boolean align, int maxNodeSize) {
 		double radius = 0;
-		double[] center = new double[dims];
-		for (int i = 0; i < dims; i++) {
+		double[] center = new double[min.length];
+		for (int i = 0; i < center.length; i++) {
 			center[i] = (max[i]+min[i])/2.0;
 			if (max[i]-min[i]>radius) {
 				radius = max[i]-min[i];
 			}
 		}
-		return create(dims, maxNodeSize, center, radius);
+		return create(center, radius, align, maxNodeSize);
 	}
-
-	public static <T> QuadTreeRKD<T> create(int dims, int maxNodeSize, 
-			double[] center, double radius) {
-		QuadTreeRKD<T> t = new QuadTreeRKD<>(dims, maxNodeSize);
+
+	/**
+	 * WARNING: Unaligned center and radius can cause precision problems.
+	 * @param dims dimensions, usually 2 or 3
+	 * @param maxNodeSize maximum entries per node, default is 10
+	 * @param center center of initial root node
+	 * @param radius radius of initial root node
+	 * @return New quadtree
+	 * @param <T> Value type
+	 * @deprecated Please use {@link #create(double[], double, boolean, int)}
+	 */
+	@Deprecated
+	public static <T> QuadTreeRKD<T> create(int dims, int maxNodeSize, double[] center, double radius) {
+		return create(center, radius, false, maxNodeSize);
+	}
+
+	/**
+	 * Note: This will align center and radius to a power of two before creating a tree.
+	 * @param center center of initial root node
+	 * @param radius radius of initial root node
+	 * @param maxNodeSize maximum entries per node, default is 10
+	 * @param align Whether center and radius should be aligned to powers of two. Aligning considerably
+	 *              reduces risk of precision problems. Recommended: "true".
+	 * @return New quadtree
+	 * @param <T> Value type
+	 */
+	public static <T> QuadTreeRKD<T> create(double[] center, double radius, boolean align, int maxNodeSize) {
+		QuadTreeRKD<T> t = new QuadTreeRKD<>(center.length, maxNodeSize);
 		if (radius <= 0) {
 			throw new IllegalArgumentException("Radius must be > 0 but was " + radius);
 		}
+		if (align) {
+			center = MathTools.floorPowerOfTwoCopy(center);
+			radius = MathTools.ceilPowerOfTwo(radius);
+		}
 		t.root = new QRNode<>(Arrays.copyOf(center, center.length), radius);
 		return t;
 	}
@@ -114,15 +156,12 @@ public void insert(double[] keyL, double[] keyU, T value) {
 	}
 
 	private void initializeRoot(double[] keyL, double[] keyU) {
-		double[] center = new double[dims];
-		double radius = 0;
-		for (int d = 0; d < dims; d++) {
-			center[d] = (keyU[d] + keyL[d]) / 2.0;
-			if (keyU[d]-keyL[d] > radius) {
-				radius = keyU[d] -keyL[d];
-			}
-		}			
-		radius *= 5; //for good measure
+		// Find a power-of-2 center, potentially inside the box entry.
+		double[] center = MathTools.floorPowerOfTwoCopy(keyU);
+		// Get minimum required radius.
+		double dMax = Math.max(MathTools.maxDelta(center, keyU), MathTools.maxDelta(center, keyL));
+		// Make radius a power of two.
+		double radius = MathTools.ceilPowerOfTwo(dMax + QUtil.EPS_MUL);
 		root = new QRNode<>(center, radius);
 	}