natas_hw
/
jidt
mirror of https://github.com/jlizier/jidt
0
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
jidt/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiV...

911 lines
37 KiB
Java
Executable File
Raw Blame History

/*
* 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 <http://www.gnu.org/licenses/>.
*/
package infodynamics.measures.continuous.kraskov;
import infodynamics.utils.ArrayFileReader;
import infodynamics.utils.MathsUtils;
import infodynamics.utils.MatrixUtils;
import infodynamics.utils.RandomGenerator;
public class ConditionalMutualInfoMultiVariateTester
extends infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateAbstractTester {
protected String NUM_THREADS_TO_USE_DEFAULT = ConditionalMutualInfoCalculatorMultiVariateKraskov.USE_ALL_THREADS;
protected String NUM_THREADS_TO_USE = NUM_THREADS_TO_USE_DEFAULT;
/**
* Utility function to create a calculator for the given algorithm number
*
* @param algNumber
* @return
*/
public ConditionalMutualInfoCalculatorMultiVariateKraskov getNewCalc(int algNumber) {
ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = null;
if (algNumber == 1) {
condMiCalc = new ConditionalMutualInfoCalculatorMultiVariateKraskov1();
} else if (algNumber == 2) {
condMiCalc = new ConditionalMutualInfoCalculatorMultiVariateKraskov2();
}
return condMiCalc;
}
/**
* Confirm that the local values average correctly back to the average value
*
*/
public void checkLocalsAverageCorrectly(int algNumber) throws Exception {
ConditionalMutualInfoCalculatorMultiVariateKraskov miCalc = getNewCalc(algNumber);
String kraskov_K = "4";
miCalc.setProperty(
MutualInfoCalculatorMultiVariateKraskov.PROP_K,
kraskov_K);
super.testLocalsAverageCorrectly(miCalc, 2, 100);
}
public void testLocalsAverageCorrectly() throws Exception {
checkLocalsAverageCorrectly(1);
checkLocalsAverageCorrectly(2);
}
/**
* Confirm that significance testing doesn't alter the average that
* would be returned.
*
* @throws Exception
*/
public void checkComputeSignificanceDoesntAlterAverage(int algNumber) throws Exception {
ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = getNewCalc(algNumber);
String kraskov_K = "4";
condMiCalc.setProperty(
MutualInfoCalculatorMultiVariateKraskov.PROP_K,
kraskov_K);
super.testComputeSignificanceDoesntAlterAverage(condMiCalc, 2, 100);
}
public void testComputeSignificanceDoesntAlterAverage() throws Exception {
checkComputeSignificanceDoesntAlterAverage(1);
checkComputeSignificanceDoesntAlterAverage(2);
}
/**
* Utility function to run Kraskov conditional MI algorithm 1
* as transfer entropy for data with known results
* from TRENTOOL. (with default parameter settings k=1, l=1)
*
* @param var1 source multivariate data set
* @param var2 dest multivariate data set
* @param kNNs array of Kraskov k nearest neighbours parameter to check
* @param expectedResults array of expected results for each k
* @return errors of the computed values against expectedResults
*/
protected double[] checkTEForGivenData(double[][] var1, double[][] var2,
int[] kNNs, double[] expectedResults) throws Exception {
return checkTEForGivenData(var1, var2, 1, 1, kNNs, expectedResults);
}
/**
* Utility function to run Kraskov conditional MI algorithm 1
* as transfer entropy for data with known results
* from TRENTOOL.
*
* @param var1 source multivariate data set
* @param var2 dest multivariate data set
* @param historyK history length k of destination
* @param historyL history length l of source
* @param kNNs array of Kraskov k nearest neighbours parameter to check
* @param expectedResults array of expected results for each k
* @return errors of the computed values against expectedResults
*/
protected double[] checkTEForGivenData(double[][] var1, double[][] var2,
int historyK, int historyL, int[] kNNs, double[] expectedResults) throws Exception {
return checkTEForGivenData(var1, var2, historyK, historyL, kNNs, expectedResults,
0, "NONE", 0.000001);
}
/**
* Utility function to run Kraskov conditional MI algorithm 1
* as transfer entropy for data with known results
* from TRENTOOL.
*
* @param var1 source multivariate data set
* @param var2 dest multivariate data set
* @param historyK history length k of destination
* @param historyL history length l of source
* @param kNNs array of Kraskov k nearest neighbours parameter to check
* @param expectedResults array of expected results for each k
* @param noiseLevel noise to add to the data - set to 0 if we
* need to exactly reproduce calculations (most cases in unit tests for consistency)
* @param noiseSeed seed for the random noise generator (either "NONE" or Long string)
* @param tolerance tolerance to accept the calculation
* @return errors of the computed values against expectedResults
*/
protected double[] checkTEForGivenData(double[][] var1, double[][] var2,
int historyK, int historyL, int[] kNNs, double[] expectedResults,
double noiseLevel, String noiseSeed, double tolerance) throws Exception {
ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = getNewCalc(1);
double[] errors = new double[expectedResults.length];
// Which is the first time index for the dest next state?
// It depends on the values of k and l for embedding the past state
// of destination and source.
int firstDestTimeIndex = Math.max(historyK, historyL);
// Normalise the data ourselves rather than letting the calculator do it -
// this ensures the extra values in the time series (e.g. last value in source)
// are taken into account, in line with TRENTOOL
var1 = MatrixUtils.normaliseIntoNewArray(var1);
var2 = MatrixUtils.normaliseIntoNewArray(var2);
for (int kIndex = 0; kIndex < kNNs.length; kIndex++) {
int k = kNNs[kIndex];
condMiCalc.setProperty(
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_K,
Integer.toString(k));
// We already normalised above, and this will do a different
// normalisation without taking the extra values in to account if we did it
condMiCalc.setProperty(
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NORMALISE,
Boolean.toString(false));
// No longer need to set this property as it's set by default:
//condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE,
// EuclideanUtils.NORM_MAX_NORM_STRING);
condMiCalc.setProperty(
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
NUM_THREADS_TO_USE);
condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE,
Double.toString(noiseLevel));
condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NOISE_SEED, noiseSeed);
condMiCalc.initialise(var1[0].length * historyL,
var2[0].length, var2[0].length * historyK);
// Construct the joint vectors of the source states
double[][] sources = null;
if (historyL == 1) {
sources = MatrixUtils.selectRows(var1, firstDestTimeIndex - historyL,
var1.length - firstDestTimeIndex);
} else {
// Build the storage for the source states
int sourceVars = var1[0].length;
sources = new double[var1.length - firstDestTimeIndex][sourceVars * historyL];
for (int t = 0; t < historyL; t++) {
MatrixUtils.arrayCopy(
var1, firstDestTimeIndex - historyL + t, 0,
sources, 0, t*sourceVars,
var1.length - firstDestTimeIndex, sourceVars);
}
}
// Construct the joint vectors of the conditionals
double[][] conditionals = null;
if (historyK == 1) {
conditionals = MatrixUtils.selectRows(var2, firstDestTimeIndex - historyK,
var2.length - firstDestTimeIndex);
} else {
// Build the storage for the conditional observations
int destVars = var2[0].length;
conditionals = new double[var2.length - firstDestTimeIndex][destVars * historyK];
for (int t = 0; t < historyK; t++) {
MatrixUtils.arrayCopy(
var2, firstDestTimeIndex - historyK + t, 0,
conditionals, 0, t*destVars,
var2.length - firstDestTimeIndex, destVars);
}
}
// And set the observations using these
condMiCalc.setObservations(sources,
MatrixUtils.selectRows(var2, firstDestTimeIndex, var2.length - firstDestTimeIndex),
conditionals);
double condMi = condMiCalc.computeAverageLocalOfObservations();
//miCalc.setDebug(false);
System.out.printf("k=%d: Average MI %.8f (expected %.8f)\n",
k, condMi, expectedResults[kIndex]);
// 6 decimal places is Matlab accuracy
assertEquals(expectedResults[kIndex], condMi, tolerance);
errors[kIndex] = condMi - expectedResults[kIndex];
}
return errors;
}
/**
* Test the computed univariate TE as a conditional MI
* against that calculated by Wibral et al.'s TRENTOOL
* on the same data.
*
* To run TRENTOOL (http://www.trentool.de/) for this
* data, run its TEvalues.m matlab script on the multivariate source
* and dest data sets as:
* TEvalues(source, dest, 1, 1, 1, kraskovK, 0)
* with these values ensuring source-dest lag 1, history k=1,
* embedding lag 1, no dynamic correlation exclusion
*
* @throws Exception if file not found
*
*/
public void testUnivariateTEforCoupledVariablesFromFile() throws Exception {
// Test set 1:
ArrayFileReader afr = new ArrayFileReader("demos/data/2coupledRandomCols-1.txt");
double[][] data = afr.getDouble2DMatrix();
// Use various Kraskov k nearest neighbours parameter
int[] kNNs = {4};
// Expected values from TRENTOOL:
double[] expectedFromTRENTOOL = {0.3058006};
System.out.println("Kraskov Cond MI as TE comparison 1 - univariate coupled data 1");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
// And now in the reverse direction:
expectedFromTRENTOOL = new double[] {-0.0029744};
System.out.println(" reverse direction:");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {0}),
kNNs, expectedFromTRENTOOL);
}
/**
* Test the computed univariate TE as a conditional MI
* against that calculated by Wibral et al.'s TRENTOOL
* on the same data.
*
* To run TRENTOOL (http://www.trentool.de/) for this
* data, run its TEvalues.m matlab script on the multivariate source
* and dest data sets as:
* TEvalues(source, dest, 1, 1, 1, kraskovK, 0)
* with these values ensuring source-dest lag 1, history k=1,
* embedding lag 1, no dynamic correlation exclusion
*
* @throws Exception if file not found
*
*/
public void testUnivariateTEforCoupledLogisticMapFromFile() throws Exception {
// Test set 1:
ArrayFileReader afr = new ArrayFileReader("demos/data/coupledLogisticMapXY.txt");
double[][] data = afr.getDouble2DMatrix();
// Use various Kraskov k nearest neighbours parameter
int[] kNNs = {4};
// Expected values from TRENTOOL:
double[] expectedFromTRENTOOL = {0.508417};
System.out.println("Kraskov Cond MI as TE comparison 1 - univariate coupled logistic map data 1");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
// And now in the reverse direction:
expectedFromTRENTOOL = new double[] {0.016257};
System.out.println(" reverse direction:");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {0}),
kNNs, expectedFromTRENTOOL);
}
/**
* Test the computed univariate TE as a conditional MI
* against that calculated by Wibral et al.'s TRENTOOL
* on the same data.
*
* To run TRENTOOL (http://www.trentool.de/) for this
* data, run its TEvalues.m matlab script on the multivariate source
* and dest data sets as:
* TEvalues(source, dest, 1, 1, 1, kraskovK, 0)
* with these values ensuring source-dest lag 1, history k=1,
* embedding lag 1, no dynamic correlation exclusion
*
* @throws Exception if file not found
*
*/
public void testUnivariateTEforRandomDataFromFile() throws Exception {
// Test set 2:
ArrayFileReader afr = new ArrayFileReader("demos/data/4randomCols-1.txt");
double[][] data = afr.getDouble2DMatrix();
// Use various Kraskov k nearest neighbours parameter
int[] kNNs = {4};
// Expected values from TRENTOOL:
double[] expectedFromTRENTOOL = {-0.0096556};
System.out.println("Kraskov Cond MI as TE comparison 2 - univariate random data 1 (col 0->1)");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
// And now for other columns
expectedFromTRENTOOL = new double[] {0.0175389};
System.out.println(" (col 1->2):");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {2}),
kNNs, expectedFromTRENTOOL);
// And now for other columns
expectedFromTRENTOOL = new double[] {0.0026367};
System.out.println(" (col 1->0):");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {0}),
kNNs, expectedFromTRENTOOL);
// And now for other columns
expectedFromTRENTOOL = new double[] {-0.00012474};
System.out.println(" (col 0->2):");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {2}),
kNNs, expectedFromTRENTOOL);
// And now for other columns
expectedFromTRENTOOL = new double[] {-5.4437e-03};
System.out.println(" (col 2->0):");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {2}),
MatrixUtils.selectColumns(data, new int[] {0}),
kNNs, expectedFromTRENTOOL);
}
/**
* Test the computed multivariate conditional MI
* against that calculated by Wibral et al.'s TRENTOOL
* on the same data.
*
* It's multivariate because we use embedding dimension 2 on both source
* and destination.
*
* To run TRENTOOL (http://www.trentool.de/) for this
* data, run its TEvalues.m matlab script on the multivariate source
* and dest data sets as:
* TEvalues(source, dest, 2, 2, 1, kraskovK, 0)
* with these values ensuring source-dest lag 1, history k=2,
* history embedding dimension l=2 on source as well.
* embedding lag 1, no dynamic correlation exclusion
*
* @throws Exception if file not found
*
*/
public void testMultivariateCondMIforCoupledDataFromFile() throws Exception {
// Test set 3:
ArrayFileReader afr = new ArrayFileReader("demos/data/4ColsPairedOneStepNoisyDependence-1.txt");
double[][] data = afr.getDouble2DMatrix();
// Use various Kraskov k nearest neighbours parameter
int[] kNNs = {4};
// Expected values from TRENTOOL:
double[] expectedFromTRENTOOL = {0.1400645};
System.out.println("Kraskov Cond MI as TE - multivariate coupled data 1, k=2,l=2");
System.out.println(" (0->2)");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {2}),
2, 2,
kNNs, expectedFromTRENTOOL);
// And now for reverse direction:
expectedFromTRENTOOL = new double[] {-0.0181459};
System.out.println(" (2->0):");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {2}),
MatrixUtils.selectColumns(data, new int[] {0}),
2, 2,
kNNs, expectedFromTRENTOOL);
// And now for other columns:
expectedFromTRENTOOL = new double[] {0.1639186};
System.out.println(" (1->3):");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {3}),
2, 2,
kNNs, expectedFromTRENTOOL);
// And in reverse:
expectedFromTRENTOOL = new double[] {0.0036976};
System.out.println(" (3->1):");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {3}),
MatrixUtils.selectColumns(data, new int[] {1}),
2, 2,
kNNs, expectedFromTRENTOOL);
// -------------
// And finally, confirm that we get different results for k=1,l=1,
// which match TRENTOOL
expectedFromTRENTOOL = new double[] {0.0072169};
System.out.println(" (0->1) but with k=1,l=1:");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
1, 1,
kNNs, expectedFromTRENTOOL);
// And in reverse
expectedFromTRENTOOL = new double[] {0.0011738};
System.out.println(" (1->2) but with k=1,l=1:");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {2}),
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(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NORMALISE));
condMiCalc.initialise(var1[0].length, var2[0].length, condVar[0].length);
condMiCalc.setObservations(var1, var2, condVar);
@SuppressWarnings("unused")
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(ConditionalMutualInfoCalculatorMultiVariateKraskov.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);
@SuppressWarnings("unused")
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.
*
* Test the computed univariate TE as a conditional MI
* against that calculated by Wibral et al.'s TRENTOOL
* on the same data.
*
* To run TRENTOOL (http://www.trentool.de/) for this
* data, run its TEvalues.m matlab script on the multivariate source
* and dest data sets as:
* TEvalues(source, dest, 1, 1, 1, kraskovK, 0)
* with these values ensuring source-dest lag 1, history k=1,
* embedding lag 1, no dynamic correlation exclusion
*
* @throws Exception if file not found
*
*/
public void testUnivariateTEVariousNumberThreads() throws Exception {
ArrayFileReader afr = new ArrayFileReader("demos/data/4randomCols-1.txt");
double[][] data = afr.getDouble2DMatrix();
// Use various Kraskov k nearest neighbours parameter
int[] kNNs = {4};
// Expected values from TRENTOOL:
double[] expectedFromTRENTOOL = {-0.0096556};
System.out.println("Kraskov Cond MI as TE - multivariate coupled data 1, k=2,l=2 (0->2)");
System.out.println(" with various numbers of threads:");
System.out.println(" -- 1 thread:");
NUM_THREADS_TO_USE = "1";
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
System.out.println(" -- 2 threads:");
NUM_THREADS_TO_USE = "2";
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
System.out.println(" -- 3 threads:");
NUM_THREADS_TO_USE = "3";
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
System.out.println(" -- all threads:");
NUM_THREADS_TO_USE = ConditionalMutualInfoCalculatorMultiVariateKraskov.USE_ALL_THREADS;
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
// And finally test that multithreading is still ok if we have
// an imbalanced number of data between each thread.
// Expected value is only generated from our own code; we're not so
// interested in checking for this precise value, as we are in
// checking that the value is stable when we change the number of
// threads and have an uneven amount of data in each thread.
double[] expectedValue = new double[] {0.026517704};
NUM_THREADS_TO_USE = "2";
checkTEForGivenData(
MatrixUtils.selectRows(
MatrixUtils.selectColumns(data, new int[] {0}),
0, 501),
MatrixUtils.selectRows(
MatrixUtils.selectColumns(data, new int[] {1}),
0, 501),
kNNs, expectedValue);
NUM_THREADS_TO_USE = "3";
checkTEForGivenData(
MatrixUtils.selectRows(
MatrixUtils.selectColumns(data, new int[] {0}),
0, 501),
MatrixUtils.selectRows(
MatrixUtils.selectColumns(data, new int[] {1}),
0, 501),
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 <dataFile> 1 1 3000 <kNearestNeighbours> 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, (double[][]) 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);
}
}
/**
* Test that observationSetIndices and observationStartTimePoints are written properly
*
* @throws Exception
*/
public void testObservationSetIndices() throws Exception {
int dimensions = 1;
int timeSteps = 100;
ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = getNewCalc(1);
condMiCalc.initialise(dimensions, dimensions, dimensions);
// generate some random data
RandomGenerator rg = new RandomGenerator();
double[][] sourceData = rg.generateNormalData(timeSteps, dimensions,
0, 1);
double[][] destData = rg.generateNormalData(timeSteps, dimensions,
0, 1);
double[][] condData = rg.generateNormalData(timeSteps, dimensions,
0, 1);
// First check that for a simple single observation set everything works:
condMiCalc.setObservations(sourceData, destData, condData);
int[] observationSetIds = condMiCalc.getObservationSetIndices();
int[] timeSeriesIndices = condMiCalc.getObservationTimePoints();
assert(observationSetIds.length == timeSteps);
for (int t = 0; t < timeSteps; t++) {
assertEquals(observationSetIds[t], 0);
assertEquals(timeSeriesIndices[t], t);
}
// Now add the same one twice:
condMiCalc.initialise(dimensions, dimensions, dimensions);
condMiCalc.startAddObservations();
condMiCalc.addObservations(sourceData, destData, condData);
condMiCalc.addObservations(sourceData, destData, condData);
condMiCalc.finaliseAddObservations();
observationSetIds = condMiCalc.getObservationSetIndices();
timeSeriesIndices = condMiCalc.getObservationTimePoints();
assert(observationSetIds.length == 2*timeSteps);
for (int t = 0; t < timeSteps; t++) {
assertEquals(observationSetIds[t], 0);
assertEquals(timeSeriesIndices[t], t);
}
for (int t = 0; t < timeSteps; t++) {
assertEquals(observationSetIds[timeSteps + t], 1);
assertEquals(timeSeriesIndices[timeSteps + t], t);
}
// Now add NUM_SEGMENTS randomly chosen segments:
int NUM_SEGMENTS = 10;
int maxLength = 10;
int[] startPoints = rg.generateRandomInts(NUM_SEGMENTS, timeSteps - maxLength);
int[] lengthsMinus1 = rg.generateRandomInts(NUM_SEGMENTS, maxLength - 1); // ensures we don't add segments of length 0
condMiCalc.initialise(dimensions, dimensions, dimensions);
condMiCalc.startAddObservations();
for (int r = 0; r < NUM_SEGMENTS; r++) {
condMiCalc.addObservations(sourceData, destData, condData, startPoints[r], lengthsMinus1[r]+1);
}
condMiCalc.finaliseAddObservations();
observationSetIds = condMiCalc.getObservationSetIndices();
timeSeriesIndices = condMiCalc.getObservationTimePoints();
assert(observationSetIds.length == MatrixUtils.sum(lengthsMinus1) + NUM_SEGMENTS);
int t = 0;
for (int r = 0; r < NUM_SEGMENTS; r++) {
for (int i = 0; i < lengthsMinus1[r]+1; i++) {
assertEquals(observationSetIds[t], r);
assertEquals(timeSeriesIndices[t], startPoints[r] + i);
t++;
}
}
}
/**
* Extends testUnivariateTEforCoupledVariablesFromFile to test
* using seed for random number generator
*
* @throws Exception if file not found
*
*/
public void testWithSeed() throws Exception {
// Test set 1:
ArrayFileReader afr = new ArrayFileReader("demos/data/2coupledRandomCols-1.txt");
double[][] data = afr.getDouble2DMatrix();
// Use various Kraskov k nearest neighbours parameter
int[] kNNs = {4};
// Expected values from TRENTOOL:
double[] expectedFromTRENTOOL = {0.3058006};
System.out.println("Kraskov Cond MI as TE comparison 1 - univariate coupled data 1");
double[] noNoiseError = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
kNNs, expectedFromTRENTOOL);
double noNoiseResult = expectedFromTRENTOOL[0] + noNoiseError[0];
// And now in the reverse direction:
double[] expectedFromTRENTOOLRev = new double[] {-0.0029744};
System.out.println(" reverse direction:");
double[] noNoiseErrorRev = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {0}),
kNNs, expectedFromTRENTOOLRev);
double noNoiseResultRev = expectedFromTRENTOOLRev[0] + noNoiseErrorRev[0];
// Check that changing the random number generator still returns close to those with no noise
// results, but with larger tolerance
System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1");
double[] withSeed1Error = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
1, 1, kNNs, expectedFromTRENTOOL,
1e-8, "1", 0.01);
double[] withSeed1Results = new double[] {expectedFromTRENTOOL[0] + withSeed1Error[0]};
// Now check that the results are exact when we repeat with the same seed:
System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1 repeat want exact");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
MatrixUtils.selectColumns(data, new int[] {1}),
1, 1, kNNs, withSeed1Results,
1e-8, "1", 1e-10);
// And in reverse direction:
System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1 reverse");
double[] withSeed1ErrorRev = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {0}),
1, 1, kNNs, expectedFromTRENTOOLRev,
1e-8, "1", 0.01);
double[] withSeed1ResultsRev = new double[] {expectedFromTRENTOOLRev[0] + withSeed1ErrorRev[0]};
// Now check that the results are exact when we repeat with the same seed:
System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1 reverse repeat want exact");
checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
MatrixUtils.selectColumns(data, new int[] {0}),
1, 1, kNNs, withSeed1ResultsRev,
1e-8, "1", 1e-10);
}
}
Reference in New Issue View Git Blame Copy Permalink