diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java index 758cd86..1b99c94 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java @@ -1,5 +1,6 @@ package infodynamics.measures.continuous.kraskov; +import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariate; import infodynamics.utils.EuclideanUtils; import infodynamics.utils.MatrixUtils; import infodynamics.utils.EmpiricalMeasurementDistribution; @@ -20,7 +21,8 @@ import infodynamics.utils.RandomGenerator; * * @author Joseph Lizier */ -public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov { +public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov + implements ConditionalMutualInfoCalculatorMultiVariate { /** * we compute distances to the kth neighbour in the joint space @@ -185,44 +187,57 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov { * The user should ensure that all values 0..N-1 are represented exactly once in the * array reordering and that no other values are included here. * + * @param variableToReorder 1 for variable 1, 2 for variable 2 * @param reordering * @return * @throws Exception */ - public abstract double computeAverageLocalOfObservations(int[] reordering) throws Exception; + public abstract double computeAverageLocalOfObservations(int variableToReorder, + int[] reordering) throws Exception; /** * Compute the significance of the mutual information of the previously supplied observations. - * We destroy the p(x,y,z) correlations, while retaining the p(x,z), p(y) marginals, to check how + * We destroy the p(x,y,z) correlations, by permuting the given variable, + * while retaining the joint distribution of the other variable + * and the conditional, and the marginal distribution of the + * permuted variable. This checks how * significant this conditional mutual information actually was. * * This is in the spirit of Chavez et. al., "Statistical assessment of nonlinear causality: * application to epileptic EEG signals", Journal of Neuroscience Methods 124 (2003) 113-128 * which was performed for Transfer entropy. * + * @param variableToReorder 1 for variable 1, 2 for variable 2 * @param numPermutationsToCheck * @return the proportion of MI scores from the distribution which have higher or equal MIs to ours. */ - public synchronized EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck) throws Exception { + public synchronized EmpiricalMeasurementDistribution computeSignificance( + int variableToReorder, int numPermutationsToCheck) throws Exception { // Generate the re-ordered indices: RandomGenerator rg = new RandomGenerator(); int[][] newOrderings = rg.generateDistinctRandomPerturbations(data1.length, numPermutationsToCheck); - return computeSignificance(newOrderings); + return computeSignificance(variableToReorder, newOrderings); } /** * Compute the significance of the mutual information of the previously supplied observations. - * We destroy the p(x,y,z) correlations, while retaining the p(x,z), p(y) marginals, to check how - * significant this mutual information actually was. + * We destroy the p(x,y,z) correlations, by permuting the given variable, + * while retaining the joint distribution of the other variable + * and the conditional, and the marginal distribution of the + * permuted variable. This checks how + * significant this conditional mutual information actually was. * * This is in the spirit of Chavez et. al., "Statistical assessment of nonlinear causality: * application to epileptic EEG signals", Journal of Neuroscience Methods 124 (2003) 113-128 * which was performed for Transfer entropy. * + * @param variableToReorder 1 for variable 1, 2 for variable 2 * @param newOrderings the specific new orderings to use * @return the proportion of conditional MI scores from the distribution which have higher or equal MIs to ours. */ - public EmpiricalMeasurementDistribution computeSignificance(int[][] newOrderings) throws Exception { + public EmpiricalMeasurementDistribution computeSignificance( + int variableToReorder, int[][] newOrderings) throws Exception { + int numPermutationsToCheck = newOrderings.length; if (!condMiComputed) { computeAverageLocalOfObservations(); @@ -235,7 +250,8 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov { int countWhereMiIsMoreSignificantThanOriginal = 0; for (int i = 0; i < numPermutationsToCheck; i++) { // Compute the MI under this reordering - double newMI = computeAverageLocalOfObservations(newOrderings[i]); + double newMI = computeAverageLocalOfObservations( + variableToReorder, newOrderings[i]); measDistribution.distribution[i] = newMI; if (debug){ System.out.println("New MI was " + newMI); @@ -256,7 +272,8 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov { public abstract double[] computeLocalOfPreviousObservations() throws Exception; - public double[] computeLocalUsingPreviousObservations(double[][] states1, double[][] states2) throws Exception { + public double[] computeLocalUsingPreviousObservations(double[][] states1, + double[][] states2, double[][] condStates) throws Exception { // If implemented, will need to incorporate any normalisation here // (normalising the incoming data the same way the previously // supplied observations were normalised). diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java index 4569b2d..d308970 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java @@ -28,32 +28,49 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1 * Compute the average conditional MI from the previously set observations */ public double computeAverageLocalOfObservations() throws Exception { - return computeAverageLocalOfObservations(null); + return computeAverageLocalOfObservations(1, null); } /** - * Compute what the average conditional MI would look like were the second time series reordered + * Compute what the average conditional MI would look like were the given + * time series reordered * as per the array of time indices in reordering. * The user should ensure that all values 0..N-1 are represented exactly once in the - * array reordering and that no other values are included here. + * array reordering and that no other values are included here. + * * If reordering is null, it is assumed there is no reordering of - * the y variable. + * the given variable. * - * @param reordering the reordered time steps of the y variable + * @param variableToReorder 1 for variable 1, 2 for variable 2 + * @param reordering the reordered time steps of the given variable * @return * @throws Exception */ - public double computeAverageLocalOfObservations(int[] reordering) throws Exception { + public double computeAverageLocalOfObservations(int variableToReorder, + int[] reordering) throws Exception { if (!tryKeepAllPairsNorms || (data1.length > MAX_DATA_SIZE_FOR_KEEP_ALL_PAIRS_NORM)) { - double[][] originalData2 = data2; + double[][] originalData; + if (variableToReorder == 1) { + originalData = data1; + } else { + originalData = data2; + } if (reordering != null) { - // Generate a new re-ordered data2 - data2 = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData2, reordering); + // Generate a new re-ordered data array + if (variableToReorder == 1) { + data1 = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData, reordering); + } else { + data2 = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData, reordering); + } } // Compute the MI double newMI = computeAverageLocalOfObservationsWhileComputingDistances(); - // restore data2 - data2 = originalData2; + // restore original data + if (variableToReorder == 1) { + data1 = originalData; + } else { + data2 = originalData; + } return newMI; } @@ -76,13 +93,21 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1 // using x, y and z norms to all neighbours // (note that norm of point t to itself will be set to infinity). - int tForY = (reordering == null) ? t : reordering[t]; + int tForReorderedVar = (reordering == null) ? t : reordering[t]; double[] jointNorm = new double[N]; for (int t2 = 0; t2 < N; t2++) { - int t2ForY = (reordering == null) ? t2 : reordering[t2]; + int t2ForReorderedVar = (reordering == null) ? t2 : reordering[t2]; // Joint norm is the max of all three marginals - jointNorm[t2] = Math.max(xNorms[t][t2], Math.max(yNorms[tForY][t2ForY], zNorms[t][t2])); + if (variableToReorder == 1) { + jointNorm[t2] = Math.max(xNorms[tForReorderedVar][t2ForReorderedVar], + Math.max(yNorms[t][t2], + zNorms[t][t2])); + } else { + jointNorm[t2] = Math.max(xNorms[t][t2], + Math.max(yNorms[tForReorderedVar][t2ForReorderedVar], + zNorms[t][t2])); + } } // Then find the kth closest neighbour, using a heuristic to // select whether to keep the k mins only or to do a sort. @@ -107,12 +132,21 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1 for (int t2 = 0; t2 < N; t2++) { if (zNorms[t][t2] < epsilon) { n_z++; - if (xNorms[t][t2] < epsilon) { - n_xz++; - } - int t2ForY = (reordering == null) ? t2 : reordering[t2]; - if (yNorms[tForY][t2ForY] < epsilon) { - n_yz++; + int t2ForReorderedVar = (reordering == null) ? t2 : reordering[t2]; + if (variableToReorder == 1) { + if (xNorms[tForReorderedVar][t2ForReorderedVar] < epsilon) { + n_xz++; + } + if (yNorms[t][t2] < epsilon) { + n_yz++; + } + } else { + if (xNorms[t][t2] < epsilon) { + n_xz++; + } + if (yNorms[tForReorderedVar][t2ForReorderedVar] < epsilon) { + n_yz++; + } } } } diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java index 14778be..94a3eab 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java @@ -38,24 +38,40 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 protected static final double CUTOFF_MULTIPLIER = 1.5; /** - * Compute what the average conditional MI would look like were the second time series reordered + * Compute what the average conditional MI would look like were the given + * time series reordered * as per the array of time indices in reordering. * The user should ensure that all values 0..N-1 are represented exactly once in the * array reordering and that no other values are included here. * - * @param reordering + * @param variableToReorder 1 for variable 1, 2 for variable 2 + * @param reordering the reordered time steps of the given variable * @return * @throws Exception */ - public double computeAverageLocalOfObservations(int[] reordering) throws Exception { + public double computeAverageLocalOfObservations(int variableToReorder, + int[] reordering) throws Exception { if (!tryKeepAllPairsNorms || (data1.length > MAX_DATA_SIZE_FOR_KEEP_ALL_PAIRS_NORM)) { - double[][] originalData2 = data2; - // Generate a new re-ordered data2 - data2 = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData2, reordering); + double[][] originalData; + if (variableToReorder == 1) { + originalData = data1; + } else { + originalData = data2; + } + // Generate a new re-ordered data array + if (variableToReorder == 1) { + data1 = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData, reordering); + } else { + data2 = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData, reordering); + } // Compute the MI double newMI = computeAverageLocalOfObservationsWhileComputingDistances(); - // restore data2 - data2 = originalData2; + // restore original data + if (variableToReorder == 1) { + data1 = originalData; + } else { + data2 = originalData; + } return newMI; } @@ -80,13 +96,19 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 // First get x and y and z norms to all neighbours // (note that norm of point t to itself will be set to infinity). - int tForY = reordering[t]; + int tForReorderedVar = reordering[t]; double[][] jointNorm = new double[N][2]; for (int t2 = 0; t2 < N; t2++) { - int t2ForY = reordering[t2]; - jointNorm[t2][JOINT_NORM_VAL_COLUMN] = Math.max(xNorms[t][t2], - Math.max(yNorms[tForY][t2ForY], zNorms[t][t2])); + int t2ForReorderedVar = reordering[t2]; + if (variableToReorder == 1) { + jointNorm[t2][JOINT_NORM_VAL_COLUMN] = Math.max( + xNorms[tForReorderedVar][t2ForReorderedVar], + Math.max(yNorms[t][t2], zNorms[t][t2])); + } else { + jointNorm[t2][JOINT_NORM_VAL_COLUMN] = Math.max(xNorms[t][t2], + Math.max(yNorms[tForReorderedVar][t2ForReorderedVar], zNorms[t][t2])); + } // And store the time step for back reference after the // array is sorted. jointNorm[t2][JOINT_NORM_TIMESTEP_COLUMN] = t2; @@ -112,11 +134,20 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 // Find eps_{x,y,z} as the maximum x and y and z norms amongst this set: for (int j = 0; j < k; j++) { int timeStepOfJthPoint = timeStepsOfKthMins[j]; - if (xNorms[t][timeStepOfJthPoint] > eps_x) { - eps_x = xNorms[t][timeStepOfJthPoint]; - } - if (yNorms[tForY][reordering[timeStepOfJthPoint]] > eps_y) { - eps_y = yNorms[tForY][reordering[timeStepOfJthPoint]]; + if (variableToReorder == 1) { + if (xNorms[tForReorderedVar][reordering[timeStepOfJthPoint]] > eps_x) { + eps_x = xNorms[tForReorderedVar][reordering[timeStepOfJthPoint]]; + } + if (yNorms[t][timeStepOfJthPoint] > eps_y) { + eps_y = yNorms[t][timeStepOfJthPoint]; + } + } else { + if (xNorms[t][timeStepOfJthPoint] > eps_x) { + eps_x = xNorms[t][timeStepOfJthPoint]; + } + if (yNorms[tForReorderedVar][reordering[timeStepOfJthPoint]] > eps_y) { + eps_y = yNorms[tForReorderedVar][reordering[timeStepOfJthPoint]]; + } } if (zNorms[t][timeStepOfJthPoint] > eps_z) { eps_z = zNorms[t][timeStepOfJthPoint]; @@ -133,11 +164,20 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 for (int t2 = 0; t2 < N; t2++) { if (zNorms[t][t2] <= eps_z) { n_z++; - if (xNorms[t][t2] <= eps_x) { - n_xz++; - } - if (yNorms[tForY][reordering[t2]] <= eps_y) { - n_yz++; + if (variableToReorder == 1) { + if (xNorms[tForReorderedVar][reordering[t2]] <= eps_x) { + n_xz++; + } + if (yNorms[t][t2] <= eps_y) { + n_yz++; + } + } else { + if (xNorms[t][t2] <= eps_x) { + n_xz++; + } + if (yNorms[tForReorderedVar][reordering[t2]] <= eps_y) { + n_yz++; + } } } }