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)

This commit is contained in:
joseph.lizier 2014-03-26 04:07:13 +00:00
parent 2571c3ae9e
commit 3e017d1be4
5 changed files with 627 additions and 0 deletions

View File

@ -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<double[]> 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<double[]>();
}
/**
* 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<int[]> 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<int[]> computeStartAndEndTimePairs(boolean[] valid) {
// Scan along the data avoiding invalid values
int startTime = 0;
int endTime = 0;
boolean lookingForStart = true;
Vector<int[]> startAndEndTimePairs = new Vector<int[]>();
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;
}
}

View File

@ -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;
/**
*
* <p>
* Implements an active information storage calculator using kernel estimation.
* </p>
*
* <p>
* Usage:
* <ol>
* <li>Construct</li>
* <li>SetProperty() for each property</li>
* <li>intialise()</li>
* <li>setObservations(), or [startAddObservations(), addObservations()*, finaliseAddObservations()]
* Note: If not using setObservations(), the results from computeLocal or getSignificance
* are not likely to be particularly sensible.</li>
* <li>computeAverageLocalOfObservations() or ComputeLocalOfPreviousObservations()</li>
* </ol>
* </p>
*
* <p>
* TODO Use only a single kernel estimator class for the joint space, and compute other
* probabilities from this. This will save much time.
* </p>
*
* <p>
* Matched against Oliver Obst's c++ implementation on 23/4/2010.
* </p>
*
* @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;
}
/**
* <p>Computes the average Active Info Storage for the previously supplied observations</p>
*
*/
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;
}
}

View File

@ -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);
}
}
}
}

View File

@ -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 {

View File

@ -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