Unknown property values are ignored.
@@ -260,6 +270,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements */ @Override public void setProperty(String propertyName, String propertyValue) { + + boolean propertySet = true; if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) { normalise = Boolean.parseBoolean(propertyValue); } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) { @@ -271,6 +283,20 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements addNoise = true; noiseLevel = Double.parseDouble(propertyValue); } + } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) { + if (propertyValue.equals(NOISE_NO_SEED_VALUE)) { + noiseSeedSet = false; + } else { + noiseSeedSet = true; + noiseSeed = Long.parseLong(propertyValue); + } + } else { + // No property was set here + propertySet = false; + } + if (debug && propertySet) { + System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName + + " to " + propertyValue); } } @@ -280,6 +306,12 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements return Boolean.toString(normalise); } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) { return Double.toString(noiseLevel); + } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) { + if (noiseSeedSet) { + return Long.toString(noiseSeed); + } else { + return NOISE_NO_SEED_VALUE; + } } else { // No property matches for this class return null; @@ -733,6 +765,9 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements // Add Gaussian noise of std dev noiseLevel to the data if required if (addNoise) { Random random = new Random(); + if (noiseSeedSet) { + random.setSeed(noiseSeed); + } for (int r = 0; r < var1Observations.length; r++) { for (int c = 0; c < dimensionsVar1; c++) { var1Observations[r][c] += diff --git a/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariate.java b/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariate.java index 8818840..49e5b43 100755 --- a/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariate.java +++ b/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariate.java @@ -91,6 +91,16 @@ public interface MutualInfoCalculatorMultiVariate * if the data is to be normalised, that will be done before adding this noise. */ public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD"; + /** + * Property name for the seed for the random number generator for noise to be + * added to the data (default is no seed) + */ + public static final String PROP_NOISE_SEED = "NOISE_SEED"; + /** + * Property value to indicate no seed for the random number generator for noise to be + * added to the data + */ + public static final String NOISE_NO_SEED_VALUE = "NONE"; /** *As per {@link #addObservations(double[][], double[][])}; diff --git a/java/source/infodynamics/measures/continuous/MutualInfoMultiVariateCommon.java b/java/source/infodynamics/measures/continuous/MutualInfoMultiVariateCommon.java index 3229b1d..1be3ea7 100755 --- a/java/source/infodynamics/measures/continuous/MutualInfoMultiVariateCommon.java +++ b/java/source/infodynamics/measures/continuous/MutualInfoMultiVariateCommon.java @@ -167,9 +167,19 @@ public abstract class MutualInfoMultiVariateCommon implements */ protected boolean addNoise = false; /** - * Amount of random Gaussian noise to add to the incoming data + * Amount of random Gaussian noise to add to the incoming data. + * 0 by default except for KSG estimators (where it is recommended + * and 1e-8 is used to match MILCA toolkit) */ protected double noiseLevel = (double) 0.0; + /** + * Has the user set a seed for the random noise + */ + protected boolean noiseSeedSet = false; + /** + * Seed that the user set for the random noise + */ + protected long noiseSeed = 0; /* (non-Javadoc) * @see infodynamics.measures.continuous.ChannelCalculatorCommon#initialise() @@ -222,6 +232,8 @@ public abstract class MutualInfoMultiVariateCommon implements * by Kraskov for the KSG method though, so for that estimator we * use 1e-8 to match the MILCA toolkit (though note it adds in * a random amount of noise in [0,noiseLevel) ). + *
Unknown property values are ignored.
@@ -250,6 +262,13 @@ public abstract class MutualInfoMultiVariateCommon implements addNoise = true; noiseLevel = Double.parseDouble(propertyValue); } + } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) { + if (propertyValue.equals(NOISE_NO_SEED_VALUE)) { + noiseSeedSet = false; + } else { + noiseSeedSet = true; + noiseSeed = Long.parseLong(propertyValue); + } } else { // No property was set here propertySet = false; @@ -270,6 +289,12 @@ public abstract class MutualInfoMultiVariateCommon implements return Boolean.toString(normalise); } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) { return Double.toString(noiseLevel); + } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) { + if (noiseSeedSet) { + return Long.toString(noiseSeed); + } else { + return NOISE_NO_SEED_VALUE; + } } else { // No property was recognised here return null; @@ -601,6 +626,9 @@ public abstract class MutualInfoMultiVariateCommon implements // Add Gaussian noise of std dev noiseLevel to the data if required if (addNoise) { Random random = new Random(); + if (noiseSeedSet) { + random.setSeed(noiseSeed); + } for (int r = 0; r < sourceObservations.length; r++) { for (int c = 0; c < dimensionsSource; c++) { sourceObservations[r][c] += diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java index cf0be00..b49055d 100755 --- a/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java +++ b/java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java @@ -98,10 +98,11 @@ public class ConditionalMutualInfoMultiVariateTester * @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 void checkTEForGivenData(double[][] var1, double[][] var2, + protected double[] checkTEForGivenData(double[][] var1, double[][] var2, int[] kNNs, double[] expectedResults) throws Exception { - checkTEForGivenData(var1, var2, 1, 1, kNNs, expectedResults); + return checkTEForGivenData(var1, var2, 1, 1, kNNs, expectedResults); } /** @@ -115,12 +116,38 @@ public class ConditionalMutualInfoMultiVariateTester * @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 void checkTEForGivenData(double[][] var1, double[][] var2, + 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. @@ -148,7 +175,9 @@ public class ConditionalMutualInfoMultiVariateTester condMiCalc.setProperty( ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS, NUM_THREADS_TO_USE); - condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency of results for unit test + 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 @@ -193,8 +222,10 @@ public class ConditionalMutualInfoMultiVariateTester 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, 0.000001); + assertEquals(expectedResults[kIndex], condMi, tolerance); + errors[kIndex] = condMi - expectedResults[kIndex]; } + return errors; } /** @@ -812,4 +843,68 @@ public class ConditionalMutualInfoMultiVariateTester } } } + + /** + * 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); +} } diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java index 1df1ed9..9f67033 100644 --- a/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java +++ b/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java @@ -96,6 +96,23 @@ public class MutualInfoMultiVariateTester checkComputeSignificanceDoesntAlterAverage(2); } + /** + * Utility function to run Kraskov MI for data with known results. + * Sets to use no noise in the calculation and a tolerance of 0.0000001 + * + * @param var1 + * @param var2 + * @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[] checkMIForGivenData(double[][] var1, double[][] var2, + int[] kNNs, double[] expectedResults) throws Exception { + // Dropping required accuracy by one order of magnitude, due + // to faster but slightly less accurate digamma estimator change + return checkMIForGivenData(var1, var2, kNNs, expectedResults, 0, "NONE", 0.0000001); + } + /** * Utility function to run Kraskov MI for data with known results * @@ -103,13 +120,20 @@ public class MutualInfoMultiVariateTester * @param var2 * @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 void checkMIForGivenData(double[][] var1, double[][] var2, - int[] kNNs, double[] expectedResults) throws Exception { - + protected double[] checkMIForGivenData(double[][] var1, double[][] var2, + int[] kNNs, double[] expectedResults, + double noiseLevel, String noiseSeed, double tolerance) throws Exception { + // The Kraskov MILCA toolkit MIhigherdim executable // uses algorithm 2 by default (this is what it means by rectangular): MutualInfoCalculatorMultiVariateKraskov miCalc = getNewCalc(2); + double[] errors = new double[expectedResults.length]; for (int kIndex = 0; kIndex < kNNs.length; kIndex++) { int k = kNNs[kIndex]; @@ -122,7 +146,9 @@ public class MutualInfoMultiVariateTester // No longer need to set this property as it's set by default: //miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE, // EuclideanUtils.NORM_MAX_NORM_STRING); - miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests + miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, + Double.toString(noiseLevel)); + miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_NOISE_SEED, noiseSeed); miCalc.initialise(var1[0].length, var2[0].length); miCalc.setObservations(var1, var2); miCalc.setDebug(true); @@ -131,10 +157,10 @@ public class MutualInfoMultiVariateTester System.out.printf("k=%d: Average MI %.8f (expected %.8f)\n", k, mi, expectedResults[kIndex]); - // Dropping required accuracy by one order of magnitude, due - // to faster but slightly less accurate digamma estimator change - assertEquals(expectedResults[kIndex], mi, 0.0000001); + assertEquals(expectedResults[kIndex], mi, tolerance); + errors[kIndex] = mi - expectedResults[kIndex]; } + return errors; } /** @@ -801,5 +827,47 @@ public class MutualInfoMultiVariateTester double miWithDynCorrExclAndSeperateSets = miCalc.computeAverageLocalOfObservations(); assertEquals(expectedMIWithoutDynCorrExcl, miWithDynCorrExclAndSeperateSets, 0.00001); } + + /** + * Extends the tests from testUnivariateMIforRandomVariablesFromFile to use + * a random seed for noise and check that results are repeatable. + * + * @throws Exception if file not found + * + */ + public void testUnivariateMIWithSeed() 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 - no seed"); + double[] noNoiseErrors = checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}), + MatrixUtils.selectColumns(data, new int[] {1}), + kNNs, expectedFromMILCA); + double[] noNoiseResults = MatrixUtils.add(expectedFromMILCA, noNoiseErrors); + + // 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 comparison 1 - univariate random data 1 - seed 1"); + double[] withSeed1Errors = checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}), + MatrixUtils.selectColumns(data, new int[] {1}), + kNNs, noNoiseResults, + 1e-8, "1", 0.01); + double[] withSeed1Results = MatrixUtils.add(noNoiseResults, withSeed1Errors); + // Now check that the results are exact when we repeat with the same seed: + System.out.println("\n Kraskov comparison 1 - univariate random data 1 - seed 1 repeat"); + checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}), + MatrixUtils.selectColumns(data, new int[] {1}), + kNNs, withSeed1Results, + 1e-8, "1", 1e-10); + } }