From 3e017d1be4a3c7cceed7caa2d0e5be40809d8c6f Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Wed, 26 Mar 2014 04:07:13 +0000 Subject: [PATCH] Added unit tests for Rearchitecting Active Info Storage calculators to use a common parent class for data collection. Mainly includes testing kernel AIS calculator against the older direct implementation of this (which is now moved into the unit test area only) --- ...StorageCalculatorCorrelationIntegrals.java | 278 ++++++++++++++++++ ...tiveInfoStorageCalculatorKernelDirect.java | 247 ++++++++++++++++ .../kernel/ActiveInfoStorageTester.java | 63 ++++ .../kraskov/MutualInfoMultiVariateTester.java | 3 + .../infodynamics/utils/MatrixUtilsTest.java | 36 +++ 5 files changed, 627 insertions(+) create mode 100755 java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorCorrelationIntegrals.java create mode 100755 java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java create mode 100755 java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageTester.java diff --git a/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorCorrelationIntegrals.java b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorCorrelationIntegrals.java new file mode 100755 index 0000000..300ad8f --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorCorrelationIntegrals.java @@ -0,0 +1,278 @@ +package infodynamics.measures.continuous.kernel; + +import infodynamics.measures.continuous.ActiveInfoStorageCalculator; +import infodynamics.utils.MatrixUtils; + +import java.util.Vector; + + +public abstract class ActiveInfoStorageCalculatorCorrelationIntegrals + implements ActiveInfoStorageCalculator { + + /** + * Length of past history to consider + */ + protected int k = 1; + protected int totalObservations = 0; + protected boolean debug = false; + protected double lastAverage; + + /** + * Storage for source observations for addObservsations + */ + protected Vector vectorOfObservations; + + protected boolean addedMoreThanOneObservationSet; + + public ActiveInfoStorageCalculatorCorrelationIntegrals() { + } + + /** + * Initialise the calculator using the existing value of k + * + */ + public void initialise() throws Exception { + initialise(k); + } + + /** + * Initialise the calculator + * + * @param k Length of past history to consider + */ + public void initialise(int k) throws Exception { + this.k = k; + addedMoreThanOneObservationSet = false; + } + + /** + * Set the observations to compute the probabilities from + * + * @param observations + */ + public void setObservations(double[] observations) throws Exception { + startAddObservations(); + addObservations(observations); + finaliseAddObservations(); + } + + /** + * Elect to add in the observations from several disjoint time series. + * + */ + public void startAddObservations() { + vectorOfObservations = new Vector(); + } + + /** + * Add some more observations. + * Note that the arrays must not be over-written by the user + * until after finaliseAddObservations() has been called. + * + * @param observations + */ + public void addObservations(double[] observations) throws Exception { + if (vectorOfObservations == null) { + // startAddObservations was not called first + throw new RuntimeException("User did not call startAddObservations before addObservations"); + } + if (observations.length <= k) { + // we won't be taking any observations here + return; + } + vectorOfObservations.add(observations); + } + + /** + * Add some more observations. + * + * @param observations + * @param startTime first time index to take observations on + * @param numTimeSteps number of time steps to use + */ + public void addObservations(double[] observations, + int startTime, int numTimeSteps) throws Exception { + if (vectorOfObservations == null) { + // startAddObservations was not called first + throw new RuntimeException("User did not call startAddObservations before addObservations"); + } + if (numTimeSteps <= k) { + // We won't be taking any observations here + return; + } + double[] obsToAdd = new double[numTimeSteps]; + System.arraycopy(observations, startTime, obsToAdd, 0, numTimeSteps); + vectorOfObservations.add(obsToAdd); + } + + /** + * Sets the observations to compute the PDFs from. + * Cannot be called in conjunction with start/add/finaliseAddObservations. + * valid is a time series (with time indices the same as observations) + * indicating whether the observation at that point is valid. + * sourceValid is the same for the source + * + * @param observation observations for the source variable + * @param valid + */ + public void setObservations(double[] observations, + boolean[] valid) throws Exception { + + Vector startAndEndTimePairs = computeStartAndEndTimePairs(valid); + + // We've found the set of start and end times for this pair + startAddObservations(); + for (int[] timePair : startAndEndTimePairs) { + int startTime = timePair[0]; + int endTime = timePair[1]; + addObservations(observations, startTime, endTime - startTime + 1); + } + finaliseAddObservations(); + } + + /** + * Compute a vector of start and end pairs of time points, between which we have + * valid series of the observations. + * + * Made public so it can be used if one wants to compute the number of + * observations prior to setting the observations. + * + * @param valid + * @return + */ + public Vector computeStartAndEndTimePairs(boolean[] valid) { + // Scan along the data avoiding invalid values + int startTime = 0; + int endTime = 0; + boolean lookingForStart = true; + Vector startAndEndTimePairs = new Vector(); + for (int t = 0; t < valid.length; t++) { + if (lookingForStart) { + // Precondition: startTime holds a candidate start time + if (valid[t]) { + // This point is OK at the destination + if (t - startTime < k) { + // We're still checking the past history only, so + continue; + } else { + // We've got the full past history ok + // set a candidate endTime + endTime = t; + lookingForStart = false; + if (t == valid.length - 1) { + // we need to terminate now + int[] timePair = new int[2]; + timePair[0] = startTime; + timePair[1] = endTime; + startAndEndTimePairs.add(timePair); + // System.out.printf("t_s=%d, t_e=%d\n", startTime, endTime); + } + } + } else { + // We need to keep looking. + // Move the potential start time to the next point + startTime = t + 1; + } + } else { + // Precondition: startTime holds the start time for this set, + // endTime holds a candidate end time + // Check if we can include the current time step + boolean terminateSequence = false; + if (valid[t]) { + // We can extend + endTime = t; + } else { + terminateSequence = true; + } + if (t == valid.length - 1) { + // we need to terminate the sequence anyway + terminateSequence = true; + } + if (terminateSequence) { + // This section is done + int[] timePair = new int[2]; + timePair[0] = startTime; + timePair[1] = endTime; + startAndEndTimePairs.add(timePair); + // System.out.printf("t_s=%d, t_e=%d\n", startTime, endTime); + lookingForStart = true; + startTime = t + 1; + } + } + } + return startAndEndTimePairs; + } + + /** + * Generate a vector for each time step, containing the past k states of the destination. + * Does not include a vector for the first k time steps. + * + * @param destination + * @return array of vectors for each time step + */ + protected double[][] makeJointVectorForPast(double[] destination) { + try { + // We want one less delay vector here - we don't need the last k point, + // because there is no next state for these. + return MatrixUtils.makeDelayEmbeddingVector(destination, k, k-1, destination.length - k); + } catch (Exception e) { + // The parameters for the above call should be fine, so we don't expect to + // throw an Exception here - embed in a RuntimeException if it occurs + throw new RuntimeException(e); + } + } + + /** + * Compute the next values into a vector + * + * @param destination + * @return + */ + protected double[][] makeJointVectorForNext(double[] destination) { + double[][] destNextVectors = new double[destination.length - k][1]; + for (int t = k; t < destination.length; t++) { + destNextVectors[t - k][0] = destination[t]; + } + return destNextVectors; + } + + /** + * Generate a vector for each time step, containing the past k states of + * the destination, and the current state. + * Does not include a vector for the first k time steps. + * + * @param destination + * @return + */ + protected double[][] makeJointVectorForNextPast(double[] destination) { + // We want all delay vectors here + return MatrixUtils.makeDelayEmbeddingVector(destination, k+1); + } + + + public double getLastAverage() { + return lastAverage; + } + + public int getNumObservations() { + return totalObservations; + } + + public void setProperty(String propertyName, String propertyValue) throws Exception { + boolean propertySet = true; + if (propertyName.equalsIgnoreCase(K_PROP_NAME)) { + k = Integer.parseInt(propertyValue); + } else { + // No property was set + propertySet = false; + } + if (debug && propertySet) { + System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName + + " to " + propertyValue); + } + } + + public void setDebug(boolean debug) { + this.debug = debug; + } +} diff --git a/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java new file mode 100755 index 0000000..0f8c5db --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageCalculatorKernelDirect.java @@ -0,0 +1,247 @@ +package infodynamics.measures.continuous.kernel; + +import infodynamics.measures.continuous.ActiveInfoStorageCalculator; +import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.EmpiricalMeasurementDistribution; + +/** + * + *

+ * Implements an active information storage calculator using kernel estimation. + *

+ * + *

+ * Usage: + *

    + *
  1. Construct
    2. + *
  2. SetProperty() for each property
    3. + *
  3. intialise()
    4. + *
  4. setObservations(), or [startAddObservations(), addObservations()*, finaliseAddObservations()] + * Note: If not using setObservations(), the results from computeLocal or getSignificance + * are not likely to be particularly sensible.
    5. + *
  5. computeAverageLocalOfObservations() or ComputeLocalOfPreviousObservations()
    6. + *
+ *

+ * + *

+ * TODO Use only a single kernel estimator class for the joint space, and compute other + * probabilities from this. This will save much time. + *

+ * + *

+ * Matched against Oliver Obst's c++ implementation on 23/4/2010. + *

+ * + * @author Joseph Lizier + * @see ActiveInfoStorageCalculator + * @see ActiveInfoStorageCalculatorCorrelationIntegrals + * + */ +public class ActiveInfoStorageCalculatorKernelDirect + extends ActiveInfoStorageCalculatorCorrelationIntegrals { + + protected MutualInfoCalculatorMultiVariateKernel miKernel = null; + + // Keep joint vectors so we don't need to regenerate them + protected double[][] destNextVectors; + protected double[][] destPastVectors; + + /** + * Default value for epsilon + */ + public static final double DEFAULT_EPSILON = 0.25; + /** + * Kernel width + */ + private double epsilon = DEFAULT_EPSILON; + + /** + * Creates a new instance of the kernel-estimate style transfer entropy calculator + * + */ + public ActiveInfoStorageCalculatorKernelDirect() { + super(); + miKernel = new MutualInfoCalculatorMultiVariateKernel(); + } + + /** + * Initialises the calculator with the existing values for k and epsilon + * + */ + public void initialise() throws Exception { + initialise(k, epsilon); + } + + /** + * Initialises the calculator with the existing value for epsilon + * + * @param k history length + */ + public void initialise(int k) throws Exception { + this.k = k; + initialise(k, epsilon); + } + + /** + * Initialises the calculator + * + * @param k history length + * @param epsilon kernel width + */ + public void initialise(int k, double epsilon) throws Exception { + super.initialise(k); + this.epsilon = epsilon; + miKernel.initialise(k, 1, epsilon); + destPastVectors = null; + destNextVectors = null; + } + + /** + * Set properties for the calculator. + * Can include any of the accepted values for + * {@link MutualInfoCalculatorMultiVariateKernel#setProperty(String, String)} + * except {@link infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate#PROP_TIME_DIFF} + * + * @param propertyName + * @param propertyValue + * @throws Exception + */ + public void setProperty(String propertyName, String propertyValue) throws Exception { + super.setProperty(propertyName, propertyValue); + if (propertyName.equalsIgnoreCase(MutualInfoCalculatorMultiVariateKernel.EPSILON_PROP_NAME)) { + // Grab epsilon here - we need it for the initialisation: + epsilon = Double.parseDouble(propertyValue); + // It will get passed onto miKernel in our initialisation routine + if (debug) { + System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName + + " to " + propertyValue); + } + } if (propertyName.equalsIgnoreCase(MutualInfoCalculatorMultiVariate.PROP_TIME_DIFF)) { + throw new Exception("Cannot set " + MutualInfoCalculatorMultiVariate.PROP_TIME_DIFF + + " property on the ActiveInfoStorageCalculator"); + } else { + // Pass any other property through to the miKernel object + miKernel.setProperty(propertyName, propertyValue); + } + } + + /** + * Flag that the observations are complete, probability distribution functions can now be built. + * + */ + public void finaliseAddObservations() { + // First work out the size to allocate the joint vectors, and do the allocation: + totalObservations = 0; + for (double[] currentObservation : vectorOfObservations) { + totalObservations += currentObservation.length - k; + } + destPastVectors = new double[totalObservations][k]; + destNextVectors = new double[totalObservations][1]; + + // Construct the joint vectors from the given observations + int startObservation = 0; + for (double[] currentObservation : vectorOfObservations) { + double[][] currentDestPastVectors = makeJointVectorForPast(currentObservation); + MatrixUtils.arrayCopy(currentDestPastVectors, 0, 0, + destPastVectors, startObservation, 0, currentDestPastVectors.length, k); + double[][] currentDestNextVectors = makeJointVectorForNext(currentObservation); + MatrixUtils.arrayCopy(currentDestNextVectors, 0, 0, + destNextVectors, startObservation, 0, currentDestNextVectors.length, 1); + startObservation += currentObservation.length - k; + } + + // Now set the joint vectors in the kernel estimator + try { + miKernel.setObservations(destPastVectors, destNextVectors); + } catch (Exception e) { + // Should only throw where our vector sizes don't match + throw new RuntimeException(e); + } + + // Store whether there was more than one observation set: + addedMoreThanOneObservationSet = vectorOfObservations.size() > 1; + + // And clear the vector of observations + vectorOfObservations = null; + } + + /** + *

Computes the average Active Info Storage for the previously supplied observations

+ * + */ + public double computeAverageLocalOfObservations() throws Exception { + lastAverage = miKernel.computeAverageLocalOfObservations(); + return lastAverage; + } + + /** + * Computes the local active information storage for the previous supplied observations. + * + * Where more than one time series has been added, the array + * contains the local values for each tuple in the order in + * which they were added. + * + * If there was only a single time series added, the array + * contains k zero values before the local values. + * (This means the length of the return array is the same + * as the length of the input time series). + * + */ + public double[] computeLocalOfPreviousObservations() throws Exception { + double[] local = miKernel.computeLocalOfPreviousObservations(); + lastAverage = miKernel.getLastAverage(); + if (!addedMoreThanOneObservationSet) { + double[] localsToReturn = new double[local.length + k]; + System.arraycopy(local, 0, localsToReturn, k, local.length); + return localsToReturn; + } else { + return local; + } + } + + /** + * Compute the significance of obtaining the given average TE from the given observations + * + * This is as per Chavez et. al., "Statistical assessment of nonlinear causality: + * application to epileptic EEG signals", Journal of Neuroscience Methods 124 (2003) 113-128. + * + * Basically, we shuffle the source observations against the destination tuples. + * This keeps the marginal PDFs the same (including the entropy rate of the destination) + * but destroys any correlation between the source and state change of the destination. + * + * @param numPermutationsToCheck number of new orderings of the source values to compare against + * @return + */ + public EmpiricalMeasurementDistribution computeSignificance( + int numPermutationsToCheck) throws Exception { + return miKernel.computeSignificance(numPermutationsToCheck); + } + + /** + * As per {@link computeSignificance(int) computeSignificance()} but supplies + * the re-orderings of the observations of the source variables. + * + * + * @param newOrderings first index is permutation number, i.e. newOrderings[i] + * is an array of 1 permutation of 0..n-1, where there were n observations. + * @return + * @throws Exception + */ + public EmpiricalMeasurementDistribution computeSignificance( + int[][] newOrderings) throws Exception { + + return miKernel.computeSignificance(newOrderings); + } + + public void setDebug(boolean debug) { + super.setDebug(debug); + miKernel.setDebug(debug); + } + + public double[] computeLocalUsingPreviousObservations( + double[] newObservations) throws Exception { + // TODO Auto-generated method stub + return null; + } +} diff --git a/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageTester.java b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageTester.java new file mode 100755 index 0000000..0d1181c --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/kernel/ActiveInfoStorageTester.java @@ -0,0 +1,63 @@ +package infodynamics.measures.continuous.kernel; + +import infodynamics.utils.ArrayFileReader; +import infodynamics.utils.MatrixUtils; +import junit.framework.TestCase; + +public class ActiveInfoStorageTester extends TestCase { + + public void testEmbedding() throws Exception { + ArrayFileReader afr = new ArrayFileReader("demos/data/2coupledRandomCols-1.txt"); + double[][] data = afr.getDouble2DMatrix(); + + int[] ks = {1,2,3,4}; + int[] taus = {1,2,3,4}; + + for (int kIndex = 0; kIndex < ks.length; kIndex++) { + int k = ks[kIndex]; + for (int tauIndex = 0; tauIndex < taus.length; tauIndex++) { + int tau = taus[tauIndex]; + + ActiveInfoStorageCalculatorKernel ais = new ActiveInfoStorageCalculatorKernel(); + ais.initialise(k, tau, 0.5); + ais.setObservations(MatrixUtils.selectColumn(data, 0)); + + assertEquals(data.length - (k-1)*tau - 1, + ais.getNumObservations()); + } + } + } + + + public void testAgainstDirectAISCalculator() throws Exception { + + int[] ks = {1,2,3,4}; + double[] kernelwidths = {0.2, 0.5, 1}; + + ArrayFileReader afr = new ArrayFileReader("demos/data/2coupledRandomCols-1.txt"); + double[][] data = afr.getDouble2DMatrix(); + + for (int kIndex = 0; kIndex < ks.length; kIndex++) { + int k = ks[kIndex]; + for (int kwIndex = 0; kwIndex < kernelwidths.length; kwIndex++) { + double kernelwidth = kernelwidths[kwIndex]; + + ActiveInfoStorageCalculatorKernelDirect aisOld = new ActiveInfoStorageCalculatorKernelDirect(); + aisOld.initialise(k, kernelwidth); + aisOld.setObservations(MatrixUtils.selectColumn(data, 0)); + double aisOldValue = aisOld.computeAverageLocalOfObservations(); + + ActiveInfoStorageCalculatorKernel aisNew = new ActiveInfoStorageCalculatorKernel(); + aisNew.initialise(k, kernelwidth); + aisNew.setObservations(MatrixUtils.selectColumn(data, 0)); + double aisNewValue = aisNew.computeAverageLocalOfObservations(); + + + System.out.printf("k=%d, kw=%.2f, Kernel calculator: %.5f bits; Direct kernel calculator: %.5f bits\n", + k, kernelwidth, aisNewValue, aisOldValue); + assertEquals(aisOldValue, aisNewValue, 0.00000001); + } + } + + } +} diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java index d4161a7..97d7049 100755 --- a/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java +++ b/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java @@ -25,6 +25,9 @@ public class MutualInfoMultiVariateTester /** * Confirm that the local values average correctly back to the average value * + * TODO Add a test with say 10000 time steps, after we introduce fast nearest + * neighbour searching. + * */ public void checkLocalsAverageCorrectly(int algNumber) throws Exception { diff --git a/java/unittests/infodynamics/utils/MatrixUtilsTest.java b/java/unittests/infodynamics/utils/MatrixUtilsTest.java index b4b9494..f6c6f53 100755 --- a/java/unittests/infodynamics/utils/MatrixUtilsTest.java +++ b/java/unittests/infodynamics/utils/MatrixUtilsTest.java @@ -252,6 +252,42 @@ public class MatrixUtilsTest extends TestCase { checkArray(new int[] {1, 4, 0, 3, 2}, MatrixUtils.sortIndices(array3)); } + public void testDelayEmbeddings() throws Exception { + double[] array1 = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; + + // Do a standard delay embedding with tau 1 + checkMatrix(new double[][] { {4, 3, 2, 1, 0}, {5, 4, 3, 2, 1}, + {6, 5, 4, 3, 2}, {7, 6, 5, 4, 3}, + {8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} }, + MatrixUtils.makeDelayEmbeddingVector(array1, 5, 4, 6), + 0.00001); + // Now specify tau explicitly + checkMatrix(new double[][] { {4, 3, 2, 1, 0}, {5, 4, 3, 2, 1}, + {6, 5, 4, 3, 2}, {7, 6, 5, 4, 3}, + {8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} }, + MatrixUtils.makeDelayEmbeddingVector(array1, 5, 1, 4, 6), + 0.00001); + + // Do same standard delay embedding but starting at an offset + checkMatrix(new double[][] { {8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} }, + MatrixUtils.makeDelayEmbeddingVector(array1, 5, 8, 2), + 0.00001); + // Now specify tau explicitly + checkMatrix(new double[][] { {8, 7, 6, 5, 4}, {9, 8, 7, 6, 5} }, + MatrixUtils.makeDelayEmbeddingVector(array1, 5, 1, 8, 2), + 0.00001); + + // Try with tau 2 + checkMatrix(new double[][] { {8, 6, 4, 2, 0}, {9, 7, 5, 3, 1} }, + MatrixUtils.makeDelayEmbeddingVector(array1, 5, 2, 8, 2), + 0.00001); + + // Try with tau 3 + checkMatrix(new double[][] { {6, 3, 0}, {7, 4, 1}, {8, 5, 2}, {9, 6, 3} }, + MatrixUtils.makeDelayEmbeddingVector(array1, 3, 3, 6, 4), + 0.00001); + } + /** * Check that all entries in the given matrix match those of the expected * matrix