=%.3f - $d/k=%.3f = %.3f" +
" (<1/n_yz>=%.3f, <1/n_xz>=%.3f)\n",
MathsUtils.digamma(k), averageDiGammas,
averageInverseCountInJointXZ + averageInverseCountInJointYZ,
- (2.0 / (double)k), averageMeasure,
+ dimensionsCond > 0 ? 2 : 1,
+ inverseKTerm, averageMeasure,
averageInverseCountInJointYZ, averageInverseCountInJointXZ);
}
double[] result = new double[1];
@@ -619,6 +662,40 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
protected abstract double[] partialComputeFromObservations(
int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception;
+ /**
+ * Protected method to be used internally for threaded implementations.
+ * This method implements the guts of each Kraskov algorithm, computing the number of
+ * nearest neighbours in each dimension for a sub-set of the data points.
+ * It is intended to be called by one thread to work on that specific
+ * sub-set of the data.
+ * In particular, this method differs from {@link #partialComputeFromObservations(int, int, boolean)}
+ * because it operates on a new set of observations (using the old set of observations for
+ * constructing the search spaces and PDFs)
+ *
+ * The method returns:
+ * - for average conditional MIs (returnLocals == false), the relevant sums of digamma(n_{xz}), digamma(n_{yz})
+ * and digamma(n_z)
+ * for a partial set of the observations
+ * - for local conditional MIs (returnLocals == true), the array of local conditional MI values
+ *
+ *
+ * @param startTimePoint start time for the partial set we examine
+ * @param numTimePoints number of time points (including startTimePoint to examine)
+ * @param newVar1Observations new time series of observations for variable 1
+ * @param newVar2Observations new time series of observations for variable 1
+ * @param newCondObservations new time series of observations for variable 1
+ * @param returnLocals whether to return an array or local values, or else
+ * sums of these values
+ * @return an array of the relevant sum of digamma(n_xz+1) and digamma(n_yz+1) and digamma(n_z), then
+ * sum of n_xz, n_yz, n_z and for algorithm 2 only, sum of 1/n_xz and 1/n_yz
+ * (these latter five are only for debugging purposes).
+ * @throws Exception
+ */
+ protected abstract double[] partialComputeFromNewObservations(
+ int startTimePoint, int numTimePoints,
+ double[][] newVar1Observations, double[][] newVar2Observations, double[][] newCondObservations,
+ boolean returnLocals) throws Exception;
+
/**
* Private class to handle multi-threading of the Kraskov algorithms.
* Each instance calls partialComputeFromObservations()
@@ -633,6 +710,7 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
protected ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc;
protected int myStartTimePoint;
protected int numberOfTimePoints;
+ protected double[][][] newObservations;
protected boolean computeLocals;
protected double[] returnValues = null;
@@ -649,11 +727,13 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
public CondMiKraskovThreadRunner(
ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc,
int myStartTimePoint, int numberOfTimePoints,
+ double[][][] newObservations,
boolean computeLocals) {
this.condMiCalc = condMiCalc;
this.myStartTimePoint = myStartTimePoint;
this.numberOfTimePoints = numberOfTimePoints;
this.computeLocals = computeLocals;
+ this.newObservations = newObservations;
}
/**
@@ -676,7 +756,16 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
*/
public void run() {
try {
- returnValues = condMiCalc.partialComputeFromObservations(myStartTimePoint, numberOfTimePoints, computeLocals);
+ if (newObservations == null) {
+ // Computing on existing observations
+ returnValues = condMiCalc.partialComputeFromObservations(myStartTimePoint, numberOfTimePoints, computeLocals);
+ } else {
+ // Computing on new observations
+ returnValues = condMiCalc.partialComputeFromNewObservations(
+ myStartTimePoint, numberOfTimePoints,
+ newObservations[0], newObservations[1],
+ newObservations[2], computeLocals);
+ }
} catch (Exception e) {
// Store the exception for later retrieval
problem = e;
diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java
index bfbfa5a..677bc43 100755
--- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java
+++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java
@@ -18,12 +18,18 @@
package infodynamics.measures.continuous.kraskov;
+import java.util.Arrays;
import java.util.Calendar;
import java.util.PriorityQueue;
import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariate;
+import infodynamics.utils.EmpiricalMeasurementDistribution;
+import infodynamics.utils.FirstIndexComparatorDouble;
+import infodynamics.utils.KdTree;
import infodynamics.utils.MathsUtils;
+import infodynamics.utils.MatrixUtils;
import infodynamics.utils.NeighbourNodeData;
+import infodynamics.utils.UnivariateNearestNeighbourSearcher;
/**
* Computes the differential conditional mutual information of two multivariate
@@ -99,7 +105,20 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
// First element in the PQ is the kth NN,
// and epsilon = kthNnData.distance
NeighbourNodeData kthNnData = nnPQ.poll();
-
+ /*
+ if (debug) {
+ System.out.print("t = " + t + " : for data point: ");
+ MatrixUtils.printArray(System.out, var1Observations[t]);
+ MatrixUtils.printArray(System.out, var2Observations[t]);
+ if (dimensionsCond > 0) {
+ MatrixUtils.printArray(System.out, condObservations[t]);
+ } else {
+ System.out.print("[]");
+ }
+ System.out.printf("t=%d: K=%d NNs found at range %.5f (point %d)\n", t, k, kthNnData.distance, kthNnData.sampleIndex);
+ }
+ */
+
// Now count the points in the conditional space, and
// the var1-conditional and var2-conditional spaces.
// We have 3 coded options for how to do this:
@@ -112,7 +131,7 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
int n_yz = kdTreeVar2Conditional.countPointsStrictlyWithinR(
t, kthNnData.distance, dynCorrExclTime);
int n_z = nnSearcherConditional.countPointsStrictlyWithinR(
- t, kthNnData.distance, dynCorrExclTime);
+ t, kthNnData.distance, dynCorrExclTime);
*/ // end option A
/* Option B --
@@ -151,8 +170,10 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
if (debug) {
methodStartTime = Calendar.getInstance().getTimeInMillis();
}
- nnSearcherConditional.findPointsWithinR(t, kthNnData.distance, dynCorrExclTime,
+ if (dimensionsCond > 0) {
+ nnSearcherConditional.findPointsWithinR(t, kthNnData.distance, dynCorrExclTime,
false, isWithinRForConditionals, indicesWithinRForConditionals);
+ }
if (debug) {
conditionalTime += Calendar.getInstance().getTimeInMillis() -
methodStartTime;
@@ -160,15 +181,26 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
}
// 2. Then compute n_xz and n_yz harnessing our knowledge of
// which points qualified for the conditional already:
- // Don't need to supply dynCorrExclTime in the following, because only
+ // Don't need to supply dynCorrExclTime in most of the following, because only
// points outside of it have been included in isWithinRForConditionals
int n_xz;
- if (dimensionsVar1 > 1) {
- n_xz = kdTreeVar1Conditional.countPointsWithinR(t, kthNnData.distance,
- false, 1, isWithinRForConditionals);
- } else { // Generally faster to search only the marginal space if it is univariate
- n_xz = uniNNSearcherVar1.countPointsWithinR(t, kthNnData.distance,
- false, isWithinRForConditionals);
+ if (dimensionsCond == 0) {
+ // We're really only counting whether the x space qualifies
+ if (dimensionsVar1 > 1) {
+ n_xz = kdTreeVar1Conditional.countPointsStrictlyWithinR(
+ t, kthNnData.distance, dynCorrExclTime);
+ } else {
+ n_xz = uniNNSearcherVar1.countPointsWithinR(t, kthNnData.distance,
+ dynCorrExclTime, false);
+ }
+ } else {
+ if (dimensionsVar1 > 1) {
+ n_xz = kdTreeVar1Conditional.countPointsWithinR(t, kthNnData.distance,
+ false, 1, isWithinRForConditionals);
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_xz = uniNNSearcherVar1.countPointsWithinR(t, kthNnData.distance,
+ false, isWithinRForConditionals);
+ }
}
if (debug) {
conditionalXTime += Calendar.getInstance().getTimeInMillis() -
@@ -176,12 +208,23 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
methodStartTime = Calendar.getInstance().getTimeInMillis();
}
int n_yz;
- if (dimensionsVar2 > 1) {
- n_yz = kdTreeVar2Conditional.countPointsWithinR(t, kthNnData.distance,
- false, 1, isWithinRForConditionals);
- } else { // Generally faster to search only the marginal space if it is univariate
- n_yz = uniNNSearcherVar2.countPointsWithinR(t, kthNnData.distance,
- false, isWithinRForConditionals);
+ if (dimensionsCond == 0) {
+ // We're really only counting whether the y space qualifies
+ if (dimensionsVar2 > 1) {
+ n_yz = kdTreeVar2Conditional.countPointsStrictlyWithinR(
+ t, kthNnData.distance, dynCorrExclTime);
+ } else {
+ n_yz = uniNNSearcherVar2.countPointsWithinR(t, kthNnData.distance,
+ dynCorrExclTime, false);
+ }
+ } else {
+ if (dimensionsVar2 > 1) {
+ n_yz = kdTreeVar2Conditional.countPointsWithinR(t, kthNnData.distance,
+ false, 1, isWithinRForConditionals);
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_yz = uniNNSearcherVar2.countPointsWithinR(t, kthNnData.distance,
+ false, isWithinRForConditionals);
+ }
}
if (debug) {
conditionalYTime += Calendar.getInstance().getTimeInMillis() -
@@ -189,8 +232,15 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
}
// 3. Finally, reset our boolean array for its next use while we count n_z:
int n_z;
- for (n_z = 0; indicesWithinRForConditionals[n_z] != -1; n_z++) {
- isWithinRForConditionals[indicesWithinRForConditionals[n_z]] = false;
+ if (dimensionsCond == 0) {
+ n_z = totalObservations - 1; // - 1 to remove the point itself.
+ // Note: This doesn't respect the dynamic correlation exclusion, but
+ // does align with a mutual information calculation here (to give the digamma(N) term).
+ // No need to reset the boolean array
+ } else {
+ for (n_z = 0; indicesWithinRForConditionals[n_z] != -1; n_z++) {
+ isWithinRForConditionals[indicesWithinRForConditionals[n_z]] = false;
+ }
}
// end option C
@@ -206,11 +256,16 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
if (returnLocals) {
localCondMi[t-startTimePoint] = digammaK - digammaNxzPlusOne - digammaNyzPlusOne + digammaNzPlusOne;
if (debug) {
- System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, local=%.4f," +
+ System.out.printf("t=%d, r=%.5f, n_xz=%d, n_yz=%d, n_z=%d, local=%.4f," +
" digamma(n_xz+1)=%.5f, digamma(n_yz+1)=%.5f, digamma(n_z+1)=%.5f, \n",
- t, n_xz, n_yz, n_z, localCondMi[t-startTimePoint],
+ t, kthNnData.distance, n_xz, n_yz, n_z, localCondMi[t-startTimePoint],
digammaNxzPlusOne, digammaNyzPlusOne, digammaNzPlusOne);
}
+ } else if (debug) {
+ System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d," +
+ " digamma(n_xz+1)=%.5f, digamma(n_yz+1)=%.5f, digamma(n_z+1)=%.5f, \n",
+ t, n_xz, n_yz, n_z,
+ digammaNxzPlusOne, digammaNyzPlusOne, digammaNzPlusOne);
}
}
@@ -244,4 +299,214 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
return results;
}
}
+
+ /**
+ * Compute the local conditional mutual information values for a new set of
+ * observations, where the PDFs are constructed from the observations added earlier.
+ *
+ * @param startTimePoint
+ * @param numTimePoints
+ * @param newStates1
+ * @param newStates2
+ * @param newCondStates
+ * @return
+ * @throws Exception
+ */
+ @Override
+ protected double[] partialComputeFromNewObservations(int startTimePoint,
+ int numTimePoints, double[][] newStates1,
+ double[][] newStates2, double[][] newCondStates, boolean returnLocals) throws Exception {
+
+ double startTime = Calendar.getInstance().getTimeInMillis();
+
+ double[] localCondMi = null;
+ if (returnLocals) {
+ localCondMi = new double[numTimePoints];
+ }
+
+ // Count the average number of points within eps_xz and eps_yz and eps_z of each point
+ double sumDiGammas = 0;
+ double sumNxz = 0;
+ double sumNyz = 0;
+ double sumNz = 0;
+
+ long knnTime = 0, conditionalTime = 0,
+ conditionalXTime = 0, conditionalYTime = 0;
+
+ // Arrays used for fast searching on conditionals with a marginal:
+ boolean[] isWithinRForConditionals = new boolean[totalObservations];
+ int[] indicesWithinRForConditionals = new int[totalObservations+1];
+
+ for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
+ // Compute eps for this time step by
+ // finding the kth closest neighbour for point t:
+ long methodStartTime = Calendar.getInstance().getTimeInMillis();
+ PriorityQueue nnPQ =
+ kdTreeJoint.findKNearestNeighbours(k,
+ new double[][] {newStates1[t], newStates2[t], newCondStates[t]});
+ knnTime += Calendar.getInstance().getTimeInMillis() -
+ methodStartTime;
+ // First element in the PQ is the kth NN,
+ // and epsilon = kthNnData.distance
+ NeighbourNodeData kthNnData = nnPQ.poll();
+ if (debug) {
+ System.out.print("t = " + t + " : for data point: ");
+ MatrixUtils.printArray(System.out, newStates1[t]);
+ MatrixUtils.printArray(System.out, newStates2[t]);
+ if (dimensionsCond > 0) {
+ MatrixUtils.printArray(System.out, newCondStates[t]);
+ } else {
+ System.out.print("[]");
+ }
+ System.out.printf("t=%d : K=%d NNs found at range %.5f (point %d)\n", t, k, kthNnData.distance, kthNnData.sampleIndex);
+ }
+
+ // Now count the points in the conditional space, and
+ // the var1-conditional and var2-conditional spaces.
+ // We use Option C determined above to do this --
+ // Identify the points satisfying the conditional criteria, then use
+ // the knowledge of which points made this cut to speed up the searching
+ // in the conditional-marginal spaces:
+ // 1. Identify the n_z points within the conditional boundaries:
+ if (debug) {
+ methodStartTime = Calendar.getInstance().getTimeInMillis();
+ }
+ if (dimensionsCond > 0) {
+ nnSearcherConditional.findPointsWithinR(kthNnData.distance,
+ new double[][] {newCondStates[t]},
+ false, isWithinRForConditionals, indicesWithinRForConditionals);
+ }
+ if (debug) {
+ conditionalTime += Calendar.getInstance().getTimeInMillis() -
+ methodStartTime;
+ methodStartTime = Calendar.getInstance().getTimeInMillis();
+ }
+ // 2. Then compute n_xz and n_yz harnessing our knowledge of
+ // which points qualified for the conditional already:
+ // Don't need to supply dynCorrExclTime in the following, because
+ // it's not relevant for the new points
+ int n_xz;
+ if (dimensionsCond == 0) {
+ // We're really only counting whether the x space qualifies
+ if (dimensionsVar1 > 1) {
+ n_xz = kdTreeVar1Conditional.countPointsWithinR(
+ new double[][] {newStates1[t], newCondStates[t]},
+ kthNnData.distance, false);
+ } else {
+ n_xz = uniNNSearcherVar1.countPointsWithinR(
+ new double[][] {newStates1[t]}, kthNnData.distance, false);
+ }
+ } else {
+ if (dimensionsVar1 > 1) {
+ n_xz = kdTreeVar1Conditional.countPointsWithinR(kthNnData.distance,
+ new double[][] {newStates1[t], newCondStates[t]},
+ false, 1, isWithinRForConditionals);
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_xz = uniNNSearcherVar1.countPointsWithinR(
+ new double[][] {newStates1[t]}, kthNnData.distance,
+ false, isWithinRForConditionals);
+ }
+ }
+
+ if (debug) {
+ conditionalXTime += Calendar.getInstance().getTimeInMillis() -
+ methodStartTime;
+ methodStartTime = Calendar.getInstance().getTimeInMillis();
+ }
+ int n_yz;
+ if (dimensionsCond == 0) {
+ // We're really only counting whether the y space qualifies
+ if (dimensionsVar2 > 1) {
+ n_yz = kdTreeVar2Conditional.countPointsWithinR(
+ new double[][] {newStates2[t], newCondStates[t]},
+ kthNnData.distance, false);
+ } else {
+ n_yz = uniNNSearcherVar2.countPointsWithinR(
+ new double[][] {newStates2[t]},
+ kthNnData.distance, false);
+ }
+ } else {
+ if (dimensionsVar2 > 1) {
+ n_yz = kdTreeVar2Conditional.countPointsWithinR(kthNnData.distance,
+ new double[][] {newStates2[t], newCondStates[t]},
+ false, 1, isWithinRForConditionals);
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_yz = uniNNSearcherVar2.countPointsWithinR(
+ new double[][] {newStates2[t]}, kthNnData.distance,
+ false, isWithinRForConditionals);
+ }
+ }
+ if (debug) {
+ conditionalYTime += Calendar.getInstance().getTimeInMillis() -
+ methodStartTime;
+ }
+ // 3. Finally, reset our boolean array for its next use while we count n_z:
+ int n_z;
+ if (dimensionsCond == 0) {
+ n_z = totalObservations - 1; // - 1 to remove the point itself.
+ // no need to reset the boolean array
+ } else {
+ for (n_z = 0; indicesWithinRForConditionals[n_z] != -1; n_z++) {
+ isWithinRForConditionals[indicesWithinRForConditionals[n_z]] = false;
+ }
+ }
+ // end option C
+
+ sumNxz += n_xz;
+ sumNyz += n_yz;
+ sumNz += n_z;
+ // And take the digammas:
+ double digammaNxzPlusOne = MathsUtils.digamma(n_xz+1);
+ double digammaNyzPlusOne = MathsUtils.digamma(n_yz+1);
+ double digammaNzPlusOne = MathsUtils.digamma(n_z+1);
+ sumDiGammas += digammaNzPlusOne - digammaNxzPlusOne - digammaNyzPlusOne;
+
+ if (returnLocals) {
+ localCondMi[t-startTimePoint] = digammaK - digammaNxzPlusOne - digammaNyzPlusOne + digammaNzPlusOne;
+ if (debug) {
+ System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, local=%.4f," +
+ " digamma(n_xz+1)=%.5f, digamma(n_yz+1)=%.5f, digamma(n_z+1)=%.5f, \n",
+ t, n_xz, n_yz, n_z, localCondMi[t-startTimePoint],
+ digammaNxzPlusOne, digammaNyzPlusOne, digammaNzPlusOne);
+ }
+ } else if (debug) {
+ double localValue = digammaK - digammaNxzPlusOne - digammaNyzPlusOne + digammaNzPlusOne;
+ System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, local=%.4f," +
+ " digamma(n_xz+1)=%.5f, digamma(n_yz+1)=%.5f, digamma(n_z+1)=%.5f, \n",
+ t, n_xz, n_yz, n_z, localValue,
+ digammaNxzPlusOne, digammaNyzPlusOne, digammaNzPlusOne);
+ }
+ }
+
+ if (debug) {
+ Calendar rightNow2 = Calendar.getInstance();
+ long endTime = rightNow2.getTimeInMillis();
+ System.out.println("Subset " + startTimePoint + ":" +
+ (startTimePoint + numTimePoints) + " Calculation time: " +
+ ((endTime - startTime)/1000.0) + " sec" );
+ System.out.println("Total exec times for: ");
+ System.out.println("\tknn search: " + (knnTime/1000.0));
+ System.out.println("\tz search: " + (conditionalTime/1000.0));
+ System.out.println("\tzx search: " + (conditionalXTime/1000.0));
+ System.out.println("\tzy search: " + (conditionalYTime/1000.0));
+ System.out.printf("%d:%d -- Returning: %.4f, %.4f, %.4f, %.4f\n",
+ startTimePoint, (startTimePoint + numTimePoints),
+ sumDiGammas, sumNxz, sumNyz, sumNz);
+ }
+
+ // Select what to return:
+ if (returnLocals) {
+ return localCondMi;
+ } else {
+ // Pad return array with two values, to allow compatibility in
+ // return length with algorithm 2
+ double[] results = new double[6];
+ results[0] = sumDiGammas;
+ results[1] = sumNxz;
+ results[2] = sumNyz;
+ results[3] = sumNz;
+ return results;
+ }
+ }
+
}
diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java
index cc3da9b..e508f37 100755
--- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java
+++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java
@@ -79,7 +79,15 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2
}
// Constants:
- double twoOnK = 2.0 / (double) k;
+ double inverseKTerm;
+ if (dimensionsCond > 0) {
+ // We will add in the 2/K term as usual
+ inverseKTerm = 2.0 / (double) k;
+ } else {
+ // We will only add in 1/K term since without a conditional there is
+ // technically one less variable in the full joint space
+ inverseKTerm = 1.0 / (double) k;
+ }
// Count the average number of points within eps_xz, eps_yz and eps_z
double sumDiGammas = 0;
@@ -136,36 +144,68 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2
// the knowledge of which points made this cut to speed up the searching
// in the conditional-marginal spaces:
// 1. Identify the n_z points within the conditional boundaries:
- nnSearcherConditional.findPointsWithinR(t, eps_z, dynCorrExclTime,
+ if (dimensionsCond > 0) {
+ nnSearcherConditional.findPointsWithinR(t, eps_z, dynCorrExclTime,
true, isWithinRForConditionals, indicesWithinRForConditionals);
+ }
// 2. Then compute n_xz and n_yz harnessing our knowledge of
// which points qualified for the conditional already:
// Don't need to supply dynCorrExclTime in the following, because only
// points outside of it have been included in isWithinRForConditionals
int n_xz;
- if (dimensionsVar1 > 1) {
- // Check only the x variable against eps_x, use existing results for z
- n_xz = kdTreeVar1Conditional.countPointsWithinRs(t,
- new double[] {eps_x, eps_z},
- true, 1, isWithinRForConditionals);
- } else { // Generally faster to search only the marginal space if it is univariate
- n_xz = uniNNSearcherVar1.countPointsWithinR(t, eps_x,
- true, isWithinRForConditionals);
+ if (dimensionsCond == 0) {
+ // We're really only counting whether the x space qualifies
+ if (dimensionsVar1 > 1) {
+ n_xz = kdTreeVar1Conditional.countPointsWithinOrOnR(
+ t, eps_x, dynCorrExclTime); // No need to pass eps_z in
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_xz = uniNNSearcherVar1.countPointsWithinOrOnR(t, eps_x,
+ dynCorrExclTime);
+ }
+ } else {
+ if (dimensionsVar1 > 1) {
+ // Check only the x variable against eps_x, use existing results for z
+ n_xz = kdTreeVar1Conditional.countPointsWithinRs(t,
+ new double[] {eps_x, eps_z},
+ true, 1, isWithinRForConditionals);
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_xz = uniNNSearcherVar1.countPointsWithinR(t, eps_x,
+ true, isWithinRForConditionals);
+ }
}
int n_yz;
- if (dimensionsVar2 > 1) {
- // Check only the y variable against eps_y, use existing results for z
- n_yz = kdTreeVar2Conditional.countPointsWithinRs(t,
- new double[] {eps_y, eps_z},
- true, 1, isWithinRForConditionals);
- } else { // Generally faster to search only the marginal space if it is univariate
- n_yz = uniNNSearcherVar2.countPointsWithinR(t, eps_y,
- true, isWithinRForConditionals);
+ if (dimensionsCond == 0) {
+ if (dimensionsVar2 > 1) {
+ // Check only the y variable against eps_y, use existing results for z
+ n_yz = kdTreeVar2Conditional.countPointsWithinOrOnR(
+ t, eps_y, dynCorrExclTime); // No need to pass eps_z in
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_yz = uniNNSearcherVar2.countPointsWithinOrOnR(t, eps_y,
+ dynCorrExclTime);
+ }
+ } else {
+ // We're really only counting whether the y space qualifies
+ if (dimensionsVar2 > 1) {
+ // Check only the y variable against eps_y, use existing results for z
+ n_yz = kdTreeVar2Conditional.countPointsWithinRs(t,
+ new double[] {eps_y, eps_z},
+ true, 1, isWithinRForConditionals);
+ } else { // Generally faster to search only the marginal space if it is univariate
+ n_yz = uniNNSearcherVar2.countPointsWithinR(t, eps_y,
+ true, isWithinRForConditionals);
+ }
}
// 3. Finally, reset our boolean array for its next use while we count n_z:
int n_z;
- for (n_z = 0; indicesWithinRForConditionals[n_z] != -1; n_z++) {
- isWithinRForConditionals[indicesWithinRForConditionals[n_z]] = false;
+ if (dimensionsCond == 0) {
+ n_z = totalObservations; // Dropping - 1 to remove the point itself, and
+ // Note: This doesn't respect the dynamic correlation exclusion, but this
+ // does align with a mutual information calculation here (to give the digamma(N) term).
+ // No need to reset the boolean array
+ } else {
+ for (n_z = 0; indicesWithinRForConditionals[n_z] != -1; n_z++) {
+ isWithinRForConditionals[indicesWithinRForConditionals[n_z]] = false;
+ }
}
// end option C
@@ -177,15 +217,25 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2
double digammaNxz = MathsUtils.digamma(n_xz);
double digammaNyz = MathsUtils.digamma(n_yz);
double digammaNz = MathsUtils.digamma(n_z);
- double invN_xz = 1.0/(double) n_xz;
- double invN_yz = 1.0/(double) n_yz;
+ double invN_xz, invN_yz;
+ if (dimensionsCond > 0) {
+ // We will add in the inverse neighbour counts as usual
+ invN_xz = 1.0/(double) n_xz;
+ invN_yz = 1.0/(double) n_yz;
+ } else {
+ // We do not add in the inverse neighbour counts, because they
+ // are technically coming from a single variable neighbour search,
+ // and so there is a 0/n_xz and 0/n_yz contribution
+ invN_xz = 0.0;
+ invN_yz = 0.0;
+ }
sumInverseCountInJointXZ += invN_xz;
sumInverseCountInJointYZ += invN_yz;
double contributionDigammas = digammaNz - digammaNxz - digammaNyz;
sumDiGammas += contributionDigammas;
if (returnLocals) {
- localCondMi[t-startTimePoint] = digammaK - twoOnK + contributionDigammas + invN_xz + invN_yz;
+ localCondMi[t-startTimePoint] = digammaK - inverseKTerm + contributionDigammas + invN_xz + invN_yz;
if (debug) {
// Only tracking this for debugging purposes:
System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, 1/n_yz=%.3f, 1/n_xz=%.3f, local=%.4f\n",
@@ -216,4 +266,13 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2
return returnValues;
}
}
+
+ @Override
+ protected double[] partialComputeFromNewObservations(int startTimePoint,
+ int numTimePoints, double[][] newVar1Observations,
+ double[][] newVar2Observations, double[][] newCondObservations,
+ boolean returnLocals) throws Exception {
+ // TODO implement me
+ throw new RuntimeException("Not implemented yet");
+ }
}
From e2182c40c2c8f71428adac57610aa4bd5c30f3b1 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Tue, 18 Oct 2016 15:04:42 +1100
Subject: [PATCH 10/34] Update to AIS kernel direct class in unit tests to
align with update to AIS interface definition
---
.../kernel/ActiveInfoStorageCalculatorKernelDirect.java | 6 ++++++
1 file changed, 6 insertions(+)
diff --git a/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java
index 09630bc..0fb246f 100755
--- a/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java
+++ b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java
@@ -272,4 +272,10 @@ public class ActiveInfoStorageCalculatorKernelDirect
double[] newObservations) throws Exception {
throw new Exception("Not yet implemented");
}
+
+ public void addObservations(double[] observations, boolean[] valid)
+ throws Exception {
+ // Probably won't need to implement this, as we're only using this class for testing
+ throw new Exception("Not implemented yet");
+ }
}
From 4cac043affd63fb4e3b1831a51c1feadb523ae05 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Tue, 18 Oct 2016 15:05:15 +1100
Subject: [PATCH 11/34] Added unit test for TE continuous
getSeparateNumObservations (in Kraskov class)
---
.../kraskov/TransferEntropyTester.java | 66 +++++++++++++++++++
1 file changed, 66 insertions(+)
diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/TransferEntropyTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/TransferEntropyTester.java
index f227eb8..92bb81c 100755
--- a/java/unittests/infodynamics/measures/continuous/kraskov/TransferEntropyTester.java
+++ b/java/unittests/infodynamics/measures/continuous/kraskov/TransferEntropyTester.java
@@ -537,4 +537,70 @@ public class TransferEntropyTester
assertEquals(teOptimisedSingleThread, teOptimisedWithValidity, 0.00000001);
System.out.println("Answer unchanged by setting validity");
}
+
+ public void testGetSeparateNumObservations() throws Exception {
+ ArrayFileReader afr = new ArrayFileReader("demos/data/SFI-heartRate_breathVol_bloodOx.txt");
+ double[][] data = afr.getDouble2DMatrix();
+
+ TransferEntropyCalculatorKraskov teCalc = new TransferEntropyCalculatorKraskov();
+
+ teCalc.initialise();
+ teCalc.startAddObservations();
+ int timeStepsPerCall = 100;
+ int calls = 10;
+ for (int i = 0; i < calls; i++) {
+ // Add more samples
+ teCalc.addObservations(MatrixUtils.selectColumn(data, 0, i*timeStepsPerCall, timeStepsPerCall),
+ MatrixUtils.selectColumn(data, 1, i*timeStepsPerCall, timeStepsPerCall));
+ }
+ teCalc.finaliseAddObservations();
+ @SuppressWarnings("unused")
+ double result = teCalc.computeAverageLocalOfObservations();
+
+ // Now we want to check how many observations were added at each call:
+ int[] samplesPerCall = teCalc.getSeparateNumObservations();
+ assertEquals(calls, samplesPerCall.length);
+ for (int i = 0; i < calls; i++) {
+ // For k = l = 1, we should have timeStepsPerCall - 1 samples per addObservations() call:
+ assertEquals(timeStepsPerCall - 1, samplesPerCall[i]);
+ }
+
+ // =====================
+ // Now run it again with different k and l and embedding lags, etc:
+ teCalc.initialise();
+ teCalc.startAddObservations();
+ // auto embed destination only
+ teCalc.setProperty(TransferEntropyCalculatorKraskov.PROP_AUTO_EMBED_METHOD,
+ TransferEntropyCalculatorKraskov.AUTO_EMBED_METHOD_RAGWITZ_DEST_ONLY);
+ teCalc.setProperty(TransferEntropyCalculatorKraskov.PROP_K_SEARCH_MAX, "5");
+ teCalc.setProperty(TransferEntropyCalculatorKraskov.PROP_TAU_SEARCH_MAX, "5");
+ // Explicitly set the source embedding params
+ teCalc.setProperty(TransferEntropyCalculatorKraskov.L_PROP_NAME, "1");
+ teCalc.setProperty(TransferEntropyCalculatorKraskov.L_TAU_PROP_NAME, "1");
+ int[] timeStepsPerCallArray = new int[] {100, 200, 150, 300, 99, 54};
+ int startTime = 0;
+ for (int i = 0; i < timeStepsPerCallArray.length; i++) {
+ // Add more samples
+ teCalc.addObservations(MatrixUtils.selectColumn(data, 0, startTime, timeStepsPerCallArray[i]),
+ MatrixUtils.selectColumn(data, 1, startTime, timeStepsPerCallArray[i]));
+ startTime += timeStepsPerCallArray[i];
+ }
+ teCalc.finaliseAddObservations();
+ result = teCalc.computeAverageLocalOfObservations();
+ int optimisedK = Integer.parseInt(teCalc.getProperty(TransferEntropyCalculatorKraskov.K_PROP_NAME));
+ int optimisedKTau = Integer.parseInt(teCalc.getProperty(TransferEntropyCalculatorKraskov.K_TAU_PROP_NAME));
+ int timeOfFirstObservationPerSet = (optimisedK - 1)*optimisedKTau + 1;
+ System.out.printf("In testing tracking of observations per addObservations() call" +
+ " we have auto-embedding dimension %d and lag %d, timeOfFirstObservationPerSet %d\n",
+ optimisedK, optimisedKTau, timeOfFirstObservationPerSet);
+
+ // Now we want to check how many observations were added at each call:
+ samplesPerCall = teCalc.getSeparateNumObservations();
+ assertEquals(timeStepsPerCallArray.length, samplesPerCall.length);
+ for (int i = 0; i < timeStepsPerCallArray.length; i++) {
+ // For timeOfFirstObservationPerSet, we should have timeStepsPerCall - timeOfFirstObservationPerSet
+ // samples per addObservations() call:
+ assertEquals(timeStepsPerCallArray[i] - timeOfFirstObservationPerSet, samplesPerCall[i]);
+ }
+ }
}
From 05beb63e58b0e25eaccb447a35995d7c4e95e5ec Mon Sep 17 00:00:00 2001
From: jlizier
Date: Tue, 18 Oct 2016 15:05:28 +1100
Subject: [PATCH 12/34] Added unit test for AND function
---
.../discrete/MutualInformationTester.java | 15 +++++++++++++++
1 file changed, 15 insertions(+)
diff --git a/java/unittests/infodynamics/measures/discrete/MutualInformationTester.java b/java/unittests/infodynamics/measures/discrete/MutualInformationTester.java
index ca69378..7306cb7 100755
--- a/java/unittests/infodynamics/measures/discrete/MutualInformationTester.java
+++ b/java/unittests/infodynamics/measures/discrete/MutualInformationTester.java
@@ -50,6 +50,21 @@ public class MutualInformationTester extends TestCase {
assertEquals(0.0, miRand, 0.000001);
}
+ public void testAnd() throws Exception {
+ MutualInformationCalculatorDiscrete miCalc = new MutualInformationCalculatorDiscrete(2, 0);
+
+ int[] X1 = new int[] {0, 0, 1, 1};
+ int[] X2 = new int[] {0, 1, 0, 1};
+ int[] Y = new int[] {0, 0, 0, 1};
+
+ // Y is dependent on X1 - MI should be 0.311 bits
+ miCalc.initialise();
+ miCalc.addObservations(X1, Y);
+ double miX1Y = miCalc.computeAverageLocalOfObservations();
+ assertEquals(0.311, miX1Y, 0.001);
+ assertEquals(4, miCalc.getNumObservations());
+ }
+
public void testXor() throws Exception {
MutualInformationCalculatorDiscrete miCalc = new MutualInformationCalculatorDiscrete(2, 0);
From 0386991f003105d74337f0a933080865b0d988a8 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Tue, 18 Oct 2016 15:07:15 +1100
Subject: [PATCH 13/34] Altering use of class methods for python demo 6 (I
think this was required for python 3 compliance?)
---
demos/python/example6DynamicCallingMutualInfo.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
mode change 100755 => 100644 demos/python/example6DynamicCallingMutualInfo.py
diff --git a/demos/python/example6DynamicCallingMutualInfo.py b/demos/python/example6DynamicCallingMutualInfo.py
old mode 100755
new mode 100644
index a0fa1ce..5da8807
--- a/demos/python/example6DynamicCallingMutualInfo.py
+++ b/demos/python/example6DynamicCallingMutualInfo.py
@@ -75,7 +75,7 @@ jointVariable2 = A[:,jointVariable2Columns]
# (in fact, all java object creation in python is dynamic - it has to be,
# since the languages are interpreted. This makes our life slightly easier at this
# point than it is in demos/java/example6 where we have to handle this manually)
-indexOfLastDot = string.rfind(implementingClass, ".")
+indexOfLastDot = implementingClass.rfind(".")
implementingPackage = implementingClass[:indexOfLastDot]
implementingBaseName = implementingClass[indexOfLastDot+1:]
miCalcClass = eval('JPackage(\'%s\').%s' % (implementingPackage, implementingBaseName))
From d74749f17dbb208ad3fc7910efd15804b4a87656 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Tue, 18 Oct 2016 15:08:15 +1100
Subject: [PATCH 14/34] Showing how to use numpy arrays of ints with JIDT in
python 3 with jpype1 in python demo 1
---
demos/python/example1TeBinaryData.py | 25 ++++++++++++++++++++-----
1 file changed, 20 insertions(+), 5 deletions(-)
diff --git a/demos/python/example1TeBinaryData.py b/demos/python/example1TeBinaryData.py
index 8e4cded..92dfa47 100755
--- a/demos/python/example1TeBinaryData.py
+++ b/demos/python/example1TeBinaryData.py
@@ -20,13 +20,14 @@
# Simple transfer entropy (TE) calculation on binary data using the discrete TE calculator:
-from jpype import *
+import jpype
import random
+import numpy
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
-startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
+jpype.startJVM(jpype.getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
# Generate some random binary data.
sourceArray = [random.randint(0,1) for r in range(100)]
@@ -34,15 +35,29 @@ destArray = [0] + sourceArray[0:99]
sourceArray2 = [random.randint(0,1) for r in range(100)]
# Create a TE calculator and run it:
-teCalcClass = JPackage("infodynamics.measures.discrete").TransferEntropyCalculatorDiscrete
+teCalcClass = jpype.JPackage("infodynamics.measures.discrete").TransferEntropyCalculatorDiscrete
teCalc = teCalcClass(2,1)
teCalc.initialise()
-# Since we have simple arrays of ints, we can directly pass these in:
+
+# First use simple arrays of ints, which we can directly pass in:
teCalc.addObservations(sourceArray, destArray)
print("For copied source, result should be close to 1 bit : %.4f" % teCalc.computeAverageLocalOfObservations())
teCalc.initialise()
teCalc.addObservations(sourceArray2, destArray)
print("For random source, result should be close to 0 bits: %.4f" % teCalc.computeAverageLocalOfObservations())
-shutdownJVM()
+# Next, demonstrate how to do this with a numpy array
+teCalc.initialise()
+# Create the numpy arrays:
+sourceNumpy = numpy.array(sourceArray, dtype=numpy.int)
+destNumpy = numpy.array(destArray, dtype=numpy.int)
+# The above can be passed straight through to JIDT in python 2:
+# teCalc.addObservations(sourceNumpy, destNumpy)
+# But you need to do this in python 3:
+sourceNumpyJArray = jpype.JArray(jpype.JInt, 1)(sourceNumpy.tolist())
+destNumpyJArray = jpype.JArray(jpype.JInt, 1)(destNumpy.tolist())
+teCalc.addObservations(sourceNumpyJArray, destNumpyJArray)
+print("Using numpy array for copied source, result confirmed as: %.4f" % teCalc.computeAverageLocalOfObservations())
+
+jpype.shutdownJVM()
From 2d4e30aeb47e18e74de51dd632d01731315e211c Mon Sep 17 00:00:00 2001
From: jlizier
Date: Tue, 18 Oct 2016 15:09:59 +1100
Subject: [PATCH 15/34] Minor bug patch for readFloatsFile for python demos
---
demos/python/readFloatsFile.py | 9 +++++----
1 file changed, 5 insertions(+), 4 deletions(-)
mode change 100755 => 100644 demos/python/readFloatsFile.py
diff --git a/demos/python/readFloatsFile.py b/demos/python/readFloatsFile.py
old mode 100755
new mode 100644
index ff37da1..2d0e6b7
--- a/demos/python/readFloatsFile.py
+++ b/demos/python/readFloatsFile.py
@@ -19,15 +19,16 @@
def readFloatsFile(filename):
"Read a 2D array of floats from a given file"
with open(filename) as f:
- # Space separate numbers, one time step per line, each column is a variable
- array = []
- for line in f: # read all lines
+ # Space separate numbers, one time step per line, each column is a variable
+ array = []
+ for line in f:
+ # read all lines
if (line.startswith("%") or line.startswith("#")):
# Assume this is a comment line
continue
if (len(line.split()) == 0):
# Line is empty
continue
- array.append([float(x) for x in line.split()])
+ array.append([float(x) for x in line.split()])
return array
From 9a7ff306c834d6eb8ec7f01904827edf38b2e75c Mon Sep 17 00:00:00 2001
From: jlizier
Date: Thu, 20 Oct 2016 09:26:12 +1100
Subject: [PATCH 16/34] Adding utilities for fast nearest neighbour searchers
to handle calculations on new samples and prepare for conditional MI fast
surrogates computation
---
java/source/infodynamics/utils/KdTree.java | 1321 ++++++++++++++++-
.../utils/NearestNeighbourSearcher.java | 93 +-
.../UnivariateNearestNeighbourSearcher.java | 480 ++++++
.../infodynamics/utils/KdTreeTest.java | 288 ++++
.../utils/UnivariateNearestNeighbourTest.java | 247 +++
5 files changed, 2422 insertions(+), 7 deletions(-)
diff --git a/java/source/infodynamics/utils/KdTree.java b/java/source/infodynamics/utils/KdTree.java
index ba2fa53..b049a2d 100755
--- a/java/source/infodynamics/utils/KdTree.java
+++ b/java/source/infodynamics/utils/KdTree.java
@@ -206,10 +206,6 @@ public class KdTree extends NearestNeighbourSearcher {
// Precondition: sortedArrayIndices[][] are currently sorted for all
// dimensions
- // Point to the correct array for the data
- double[][] data = dimensionToArray[currentDim];
- int actualDim = dimensionToArrayIndex[currentDim];
-
// Handle non-recursive solutions:
if (numPoints == 0) {
return null;
@@ -225,6 +221,10 @@ public class KdTree extends NearestNeighbourSearcher {
new KdTreeNode(sortedArrayIndices[currentDim][startPoint+1], null, null));
}
+ // Point to the correct array for the data
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+
// Identify the point on which to split here
int candidateSplitPoint = startPoint + numPoints/2;
while ((candidateSplitPoint > startPoint) &&
@@ -776,6 +776,137 @@ public class KdTree extends NearestNeighbourSearcher {
}
}
+ public PriorityQueue
+ findKNearestNeighbours(int K, double[][] sampleVectors) throws Exception {
+
+ if (originalDataSets[0].length <= K) {
+ throw new Exception("Not enough data points for a K nearest neighbours search");
+ }
+
+ PriorityQueue pq = new PriorityQueue(K);
+ if (rootNode == null) {
+ return pq;
+ }
+ findKNearestNeighbours(K, sampleVectors, rootNode, 0, pq);
+ return pq;
+ }
+
+ /**
+ * Protected method to Update the k nearest neighbours to a given sample (samplePoint), from the tree
+ * rooted at node (which is at the specified level in the tree).
+ * Incorporate neighbours found in this sub-tree into the PriorityQueue
+ * of the current K closest.
+ * @param K the number of nearest neighbour
+ * @param sampleVectors sample vectors (not necessarily in the tree) 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 currentKBest a PriorityQueue of NeighbourNodeData objects
+ * capturing the current K closest neighbours and their distances.
+ * Assumed not to be null, but may be empty or with less than K elements
+ * so far. It must be sorted from furthest away first to nearest last.
+ */
+ protected void findKNearestNeighbours(int K,
+ double[][] sampleVectors,
+ KdTreeNode node, int level,
+ PriorityQueue currentKBest) {
+
+ // Point to the correct array for the data at this level
+ int currentDim = level % totalDimensions;
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+ int variableNumber = dimensionToVariableNumber[currentDim];
+
+ // Check the distance on this particular dimension
+ double distOnThisDim = sampleVectors[variableNumber][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 (this saves taking square roots anywhere)
+ absDistOnThisDim = distOnThisDim * distOnThisDim;
+ }
+
+ // Grab the current furthest nearest neighbour in our cached list
+ // (will not throw an Exception if the PQ is empty)
+ NeighbourNodeData furthestCached = currentKBest.peek();
+
+ if ((currentKBest.size() < K) || (absDistOnThisDim < furthestCached.distance)) {
+ // Preliminary check says we need to compute the full distance
+ // to use or at least to check if it should be
+ // added to our currentKBest properly.
+ double maxNorm = 0;
+ double[] norms = new double[originalDataSets.length];
+ for (int v = 0; v < originalDataSets.length; v++) {
+ // For each of our separate (multivariate) variables,
+ // compute the (specified) norm in that variable's space:
+ if (currentKBest.size() < K) {
+ norms[v] = norm(
+ sampleVectors[v],
+ originalDataSets[v][node.indexOfThisPoint], normTypeToUse);
+ } else {
+ // Distance calculation terminates early with Double.POSITIVE_INFINITY
+ // if it is clearly larger than currentBest.distance:
+ norms[v] = normWithAbort(
+ sampleVectors[v],
+ originalDataSets[v][node.indexOfThisPoint],
+ furthestCached.distance, normTypeToUse);
+ }
+ if (norms[v] > maxNorm) {
+ maxNorm = norms[v];
+ if (Double.isInfinite(maxNorm)) {
+ // we've aborted the norm check early;
+ // no point checking the other variables.
+ break;
+ }
+ }
+ }
+ if ((currentKBest.size() < K) ||
+ (maxNorm < furthestCached.distance)) {
+ // We add this to our cache of K nearest neighbours:
+ if (currentKBest.size() == K) {
+ // Remove the current Kth nearest neighbour
+ // as it is about to be replaced.
+ currentKBest.poll();
+ }
+ currentKBest.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) {
+ findKNearestNeighbours(K, sampleVectors,
+ closestSubTree, level + 1, currentKBest);
+ }
+ // Grab the current furthest nearest neighbour in our cached list again
+ // (will not throw an Exception if the PQ is empty)
+ // as it may have been changed above:
+ furthestCached = currentKBest.peek();
+ if ((currentKBest.size() < K) || (absDistOnThisDim < furthestCached.distance)) {
+ // It's possible we could have a closer node than the current best
+ // in the other branch as well, so search there too:
+ if (furthestSubTree != null) {
+ findKNearestNeighbours(K, sampleVectors,
+ furthestSubTree, level + 1, currentKBest);
+ }
+ }
+ }
+
@Override
public int countPointsWithinR(int sampleIndex, double r, boolean allowEqualToR) {
if (rootNode == null) {
@@ -1009,6 +1140,122 @@ public class KdTree extends NearestNeighbourSearcher {
return count;
}
+ @Override
+ public int countPointsWithinR(double[][] sampleVectors, double r, boolean allowEqualToR) {
+ if (rootNode == null) {
+ return 0;
+ }
+ return countPointsWithinR(sampleVectors, rootNode, 0, r, allowEqualToR);
+ }
+
+ /**
+ * Count the number of points within radius r of a given sample (sampleVector),
+ * 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 sampleVectors sample vectors to find the neighbours within r
+ * 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
+ * @return count of points within r
+ */
+ protected int countPointsWithinR(double[][] sampleVectors,
+ KdTreeNode node, int level, double r, boolean allowEqualToR) {
+
+ int count = 0;
+
+ // Point to the correct array for the data at this level
+ int currentDim = level % totalDimensions;
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+ int variableNumber = dimensionToVariableNumber[currentDim];
+
+ // Check the distance on this particular dimension
+ double distOnThisDim = sampleVectors[variableNumber][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 (((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;
+ 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(
+ sampleVectors[v],
+ originalDataSets[v][node.indexOfThisPoint],
+ r, normTypeToUse);
+ if ((distForVariableV >= r) &&
+ !(allowEqualToR && (distForVariableV == r))) {
+ // We don't fit on this dimension, no point
+ // checking the others:
+ withinBounds = false;
+ break;
+ }
+ }
+ if (withinBounds) {
+ // This node gets counted
+ count++;
+ }
+ }
+
+ 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) {
+ count += countPointsWithinR(sampleVectors, closestSubTree,
+ level + 1, r, allowEqualToR);
+ }
+ 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) {
+ count += countPointsWithinR(sampleVectors, furthestSubTree,
+ level + 1, r, allowEqualToR);
+ }
+ }
+
+ return count;
+ }
+
@Override
public Collection findPointsWithinR(int sampleIndex,
double r, boolean allowEqualToR) {
@@ -1300,7 +1547,7 @@ public class KdTree extends NearestNeighbourSearcher {
/**
* As per {@link #findPointsWithinR(int, KdTreeNode, int, double, boolean, boolean[], int[], int)}
- * however with dynamic correaltion time specified
+ * however with dynamic correlation time specified
*
* @param sampleIndex sample index in the data to find a nearest neighbour
* for
@@ -1420,6 +1667,487 @@ public class KdTree extends NearestNeighbourSearcher {
return nextIndexInIndicesWithinR;
}
+ @Override
+ public int findPointsWithinR(int sampleIndex, double r,
+ int dynCorrExclTime, boolean allowEqualToR,
+ boolean[] isWithinR,
+ double[] distancesWithinRInOrder,
+ double[][] distancesAndIndicesWithinR) {
+ if (rootNode == null) {
+ distancesAndIndicesWithinR[0][1] = -1;
+ return 0;
+ }
+ int currentIndexInIndicesWithinR =
+ findPointsWithinR(sampleIndex,
+ rootNode, 0, r, dynCorrExclTime, allowEqualToR, isWithinR,
+ distancesWithinRInOrder, distancesAndIndicesWithinR, 0);
+ distancesAndIndicesWithinR[currentIndexInIndicesWithinR][1] = -1;
+ return currentIndexInIndicesWithinR;
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, KdTreeNode, int, double, boolean, boolean[], int[], int)}
+ * however with dynamic correlation time specified, and
+ * returning distances to neighbours as well as neighbour indices.
+ *
+ * @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 dynCorrExclTime time window within which to exclude
+ * points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param isWithinR the array should be passed in with all points set to
+ * false initially, and is return indicating whether each sample was
+ * found to be within r of that at sampleIndex.
+ * @param distancesWithinRInOrder the array must be passed in and is
+ * returned with distances for each point found to be within r. Values
+ * at other indices are not defined.
+ * @param distancesAndIndicesWithinR a list of distances (in column index 0)
+ * and array indices (in column index 1)
+ * for points marked as true in isWithinR, terminated with a -1 value on the index.
+ * @param nextIndexInIndicesWithinR the next available index in indicesWithinR
+ * before the method.
+ * @return the next available index in distancesAndIndicesWithinR after the method is complete
+ */
+ protected int findPointsWithinR(int sampleIndex,
+ KdTreeNode node, int level, double r,
+ int dynCorrExclTime, boolean allowEqualToR,
+ boolean[] isWithinR,
+ double[] distancesWithinRInOrder,
+ double[][] distancesAndIndicesWithinR,
+ int nextIndexInIndicesWithinR) {
+
+ // 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 ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+ ((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
+ isWithinR[node.indexOfThisPoint] = true;
+ distancesWithinRInOrder[node.indexOfThisPoint] = maxNorm;
+ distancesAndIndicesWithinR[nextIndexInIndicesWithinR][0] =
+ maxNorm;
+ distancesAndIndicesWithinR[nextIndexInIndicesWithinR++][1] =
+ node.indexOfThisPoint;
+ }
+ }
+
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinR(sampleIndex, closestSubTree,
+ level + 1, r, dynCorrExclTime, allowEqualToR, isWithinR,
+ distancesWithinRInOrder,
+ distancesAndIndicesWithinR, nextIndexInIndicesWithinR);
+ }
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinR(sampleIndex, furthestSubTree,
+ level + 1, r, dynCorrExclTime, allowEqualToR, isWithinR,
+ distancesWithinRInOrder,
+ distancesAndIndicesWithinR, nextIndexInIndicesWithinR);
+ }
+ }
+ return nextIndexInIndicesWithinR;
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, double, int, boolean, boolean[], double[], double[][])}
+ * only we have a sub-variable already tested, and we incorporate its results here.
+ *
+ * @param sampleIndex sample index in the data to find a nearest neighbour
+ * for
+ * @param r radius within which to count points
+ * @param dynCorrExclTime time window within which to exclude
+ * points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param variableAlreadyTested
+ * @param testResultsForGivenVariable an array indicating whether each sample for
+ * variableAlreadyTested was within radius for this search point
+ * @param distancesWithinRForGivenVariable an array of distances for
+ * variableAlreadyTested for each sample
+ * found to be within r. Values at other indices are not defined.
+ * @param distancesAndIndicesWithinR a list of distances (in column index 0)
+ * and array indices (in column index 1)
+ * for points found to be within r, terminated with a -1 value on the index.
+ * @return the next available index in distancesAndIndicesWithinR after the method is complete
+ */
+ public int findPointsWithinRAndTakeRadiusMax(int sampleIndex, double r,
+ int dynCorrExclTime, boolean allowEqualToR,
+ int variableAlreadyTested, boolean[] testResultsForGivenVariable,
+ double[] distancesWithinRForGivenVariable,
+ double[][] distancesAndIndicesWithinR) {
+ if (rootNode == null) {
+ distancesAndIndicesWithinR[0][1] = -1;
+ return 0;
+ }
+ int currentIndexInIndicesWithinR =
+ findPointsWithinRAndTakeRadiusMax(sampleIndex,
+ rootNode, 0, r, dynCorrExclTime, allowEqualToR,
+ variableAlreadyTested, testResultsForGivenVariable,
+ distancesWithinRForGivenVariable, distancesAndIndicesWithinR, 0);
+ distancesAndIndicesWithinR[currentIndexInIndicesWithinR][1] = -1;
+ return currentIndexInIndicesWithinR;
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, KdTreeNode, int, double, int, boolean, boolean[], double[], double[][], int)}
+ * only we have a sub-variable already tested, and we incorporate its results here.
+ *
+ * @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 dynCorrExclTime time window within which to exclude
+ * points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param variableAlreadyTested
+ * @param testResultsForGivenVariable an array indicating whether each sample for
+ * variableAlreadyTested was within radius for this search point
+ * @param distancesWithinRForGivenVariable an array of distances for
+ * variableAlreadyTested for each sample
+ * found to be within r. Values at other indices are not defined.
+ * @param distancesAndIndicesWithinR a list of distances (in column index 0)
+ * and array indices (in column index 1)
+ * for points found to be within r, terminated with a -1 value on the index.
+ * @param nextIndexInIndicesWithinR the next available index in indicesWithinR
+ * before the method.
+ * @return the next available index in distancesAndIndicesWithinR after the method is complete
+ */
+ protected int findPointsWithinRAndTakeRadiusMax(int sampleIndex,
+ KdTreeNode node, int level, double r,
+ int dynCorrExclTime, boolean allowEqualToR,
+ int variableAlreadyTested, boolean[] testResultsForGivenVariable,
+ double[] distancesWithinRForGivenVariable,
+ double[][] distancesAndIndicesWithinR,
+ int nextIndexInIndicesWithinR) {
+
+ // 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 (testResultsForGivenVariable[node.indexOfThisPoint] &&
+ (Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+ ((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++) {
+ if (v == variableAlreadyTested) {
+ // This variable has already been checked to be ok above
+ continue;
+ }
+ // 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
+ if (maxNorm < distancesWithinRForGivenVariable[node.indexOfThisPoint]) {
+ // The max norm actually comes from the pre-tested variable:
+ maxNorm = distancesWithinRForGivenVariable[node.indexOfThisPoint];
+ }
+ distancesAndIndicesWithinR[nextIndexInIndicesWithinR][0] =
+ maxNorm;
+ distancesAndIndicesWithinR[nextIndexInIndicesWithinR++][1] =
+ node.indexOfThisPoint;
+ }
+ }
+
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinRAndTakeRadiusMax(sampleIndex, closestSubTree,
+ level + 1, r, dynCorrExclTime, allowEqualToR,
+ variableAlreadyTested, testResultsForGivenVariable,
+ distancesWithinRForGivenVariable,
+ distancesAndIndicesWithinR, nextIndexInIndicesWithinR);
+ }
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinRAndTakeRadiusMax(sampleIndex, furthestSubTree,
+ level + 1, r, dynCorrExclTime, allowEqualToR,
+ variableAlreadyTested, testResultsForGivenVariable,
+ distancesWithinRForGivenVariable,
+ distancesAndIndicesWithinR, nextIndexInIndicesWithinR);
+ }
+ }
+ return nextIndexInIndicesWithinR;
+ }
+
+ @Override
+ public void findPointsWithinR(
+ double r, double[][] samplePoint,
+ boolean allowEqualToR, boolean[] isWithinR,
+ int[] indicesWithinR) {
+ if (rootNode == null) {
+ indicesWithinR[0] = -1;
+ return;
+ }
+ int currentIndexInIndicesWithinR =
+ findPointsWithinR(samplePoint,
+ rootNode, 0, r, allowEqualToR, isWithinR,
+ indicesWithinR, 0);
+ indicesWithinR[currentIndexInIndicesWithinR] = -1;
+ return;
+ }
+
+ /**
+ * Record the collection of points within radius r of a given sample (sampleVectors),
+ * 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).
+ *
+ * The recording of nearest neighbours is made within the isWithinR
+ * and indicesWithinR arrays, which must be constructed before
+ * calling this method, with length at or exceeding the total
+ * number of data points. indicesWithinR is
+ *
+ *
+ * @param sampleVectors sample vectors to find the neighbours within r
+ * 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 isWithinR the array should be passed in with all points set to
+ * false initially, and is return indicating whether each sample was
+ * found to be within r of that at sampleIndex.
+ * @param indicesWithinR a list of array indices
+ * for points marked as true in isWithinR, terminated with a -1 value.
+ * @param nextIndexInIndicesWithinR the next available index in indicesWithinR
+ * before the method.
+ * @return the next available index in indicesWithinR after the method is complete
+ */
+ protected int findPointsWithinR(double[][] sampleVectors,
+ KdTreeNode node, int level, double r, boolean allowEqualToR,
+ boolean[] isWithinR, int[] indicesWithinR,
+ int nextIndexInIndicesWithinR) {
+
+ // Point to the correct array for the data at this level
+ int currentDim = level % totalDimensions;
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+ int variableNumber = dimensionToVariableNumber[currentDim];
+
+ // Check the distance on this particular dimension
+ double distOnThisDim = sampleVectors[variableNumber][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 ((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(
+ sampleVectors[v],
+ 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
+ isWithinR[node.indexOfThisPoint] = true;
+ indicesWithinR[nextIndexInIndicesWithinR++] =
+ node.indexOfThisPoint;
+ }
+ }
+
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinR(sampleVectors, closestSubTree,
+ level + 1, r, allowEqualToR, isWithinR, indicesWithinR,
+ nextIndexInIndicesWithinR);
+ }
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinR(sampleVectors, furthestSubTree,
+ level + 1, r, allowEqualToR, isWithinR, indicesWithinR,
+ nextIndexInIndicesWithinR);
+ }
+ }
+ return nextIndexInIndicesWithinR;
+ }
+
/**
* As per {@link #countPointsWithinR(int, double, boolean)}
* with an additional factor as follows.
@@ -1616,6 +2344,466 @@ public class KdTree extends NearestNeighbourSearcher {
return count;
}
+ /**
+ * As per {@link #countPointsWithinR(int, double, boolean)}, except searches for a specified
+ * sample point rather than a point associated with an index in the existing sample set, and
+ * with an additional factor as follows.
+ *
+ * For this method, one of the high-level variables has already had their
+ * norms tested against r; this variable is indicated by variableAlreadyTested
+ * and the results are contained in the testResultsForGivenVariable array.
+ *
+ *
+ * @param r radius within which to count points
+ * @param sampleVectors sample vectors (not necessarily in the tree) to find a nearest neighbour
+ * for
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param variableAlreadyTested the variable for which the test results
+ * are provided in testResultsForGivenVariable
+ * @param testResultsForGivenVariable array of booleans indicating whether
+ * the variable identified by variableAlreadyTested was within r for each
+ * sample index.
+ * @return the count of points within r.
+ */
+ public int countPointsWithinR(double r, double[][] sampleVectors, boolean allowEqualToR,
+ int variableAlreadyTested, boolean[] testResultsForGivenVariable) {
+ if (rootNode == null) {
+ return 0;
+ }
+ return countPointsWithinR(sampleVectors, rootNode, 0, r, allowEqualToR,
+ variableAlreadyTested, testResultsForGivenVariable);
+ }
+
+ /**
+ * Count the number of points within radius r of a given sample (sampleVectors),
+ * in the tree rooted at node (which is at the specified level in the tree),
+ * when one of the variables has already been tested.
+ * 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 sampleVectors sample vectors (not necessarily in the tree) 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 variableAlreadyTested the variable for which the test results
+ * are provided in testResultsForGivenVariable
+ * @param testResultsForGivenVariable array of booleans indicating whether
+ * the variable identified by variableAlreadyTested was within r for each
+ * sample index.
+ * @return count of points within r
+ */
+ protected int countPointsWithinR(double[][] sampleVectors,
+ KdTreeNode node, int level, double r, boolean allowEqualToR,
+ int variableAlreadyTested, boolean[] testResultsForGivenVariable) {
+
+ int count = 0;
+
+ // Point to the correct array for the data at this level
+ int currentDim = level % totalDimensions;
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+ int variableNumber = dimensionToVariableNumber[currentDim];
+
+ // Check the distance on this particular dimension
+ double distOnThisDim = sampleVectors[variableNumber][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 (testResultsForGivenVariable[node.indexOfThisPoint] &&
+ ((variableNumber == variableAlreadyTested) ||
+ (absDistOnThisDim < r) ||
+ (allowEqualToR && (absDistOnThisDim == r)))) {
+ // Preliminary check (including that the dimension already tested was ok)
+ // says we need to compute the full distance
+ // to use or at least to check if it should be counted.
+
+ boolean withinBounds = true;
+ for (int v = 0; v < originalDataSets.length; v++) {
+ if (v == variableAlreadyTested) {
+ // This variable has already been checked to be ok above
+ continue;
+ }
+ // 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(
+ sampleVectors[v],
+ originalDataSets[v][node.indexOfThisPoint],
+ r, normTypeToUse);
+ if ((distForVariableV >= r) &&
+ !(allowEqualToR && (distForVariableV == r))) {
+ // We don't fit on this dimension, no point
+ // checking the others:
+ withinBounds = false;
+ break;
+ }
+ }
+ if (withinBounds) {
+ // This node gets counted
+ count++;
+ }
+ }
+
+ 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) {
+ count += countPointsWithinR(sampleVectors, closestSubTree,
+ level + 1, r, allowEqualToR,
+ variableAlreadyTested, testResultsForGivenVariable);
+ }
+ 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) {
+ count += countPointsWithinR(sampleVectors, furthestSubTree,
+ level + 1, r, allowEqualToR,
+ variableAlreadyTested, testResultsForGivenVariable);
+ }
+ }
+
+ return count;
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, double, int, boolean, boolean[], int[])}
+ * however using multiple radii for each variable.
+ *
+ * @param sampleIndex sample index in the data to find a nearest neighbour for
+ * @param rs radii within which to count points (one for each variable)
+ * @param dynCorrExclTime time window within which to exclude points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also, otherwise only those strictly within r
+ * @param isWithinR the array MUST be passed in with all points set to false initially, and is returned indicating whether each sample was found to be within r of that at sampleIndex.
+ * @param indicesWithinR a list of array indices for points marked as true in isWithinR, terminated with a -1 value.
+ */
+ public void findPointsWithinRs(int sampleIndex,
+ double[] rs, int dynCorrExclTime, boolean allowEqualToR,
+ boolean[] isWithinR, int[] indicesWithinR) {
+ if (rootNode == null) {
+ indicesWithinR[0] = -1;
+ return;
+ }
+ int currentIndexInIndicesWithinR =
+ findPointsWithinRs(sampleIndex,
+ rootNode, 0, rs, dynCorrExclTime, allowEqualToR, isWithinR,
+ indicesWithinR, 0);
+ indicesWithinR[currentIndexInIndicesWithinR] = -1;
+ return;
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, KdTreeNode, int, double, int, boolean, boolean[], int[], int)}
+ * however with dynamic correlation time specified
+ *
+ * @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 rs radii within which to count points (one for each variable)
+ * @param dynCorrExclTime time window within which to exclude
+ * points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param isWithinR the array should be passed in with all points set to
+ * false initially, and is return indicating whether each sample was
+ * found to be within r of that at sampleIndex.
+ * @param indicesWithinR a list of array indices
+ * for points marked as true in isWithinR, terminated with a -1 value.
+ * @param nextIndexInIndicesWithinR the next available index in indicesWithinR
+ * before the method.
+ * @return the next available index in indicesWithinR after the method is complete
+ */
+ protected int findPointsWithinRs(int sampleIndex,
+ KdTreeNode node, int level, double[] rs,
+ int dynCorrExclTime, boolean allowEqualToR,
+ boolean[] isWithinR, int[] indicesWithinR,
+ int nextIndexInIndicesWithinR) {
+
+ // Point to the correct array for the data at this level
+ int currentDim = level % totalDimensions;
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+ int variableNumber = dimensionToVariableNumber[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 ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+ ((absDistOnThisDim < rs[variableNumber]) ||
+ ( allowEqualToR && (absDistOnThisDim == rs[variableNumber])))) {
+ // 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],
+ rs[v], normTypeToUse);
+ if ((distForVariableV >= rs[v]) &&
+ !(allowEqualToR && (distForVariableV == rs[v]))) {
+ // 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
+ isWithinR[node.indexOfThisPoint] = true;
+ indicesWithinR[nextIndexInIndicesWithinR++] =
+ node.indexOfThisPoint;
+ }
+ }
+
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinRs(sampleIndex, closestSubTree,
+ level + 1, rs, dynCorrExclTime, allowEqualToR, isWithinR,
+ indicesWithinR, nextIndexInIndicesWithinR);
+ }
+ if ((absDistOnThisDim < rs[variableNumber]) ||
+ ( allowEqualToR && (distOnThisDim < 0) && (absDistOnThisDim == rs[variableNumber]))) {
+ // 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) {
+ nextIndexInIndicesWithinR = findPointsWithinRs(sampleIndex, furthestSubTree,
+ level + 1, rs, dynCorrExclTime, allowEqualToR, isWithinR,
+ indicesWithinR, nextIndexInIndicesWithinR);
+ }
+ }
+ return nextIndexInIndicesWithinR;
+ }
+ /**
+ * As per {@link #findPointsWithinR(double[], double[][], boolean, boolean[], int[]))}
+ * however using multiple radii for each variable.
+ *
+ * @param rs radii within which to count points (one for each variable)
+ * @param samplePoint sample vectors to find the neighbours within rs for
+ * @param allowEqualToR if true, then count points at radius r also, otherwise only those strictly within r
+ * @param isWithinR the array MUST be passed in with all points set to false initially, and is returned indicating whether each sample was found to be within r of that at sampleIndex.
+ * @param indicesWithinR a list of array indices for points marked as true in isWithinR, terminated with a -1 value.
+ */
+ public void findPointsWithinRs(
+ double[] rs, double[][] samplePoint,
+ boolean allowEqualToR, boolean[] isWithinR,
+ int[] indicesWithinR) {
+ if (rootNode == null) {
+ indicesWithinR[0] = -1;
+ return;
+ }
+ int currentIndexInIndicesWithinR =
+ findPointsWithinRs(samplePoint,
+ rootNode, 0, rs, allowEqualToR, isWithinR,
+ indicesWithinR, 0);
+ indicesWithinR[currentIndexInIndicesWithinR] = -1;
+ return;
+ }
+
+ /**
+ * Record the collection of points within radii rs of a given sample (sampleVectors),
+ * in the tree rooted at node (which is at the specified level in the tree).
+ * Nearest neighbour function to compare to rs 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).
+ *
+ * The recording of nearest neighbours is made within the isWithinR
+ * and indicesWithinR arrays, which must be constructed before
+ * calling this method, with length at or exceeding the total
+ * number of data points. indicesWithinR is
+ *
+ *
+ * @param sampleVectors sample vectors to find the neighbours within r
+ * 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 rs radii within which to count points (one for each variable)
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param isWithinR the array should be passed in with all points set to
+ * false initially, and is return indicating whether each sample was
+ * found to be within r of that at sampleIndex.
+ * @param indicesWithinR a list of array indices
+ * for points marked as true in isWithinR, terminated with a -1 value.
+ * @param nextIndexInIndicesWithinR the next available index in indicesWithinR
+ * before the method.
+ * @return the next available index in indicesWithinR after the method is complete
+ */
+ protected int findPointsWithinRs(double[][] sampleVectors,
+ KdTreeNode node, int level, double[] rs, boolean allowEqualToR,
+ boolean[] isWithinR, int[] indicesWithinR,
+ int nextIndexInIndicesWithinR) {
+
+ // Point to the correct array for the data at this level
+ int currentDim = level % totalDimensions;
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+ int variableNumber = dimensionToVariableNumber[currentDim];
+
+ // Check the distance on this particular dimension
+ double distOnThisDim = sampleVectors[variableNumber][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 ((absDistOnThisDim < rs[variableNumber]) ||
+ ( allowEqualToR && (absDistOnThisDim == rs[variableNumber]))) {
+ // 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(
+ sampleVectors[v],
+ originalDataSets[v][node.indexOfThisPoint],
+ rs[v], normTypeToUse);
+ if ((distForVariableV >= rs[v]) &&
+ !(allowEqualToR && (distForVariableV == rs[v]))) {
+ // 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
+ isWithinR[node.indexOfThisPoint] = true;
+ indicesWithinR[nextIndexInIndicesWithinR++] =
+ node.indexOfThisPoint;
+ }
+ }
+
+ 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) {
+ nextIndexInIndicesWithinR = findPointsWithinRs(sampleVectors, closestSubTree,
+ level + 1, rs, allowEqualToR, isWithinR, indicesWithinR,
+ nextIndexInIndicesWithinR);
+ }
+ if ((absDistOnThisDim < rs[variableNumber]) ||
+ ( allowEqualToR && (distOnThisDim < 0) && (absDistOnThisDim == rs[variableNumber]))) {
+ // 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) {
+ nextIndexInIndicesWithinR = findPointsWithinRs(sampleVectors, furthestSubTree,
+ level + 1, rs, allowEqualToR, isWithinR, indicesWithinR,
+ nextIndexInIndicesWithinR);
+ }
+ }
+ return nextIndexInIndicesWithinR;
+ }
+
/**
* As per {@link #countPointsWithinRs(int, double[], boolean)}
* with an additional factor as follows.
@@ -1940,6 +3128,129 @@ public class KdTree extends NearestNeighbourSearcher {
return count;
}
+ @Override
+ public int countPointsWithinR(double[][] sampleVectors, double r,
+ boolean allowEqualToR, boolean[] additionalCriteria) {
+ if (rootNode == null) {
+ return 0;
+ }
+ return countPointsWithinR(sampleVectors, rootNode, 0, r, allowEqualToR,
+ additionalCriteria);
+ }
+
+ /**
+ * Count the number of points within radius r of a given sample (sampleVectors -- from outside the tree),
+ * in the tree rooted at node (which is at the specified level in the tree),
+ * subject to an additional criteria.
+ * 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 sampleVectors sample vectors to find the neighbours within r
+ * 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 additionalCriteria array of booleans. Only count a point if it
+ * is within r and is true in additionalCrtieria.
+ * @return count of points within r
+ */
+ protected int countPointsWithinR(double[][] sampleVectors,
+ KdTreeNode node, int level, double r, boolean allowEqualToR,
+ boolean[] additionalCriteria) {
+
+ int count = 0;
+
+ // Point to the correct array for the data at this level
+ int currentDim = level % totalDimensions;
+ double[][] data = dimensionToArray[currentDim];
+ int actualDim = dimensionToArrayIndex[currentDim];
+ int variableNumber = dimensionToVariableNumber[currentDim];
+
+ // Check the distance on this particular dimension
+ double distOnThisDim = sampleVectors[variableNumber][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 (additionalCriteria[node.indexOfThisPoint] &&
+ ((absDistOnThisDim < r) ||
+ ( allowEqualToR && (absDistOnThisDim == r)))) {
+ // Preliminary check (including additional criteria) says we need to compute the full distance
+ // to use or at least to check if it should be counted.
+ boolean withinBounds = true;
+ 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(
+ sampleVectors[v],
+ originalDataSets[v][node.indexOfThisPoint],
+ r, normTypeToUse);
+ if ((distForVariableV >= r) &&
+ !(allowEqualToR && (distForVariableV == r))) {
+ // We don't fit on this dimension, no point
+ // checking the others:
+ withinBounds = false;
+ break;
+ }
+ }
+ if (withinBounds) {
+ // This node gets counted
+ count++;
+ }
+ }
+
+ 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) {
+ count += countPointsWithinR(sampleVectors, closestSubTree,
+ level + 1, r, allowEqualToR, additionalCriteria);
+ }
+ 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) {
+ count += countPointsWithinR(sampleVectors, furthestSubTree,
+ level + 1, r, allowEqualToR, additionalCriteria);
+ }
+ }
+
+ return count;
+ }
+
/**
* Count the number of points within norms {r1,r2,etc} for each high-level
* variable, for a given
diff --git a/java/source/infodynamics/utils/NearestNeighbourSearcher.java b/java/source/infodynamics/utils/NearestNeighbourSearcher.java
index b70a561..9d9f216 100755
--- a/java/source/infodynamics/utils/NearestNeighbourSearcher.java
+++ b/java/source/infodynamics/utils/NearestNeighbourSearcher.java
@@ -49,7 +49,10 @@ public abstract class NearestNeighbourSearcher {
public static NearestNeighbourSearcher create(double[][] data)
throws Exception {
- if (data[0].length == 1) {
+ if ((data == null) || (data[0].length == 0)) {
+ // We have null data:
+ return null;
+ } else if (data[0].length == 1) {
// We have univariate data:
return new UnivariateNearestNeighbourSearcher(MatrixUtils.selectColumn(data, 0));
} else {
@@ -271,6 +274,23 @@ public abstract class NearestNeighbourSearcher {
public abstract int countPointsWithinR(int sampleIndex, double r,
int dynCorrExclTime, boolean allowEqualToR);
+ /**
+ * As per {@link #countPointsWithinR(int, double, boolean)}
+ * however the search is to match a specified sample point (not a sample
+ * point within the search space itself).
+ *
+ * @param sampleVectors sample vectors to find the neighbours within r
+ * for
+ * @param r radius within which to count points
+ * @param dynCorrExclTime time window within which to exclude
+ * points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @return the count of points within r.
+ */
+ public abstract int countPointsWithinR(double[][] sampleVectors, double r,
+ boolean allowEqualToR);
+
/**
* As per {@link #countPointsWithinR(int, double, boolean)}
* however returns a collection rather than a count.
@@ -331,6 +351,60 @@ public abstract class NearestNeighbourSearcher {
int sampleIndex, double r, int dynCorrExclTime,
boolean allowEqualToR, boolean[] isWithinR, int[] indicesWithinR);
+ /**
+ * As per {@link #findPointsWithinR(int, double, boolean, boolean[], int[])}
+ * however incorporates dynamic correlation exclusion.
+ *
+ *
+ * @param sampleIndex sample index in the data to find a nearest neighbour
+ * for
+ * @param r radius within which to count points
+ * @param dynCorrExclTime time window within which to exclude
+ * points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param isWithinR the array MUST be passed in with all points set to
+ * false initially, and is returned indicating whether each sample was
+ * found to be within r of that at sampleIndex.
+ * @param distancesWithinRInOrder the array must be passed in and is
+ * returned with distances for each point found to be within r. Values
+ * at other indices are not defined.
+ * @param distancesAndIndicesWithinR is returned as
+ * a list of distances (in column index 0)
+ * and array indices (in column index 1)
+ * for points marked as true in isWithinR, terminated with a -1 value on the index.
+ * @return the point count
+ */
+ public abstract int findPointsWithinR(
+ int sampleIndex, double r, int dynCorrExclTime,
+ boolean allowEqualToR, boolean[] isWithinR,
+ double[] distancesWithinRInOrder,
+ double[][] distancesAndIndicesWithinR);
+
+ /**
+ * As per {@link #countPointsWithinR(int, double, boolean)}
+ * however records the nearest neighbours for a sample data point
+ * (which may not be in the search tree), within the isWithinR
+ * and indicesWithinR arrays, which must be constructed before
+ * calling this method, with length at or exceeding the total
+ * number of data points. indicesWithinR is
+ *
+ *
+ * @param r radius within which to count points
+ * @param sampleVectors sample vectors to find the neighbours within r
+ * for
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param isWithinR the array MUST be passed in with all points set to
+ * false initially, and is returned indicating whether each sample was
+ * found to be within r of that at sampleIndex.
+ * @param indicesWithinR a list of array indices
+ * for points marked as true in isWithinR, terminated with a -1 value.
+ */
+ public abstract void findPointsWithinR(
+ double r, double[][] sampleVectors,
+ boolean allowEqualToR, boolean[] isWithinR, int[] indicesWithinR);
+
/**
* As per {@link #countPointsWithinR(int, double, int, boolean)}
* with allowEqualToR == false
@@ -472,5 +546,20 @@ public abstract class NearestNeighbourSearcher {
public abstract int countPointsWithinR(int sampleIndex, double r,
boolean allowEqualToR, boolean[] additionalCriteria);
-
+ /**
+ * As per {@link #countPointsWithinR(double[][], double, boolean)}
+ * however each point is subject to also meeting the additional
+ * criteria of being true in additionalCriteria.
+ *
+ * @param sampleVectors sample vectors to find the neighbours within r
+ * 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
+ * @param additionalCriteria array of booleans. Only count a point if it
+ * is within r and is true in additionalCrtieria.
+ * @return the count of points within r.
+ */
+ public abstract int countPointsWithinR(double[][] sampleVectors, double r,
+ boolean allowEqualToR, boolean[] additionalCriteria);
}
diff --git a/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java b/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java
index d56d2be..73e7b58 100755
--- a/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java
+++ b/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java
@@ -18,6 +18,7 @@
package infodynamics.utils;
+import java.util.Arrays;
import java.util.Collection;
import java.util.PriorityQueue;
import java.util.Vector;
@@ -43,6 +44,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
* Number of samples in the data
*/
protected int numObservations = 0;
+
/**
* An array of indices to the data in originalDataSet,
* sorted in order (min to max).
@@ -55,6 +57,11 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
*/
protected int[] indicesInSortedArray = null;
+ /**
+ * An array of the original values sorted into ascending order.
+ */
+ protected double[] sortedValues = null;
+
public UnivariateNearestNeighbourSearcher(double[][] data) throws Exception {
// Ideally we would not call the constructor until after the following check,
// but the constructor must come first in Java.
@@ -88,6 +95,13 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
sortedArrayIndices[t] = (int) dataWithIndices[t][1];
indicesInSortedArray[(int) dataWithIndices[t][1]] = t;
}
+ // Finally, create the sorted array of the values:
+ // (these are only used by a few methods, but better to create the array now
+ // rather than hit concurrency issues trying to create it when required
+ sortedValues = new double[numObservations];
+ for (int i = 0; i < numObservations; i++) {
+ sortedValues[i] = originalDataSet[sortedArrayIndices[i]];
+ }
}
/**
@@ -359,6 +373,47 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
return count;
}
+ /**
+ * Count the number of points within norms r_upper and r_lower for a given
+ * sample index in the data set. The node itself is
+ * excluded from the search.
+ * Nearest neighbour function to compare to radii 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 countPointsWithinRs(int sampleIndex, double r_upper,
+ double r_lower, boolean allowEqualToR) {
+ int count = 0;
+ // 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_lower)) ||
+ (!allowEqualToR && (theNorm < r_lower))) {
+ count++;
+ 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_upper)) ||
+ (!allowEqualToR && (theNorm < r_upper))) {
+ count++;
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ return count;
+ }
+
/**
* Count the number of points within norm r for a given
* sample index in the data set. Nodes within dynCorrExclTime
@@ -408,6 +463,124 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
return count;
}
+ /**
+ * Count the number of points within norm r for a given
+ * sample index in the data set, or larger than r. 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 countPointsWithinROrLarger(int sampleIndex, double r, boolean allowEqualToR) {
+ int count = 0;
+ // 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))) {
+ count++;
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ count += numObservations - indexInSortedArray - 1;
+ return count;
+ }
+
+ /**
+ * Count the number of points less than the sample at the given index,
+ * and outside norm r. allowEqualToR means to be considered same as the sample.
+ * 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 countPointsSmallerAndOutsideR(int sampleIndex, double r, boolean allowEqualToR) {
+ // Find where this node sits in the sorted array:
+ int indexInSortedArray = indicesInSortedArray[sampleIndex];
+ // Check the points with smaller data values first:
+ int i;
+ for (i = indexInSortedArray - 1; i >= 0; i--) {
+ double theNorm = norm(originalDataSet[sampleIndex],
+ originalDataSet[sortedArrayIndices[i]], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ return i + 1;
+ }
+
+ @Override
+ public int countPointsWithinR(double[][] sampleVectors, double r, boolean allowEqualToR) {
+ if (sampleVectors.length > 1) {
+ throw new RuntimeException("sampleVectors[][] argument must have length 1 on each dimension when calling a UnivariateNearstNeighbourSearcher");
+ }
+ if (sampleVectors[0].length > 1) {
+ throw new RuntimeException("sampleVectors[][] argument must have length 1 on each dimension when calling a UnivariateNearstNeighbourSearcher");
+ }
+
+ return countPointsWithinR(sampleVectors[0][0], r, allowEqualToR);
+ }
+
+ /**
+ * As per {@link #countPointsWithinR(double, double, boolean)} however
+ * for a specific sample value (not one within the data set itself).
+ *
+ * @param sampleValue the specific value to search for
+ */
+ public int countPointsWithinR(double sampleValue, double r, boolean allowEqualToR) {
+
+ int count = 0;
+
+ // Find where this value would sit in the sorted array.
+ int potentialArrayIndex = Arrays.binarySearch(sortedValues, sampleValue);
+ if (potentialArrayIndex < 0) {
+ // The sampleValue was not found because it wasn't in the original set of values:
+ // invert the return value of binarySearch to find where the sampleValue would be inserted into the array:
+ potentialArrayIndex = -potentialArrayIndex - 1;
+ // This means at potentialArrayIndex and higher are values >= sampleValue, while
+ // at potentialArrayIndex-1 and lower are values <= sampleValue.
+ // potentialArrayIndex will be array.length if all values are smaller than it.
+ } // else potentialArrayIndex is pointing to one instance of this value in the sorted array
+
+ // Check the points with smaller data values first
+ for (int i = potentialArrayIndex - 1; i >= 0; i--) {
+ double theNorm = norm(sampleValue,
+ sortedValues[i], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ count++;
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Next check the points with larger data values: (and including potentialArrayIndex, since we're
+ // not excluding it since it's not nominally part of the underlying data set)
+ for (int i = potentialArrayIndex; i < numObservations; i++) {
+ double theNorm = norm(sampleValue,
+ sortedValues[i], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ count++;
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ return count;
+ }
+
/**
* Find the collection of points within norm r for a given
* sample index in the data set. The node itself is
@@ -565,6 +738,214 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
indicesWithinR[indexInIndicesWithinR++] = -1;
}
+ @Override
+ public int findPointsWithinR(int sampleIndex, double r,
+ int dynCorrExclTime, boolean allowEqualToR, boolean[] isWithinR,
+ double[] distancesWithinRInOrder,
+ double[][] distancesAndIndicesWithinR) {
+ int indexInIndicesWithinR = 0;
+ // 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--) {
+ if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+ // Can't count this point, but keep checking:
+ continue;
+ }
+ double theNorm = norm(originalDataSet[sampleIndex],
+ originalDataSet[sortedArrayIndices[i]], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ isWithinR[sortedArrayIndices[i]] = true;
+ distancesWithinRInOrder[sortedArrayIndices[i]] = theNorm;
+ distancesAndIndicesWithinR[indexInIndicesWithinR][0] = theNorm;
+ distancesAndIndicesWithinR[indexInIndicesWithinR++][1] = sortedArrayIndices[i];
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Next check the points with larger data values:
+ for (int i = indexInSortedArray + 1; i < numObservations; i++) {
+ if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+ // Can't count this point, but keep checking:
+ continue;
+ }
+ double theNorm = norm(originalDataSet[sampleIndex],
+ originalDataSet[sortedArrayIndices[i]], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ isWithinR[sortedArrayIndices[i]] = true;
+ distancesWithinRInOrder[sortedArrayIndices[i]] = theNorm;
+ distancesAndIndicesWithinR[indexInIndicesWithinR][0] = theNorm;
+ distancesAndIndicesWithinR[indexInIndicesWithinR++][1] = sortedArrayIndices[i];
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Write the terminating integer into the indicesWithinR array:
+ distancesAndIndicesWithinR[indexInIndicesWithinR++][1] = -1;
+ return indexInIndicesWithinR - 1;
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, double, int, boolean, boolean[], double[], double[][])}
+ * only we have additional criteria, and we report max values taking that additional
+ * criteria into account.
+ *
+ * @param sampleIndex sample index in the data to find a nearest neighbour
+ * for
+ * @param r radius within which to count points
+ * @param dynCorrExclTime time window within which to exclude
+ * points to be counted. Is >= 0. 0 means only exclude sampleIndex.
+ * @param allowEqualToR if true, then count points at radius r also,
+ * otherwise only those strictly within r
+ * @param additionalCriteria an array indicating whether each sample for
+ * should be tested for this search point
+ * @param reportDistancesAsMaxWith an array of distances for
+ * points flagged in additionalCriteria: distancesAndIndicesWithinR should
+ * return the max of this value or the distances within r computed here.
+ * Values at indices where distances are not within r are not defined.
+ * @param distancesAndIndicesWithinR a list of distances (in column index 0)
+ * and array indices (in column index 1)
+ * for points found to be within r, terminated with a -1 value on the index.
+ * @return the next available index in distancesAndIndicesWithinR after the method is complete
+ */
+ public int findPointsWithinRAndTakeRadiusMax(int sampleIndex, double r,
+ int dynCorrExclTime, boolean allowEqualToR,
+ boolean[] additionalCriteria,
+ double[] reportDistancesAsMaxWith,
+ double[][] distancesAndIndicesWithinR) {
+ int indexInIndicesWithinR = 0;
+ // 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--) {
+ if (!additionalCriteria[sortedArrayIndices[i]] ||
+ (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime)) {
+ // Can't count this point, but keep checking:
+ continue;
+ }
+ double theNorm = norm(originalDataSet[sampleIndex],
+ originalDataSet[sortedArrayIndices[i]], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ if (theNorm < reportDistancesAsMaxWith[sortedArrayIndices[i]]) {
+ theNorm = reportDistancesAsMaxWith[sortedArrayIndices[i]];
+ }
+ distancesAndIndicesWithinR[indexInIndicesWithinR][0] = theNorm;
+ distancesAndIndicesWithinR[indexInIndicesWithinR++][1] = sortedArrayIndices[i];
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Next check the points with larger data values:
+ for (int i = indexInSortedArray + 1; i < numObservations; i++) {
+ if (!additionalCriteria[sortedArrayIndices[i]] ||
+ (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime)) {
+ // Can't count this point, but keep checking:
+ continue;
+ }
+ double theNorm = norm(originalDataSet[sampleIndex],
+ originalDataSet[sortedArrayIndices[i]], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ if (theNorm < reportDistancesAsMaxWith[sortedArrayIndices[i]]) {
+ theNorm = reportDistancesAsMaxWith[sortedArrayIndices[i]];
+ }
+ distancesAndIndicesWithinR[indexInIndicesWithinR][0] = theNorm;
+ distancesAndIndicesWithinR[indexInIndicesWithinR++][1] = sortedArrayIndices[i];
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Write the terminating integer into the indicesWithinR array:
+ distancesAndIndicesWithinR[indexInIndicesWithinR++][1] = -1;
+ return indexInIndicesWithinR - 1;
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, double, int, boolean, boolean[], int[])}
+ * only for a specific data point.
+ *
+ * Note that the method signature takes a 2D sample vector, however since this is class
+ * only searches univariate spaces, then there must only be one value in #sampleVector
+ * (i.e. length one on both array indices).
+ *
+ *
+ * @see super{@link #findPointsWithinR(double, double[][], boolean, boolean[], int[])}
+ */
+ @Override
+ public void findPointsWithinR(double r, double[][] sampleVectors,
+ boolean allowEqualToR, boolean[] isWithinR, int[] indicesWithinR) {
+
+ if (sampleVectors.length > 1) {
+ throw new RuntimeException("sampleVectors[][] argument must have length 1 on each dimension when calling a UnivariateNearstNeighbourSearcher");
+ }
+ if (sampleVectors[0].length > 1) {
+ throw new RuntimeException("sampleVectors[][] argument must have length 1 on each dimension when calling a UnivariateNearstNeighbourSearcher");
+ }
+
+ findPointsWithinR(r, sampleVectors[0][0], allowEqualToR, isWithinR, indicesWithinR);
+ }
+
+ /**
+ * As per {@link #findPointsWithinR(int, double, int, boolean, boolean[], int[])}
+ * only for a specific data point.
+ *
+ * @see super{@link #findPointsWithinR(double, double[][], boolean, boolean[], int[])}
+ *
+ * @param sampleValue the specific point being searched (not a data point in the search structure)
+ */
+ public void findPointsWithinR(double r, double sampleValue,
+ boolean allowEqualToR, boolean[] isWithinR, int[] indicesWithinR) {
+ int indexInIndicesWithinR = 0;
+
+ // Find where this value would sit in the sorted array.
+ int potentialArrayIndex = Arrays.binarySearch(sortedValues, sampleValue);
+ if (potentialArrayIndex < 0) {
+ // The sampleValue was not found because it wasn't in the original set of values:
+ // invert the return value of binarySearch to find where the sampleValue would be inserted into the array:
+ potentialArrayIndex = -potentialArrayIndex - 1;
+ // This means at potentialArrayIndex and higher are values >= sampleValue, while
+ // at potentialArrayIndex-1 and lower are values <= sampleValue.
+ // potentialArrayIndex will be array.length if all values are smaller than it.
+ } // else potentialArrayIndex is pointing to one instance of this value in the sorted array
+
+ // Check the points with smaller data values first:
+ for (int i = potentialArrayIndex - 1; i >= 0; i--) {
+ double theNorm = norm(sampleValue,
+ sortedValues[i], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ isWithinR[sortedArrayIndices[i]] = true;
+ indicesWithinR[indexInIndicesWithinR++] = sortedArrayIndices[i];
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Next check the points with larger data values: (and including potentialArrayIndex, since we're
+ // not excluding it since it's not nominally part of the underlying data set)
+ for (int i = potentialArrayIndex; i < numObservations; i++) {
+ double theNorm = norm(sampleValue,
+ sortedValues[i], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ isWithinR[sortedArrayIndices[i]] = true;
+ indicesWithinR[indexInIndicesWithinR++] = sortedArrayIndices[i];
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Write the terminating integer into the indicesWithinR array:
+ indicesWithinR[indexInIndicesWithinR++] = -1;
+ }
+
@Override
public int countPointsWithinR(int sampleIndex, double r, boolean allowEqualToR,
boolean[] additionalCriteria) {
@@ -618,4 +999,103 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
}
return count;
}
+
+ /**
+ * As per {@link #countPointsWithinR(int, double, boolean, boolean[])}
+ * only for a specific data point.
+ *
+ * Note that the method signature takes a 2D sample vector, however since this is class
+ * only searches univariate spaces, then there must only be one value in #sampleVector
+ * (i.e. length one on both array indices).
+ *
+ *
+ * @see super{@link #findPointsWithinR(double, double[][], boolean, boolean[], int[])}
+ */
+ @Override
+ public int countPointsWithinR(double[][] sampleVectors, double r, boolean allowEqualToR,
+ boolean[] additionalCriteria) {
+
+ if (sampleVectors.length > 1) {
+ throw new RuntimeException("sampleVectors[][] argument must have length 1 on each dimension when calling a UnivariateNearstNeighbourSearcher");
+ }
+ if (sampleVectors[0].length > 1) {
+ throw new RuntimeException("sampleVectors[][] argument must have length 1 on each dimension when calling a UnivariateNearstNeighbourSearcher");
+ }
+
+ return countPointsWithinR(sampleVectors[0][0], r, allowEqualToR, additionalCriteria);
+ }
+
+ /**
+ * As per {@link #countPointsWithinR(int, double, boolean, boolean[])}
+ * only for a specific data point.
+ *
+ * @param sampleValue the specific point being searched (not a data point in the search structure)
+ */
+ public int countPointsWithinR(double sampleValue, double r, boolean allowEqualToR,
+ boolean[] additionalCriteria) {
+
+ int count = 0;
+
+ // Find where this value would sit in the sorted array.
+ int potentialArrayIndex = Arrays.binarySearch(sortedValues, sampleValue);
+ if (potentialArrayIndex < 0) {
+ // The sampleValue was not found because it wasn't in the original set of values:
+ // invert the return value of binarySearch to find where the sampleValue would be inserted into the array:
+ potentialArrayIndex = -potentialArrayIndex - 1;
+ // This means at potentialArrayIndex and higher are values >= sampleValue, while
+ // at potentialArrayIndex-1 and lower are values <= sampleValue.
+ // potentialArrayIndex will be array.length if all values are smaller than it.
+ } // else potentialArrayIndex is pointing to one instance of this value in the sorted array
+
+ // Check the points with smaller data values first:
+ for (int i = potentialArrayIndex - 1; i >= 0; i--) {
+ if (!additionalCriteria[sortedArrayIndices[i]]) {
+ // This point failed the additional criteria,
+ // so skip checking it.
+ // We check this first, even though it's the norm
+ // that really determines whether we need to continue
+ // checking or not. In some circumstances this may be
+ // slower, but for our main application - KSG conditional MI -
+ // this should be faster.
+ continue;
+ }
+ double theNorm = norm(sampleValue,
+ sortedValues[i], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ count++;
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ // Next check the points with larger data values: (and including potentialArrayIndex, since we're
+ // not excluding it since it's not nominally part of the underlying data set)
+ for (int i = potentialArrayIndex; i < numObservations; i++) {
+ if (!additionalCriteria[sortedArrayIndices[i]]) {
+ // This point failed the additional criteria,
+ // so skip checking it.
+ // We check this first, even though it's the norm
+ // that really determines whether we need to continue
+ // checking or not. In some circumstances this may be
+ // slower, but for our main application - KSG conditional MI -
+ // this should be faster.
+ continue;
+ }
+ double theNorm = norm(sampleValue,
+ sortedValues[i], normTypeToUse);
+ if ((allowEqualToR && (theNorm <= r)) ||
+ (!allowEqualToR && (theNorm < r))) {
+ count++;
+ continue;
+ }
+ // Else no point checking further points
+ break;
+ }
+ return count;
+ }
+
+ public int getNumObservations() {
+ return numObservations;
+ }
}
diff --git a/java/unittests/infodynamics/utils/KdTreeTest.java b/java/unittests/infodynamics/utils/KdTreeTest.java
index baa3bfb..e4c28f7 100755
--- a/java/unittests/infodynamics/utils/KdTreeTest.java
+++ b/java/unittests/infodynamics/utils/KdTreeTest.java
@@ -502,6 +502,48 @@ public class KdTreeTest extends TestCase {
long endTimeValidate = Calendar.getInstance().getTimeInMillis();
System.out.printf("All %d nearest neighbours found in: %.3f sec\n",
K, ((double) (endTimeValidate - startTime)/1000.0));
+
+ // And also test the K nearest neighbours for a few points that
+ // weren't in the original data set:
+ for (int n = 0; n < 100; n++) {
+ // Generate a sample vector
+ double[][] sampleVectors = new double[variables][];
+ for (int v = 0; v < variables; v++) {
+ int vLength = data[v][0].length;
+ sampleVectors[v] = rg.generateNormalData(vLength, 0, 1);
+ }
+ // Now pull the kNNs for it:
+ PriorityQueue nnPQ =
+ kdTree.findKNearestNeighbours(K, sampleVectors);
+ assertTrue(nnPQ.size() == K);
+ // Now find the K nearest neighbours with a naive all-pairs comparison
+ double[][] distancesAndIndices = new double[samples][2];
+ for (int t2 = 0; t2 < samples; t2++) {
+ double maxNorm = 0;
+ for (int v = 0; v < variables; v++) {
+ double normForThisVariable = normCalculator.norm(
+ sampleVectors[v], data[v][t2]);
+ if (maxNorm < normForThisVariable) {
+ maxNorm = normForThisVariable;
+ }
+ }
+ distancesAndIndices[t2][0] = maxNorm;
+ distancesAndIndices[t2][1] = t2;
+ }
+ int[] timeStepsOfKthMins =
+ MatrixUtils.kMinIndices(distancesAndIndices, 0, K);
+ for (int i = 0; i < K; i++) {
+ // Check that the ith nearest neighbour matches
+ NeighbourNodeData nnData = nnPQ.poll();
+ if (timeStepsOfKthMins[K - 1 - i] != nnData.sampleIndex) {
+ // We have an error:
+ System.out.printf("Erroneous match between indices %d (expected) " +
+ " and %d\n", timeStepsOfKthMins[K - 1 - i], nnData.sampleIndex);
+ }
+ assertEquals(timeStepsOfKthMins[K - 1 - i], nnData.sampleIndex);
+ }
+ }
+
}
public void testCountNeighboursWithinRSeparateArrays() {
@@ -534,6 +576,16 @@ public class KdTreeTest extends TestCase {
data, samples, kdTree, rs[i], true);
verifyCountNeighboursWithinRForSeparateArraysArrayArgs(EuclideanUtils.NORM_MAX_NORM,
data, samples, kdTree, rs[i], false);
+ // And test the array return variant for radii returned:
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsPlusRadiusOut(EuclideanUtils.NORM_MAX_NORM,
+ data, samples, kdTree, rs[i], true);
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsPlusRadiusOut(EuclideanUtils.NORM_MAX_NORM,
+ data, samples, kdTree, rs[i], false);
+ // And test the array return variant for new data points too:
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsNewPoints(EuclideanUtils.NORM_MAX_NORM,
+ data, samples, kdTree, rs[i], true);
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsNewPoints(EuclideanUtils.NORM_MAX_NORM,
+ data, samples, kdTree, rs[i], false);
// Must test with EuclideanUtils.NORM_EUCLIDEAN_SQUARED instead
// of EuclideanUtils.NORM_EUCLIDEAN if we want the min distances
// to match when tested inside verifyNearestNeighbourForSeparateArrays()
@@ -545,6 +597,14 @@ public class KdTreeTest extends TestCase {
data, samples, kdTree, rs[i], true);
verifyCountNeighboursWithinRForSeparateArraysArrayArgs(EuclideanUtils.NORM_EUCLIDEAN_SQUARED,
data, samples, kdTree, rs[i], false);
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsPlusRadiusOut(EuclideanUtils.NORM_EUCLIDEAN_SQUARED,
+ data, samples, kdTree, rs[i], true);
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsPlusRadiusOut(EuclideanUtils.NORM_EUCLIDEAN_SQUARED,
+ data, samples, kdTree, rs[i], false);
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsNewPoints(EuclideanUtils.NORM_EUCLIDEAN_SQUARED,
+ data, samples, kdTree, rs[i], true);
+ verifyCountNeighboursWithinRForSeparateArraysArrayArgsNewPoints(EuclideanUtils.NORM_EUCLIDEAN_SQUARED,
+ data, samples, kdTree, rs[i], false);
}
}
@@ -694,6 +754,198 @@ public class KdTreeTest extends TestCase {
((double) (endTimeValidate - endTimeCount)/1000.0));
}
+ private void verifyCountNeighboursWithinRForSeparateArraysArrayArgsPlusRadiusOut(int normType,
+ double[][][] data, int samples, KdTree kdTree, double r,
+ boolean allowEqualToR) {
+
+ int variables = data.length;
+
+ // Set up the given norm type:
+ EuclideanUtils normCalculator = new EuclideanUtils(normType);
+ kdTree.setNormType(normType);
+
+ long startTime = Calendar.getInstance().getTimeInMillis();
+ int totalCount = 0;
+ int[] counts = new int[samples];
+ boolean[] withinR = new boolean[samples];
+ double[] distancesWithinRInOrder = new double[samples];
+ double[][] distancesAndIndicesWithinR = new double[samples][2];
+ for (int t = 0; t < samples; t++) {
+ // int count = kdTree.countPointsWithinR(t, r, allowEqualToR);
+ kdTree.findPointsWithinR(t, r, 0, allowEqualToR, withinR,
+ distancesWithinRInOrder, distancesAndIndicesWithinR);
+ int count = 0;
+ // Run through the list of points returned as within r and check them
+ for (int t2 = 0; distancesAndIndicesWithinR[t2][1] >= 0; t2++) {
+ count++;
+ assertTrue(withinR[(int) distancesAndIndicesWithinR[t2][1]]);
+ // Reset this marker
+ withinR[(int) distancesAndIndicesWithinR[t2][1]] = false;
+ // Check that it's really within r
+ double maxNorm = 0;
+ for (int v = 0; v < variables; v++) {
+ double normForThisVariable = normCalculator.norm(
+ data[v][t], data[v][(int) distancesAndIndicesWithinR[t2][1]]);
+ if (maxNorm < normForThisVariable) {
+ maxNorm = normForThisVariable;
+ }
+ }
+
+ assertTrue((!allowEqualToR && (maxNorm < r)) ||
+ ( allowEqualToR && (maxNorm <= r)));
+ // Check that its distance was properly returned in the other arrays:
+ assertEquals(maxNorm, distancesAndIndicesWithinR[t2][0], 0.000001);
+ assertEquals(distancesAndIndicesWithinR[t2][0],
+ distancesWithinRInOrder[(int) distancesAndIndicesWithinR[t2][1]],
+ 0.00000001);
+ }
+ // Check that there were no other points marked as within r:
+ int count2 = 0;
+ for (int t2 = 0; t2 < samples; t2++) {
+ if (withinR[t2]) {
+ count2++;
+ }
+ }
+ assertEquals(0, count2); // We set all the ones that should have been true to false already
+ // Now reset the boolean array
+ counts[t] = count;
+ totalCount += count;
+ }
+ long endTimeCount = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("All neighbours within %.3f (average of %.3f) found in: %.3f sec\n",
+ r, (double) totalCount / (double) samples,
+ ((double) (endTimeCount - startTime)/1000.0));
+ // Now find the neighbour count with a naive all-pairs comparison
+ for (int t = 0; t < samples; t++) {
+ int naiveCount = 0;
+ for (int t2 = 0; t2 < samples; t2++) {
+ double norm = Double.POSITIVE_INFINITY;
+ if (t2 != t) {
+ double maxNorm = 0;
+ for (int v = 0; v < variables; v++) {
+ double normForThisVariable = normCalculator.norm(
+ data[v][t], data[v][t2]);
+ if (maxNorm < normForThisVariable) {
+ maxNorm = normForThisVariable;
+ }
+ }
+ norm = maxNorm;
+ }
+ if ((!allowEqualToR && (norm < r)) ||
+ ( allowEqualToR && (norm <= r))) {
+ naiveCount++;
+ }
+ }
+ if (naiveCount != counts[t]) {
+ System.out.printf("All pairs : count %d\n", naiveCount);
+ System.out.printf("kdTree search: count %d\n", counts[t]);
+ }
+ assertEquals(naiveCount, counts[t]);
+ }
+ long endTimeValidate = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("All neighbours within %.3f (average of %.3f) validated in: %.3f sec\n",
+ r, (double) totalCount / (double) samples,
+ ((double) (endTimeValidate - endTimeCount)/1000.0));
+ }
+
+ private void verifyCountNeighboursWithinRForSeparateArraysArrayArgsNewPoints(int normType,
+ double[][][] data, int samples, KdTree kdTree, double r,
+ boolean allowEqualToR) {
+
+ int variables = data.length;
+
+ int numNewSamplesToCheck = 100;
+ double[][][] allSampleVectors = new double[numNewSamplesToCheck][][];
+
+ // Set up the given norm type:
+ EuclideanUtils normCalculator = new EuclideanUtils(normType);
+ kdTree.setNormType(normType);
+
+ long startTime = Calendar.getInstance().getTimeInMillis();
+ int totalCount = 0;
+ int[] counts = new int[samples];
+ boolean[] withinR = new boolean[samples];
+ int[] indicesWithinR = new int[samples];
+
+ // Test the range search for a few points that
+ // weren't in the original data set:
+ for (int n = 0; n < numNewSamplesToCheck; n++) {
+ // Generate a sample vector
+ double[][] sampleVectors = new double[variables][];
+ for (int v = 0; v < variables; v++) {
+ int vLength = data[v][0].length;
+ sampleVectors[v] = rg.generateNormalData(vLength, 0, 1);
+ }
+ allSampleVectors[n] = sampleVectors;
+
+ kdTree.findPointsWithinR(r, sampleVectors, allowEqualToR, withinR, indicesWithinR);
+ int count = 0;
+ // Run through the list of points returned as within r and check them
+ for (int t2 = 0; indicesWithinR[t2] != -1; t2++) {
+ count++;
+ assertTrue(withinR[indicesWithinR[t2]]);
+ // Reset this marker
+ withinR[indicesWithinR[t2]] = false;
+ // Check that it's really within r
+ double maxNorm = 0;
+ for (int v = 0; v < variables; v++) {
+ double normForThisVariable = normCalculator.norm(
+ sampleVectors[v], data[v][indicesWithinR[t2]]);
+ if (maxNorm < normForThisVariable) {
+ maxNorm = normForThisVariable;
+ }
+ }
+
+ assertTrue((!allowEqualToR && (maxNorm < r)) ||
+ ( allowEqualToR && (maxNorm <= r)));
+ }
+ // Check that there were no other points marked as within r:
+ int count2 = 0;
+ for (int t2 = 0; t2 < samples; t2++) {
+ if (withinR[t2]) {
+ count2++;
+ }
+ }
+ assertEquals(0, count2); // We set all the ones that should have been true to false already
+ // Now reset the boolean array
+ counts[n] = count;
+ totalCount += count;
+ }
+ long endTimeCount = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("All neighbours within %.3f (average of %.3f) found in: %.3f sec\n",
+ r, (double) totalCount / (double) samples,
+ ((double) (endTimeCount - startTime)/1000.0));
+ // Now find the neighbour count with a naive all-pairs comparison
+ for (int n = 0; n < numNewSamplesToCheck; n++) {
+ int naiveCount = 0;
+ for (int t2 = 0; t2 < samples; t2++) {
+ double norm = Double.POSITIVE_INFINITY;
+ double maxNorm = 0;
+ for (int v = 0; v < variables; v++) {
+ double normForThisVariable = normCalculator.norm(
+ allSampleVectors[n][v], data[v][t2]);
+ if (maxNorm < normForThisVariable) {
+ maxNorm = normForThisVariable;
+ }
+ }
+ norm = maxNorm;
+ if ((!allowEqualToR && (norm < r)) ||
+ ( allowEqualToR && (norm <= r))) {
+ naiveCount++;
+ }
+ }
+ if (naiveCount != counts[n]) {
+ System.out.printf("All pairs : count %d\n", naiveCount);
+ System.out.printf("kdTree search: count %d\n", counts[n]);
+ }
+ assertEquals(naiveCount, counts[n]);
+ }
+ long endTimeValidate = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("All neighbours within %.3f (average of %.3f) validated in: %.3f sec\n",
+ r, (double) totalCount / (double) samples,
+ ((double) (endTimeValidate - endTimeCount)/1000.0));
+ }
+
public void testCountNeighboursWithinRExclusionWindow() {
int dimensionsPerVariable = 3;
int samples = 1000;
@@ -1106,6 +1358,42 @@ public class KdTreeTest extends TestCase {
}
}
}
+
+ // And now test some fresh samples:
+ int numNewSamples = 1000;
+ for (int i = 0; i < r1.length; i++) {
+
+ for (int n = 0; n < numNewSamples; n++) {
+ double[][] newSample = rg.generateNormalData(variables, dimensionsPerVariable, 0, 1);
+ double[] variable1 = newSample[0];
+ double[][] newSampleSansVariable1 = new double[variables - 1][];
+ for (int v = 1; v < variables; v++) {
+ newSampleSansVariable1[v-1] = newSample[v];
+ }
+
+ // Establish our ground truth by searching all variables:
+ int count = kdTree.countPointsWithinR(newSample, r1[i], allowEqualsR);
+
+ // Now search in variable 1 first, then use that
+ // to help the others:
+ nnSearcherVar1.findPointsWithinR(r1[i], new double[][] {variable1}, allowEqualsR,
+ isWithinR, indicesWithinR);
+ int countUsingVar1 = kdTree.countPointsWithinR(r1[i], newSample,
+ allowEqualsR, 0, isWithinR);
+ assertEquals(count, countUsingVar1);
+
+ // And then search the kdTree which does not have variable 1:
+ int countAdditionalCriteria =
+ nnSearcherExceptVar1.countPointsWithinR(newSampleSansVariable1, r1[i],
+ allowEqualsR, isWithinR);
+ assertEquals(count, countAdditionalCriteria);
+
+ // Finally, reset our boolean array:
+ for (int t2 = 0; indicesWithinR[t2] != -1; t2++) {
+ isWithinR[indicesWithinR[t2]] = false;
+ }
+ }
+ }
}
public void testPreviouslyTestedVariableMethodMultipleRs() throws Exception {
diff --git a/java/unittests/infodynamics/utils/UnivariateNearestNeighbourTest.java b/java/unittests/infodynamics/utils/UnivariateNearestNeighbourTest.java
index 03b38a8..bcd1b74 100755
--- a/java/unittests/infodynamics/utils/UnivariateNearestNeighbourTest.java
+++ b/java/unittests/infodynamics/utils/UnivariateNearestNeighbourTest.java
@@ -117,6 +117,21 @@ public class UnivariateNearestNeighbourTest extends TestCase {
checkRangeFinderArrayCalls(false, 100);
}
+ public void testRangeFinderArrayCallsRaidusOut() throws Exception {
+ checkRangeFinderArrayCallsPlusRadiusOut(true, 100);
+ checkRangeFinderArrayCallsPlusRadiusOut(false, 100);
+ }
+
+ public void testRangeFinderArrayCallsMaxRadiusOut() throws Exception {
+ checkRangeFinderArrayCallsPlusMaxRadiusOut(true, 100);
+ checkRangeFinderArrayCallsPlusMaxRadiusOut(false, 100);
+ }
+
+ public void testRangeFinderArrayCallsNewData() throws Exception {
+ checkRangeFinderArrayCallsNewData(true);
+ checkRangeFinderArrayCallsNewData(false);
+ }
+
public void checkRangeFinder(boolean strict, int exclWindow) throws Exception {
int numTimeStepsInitial = 10000;
int duplicateSteps = 100;
@@ -242,6 +257,238 @@ public class UnivariateNearestNeighbourTest extends TestCase {
}
}
+ public void checkRangeFinderArrayCallsPlusRadiusOut(boolean strict, int exclWindow) throws Exception {
+ int numTimeStepsInitial = 10000;
+ int duplicateSteps = 100;
+ int numTimeSteps = numTimeStepsInitial+duplicateSteps;
+ double[] dataRaw = rg.generateNormalData(numTimeStepsInitial, 0, 1);
+ double[] data = new double[numTimeSteps];
+ // Now duplicate some of the time steps as a test:
+ System.arraycopy(dataRaw, 0, data, 0, numTimeStepsInitial);
+ System.arraycopy(dataRaw, 0, data, numTimeStepsInitial, duplicateSteps);
+
+ UnivariateNearestNeighbourSearcher searcher =
+ new UnivariateNearestNeighbourSearcher(data);
+ int[] counts = new int[data.length];
+ boolean[] withinR = new boolean[data.length];
+ double[] distancesWithinRInOrder = new double[data.length];
+ double[][] distancesAndIndicesWithinR = new double[data.length][2];
+ double r = 0.2;
+ long startTime = Calendar.getInstance().getTimeInMillis();
+ for (int t = 0; t < data.length; t++) {
+ searcher.findPointsWithinR(t, r, exclWindow, !strict,
+ withinR, distancesWithinRInOrder, distancesAndIndicesWithinR);
+ int count = 0;
+ // Run through the list of points returned as within r and check them
+ for (int t2 = 0; distancesAndIndicesWithinR[t2][1] >= 0; t2++) {
+ count++;
+ assertTrue(withinR[(int) distancesAndIndicesWithinR[t2][1]]);
+ // Reset this marker
+ withinR[(int) distancesAndIndicesWithinR[t2][1]] = false;
+ // Check that it's really within r
+ double maxNorm = Math.abs(data[t] - data[(int) distancesAndIndicesWithinR[t2][1]]);
+ assertTrue((strict && (maxNorm < r)) ||
+ ( !strict && (maxNorm <= r)));
+ // Check that its distance was properly returned in the other arrays:
+ assertEquals(maxNorm, distancesAndIndicesWithinR[t2][0], 0.000001);
+ assertEquals(distancesAndIndicesWithinR[t2][0],
+ distancesWithinRInOrder[(int) distancesAndIndicesWithinR[t2][1]],
+ 0.00000001);
+ }
+ // Check that there were no other points marked as within r in withinR:
+ int count2 = 0;
+ for (int t2 = 0; t2 < data.length; t2++) {
+ if (withinR[t2]) {
+ count2++;
+ }
+ }
+ assertEquals(0, count2); // We set all the ones that should have been true to false already
+ // Now reset the boolean array
+ counts[t] = count;
+ }
+ long endTimeNNs = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("Found all neighbours within %.3f (mean = %.3f) in: %.3f sec\n",
+ r, MatrixUtils.mean(counts),
+ ((double) (endTimeNNs - startTime)/1000.0));
+
+ // Now do brute force:
+ int[] bruteForceCounts = new int[data.length];
+ for (int t = 0; t < data.length; t++) {
+ int count = 0;
+ for (int t2 = 0; t2 < data.length; t2++) {
+ if (Math.abs(t2 - t) <= exclWindow) {
+ continue;
+ }
+ double norm = Math.abs(data[t] - data[t2]);
+ if ((strict && (norm < r) ) ||
+ (!strict && (norm <= r))) {
+ count++;
+ }
+ }
+ bruteForceCounts[t] = count;
+ }
+ long endTimeBruteForce = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("Found all neighbours within %.3f (mean = %.3f) by brute force in: %.3f sec\n",
+ r, MatrixUtils.mean(bruteForceCounts),
+ ((double) (endTimeBruteForce - endTimeNNs)/1000.0));
+
+ // Now validate the nearest neighbour counts:
+ for (int t = 0; t < data.length; t++) {
+ assertEquals(bruteForceCounts[t], counts[t]);
+ }
+ }
+
+ public void checkRangeFinderArrayCallsPlusMaxRadiusOut(boolean strict, int exclWindow) throws Exception {
+ int numTimeStepsInitial = 10000;
+ int duplicateSteps = 100;
+ int numTimeSteps = numTimeStepsInitial+duplicateSteps;
+ double[] dataRaw = rg.generateNormalData(numTimeStepsInitial, 0, 1);
+ double[] data2 = new double[numTimeSteps];
+ // Now duplicate some of the time steps as a test:
+ System.arraycopy(dataRaw, 0, data2, 0, numTimeStepsInitial);
+ System.arraycopy(dataRaw, 0, data2, numTimeStepsInitial, duplicateSteps);
+ // And generate a first set of data:
+ dataRaw = rg.generateNormalData(numTimeStepsInitial, 0, 1);
+ double[] data1 = new double[numTimeSteps];
+ // Now duplicate some of the time steps as a test:
+ System.arraycopy(dataRaw, 0, data1, 0, numTimeStepsInitial);
+ System.arraycopy(dataRaw, 0, data1, numTimeStepsInitial, duplicateSteps);
+
+ UnivariateNearestNeighbourSearcher searcher =
+ new UnivariateNearestNeighbourSearcher(data2);
+ boolean[] withinR = new boolean[data2.length];
+ double[] distancesForFirstVariable = new double[data1.length];
+ double[][] distancesAndIndicesWithinR = new double[data1.length][2];
+ double r = 0.2;
+ for (int t = 0; t < data2.length; t++) {
+ // Compute the distances for the first variable:
+ for (int t2 = 0; t2 < data1.length; t2++) {
+ if (Math.abs(t2 - t) <= exclWindow) {
+ // Within the exclusion window (always includes the point itself for exclWindow==0)
+ withinR[t2] = false;
+ distancesForFirstVariable[t2] = Double.POSITIVE_INFINITY;
+ continue;
+ }
+ distancesForFirstVariable[t2] = Math.abs(data1[t] - data1[t2]);
+ if ((strict && (distancesForFirstVariable[t2] < r)) ||
+ (!strict && (distancesForFirstVariable[t2] <= r))) {
+ withinR[t2] = true;
+ } else {
+ withinR[t2] = false;
+ }
+ }
+
+ // Now compute the max distances across both variables:
+ searcher.findPointsWithinRAndTakeRadiusMax(t, r, exclWindow, !strict,
+ withinR, distancesForFirstVariable, distancesAndIndicesWithinR);
+
+ // Now check all of the results -- first those said to be within r
+ for (int t2 = 0; distancesAndIndicesWithinR[t2][1] >= 0; t2++) {
+ // Check that the first variable was within r
+ assertTrue(withinR[(int) distancesAndIndicesWithinR[t2][1]]);
+ // Reset this marker
+ withinR[(int) distancesAndIndicesWithinR[t2][1]] = false;
+ // And check that the max norm was reported
+ double maxNorm = (Math.abs(data2[t] - data2[(int) distancesAndIndicesWithinR[t2][1]]));
+ if (maxNorm < distancesForFirstVariable[(int) distancesAndIndicesWithinR[t2][1]]) {
+ maxNorm = distancesForFirstVariable[(int) distancesAndIndicesWithinR[t2][1]];
+ }
+ assertTrue((strict && (maxNorm < r)) ||
+ ( !strict && (maxNorm <= r)));
+ // Check that its distance was properly returned in the array:
+ assertEquals(maxNorm, distancesAndIndicesWithinR[t2][0], 0.000001);
+ }
+
+ // Check that there were no other points that should have been returned:
+ int count2 = 0;
+ for (int t2 = 0; t2 < data2.length; t2++) {
+ if (withinR[t2]) {
+ // This won't count points that came through above
+ // since we already disabled withinR for them,
+ // nor the points within the exclusion window
+ double maxNorm = Math.abs(data2[t] - data2[t2]);
+ if ((strict && (maxNorm < r)) ||
+ (!strict && (maxNorm <= r))) {
+ count2++;
+ }
+ }
+ }
+ assertEquals(0, count2); // We set all the ones that should have been true to false already
+ }
+ }
+
+ public void checkRangeFinderArrayCallsNewData(boolean strict) throws Exception {
+ int numTimeSteps = 10000;
+ double[] data = rg.generateNormalData(numTimeSteps, 0, 1);
+
+ int newSamplesToCheck = 1000;
+ double[] samplePoints = rg.generateNormalData(newSamplesToCheck, 0, 1);
+
+ UnivariateNearestNeighbourSearcher searcher =
+ new UnivariateNearestNeighbourSearcher(data);
+ int[] counts = new int[newSamplesToCheck];
+ boolean[] withinR = new boolean[data.length];
+ int[] indicesWithinR = new int[data.length];
+ double r = 0.2;
+ long startTime = Calendar.getInstance().getTimeInMillis();
+ for (int n = 0; n < newSamplesToCheck; n++) {
+ int simpleCount = searcher.countPointsWithinR(samplePoints[n], r, !strict);
+
+ searcher.findPointsWithinR(r, samplePoints[n], !strict,
+ withinR, indicesWithinR);
+ int count = 0;
+ // Run through the list of points returned as within r and check them
+ for (int t2 = 0; indicesWithinR[t2] != -1; t2++) {
+ count++;
+ assertTrue(withinR[indicesWithinR[t2]]);
+ // Reset this marker
+ withinR[indicesWithinR[t2]] = false;
+ // Check that it's really within r
+ double maxNorm = Math.abs(samplePoints[n] - data[indicesWithinR[t2]]);
+ assertTrue((strict && (maxNorm < r)) ||
+ ( !strict && (maxNorm <= r)));
+ }
+ // Check that there were no other points marked as within r:
+ int count2 = 0;
+ for (int t2 = 0; t2 < data.length; t2++) {
+ if (withinR[t2]) {
+ count2++;
+ }
+ }
+ assertEquals(0, count2); // We set all the ones that should have been true to false already
+ assertEquals(simpleCount, count);
+ // Now reset the boolean array
+ counts[n] = count;
+ }
+ long endTimeNNs = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("Found all neighbours within %.3f (mean = %.3f) in: %.3f sec\n",
+ r, MatrixUtils.mean(counts),
+ ((double) (endTimeNNs - startTime)/1000.0));
+
+ // Now do brute force:
+ int[] bruteForceCounts = new int[newSamplesToCheck];
+ for (int n = 0; n < newSamplesToCheck; n++) {
+ int count = 0;
+ for (int t2 = 0; t2 < data.length; t2++) {
+ double norm = Math.abs(samplePoints[n] - data[t2]);
+ if ((strict && (norm < r) ) ||
+ (!strict && (norm <= r))) {
+ count++;
+ }
+ }
+ bruteForceCounts[n] = count;
+ }
+ long endTimeBruteForce = Calendar.getInstance().getTimeInMillis();
+ System.out.printf("Found all neighbours within %.3f (mean = %.3f) by brute force in: %.3f sec\n",
+ r, MatrixUtils.mean(bruteForceCounts),
+ ((double) (endTimeBruteForce - endTimeNNs)/1000.0));
+
+ // Now validate the nearest neighbour counts:
+ for (int n = 0; n < newSamplesToCheck; n++) {
+ assertEquals(bruteForceCounts[n], counts[n]);
+ }
+ }
+
public void testFindKNearestNeighbours() throws Exception {
for (int K = 1; K < 5; K++) {
double[] data = rg.generateNormalData(1000, 0, 1);
From b0d48bc96c2a86c723727c3ea5aca57420f6e5a9 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Thu, 20 Oct 2016 09:29:56 +1100
Subject: [PATCH 17/34] Unit tests for conditional MI on new observations, and
for no conditionals
---
...nditionalMutualInfoMultiVariateTester.java | 238 +++++++++++++++++-
1 file changed, 237 insertions(+), 1 deletion(-)
diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java
index 00e6c42..8f769b8 100755
--- a/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java
+++ b/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java
@@ -18,8 +18,13 @@
package infodynamics.measures.continuous.kraskov;
+import java.util.Calendar;
+
import infodynamics.utils.ArrayFileReader;
+import infodynamics.utils.EmpiricalMeasurementDistribution;
+import infodynamics.utils.MathsUtils;
import infodynamics.utils.MatrixUtils;
+import infodynamics.utils.RandomGenerator;
public class ConditionalMutualInfoMultiVariateTester
extends infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateAbstractTester {
@@ -423,7 +428,109 @@ public class ConditionalMutualInfoMultiVariateTester
1, 1,
kNNs, expectedFromTRENTOOL);
}
-
+
+ /**
+ * Unit test for conditional MI on new observations.
+ * We can test this against the calculator itself. If we send in the original data set as new observations,
+ * we can recreate the neighbour counts (plus one) by setting K to 1 larger (to account for the data point itself), and
+ * account for the change in bias.
+ *
+ * @throws Exception
+ */
+ public void testMultivariateCondMIForNewObservations() throws Exception {
+
+ ArrayFileReader afr = new ArrayFileReader("demos/data/4ColsPairedOneStepNoisyDependence-1.txt");
+ double[][] data = afr.getDouble2DMatrix();
+ // RandomGenerator rg = new RandomGenerator();
+ //double[][] data = rg.generateNormalData(50, 4, 0, 1);
+
+ // Use various Kraskov k nearest neighbours parameter
+ int[] kNNs = {4};
+
+ System.out.println("Kraskov Cond MI testing new Observations:");
+
+ for (int ki = 0; ki < kNNs.length; ki++) {
+ ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = getNewCalc(1);
+ ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalcForNew = getNewCalc(1);
+ // Let it normalise by default
+ // And no noise addition to protect the integrity of our neighbour counts under both techniques here:
+ condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0");
+ condMiCalcForNew.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0");
+ double[][] var1 = MatrixUtils.selectColumns(data, new int[] {0});
+ double[][] var2 = MatrixUtils.selectColumns(data, new int[] {1});
+ double[][] condVar = MatrixUtils.selectColumns(data, new int[] {2, 3});
+
+ // Compute MI(0;1|2,3) :
+ condMiCalc.setProperty(
+ ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_K,
+ Integer.toString(kNNs[ki]));
+ System.out.println("Main calc normalisation is " + condMiCalc.getProperty(condMiCalc.PROP_NORMALISE));
+ condMiCalc.initialise(var1[0].length, var2[0].length, condVar[0].length);
+ condMiCalc.setObservations(var1, var2, condVar);
+ double miAverage = condMiCalc.computeAverageLocalOfObservations();
+ // Now compute as new observations:
+ condMiCalcForNew.setProperty(
+ ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_K,
+ Integer.toString(kNNs[ki] + 1)); // Using K = K + 1
+ // condMiCalcForNew.setProperty(
+ // ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
+ // "1");
+ System.out.println("New obs calc normalisation is " + condMiCalcForNew.getProperty(condMiCalc.PROP_NORMALISE));
+ condMiCalcForNew.initialise(var1[0].length, var2[0].length, condVar[0].length);
+ condMiCalcForNew.setObservations(var1, var2, condVar);
+ // condMiCalc.setDebug(true);
+ //condMiCalcForNew.setDebug(true);
+ double[] newLocals = condMiCalcForNew.computeLocalUsingPreviousObservations(var1, var2, condVar);
+ //condMiCalcForNew.setDebug(false);
+ double averageFromNewObservations = MatrixUtils.mean(newLocals);
+ // We can't check this directly, so test each point individually:
+ for (int t = 0; t < data.length; t++) {
+ double[] originalNeighbourCounts = condMiCalc.partialComputeFromObservations(t, 1, false);
+ // Need to normalise the data before passing it in here -- this is
+ // what is happening inside computeLocalUsingPreviousObservations above
+ double[] newObsNeighbourCounts = condMiCalcForNew.partialComputeFromNewObservations(
+ t, 1,
+ MatrixUtils.normaliseIntoNewArray(var1),
+ MatrixUtils.normaliseIntoNewArray(var2),
+ MatrixUtils.normaliseIntoNewArray(condVar), false);
+ // Now check each return count in the array:
+ if (originalNeighbourCounts[1] != newObsNeighbourCounts[1] - 1) {
+ System.out.println("Assertion failure for t=" + t + ": expected " + originalNeighbourCounts[1] +
+ " from original, plus 1, but got " + newObsNeighbourCounts[1]);
+ System.out.print("Actual raw data was: ");
+ MatrixUtils.printArray(System.out, data[0]);
+ }
+ assertEquals(originalNeighbourCounts[1], newObsNeighbourCounts[1] - 1); // Nxz should be 1 higher
+ assertEquals(originalNeighbourCounts[2], newObsNeighbourCounts[2] - 1); // Nyz should be 1 higher
+ assertEquals(originalNeighbourCounts[3], newObsNeighbourCounts[3] - 1); // Nz should be 1 higher
+ // Now check the local value at each point using these verified counts:
+ double newLocalValue = condMiCalcForNew.digammaK -
+ MathsUtils.digamma((int) newObsNeighbourCounts[1] + 1) -
+ MathsUtils.digamma((int) newObsNeighbourCounts[2] + 1) +
+ MathsUtils.digamma((int) newObsNeighbourCounts[3] + 1);
+ /* System.out.printf("n=%d, n_xz=%d, n_yz=%d, n_z=%d, local=%.4f," +
+ " digamma(n_xz+1)=%.5f, digamma(n_yz+1)=%.5f, digamma(n_z+1)=%.5f, \n",
+ t, (int) newObsNeighbourCounts[1], (int) newObsNeighbourCounts[2],
+ (int) newObsNeighbourCounts[3], newLocalValue,
+ MathsUtils.digamma((int) newObsNeighbourCounts[1] + 1),
+ MathsUtils.digamma((int) newObsNeighbourCounts[2] + 1),
+ MathsUtils.digamma((int) newObsNeighbourCounts[3] + 1)); */
+ if (Math.abs(newLocalValue - newLocals[t]) > 0.00000001) {
+ System.out.printf("t=%d: Assertion failed: computed local was %.5f, local from nn counts was %.5f\n",
+ t, newLocals[t], newLocalValue);
+ }
+ assertEquals(newLocalValue, newLocals[t], 0.00000001);
+ }
+
+ // And check that the two methods are equivalent, after correcting for the different K biases:
+ // CAN'T EXPECT THIS TO WORK though, because different neighbour counts (by one in each case) puts
+ // different inputs to the digamma functions, and we have no way to check that for the average here;
+ // Our above checks on the counts for each point suffice though.
+ // assertEquals(miAverage, averageFromNewObservations - MathsUtils.digamma(kNNs[ki]+1) + MathsUtils.digamma(kNNs[ki]),
+ // 0.00000001);
+ }
+ }
+
/**
* Test the computed univariate TE as a conditional MI
* using various numbers of threads.
@@ -502,4 +609,133 @@ public class ConditionalMutualInfoMultiVariateTester
kNNs, expectedValue);
NUM_THREADS_TO_USE = NUM_THREADS_TO_USE_DEFAULT;
}
+
+ /**
+ * Test the conditional MI without any conditionals
+ * against the MI calculated by Kraskov's own MILCA
+ * tool on the same data.
+ *
+ * To run Kraskov's tool (http://www.klab.caltech.edu/~kraskov/MILCA/) for this
+ * data, run:
+ * ./MIxnyn 1 1 3000 0
+ *
+ * @throws Exception if file not found
+ *
+ */
+ public void testNoConditionalsIsMIforRandomVariablesFromFile() throws Exception {
+
+ // Test set 1:
+
+ ArrayFileReader afr = new ArrayFileReader("demos/data/2randomCols-1.txt");
+ double[][] data = afr.getDouble2DMatrix();
+
+ // Use various Kraskov k nearest neighbours parameter
+ int[] kNNs = {1, 2, 3, 4, 5, 6, 10, 15};
+ // Expected values from Kraskov's MILCA toolkit:
+ double[] expectedFromMILCA = {-0.05294175, -0.03944338, -0.02190217,
+ 0.00120807, -0.00924771, -0.00316402, -0.00778205, -0.00565778};
+
+ System.out.println("Kraskov comparison 1 - univariate random data 1 --");
+ // Check for algorithm 1 first (no expected results from MILCA included):
+ System.out.println("Algorithm 1:");
+ checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+ MatrixUtils.selectColumns(data, new int[] {1}),
+ kNNs, null);
+ // then check for algorithm 2 with expected results from MILCA
+ System.out.println("Algorithm 2:");
+ checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+ MatrixUtils.selectColumns(data, new int[] {1}),
+ kNNs, expectedFromMILCA);
+
+ //------------------
+ // Test set 2:
+
+ // We'll just take the first two columns from this data set
+ afr = new ArrayFileReader("demos/data/4randomCols-1.txt");
+ data = afr.getDouble2DMatrix();
+
+ // Expected values from Kraskov's MILCA toolkit:
+ double[] expectedFromMILCA_2 = {-0.04614525, -0.00861460, -0.00164540,
+ -0.01130354, -0.01339670, -0.00964035, -0.00237072, -0.00096891};
+
+ System.out.println("Kraskov comparison 2 - univariate random data 2 --");
+ // Check for algorithm 1 first (no expected results from MILCA included):
+ System.out.println("Algorithm 1:");
+ checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+ MatrixUtils.selectColumns(data, new int[] {1}),
+ kNNs, null);
+ // then check for algorithm 2 with expected results from MILCA
+ System.out.println("Algorithm 2:");
+ checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+ MatrixUtils.selectColumns(data, new int[] {1}),
+ kNNs, expectedFromMILCA_2);
+
+ }
+
+ /**
+ * Utility function to run Kraskov MI for data with known results
+ *
+ * @param var1
+ * @param var2
+ * @param kNNs array of Kraskov k nearest neighbours parameter to check
+ * @param expectedResults array of expected results for each k
+ */
+ protected void checkMIForGivenData(double[][] var1, double[][] var2,
+ int[] kNNs, double[] expectedResults) throws Exception {
+
+ // The Kraskov MILCA toolkit MIhigherdim executable
+ // uses algorithm 2 by default (this is what it means by rectangular);
+ // These are the only ones we have expected results for
+ ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = getNewCalc(expectedResults == null ? 1 : 2);
+
+ for (int kIndex = 0; kIndex < kNNs.length; kIndex++) {
+ int k = kNNs[kIndex];
+ condMiCalc.setProperty(
+ MutualInfoCalculatorMultiVariateKraskov.PROP_K,
+ Integer.toString(k));
+ condMiCalc.setProperty(
+ MutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
+ NUM_THREADS_TO_USE);
+ // No longer need to set this property as it's set by default:
+ //miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE,
+ // EuclideanUtils.NORM_MAX_NORM_STRING);
+ condMiCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests
+ condMiCalc.initialise(var1[0].length, var2[0].length, 0);
+ condMiCalc.setObservations(var1, var2, null);
+ // condMiCalc.setDebug(true);
+ double mi = condMiCalc.computeAverageLocalOfObservations();
+ // condMiCalc.setDebug(false);
+
+ // And compute also passing in an array of vectors of length 0
+ condMiCalc.initialise(var1[0].length, var2[0].length, 0);
+ condMiCalc.setObservations(var1, var2, new double[var1.length][]);
+ // condMiCalc.setDebug(true);
+ double miFromEmptyVectors = condMiCalc.computeAverageLocalOfObservations();
+
+ double expectedMi = 0.0;
+ if (expectedResults == null) {
+ // We don't have expected results from MILCA toolkit for algorithm 1.
+ // In this case, we'll generate our own expectations:
+ MutualInfoCalculatorMultiVariateKraskov1 miCalc = new MutualInfoCalculatorMultiVariateKraskov1();
+ miCalc.setProperty(
+ MutualInfoCalculatorMultiVariateKraskov.PROP_K,
+ Integer.toString(k));
+ miCalc.setProperty(
+ MutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
+ NUM_THREADS_TO_USE);
+ miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests
+ miCalc.initialise(var1[0].length, var2[0].length);
+ miCalc.setObservations(var1, var2);
+ expectedMi = miCalc.computeAverageLocalOfObservations();
+ } else {
+ expectedMi = expectedResults[kIndex];
+ }
+ System.out.printf("k=%d: Average MI %.8f (expected %.8f)\n",
+ k, mi, expectedMi);
+ // Dropping required accuracy by one order of magnitude, due
+ // to faster but slightly less accurate digamma estimator change
+ assertEquals(expectedMi, mi, 0.0000001);
+ assertEquals(expectedMi, miFromEmptyVectors, 0.0000001);
+ }
+ }
}
From f761a291e8e7a81ae24b94515dce8ad914139a32 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Thu, 20 Oct 2016 10:14:44 +1100
Subject: [PATCH 18/34] Adding Matlab/Octave example 10 for conditional TE via
Kraskov (KSG) algorithm
---
demos/octave/example10ConditionalTeKraskov.m | 76 ++++++++++++++++++++
1 file changed, 76 insertions(+)
create mode 100644 demos/octave/example10ConditionalTeKraskov.m
diff --git a/demos/octave/example10ConditionalTeKraskov.m b/demos/octave/example10ConditionalTeKraskov.m
new file mode 100644
index 0000000..dd6640c
--- /dev/null
+++ b/demos/octave/example10ConditionalTeKraskov.m
@@ -0,0 +1,76 @@
+%%
+%% Java Information Dynamics Toolkit (JIDT)
+%% Copyright (C) 2016 Joseph T. Lizier
+%%
+%% This program is free software: you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published by
+%% the Free Software Foundation, either version 3 of the License, or
+%% (at your option) any later version.
+%%
+%% This program is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty of
+%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+%% GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public License
+%% along with this program. If not, see .
+%%
+
+% = Example 10 - Conditional Transfer entropy on continuous multivariate data using Kraskov estimators =
+
+% Conditional Transfer entropy (TE) calculation on multivariate continuous-valued data using the Kraskov-estimator TE calculator.
+
+% Change location of jar to match yours:
+javaaddpath('../../infodynamics.jar')
+
+% Generate some random normalised data.
+numObservations = 100000;
+% Keep the sum of squares of these covariances below 1 to allow proper calculation of noise term.
+covarianceToSource = 0.4;
+covarianceToConds = [0.3,0.3];
+if (sumsq([covarianceToSource, covarianceToConds]) >= 1)
+ error('Sum of squares of the covariances must be < 1 here');
+end
+noiseCovar = sqrt(1 - sumsq([covarianceToSource, covarianceToConds]));
+
+% Generate the random variables
+sourceArray = randn(numObservations, 1);
+condArray = randn(numObservations, length(covarianceToConds));
+destArray = [0; covarianceToSource*sourceArray(1:numObservations-1,:) + covarianceToConds(1)*condArray(1:numObservations-1,1) + covarianceToConds(2)*condArray(1:numObservations-1,2) + noiseCovar*randn(numObservations-1, 1)];
+
+% Expected results:
+expectedConditional = -0.5 * log(noiseCovar .* noiseCovar ./ (1 - sumsq(covarianceToConds)));
+expectedPairwise = -0.5 * log(1-covarianceToSource.^2);
+
+% Create a conditional TE calculator and run it:
+teCalc=javaObject('infodynamics.measures.continuous.kraskov.ConditionalTransferEntropyCalculatorKraskov');
+teCalc.initialise(1,1, ... % Destination embedding length (Schreiber k=1) and delays
+ 1,1, ... % Source embedding length (Schreiber l=1) and delays
+ 1, ... % Source-destination delay of 1 (default)
+ octaveToJavaDoubleArray([1,1]), ... % Embedding lengths for each conditional variable
+ octaveToJavaDoubleArray([1,1]), ... % Embedding delays for each conditional variable
+ octaveToJavaDoubleArray([1,1]) ... % conditional-destination delays for each conditional variable
+);
+teCalc.setObservations(octaveToJavaDoubleArray(sourceArray), ...
+ octaveToJavaDoubleArray(destArray), ...
+ octaveToJavaDoubleMatrix(condArray));
+% Perform calculation with correlated source, but no conditioning on other sources:
+conditionalResult = teCalc.computeAverageLocalOfObservations();
+
+% Create a pairwise TE calculator and run it:
+teCalc=javaObject('infodynamics.measures.continuous.kraskov.TransferEntropyCalculatorKraskov');
+teCalc.initialise(); % Use default embeddings of 1 (e.g. Schreiber k=1 and l=1)
+teCalc.setObservations(octaveToJavaDoubleArray(sourceArray), octaveToJavaDoubleArray(destArray));
+% Perform calculation with correlated source, but no conditioning on other sources:
+pairwiseResult = teCalc.computeAverageLocalOfObservations();
+
+% Note that the calculation is a random variable (because the generated
+% data is a set of random variables) - the result will be of the order
+% of what we expect, but not exactly equal to it; in fact, there will
+% be some variance around it. It will probably be biased down here
+% due to small correlations between the supposedly uncorrelated variables.
+fprintf('From %d samples:\nTE result conditional result = %.4f nats, pairwise = %.4f nats;\nexpected around %.4f nats (conditional) and %.4f nats (pairwise) for the correlated Gaussians\n', ...
+ numObservations, conditionalResult, pairwiseResult, expectedConditional, expectedPairwise);
+
+clear teCalc
+
From bea2c77aa35dd912927f81188f5c04a733767f35 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Thu, 20 Oct 2016 10:20:20 +1100
Subject: [PATCH 19/34] Added python examples 7 (ensemble method with TE
Kraskov) and 9 (TE Kraskov with auto-embedding Ragwitz)
---
...e7EnsembleMethodTeContinuousDataKraskov.py | 83 +++++++++++++++++
.../python/example9TeKraskovAutoEmbedding.py | 93 +++++++++++++++++++
2 files changed, 176 insertions(+)
create mode 100644 demos/python/example7EnsembleMethodTeContinuousDataKraskov.py
create mode 100644 demos/python/example9TeKraskovAutoEmbedding.py
diff --git a/demos/python/example7EnsembleMethodTeContinuousDataKraskov.py b/demos/python/example7EnsembleMethodTeContinuousDataKraskov.py
new file mode 100644
index 0000000..f2f457d
--- /dev/null
+++ b/demos/python/example7EnsembleMethodTeContinuousDataKraskov.py
@@ -0,0 +1,83 @@
+##
+## Java Information Dynamics Toolkit (JIDT)
+## Copyright (C) 2012, Joseph T. Lizier
+##
+## This program is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see .
+##
+
+# = Example 7 - Ensemble method with transfer entropy on continuous data using Kraskov estimators =
+
+# Calculation of transfer entropy (TE) by supplying an ensemble of samples from multiple time series.
+# We use continuous-valued data using the Kraskov-estimator TE calculator here.
+
+from jpype import *
+import random
+import math
+
+# Change location of jar to match yours:
+jarLocation = "../../infodynamics.jar"
+# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
+startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
+
+# Generate some random normalised data.
+numObservations = 1000
+covariance=0.4
+numTrials=10
+kHistoryLength=1
+
+# Create a TE calculator and run it:
+teCalcClass = JPackage("infodynamics.measures.continuous.kraskov").TransferEntropyCalculatorKraskov
+teCalc = teCalcClass()
+teCalc.setProperty("k", "4") # Use Kraskov parameter K=4 for 4 nearest points
+teCalc.initialise(kHistoryLength) # Use target history length of kHistoryLength (Schreiber k)
+teCalc.startAddObservations()
+
+for trial in range(0,numTrials):
+ # Create a new trial, with destArray correlated to
+ # previous value of sourceArray:
+ sourceArray = [random.normalvariate(0,1) for r in range(numObservations)]
+ destArray = [0] + [sum(pair) for pair in zip([covariance*y for y in sourceArray[0:numObservations-1]], \
+ [(1-covariance)*y for y in [random.normalvariate(0,1) for r in range(numObservations-1)]] ) ]
+
+ # Add observations for this trial:
+ print("Adding samples from trial %d ..." % trial)
+ teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray))
+
+# We've finished adding trials:
+print("Finished adding trials")
+teCalc.finaliseAddObservations()
+
+# Compute the result:
+print("Computing TE ...")
+result = teCalc.computeAverageLocalOfObservations()
+# Note that the calculation is a random variable (because the generated
+# data is a set of random variables) - the result will be of the order
+# of what we expect, but not exactly equal to it; in fact, there will
+# be some variance around it (smaller than example 4 since we have more samples).
+print("TE result %.4f nats; expected to be close to %.4f nats for these correlated Gaussians " % \
+ (result, math.log(1.0/(1-math.pow(covariance,2)))))
+
+# And here's how to pull the local TEs out corresponding to each input time series.
+# Normally you would need to track how to split these up yourself -- here
+# it's easy because our input time series are all of the same length
+localTEs=teCalc.computeLocalOfPreviousObservations()
+localValuesPerTrial = int(len(localTEs)/numTrials) # Need to convert to int for indices later
+for trial in range(0,numTrials):
+ startIndex = localValuesPerTrial*trial
+ endIndex = localValuesPerTrial*(trial+1)-1
+ print("Local TEs for trial %d go from array index %d to %d" % (trial, startIndex, endIndex))
+ print(" corresponding to time points %d:%d (indexed from 0) of that trial" % (kHistoryLength, numObservations-1))
+ # Access the local TEs for this trial as:
+ localTEForThisTrial = localTEs[startIndex:endIndex]
+
diff --git a/demos/python/example9TeKraskovAutoEmbedding.py b/demos/python/example9TeKraskovAutoEmbedding.py
new file mode 100644
index 0000000..050aa06
--- /dev/null
+++ b/demos/python/example9TeKraskovAutoEmbedding.py
@@ -0,0 +1,93 @@
+##
+## Java Information Dynamics Toolkit (JIDT)
+## Copyright (C) 2012, Joseph T. Lizier
+##
+## This program is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see .
+##
+
+# = Example 9 - Transfer entropy on continuous data using Kraskov estimators with auto-embedding =
+
+# Transfer entropy (TE) calculation on continuous-valued data using the Kraskov-estimator TE calculator,
+# with automatic selection of embedding parameters
+
+
+from jpype import *
+import random
+import math
+import numpy
+import readFloatsFile
+
+# Change location of jar to match yours:
+jarLocation = "../../infodynamics.jar"
+# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
+startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
+
+# Examine the heart-breath interaction that Schreiber originally looked at:
+datafile = '../data/SFI-heartRate_breathVol_bloodOx.txt'
+data = readFloatsFile.readFloatsFile(datafile)
+# As numpy array:
+A = numpy.array(data)
+# Select data points 2350:3550, pulling out the relevant columns:
+breathRate = A[2350:3551,1];
+heartRate = A[2350:3551,0];
+
+# Create a Kraskov TE calculator:
+teCalcClass = JPackage("infodynamics.measures.continuous.kraskov").TransferEntropyCalculatorKraskov
+teCalc = teCalcClass()
+
+# Set properties for auto-embedding of both source and destination
+# using the Ragwitz criteria:
+# a. Auto-embedding method
+teCalc.setProperty(teCalcClass.PROP_AUTO_EMBED_METHOD,
+ teCalcClass.AUTO_EMBED_METHOD_RAGWITZ)
+# b. Search range for embedding dimension (k) and delay (tau)
+teCalc.setProperty(teCalcClass.PROP_K_SEARCH_MAX, "6")
+teCalc.setProperty(teCalcClass.PROP_TAU_SEARCH_MAX, "6")
+# Since we're auto-embedding, no need to supply k, l, k_tau, l_tau here:
+teCalc.initialise()
+# Compute TE from breath (column 1) to heart (column 0)
+teCalc.setObservations(breathRate, heartRate)
+teBreathToHeart = teCalc.computeAverageLocalOfObservations()
+
+# Check the auto-selected parameters and print out the result:
+optimisedK = int(teCalc.getProperty(teCalcClass.K_PROP_NAME))
+optimisedKTau = int(teCalc.getProperty(teCalcClass.K_TAU_PROP_NAME))
+optimisedL = int(teCalc.getProperty(teCalcClass.L_PROP_NAME))
+optimisedLTau = int(teCalc.getProperty(teCalcClass.L_TAU_PROP_NAME))
+print(("TE(breath->heart) was %.3f nats for (heart embedding:) k=%d," + \
+ "k_tau=%d, (breath embedding:) l=%d,l_tau=%d optimised via Ragwitz criteria") % \
+ (teBreathToHeart, optimisedK, optimisedKTau, optimisedL, optimisedLTau))
+
+# Next, embed the destination only using the Ragwitz criteria:
+teCalc.setProperty(teCalcClass.PROP_AUTO_EMBED_METHOD,
+ teCalcClass.AUTO_EMBED_METHOD_RAGWITZ_DEST_ONLY)
+teCalc.setProperty(teCalcClass.PROP_K_SEARCH_MAX, "6")
+teCalc.setProperty(teCalcClass.PROP_TAU_SEARCH_MAX, "6")
+# Since we're only auto-embedding the destination, we supply
+# source embedding here (to overwrite the auto embeddings from above):
+teCalc.setProperty(teCalcClass.L_PROP_NAME, "1")
+teCalc.setProperty(teCalcClass.L_TAU_PROP_NAME, "1")
+# Since we're auto-embedding, no need to supply k and k_tau here:
+teCalc.initialise()
+# Compute TE from breath (column 1) to heart (column 0)
+teCalc.setObservations(breathRate, heartRate)
+teBreathToHeartDestEmbedding = teCalc.computeAverageLocalOfObservations()
+
+# Check the auto-selected parameters and print out the result:
+optimisedK = int(teCalc.getProperty(teCalcClass.K_PROP_NAME))
+optimisedKTau = int(teCalc.getProperty(teCalcClass.K_TAU_PROP_NAME))
+print(("TE(breath->heart) was %.3f nats for (heart embedding:) k=%d," + \
+ "k_tau=%d, optimised via Ragwitz criteria, plus (breath embedding:) l=1,l_tau=1") % \
+ (teBreathToHeartDestEmbedding, optimisedK, optimisedKTau))
+
From fbe95040f6a153a62d625d462a5df4e43d07dda9 Mon Sep 17 00:00:00 2001
From: jlizier
Date: Fri, 21 Oct 2016 00:13:26 +1100
Subject: [PATCH 20/34] Updating PDFs of all demos before making new release
---
.../AutoAnalyser/README-AutoAnalyserDemo.pdf | Bin 237741 -> 231451 bytes
demos/java/README-SimpleJavaDemos.pdf | Bin 181776 -> 158817 bytes
demos/octave/README-OctaveMatlabDemos.pdf | Bin 151287 -> 141877 bytes
demos/python/README-PythonDemos.pdf | Bin 156017 -> 155811 bytes
demos/python/README-UseInPython.pdf | Bin 83238 -> 94314 bytes
5 files changed, 0 insertions(+), 0 deletions(-)
diff --git a/demos/AutoAnalyser/README-AutoAnalyserDemo.pdf b/demos/AutoAnalyser/README-AutoAnalyserDemo.pdf
index 7032b39714e76ac049f7a61a97c85b89b7aaaa0e..a80ab622085e448479968ce869ef4aa8a258665b 100644
GIT binary patch
literal 231451
zcmc$`1z1(v_C5?qw{(bXknV1zySuxkySuwXQc94N4(aahP60u>%Riu=sd
z&u2Zb*4%ULHAlST9doQ9gwp(il(bZg(1cZAOKza)0W<(h-8ax28~|!@0}CTNV*ta$
zkSsI+0H78$F|#wUdH8LnV`so`pl7LX0R8$kw5^?ufsQ$}Q|edF`49?Ln}aLVE3@qv
z8gfW4u3NQXVZi6vR{{-ihd%5t>3pbgqFz~cT*w_LDkv#93ez$S;LFR(k>)3{sLwer
z9rqmYOw(1?Ewr^a@1+`Gt?#U1YL0u62+ei-)&0KoiaIvq=v(uh$^B|!%C}QT^Xlr_
zOKF!1#soqjpk()j4urE~;+`rA728rz?A?erfnezKZf=*uo|8+rw_LF)o)F*eIqp^?
zrP~KBV;xVL2Xo`=;y&hF-MsMSNZ(q1wuE-ne02%rBnEf#1)EV!;ANQ#m5wCkHo|QX
z7eI^9*Nn*^ai#*naCOaD>e`me8wTiHt7&r7P~$
z%XqzArY$0Y4`6VIqnEy0LoaE8n4(jv33jLE$;;kC3Mnd-i(BfVQ
zwAfjw)VHkjx!45C1V*@wk2?#^By6x02ruOK^l-hB?uE5sH}`v^N65PX+wl(HSvuO9
zquW{$Rp!WCg<$>dk!Lxwlmn6oH-$IFu#pQ@F0J$EOKa|W1sS#F-HA+-M(JQW)htwgw{WmI<*rM8xei(pgRPWv7rX^$Am7gE2446ky$tZn?%=q+Uw7*i||*^+obzLFfyki&~Ny+?P!?y#idKDmSFe5
zlqHSV)W^1N?wMPqvi_dQ?|Q|QpC5ibi{!euxK6|sF@B9YG9t(qQR$R)S`FRN-cd5F
z^8uthe#C9eyBN0#JPyky5K9tv%@<;ncmwuxvb{#1J&{z~K#Fv8li}@2L__515^E2w
z=@MA7vk5>y6#=DN9Y`uW7!wu-$qePDza1ZC{;;75afjUI^g5;q5vC0lK9Qp)7&`JK
z!G}x*8sR7*%P7(|wzTwF_!;*N@#Hkogi@6IqQ&NkkWR
zmwpXApj;_{x2ZMC5hgU{CHVU_%J&>(F4*}+x{z-%Pxkz_X+t$oMsX?7M_ilx5CCx$
z9%0(EspxeDjq`>0rT9+o?$?oy8AF5KxC!$%RfBbrc{6{a#7990FAf3SK9_L9u=BZ`
znLN_=yq`xl1M30O-{|Kv^Jk=jhx>~BK5k?@)&{}!IS8#C7q9k_SU4ibm{GwMQ$C`f
z@U@z9=k2Etgn_*}Nn^1NpRLf+Id6Lk+oo%w2E;JTD42y8`&xu+G1Xlb5ccpu(;={n
zK#g3--f=$K`KULxIVPEa9Lx?clZ24bhTl=CZ(4CK+9cjcM3ii
z&$+KmCo2M91(x+~Fejq(+wxdm!$jY7Pze}`J_6~Ejfp$_tfBY#buF@3YwYX>*0x1r
z2je3V0X(OakUU>5o*@Cw+7ub)rYqx@4F}T4q$lWtK2dRIXA&-qNy~u88r5cI;&vUw
z^b0a*6=;ds1F7*n8L>
kz$I
zSErXK=|%uR2hR3K@{B2JRcGhf#pAiS$x$x_Xeacdl#_L1loP=5u#LZi(taLO+hB!=
z+mR>o@$+&B{+sb3k*E#`f!>PH-ku|m)L{L=d;!6^(Gkmxgl5q5g_MCi5jZeP-))pd
zdaRI0B}N96yc3vMH#Ev}QbDBz=5e`~Cmqp_bL{xKQ^<3JWCFHm=cL=l<)`mOuc>;4
zF!-HEJVp;&GqSEAQqkQS-iw2Jk6W<_uD%6*A@EXUWB*>6CTL+9>8gt~pSHMhfX6J#
z0B57SYlAc~WHCC)pb<_mBH=E6vrUkzj%~tuSJ$awKIfa@oGRdabVWKMV?|O
z+O2LQX)nVyDc-gFN3p3jVxlINRPn0jls5;{okw5#z4s&pg^?F@16GA_1@%ujd~N59
zI07J4=&$ZQ&d%!EV2|rtWPl|j^AwnG{HSm!l^6Gw_~60_FAzjgO?mIf4BbjhM(gj*
zC~fBRi#(zR$sBR>@1@$pmanB4cgeDEzN&x`MrVTB)VcckhSrnWM)f%K`NXg^y67rw
zJ#$y=K1wAfkfN|4z~K<4JRrhR;+yhzqIuXtP~3
zj>%o38aSjPm!fa3!_MS&*y?ZY*f&7LkZ>AUUw3P>(=cqPSo%iFyfr2wJxU9lGv5l@
zHC-JLtarh+$-GFHQy(rjknyJv4!^f*oIKO-UmLUYb{gejtfzH#WIjpW`Jle0+26p%
zrscGGqW#t>$l|Lly-2BP_t{BQLL%0T`5u_H7503U9CMK-+zqVz
zSDyCgbv+k4wEb=a0H0Fdg{2P2rB4fu0zAQ!SadoG6N40UA&bcOMiThVNe7NxrW#mp
z`y?pmD9Rd}L%uotbbYo9Yk$UHons`lkFhSe6Gy?rqIOVMdHjJ2$dibxcfGHEsoGR<
zkBDK@*}~-62TR*VE;c)Igfm5W%sZyJ@cuel?XjAa{C$7#lzXXS_r+b{(ukBC9&Ci9
zwHMHY3-dUcTiSf@a!cw^?M8U?tpU{y4ysGmh-
z_7#0Ewv>+HD~gq%Ron+<)4<3y?%L&k-)a-;IBix91ZkX4O+ZlJH;72
zqouwe>_+VP+XC8(BNJYN0<-D0NdDJcEH$!liJO)fLik
zXNjzxs*l-RwMeWuZC^isZ(_HEnJRNy6BPYQZe8ZJ)@~a}W2L#o2?GajEUE&vEg
zR^Q}Peb%;6-4u3!(FXByG8gmCEYiLM8|XtHw%$QJShfa3|8J}R%)sgNU|b^$m>L$X%Y6_
zEtDXsvtlt88rDs#HPob(gk~sm#o_)vC7NQQA9Bikfc3&V9K&=ML77Nrb@t*i
zwHQuiS)fJ!gaHsY4Ox~MS(3M!oBMM}N4He++g&I}DDMxtP*SvXg%^|9xWjt^hT=GF
z1(w2-Nkv-pMA?Hy!JnU(Fw+r%nQYC$7WDH_7Vssx@f8SGCPHCN%19;_eSy#9(pU15
zf`wVT4aCt;C_vH)+&&uh+klS5)W#nZ9|3T!=yI
z5(o^R!Vo76_o2MFsVTJG+4UaCPwHw`0R)%Z0H&BoOLwiyne+%!&qw(1qE$uEqvE}k
z2{Sc^LD`gY;-qiPGhcLk+BprX#yQ6kc8^>UAZ>G}b)YT?Z@-XXuq>&Zt=?|?;6fp}
z&J134XxMudQJr2=M}kXPU5cqt6e-(0er}`V=d|ZcT`?(cnwPy&ygdEI8>jx+FFtPX@@X~6ICxNe}4k_jTJa;_AVk55+9i$
zgG#fZi(rb1<^Vb51W;Gj${t6`rafu?sy(D!DT{xdFaLIOgl`-b=F%3DXU{&qXy9H1
z{^Xq@x~>Z!!vg7SBS21%A!Sf?`eW%
zCbFdtNk`dxRN?e2Rd~(8Tjgtfpj>>l7Y3f20x>7mIRXu6xV>{z0XW6s=CdFm_RP$`
z5mr;fKjkv1kqvlTDNO_vlm!LLSqx1B4yLoEe9du2#yjooMY9tkr^4Rj%X`Elh|YJk
z1>#A8lE?j)D`4$fnMAU1-}5EP{xm|T2geLf26V5A`QH6);$%2l4Hqq2_<%{x7XGYj
z+@z;wmt*xM8Z!;r66qMWOup
z?5Iclp5Uo4bFG8LZX0XqJ0pjHlh5v?e8X;O;7F(bn~pWC4r>P4321CfnqMzt_B&R?
z3^-1ltVIKqPU{9+cv8tSG#K)Wx3gN8)>^SwJ*MI@ohAXA?@5&4j$bWe)`F;Lmc|=j
znA5x);&m6ColZDJ!Hn$bneYrzK}6;qmu1okhGu{Sbb4YSAC)XEuhVr2QN(+Qwbvu
z2{_I90OYojQC1V>LwqP6t;1UeKAkL=bWMQw{6w@G4>=wI1KII!&f8#v+F@-r2
z+V}uyF|u6drqdI7?9SkPho#
zrBXqGA|t#M%H4n*6~v*#)Wq1}IIYh{hS(Rr$@RWbhmTyE#@2-o=WXcbK(L_>KH}tr
zQ5vteCc~IkVa)gMEy_|0a|#&L*0!50~{-z_jRQ>XE7%5KA{7XkBP2
z!M7N@C_GAT?r{IUOzE<-rAwm5VA3Nu;AN?C;hFB@(DLA5lA#m!i|E`Sf{(99>|GV}
z(}O$752BQ(KYoKdBVL!r%sCh9B#+>!uM85o^HJB$-4b!ZXT&nD7vX}Xe!n;|-vX8(
zU1g5KwVJW(k+`KjKo0+2*woz%ByXwftO?^~<)HK~JylONc&f=x`ynX=mHmq&-o^k|Coi)tmJiG6p#QiyPPc-Y
z4uNl)tLj;@g!=QNX~QGZ*ut6{v0}Oip`C0ehbc4A?@G5T}e*OMRz)Q~~BxS-_v5$Kz_wP1(=d
z6~AxNIr)Z!DML7TDz%mkgvcMTODPLBiiKq-UTqSmeTtMw&d?55m~@3U6O^U0BORQO
z`8rJKqfLO!G2d=iES~^8&&()l{S52UF!*QG!snmJdEJ9zn@+0x@^Jl5p!y<8t`tVK
z(RM0&EjovNonp*d0Qmy(4M0R#SWe2;oBpM!{?^Ylb@SJ26c4z}mF^2rQD0>@&J5T(fj4=~HxF4%5=Zfpx9`gpcZ%%ymaW=sXp
zp-euZ0b)mwlD8UFA0A@N)OzhieiX9pMC)ywq7cz2N6}CjoB)<(;Xd`OuMj5g2Er~n
zVtzHZk)y)+)*1mgmHk~0MN2ZtN?tE_{OfI6?GMNPk-#n9j>c_X5&kb`g9Lz{B@lIB
zqJhCG-$tMTrwX+spgo(>ez&1Q1ue|Bi3J3ob`ioX$Q#XN17g?^Q!05qip(u}9d+7I
z4t}Wd5iAz7Uc1PSH&3$@^qnxbn+FU@Wn(VTXsAf9Jf{d7AAKhj+Y-)*6O7Mcu&%^h
zW}H;t2aq8mC}7waemB|FO1?2Sp6GJpChK@c2|4Zs(_9cDUxHk|O;?~9#LqY}k|E*j
zf}O=~pi$)!dU;J70l0j>AXpn9H=%si1DWU}?0FF{M_iZo^n=YQS
z4<{)h2P*d{Mq&b8$(>-gEeQA}Tq!+f|%f^{WgbQ^-;%AT0wp}
z%*Y+
zk3Bb3yqT4_AGX`Q{i0I6^ngFvDCiCZ%TbzCj9?|YJ
zWPFdvtxD^FA|0X6HwVQlE0Gq~0`pe%GmnTdoZEJtAAiU*(i2r!-s~oqbN-f1Wv1t(
zkVbFQWmZ5S&%;VhRfyGJ2Oop^BAQp`n37Z$HHIO+=aZpHgHTLrEG}_hHrDDl26sB<
zN?G*F0MO
zK~I74gq=~YO}xKNu%Bq%T)LuGhCXgcuwPmgoAOn(#RO~$>>ELlMPrpfni?a)HIfFN
zJ&UGDXsWhEtHb9U3D((2%UjRaUZDzwW6T=(1E4+{zCmr4vuKvvhu7(RA)Rj%=xfs*
zZsSM?C$K6+;uVIwr
zlFK)1@Zc8ehEmL0F@}cstNDTGAMCygga^^ip5u3;`34=O8<`p)LiW)b!owy&Vd`_|
zBEIvd+7)RC^6f^%!7=7(n6>E$632gLsp!%P4%Cq9LK6=o*j$L8*DG
z`*=kKQ|K{i%vHNK#gJqnig;7lA$!02ZJhJkxZ**XXb5Lo_%p<*2tVUKlK~_EA`x?%
z5@Qj|;~`;4#c=Gmz;|Hb5p}IdV3s(_eJkVPx{Z9ocrp9iv}b$HCJJZF?EGkaJ2@f|
z@$aBl7|2aPFp{44kN~|rPUEK%+bK(g$x4JPzn$sko;g2YsmL;-&R##suyGY^kM$}&
zm!;HTpl`UT;iOe7@Js3!QaEfIK5}&HL~-0?1NBoQZ=}yfbc?4p#wa`vAqCzW5Dq;2
zfW(KFCs`EhA>T<-mQ5CYZg9Pdx-Jl{IvY`7_n44o<{ZorGO`4%tf@?3`?!glTfuyJ
zm{_((=lM<|?k2LJNAy|7;fw6%{dr9nEX90lHX%;MvIvx0aGRF^wz|Dyv-+!XuQ#Yr
z7XkGv3|YDmj9JI4UZQuNmihP1d=r7^itXim(1qLDMl=?HEnr<()xzjj;BCkJFdSdI
zglQaoUozbMK8S)b?*ax{1*RRsEV2Yja=6a1O>oaL@f)QzWGJXmNgN|+7FJSU2q`<>
zOxSyjK!C*;@ylve45SnKMj|EKKyB6=Z#444GD+}q8YW(=`R*uD;c|-`-jL@;)&m&Z
z#__$G5s+HIf}%5Y$&w?Oqe&5ITMs*%IqvP8oGsohJd}G^
z=?!&5YP#ZFYW0#SQr@k71)WP73*jC@R7RJ59ZkEU0ow@f&H-FA0OgG$G_{rhF7S?Fzcwj}ahO(Sy(qL#5s#vu<4p?umS0W2fsz2JmhUail^JIl8l(%)!g*l&S^Wd^j`WR8On~4
zB{8v5qU}{-1^MY+uLaSBTF9MMcUM5?j^@koh#Jn&OYE$&!U#=f%b~f(*3xsABd5;O
z^Y>-w$dfQjT7J*gWDK=7*fL|SE?b9Nb!5FE$Z<4-W;U@onitLFn@v^HvjK-TcpFxO
zc?Y=dDfj5(*Jhy51{V51Uh+SjJ+=ot_5?ise&`Bdp=D%vobY%__gnYCQ@??4;Nc&&4$Z>_fSQj7KntMO(|P#$&_VF?Guoe@(NNJaGceMyvOLZC{py#w
zbidDKq-UmMWc)Gp$Hn)V-+xj|>)1UsE?7L>VZ{hweC%idP+Ktpn0}lw1DJoDu>e?p
zoUsB}f1J_M03MF-_t4S;9*!U9bO72%k!63#KrLu#VfXNmpYj+|(ojC!`q)$P)Ijkl
zzab^<<=R_4~c_k7G{@{C$k!aqMZSk7EzgJsdxDPy4v;pRF76&Q=D1
z2eFL|9_Kxr*cd#l>}lOk0}=-MCOW*9P5`xsn`i)xtc+Cj02T&XDpr8T(*(W;Ngi&o
z1^l*nY9Sj-d#hih`~Hll(Z^dJpCbpLmbcNdu(f)8o}Tm1VZMiPegg*+Jp)-G-p9co
zW|1|pwY0aM{+{H(pMZK$iK3j_pm1ElNZxEkNrIWI~)7&-tz~-KdwzqLdVwhL6ippWdBARdfEV3ARq!D
z;6G5p&j|7G1+=t3mGJ2Jf3Jhz{q6^NpKMF}WGLPz&wYF_J%Hi6zbebw8CX5Qih=QI
zhRWYw%A1&(SQ!2C4bnO$7Iwc}OX%3ym^`B4(Z?Pj24H3cFtalK28Z80^)Yz;aImLv
z;zY{|csRZXfggzc5i0@KPw1g}>Y$jO5KWl62J_7=o{6i!>F8?2pK=+eVerY|U`z_S{&jkH9knlZz
zKg71*8sh#L5*|hSp@iR&@KZCt>fniwp1k`17bMU#(EL>>JecswijREw{|$wQsecNE
zpIr7gP@w+}
zzbB6WLTF(5tMo|E3V1mFBvbsJcs_d0kI?Wi_0P2OBs@NoQzTD~jfyK&7^{!1ug
zX8B9*!U$k^3b22qmVZFe!_+@@mtWZD_x%>b&v5)Jihhs0|A7tp_dV>7%ENc4{p(34
z@`k}<$pI(_i1;6KOom^gCoTQY9pzID`g_Aa#ib{x|F0zc*ZK9yME{Uu0s{dd{85g{
z@RKe6hF=-~U@yV=6ovj%O89QsA7=h;^Z$~WpVs{&?)`4&rt=p3(O)+6lW0Gc@Y_~`@n=ot7c?{e9{m4TO87sd%j>A1KtQ5@gyyGp|CI1IG(STPyA>K+vkar%Z@G9ybNH`hG&Z#jX_A;I$5A?&ckU(b^e!b>
z+3TILh&Us8=J8lKlbv5K$y%+0GJ7<`eS6`dVA#^!{Eg%Gv?MU&o71H-C=kba^wG@(
zHA5yvC2hpe%&XhoZnN!|S5E%t%K;hp_xD>FHOrnj30w~s*k2F3k9SWjWtV_0=!bkS
zpe%?4iy&_2gYNs4z!30Y@1Hejw&N~TU=POmM6^#Ie94Psz5^!@D8>DvrB&&mq#Q@Z
zv@(3z!Ol^kz67hj3-ql?qo^yxs2C0(b~C-*H1r#l(+3RA>BjMs19bjk%0)>qGhL73
z-hAtw_D32`w2EE@JqSllf4@D%0>9-iA`M&k<_Eg@HoQ30fIEEoYGhe
zd%&2)nAErhuqm%~P)3?$v}`*#z_i2yt8>n*aRiI4){J~2Gt2PV2$zq~*JDa1LV+`L
zSzexJP%z~qTIH{sGgcP(!2L}otSS#@Xg#HGmO>T1C*1GAmJyQIh~e_0i#Wnw4F;Wm
z8wm^7F)T@w1C7{cvOaa{
z*ooEqsDDX`?UJ@BiKV}*_dMVwNji&CkaeSF^#2iCes5!$&}
zTHOR+fYQ3w>3E+*aWv_Ym~7xG1+bKhK{-I}fKXi|^dg
z*in!doi-yj)LA&L6GcmcTF87BZ(~w9g<$=XbAK#HhFI(CVUZAfmBKOvcAIhdhbVUL
zmM`E+cud?0NY(;=@z&CCKnFyT5(QFHx$p@O_MRZ9{x#}t4rTE*#|
zE;_vG*?DKGQhZZ!V&Z)*>)h@dZt?qt2;Kr8(7OsRUi6?4{mmAaqx-
zVR)WdsU5d%Umx+JPGU^H;XsNnGU=uQ28Uw(UW1s23Jy1;*9Hgcq&nv5(hb8dv3t#V
z2;7TXE}1wF6i&KS{17zz1mo3Ifl6S`qf<9wQi?3_;J_6mFINbMD0w%A=@tJJb6R55
z7c7NFN&-Tg7Hg27D6|YO@p0E`#foDAtMHMbs2cEo5PiKA%gC9r6hv13?{c^`%)Kh8
zhfAo>O#RlE=fjMemC3nxVM~TPY`8<{oOAgVL2={aieB7=U^7O{BKOD2Ld4J`IJh(^
zPUQ?`5^^QfOS1a8(O#6gmBh5e?Hg7p%r1^34`pcB;o@(8LxrJ5+HlbAF>7ad{z8dY
z2`aXmw-s1esjELYGKhRMIZ3{q$eH8i3sg4j;#Qf?R^8FjTKlLKAIwe%hOeFWMDig)K!
zw(=tOvSeh3KGj63=&dbm+!S#|VmO`J2b%q&qJBpm6Bx
z#Mk|3W)(Vj8llh>_Xe0(E#(Qsx`WWXNup&Mph*dL7_Ti
zq3IDLHcCqwS}}-YntIy{n{Mz;zjckIP)D7}A!Xp$yY^GvEHM!8$vwEndH(5OqM=Xw
zHb7MH*3znM4GPSnEg#0d2;x>A7;)z03C49ir;4vg5HUr{u5+74k1@R}nTrWYl(B-h
zz!tfzr~2kVm$(Y8VzG1E;xVREu_|Uwvg;u`)q;0Pd@V90xl+bsaz5~_Lt$Rwm3_^o
zm7*1W@i>fxe3V*_lbG8wZ(4ur#?(%+l#^wfmefXnt{1h@IUV*K41f8srC|WUJ39s3
zG|)WNDXytFjo}^vEGp)~VkSIW4hL_Dafsb#-EYR#&bP+ji%g@6hv1?+)rl{kzhgqi
zY#QzLzr?cLGc64kUCEbNS}Zm*qyk0f$m^@1q1v_Rt;2toe#kL1n@&pS70MbH(S8*Z
zT6AsMEv}lCZp14}osAGQ)tU5$Lcl$UEv6*3e(VM(XW^hcMX<1+z2vXwW
z=Fo9PCXHOn*j%y_;%~S!iUr`VOY0dehf|?9C_37=HS;(1R6_(R8Kqp6j-ub{zwtX9
zjp5bTT*Z%Eor#49kqmOsUO)NzI`1y&)A>=RvfF2X{_=ipePQ*cMkDh1UJ$Sv2-C|d
zL+2T)uco0QK`RQXG`WTI3CM_ag4u9sbiF-tdrYYxZd5OiXzY$niJ(<4KEEVtRGZx5
zfnyhQ?6UP&GRc6NYOpqXH4za_v$*r7ghVWH7B#@!p}@T82s6_dykU@Fn@u;#v}Z_z
zXvrxNdT%;GjE%~%)jYIrCfQ!A=y=$sN%~wRJ63|A?G=WHZH5HzS)K-2yK|q2inz0;
zqswa-l=7g-iB~EchQMm2;e`gyW`P0~3?2HT$aV*f;M2WcaLAArH4t*S^=YS)!TCeY
zK^WYU4}a=nwqQ`yl+K=3ti^_(vUkJ+zF${5TZ-puqIni9B%6%1tu4MftVmHVbyE>`
zIf6SKIp0eGmSQa(kpw*=Eu2BI&ExsBM%pl%Ljl6{q{r?D43ay)Dq?~_olw&+CeAeK
zGq6Y$SlZOK)nBYccV*m8gBLZAve>&fPIz<~ncA-Snc2*l21`+Ffy0#CaXC`ME^Pf5
zPJ~7LPGj*{Adt3LFVM2yLQDc3_QeI{Qd%2|b9ZQP14bds@Tvet++GEbXqKOYwCSE;
zRyxIyy$}|jf4|cxrp8;j2$mD26+L#L+ug3YX!KqLp!cTlE>8LywclezbY&Rup*(
zx@@mWwz5PvDO4keg``XUvJbT*67B0I=^gzmx+T_fq^_pJw~#^5eKl7d;yoj})S(tr
znTn^M73ikJfyC?M-!1tIC`&+O4Bd|CuQj0Y`x|{eDhht}u`BP+(lG2}LgBlwm+c7&
z9f_;3t1#2vS{a-68e=w}Y|9EH(N!9uhh6CFSvl+YIfq*NIiS;$QtgY0>e$!lMdB~)
z3~j%Wwi^$=IAhQj%J2xfB*r09eb+r6mogHGKu3;iE}9U`tOaUbC{Bzr0Y0bY)YQDy
zdM&6|Veedviaos~0ZJ@oSACd*h`=;L5cIme%TeYaKuJez5e=6Q(JavS=;%G~2+Y8?
zJ_AM%h8R}^kP}%9p5eJ3EGERW7x`=Q)@v+TLa>{;o8=BweFk4m;>G6B_0}TL8Vi^z
zlagGYq1ZrS9vF54&5pbgcODsWX<1fcSDlJO@X$PwE6z3e`~jVyhO(MD+hCU0jXcD{
zFcW>5pBXp>eL+A87c>W~e9NM&I=b+x?OFYKj0dG)fYEE!3QeQ%R^G*vlRyH`8R2;=
zlX>An9SYXsc)3&o)g}v*UX(L33=hLgM2Q-;4*XJ++Pm&}@knUp+U`h&-U;6x2_3y@
zgx5-p&rxjHUbND?Pk2)mW5d^4U(2Wy=NKz*f1Ji?%CpC?+=~lLhx4-=>cUHsv6R4b
z>V1C0+wjHpJ_+*E1Z>?_hIlGjM~<_f!l5Tt8V|iZaCATu87RnXRQE7eM+DOVoC$;T
zocQDg5>OP4IVdo}3OB4@0PF`>f*kYl)1_QWyt#c4*`88i6Ibl?4eyHw&x2Uv+PpG-
zcHq$F{mQQVyvp0r$(zwx(e0JDd%u<>q=rA|#@yap9TQ)-KDM>b=&xJ)L1IkTGK>e2
z&A+jcYRc!x@A7*k`6@DL3E_N`l~s#>*k&9T*m0?A^+VT}B7Jqkn^Rz|Q*?TbPRf|`
zHOLljsB?Nn9B<2{DV)h(Z^eEfNaFh3)W*FlGZQtuj@Uf?AQ`d(nzqb3a9oGIs(gx{
zC_}z^FPGKY;aqh(qFnU?4Iz}09ciM8E}I(IKxfbiqs=jhV+$Y{uK
z${i^JdyEhHC#5hhj`Q|~TSErHDMYxz>>(mLz%qoQoF^|hdQ!q3vtl4QE7X`R+gByr
zT0jXk!1xlbijayr$4Mr>wWHe;QM@eS?YRX_i9jr{5|z(E12chxt5E+K$)xQgrmw$F
z`M!S8=NaHMe)JH1v=o-FySV($2_eQd%?$VU#8yv$CS=^z7&qM(KwO**^R|C@+Qd(R
z-r#m}U-DC3WGIYTrYcwsGl5LbU8k@;ACe*^LR`}e?GnRaBJWGyDolC8Pc(&7O188q
znrH~+lbTv5F7DGpQ&IW{`gi`}sp;{>8|_7Wsh4loa>8-VUQbeX<|HF{sJpFf-;I#x
zOD*^mY2(<(2*(5N*h#~{KV3&%q#+6lV?-K6zcF<;wln|K-ws%IpU^rIqf$^Nb=czr%L@kDywhEtd4teI;Lj
zl(Xu}SgL#DS{mu7!_`VECX4FzyDb{g-=EQ^erBEm+7swU~cUVnD8RnWQA*y
zo%*!e#)Nl9JgywHRjqt3XQ#4Asr>R~vWJ6TOB@!GEypm!a=}8uH(-9^)%$;ZS@2Xy`B7eEpl4$G`MQAd
zXZ`$pUn}+hrK$I!`tnf4`uUFSd;2XD9pK^k{<`mbiSn;B&He`(WSM@g1wXwVepuMQ
z-yr)`l#+dtR}eN0j%v6T(r6@h`X
zi?6KAjmsR!H1JOVGFL$n?MDG10?Yv;BZx*3zBr+b2aQvjrWy8HA`!3V%P9g0sy7d|
zu3iThIbMsD=UnV6@d*6|2_@U=bc|+sNQ5BVTBRuaE20-ql1Wnz;jR@#HT~L&3>%{
zJgX1dLEb>o8Sl5@7~`WhR(@NuP%=X_?(b6(K?x)K3^$M!6j2fq>wY{K(kDhJm(ZD?
zUMvlXjJ&shu*@*|4BS0HGm+Km@WKC{zM
zIGUWxy*rq?bSNhD>^)Ue5Fp}A%P@*1F&-IjPGY_Zp9?=|z8Kn%Tzd8JF?9i1xSSlK
zkE{m0D73B3SEj-bWu1MS>u};rDRzIjf4tm$H+E1%3BCO_
zu_yJMQcSdBppGC?xpj|XzPg^7>-OFK;`LYdr4Q=@?|W0t?Jf0th{=X
z7(+G-wd_@0^;?H=tHrftix1c4=D9)D^3fBrDWuCf2;{;XU|Wn!F0MBz(fB7KZxR;S
zRKN^+-kLDQ-n#`zYRT-8kdZMNeZ0GQky28bWOo@V;es>0yd_;|?^4`uk*cx*Fc1q-
zf9s-32wqHB1(qE3*4#!xUS5BQLiJPW&VF`aqZF>gviX9lZt7;GWPU*l8bsuxrDb^FWkjN1i3@<$ftw57)c0+r^rc&p1Q$rqdnz;c=jVysYx6Wx
zt%uJ|S5Gz`O-)T_o7_}0FXb<4Cx!XY#aYrBDeQj|xjt!JZu&yAEgleap*pkg-hb}H
zQMhDA<#JM}wkS2rTIBJj&V6;J8+$9t%aFB#9&UDWt2?I6&Kj;)PIH5V%C4`Mh+SKS
zjEbN?p40LhTvPnM(8VWqETy^13YfhkMFn|f{g-OQnb*e(pFFDjiB!p!8QYE3oV^Fn
z7uf}_`)&t%HMp;I-G$D-XmRQ2c}&o7b>~x3h8>+H6|E4D*fG>W8$m9L3Z&D@>>__R
zu8910AJLkSvR_|sUs0OQNwXrIYMKzkN4kqCU3M-B{`b8vJd!tN<#h*phhWi)NT>^R_eMUpwNP&}`Ew5ro%Ulg4r!zJe;3=7X@uE+9_pDwP;uDKMm_O0l9fNmja#)!9
zrzVDvDDbA$e7Ln6ViKlR23OAD>qYHcQABe?MIUVp7;H}{1pKewD#8>i+A@wft9WBU
z!$8vw5)}i&RKQd^y>oc2YsRc?sx-9r`%uHlmQ{2=me!+eAC9Z1>6>o^G4W@WD^sY3
zCx6lnA1cEXg|P(*HKQBxw&G7E17jnK2kEvrn$Q=b8M^zS<}>R}B2Su9c+b+|0CH&!um*Z8}+3lG&mLZaWi^McJ$sJsnzWx#RKSRkGWC
zp$QVj)XJZ3E-xr+u2xerN<3($jWs{s^xq)N)oJc8PTArXXXsE3R$)7A(LAAoPPMnGxNmk~kUUll>
z+|z0-ZI77@CX59I(7^_z`p?5_G%XKJTAgoyR(SwyaK8xIvdis2v8!XO<5}AHC`txFUWTwN_JF{qb!1
z<{d+~cp}OMOV=xX2x(Xs)5+C&0=}(avJ@V2v<8a+`ehQdf{y~M5}(w;MyJ$*`*r&4
zD8ep5;x-I-z$^sv`Xj2rZK|#tf~fgaJR)Sy-zPR6d}kU$*%CMQ=Ui?0ahiEoSEx5V
zV)fCTs#!+lUyHFnZy0~GuMrctuP(_rG4l_9Ffp4}*^WJCt4|73$>OEzpec4iuE4R`
zoHBAT;D|`)%Gb6Hx=*@fR3F!o++~i$j{GxXO^i{dMh=XG(Ga_pY0pK+RY?mwKZ-py;cDP0Y@q6(%_C7o>D
zq5Qtb;PjyYF$yN_#GxoFy+V!;Bs**_g@Jyn#MD0^d;x?E@>%?Gb
zGnotx3l6K0vUGq=m|&(_FQhGbTcq<_O#dMYTSf1yO~L)Z?b4~)DpYt|Z~AGHk9bTt
z8Wc;4PE1N3(JK(PX`;nbSmA_M9TOwTu#JT5(<)_GuGN_^YH4h$NwGto>uUQJLD{xN
zY-wd~Tkcfp(P#NhAhWmX=-Uie0-+#b`fW8=FO>@a)orq%&FF~k;|4=Un=W^
zXbd+Xm>TJqO^`kz6`4D#ND8B~5s-cIo`Ruzj}N>JhgovlcDTv7eHHJ;)vT;4e(xc7JluLkFr0l=shJ0bR-{?c@ZQDwEZKnH
zvZ;<%)(t7q$oVl|9wPIAZTX?bWxnvGwKEgovXHU~R(tL?KyQSyVT5U3GN2?`|E%Wi
zuHolojjOr3`7M_UfSOc
z%D3wq9#lXY5hRpO8M+&kMp_w=p&RM$5&;PbDQQWmpB;UN@(vDOo@@g<&@kBB)O88RiS;{T;oT?F%r`*jTiR`cM
zN{EfOptWB`v>I`I80dS4(|`NpXZM|?W970Bg7bFeu-jmuo9QrKUKe~JEfTzZ;hr;T
z^cw|(MU&6wBtN?V!V>kYOWa-xEK@J5aXX(}@)57Qzptm182yEF!%u-VJq(vVAE(rC
zo=?@WkrCMCbi5tuTQoPgDClUV!CcL#AEC(uE6FUz+T$j7~F
zE+oL_$783j%lx#KqD;52PHWPp(f2Y?+o>G`PMI34+S*Fd1FYK^pVb@XRQ4C~Sd+T8lfwBzsxlPD
zJoUXC-AcyKHFE)RXpiywGK42(Hs)>h4s=)gyt0NQc2fDIAbi5QEm1o$DNn5&2Z>>7
zu(#bg=s5UGOJe1iW&G@k6sJpT%XpIMH}1l<-fIg`B-
zl8lquok`3T^FiWMljy1(<1x%iMsW*<-aZ_Q&yKJYzt9*>CAFh?^%VW9(Ud05R|30!
zz0A_9BqNNlBko*(TLkbC~5RzfqJf%F-tu;lF_@2bf+kx--h+OJrN_
z+qrKCN;jN$=}p^ryR3I<{%UvKinb($mv@~noZVMNwbBlhP`sq9`(0l(D3rGP)vrgOx)xQXN
zdO9D;Kh;h)T>a296flYV<6cts4q7c+wjov(gK;*gO-=q58B3SygsfVdFY0h*sOz0&
zSdWoOmI}2+MZ6|PowHY@ILGa74z@;5F%4zZbjD&b7l?s+PkXKm~UmN~U22c&VKZocr4
za_qU{EdS!T{|c!=X-ZJ=fj^-rS#CGG*`dR>K!exp&*pQEjP^nwf(}MrK?eW(p&c9umGHHcO$rT3Ld5Vkrrp(?^tAw-t?lMk(es4h>AS
zP+eR3HR^~$YTSTW4>QyC#q|0)HGVDyB_T9@Hdzs6(QM!l;!RvJS%_xo{!f4iagIx1
z;O)SBnKY_l@SQC7w#h5|F4LynhfrnXjv`O;(ubPA6+Ka*RuZG&8Ec~StyD;w4bLWh
zfgV?%*4a<>ap<%sYURqvt$HYMQ>X>i$rmA90*3lKh}9oISuQi5_vwJ(vHxJRD77t?*R@*2nv}Enn4H)g!N0yL}k<4CNet@ub~#
z0`buaFUMwFBOq0VNG(p$u&0}InR({uq>`l(_GWLOr}uMMNP+0tiVI0%P0*I}oah;1
ztfuBNtbY6YJJii{jSkusy&8I@9{a5#0FjfNEME>2Q?Fj3gI?UcE!IGXeRwtTD{-Yz
zpr>)E{pLh2h^4iyk5^MW@bGMs&tY;OS)0ZtCxX$)=r`gr*<_(ovBZ1P{OK&EPlnO
z`(7C+=?#njv4gbbW!rRrawiE1wN%1R+4JP?Z*G+!ZFkHRr%QKEHS^>lm@#dIgv@3f
zn~bJ?$#J2t{>d3kzI|W0gyH0gf9?OmMxW6}*)a+=jM(GkyHtuz
z*JnJbobT`SQa(kc<$e7eGGlc50>)#wUFo0LV(D)tjzBo440giBUYC92x&-XyP34<0Ce~EZN>P7?`4`Tq7s)myt49S@%`GT;
zJ>6rQPMh`9nmW@c3{O*BScgTj`^;~>$m*N
z`F|$qWkprrzyzO0YVxdqp
zBmE$fQc?-ZUy%OU5<;yKH^W3}V9L?w?oTQHW7N!?`@-V)uu=F7of;!0a{9>E(Ouly
z*cN$i(K+jZ3@k|-@4ne5b64l~!&Tju4FVCQsPdk@!_Ph*n848jC^=I69<+xEeEtV<^7PA
z8xn{Wd&x$$++(SSv%Ie&gy%hJKy0Xm0r}7mzfafgv5_Buwq1Fausc<51(lJZr!Ldf
z-h1u~5C8`#^zPaCj)PH8Kd)
z{A(R571dgK#6z9ft@ZosJ2}au7X&?3ydp*q5G~k;q(A10^XwBUKfnz2XY%?xF}aRy=cPVbOW`KV%bFg2z~k
zE2~`GS{CuWpHe)QBaoa5_2-18cQ}j^ao|v71oc&83$jN9ap0I$K9P1EEDu~qNH1nWLr5UiC;ZmgHpm{IZKrWvz=t>#Per8Koc3dK`41nnx)TxpZwXVnI~!il6&a+V
zRb`Zp?wsHVd)bzk7JXCFY$a{LamI`YueNz%DVSv
ztYg&K>SqAiD&w8w+DTQn(4@=@o+W|pUrsXdxfcjUf$jhM@u=Cv)g!F5AyHNa
zl3v
zRSX6tz@2R3@gkZ%aC7^ahnaAq*vTA
z=5Mr2dDZYGq9rW!gvOt5Lf)-vnpO61GUzBbsF`FQ`X`Gax;VJ;woLDwJ~sKi~(6A
ze^h;4553q)yX8RodCwdfXwER~m^o91Y|ToFE4Jo4ZC*?W4?km^m5pa7TL!mFvUgKX
zxA9nhWq5$o4*_b<$D)C>!90`#*e(ATr$)A67af^@-SIwDN=PSqVL)nO$EYc_|E4=DhHeS1WbDQiqD)0}3t8Cc~YlQ$PzxP6crE^JF)k=%-
zh9sYNquojiM<^t7$gDDFdV1)B^2GD#Wc8#c^stad8^)Pb{sw(*)b@m&z4U~zeqa>m
zRWRNjGo(-N67H&78nsxa)+x;5k8>1)C0U5KK*-J9G@yt0{
zkIP`a7dan~%tQQ5_aQg_WlthOy!nJSji
z5L~ZTvvlj#LhA?)*&x6mlt9AB_LT%+kLd6Q#9en`7ivym2Yz0ideK-XtRvdWwrdtH
z4x~T(oXsAJ^Zett_D2H4U|MRDVX$9y-rrEmC*&Q6
zZbHJhXT>lHjWwz^jyzhjD6V^#$cMDvk-#c*9ovG)zkCsSA$r|E&G91&QeakK*xz1x
z)dKr<;DF^&IXzY=Jm6jbTT{2LmaQWfV(y3_aNzYBv|S*x!$ooK*+w66TNp1)``&w!
zZusJaq>5}@{mg*nV*-?Ns_6+kY0xv<@yVTMXVe`{i!2l)T(tE~=-K6}i;c0%rHrrP
z8CfiLw(ny~mAf_2L1JGy8wMU-ITLuaUWBwm;Nzml$3OepYYjxqvBsYWPp}IKEH>IT
z6Zgphj`9rV`JaxGskhxj+*@~`Pdm`em^LwY@oVEr*wuqQmeXf=M^VlZbHm;f5;;3T
zo>5NSx}um%2oWAgnb|kA9`os`F
zkkcFc;K0Jw7MNgv#M$fxjr&fQlDDwPQJD#v)P(Re!qb{(E38ie803A>ypLCME9J1aK6KsbC{Sa?ZVE<%^m
z%G`3NG4!>oUl+k#Sf@j49usIOb}Nfmw8bTocFft3C?)pLHju7p7B3Md@;&DhT)0tV
zNFlh8XN#p7zLPJ>Dr0;F=NP@w=MOrT{iK_0zPsF`)s9w(4w&|4F1D02z!&b#SG!wo
zm4#S7g=Hx4jMsBR&o6GF9{D{m$Ms633@azD26f@_t~eK_$#&tVVKwmL+G(i*gBSpO
z>@|nrg-Y=S=?AiSRm9nJEGIUZ`;_x`b9NgQkyrH}Ju4c}LJWMWBOnXuPk^5|&Yg4q
zyx!3rDMov=0DPI@-L|`u(M^=KWE=C+v_GsR6)#Pdyi#k4xvVHQ3ZPvE|`BYp>w5O^5_x
zy8c=F>!tJ6Oa_MwL6>tH()hdefPKHL(rk-`!(p_&;@2BYgJ&=aGJdz5J1;Oh44r&}
ztBzo!zhvI*``G)QxREcS@~X9;kCR_Y_NQ$41!#`aCL=0Cz7lVKinn`E@BnacY*_Qb
z@xv%b+gPWM-Yq`s_dIxNoCQVPof_;F;
zbxy;G*ukg5@aG%yIRc1*A=oER$5t{N)Q#0PK=;V|q<
zua+@?J}Xh+j>a07>z3FO%#HzFgx7Xx
zt;?7#ecxz{90_)jDcEIXkK6%OO{E?3nPjX)HvHlC)AYQ}u}+q0?EYu`)Vt7Ntw-tM
zx_sT0c-0a+&6ZVoNAxPw-Rhi(-p!Y8uQ?ud;qO5%rSFJdrmH>ai~MB0a{BW2L`_q?
zzc7C?TG#1r*VzSpmHBDYh(i2~`gDx5dv@YRaPT^M9gtR6ycABgf;=OU!@wMZYdJlGkpUaMRepkID)#JKK|d5KQ3SpT;3KFzJCD%
zxxSi*Zz9j__+m~$PkE<)e{=Hw48T`h-Vtrcxg
z0G|`9=aqL0GgJ>xwh4;)v0K|5)jgzt1R63{{{Tuc|qx4SsXuq#$&}P8O
z$fZ4g#%Q(XwX<%4S2rU1`LqpwfD+B2QUKJAWr;vNDt7&!FnCu@>z!~K{!tq=BReZ<
z5AuyQ_@dmq3%9XCHD=`v9h2*Y4Q^<`_w}nL9V0ssXdc-QD(gfu3m>wz^dLKHu4v(9
zm8gvkd}tH9mN1&md$=!)dOE^d$i#*FJe_^22QvBsgZy0YXhf}Ktx|o!t6Tf~$yhC|
z8;`$c_TVJrs^teEIyO3J%@1#t#LNusYQUsz56;19Uw7B-
zp%+OB-a0K!5@5W5oYojs3|cdPk!aJQ#N=lH14bfle&^+u@>H9sBI~`nO$mc$VjsKC}}Cr*b9VimO8kv-sI?Idk!Z
zrbSjm;Sz6Tp9VqoQ<*4(iYl$oG!Jbyp3+Uj9tB^ZlvbmMG~R7VJ_DeX4sLbqJJVsV
z#vbnwN28mOrXKHBSYGankISJh;0RP7K_uEFXw8jR>!#$2b2P+uJSotDN~FqP0aiUW
z)7EfuQfYF_QNnOR-v+a-%rmpt=&`xk{5cW(X6T-xU*ziEUAWe<`P=1Tc?x{y-S^I8_
z|CcQu0z15
zl5Es|)8f+Ub1;TqLlAsn*YHZ-7My2t7O
zCn)R*$N>&5v;WDGf6-fUy`AKX2S9jH(wiJ}v?>2e?GSTzIu-Vr5Uz~kYB@gity!=o
zpyXX>03-kPuIckGK!{^70IHeha-RC*Yb4zsp1eX+&3Hb2CGj=>{rIurt`%n$-
zKBi;vm$d_^U2Zh^6Kj8O!fniJin0=IWe)HRRYs}kS^xd45Pux%*+%Q7q1?}@O8|rtMjn6
z#f8J``R*h3??UqsJnU9BYHI3+hKAr>b|erUCO@{T5ZUw;hyK5xfIdV?yx1(>Ko=;N
zE$8@nNoX=JJ5Ed95x=WHApYe}e@AdQ74m$L=Hn9PgFoSZH}&WYdh5B1$|vnx!n*w=
zN3T{o2~kD#dL@-q*#wg-w$jxWIR&|1ZEZpfhP8NYZB`w%5Gn_5XQ!0@?2tyKB?1E0
zdMiKED4aisjEHdoH}byxFeeh?KFrGx&r=K#M$34Y))F0JZLGR}DX+ln!?!V-gL^ps
z;R;XA0neJ&dj7r9ze(xd{~&k6-lD-i_TfJb(>(AJI^XgrcjP%9XkG#TA)XuDahPsup7f2ysm
zKyy3nw-NF3gNFuc;Jc~o!F>CIxDf%<;nJ$ov_hZSv(c-b9(zX;0xK4s!~Tn+B5Ulg
zckQW_zjC-FPUly%MntW;J@z`1sLm6s`
z7i@3wcCx`Z!zi&Gl{E^*6qgxS7mSNH-WvIk>EOK+CCI)N)xUpi-a$AOAEQz6!EmrU
z!sau|%iQIlz)3en>Z#NE$w2Z)u!M))I!Jf9Q$d|>xtp~ZKN`_>F&ASy$q
ze5_iZv2%+@z7i=0Z((aVdbvXzxDSg^L$T@BQVqpgww?bS|4ZVX}eF&d;UCeFISM&7-h$p&4dW
z@i>+8l_Y6#zcp;xLiDNG=etpZe$U^$%pYmIi^7YLLr`I3)H&Egad*;VEMb^lp4-
z^Id1_skRdy;*b;fLR)ov5;KggC
ziJERXUcc9#TZ)!{%$Pfnv2mdceYx2<|60Piu2AI47q_H#cz
zRc#3#25mzy?Xy>oJ6+Sq0|>8E7D?vtFwqqXR8oT%Qqt01B5C7Wv{$5vtItX6TU$x0
zg5r$yShVgN(S+3m$l_J=b4Z{NsRZln^02b1{b6QEeql5xU8>u=L5c)3xaymrc;&Gm
zquK$i
zQ8Dp$=aeKRPP-$OWHzQJ>&wyFFxtzn@hPa#iL~;BTpD;r$*wt6uIER0K4T8U!mqq3
zc#bnenNs?#PqIh(jS4*#Nl29!uj;!cSK6kGT!o*Q8B^My<1*^gttFL5&z;er`)2x~
zMUMs(UP}-{4AQfBOMg*4-ioLimwD(GkB9?_$;ejgo(=A9^Sq5o48jQTT(s=D5&Gby
zKAxh*SPUy*j_9r%KRM#y5I&{iU-bRR{y2$_j)9HTi#T9_>T?Zah60p9PMin_N?jX5
zQzq#_^I&*Zqsdvttx6MGZ?Hz2*1MWv%!a}9WgSmAqLI`e%&F;puTIvv&|pZc?a#~r
z=2|=nQ}id=8vGcv0?eh4Gc3ecb9{=vm&R68TrnXWP*G6m&aC7B&*BO*0(=`f3J^E8
zWeNc|x-Tn;uHGc|I4`lr!KE(vQf*FNt|Kv?hm_l6|FN<4+8e8!+_PfgP2$7wui?~3
z3}tna9}a>}g<9!aM65G%b3T_66pw?iUvVcw5i;E&!sHwIQ3Cto?AZ-1id(twwbPVW
zt2(*GHu?EykGtk^=Go&_ym@daqQ0xo8Mh
zsj^vTlc1Kg!v6x{I){GPq_vSr$(W)ow5%xUKo2Z@6r3Q$(GK}V$e0LyqW!3>GBh&dsyaYujBbLnA^jO34&blBNjm9OywkEPX6xi-?4jbdwJrCrpQfjU9NxkWmxc
zKZ6(vX`QQb7TVafw7I=(rXv(>+5%Ua(AjZ3{B`oFjr`d|wWAa9s*NDsapmLiCkZ_a
zqDku&donC1n`|}pA1rUpKQHe!+{x~DSA08>=$W4}
z{vPKwIJR}2kn4hTZUq=I{_2w%7JcLY
zr91?Eu+JU_0lEAoq8Qk5N#igk%#Kf
zBxDQ5d*nU#*21j3p~QY;%cK8#KK`UK!MA+=ea}Tq^v%-v!}}5gjdCiKKRshCa5khx
z2&$B)D75A+MJMTuTGT$re4RewlILw6mEeU5`T22jiRvM`FnPoBNSrX=aZX*S
zvo#NFXp&J^__v%&TxX*M7kp@R^rElG$svUf>!o(1{aOMKM|pAaz|sYm=CwY%5MO=K
zP>!0p!>ZLq3$q9Dis{o76C&N7^vM8Zld|8Z2#YhHy=VXMvbKLu^Q|?d1$AA9`x$pr
zuX~2E-nXAed(5omzdsE3o?t~t_~XA@yv+!|Ay77SWX=~P&`KzqMjzIcBhY$=|DSve
z9HT2Mv48egEpx8Q$JbUjN0tT4T`iWB22T^$q09D2YLg@uznB%K(Q)jE0W@b+?DML1
zXlm+ty%bl9A2f=>4A-+&+iFdOA)BGQs!}p=Af*JVVrMORJfwS{BcXY@08u3=9-o?N
z2+EjEw)HX8)1xR7G$ojedP6C6SrnTfpHvjp%v3~5zLej_m`YdJ$e#Xh=T#c%RD4*Lv;M+=iqx#S=q@F|dzYW%5
z*)YE>FKhZ#{mfEiX!GfH1zs--{DJSHA2dP0WT)o3WPHtT?#}ULwtK@R6KRUrMLy_nF8zv;_fSV=o@FlmmF|cZY|4WpI5g+PzVI2sI7HI)hoxGqv8Q(N7ajcg^~#Oz!C%t?@5_
z@sHgAK^TaEtU6XU$PtBZzlc}W$mkI-$wSccArNVxz4%aDS2z5-u^k42PRrfP6{
z7(&V-M4wt>nG$3Ft9FY%BS8><44bD+#!HZdhFvoxpZe%XC=Il(=FD)r
zMIrvULqJ{_8z}Pgl4!2X4?ubkBA42K|Ni|5Gjx1Q!Et*S9gw9V;rl%E%s#b)!}ggU
zKKS)EZ9W!CPBpRE4Gf`5n}Q`7zcSDc2nRW(s~9d-`1;sxPgys|dz})^2(0Yy%c9N=
z^c??UJ>CKtR#$-1;F%&;H=T69{Y$Na8Bom5HW9J_41@KN|3Cd4#-JjhKX_
zB@!K^r%FUd7QS^ZE|XTEv;e`!{6~5PNUl
zai+UN#oZsc+S0TcDzbsnM<>1tCzPDIo^stfYX
z^M9eg=tUp8xR>e=Rbey2cnyTdl2bRdBHfgvfXc3$X({2JxC#Yh6N4h3{c>{4UnIA{
z1i(cnN&tn-0p85b0>ohP>}ts!*+@Xk8VMvrw#dO`4Tt`;shxWS&(P{b*&)B7e)J0c
z*>v|WKoN_JmAsWGU}VqGK{9xY?uG6)HpQyNZv*Z(_jgad6Z~f>5fM=&0!gNouoMaO
zHK-c%jjpb4#!?F2Jirdb_spg%&c6#EMcy5nwdrDEgJirrRb*~~+?Bt7aJb8%-=8Pm
zo+|nHQONxK2n{AI9>*wn0Fq%|bf@`zSfurL%iE0sAK)awf+(C;%A>FdjAjr>HsHDg
z7H#PrE)gVT{wGFwb$3*Jhx{t?BQl&S%^N#GRB;ha59Dhkep_I9#fL|_mxTl{<@^iA
z5E6lOMkWXiWGzW&Cyar=r-lTgq2F+0xMwdBR95n4MRYRUGIft*ka1G4izV@9h=f1e
z70yxmOgD=u#`Y(y6ZkrycemnBfq{W|q^uy2gt4ouMx(w;d!Nfas5ygWW$Y5D{1V*u
z?!CU9yW^X7J1S=mu$MasH#A^wlzGpv02~&rmF;lu0b*`0ZQT_KM1gXc_ERz38{HKR
z15mo3oIh1n+<>ZT(e1%Nlx)+*QWKoJvD_#};YUtiK&VB8=zfIp5~*~rr9?#2?9FEV
z>vBa{pNIKW{7t?xUgNoPJ+MRY^%?>AeC7iP6hL`*Q-$uD!+!thD~PH}vAuNw2dmcA
z=`FRXO~vV*G~32|L&`0Ns=1BRyVZrR2OJwh7ki#~d)x*lN2x_d;NpuMtM!NhscFK~
zcY(11s*1MYp|J4qZlQD+7nkVX;BZ0^$mq#9a#5bQ@cpg%S-1UG#ec4;>9-s*svtvm
zWB!JV=MXeG>LTXiX#$x^ZOurECte!ANfo>%YoI#o`zPb{n5WYrT&|CaKw^2FO-mx9
zX*ei<)c9UKW9&DNE;OM<$$W<@3P`;DuVD*o&|Eu0&Q)VaPrN&J6CRPgo-DA?X
zRzhpCsMBTt^?(H4@wypn+6E~sP^3twPsn5b-`TrfZh(LTEq?=;DJ#{Ny
zpHnC7%RQZ`w%XXQ5Rts2VUEM7S?|j1HD1w^W8JIWUcrK+bDCMZU%g{)8Obo)(Q);<
zZ;EHpoONrF3FUkkM>SV)G+k%22Y}E5JC>p@qL<1|2k=OlfNnV4v;q1fjriIK)1bfC
z1;C*UmB*TWA@AF;LuCtc%(PE?k&~YHpWOs1xeM(l?L{=zWeu5+QQ%2*!Q})abD;;8Q=Q2}=G-0;=TscUOH(9Go8h
z_dnL#ra&MHzh3C6>h6Nx-$kn(0~C){OXeY;*f6P0^#Yb<2T_-RhL1-asPVF?BkpH_
zR{^x+SgzyopAAEB?-7mwkm*##E_-#Ke2{RtJ9Y$f%smcvDFlIHLV8__F>!vV()&F9
zlWq$C^G7~kKLU`iD!@y7Z2x7mgx9U$zbUh!->^>488JqYhNmFVGWF5R-Q{ly>oO}q
zy2+PH#ZF%PrtxU%Sxp#8$RhL!5AUdgk}0@r*t$L|9lEi%7w_;IOO;eLuCj_`V#Ly_+c64s{7O9%EYw9J3c$F}gD|67JfA;Hf!p=7-SC0L^h%
zTlXgLC7@6IH!s++-sIQ2{&_$qiQ#+G?rH*pHORGrdP56L-bXDeD+wCqKs>=2E+GwM
zP~V0}P_L>CE*;PasJ;grRHW^>jzwld2Id#fJKc#0<3DItQt}y)2SW#tnQIUI>Jy@z
zJNcOhLE)-sL`n_cje~67RdlJ8T6!#rzXAh|V<=2aL{^po`civ5^_5w@`EQ3!#y|E+
zYqt5@2v@t=m1_S`*(d|6@9=py(3{v@z#RW~Yl~#OW(LXzL66Bk`Hz>cN2!Y*19`eq
z@-w;d7gj4u&o)*emzS4WO5MW3!olHBfF`Bb*;$>9htyq>>qN##@CRZo3s(ib!|Gjg
z*`fh;1apI0999C;BHL#bl&tX
z9#IrjydW-LmtHvjFUW{b>@xf7%gA>9%IRfbYDMamQ3I!mT>|VJ=z@{V6`2y5
ze@^P|XzmiQuKH$#qQSszH@4*9`GCFk%3J5eeA7kEPqOqdKwzfxS=~yXf&25m(=FFS
zPKf%+y}X}pSKbN)WjLmN;Rq`shjyGS7bAPan_5mk<<)k|eV{H2KK(@8-Z5IjRcNjP
zw2-vFxKQD{;qmc!Ly_O1WhF>uQJ{xPU8t*LVYP?8p(=+TkB}3NKXzO9Rb0NO$Fz+@
zJgJ82PZtZz-80dr<5cnj<|~lpTz5uhDp{Ishu7Sa{?ClfQ#-91LaR
z!VMN`>t@AV#(SJxU{SH(sLPNo^fV)^-)p9t=3Gd-+Hsq;Bw|xVz2{o(jA$_j6WAjW
zs3)dKWpMQT$a@D*nGhu!8S?GV$2=
zJE(x`?LP!U6Y!UX3OECF?z`EdKPS^i(%d=B0115=z7`$u(uK0{sfCbeM_p)xdW&b)
z%o3nSJW=LXlu-|R%2OOQs`Rids?q=5$cP@DU!7$myyRsb+&gICO*R8N7iFa$J`V3u
z67o*_v8L}Zyqj3xIf>3gUN^FOdM+iL@+BFHz$>rrgCbvmDHh
z%Y9|ghOSREa9WB5H`9RhMd^2HwMLOs%aX&X*0hpe7h9qL)>Wc^cqaO;N(3-#c#0YW
zN3+woK+@&Dqfl_%JVULbi`!@T)}$EdNN`W8>K_;k5vuL}i>Y%YzFuM%GL)JbnZ$kO
zPlW0d`CJj6SUeaMEXdyytz+R>pLSSTuQFcYfyu!t30RASC)tH=^}J`3lq!tP^5*+|q%R4p`FOVKiRT-ObIu$HCQ2O`?PB4?toq%FVCv
z?i=RlkXI4U3!yLV=ml`&j_tqzXuXsO>jB{3MlawXvv9rwab$2fJIcM4-Cw*f*nfSa
z65|Y2SJiTip|h6kR{*~}@Kgl!mWn7a*ezN95ULIY89=dCpD02uTNP?LnkpN2aLi=~
zX1CF_tBW?s#Q?r)NTwH+p6&wnFXR0S_N$&O$wCOJw{ppNDI-=8QtdL0_H6Go078pYMRUa
zNyOZ;{Gl6LThSsH?hUgN9ghKM>^HC$EuwdNdP*nIp8Ic#@Q(cSs@uqbioU_DSz2!d
zG7#DVg`;#=_}3v*9e_pN2DdFsXID>2Gyf^{0G$U0I-|dbhojet0E8MWGgDxCAL}5#
z10EMT=Tuh1kEUUf(SS_OBJi1H}INGe3<&ctLQu$_pTk4ESE(GFekN2OAqu
zxCMa@0slV~)R*P8fB0tr}HdT`5CYWcML*4rPs+wqv9FRfZB!&>Mtld^C
z06kCw&{mF|r^}9ekUb>2*LnqlBw2e#UUhPHO)V^{BT$8iy?+o;-ybu`zdjzGyW6p8
z@^P^-hECx3i8oh;S|>)m&1V!c}dW(|(i)2zt%zhDd3R@8m$eIdMHa5VZqGXYel!h>Z)erVKD*k$|4
zQ{!jsqB~T>4XD*r&`)lTf=M$vcY0eETB=?~YOmmBz8RMtf_z^1m26K)!JuXuKojm2
zfCI?$pi4Cs>sn^FqdsAPI=!72SQhB-)LIH;CR#<0*X*FT8?M!(jgtB<$)>D2FnZ*+
z{<~{xx;NUdo7{n5_cNwaMWLe}khN#*<+66fHTx
z@%={g>zR;uk={jtog1A(`fUMvmdGIE=zi;fWHmsegRz0z{xy7e*E+bBPFBIucaf$B
ziq^tLSVPtreW&r908nUgb(Xkmb#u30Xhy#I`nZD4bi7=#hIwXQ^!xX`dq*
zA4(ZxDCgI~J$1CaaC*NmnknsR8Dw#wZ`+gSf?luFaM*OhqHnyzz_
z(&4Uw4I2A@^w2ZddOWtF;+cJZc;#-r%OO)oaQE_dHkx$`)hn}AR2w0hlyc!loFKbb
zef?|c%*~yQ=TAHRJ=Z1lTE>h~_>tgL$z(os6Dh!*=^Ylm-m&pY!4wGJdyDm-O8_A0TjmG7
zZ*&J;G{bz2{O7(evp#wyM=<@zGju+(My8<*zF8myTH`+NUfF{G$`iYqAkvcNT3z??
z#~e+-h*rl+px#?;)M<1X9-f6pAWhWgV1510@hk``>*x!x#PR;=4*`d0T#k13mq>s#
zv9kZ7K%OOIFGZLLx%a)x(SgI=kQrSpbG~+7N$X(W7zxT#3skqt>g-3um&QFlRw4&mzWyecS8;l{>E+Uq^2Fg)Pu~xP@DSsv=j8TQZH3YLp3FfX
z2npT;X}gyiB^UEBy~v?Hv<>>d!-tAuiO6&TRE|qyTvC9cu6~SRczQ^l>i;IoiaHN
z;g`!?gumwA_0y^Sk41c~SozaDMD-1mQ1zWle8ilAQdfZ1}P4MSy)=o`++3cxfHuMG3ZHepL0C6S0gaS
zm8c*i8!RPA7`1;$n;#`sl{DBA$I1NUfA5iM-L|%2EqJC+yx~68k%qQA!
zD>XH!Pod7@@I32dUk6KNjc}uil4+)#EwvQEc-P%ns|3%1R+&R%AlS9U(ADvm*OKqs
zZ#3FJTnIhh24Be3b?ljoT2c+4ic|4^lkmn|6m6S>9~xhMw{blG0lqTrDUNdy`1d>#
z-5~7ok3R^-_79GvNaJz@<7el5;ZbkvdG3;StFqrn5O@IhQE4(M8kQ(9i!RTvz2)3R
z=zjX|QJpHWyvB!MIi^H3KoIcLck$4T-xdSQBVB(0{PWL)f5H1iYaNv+`i7Z%abj^x
zsTp)WL9Z;|cz36O1wP0aRq~?+mit=PXb&C4=Jqc;Hq;ICjREbIq5i*^wm-(_Un&D!
z`Yx2FS6zSsi(&l7Y6HDiGOzKjT3JDQB)BtY3lKa~3ZU9XdJoJ&UmyL4jRWv84U#qO
zAQFK1h5!o&01#Q<|1osHav30l{Qt+dXGK6M5kMQ+cO&HrL&0gejE9dvzFE|;>29DO
z!WrCWuh{y_ZKIzM0?
zo2XrrVmwClM_c0mGRe$y6RZ<0$uTT0uVxT_r8I%W4O9DLUpWdC6_T_ufDt?bW&*S&
zQN_2P@};%B1Nrj|M0dYb;rorLX?nPwG}bMN!Z-;lrak`@{k*Q#UG!i6pN8BY^=-$e
zneL;JYZ6ZR&@a7jxx>9I;pJ|Y%$`rWSZ$Xbf1xMLyAb^kx8#x8>)hJ4d8cQ^=Gp-O
z+sjUBO@8Vzi9_7w*sr21v-)SURYbO~r-%Mas-Djl!5}-)Q9jFJQX}H9*11
z372&NcinM&|5`z;tQ`vV>L0}jtwv`xHQ#%K(4t5*pP>Jg9q#rh+wDToH7!@D;3rw;
zk7@cqLSEd7dH?AFL6QGE`d&|3nhh+G4XdJ&F0mm)7IWzex|UkE-rpuHcS+)y4&Mq)
zLq+qU=Ga{9odkmtCWe11j~e-3^v|h2-DBOp+d^14vCuaGoq4&EcQ+)Q4HyWa%0hnk
zWxXS&%`ctQr^wiV0axVm{V&P@tvWXQrw7+>!ii4kS<%!ImZL|#*2@Am2eYpUO29W7
ze0o%;76nX;nh%w>z1iDaLoa&T@Hv2W1*c+{5YqDRDcJUw$0Wr1d7T^Ie){JvA^$(L
zePvu!UHk4Bgs3zE0wOI8-Hj3w5>lfwbT<;xA>HlJp-3u2J9J1l(hbtx-Dd;p-5Lw^kvyXMR`F4AW$Zi~vBJW5#&ia&ICcM~e^1v^#g(
z@1uUXiR@d-1DpkU9oYQtmf$ht)%SPotZV1aa{XNn?|NKRnz?*-r#$w28r^>3KzL%=
zG)MEOx~q35uhhs+M!ksoOG#3H@uTko+d&2=L!otnJ`GN)*yeioU+rr)Z7%wqG~O_%
zFc-UAl5`eSas)5~@heM*=4831cJFqkw*t{-GnHi3xeN1*b*aw(xm8|qg*k0!j
zlk@bm4)tzd+HRg5MAABuTTI0D9H@5Zm*wNGoC8sl>-iLFQ$g#A+ILfQC!bFR>8Bv4
zzw$|-7I~pqd=%;wxE7>D#XOZeW5-Pa!5QSZU$-L_REH6`q1kPJb3FPal{4pbty!c7^3B|o|`CjTZ8T|
zor($|Ygo{4Dp~K-MhdD9%+?0r<>O{AMsw5D%72)DQ~vrID{Wk5t$ay4&|+kIy8Ova
z;P8hYI*;(MBXw*E>X!CI(pXB7`6F&?n-c4V&cyhZ%7gG&J}!zHh$fKWKi+vU+WJ&c
zCdayK=~)Z4axfeh>u8o}t^~KjoNfvC%;lSxRHXuoi)hZK${_Lx+XEG%Q4u+t7E(0O
z$9tv?_kRCM+ekq5St3W#m-Q1AE1543p+@hzeJGOFr%u@D6f((=S@fMya9zhYhmd6)mkMea+bQKp
z%a$Fnli487`c0?vR`|(7Cb@=;$LEFbpFo1s{m`&O#fkEJXA}%_~IR*84tBUz4}W7@Oou1ewh#TBjj6=6p0zF
zv={VDNE&K_Yo^N&lU-2rVWtZG
zL#>B%a_Eq3OCHmZn+J{f6P4Vb!fHP(RF1C%BybvZmyK^YrR(ZUz^TRAXh=Zs(Xkhx
z-?PC>D2SVvJ)q3%9{?CRa2xiX@IAcd5|_lu1E?^ur&OBAnqc5h(2v`uVquPNERR*p
zK~Ky6elA+*cMqhy59sw@UtjssKh2b?u=>eO5bIsK8-L#C$Bi)WzX5u_tzGL3j_nBLaR0M7LG17h3``Do
zOp=&+F%fZ8P$8C`1(HTj(C__2Zw9_4!%^BPkb%N;-$dI@zrbg;!{cqEEJtL5;Uek{
zP0qh={CWU#y0(%1N%`HR&`+-HJ7H@f%}P^w6a_N{pE0}%*geQ->kV23(f2>4Ug+MZQLtIf
z2Go-k>K1O0I@z#9{X+_Ljr_QPwm1aqm}>k&reJr*SYG3Sflq65KcV`Xjr?$!G`o!#
zh?(^2owLBe@!;5Qc_w9U!{R@-w2F%|cvxqFk-mrxibiL5~^rjc0t|sjuzMg6HdghpIFffaK
zn!7Sh$yuNx?!|2AaF^Vybcz59hydfty}9q3r<(#AiI{gs_tyMVi6M||Ph->K1enNx
z4aFN2LA=gd$C2#F*riCt_i98gE-upQugw)1%vH3`cBJt`O9_aHl@!efitwhU=7jLO
z5nVobL<2+T(&A*{hYVVx-Xx#ypZJF~t=a_{R(+UI0SE&aqL;5;1Tgb@!QmgP}r
z-Uh;&ckj>iNyuxV?>}6{^$7?GpYb$n^|^Ru(~}RNpbL@dOJLFV3=T#wv2jvxSxhL=
z^A4ebYRErgGvCU?ll6>ug|>tkIQ3P(c~ai?3>)@22Cga!aIVIWo(mJEqe9LD#89AE
z*pKdhk0b->QeExC>!o^g8}f!ji`g}3-n>DNu&-2Cr8AZ|r&G8svPvbdxbdbHlBWs8
zgh0XvLjmiM(WijN++VZN4Wev{m@KuX=`{!wc64&G%*vdd?}(>O%%NHdmg5Q#GoXdG
ztTF|F-nV;Pg$ZEW#hK}unzsLVB;F{|11Ts6PXLc|nierGM(N{%66u!rGDJE#!9m^`
z|27a~HL+gaxeT`%-(@Z$HSbzB?u%|x5}-{q6$uwRizwD0_59zTP)F}&DA88{97>IO)_
zIeA71*lg27-MJNqI)XjKS9(Tro*>6Wm6+8>k6@H}uoAV1Vpwi;_wfvg=Rvw>?}2Yw
zm2u{~9fiZamt+E?eH`=2TcLIw$_i7TG-h;=)1PaBf3MwU1Z|~}ax$yg;|(hni|%qe
zbGsv}3nyVs4Zinx2=^0S4pXtOS2=xOPU8DYdW3BEnv9pu`CrfdNa)gbhpWp4_dB
zv+a&6W@oQYc2h03>|o6hFbH>A=@cQ!YL$~mf|DW}R^>qiSR4XKnEQb?qUm8_2wN5K;hwNI%S
zJyr~*Js32UMfZCw@jB+o4CE&wIP|*C*M2(qM|5xTET8uUNr}slpVrItcguW(rxqmy
zOTCmiW_VSl{(QYS=U(Qd7^{kQOOdnNc|!nA)a=5D6XBIsaE
z3v(enQz&*Wjg3mWlh{d!FL3-F7PH|rV%kuXTm>YMLes|%5xG)@!dO}Ei3Mc29EkZ2
zT38Ka+EyO;#^dzC_PrliVxfkvjz@E2`q
zLOe@L(XyvQyCQDA)r?gtCx}k+M)DBf+8($hpQsnozVv(fo$1<3{ke52FO6=Np~X
z)s;zB=v`_19pA-@^e$Sy0-q)2cm1NpFt>3`T@{KO{bBamq;2n)URK0*-G={!v-{!_
zsWy3bicnHgby-aJDeVQ;_;}~r)2RKb-sq}ap&~b0xc!t{ky`O{7wN|INawhy0z5Wd
z1(s*Zh;4Orv!#JJDh|Z+I|nev*u7H68*
zQ4$cX8}rl2W1{A`GH!@T?@hO5g+mB`H{8G
zxp>=(aZAR1l;Yv2MX4_6iR}m6`4M8VP^F5C{6#yVElVDOfchEg^AGsZDs8ZwS9Z}%
zZ$9W~C>1y!gwLdSy4vS`Df(g5rgpq&T{s2oE;U*3MCiV2u89
zN9`Fu&W_#8G~~lfpu3BW)y`cT_M&b9AEY?9*m)MH1kdVq5_yhWBN7tb%IeEIu+~gf
zm&jGy&rEq3&fJrMHwmqB;~L!Gmdw)dj87sT9lgb1F}##TG_BTI
z@w9D0C@#v)EMu*j(g++rXR3q6P|2Tcy|e<^@(4X~rmE(k6Ef#WF4sK{n0#zGZQLob
zRaeB98nJL$Dl>RH&C!6-e%zw8@3d?aOt+ncBbc
zp}D5hS7oLiR-|gE!~NF?HN?%BA=1Q
zPGtPz%h7{+JPCWM(8ZEdjF*EUI#!Pcb-R#gvZ(4j5#
zd3ZYR@YIJxQe%wBBT}fO@>n$9RoGTdK~T%cav_#Y`f&N(9n_0my^V7cbCg)M15S8hiS}D^mB&b5i_pq
zCT+om>;9)jpuBg`sap7fh*eG6bm)KatutZ!PHw_*v^bH_>3hABHMdDehhSCn83=Bq
zGC$MSQWMFFy^DA{hddL^CEdUJ9@mm2F?iRj_vV-ADV`P#_SNB;P-%3&8xZX^~fLJp>PyIW6hX#juo
zxq86|#mOmM&3Nv(s>kE_cUckf;FrE_=vvtgeMJb#VT8;(#c{%5s{v0B=C2)M1Bug(1WIylU%ZZUl+8yWnj@gA-i9H_9A;3o&YVzW}o9g)+ez78vM#20~F
zl~ylzXP@6rS>{@R(d0O*Pnfb4kNXjmFrx{iK&|AlIS4J
zD&3^47N7Q3#B;_)GYDuM5K7R%*O2$ufXaHr``2$#ERf^p96IEaV2IWUO^>22uyar)
z*(4ZjQRKKBTb#pU&j*8Z`f{qH88Fu+HNTU~bT-8@ayYC?ULAROv}@$5kl@LA!7zy^
zh7xawl3R7y<*TSbUWU_93rz~X2&p?sEZSEHn#mL@(me%07PIg!`Ry|0il013-N%;V
z=d;f10i-G0-!i6yl+^A0Cn>gSv(PB5!JS`byXUZj+Tf}zQX>Iz+Y8Z!alJ~B+!uAQ
z_B!K~x_sC&!nsx8&{|s>SOHVt8@Urcw>Kh3$yN$4vx%`+B8O}CDNI+suCYW1yB|Qz
zCn{p2WJS?RBXOmA$N|%a5M>JBK-0t`~7kNcFGlrWrY=>7A*YYwai-kJT
zS~3mm>c~Y}t4E~a5BKdJ#Ho(XN1vzdw^+AqM(sx)?dl%#oMJ@R3&d4)5wlythVLcEq-xT;6X7&N?nNKw>H!No$%B$y_e
zyBZKHiBtwZDkm;{(RN^vO#87v_acTMHWK&ZZok!>kfmq$bV^n>as8hE*pv#R1-X+n
z7#S*wJQPrMZ+CZ>G8sh>XjzRIQ;-b6$Lf9Ahuw9M_V?(}H+5wT*tk3UH%w()knx!3ey;xepm^@iV4z%7@T{&Nn>gTFTF#PI=?pR^f9hcPg
ziDrS+#-L-|yt@PjZxTAqsXi)ia~qc8VS1iC=Sq=`5RzND5STP8v8gd~(!VeeN;KcF
zI9^tmO$0Wx^?WsccVpuM6782z(CL9d*&S>1u2InVT@%}N2+MVy3k%S+v;ZvuG+6_6JDy-gBH_jE>ZPF{67V8>&)2vn;fnPD!|Ky-F
z^R$FwidWFCw5Glyi-!V6$Q`{YfX!t5{dg>#s_kaY^Gj`ZiFg(~fORsWXA60$!aVW&i>X2!Lls0eF
z>Mlg^6hA;T+PjVxL_
z9_!Ukysswbgbq&{9ypeUIqkm_cmaX9QM$Itv{kph!B;q)p6j|ens021njf)#05csw
zGWTuWHY(k2(>M^U1S&j)q4JT&!h~0^)zs#EQxvcM!AKE&U20n0s--cJdZA)gP;u)x
zHuQoo88Xn^UK~O2L_TlQ{BuHR?Meob#%P|EWe5Y8gut0$17LGw{o{(GTjCw(L
z;}FUWbmU9RzPsx5^tO^UUW3{@-8MC)-J!VQBS?K{u7t+l?=tdeZ0;rw&MWKZbZl>B
zmz`2yry2v5m8bO1cu9S$&xmhRBA
z81e26H$O7u>0eGD#sq-VDXNufafVmR{_31&(}mWqH#&<|LZ0v~ta{3LT1`=@tp^d7eqjPq4*0Ml{0#85!Az
z_{D*zL2LmYWR*80?a;$;C;$M2>lv&(I*Zud}HvYKRO;vDQHDQ6Aw-y-1mH>CVGLzeG81uM-8aA$=Xn
zzeI@cMrFUT*$dSK7uG<6#FD8Al9XK(dlr+P(NPr)}W#;8=y;xp)y3)5lkdH?qb|1)AVn6i2ksP`f54_lHzr@apiY%aj
z*ls~tuF$o+0F7Z&2T|}O8_XbplPYRxB!_%n6rA+d0|#h04@r^V!E%moS&CKL;x;Oj
z<9x^~1d4fhW$m{QRRPK%jR+CJdMrYj+{L^)`emsR<>d7{V&Q`p{T~n;;`VH0P;azSWD!m{t3BhEnA&%VhMPTvik5od_xY1D&5fZ9iD-B8`
zTo>XNw0IgAx%1dZI{P~~Jdx%tH$Yi8uK@1^WPAXeU-E&c5|9i$_MN`Bs%;arJIE3x
zbQ8daQ@bI&X#VvUcPQvqQO$*(%YRl`XPc7%A@yrgd1sVgK9ME$4I5cib>lz61o*`Nt?-v08S&vA_!C`
znP7obn3GC+i5z+Z1Vfn?IIF4V&U$}!(`#Qr%AI!r+u^Qz2{6k11{nXjEPsF5H)pY=
zxZtUoSW|i27Uk3hCcO9pRc?RvaOpj?(uK)rE
z006So-GqkU{`H2u-QC@1fv=#<5r?m?_Plv~Y2-SOn4gpmyk*w$3hMQdnjlc{Z~XVy
z4(;(Ei#CY~569vV?~@l@_T8k30Pev;?{`NwBI7e(G3_OJ_5pWiTHk>H_MIW6)_<9NZFs|*e1Y9$NJOZA;_xqUj
zA)+*&HOUCttwfZ0G{#^EwpQ%n+g7Do0ju+lE#H40Ctv#MfK_MHbFV%#boD}wR{5Sk
z0z>)bQgyKT@(mEL1B+G_tbM*It_H$=pAn$>dhSCd@p;v0>HUwR2|~XnEeISj3G18b
zsqQsDj*{cJHE1jh8gq$=CztS%8wjg>gV{NhrO%2z^v!)g$a|o1{|AnmB-bji2d@dW
zKl^W~HZz)F#lCB^xf#`(ILcjuCc))kxS=aTfk$5hi){$avjz3!HC5M#egfpA5VBGZ
z;#+27m?QqSt)|$sbuNDb3llcK3#DBP)&4~T`9k=TS0V@$drf5h^&ljQIWJONJtQOn
z5-GNO*D9tSQXk~~qnOE53^U#DfjDOtNinO^^PJ92;sVnU{}E<5
zO*W#tgp2S#Q~R+mh}{KZDRL;(p&ndHD=ch(Nf$l1C$ELVBnW!E^N*{v4{!AX>?D^F
z0*@rJ?L`kT@X%@c^)?W`g7RQBw~`Msb}Vwzdm5dDSJUmM>zH^AM=hq``=eQ@eh9??
zG2i(2`hQNKB}*qiyEjk-8!jiaHPG5BBVi&ruO%{j&NBn4PbiRIQ!>4^I{5;jIBlaa
zGIk4O6UIG>rmrR`bjhm#t2Sto1NhT%&9MOY%Vo|NpnwdKuDI1-!F@ElQbRBgukglI
zwgA0{{6}Ev7&)0M#579-25OG@?P$)BXK?@P!hT&aki4@4n|bZC3jsn9N7fylSH*Z&
z&Xv&=h+(xg5SzlM3xFf@={~&NTrP+!lJ<3b9ze3nTY#nSA9ir{J{T~QT`Hy2mz0Zs
z^5802sO1GR9PzUOOeJq292-bT_E-l1AN=sy6%(9*{3;Eexk^~DFMwyQ+jawb=N7bE
zCeDn?N}N>^s2_S_%>WSll5q!p#8bKn5W4`02awZ16EB3tva~Ax{c2^A`cTO#q2%!E
zw{*!2!+%GmuZ<976qNo0t{xs!QUo?KA|;;P?j*A7Ex|ra+5n+^6R=pzQ0uBLSgeCB
zWk8XC%A