Unknown property values are ignored.
@@ -208,7 +248,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements // startAddObservations was not called first throw new RuntimeException("User did not call startAddObservations before addObservations"); } - if ((var1.length != var2.length) || (var1.length != cond.length)) { + if ((var1.length != var2.length) || + ((dimensionsCond != 0) && (var1.length != cond.length))) { throw new Exception(String.format("Observation vector lengths (%d, %d and %d) must match!", var1.length, var2.length, cond.length)); } @@ -220,7 +261,7 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements throw new Exception("Number of joint variables in var2 data " + "does not match the initialised value"); } - if (cond[0].length != dimensionsCond) { + if ((dimensionsCond != 0) && (cond[0].length != dimensionsCond)) { throw new Exception("Number of joint variables in cond data " + "does not match the initialised value"); } @@ -244,8 +285,11 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements System.arraycopy(var1, startTime, var1ToAdd, 0, numTimeSteps); double[][] var2ToAdd = new double[numTimeSteps][]; System.arraycopy(var2, startTime, var2ToAdd, 0, numTimeSteps); - double[][] condToAdd = new double[numTimeSteps][]; - System.arraycopy(cond, startTime, condToAdd, 0, numTimeSteps); + double[][] condToAdd = null; + if (dimensionsCond != 0) { + condToAdd = new double[numTimeSteps][]; + System.arraycopy(cond, startTime, condToAdd, 0, numTimeSteps); + } addObservations(var1ToAdd, var2ToAdd, condToAdd); } @@ -317,17 +361,27 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements MatrixUtils.arrayCopy(var2, 0, 0, var2Observations, startObservation, 0, var2.length, dimensionsVar2); - MatrixUtils.arrayCopy(cond, 0, 0, + if (dimensionsCond != 0) { + MatrixUtils.arrayCopy(cond, 0, 0, condObservations, startObservation, 0, cond.length, dimensionsCond); + } // else we can do nothing there startObservation += var2.length; } // Normalise the data if required if (normalise) { - MatrixUtils.normalise(var1Observations); - MatrixUtils.normalise(var2Observations); - MatrixUtils.normalise(condObservations); + var1Means = MatrixUtils.means(var1Observations); + var1Stds = MatrixUtils.stdDevs(var1Observations, var1Means); + MatrixUtils.normalise(var1Observations, var1Means, var1Stds); + var2Means = MatrixUtils.means(var2Observations); + var2Stds = MatrixUtils.stdDevs(var2Observations, var2Means); + MatrixUtils.normalise(var2Observations, var2Means, var2Stds); + if (dimensionsCond != 0) { + condMeans = MatrixUtils.means(condObservations); + condStds = MatrixUtils.stdDevs(condObservations, condMeans); + MatrixUtils.normalise(condObservations, condMeans, condStds); + } } // We don't need to keep the vectors of observation sets anymore: diff --git a/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java index c108cb9..8a7aed7 100755 --- a/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java +++ b/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java @@ -158,6 +158,22 @@ public class ConditionalMutualInfoCalculatorMultiVariateGaussian var2IndicesInCovariance = null; } + @Override + public void setProperty(String propertyName, String propertyValue) { + if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) { + // We do not allow the user to alter this property away from false, + // otherwise we would need to alter the functionality of the + // local calculations, etc. + return; + } + super.setProperty(propertyName, propertyValue); + } + + @Override + protected Object clone() throws CloneNotSupportedException { + return super.clone(); + } + /** * @throws Exception if the observation variables are not linearly independent * (leading to a non-positive definite covariance matrix). diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java index 99f48e9..7a1644e 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java @@ -46,7 +46,8 @@ import infodynamics.utils.UnivariateNearestNeighbourSearcher; * actually implement the two KSG algorithms. * *Crucially, the calculation is performed by examining - * neighbours in the full joint space (as specified by Frenzel and Pompe) + * neighbours in the full joint space (as specified by Frenzel and Pompe, + * and later by Vejmelka and Paluš, and Vlachos and Kugiumtzis) * rather than two MI calculators.
* *Usage is as per the paradigm outlined for {@link ConditionalMutualInfoCalculatorMultiVariate}, @@ -69,6 +70,12 @@ import infodynamics.utils.UnivariateNearestNeighbourSearcher; *
The method returns:
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