diff --git a/java/source/infodynamics/utils/KdTree.java b/java/source/infodynamics/utils/KdTree.java index d1d7908..dea14a2 100755 --- a/java/source/infodynamics/utils/KdTree.java +++ b/java/source/infodynamics/utils/KdTree.java @@ -18,7 +18,9 @@ package infodynamics.utils; +import java.util.Collection; import java.util.PriorityQueue; +import java.util.Vector; /** * K-d tree implementation to be used for fast neighbour searching @@ -82,13 +84,6 @@ public class KdTree extends NearestNeighbourSearcher { */ protected KdTreeNode rootNode = null; - /** - * Calculator for computing the norms for each variable; defaults - * to a max norm. - */ - protected EuclideanUtils normCalculator; - - /** * Protected class to implement nodes of a k-d tree * @@ -131,9 +126,7 @@ public class KdTree extends NearestNeighbourSearcher { * within this data set) */ public KdTree(int[] dimensions, double[][][] data) { - - normCalculator = new EuclideanUtils(normTypeToUse); - + this.originalDataSets = data; int numObservations = data[0].length; @@ -304,15 +297,100 @@ public class KdTree extends NearestNeighbourSearcher { @Override public void setNormType(int normType) { super.setNormType(normType); - normCalculator.setNormToUse(normTypeToUse); } @Override public void setNormType(String normTypeString) { super.setNormType(normTypeString); - normCalculator.setNormToUse(normTypeToUse); } - + + /** + * Computing the configured norm between vectors x1 and x2. + * Adding here instead of using {@link EuclideanUtils#norm(double[], double[])} + * to attempt speed-up. + * Also hoping this method is inlined by the JVM, but haven't checked this. + * + * @param x1 vector of doubles + * @param x2 vector of doubles + * @return the selected norm + */ + public final static double norm(double[] x1, double[] x2, int normToUse) { + double distance = 0.0; + switch (normToUse) { + case EuclideanUtils.NORM_MAX_NORM: + // Inlined from {@link EuclideanUtils}: + for (int d = 0; d < x1.length; d++) { + double difference = x1[d] - x2[d]; + // Take the abs + if (difference < 0) { + difference = -difference; + } + if (difference > distance) { + distance = difference; + } + } + return distance; + // case EuclideanUtils.NORM_EUCLIDEAN_SQUARED: + default: + // Inlined from {@link EuclideanUtils}: + for (int d = 0; d < x1.length; d++) { + double difference = x1[d] - x2[d]; + distance += difference * difference; + } + return distance; + } + } + + /** + * Computing the configured norm between vectors x1 and x2; if + * it becomes clear that norm will be larger than limit, + * then return Double.POSITIVE_INFINITY immediately. + * + *

Adding here instead of using {@link EuclideanUtils#normWithAbort(double[], double[], double)} + * to attempt speed-up. + * Also hoping this method is inlined by the JVM, but haven't checked this.

+ * + * @param x1 vector 1 of doubles + * @param x2 vector 2 of doubles + * @param limit if it becomes clear that norm will be larger than limit, + * then return Double.POSITIVE_INFINITY immediately. + * @param normToUse which norm to use, as defined by {@link #setNormType(int)} + * @return the selected norm + */ + public final static double normWithAbort(double[] x1, double[] x2, + double limit, int normToUse) { + double distance = 0.0; + switch (normToUse) { + case EuclideanUtils.NORM_MAX_NORM: + // Inlined from {@link EuclideanUtils}: + for (int d = 0; d < x1.length; d++) { + double difference = x1[d] - x2[d]; + // Take the abs + if (difference < 0) { + difference = -difference; + } + if (difference > distance) { + if (difference > limit) { + return Double.POSITIVE_INFINITY; + } + distance = difference; + } + } + return distance; + // case EuclideanUtils.NORM_EUCLIDEAN_SQUARED: + default: + // Inlined from {@link EuclideanUtils}: + for (int d = 0; d < x1.length; d++) { + double difference = x1[d] - x2[d]; + distance += difference * difference; + if (distance > limit) { + return Double.POSITIVE_INFINITY; + } + } + return distance; + } + } + @Override public NeighbourNodeData findNearestNeighbour(int sampleIndex) { if (rootNode == null) { @@ -350,7 +428,7 @@ public class KdTree extends NearestNeighbourSearcher { data[node.indexOfThisPoint][actualDim]; double absDistOnThisDim; - if (normCalculator.getNormInUse() == EuclideanUtils.NORM_MAX_NORM) { + if (normTypeToUse == EuclideanUtils.NORM_MAX_NORM) { absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; } else { // norm type is EuclideanUtils#NORM_EUCLIDEAN_SQUARED @@ -369,16 +447,16 @@ public class KdTree extends NearestNeighbourSearcher { // For each of our separate (multivariate) variables, // compute the (specified) norm in that variable's space: if (currentBest == null) { - norms[v] = normCalculator.norm( + norms[v] = norm( originalDataSets[v][sampleIndex], - originalDataSets[v][node.indexOfThisPoint]); + originalDataSets[v][node.indexOfThisPoint], normTypeToUse); } else { // Distance calculation terminates early with Double.POSITIVE_INFINITY // if it is clearly larger than currentBest.distance: - norms[v] = normCalculator.normWithAbort( + norms[v] = normWithAbort( originalDataSets[v][sampleIndex], originalDataSets[v][node.indexOfThisPoint], - currentBest.distance); + currentBest.distance, normTypeToUse); } if (norms[v] > maxNorm) { maxNorm = norms[v]; @@ -472,7 +550,7 @@ public class KdTree extends NearestNeighbourSearcher { double distOnThisDim = data[sampleIndex][actualDim] - data[node.indexOfThisPoint][actualDim]; double absDistOnThisDim; - if (normCalculator.getNormInUse() == EuclideanUtils.NORM_MAX_NORM) { + if (normTypeToUse == EuclideanUtils.NORM_MAX_NORM) { absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; } else { // norm type is EuclideanUtils#NORM_EUCLIDEAN_SQUARED @@ -495,16 +573,16 @@ public class KdTree extends NearestNeighbourSearcher { // For each of our separate (multivariate) variables, // compute the (specified) norm in that variable's space: if (currentKBest.size() < K) { - norms[v] = normCalculator.norm( + norms[v] = norm( originalDataSets[v][sampleIndex], - originalDataSets[v][node.indexOfThisPoint]); + originalDataSets[v][node.indexOfThisPoint], normTypeToUse); } else { // Distance calculation terminates early with Double.POSITIVE_INFINITY // if it is clearly larger than currentBest.distance: - norms[v] = normCalculator.normWithAbort( + norms[v] = normWithAbort( originalDataSets[v][sampleIndex], originalDataSets[v][node.indexOfThisPoint], - furthestCached.distance); + furthestCached.distance, normTypeToUse); } if (norms[v] > maxNorm) { maxNorm = norms[v]; @@ -583,12 +661,15 @@ public class KdTree extends NearestNeighbourSearcher { } return countPointsWithinR(sampleIndex, rootNode, 0, r, true); } - + /** * Count the number of points within radius r of a given sample (sampleIndex), * in the tree rooted at node (which is at the specified level in the tree). * Nearest neighbour function to compare to r is a max norm between the * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -614,7 +695,7 @@ public class KdTree extends NearestNeighbourSearcher { data[node.indexOfThisPoint][actualDim]; double absDistOnThisDim; - if (normCalculator.getNormInUse() == EuclideanUtils.NORM_MAX_NORM) { + if (normTypeToUse == EuclideanUtils.NORM_MAX_NORM) { absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; } else { // norm type is EuclideanUtils#NORM_EUCLIDEAN_SQUARED @@ -634,10 +715,10 @@ public class KdTree extends NearestNeighbourSearcher { double distForVariableV; // Distance calculation terminates early with Double.POSITIVE_INFINITY // if it is clearly larger than r: - distForVariableV = normCalculator.normWithAbort( + distForVariableV = normWithAbort( originalDataSets[v][sampleIndex], originalDataSets[v][node.indexOfThisPoint], - r); + r, normTypeToUse); if ((distForVariableV >= r) && !(allowEqualToR && (distForVariableV == r))) { // We don't fit on this dimension, no point @@ -689,6 +770,145 @@ public class KdTree extends NearestNeighbourSearcher { return count; } + @Override + public Collection findPointsWithinR(int sampleIndex, + double r, boolean allowEqualToR) { + Vector pointsWithinR = new Vector(); + + if (rootNode == null) { + return pointsWithinR; + } + findPointsWithinR(sampleIndex, + rootNode, 0, r, allowEqualToR, pointsWithinR); + + return pointsWithinR; + } + + @Override + public Collection findPointsStrictlyWithinR( + int sampleIndex, double r) { + return findPointsWithinR(sampleIndex, r, false); + } + + @Override + public Collection findPointsWithinOrOnR(int sampleIndex, + double r) { + return findPointsWithinR(sampleIndex, r, true); + } + + /** + * Add to the collection of points within radius r of a given sample (sampleIndex), + * in the tree rooted at node (which is at the specified level in the tree). + * Nearest neighbour function to compare to r is a max norm between the + * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). + * + * @param sampleIndex sample index in the data to find a nearest neighbour + * for + * @param node node to start searching from in the kd-tree. Cannot be null + * @param level which level we're currently at in the tree + * @param r radius within which to count points + * @param allowEqualToR if true, then count points at radius r also, + * otherwise only those strictly within r + * @param pointsWithinR the collection of points to add to + */ + protected void findPointsWithinR(int sampleIndex, + KdTreeNode node, int level, double r, boolean allowEqualToR, + Collection pointsWithinR) { + + // Point to the correct array for the data at this level + int currentDim = level % totalDimensions; + double[][] data = dimensionToArray[currentDim]; + int actualDim = dimensionToArrayIndex[currentDim]; + + // Check the distance on this particular dimension + double distOnThisDim = data[sampleIndex][actualDim] - + data[node.indexOfThisPoint][actualDim]; + + double absDistOnThisDim; + if (normTypeToUse == EuclideanUtils.NORM_MAX_NORM) { + absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; + } else { + // norm type is EuclideanUtils#NORM_EUCLIDEAN_SQUARED + // Track the square distance + absDistOnThisDim = distOnThisDim * distOnThisDim; + } + + if ((node.indexOfThisPoint != sampleIndex) && + ((absDistOnThisDim < r) || + ( allowEqualToR && (absDistOnThisDim == r)))) { + // Preliminary check says we need to compute the full distance + // to use or at least to check if it should be counted. + boolean withinBounds = true; + double[] norms = new double[originalDataSets.length]; + double maxNorm = 0; + for (int v = 0; v < originalDataSets.length; v++) { + // For each of our separate (multivariate) variables, + // compute the (specified) norm in that variable's space: + double distForVariableV; + // Distance calculation terminates early with Double.POSITIVE_INFINITY + // if it is clearly larger than r: + distForVariableV = normWithAbort( + originalDataSets[v][sampleIndex], + originalDataSets[v][node.indexOfThisPoint], + r, normTypeToUse); + if ((distForVariableV >= r) && + !(allowEqualToR && (distForVariableV == r))) { + // We don't fit on this variable, no point + // checking the others: + withinBounds = false; + break; + } + norms[v] = distForVariableV; + if (distForVariableV > maxNorm) { + maxNorm = distForVariableV; + } + } + if (withinBounds) { + // This node gets counted + pointsWithinR.add( + new NeighbourNodeData(node.indexOfThisPoint, + norms, maxNorm)); + } + } + + KdTreeNode closestSubTree = null; + KdTreeNode furthestSubTree = null; + // And translate this to which subtree is closer + if (distOnThisDim < 0) { + // We need to search the left tree + closestSubTree = node.leftTree; + furthestSubTree = node.rightTree; + } else { + // We need to search the right tree + closestSubTree = node.rightTree; + furthestSubTree = node.leftTree; + } + // Update the search on that subtree + if (closestSubTree != null) { + findPointsWithinR(sampleIndex, closestSubTree, + level + 1, r, allowEqualToR, pointsWithinR); + } + if ((absDistOnThisDim < r) || + ( allowEqualToR && (distOnThisDim < 0) && (absDistOnThisDim == r))) { + // It's possible we could have a node within (or on) r + // in the other branch as well, so search there too. + // (Note: we only check furthest subtree in the == case + // when it's allowed + // *if* it's the right subtree, as only the right sub-tree + // can have node with distance in this coordinate *equal* to + // that of the current node -- left subtree must be strictly + // less than the coordinate of the current node, so + // distance to any of those points could not be equal.) + if (furthestSubTree != null) { + findPointsWithinR(sampleIndex, furthestSubTree, + level + 1, r, allowEqualToR, pointsWithinR); + } + } + } + /** * Count the number of points within norms {r1,r2,etc} for each high-level * variable, for a given @@ -696,6 +916,9 @@ public class KdTree extends NearestNeighbourSearcher { * excluded from the search. * Nearest neighbour function to compare to {r1,r2,etc} * for each variable is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -719,6 +942,9 @@ public class KdTree extends NearestNeighbourSearcher { * excluded from the search. * Nearest neighbour function to compare to {r1,r2,etc} * for each variable is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -739,6 +965,9 @@ public class KdTree extends NearestNeighbourSearcher { * excluded from the search. * Nearest neighbour function to compare to {r1,r2,etc} * for each variable is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -760,6 +989,9 @@ public class KdTree extends NearestNeighbourSearcher { * The node itself is excluded from the search. * Nearest neighbour function to compare to {r1,r2,etc} * for each variable is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -786,7 +1018,7 @@ public class KdTree extends NearestNeighbourSearcher { data[node.indexOfThisPoint][actualDim]; double absDistOnThisDim; - if (normCalculator.getNormInUse() == EuclideanUtils.NORM_MAX_NORM) { + if (normTypeToUse == EuclideanUtils.NORM_MAX_NORM) { absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; } else { // norm type is EuclideanUtils#NORM_EUCLIDEAN_SQUARED @@ -806,10 +1038,10 @@ public class KdTree extends NearestNeighbourSearcher { double distForVariableV; // Distance calculation terminates early with Double.POSITIVE_INFINITY // if it is clearly larger than rs[v]: - distForVariableV = normCalculator.normWithAbort( + distForVariableV = normWithAbort( originalDataSets[v][sampleIndex], originalDataSets[v][node.indexOfThisPoint], - rs[v]); + rs[v], normTypeToUse); if ((distForVariableV >= rs[v]) && !(allowEqualToR && (distForVariableV == rs[v]))) { // We don't fit on this dimension, no point diff --git a/java/source/infodynamics/utils/NearestNeighbourSearcher.java b/java/source/infodynamics/utils/NearestNeighbourSearcher.java index 603829f..524d35f 100755 --- a/java/source/infodynamics/utils/NearestNeighbourSearcher.java +++ b/java/source/infodynamics/utils/NearestNeighbourSearcher.java @@ -18,6 +18,7 @@ package infodynamics.utils; +import java.util.Collection; import java.util.PriorityQueue; /** @@ -186,6 +187,9 @@ public abstract class NearestNeighbourSearcher { * excluded from the search. * Nearest neighbour function to compare to r is a max norm between the * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -197,12 +201,36 @@ public abstract class NearestNeighbourSearcher { public abstract int countPointsWithinR(int sampleIndex, double r, boolean allowEqualToR); + /** + * Return the collection of points within norm r for a given + * sample index in the data set. The node itself is + * excluded from the search. + * Nearest neighbour function to compare to r is a max norm between the + * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). + * + * @param sampleIndex sample index in the data to find a nearest neighbour + * for + * @param r radius within which to count points + * @param allowEqualToR if true, then count points at radius r also, + * otherwise only those strictly within r + * @return the collection of points within r. + */ + public abstract Collection findPointsWithinR( + int sampleIndex, double r, + boolean allowEqualToR); + /** * Count the number of points strictly within norm r for a given * sample index in the data set. The node itself is * excluded from the search. * Nearest neighbour function to compare to r is a max norm between the * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -211,12 +239,32 @@ public abstract class NearestNeighbourSearcher { */ public abstract int countPointsStrictlyWithinR(int sampleIndex, double r); + /** + * Return the collection of points strictly within norm r for a given + * sample index in the data set. The node itself is + * excluded from the search. + * Nearest neighbour function to compare to r is a max norm between the + * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). + * + * @param sampleIndex sample index in the data to find a nearest neighbour + * for + * @param r radius within which to count points + * @return the collection of points within r. + */ + public abstract Collection findPointsStrictlyWithinR(int sampleIndex, double r); + /** * Count the number of points within or at norm r for a given * sample index in the data set. The node itself is * excluded from the search. * Nearest neighbour function to compare to r is a max norm between the * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * * @param sampleIndex sample index in the data to find a nearest neighbour * for @@ -224,4 +272,21 @@ public abstract class NearestNeighbourSearcher { * @return the count of points within or on r. */ public abstract int countPointsWithinOrOnR(int sampleIndex, double r); + + /** + * Return the collection of points within or at norm r for a given + * sample index in the data set. The node itself is + * excluded from the search. + * Nearest neighbour function to compare to r is a max norm between the + * high-level variables, with norm for each variable being the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). + * + * @param sampleIndex sample index in the data to find a nearest neighbour + * for + * @param r radius within which to count points + * @return the collection of points within or on r. + */ + public abstract Collection findPointsWithinOrOnR(int sampleIndex, double r); } diff --git a/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java b/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java index b26a2a8..267bc1d 100755 --- a/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java +++ b/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java @@ -18,7 +18,9 @@ package infodynamics.utils; +import java.util.Collection; import java.util.PriorityQueue; +import java.util.Vector; /** @@ -80,17 +82,18 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher } /** - * Computed the configured norm between the unidimensional variables + * Computed the configured norm between the unidimensional variables. + * Hoping this method is inlined by the JVM, but haven't checked this. * * @param x1 data point 1 * @param x2 data point 2 * @return the norm */ - protected double norm(double x1, double x2) { + protected static double norm(double x1, double x2, int normTypeToUse) { switch (normTypeToUse) { case EuclideanUtils.NORM_MAX_NORM: return Math.abs(x1-x2); - case EuclideanUtils.NORM_EUCLIDEAN_SQUARED: + // case EuclideanUtils.NORM_EUCLIDEAN_SQUARED: default: double difference = x1 - x2; return difference * difference; @@ -111,7 +114,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher // Assumes we have more than 1 data point -- this is // checked in the constructor for us. double theNorm = norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[1]]); + originalDataSet[sortedArrayIndices[1]], normTypeToUse); return new NeighbourNodeData(sortedArrayIndices[1], new double[] {theNorm}, theNorm); } else if (indexInSortedArray == numObservations - 1) { @@ -119,15 +122,18 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher // Assumes we have more than 1 data point -- this is // checked in the constructor for us. double theNorm = norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[numObservations - 2]]); + originalDataSet[sortedArrayIndices[numObservations - 2]], + normTypeToUse); return new NeighbourNodeData(sortedArrayIndices[numObservations - 2], new double[] {theNorm}, theNorm); } else { // We need to check candidates on both sides of the data point: double normAbove = norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[indexInSortedArray+1]]); + originalDataSet[sortedArrayIndices[indexInSortedArray+1]], + normTypeToUse); double normBelow = norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[indexInSortedArray-1]]); + originalDataSet[sortedArrayIndices[indexInSortedArray-1]], + normTypeToUse); if (normAbove < normBelow) { return new NeighbourNodeData(sortedArrayIndices[indexInSortedArray+1], new double[] {normAbove}, normAbove); @@ -168,11 +174,13 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher double normAbove = (upperCandidate == -1) ? Double.POSITIVE_INFINITY : norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[upperCandidate]]); + originalDataSet[sortedArrayIndices[upperCandidate]], + normTypeToUse); double normBelow = (lowerCandidate == -1) ? Double.POSITIVE_INFINITY : norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[lowerCandidate]]); + originalDataSet[sortedArrayIndices[lowerCandidate]], + normTypeToUse); NeighbourNodeData nextNearest; if (normAbove < normBelow) { nextNearest = new NeighbourNodeData(sortedArrayIndices[upperCandidate], @@ -199,6 +207,9 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher * sample index in the data set. The node itself is * excluded from the search. * Nearest neighbour function to compare to r is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). */ public int countPointsWithinR(int sampleIndex, double r, boolean allowEqualToR) { int count = 0; @@ -207,7 +218,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher // Check the points with smaller data values first: for (int i = indexInSortedArray - 1; i >= 0; i--) { double theNorm = norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[i]]); + originalDataSet[sortedArrayIndices[i]], normTypeToUse); if ((allowEqualToR && (theNorm <= r)) || (!allowEqualToR && (theNorm < r))) { count++; @@ -219,7 +230,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher // Next check the points with larger data values: for (int i = indexInSortedArray + 1; i < numObservations; i++) { double theNorm = norm(originalDataSet[sampleIndex], - originalDataSet[sortedArrayIndices[i]]); + originalDataSet[sortedArrayIndices[i]], normTypeToUse); if ((allowEqualToR && (theNorm <= r)) || (!allowEqualToR && (theNorm < r))) { count++; @@ -231,24 +242,101 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher return count; } + /** + * Count the number of points within norm r for a given + * sample index in the data set. The node itself is + * excluded from the search. + * Nearest neighbour function to compare to r is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). + */ + public Collection findPointsWithinR(int sampleIndex, double r, boolean allowEqualToR) { + Vector pointsWithinR = new Vector(); + // Find where this node sits in the sorted array: + int indexInSortedArray = indicesInSortedArray[sampleIndex]; + // Check the points with smaller data values first: + for (int i = indexInSortedArray - 1; i >= 0; i--) { + double theNorm = norm(originalDataSet[sampleIndex], + originalDataSet[sortedArrayIndices[i]], normTypeToUse); + if ((allowEqualToR && (theNorm <= r)) || + (!allowEqualToR && (theNorm < r))) { + pointsWithinR.add( + new NeighbourNodeData(sortedArrayIndices[i], + new double[] {theNorm}, theNorm)); + continue; + } + // Else no point checking further points + break; + } + // Next check the points with larger data values: + for (int i = indexInSortedArray + 1; i < numObservations; i++) { + double theNorm = norm(originalDataSet[sampleIndex], + originalDataSet[sortedArrayIndices[i]], normTypeToUse); + if ((allowEqualToR && (theNorm <= r)) || + (!allowEqualToR && (theNorm < r))) { + pointsWithinR.add( + new NeighbourNodeData(sortedArrayIndices[i], + new double[] {theNorm}, theNorm)); + continue; + } + // Else no point checking further points + break; + } + return pointsWithinR; + } + /** * Count the number of points strictly within norm r for a given * sample index in the data set. The node itself is * excluded from the search. * Nearest neighbour function to compare to r is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). * */ public int countPointsStrictlyWithinR(int sampleIndex, double r) { return countPointsWithinR(sampleIndex, r, false); } + /** + * Return a collection of points strictly within norm r for a given + * sample index in the data set. The node itself is + * excluded from the search. + * Nearest neighbour function to compare to r is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). + * + */ + public Collection findPointsStrictlyWithinR(int sampleIndex, double r) { + return findPointsWithinR(sampleIndex, r, false); + } + /** * Count the number of points within or at norm r for a given * sample index in the data set. The node itself is * excluded from the search. * Nearest neighbour function to compare to r is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). */ public int countPointsWithinOrOnR(int sampleIndex, double r) { return countPointsWithinR(sampleIndex, r, true); } + + /** + * Return a collection of points within or at norm r for a given + * sample index in the data set. The node itself is + * excluded from the search. + * Nearest neighbour function to compare to r is the specified norm. + * (If {@link EuclideanUtils#NORM_EUCLIDEAN} was selected, then the supplied + * r should be the required Euclidean norm squared, since we switch it + * to {@link EuclideanUtils#NORM_EUCLIDEAN_SQUARED} internally). + */ + public Collection findPointsWithinOrOnR(int sampleIndex, double r) { + return findPointsWithinR(sampleIndex, r, true); + } }