@@ -293,24 +295,30 @@
-
+
-
+
-
+
+
+
+
+
+
+
diff --git a/cuda/README-GPU.pdf b/cuda/README-GPU.pdf
new file mode 100644
index 0000000..819ae62
Binary files /dev/null and b/cuda/README-GPU.pdf differ
diff --git a/demos/AutoAnalyser/README-AutoAnalyserDemo.pdf b/demos/AutoAnalyser/README-AutoAnalyserDemo.pdf
index bdde033..6992bf7 100644
Binary files a/demos/AutoAnalyser/README-AutoAnalyserDemo.pdf and b/demos/AutoAnalyser/README-AutoAnalyserDemo.pdf differ
diff --git a/demos/AutoAnalyser/runAISAutoAnalyser.bat b/demos/AutoAnalyser/runAISAutoAnalyser.bat
deleted file mode 100755
index 11bdefc..0000000
--- a/demos/AutoAnalyser/runAISAutoAnalyser.bat
+++ /dev/null
@@ -1,5 +0,0 @@
-@ECHO OFF
-
-REM Run the example:
-java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserAIS
-
diff --git a/demos/AutoAnalyser/runAISAutoAnalyser.sh b/demos/AutoAnalyser/runAISAutoAnalyser.sh
deleted file mode 100755
index d0898eb..0000000
--- a/demos/AutoAnalyser/runAISAutoAnalyser.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/bash
-
-# Run the example:
-java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserAIS
-
diff --git a/demos/AutoAnalyser/runCMIAutoAnalyser.bat b/demos/AutoAnalyser/runCMIAutoAnalyser.bat
deleted file mode 100755
index 0cf7be1..0000000
--- a/demos/AutoAnalyser/runCMIAutoAnalyser.bat
+++ /dev/null
@@ -1,5 +0,0 @@
-@ECHO OFF
-
-REM Run the example:
-java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCMI
-
diff --git a/demos/AutoAnalyser/runCMIAutoAnalyser.sh b/demos/AutoAnalyser/runCMIAutoAnalyser.sh
deleted file mode 100755
index aad7edb..0000000
--- a/demos/AutoAnalyser/runCMIAutoAnalyser.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/bash
-
-# Run the example:
-java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCMI
-
diff --git a/demos/AutoAnalyser/runCTEAutoAnalyser.bat b/demos/AutoAnalyser/runCTEAutoAnalyser.bat
deleted file mode 100755
index 0e4b397..0000000
--- a/demos/AutoAnalyser/runCTEAutoAnalyser.bat
+++ /dev/null
@@ -1,5 +0,0 @@
-@ECHO OFF
-
-REM Run the example:
-java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCTE
-
diff --git a/demos/AutoAnalyser/runCTEAutoAnalyser.sh b/demos/AutoAnalyser/runCTEAutoAnalyser.sh
deleted file mode 100755
index 5d3e56a..0000000
--- a/demos/AutoAnalyser/runCTEAutoAnalyser.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/bash
-
-# Run the example:
-java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCTE
-
diff --git a/demos/AutoAnalyser/runHAutoAnalyser.bat b/demos/AutoAnalyser/runHAutoAnalyser.bat
deleted file mode 100755
index 825ded1..0000000
--- a/demos/AutoAnalyser/runHAutoAnalyser.bat
+++ /dev/null
@@ -1,5 +0,0 @@
-@ECHO OFF
-
-REM Run the example:
-java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserEntropy
-
diff --git a/demos/AutoAnalyser/runHAutoAnalyser.sh b/demos/AutoAnalyser/runHAutoAnalyser.sh
deleted file mode 100755
index 2969f73..0000000
--- a/demos/AutoAnalyser/runHAutoAnalyser.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/bash
-
-# Run the example:
-java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserEntropy
-
diff --git a/demos/AutoAnalyser/runMIAutoAnalyser.bat b/demos/AutoAnalyser/runMIAutoAnalyser.bat
deleted file mode 100755
index 383aaee..0000000
--- a/demos/AutoAnalyser/runMIAutoAnalyser.bat
+++ /dev/null
@@ -1,5 +0,0 @@
-@ECHO OFF
-
-REM Run the example:
-java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserMI
-
diff --git a/demos/AutoAnalyser/runMIAutoAnalyser.sh b/demos/AutoAnalyser/runMIAutoAnalyser.sh
deleted file mode 100755
index cb695ea..0000000
--- a/demos/AutoAnalyser/runMIAutoAnalyser.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/bash
-
-# Run the example:
-java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserMI
-
diff --git a/demos/AutoAnalyser/runTEAutoAnalyser.bat b/demos/AutoAnalyser/runTEAutoAnalyser.bat
deleted file mode 100755
index 8c5f7aa..0000000
--- a/demos/AutoAnalyser/runTEAutoAnalyser.bat
+++ /dev/null
@@ -1,5 +0,0 @@
-@ECHO OFF
-
-REM Run the example:
-java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserTE
-
diff --git a/demos/AutoAnalyser/runTEAutoAnalyser.sh b/demos/AutoAnalyser/runTEAutoAnalyser.sh
deleted file mode 100755
index be3e3fa..0000000
--- a/demos/AutoAnalyser/runTEAutoAnalyser.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/bash
-
-# Run the example:
-java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserTE
-
diff --git a/demos/README-Demos.pdf b/demos/README-Demos.pdf
index 36949ef..130e43c 100755
Binary files a/demos/README-Demos.pdf and b/demos/README-Demos.pdf differ
diff --git a/demos/java/example10GPUBenchmark.sh b/demos/java/example10GPUBenchmark.sh
old mode 100644
new mode 100755
diff --git a/demos/java/infodynamics/demos/Example10GPUBenchmark.java b/demos/java/infodynamics/demos/Example10GPUBenchmark.java
index d3d5269..d3d0a44 100755
--- a/demos/java/infodynamics/demos/Example10GPUBenchmark.java
+++ b/demos/java/infodynamics/demos/Example10GPUBenchmark.java
@@ -67,9 +67,10 @@ public class Example10GPUBenchmark {
miCalc.setProperty("k", "4");
miCalc.initialise(src[0].length, tgt[0].length);
miCalc.setObservations(src, tgt);
+ int nb_surrogates = 150;
double[] timeAndValue = new double[2];
long startTime = System.nanoTime();
- timeAndValue[1] = miCalc.computeAverageLocalOfObservations();
+ timeAndValue[1] = miCalc.computeSignificance(nb_surrogates).actualValue;
timeAndValue[0] = (System.nanoTime() - startTime)/1000000.0;
return timeAndValue;
}
diff --git a/demos/python/README-PythonDemos.pdf b/demos/python/README-PythonDemos.pdf
index db26c03..049fe85 100755
Binary files a/demos/python/README-PythonDemos.pdf and b/demos/python/README-PythonDemos.pdf differ
diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java
index e09335e..1482838 100755
--- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java
+++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java
@@ -717,9 +717,10 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
System.load(gpuLibraryPath);
}
} catch (Throwable e) {
- String errmsg = "GPU library not found. To compile GPU code set the enablegpu flag to true in build.xml";
+ String errmsg = "GPU library not found. To compile GPU code set the enablegpu flag to true in build.xml, or run `ant gpu jar`.";
+ errmsg += "\nFor more information see the JIDT GPU wiki page: https://github.com/jlizier/jidt/wiki/GPU";
if (gpuLibraryPath.length() > 0) {
- errmsg += "\nGPU library was not found in the path provided. Provide full path including library file name.";
+ errmsg += "\n\nGPU library was not found in the path provided. Provide full path including library file name.";
errmsg += "\nExample: /home/johndoe/myfolder/libKraskov.so";
}
throw new Exception(errmsg);
diff --git a/java/source/infodynamics/measures/spiking/TransferEntropyCalculatorSpiking.java b/java/source/infodynamics/measures/spiking/TransferEntropyCalculatorSpiking.java
index 5aaaf9f..630be29 100755
--- a/java/source/infodynamics/measures/spiking/TransferEntropyCalculatorSpiking.java
+++ b/java/source/infodynamics/measures/spiking/TransferEntropyCalculatorSpiking.java
@@ -282,68 +282,11 @@ public interface TransferEntropyCalculatorSpiking {
* will vary depending on the underlying implementation
*/
public SpikingLocalInformationValues computeLocalOfPreviousObservations() throws Exception;
-
- /**
- * Generate a bootstrapped distribution of what the TE would look like,
- * under a null hypothesis that the source values of our
- * samples had no temporal relation to the destination value.
- *
- * See Section II.E "Statistical significance testing" of
- * the JIDT paper below for a description of how this is done for MI and TE in general.
- *
- *
- * Note that if several disjoint time-series have been added
- * as observations using {@link #addObservations(double[])} etc.,
- * then these separate "trials" will be mixed up in the generation
- * of surrogates here.
- *
- * This method (in contrast to {@link #computeSignificance(int[][])})
- * creates random shufflings of the source embedding vectors for the surrogate
- * calculations.
- *
- * @param numPermutationsToCheck number of surrogate samples to bootstrap
- * to generate the distribution.
- * @return the distribution of TE scores under this null hypothesis.
- * @see "J.T. Lizier, 'JIDT: An information-theoretic
- * toolkit for studying the dynamics of complex systems', 2014."
- */
- public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck) throws Exception;
-
- /**
- * Generate a bootstrapped distribution of what the TE would look like,
- * under a null hypothesis that the source values of our
- * samples had no relation to the destination value.
- *
- * See Section II.E "Statistical significance testing" of
- * the JIDT paper below for a description of how this is done for MI and TE.
- *
- *
- * Note that if several disjoint time-series have been added
- * as observations using {@link #addObservations(double[])} etc.,
- * then these separate "trials" will be mixed up in the generation
- * of surrogates here.
- *
- * This method (in contrast to {@link #computeSignificance(int)})
- * allows the user to specify how to construct the surrogates,
- * such that repeatable results may be obtained.
- *
- * @param newOrderings a specification of how to shuffle the source embedding vectors
- * between all of our samples to create the surrogates to generate the distribution with. The first
- * index is the permutation number (i.e. newOrderings.length is the number
- * of surrogate samples we use to bootstrap to generate the distribution here.)
- * Each array newOrderings[i] should be an array of length N (where
- * would be the value returned by {@link #getNumObservations()}),
- * containing a permutation of the values in 0..(N-1).
- * TODO Need to think this through a little more before implementing.
- * @return the distribution of channel measure scores under this null hypothesis.
- * @see "J.T. Lizier, 'JIDT: An information-theoretic
- * toolkit for studying the dynamics of complex systems', 2014."
- * @throws Exception where the length of each permutation in newOrderings
- * is not equal to the number N samples that were previously supplied.
- */
- public EmpiricalMeasurementDistribution computeSignificance(
- int[][] newOrderings) throws Exception;
-
+
+ public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck, double estimatedValue) throws Exception;
+
+ public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck,
+ double estimatedValue, long randomSeed) throws Exception;
/**
* Set or clear debug mode for extra debug printing to stdout
*
diff --git a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java
index 690c397..89094d8 100644
--- a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java
+++ b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java
@@ -5,6 +5,15 @@ import java.util.Iterator;
import java.util.PriorityQueue;
import java.util.Vector;
+import java.util.Collections;
+import java.util.Iterator;
+import java.util.PriorityQueue;
+import java.util.Random;
+import java.util.Vector;
+import java.util.ArrayList;
+import java.lang.Math;
+
+//import infodynamics.measures.continuous.kraskov.EuclideanUtils;
import infodynamics.measures.spiking.TransferEntropyCalculatorSpiking;
import infodynamics.utils.EmpiricalMeasurementDistribution;
import infodynamics.utils.KdTree;
@@ -14,1217 +23,904 @@ import infodynamics.utils.NeighbourNodeData;
import infodynamics.utils.FirstIndexComparatorDouble;
import infodynamics.utils.UnivariateNearestNeighbourSearcher;
+import infodynamics.utils.EuclideanUtils;
+import infodynamics.utils.ParsedProperties;
+
/**
- * Computes the transfer entropy between a pair of spike trains,
- * using an integration-based measure in order to match the theoretical
- * form of TE between such spike trains.
+ * Computes the transfer entropy between a pair of spike trains, using an
+ * integration-based measure in order to match the theoretical form of TE
+ * between such spike trains.
*
- * Usage paradigm is as per the interface {@link TransferEntropyCalculatorSpiking}
+ *
+ * Usage paradigm is as per the interface
+ * {@link TransferEntropyCalculatorSpiking}
+ *
*
* @author Joseph Lizier (email,
- * www)
+ * www)
*/
-public class TransferEntropyCalculatorSpikingIntegration implements
- TransferEntropyCalculatorSpiking {
+
+/*
+ * TODO
+ * This implementation of the estimator does not implement dynamic exclusion windows. Such windows make sure
+ * that history embeddings that overlap are not considered in nearest-neighbour searches (as this breaks the
+ * independece assumption). Getting dynamic exclusion windows working will probably require modifications to the
+ * KdTree class.
+ */
+public class TransferEntropyCalculatorSpikingIntegration implements TransferEntropyCalculatorSpiking {
/**
- * Number of past destination spikes to consider (akin to embedding length)
+ * The past destination interspike intervals to consider
+ * (and associated property name and convenience length variable)
+ * The code assumes that the first interval is numbered 1, the next is numbered 2, etc.
+ * It also assumes that the intervals are sorted. The setter method performs sorting to ensure this.
*/
- protected int k = 1;
+ public static final String DEST_PAST_INTERVALS_PROP_NAME = "DEST_PAST_INTERVALS";
+ protected int[] destPastIntervals = new int[] {};
+ protected int numDestPastIntervals = 0;
/**
- * Number of past source spikes to consider (akin to embedding length)
+ * As above but for the sources
*/
- protected int l = 1;
+ public static final String SOURCE_PAST_INTERVALS_PROP_NAME = "SOURCE_PAST_INTERVALS";
+ protected int[] sourcePastIntervals = new int[] {};
+ protected int numSourcePastIntervals = 0;
+
+ /**
+ * As above but for the conditioning processes. There is no property name for this variable,
+ * due to there not currently being a method for converting strings to 2d arrays. Instead,
+ * a separate setter method is implemented.
+ */
+ protected Vector vectorOfCondPastIntervals = new Vector();
+ protected int numCondPastIntervals = 0;
/**
* Number of nearest neighbours to search for in the full joint space
*/
protected int Knns = 4;
-
+
+
/**
- * Storage for source observations supplied via {@link #addObservations(double[], double[])} etc.
+ * Storage for source observations supplied via
+ * {@link #addObservations(double[], double[])} etc.
+
*/
protected Vector vectorOfSourceSpikeTimes = null;
/**
- * Storage for destination observations supplied via {@link #addObservations(double[], double[])} etc.
+ * Storage for destination observations supplied via
+ * {@link #addObservations(double[], double[])} etc.
*/
protected Vector vectorOfDestinationSpikeTimes = null;
-
- // constants for indexing our data storage
- protected final static int NEXT_DEST = 0;
- protected final static int NEXT_SOURCE = 1;
- protected final static int NEXT_POSSIBILITIES = 2;
/**
- * Cache of the timing data for each new observed spiking event in both the source
- * and destination
+ * Storage for conditional observations supplied via
+ * {@link #addObservations(double[], double[])} etc.
*/
- Vector[] eventTimings = null;
- /**
- * Cache of the timing data for each new observed spiking event for the
- * destination only
- */
- Vector destPastAndNextTimings = null;
- /**
- * Cache of the type of event for each new observed spiking event in both the source
- * and destination (i.e. which spiked next)
- */
- Vector eventTypeLocator = null;
- /**
- * Cache for each new observed spiking event of which index it has in the vector
- * of spiking events of the same type
- */
- Vector eventIndexLocator = null;
- /**
- * Cache for each time-series of observed spiking events of how many
- * observations were in that set.
- */
- Vector numEventsPerObservationSet = null;
+ protected Vector vectorOfConditionalSpikeTimes = null;
- /**
- * KdTrees for searching the joint past spaces and time to next spike,
- * for each possibility of which spiked next
- */
- protected KdTree[] kdTreesJoint = null;
+ Vector conditioningEmbeddingsFromSpikes = null;
+ Vector jointEmbeddingsFromSpikes = null;
+ Vector conditioningEmbeddingsFromSamples = null;
+ Vector jointEmbeddingsFromSamples = null;
+ Vector processTimeLengths = null;
- /**
- * KdTrees for searching the joint past spaces,
- * for each possibility of which spiked next
- */
- protected KdTree[] kdTreesSourceDestHistories = null;
+ protected KdTree kdTreeJointAtSpikes = null;
+ protected KdTree kdTreeJointAtSamples = null;
+ protected KdTree kdTreeConditioningAtSpikes = null;
+ protected KdTree kdTreeConditioningAtSamples = null;
- /**
- * KdTrees for searching the past destination space and time to next spike
- */
- protected KdTree kdTreeDestNext = null;
-
- /**
- * KdTrees for searching the past destination space
- */
- protected KdTree kdTreeDestHistory = null;
-
- /**
- * NN searcher for the time to next spike space only, if required
- */
- protected UnivariateNearestNeighbourSearcher nnSearcherDestTimeToNextSpike = null;
-
- /**
- * Property name for the number of nearest neighbours to search
- */
public static final String KNNS_PROP_NAME = "Knns";
/**
- * Property name for adjusting the search radius for the next spike such that
- * it does not cover negative times (with respect to the previous spike, being either
- * source or destination spike)
+ * Property name for an amount of random Gaussian noise to be added to the data
+ * (default is 1e-8, matching the MILCA toolkit).
*/
- public static final String TRIM_TO_POS_PROP_NAME = "TRIM_RANGE_TO_POS_TIMES";
-
- protected boolean trimToPosNextSpikeTimes = false;
+ public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD";
+ /**
+ * Whether to add an amount of random noise to the incoming data
+ */
+ protected boolean addNoise = true;
+ /**
+ * Amount of random Gaussian noise to add to the incoming data
+ */
+ protected double noiseLevel = (double) 1e-8;
+
/**
* Stores whether we are in debug mode
*/
protected boolean debug = false;
+
+ /**
+ * Property name for the number of random sample points to use as a multiple
+ * of the number of target spikes.
+ */
+ public static final String PROP_SAMPLE_MULTIPLIER = "NUM_SAMPLES_MULTIPLIER";
+ protected double numSamplesMultiplier = 2.0;
+ /**
+ * Property name for the number of random sample points to use in the construction of the surrogates as a multiple
+ * of the number of target spikes.
+ */
+ public static final String PROP_SURROGATE_SAMPLE_MULTIPLIER = "SURROGATE_NUM_SAMPLES_MULTIPLIER";
+ protected double surrogateNumSamplesMultiplier = 2.0;
+
+ /**
+ * Property for the number of nearest neighbours to use in the construction of the surrogates
+ */
+ public static final String PROP_K_PERM = "K_PERM";
+ protected int kPerm = 10;
+
+
+ /**
+ * Property name for what type of norm to use between data points
+ * for each marginal variable -- Options are defined by
+ * {@link KdTree#setNormType(String)} and the
+ * default is {@link EuclideanUtils#NORM_EUCLIDEAN}.
+ */
+ public final static String PROP_NORM_TYPE = "NORM_TYPE";
+ protected int normType = EuclideanUtils.NORM_EUCLIDEAN;
public TransferEntropyCalculatorSpikingIntegration() {
super();
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(int)
+ /*
+ * (non-Javadoc)
+ *
+ * @see
+ * infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(
+ * int)
*/
@Override
public void initialise() throws Exception {
- initialise(k,l);
+ initialise(0, 0);
}
-
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(int)
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see
+ * infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(
+ * int)
*/
@Override
public void initialise(int k) throws Exception {
- initialise(k,this.l);
+ initialise(0, 0);
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(int, int)
+ /*
+ * (non-Javadoc)
+ *
+ * @see
+ * infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(
+ * int, int)
*/
@Override
public void initialise(int k, int l) throws Exception {
- if ((k < 1) || (l < 1)) {
- throw new Exception("Zero history length not supported");
- }
- this.k = k;
- this.l = l;
vectorOfSourceSpikeTimes = null;
vectorOfDestinationSpikeTimes = null;
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#setProperty(java.lang.String, java.lang.String)
+ /*
+ * (non-Javadoc)
+ *
+ * @see
+ * infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#setProperty(
+ * java.lang.String, java.lang.String)
*/
@Override
- public void setProperty(String propertyName, String propertyValue)
- throws Exception {
+ public void setProperty(String propertyName, String propertyValue) throws Exception {
boolean propertySet = true;
- if (propertyName.equalsIgnoreCase(K_PROP_NAME)) {
- k = Integer.parseInt(propertyValue);
- } else if (propertyName.equalsIgnoreCase(L_PROP_NAME)) {
- l = Integer.parseInt(propertyValue);
+ if (propertyName.equalsIgnoreCase(DEST_PAST_INTERVALS_PROP_NAME)) {
+ if (propertyValue.length() == 0) {
+ destPastIntervals = new int[] {};
+ } else {
+ int[] destPastIntervalsTemp = ParsedProperties.parseStringArrayOfInts(propertyValue);
+ for (int interval : destPastIntervalsTemp) {
+ if (interval < 1) {
+ throw new Exception ("Invalid interval number less than 1.");
+ }
+ }
+ destPastIntervals = destPastIntervalsTemp;
+ Arrays.sort(destPastIntervals);
+ }
+ } else if (propertyName.equalsIgnoreCase(SOURCE_PAST_INTERVALS_PROP_NAME)) {
+ if (propertyValue.length() == 0) {
+ sourcePastIntervals = new int[] {};
+ } else {
+ int[] sourcePastIntervalsTemp = ParsedProperties.parseStringArrayOfInts(propertyValue);
+ for (int interval : sourcePastIntervalsTemp) {
+ if (interval < 1) {
+ throw new Exception ("Invalid interval number less than 1.");
+ }
+ }
+ sourcePastIntervals = sourcePastIntervalsTemp;
+ Arrays.sort(sourcePastIntervals);
+ }
} else if (propertyName.equalsIgnoreCase(KNNS_PROP_NAME)) {
Knns = Integer.parseInt(propertyValue);
- } else if (propertyName.equalsIgnoreCase(TRIM_TO_POS_PROP_NAME)) {
- trimToPosNextSpikeTimes = Boolean.parseBoolean(propertyValue);
+ } else if (propertyName.equalsIgnoreCase(PROP_K_PERM)) {
+ kPerm = Integer.parseInt(propertyValue);
+ } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
+ if (propertyValue.equals("0") || propertyValue.equalsIgnoreCase("false")) {
+ addNoise = false;
+ noiseLevel = 0;
+ } else {
+ addNoise = true;
+ noiseLevel = Double.parseDouble(propertyValue);
+ }
+ } else if (propertyName.equalsIgnoreCase(PROP_SAMPLE_MULTIPLIER)) {
+ double tempNumSamplesMultiplier = Double.parseDouble(propertyValue);
+ if (tempNumSamplesMultiplier <= 0) {
+ throw new Exception ("Num samples multiplier must be greater than 0.");
+ } else {
+ numSamplesMultiplier = tempNumSamplesMultiplier;
+ }
+ } else if (propertyName.equalsIgnoreCase(PROP_NORM_TYPE)) {
+ normType = KdTree.validateNormType(propertyValue);
+ } else if (propertyName.equalsIgnoreCase(PROP_SURROGATE_SAMPLE_MULTIPLIER)) {
+ double tempSurrogateNumSamplesMultiplier = Double.parseDouble(propertyValue);
+ if (tempSurrogateNumSamplesMultiplier <= 0) {
+ throw new Exception ("Surrogate Num samples multiplier must be greater than 0.");
+ } else {
+ surrogateNumSamplesMultiplier = tempSurrogateNumSamplesMultiplier;
+ }
} else {
// No property was set on this class
propertySet = false;
}
if (debug && propertySet) {
- System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
- " to " + propertyValue);
+ System.out.println(
+ this.getClass().getSimpleName() + ": Set property " + propertyName + " to " + propertyValue);
}
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getProperty(java.lang.String)
+ /*
+ * (non-Javadoc)
+ *
+ * @see
+ * infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getProperty(
+ * java.lang.String)
*/
@Override
public String getProperty(String propertyName) throws Exception {
- if (propertyName.equalsIgnoreCase(K_PROP_NAME)) {
- return Integer.toString(k);
- } else if (propertyName.equalsIgnoreCase(L_PROP_NAME)) {
- return Integer.toString(l);
- } else if (propertyName.equalsIgnoreCase(KNNS_PROP_NAME)) {
+ if (propertyName.equalsIgnoreCase(KNNS_PROP_NAME)) {
return Integer.toString(Knns);
- } else if (propertyName.equalsIgnoreCase(TRIM_TO_POS_PROP_NAME)) {
- return Boolean.toString(trimToPosNextSpikeTimes);
+ } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
+ return Double.toString(noiseLevel);
+ } else if (propertyName.equalsIgnoreCase(PROP_SAMPLE_MULTIPLIER)) {
+ return Double.toString(numSamplesMultiplier);
} else {
// No property matches for this class
return null;
}
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#setObservations(double[], double[])
+ public void appendConditionalIntervals(int[] intervals) throws Exception{
+ for (int interval : intervals) {
+ if (interval < 1) {
+ throw new Exception ("Invalid interval number less than 1.");
+ }
+ }
+ Arrays.sort(intervals);
+ vectorOfCondPastIntervals.add(intervals);
+ }
+
+ public void clearConditionalIntervals() {
+ vectorOfCondPastIntervals = new Vector();
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#
+ * setObservations(double[], double[])
*/
@Override
- public void setObservations(double[] source, double[] destination)
- throws Exception {
+ public void setObservations(double[] source, double[] destination) throws Exception {
startAddObservations();
addObservations(source, destination);
finaliseAddObservations();
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#startAddObservations()
+ public void setObservations(double[] source, double[] destination, double[][] conditionals) throws Exception {
+ startAddObservations();
+ addObservations(source, destination, conditionals);
+ finaliseAddObservations();
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#
+ * startAddObservations()
*/
@Override
public void startAddObservations() {
vectorOfSourceSpikeTimes = new Vector();
vectorOfDestinationSpikeTimes = new Vector();
+ vectorOfConditionalSpikeTimes = new Vector();
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#addObservations(double[], double[])
+ /*
+ * (non-Javadoc)
+ *
+ * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#
+ * addObservations(double[], double[])
*/
@Override
- public void addObservations(double[] source, double[] destination)
- throws Exception {
+ public void addObservations(double[] source, double[] destination) throws Exception {
// Store these observations in our vector for now
vectorOfSourceSpikeTimes.add(source);
vectorOfDestinationSpikeTimes.add(destination);
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#finaliseAddObservations()
+ public void addObservations(double[] source, double[] destination, double[][] conditionals) throws Exception {
+ // Store these observations in our vector for now
+ vectorOfSourceSpikeTimes.add(source);
+ vectorOfDestinationSpikeTimes.add(destination);
+ vectorOfConditionalSpikeTimes.add(conditionals);
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#
+ * finaliseAddObservations()
*/
@Override
public void finaliseAddObservations() throws Exception {
- // TODO Auto embed if required
- // preFinaliseAddObservations();
- // Run through each spiking time series set and pull out the observation
- // tuples we'll store.
- // Initialise our data stores:
- eventTimings = new Vector[NEXT_POSSIBILITIES];
- for (int next = 0; next < NEXT_POSSIBILITIES; next++) {
- eventTimings[next] = new Vector();
+ // Set these conveniance variables as they are used quite a bit later on
+ numDestPastIntervals = destPastIntervals.length;
+ numSourcePastIntervals = sourcePastIntervals.length;
+ numCondPastIntervals = 0;
+ for (int[] intervals : vectorOfCondPastIntervals) {
+ numCondPastIntervals += intervals.length;
}
- destPastAndNextTimings = new Vector();
- eventTypeLocator = new Vector();
- eventIndexLocator = new Vector();
- numEventsPerObservationSet = new Vector();
-
+
+ conditioningEmbeddingsFromSpikes = new Vector();
+ jointEmbeddingsFromSpikes = new Vector();
+ conditioningEmbeddingsFromSamples = new Vector();
+ jointEmbeddingsFromSamples = new Vector();
+ processTimeLengths = new Vector();
+
// Send all of the observations through:
Iterator sourceIterator = vectorOfSourceSpikeTimes.iterator();
+ Iterator conditionalIterator = null;
+ if (vectorOfConditionalSpikeTimes.size() > 0) {
+ conditionalIterator = vectorOfConditionalSpikeTimes.iterator();
+ }
int timeSeriesIndex = 0;
for (double[] destSpikeTimes : vectorOfDestinationSpikeTimes) {
double[] sourceSpikeTimes = sourceIterator.next();
- timeSeriesIndex++;
-
- processEventsFromSpikingTimeSeries(sourceSpikeTimes, destSpikeTimes,
- timeSeriesIndex, eventTimings, destPastAndNextTimings,
- eventTypeLocator, eventIndexLocator, numEventsPerObservationSet);
- }
-
- // Now we have collected all the events.
- // Load up the search structures:
- // 1. Full joint space:
- // 2. Histories of source and dest only:
- kdTreesJoint = new KdTree[NEXT_POSSIBILITIES];
- kdTreesSourceDestHistories = new KdTree[NEXT_POSSIBILITIES];
- for (int next = 0; next < NEXT_POSSIBILITIES; next++) {
- // This line does not work:
- // double[][][] jointEventTimings = (double[][][]) eventTimings[prev][next].toArray();
- // So we'll do it manually:
- double[][] sourcePastTimings = new double[eventTimings[next].size()][];
- double[][] destPastTimings = new double[eventTimings[next].size()][];
- double[][] nextTimings = new double[eventTimings[next].size()][];
- int i = 0;
- for (double[][] timing : eventTimings[next]) {
- sourcePastTimings[i] = timing[0];
- destPastTimings[i] = timing[1];
- nextTimings[i] = timing[2];
- i++;
+ double[][] conditionalSpikeTimes = null;
+ if (vectorOfConditionalSpikeTimes.size() > 0) {
+ conditionalSpikeTimes = conditionalIterator.next();
+ } else {
+ conditionalSpikeTimes = new double[][] {};
}
- // TODO Should we normalise before we supply to the KdTree?
- // Think about this later. I'm not convinced it's the best
- // approach in this particular case.
- kdTreesJoint[next] = new KdTree(
- new int[] {l, k - 1, 1},
- new double[][][] {sourcePastTimings, destPastTimings, nextTimings});
- kdTreesSourceDestHistories[next] = new KdTree(
- new int[] {l, k - 1},
- new double[][][] {sourcePastTimings, destPastTimings});
+ processEventsFromSpikingTimeSeries(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes, conditioningEmbeddingsFromSpikes,
+ jointEmbeddingsFromSpikes, conditioningEmbeddingsFromSamples, jointEmbeddingsFromSamples,
+ numSamplesMultiplier);
}
- // 3. For the dest past and time to next spike
- // 4. For the dest past only
- double[][] destPastOnlyTimings = new double[destPastAndNextTimings.size()][];
- double[][] nextTimingsForDestPastOnly = new double[destPastAndNextTimings.size()][];
- int i = 0;
- for (double[][] timing : destPastAndNextTimings) {
- destPastOnlyTimings[i] = timing[0];
- nextTimingsForDestPastOnly[i] = timing[1];
- i++;
+
+
+
+ // Convert the vectors to arrays so that they can be put in the trees
+ double[][] arrayedTargetEmbeddingsFromSpikes = new double[conditioningEmbeddingsFromSpikes.size()][numDestPastIntervals + numCondPastIntervals];
+ double[][] arrayedJointEmbeddingsFromSpikes = new double[conditioningEmbeddingsFromSpikes.size()][numDestPastIntervals +
+ numCondPastIntervals + numSourcePastIntervals];
+ for (int i = 0; i < conditioningEmbeddingsFromSpikes.size(); i++) {
+ arrayedTargetEmbeddingsFromSpikes[i] = conditioningEmbeddingsFromSpikes.elementAt(i);
+ arrayedJointEmbeddingsFromSpikes[i] = jointEmbeddingsFromSpikes.elementAt(i);
}
- kdTreeDestNext = new KdTree(
- new int[] {k - 1, 1},
- new double[][][] {destPastOnlyTimings, nextTimingsForDestPastOnly});
-
- if (k == 1) {
- // We need an NN searcher for the time to next spike (dest only)
- nnSearcherDestTimeToNextSpike = new UnivariateNearestNeighbourSearcher(nextTimingsForDestPastOnly);
- } else {
- kdTreeDestHistory = new KdTree(destPastOnlyTimings);
+ double[][] arrayedTargetEmbeddingsFromSamples = new double[conditioningEmbeddingsFromSamples.size()][numDestPastIntervals + numCondPastIntervals];
+ double[][] arrayedJointEmbeddingsFromSamples = new double[conditioningEmbeddingsFromSamples.size()][numDestPastIntervals +
+ numCondPastIntervals + numSourcePastIntervals];
+ for (int i = 0; i < conditioningEmbeddingsFromSamples.size(); i++) {
+ arrayedTargetEmbeddingsFromSamples[i] = conditioningEmbeddingsFromSamples.elementAt(i);
+ arrayedJointEmbeddingsFromSamples[i] = jointEmbeddingsFromSamples.elementAt(i);
+ }
+
+ kdTreeJointAtSpikes = new KdTree(arrayedJointEmbeddingsFromSpikes);
+ kdTreeJointAtSamples = new KdTree(arrayedJointEmbeddingsFromSamples);
+ kdTreeConditioningAtSpikes = new KdTree(arrayedTargetEmbeddingsFromSpikes);
+ kdTreeConditioningAtSamples = new KdTree(arrayedTargetEmbeddingsFromSamples);
+
+ kdTreeJointAtSpikes.setNormType(normType);
+ kdTreeJointAtSamples.setNormType(normType);
+ kdTreeConditioningAtSpikes.setNormType(normType);
+ kdTreeConditioningAtSamples.setNormType(normType);
+ }
+
+ protected void makeEmbeddingsAtPoints(double[] pointsAtWhichToMakeEmbeddings, int indexOfFirstPointToUse,
+ double[] sourceSpikeTimes, double[] destSpikeTimes,
+ double[][] conditionalSpikeTimes,
+ Vector conditioningEmbeddings,
+ Vector jointEmbeddings) {
+
+ Random random = new Random();
+
+ // Initialise the starting points of all the tracking variables
+ int embeddingPointIndex = indexOfFirstPointToUse;
+ int mostRecentDestIndex = destPastIntervals[destPastIntervals.length - 1];
+ int mostRecentSourceIndex = sourcePastIntervals[sourcePastIntervals.length - 1];
+ int[] mostRecentConditioningIndices = new int[vectorOfCondPastIntervals.size()];
+ for (int i = 0; i < vectorOfCondPastIntervals.size(); i++) {
+ mostRecentConditioningIndices[i] = vectorOfCondPastIntervals.elementAt(i).length;
+ }
+
+
+ // Loop through the points at which embeddings need to be made
+ for (; embeddingPointIndex < pointsAtWhichToMakeEmbeddings.length; embeddingPointIndex++) {
+
+ // Advance the tracker of the most recent dest index
+ while (mostRecentDestIndex < (destSpikeTimes.length - 1)) {
+ if (destSpikeTimes[mostRecentDestIndex + 1] < pointsAtWhichToMakeEmbeddings[embeddingPointIndex]) {
+ mostRecentDestIndex++;
+ } else {
+ break;
+ }
+ }
+ // Do the same for the most recent source index
+ while (mostRecentSourceIndex < (sourceSpikeTimes.length - 1)) {
+ if (sourceSpikeTimes[mostRecentSourceIndex + 1] < pointsAtWhichToMakeEmbeddings[embeddingPointIndex]) {
+ mostRecentSourceIndex++;
+ } else {
+ break;
+ }
+ }
+ // Now advance the trackers for the most recent conditioning indices
+ for (int j = 0; j < mostRecentConditioningIndices.length; j++) {
+ while (mostRecentConditioningIndices[j] < (conditionalSpikeTimes[j].length - 1)) {
+ if (conditionalSpikeTimes[j][mostRecentConditioningIndices[j] + 1] < pointsAtWhichToMakeEmbeddings[embeddingPointIndex]) {
+ mostRecentConditioningIndices[j]++;
+ } else {
+ break;
+ }
+ }
+ }
+
+
+ double[] conditioningPast = new double[numDestPastIntervals + numCondPastIntervals];
+ double[] jointPast = new double[numDestPastIntervals + numCondPastIntervals + numSourcePastIntervals];
+
+ // Add the embedding intervals from the target process
+ for (int i = 0; i < destPastIntervals.length; i++) {
+ // Case where we are inserting an interval from an observation point back to the most recent event in the target process
+ if (destPastIntervals[i] == 1) {
+ conditioningPast[i] = pointsAtWhichToMakeEmbeddings[embeddingPointIndex] - destSpikeTimes[mostRecentDestIndex];
+ jointPast[i] = pointsAtWhichToMakeEmbeddings[embeddingPointIndex]
+ - destSpikeTimes[mostRecentDestIndex];
+ // Case where we are inserting an inter-event intervvl from the target process
+ } else {
+ conditioningPast[i] = destSpikeTimes[mostRecentDestIndex - destPastIntervals[i] + 2]
+ - destSpikeTimes[mostRecentDestIndex - destPastIntervals[i] + 1];
+ jointPast[i] = destSpikeTimes[mostRecentDestIndex - destPastIntervals[i] + 2]
+ - destSpikeTimes[mostRecentDestIndex - destPastIntervals[i] + 1];
+ }
+ }
+
+ // Add the embeding intervals from the conditional processes
+ int indexOfNextEmbeddingInterval = numDestPastIntervals;
+ for (int i = 0; i < vectorOfCondPastIntervals.size(); i++) {
+ for (int j = 0; j < vectorOfCondPastIntervals.elementAt(i).length; j++) {
+ // Case where we are inserting an interval from an observation point back to the most recent event in the conditioning process
+ if (vectorOfCondPastIntervals.elementAt(i)[j] == 1) {
+ conditioningPast[indexOfNextEmbeddingInterval] = pointsAtWhichToMakeEmbeddings[embeddingPointIndex]
+ - conditionalSpikeTimes[i][mostRecentConditioningIndices[i]];
+ jointPast[indexOfNextEmbeddingInterval] = pointsAtWhichToMakeEmbeddings[embeddingPointIndex]
+ - conditionalSpikeTimes[i][mostRecentConditioningIndices[i]];
+ // Case where we are inserting an inter-event interval from the conditioning process
+ } else {
+ // Convenience variable
+ int intervalNumber = vectorOfCondPastIntervals.elementAt(i)[j];
+ conditioningPast[indexOfNextEmbeddingInterval] = conditionalSpikeTimes[i][mostRecentConditioningIndices[i] - intervalNumber + 2]
+ - conditionalSpikeTimes[i][mostRecentConditioningIndices[i] - intervalNumber + 1];
+ jointPast[indexOfNextEmbeddingInterval] = conditionalSpikeTimes[i][mostRecentConditioningIndices[i] - intervalNumber + 2]
+ - conditionalSpikeTimes[i][mostRecentConditioningIndices[i] - intervalNumber + 1];
+ }
+ indexOfNextEmbeddingInterval++;
+ }
+ }
+
+ // Add the embedding intervals from the source process (this only gets added to the joint embeddings)
+ for (int i = 0; i < sourcePastIntervals.length; i++) {
+ // Case where we are inserting an interval from an observation point back to the most recent event in the source process
+ if (sourcePastIntervals[i] == 1) {
+ jointPast[indexOfNextEmbeddingInterval] = pointsAtWhichToMakeEmbeddings[embeddingPointIndex]
+ - sourceSpikeTimes[mostRecentSourceIndex];
+ // Case where we are inserting an inter-event interval from the source process
+ } else {
+ jointPast[indexOfNextEmbeddingInterval] = sourceSpikeTimes[mostRecentSourceIndex - sourcePastIntervals[i] + 2]
+ - sourceSpikeTimes[mostRecentSourceIndex - sourcePastIntervals[i] + 1];
+ }
+ indexOfNextEmbeddingInterval++;
+ }
+
+ // Add Gaussian noise, if necessary
+ if (addNoise) {
+ for (int i = 0; i < conditioningPast.length; i++) {
+ //conditioningPast[i] = Math.exp(-conditioningPast[i]);
+ conditioningPast[i] += random.nextGaussian() * noiseLevel;
+ }
+ for (int i = 0; i < jointPast.length; i++) {
+ //jointPast[i] = Math.exp(-jointPast[i]);
+ jointPast[i] += random.nextGaussian() * noiseLevel;
+ }
+ }
+
+ conditioningEmbeddings.add(conditioningPast);
+ jointEmbeddings.add(jointPast);
}
}
- protected void processEventsFromSpikingTimeSeries(double[] sourceSpikeTimes, double[] destSpikeTimes,
- int timeSeriesIndex, Vector[] eventTimings,
- Vector destPastAndNextTimings, Vector eventTypeLocator,
- Vector eventIndexLocator, Vector numEventsPerObservationSet) throws Exception {
- // addObservationsAfterParamsDetermined(sourceSpikeTimes, destSpikeTimes);
-
+
+
+ protected int getFirstDestIndex(double[] sourceSpikeTimes, double[] destSpikeTimes, double[][] conditionalSpikeTimes, Boolean setProcessTimeLengths)
+ throws Exception{
+
// First sort the spike times in case they were not properly in ascending order:
Arrays.sort(sourceSpikeTimes);
Arrays.sort(destSpikeTimes);
- // Scan to find the indices by which we have k and l spikes for dest and source
- // respectively
- int dest_index = k - 1;
- int source_index = l - 1;
- if (sourceSpikeTimes[source_index] > destSpikeTimes[dest_index]) {
- // Minimum required Source spikes are later than the dest.
- // Need to advance dest_index until it's the most recent before source_index
- for(;dest_index < destSpikeTimes.length; dest_index++) {
- if (destSpikeTimes[dest_index] > sourceSpikeTimes[source_index]) {
- // We've gone past the set of source spikes we have, we
- // can move back one in the dest series
- dest_index--;
- break;
- }
- }
- if (dest_index == destSpikeTimes.length) {
- // We didn't have enough spikes in this series to generate any observations
- // TODO work out how to handle this later -- I think this is ok
- numEventsPerObservationSet.add(0);
- return;
- // throw new Exception("Dest spikes stop before enough source spikes in time-series " + timeSeriesIndex);
- }
- } else {
- // Minimum required Dest spikes are later than the source.
- // Need to advance source_index until it's the most recent before dest_index
- for(;source_index < sourceSpikeTimes.length; source_index++) {
- if (sourceSpikeTimes[source_index] > destSpikeTimes[dest_index]) {
- // We've gone past the set of dest spikes we have, we
- // can move back one in the source series
- source_index--;
- break;
- }
- }
- if (source_index == sourceSpikeTimes.length) {
- // We didn't have enough spikes in this series to generate any observations
- // TODO work out how to handle this later -- I think this is ok
- numEventsPerObservationSet.add(0);
- return;
- // throw new Exception("Source spikes stop before enough dest spikes in time-series " + timeSeriesIndex);
+ int firstTargetIndexOfEmbedding = destPastIntervals[destPastIntervals.length - 1];
+
+ int furthestInterval = sourcePastIntervals[sourcePastIntervals.length - 1];
+ while (destSpikeTimes[firstTargetIndexOfEmbedding] <= sourceSpikeTimes[furthestInterval - 1]) {
+ firstTargetIndexOfEmbedding++;
+ }
+ if (conditionalSpikeTimes.length != vectorOfCondPastIntervals.size()) {
+ throw new Exception("Number of conditional embedding lengths does not match the number of conditional processes");
+ }
+ for (int i = 0; i < conditionalSpikeTimes.length; i++) {
+ furthestInterval = vectorOfCondPastIntervals.elementAt(i)[vectorOfCondPastIntervals.elementAt(i).length - 1];
+ while (destSpikeTimes[firstTargetIndexOfEmbedding] <= conditionalSpikeTimes[i][furthestInterval - 1]) {
+ firstTargetIndexOfEmbedding++;
}
}
- // Post-condition: dest_index and source_index are set correctly for the first set of pasts
+
+ // We don't want to reset these lengths when resampling for surrogates
+ if (setProcessTimeLengths) {
+ processTimeLengths.add(destSpikeTimes[destSpikeTimes.length - 1] - destSpikeTimes[firstTargetIndexOfEmbedding]);
+ }
- double timeToNextSpike;
- boolean nextIsDest = false;
- double[] spikeTimesForNextSpiker;
- double timeOfPrevDestSpike = destSpikeTimes[dest_index];
- int numEvents = 0;
- while(true) {
- // 0. Check whether we're finished
- if ((source_index == sourceSpikeTimes.length - 1) &&
- (dest_index == destSpikeTimes.length - 1)) {
- // We have no next spike so we can't take an observation here
- // and we're done
- break;
- }
- // Otherwise:
- // 1. Determine which of source / dest fires next
- if (source_index == sourceSpikeTimes.length - 1) {
- nextIsDest = true;
- } else if (dest_index == destSpikeTimes.length - 1) {
- nextIsDest = false;
- } else if (sourceSpikeTimes[source_index+1] < destSpikeTimes[dest_index+1]) {
- nextIsDest = false;
- } else {
- nextIsDest = true;
- }
- spikeTimesForNextSpiker = nextIsDest ? destSpikeTimes : sourceSpikeTimes;
- int indexForNextSpiker = nextIsDest ? dest_index : source_index;
- timeToNextSpike = spikeTimesForNextSpiker[indexForNextSpiker+1] - timeOfPrevDestSpike;
- // 2. Embed the past spikes
- double[] sourcePast = new double[l];
- double[] destPast = new double[k - 1];
- /* if (debug) {
- System.out.println("previousIsDest = " + previousIsDest + " and nextIsDest = " + nextIsDest);
- }*/
- sourcePast[0] = timeOfPrevDestSpike -
- sourceSpikeTimes[source_index];
- for (int i = 1; i < k; i++) {
- destPast[i - 1] = destSpikeTimes[dest_index - i + 1] -
- destSpikeTimes[dest_index - i];
- }
- for (int i = 1; i < l; i++) {
- sourcePast[i] = sourceSpikeTimes[source_index - i + 1] -
- sourceSpikeTimes[source_index - i];
- }
- // 3. Store these embedded observations
- double[][] observations = new double[][]{sourcePast, destPast,
- new double[] {timeToNextSpike}};
- if (debug) {
- System.out.printf("Adding event %d with: timeToNextSpike=%.4f, sourceSpikeTimes=", numEvents, timeToNextSpike);
- MatrixUtils.printArray(System.out, sourcePast, 3);
- System.out.printf(", destSpikeTimes=");
- MatrixUtils.printArray(System.out, destPast, 3);
- System.out.println();
- }
- // Add the index locator first so it gets the index correct before
- // we add the new event in:
- eventIndexLocator.add(eventTimings[nextIsDest ? NEXT_DEST : NEXT_SOURCE].size());
- eventTimings[nextIsDest ? NEXT_DEST : NEXT_SOURCE].add(observations);
- // TODO Switch eventTypeLocation to be of type Integer rather than int[]
- eventTypeLocator.add(nextIsDest ? NEXT_DEST : NEXT_SOURCE);
- // And finally store the observations for the dest only
- // search structure if required:
- if (nextIsDest) {
- double[][] destOnlyObservations;
- destOnlyObservations = new double[][] {
- destPast,
- new double[] {timeToNextSpike}
- };
- destPastAndNextTimings.add(destOnlyObservations);
- }
- // 4. Reset variables
- if (nextIsDest) {
- dest_index++;
- } else {
- source_index++;
- }
- timeOfPrevDestSpike = destSpikeTimes[dest_index];
- numEvents++;
- }
- numEventsPerObservationSet.add(numEvents);
- if (debug) {
- System.out.printf("Finished processing %d source-target events for observation set %d\n", numEvents, timeSeriesIndex);
+ return firstTargetIndexOfEmbedding;
+ }
+
+ protected double[] generateRandomSampleTimes(double[] sourceSpikeTimes, double[] destSpikeTimes, double[][] conditionalSpikeTimes,
+ double actualNumSamplesMultiplier, int firstTargetIndexOfEmbedding) {
+
+ double sampleLowerBound = destSpikeTimes[firstTargetIndexOfEmbedding];
+ double sampleUpperBound = destSpikeTimes[destSpikeTimes.length - 1];
+ int num_samples = (int) Math.round(actualNumSamplesMultiplier * (destSpikeTimes.length - firstTargetIndexOfEmbedding + 1));
+ double[] randomSampleTimes = new double[num_samples];
+ Random rand = new Random();
+ for (int i = 0; i < randomSampleTimes.length; i++) {
+ randomSampleTimes[i] = sampleLowerBound + rand.nextDouble() * (sampleUpperBound - sampleLowerBound);
}
+ Arrays.sort(randomSampleTimes);
+
+ return randomSampleTimes;
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getAddedMoreThanOneObservationSet()
- */
- @Override
- public boolean getAddedMoreThanOneObservationSet() {
- return (vectorOfDestinationSpikeTimes != null) &&
- (vectorOfDestinationSpikeTimes.size() > 1);
+ protected void processEventsFromSpikingTimeSeries(double[] sourceSpikeTimes, double[] destSpikeTimes, double[][] conditionalSpikeTimes,
+ Vector conditioningEmbeddingsFromSpikes, Vector jointEmbeddingsFromSpikes,
+ Vector conditioningEmbeddingsFromSamples, Vector jointEmbeddingsFromSamples,
+ double actualNumSamplesMultiplier)
+ throws Exception {
+
+ int firstTargetIndexOfEmbedding = getFirstDestIndex(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes, true);
+ double[] randomSampleTimes = generateRandomSampleTimes(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes,
+ actualNumSamplesMultiplier, firstTargetIndexOfEmbedding);
+
+ makeEmbeddingsAtPoints(destSpikeTimes, firstTargetIndexOfEmbedding, sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes,
+ conditioningEmbeddingsFromSpikes, jointEmbeddingsFromSpikes);
+ makeEmbeddingsAtPoints(randomSampleTimes, 0, sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes,
+ conditioningEmbeddingsFromSamples, jointEmbeddingsFromSamples);
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeAverageLocalOfObservations()
- */
- @Override
- public double computeAverageLocalOfObservations() throws Exception {
-
- int numberOfEvents = eventTypeLocator.size();
-
- double contributionFromSpikes = 0;
- double totalTimeLength = 0;
-
- double digammaK = MathsUtils.digamma(Knns);
- double twoInverseKTerm = 2.0 / (double) Knns;
- double inverseKTerm = 1.0 / (double) Knns;
-
- // Create temporary storage for arrays used in the neighbour counting:
- boolean[] isWithinR = new boolean[numberOfEvents]; // dummy, we don't really use this
- int[] indicesWithinR = new int[numberOfEvents];
-
- // Iterate over all the spiking events:
- Iterator eventIndexIterator = eventIndexLocator.iterator();
- int eventIndex = -1;
- int indexForNextIsDest = -1;
- for (Integer eventType : eventTypeLocator) {
- eventIndex++;
- int eventIndexWithinType = eventIndexIterator.next().intValue();
- double[][] thisEventTimings = eventTimings[eventType].elementAt(eventIndexWithinType);
- double timeToNextSpikeSincePreviousDestSpike = thisEventTimings[2][0];
- double timePreviousSourceSpikeBeforePreviousDestSpike = thisEventTimings[0][0];
- totalTimeLength += (timePreviousSourceSpikeBeforePreviousDestSpike < 0) ?
- // Source spike is after previous dest spike
- timeToNextSpikeSincePreviousDestSpike + timePreviousSourceSpikeBeforePreviousDestSpike :
- // Source spike is before previous dest spike
- timeToNextSpikeSincePreviousDestSpike;
- // Pull out the data for this observation:
- if (debug && (eventIndex < 10000)) {
- System.out.print("index = " + eventIndex + ", " +
- eventIndexWithinType + " for ->" +
- (eventType == NEXT_DEST ? "dst" : "src"));
- }
+ protected void processEventsFromSpikingTimeSeries(double[] sourceSpikeTimes, double[] destSpikeTimes, double[][] conditionalSpikeTimes,
+ Vector conditioningEmbeddingsFromSamples, Vector jointEmbeddingsFromSamples,
+ double actualNumSamplesMultiplier)
+ throws Exception {
- // Select only events where the destination spiked next:
- if (eventType == NEXT_DEST) {
- // Find the Knns nearest neighbour matches to this event,
- // with the same previous spiker and the next.
- // TODO Add dynamic exclusion time later
- PriorityQueue nnPQ =
- kdTreesJoint[eventType].findKNearestNeighbours(
- Knns, eventIndexWithinType);
-
- // Find eps_{x,y,z} as the maximum x, y and z norms amongst this set:
- double radius_sourcePast = 0.0;
- double radius_destPast = 0.0;
- double radius_destNext = 0.0;
- int radius_destNext_sampleIndex = -1;
- for (int j = 0; j < Knns; j++) {
- // Take the furthest remaining of the nearest neighbours from the PQ:
- NeighbourNodeData nnData = nnPQ.poll();
- if (nnData.norms[0] > radius_sourcePast) {
- radius_sourcePast = nnData.norms[0];
- }
- if (nnData.norms[1] > radius_destPast) {
- radius_destPast = nnData.norms[1];
- }
- if (nnData.norms[2] > radius_destNext) {
- radius_destNext = nnData.norms[2];
- radius_destNext_sampleIndex = nnData.sampleIndex;
- }
- }
-
- if (debug && (eventIndex < 10000)) {
- System.out.print(", timings: src: ");
- MatrixUtils.printArray(System.out, thisEventTimings[0], 5);
- System.out.print(", dest: ");
- MatrixUtils.printArray(System.out, thisEventTimings[1], 5);
- System.out.print(", time to next: ");
- MatrixUtils.printArray(System.out, thisEventTimings[2], 5);
- System.out.printf("index=%d: K=%d NNs at next_range %.5f (point %d)", eventIndexWithinType, Knns, radius_destNext, radius_destNext_sampleIndex);
- }
-
- indexForNextIsDest++;
-
- // Now find the matching samples in each sub-space;
- // first match dest history and source history, with a next spike in dest:
- kdTreesSourceDestHistories[NEXT_DEST].
- findPointsWithinRs(eventIndexWithinType,
- new double[] {radius_sourcePast, radius_destPast}, 0,
- true, isWithinR, indicesWithinR);
- // And check which of these samples had spike time in dest after ours:
- int countOfDestNextAndGreater = 0;
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the full joint space
- double[][] matchedHistoryEventTimings = eventTimings[NEXT_DEST].elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[2][0] > thisEventTimings[2][0] + radius_destNext) {
- // This sample had a matched history and next spike was a destination
- // spike with a longer interval than the current sample
- countOfDestNextAndGreater++;
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
- // And count how many samples with the matching history actually had a
- // *source* spike next, after ours.
- // Note that we now must go to the other kdTree for next source spike
- kdTreesSourceDestHistories[NEXT_SOURCE].
- findPointsWithinRs(
- new double[] {radius_sourcePast, radius_destPast}, thisEventTimings,
- true, isWithinR, indicesWithinR);
- // And check which of these samples had spike time in source at or after ours:
- int countOfSourceNextAndGreater = 0;
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the full joint space
- double[][] matchedHistoryEventTimings = eventTimings[NEXT_SOURCE].elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[2][0] >= thisEventTimings[2][0] - radius_destNext) {
- // This sample had a matched history and next spike was a source
- // spike with an interval longer than or considered equal to the current sample.
- // (The "equal to" is why we look for matches within radius_destNext here as well.)
- countOfSourceNextAndGreater++;
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
-
- // We need to correct radius_destNext to have a different value below, chopping it
- // where it pushes into negative times (i.e. *before* the previous spike).
- // TODO I think this is justified by the approach of Greg Ver Steeg et al. in
- // http://www.jmlr.org/proceedings/papers/v38/gao15.pdf -- it's not the same
- // as that approach, but a principled version where we know the space cannot
- // be explored (so we don't need to only apply the correction where it is large).
- double searchAreaRatio = 0; // Ratio of actual search area for full joint space compared to if the timings were independent
- double prevSource_timing_upper_original = 0, prevSource_timing_lower = 0;
- double destNext_timing_lower_original = timeToNextSpikeSincePreviousDestSpike - radius_destNext;
- double destNext_timing_upper = timeToNextSpikeSincePreviousDestSpike + radius_destNext;
- double radius_destNext_lower = radius_destNext;
- if (trimToPosNextSpikeTimes) {
- if (destNext_timing_lower_original < 0) {
- // This is the range that's actually being searched in the lower
- // dimensional space.
- // (We don't need to adjust margins in lower dimensional space for
- // this because it simply won't get any points below this!)
- destNext_timing_lower_original = 0;
- }
- // timePreviousSourceSpikeBeforePreviousDestSpike is negative if source is after dest.
- prevSource_timing_upper_original = -timePreviousSourceSpikeBeforePreviousDestSpike + radius_sourcePast;
- prevSource_timing_lower = -timePreviousSourceSpikeBeforePreviousDestSpike - radius_sourcePast;
- double destNext_timing_lower = Math.max(destNext_timing_lower_original,
- prevSource_timing_lower);
- double prevSource_timing_upper = Math.min(destNext_timing_upper, prevSource_timing_upper_original);
- // Now look at various cases for the corrective ratio:
- // TODO I think we still need to correct destNext_timing_lower_original to radius_destNext_lower here!
- double denominator = (destNext_timing_upper - destNext_timing_lower_original) *
- (prevSource_timing_upper_original - prevSource_timing_lower);
- if (destNext_timing_upper - destNext_timing_lower >
- prevSource_timing_upper - prevSource_timing_lower) {
- // Case 1:
- if (prevSource_timing_upper < destNext_timing_lower) {
- // Windows do not overlap, so there will be no correction here
- searchAreaRatio = (destNext_timing_upper - destNext_timing_lower) *
- (prevSource_timing_upper - prevSource_timing_lower);
- } else if (destNext_timing_upper < prevSource_timing_lower) {
- // Should never happen, because next spike can't be before previous source spike
- searchAreaRatio = 0;
- throw new RuntimeException("Encountered a next target spike *before* the previous source spike");
- } else if (destNext_timing_lower < prevSource_timing_lower) {
- // Dest next lower bound is below source previous lower bound -- definitely
- // can't get any points in the part below prevSource_timing_lower
- // for dest next.
- searchAreaRatio = (prevSource_timing_upper - prevSource_timing_lower) *
- (destNext_timing_upper - prevSource_timing_upper) +
- 0.5 * (prevSource_timing_upper - prevSource_timing_lower) *
- (prevSource_timing_upper - prevSource_timing_lower);
- } else {
- // destNext_timing_lower >= prevSource_timing_lower
- // Dest next lower bound is within the previous source window
- // but above the source lower bound.
- searchAreaRatio = (destNext_timing_upper - destNext_timing_lower) *
- (prevSource_timing_upper - prevSource_timing_lower) -
- 0.5 * (prevSource_timing_upper - destNext_timing_lower) *
- (prevSource_timing_upper - destNext_timing_lower);
- }
- } else {
- // Case 2:
- if (prevSource_timing_upper < destNext_timing_lower) {
- // Windows do not overlap, so there will be no correction here
- searchAreaRatio = (destNext_timing_upper - destNext_timing_lower) *
- (prevSource_timing_upper - prevSource_timing_lower);
- } else if (destNext_timing_upper < prevSource_timing_lower) {
- // Should never happen, because next spike can't be before previous source spike
- searchAreaRatio = 0;
- throw new RuntimeException("Encountered a next target spike *before* the previous source spike");
- } else if (destNext_timing_lower < prevSource_timing_lower) {
- searchAreaRatio = 0.5 * (destNext_timing_upper - prevSource_timing_lower) *
- (destNext_timing_upper - prevSource_timing_lower);
- } else {
- // destNext_timing_lower >= prevSource_timing_lower
- searchAreaRatio = (destNext_timing_lower - prevSource_timing_lower) *
- (destNext_timing_upper - destNext_timing_lower) +
- 0.5 * (destNext_timing_upper - destNext_timing_lower) *
- (destNext_timing_upper - destNext_timing_lower);
- }
- }
- if (denominator <= 0) {
- // We'll get NaN if we divide by it. No correction is necessary in this
- // case, the margins are too thin to cause a discrepancy anyway.
- searchAreaRatio = 1;
- } else {
- searchAreaRatio /= denominator;
- }
- // And finally fix up the lower radius on destNext if required:
- radius_destNext_lower = timeToNextSpikeSincePreviousDestSpike - destNext_timing_lower;
- }
-
- if (debug && (eventIndex < 10000)) {
- System.out.printf(" of %d + %d + %d points with matching S-D history",
- Knns, countOfSourceNextAndGreater, countOfDestNextAndGreater);
- }
-
- // Now find the matching samples in the dest history and
- // with a next spike timing.
- int countOfDestNextAndGreaterMatchedDest = 0;
- int countOfDestNextMatched = 0;
- if (k > 1) {
- // Search only the space of dest past -- no point
- // searching dest past and next, since we need to run through
- // all matches of dest past to count those with greater next spike
- // times we might as well count those with matching spike times
- // while we're at it.
- kdTreeDestHistory.findPointsWithinR(indexForNextIsDest, radius_destPast,
- true, isWithinR, indicesWithinR);
- // And check which of these samples had next spike time after ours:
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the dest history space
- double[][] matchedHistoryEventTimings = destPastAndNextTimings.elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[1][0] >= timeToNextSpikeSincePreviousDestSpike - radius_destNext_lower) {
- // This sample had a matched history and next spike was a
- // spike with an interval longer than or considered equal to the current sample.
- // (The "equal to" is why we look for matches within kthNnData.distance here as well.)
- countOfDestNextAndGreaterMatchedDest++;
- if (matchedHistoryEventTimings[1][0] <= timeToNextSpikeSincePreviousDestSpike + radius_destNext) {
- // Then we also have a match on the next spike itself
- countOfDestNextMatched++;
- }
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
- } else {
- // We don't take any past dest spike ISIs into account, so we just need to look at the proportion of next
- // spike times that match.
- countOfDestNextMatched = nnSearcherDestTimeToNextSpike.countPointsWithinRs(indexForNextIsDest,
- radius_destNext, radius_destNext_lower, true);
- // To check for points matching or larger, we just need to make this call with the
- // revised lower radius, because we don't check the upper one.
- countOfDestNextAndGreaterMatchedDest =
- nnSearcherDestTimeToNextSpike.countPointsWithinROrLarger(indexForNextIsDest, radius_destNext_lower, true);
- }
- double totalSearchTimeWindowCondDestPast = radius_destNext_lower + radius_destNext;
-
- if (debug && (eventIndex < 10000)) {
- System.out.printf(", and %d of %d points for D history only; ",
- countOfDestNextMatched, countOfDestNextAndGreaterMatchedDest);
- }
-
- //============================
- // This code section is if we wish to compute using digamma logs:
- //============================
- // With these neighbours counted, we're ready to compute the probability of the spike given the past
- // of source and dest.
- double logPGivenSourceAndDest;
- double logPGivenDest;
- if (k > 1) {
- // We're handling three variables:
- logPGivenSourceAndDest = digammaK - twoInverseKTerm
- - MathsUtils.digamma(Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater)
- + 1.0 / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater);
- logPGivenDest = MathsUtils.digamma(countOfDestNextMatched)
- - 1.0 / ((double) countOfDestNextMatched)
- - MathsUtils.digamma(countOfDestNextAndGreaterMatchedDest);
- } else {
- // We're really only handling two variables like an MI:
- logPGivenSourceAndDest = digammaK - inverseKTerm
- - MathsUtils.digamma(Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater);
- logPGivenDest = MathsUtils.digamma(countOfDestNextMatched)
- - MathsUtils.digamma(countOfDestNextAndGreaterMatchedDest);
- }
- if (trimToPosNextSpikeTimes) {
- logPGivenSourceAndDest -= Math.log(searchAreaRatio);
- }
- //============================
-
- if (debug && (eventIndex < 10000)) {
- System.out.printf(" te ~~ log (%d/%d)/(%d/%d) = %.4f wc-> %.5f bc-> %.4f (inferred rates %.4f vs %.4f, " +
- "win-cor %.5f vs %.5f, bias-corrected %.5f vs %.5f)\n", Knns,
- Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater,
- countOfDestNextMatched, countOfDestNextAndGreaterMatchedDest,
- Math.log(((double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater)) /
- ((double) (countOfDestNextMatched) / (double) (countOfDestNextAndGreaterMatchedDest))),
- // TE from Window corrected rates:
- Math.log((((double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater)) /
- ((destNext_timing_upper - destNext_timing_lower_original)*searchAreaRatio)) /
- (((double) (countOfDestNextMatched) / (double) (countOfDestNextAndGreaterMatchedDest)) /
- totalSearchTimeWindowCondDestPast)),
- // TE from Bias corrected rates:
- logPGivenSourceAndDest - logPGivenDest,
- // Inferred rates raw:
- (double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater) / (2.0*radius_destNext),
- (double) (countOfDestNextMatched) / (double) (countOfDestNextAndGreaterMatchedDest) / (2.0*radius_destNext),
- // Inferred rates with window correction:
- ((double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater)) /
- ((destNext_timing_upper - destNext_timing_lower_original)*searchAreaRatio),
- ((double) (countOfDestNextMatched) / (double) (countOfDestNextAndGreaterMatchedDest)) /
- totalSearchTimeWindowCondDestPast,
- // Transform the log likelihoods into bias corrected rates:
- Math.exp(logPGivenSourceAndDest) / (2.0*radius_destNext), Math.exp(logPGivenDest) / (2.0*radius_destNext));
- if (trimToPosNextSpikeTimes) {
- System.out.printf("Search area ratio: %.5f, correction %.5f, t_y_upper %.5f, t_y_lower %.5f\n",
- searchAreaRatio, -Math.log(searchAreaRatio), prevSource_timing_upper_original, prevSource_timing_lower);
- }
- }
- // Unexplained case:
- if (countOfDestNextMatched < Knns) {
- // TODO Should not happen, print something!
- System.out.printf("SHOULD NOT HAPPEN!\n");
- // So debug this:
- nnPQ = kdTreesJoint[eventType].findKNearestNeighbours(
- Knns, eventIndexWithinType);
- for (int j = 0; j < Knns; j++) {
- // Take the furthest remaining of the nearest neighbours from the PQ:
- NeighbourNodeData nnData = nnPQ.poll();
- System.out.printf("NN data %d norms: source %.5f, destPast %.5f, destNext %.5f, nextIndex %d\n",
- nnData.norms[0], nnData.norms[1], nnData.norms[2], nnData.sampleIndex);
- }
-
- }
-
- //======================
- // Add the contribution in:
- // a.If we were using digamma logs:
- contributionFromSpikes += logPGivenSourceAndDest - logPGivenDest;
- // b. If we are only using window corrections but actual ratios:
- //contributionFromSpikes +=
- // Math.log((((double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater)) /
- // ((destNext_timing_upper - destNext_timing_lower_original)*searchAreaRatio)) /
- // (((double) (countOfDestNextMatched) / (double) (countOfDestNextAndGreaterMatchedDest)) /
- // totalSearchTimeWindowCondDestPast));
- } else {
- if (debug && (eventIndex < 10000)) {
- System.out.println();
- }
- }
- }
- System.out.println("All done!");
- contributionFromSpikes /= totalTimeLength;
- return contributionFromSpikes;
+ int firstTargetIndexOfEmbedding = getFirstDestIndex(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes, false);
+ double[] randomSampleTimes = generateRandomSampleTimes(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes,
+ actualNumSamplesMultiplier, firstTargetIndexOfEmbedding);
+ makeEmbeddingsAtPoints(randomSampleTimes, 0, sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes,
+ conditioningEmbeddingsFromSamples, jointEmbeddingsFromSamples);
}
/*
- * This old method is not adjusted for the newer representation yet
- *
- public double computeAverageLocalOfObservationsAlg1() throws Exception {
-
- int numberOfEvents = eventTypeLocator.size();
-
- double te = 0;
- double contributionFromSpikes = 0;
- double contributionFromNonSpikes = 0;
- double contributionFromNonSpikes_destOnly = 0;
- double contributionFromNonSpikes_destAndSource = 0;
- double totalTimeLength = 0;
-
- // Create temporary storage for arrays used in the neighbour counting:
- boolean[] isWithinR = new boolean[numberOfEvents]; // dummy, we don't really use this
- int[] indicesWithinR = new int[numberOfEvents];
-
- // Iterate over all the spiking events:
- Iterator eventIndexIterator = eventIndexLocator.iterator();
- int eventIndex = -1;
- int indexForNextIsDest = -1;
- for (int[] eventType : eventTypeLocator) {
- eventIndex++;
- int eventIndexWithinType = eventIndexIterator.next().intValue();
- double[][] thisEventTimings = eventTimings[eventType[0]][eventType[1]].elementAt(eventIndexWithinType);
- totalTimeLength += thisEventTimings[2][0];
- // Find the Knns nearest neighbour matches to this event,
- // with the same previous spiker and the next.
- // TODO Add dynamic exclusion time later
- PriorityQueue nnPQ =
- kdTreesJoint[eventType[0]][eventType[1]].findKNearestNeighbours(
- Knns, eventIndexWithinType);
- // First element in the PQ is the kth NN,
- // and epsilon = kthNnData.distance
- NeighbourNodeData kthNnData = nnPQ.poll();
- double radiusToKnn = kthNnData.distance;
- if (debug && (eventIndex < 10000)) {
- // Pull out the data for this observation:
- System.out.print("index = " + eventIndex + ", " +
- eventIndexWithinType + " for " +
- (eventType[0] == PREV_DEST ? "dst" : "src") +
- "->" +
- (eventType[1] == NEXT_DEST ? "dst" : "src") +
- ", timings: src: ");
- MatrixUtils.printArray(System.out, thisEventTimings[0], 3);
- System.out.print(", dest: ");
- MatrixUtils.printArray(System.out, thisEventTimings[1], 3);
- System.out.print(", time to next: ");
- MatrixUtils.printArray(System.out, thisEventTimings[2], 3);
- System.out.printf("index=%d: K=%d NNs at range %.5f (point %d)", eventIndexWithinType, Knns, radiusToKnn, kthNnData.sampleIndex);
- }
-
- // Select only events where the destination spiked next:
- if (eventType[1] == NEXT_DEST) {
- indexForNextIsDest++;
-
- // Now find the matching samples in each sub-space;
- // first match dest history and source history, with a next spike in dest:
- kdTreesSourceDestHistories[eventType[0]][NEXT_DEST].
- findPointsWithinR(eventIndexWithinType, radiusToKnn, 0,
- false, isWithinR, indicesWithinR);
- // And check which of these samples had spike time in dest after ours:
- int countOfDestNextAndGreater = 0;
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the full joint space
- double[][] matchedHistoryEventTimings = eventTimings[eventType[0]][NEXT_DEST].elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[2][0] >= thisEventTimings[2][0] + radiusToKnn) {
- // This sample had a matched history and next spike was a destination
- // spike with a longer interval than the current sample
- countOfDestNextAndGreater++;
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
- // And count how many samples with the matching history actually had a
- // *source* spike next, after ours.
- // Note that we now must go to the other kdTree for next source spike
- kdTreesSourceDestHistories[eventType[0]][NEXT_SOURCE].
- findPointsWithinR(radiusToKnn, thisEventTimings,
- false, isWithinR, indicesWithinR);
- // And check which of these samples had spike time in source after ours:
- int countOfSourceNextAndGreater = 0;
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the full joint space
- double[][] matchedHistoryEventTimings = eventTimings[eventType[0]][NEXT_SOURCE].elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[2][0] > thisEventTimings[2][0] - radiusToKnn) {
- // This sample had a matched history and next spike was a source
- // spike with an interval longer than or considered equal to the current sample.
- // (The "equal to" is why we look for matches within kthNnData.distance here as well.)
- countOfSourceNextAndGreater++;
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
-
- if (debug && (eventIndex < 10000)) {
- System.out.printf(" of %d + %d + %d points with matching S-D history",
- Knns, countOfSourceNextAndGreater, countOfDestNextAndGreater);
- }
-
- // Now find the matching samples in the dest history and
- // with a next spike timing.
- // Construct the appropriate timings to compare to here:
- double[][] destOnlyObservations;
- double[][] destPastOnlyObservations;
- double timeToNextSpikeSincePreviousDestSpike;
- if (eventType[0] == PREV_DEST) {
- destOnlyObservations = new double[][] {
- thisEventTimings[1], // timing of past dest spikes
- thisEventTimings[2] // time to next spike
- };
- timeToNextSpikeSincePreviousDestSpike = thisEventTimings[2][0];
- destPastOnlyObservations = new double[][] {
- thisEventTimings[1] // timing of past dest spikes
- };
- } else {
- // previous is source:
- // We can take a copy of the dest past timings, removing the first entry
- // (since this only signals time the dest last fired before the source)
- // and add that entry to the timeToNextSpike (which was back to the
- // source firing).
- double[] destPastOnly = Arrays.copyOfRange(
- thisEventTimings[1], 1, thisEventTimings[1].length);
- timeToNextSpikeSincePreviousDestSpike =
- thisEventTimings[1][0] + thisEventTimings[2][0];
- destOnlyObservations = new double[][] {
- destPastOnly,
- new double[] {timeToNextSpikeSincePreviousDestSpike}
- };
- destPastOnlyObservations = new double[][] {
- destPastOnly
- };
- }
- int countOfDestNextAndGreaterMatchedDest = 0;
- int countOfDestNextMatched = 0;
- if (k > 1) {
- // Search only the space of dest past -- no point
- // searching dest past and next, since we need to run through
- // all matches of dest past to count those with greater next spike
- // times we might as well count those with matching spike times
- // while we're at it.
-
- // OLD WAY:
- // NO NO NO -- Can't search for it this way, because it's biased --
- // should search for it by giving the index of this dest past-next
- // observation, so that it doesn't match to this observation.
- // Should be able to use indexForNextIsDest here
- //kdTreeDestHistory.findPointsWithinR(radiusToKnn, destPastOnlyObservations,
- // false, isWithinR, indicesWithinR);
- // Proper way:
- kdTreeDestHistory.findPointsWithinR(indexForNextIsDest, radiusToKnn,
- false, isWithinR, indicesWithinR);
- // And check which of these samples had next spike time after ours:
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the dest history space
- double[][] matchedHistoryEventTimings = destPastAndNextTimings.elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[1][0] > timeToNextSpikeSincePreviousDestSpike - radiusToKnn) {
- // This sample had a matched history and next spike was a
- // spike with an interval longer than or considered equal to the current sample.
- // (The "equal to" is why we look for matches within kthNnData.distance here as well.)
- countOfDestNextAndGreaterMatchedDest++;
- if (matchedHistoryEventTimings[1][0] < timeToNextSpikeSincePreviousDestSpike + radiusToKnn) {
- // Then we also have a match on the next spike itself
- countOfDestNextMatched++;
- }
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
- } else {
- // We don't take any past dest spike times into account, so we just need to look at the proportion of next
- // spike times that match.
- countOfDestNextMatched = nnSearcherDestTimeToNextSpike.countPointsStrictlyWithinR(indexForNextIsDest, radiusToKnn);
- countOfDestNextAndGreaterMatchedDest = countOfDestNextMatched +
- nnSearcherDestTimeToNextSpike.countPointsWithinROrLarger(indexForNextIsDest, radiusToKnn, false);
- }
-
- if (debug && (eventIndex < 10000)) {
- System.out.printf(", and %d of %d points for D history only; ",
- countOfDestNextMatched, countOfDestNextAndGreaterMatchedDest);
- }
+ * (non-Javadoc)
+ *
+ * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#
+ * getAddedMoreThanOneObservationSet()
+ */
+ @Override
+ public boolean getAddedMoreThanOneObservationSet() {
+ return (vectorOfDestinationSpikeTimes != null) && (vectorOfDestinationSpikeTimes.size() > 1);
+ }
- // With these neighbours counted, we're ready to compute the probability of the spike given the past
- // of source and dest.
- // Digammas for algorithm 1 include the extra "+1" on all terms except
- // for the full joint space
- double logPGivenSourceAndDest = MathsUtils.digamma(Knns) -
- MathsUtils.digamma(Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater + 1);
- double logPGivenDest = MathsUtils.digamma(countOfDestNextMatched + 1) -
- MathsUtils.digamma(countOfDestNextAndGreaterMatchedDest + 1);
- if (debug && (eventIndex < 10000)) {
- System.out.printf(" te ~~ log (%d/%d)/(%d/%d) = %.4f -> %.4f (inferred rates %.4f vs %.4f)\n", Knns,
- Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater + 1,
- countOfDestNextMatched + 1, countOfDestNextAndGreaterMatchedDest + 1,
- Math.log(((double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater + 1)) /
- ((double) (countOfDestNextMatched + 1) / (double) (countOfDestNextAndGreaterMatchedDest + 1))),
- logPGivenSourceAndDest - logPGivenDest,
- (double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater + 1) / (2.0*radiusToKnn),
- (double) (countOfDestNextMatched + 1) / (double) (countOfDestNextAndGreaterMatchedDest + 1) / (2.0*radiusToKnn));
- }
- contributionFromSpikes += logPGivenSourceAndDest - logPGivenDest;
- } else {
- if (debug) {
- System.out.println();
- }
- }
-
- // Regardless of which type of event it was, we need to integrate
- // the spiking rates up until the next spiking event
- // Our first attempt at a solution uses the search width defined
- // using the history and the next spike.
-
- // Consider first the destination process only.
- // Match dest history
- double[][] destPastOnlyObservations;
- double timeToNextSpikeSincePreviousDestSpike;
- if (eventType[0] == PREV_DEST) {
- timeToNextSpikeSincePreviousDestSpike = thisEventTimings[2][0];
- destPastOnlyObservations = new double[][] {
- thisEventTimings[1] // timing of past dest spikes
- };
- } else {
- // previous is source:
- // We can take a copy of the dest past timings, removing the first entry
- // (since this only signals time the dest last fired before the source)
- // and add that entry to the timeToNextSpike (which was back to the
- // source firing).
- double[] destPastOnly = Arrays.copyOfRange(
- thisEventTimings[1], 1, thisEventTimings[1].length);
- timeToNextSpikeSincePreviousDestSpike =
- thisEventTimings[1][0] + thisEventTimings[2][0];
- destPastOnlyObservations = new double[][] {
- destPastOnly
- };
- }
- int countOfDestNextEarlier = 0;
- int countOfDestMatches = 0;
- if (k > 1) {
- kdTreeDestHistory.findPointsWithinR(radiusToKnn, destPastOnlyObservations,
- false, isWithinR, indicesWithinR);
- // And check which of these samples had next spike time before ours:
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the dest history space
- double[][] matchedHistoryEventTimings = destPastAndNextTimings.elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[1][0] < timeToNextSpikeSincePreviousDestSpike) {
- // This sample had a matched history and next spike was a
- // spike with an interval shorter than the current sample.
- countOfDestNextEarlier++;
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- countOfDestMatches++;
- }
- } else {
- // We're not using the past, so we match on everything up to the last spike
- countOfDestMatches = nnSearcherDestTimeToNextSpike.getNumObservations();
- countOfDestNextEarlier = nnSearcherDestTimeToNextSpike.countPointsSmallerAndOutsideR(
- // indexForNextIsDest must point to the next event (possibly this one) where the dest spikes next.
- (eventType[1] == NEXT_DEST) ? indexForNextIsDest : indexForNextIsDest + 1,
- radiusToKnn, false);
+ // Class to allow returning two values in the subsequent method
+ private static class distanceAndNumPoints {
+ public double distance;
+ public int numPoints;
- }
- // And include the contribution for each of these
- double integralForDestHistorySpace = 0;
- for (int hi = 0; hi < countOfDestNextEarlier; hi++) {
- integralForDestHistorySpace += (double) 1 /
- (double) (countOfDestMatches - hi);
- }
-
- // First match dest history and source history, with a next spike in dest:
- kdTreesSourceDestHistories[eventType[0]][NEXT_DEST].
- findPointsWithinR(radiusToKnn, thisEventTimings,
- false, isWithinR, indicesWithinR);
- // And store which of these samples had spike time in dest before ours.
- // Store them in a vector of double arrays, with each array holding the
- // spike time then -1 for a next dest spike and +1 for a source spike
- Vector spikesBeforeOurs = new Vector();
- int countOfSpikesAfterAndIncludingOurs = 0;
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the full joint space
- double[][] matchedHistoryEventTimings = eventTimings[eventType[0]][NEXT_DEST].elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[2][0] < thisEventTimings[2][0]) {
- // This sample had a matched history and next spike was a destination
- // spike with a shorted interval than the current sample
- spikesBeforeOurs.add(new double[] {
- matchedHistoryEventTimings[2][0], -1});
- } else {
- countOfSpikesAfterAndIncludingOurs++;
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
- // And store which of these samples had spike time in source before ours.
- // Store them in a vector of double arrays, with each array holding the
- // spike time then -1 for a next dest spike and +1 for a source spike
- // Note that we now must go to the other kdTree for next source spike
- kdTreesSourceDestHistories[eventType[0]][NEXT_SOURCE].
- findPointsWithinR(radiusToKnn, thisEventTimings,
- false, isWithinR, indicesWithinR);
- for (int nIndex = 0; indicesWithinR[nIndex] != -1; nIndex++) {
- // Pull out this matching event from the full joint space
- double[][] matchedHistoryEventTimings = eventTimings[eventType[0]][NEXT_SOURCE].elementAt(indicesWithinR[nIndex]);
- if (matchedHistoryEventTimings[2][0] < thisEventTimings[2][0]) {
- // This sample had a matched history and next spike was a source
- // spike with a shorted interval than the current sample
- spikesBeforeOurs.add(new double[] {
- matchedHistoryEventTimings[2][0], +1});
- } else {
- countOfSpikesAfterAndIncludingOurs++;
- }
- // Reset the isWithinR array while we're here
- isWithinR[indicesWithinR[nIndex]] = false;
- }
- // Now we can sort the spikes which occur before ours and process
- // them in order:
- // Next line doesn't work, so replaced with clunkier code:
- // double[][] nextSpikeTimesAndType = (double[][]) spikesBeforeOurs.toArray();
- double[][] nextSpikeTimesAndType = new double[spikesBeforeOurs.size()][];
- for (int si = 0; si < nextSpikeTimesAndType.length; si++) {
- nextSpikeTimesAndType[si] = spikesBeforeOurs.elementAt(si);
- }
- Arrays.sort(nextSpikeTimesAndType, FirstIndexComparatorDouble.getInstance());
- double integralForJointSpace = 0;
- for (int si = 0; si < nextSpikeTimesAndType.length; si++) {
- if (nextSpikeTimesAndType[si][1] < 0) {
- // We have a next spike from the dest, which is
- // earlier than our spike.
- // Integrated Prob for getting a spike here is 1 / N, where
- // N is the number of properly matched histories (i.e.
- // which don't have a next spike before this one)
- double intPNext = (double) 1 / (double)
- (nextSpikeTimesAndType.length - si + countOfSpikesAfterAndIncludingOurs);
- integralForJointSpace += intPNext;
- }
- // Ignore next spikes on the source, they simply get removed
- // from the matched histories count
- }
-
- // We now have the integral of spike rates given the dest and
- // joint histories, so subtract this out:
- contributionFromNonSpikes += integralForDestHistorySpace - integralForJointSpace;
- contributionFromNonSpikes_destAndSource += integralForJointSpace;
- contributionFromNonSpikes_destOnly += integralForDestHistorySpace;
+ public distanceAndNumPoints(double distance, int numPoints) {
+ this.distance = distance;
+ this.numPoints = numPoints;
}
- contributionFromSpikes /= totalTimeLength;
- contributionFromNonSpikes /= totalTimeLength;
- contributionFromNonSpikes_destAndSource /= totalTimeLength;
- contributionFromNonSpikes_destOnly /= totalTimeLength;
- te = contributionFromSpikes + contributionFromNonSpikes;
- System.out.printf("TE = %.4f (spikes) + %.4f (non-spikes: d:%.4f - s-d:%.4f) = %.4f\n",
- contributionFromSpikes, contributionFromNonSpikes,
- contributionFromNonSpikes_destOnly, contributionFromNonSpikes_destAndSource, te);
- return te;
}
- */
+
+ private distanceAndNumPoints findMaxDistanceAndNumPointsFromIndices(double[] point, int[] indices, Vector setOfPoints) {
+ double maxDistance = 0;
+ int i = 0;
+ for (; indices[i] != -1; i++) {
+ double distance = KdTree.norm(point, setOfPoints.elementAt(indices[i]), normType);
+ if (distance > maxDistance) {
+ maxDistance = distance;
+ }
+ }
+ return new distanceAndNumPoints(maxDistance, i);
+ }
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeLocalOfPreviousObservations()
+ @Override
+ public double computeAverageLocalOfObservations() throws Exception {
+ return computeAverageLocalOfObservations(kdTreeJointAtSpikes, jointEmbeddingsFromSpikes);
+ }
+
+
+ /*
+ * We take the actual joint tree at spikes (along with the associated embeddings) as an argument, as we will need to swap these out when
+ * computing surrogates.
+ */
+ public double computeAverageLocalOfObservations(KdTree actualKdTreeJointAtSpikes, Vector actualJointEmbeddingsFromSpikes) throws Exception {
+
+ double currentSum = 0;
+ for (int i = 0; i < conditioningEmbeddingsFromSpikes.size(); i++) {
+
+ double radiusJointSpikes = actualKdTreeJointAtSpikes.findKNearestNeighbours(Knns, i).poll().norms[0];
+ double radiusJointSamples = kdTreeJointAtSamples.findKNearestNeighbours(Knns,
+ new double[][] { actualJointEmbeddingsFromSpikes.elementAt(i) }).poll().norms[0];
+
+ /*
+ The algorithm specified in box 1 of doi.org/10.1371/journal.pcbi.1008054 specifies finding the maximum of the two radii
+ just calculated and then redoing the searches in both sets at this radius. In this implementation, however, we make use
+ of the fact that one radius is equal to the maximum, and so only one search needs to be redone.
+ */
+ double eps = 0.0;
+ // Need variables for the number of neighbours as this is now variable within the maximum radius
+ int kJointSpikes = 0;
+ int kJointSamples = 0;
+ if (radiusJointSpikes >= radiusJointSamples) {
+ /*
+ The maximum was the radius in the set of embeddings at spikes, so redo search in the set of embeddings at randomly
+ sampled points, using this larger radius.
+ */
+ kJointSpikes = Knns;
+ int[] indicesWithinR = new int[jointEmbeddingsFromSamples.size()];
+ boolean[] isWithinR = new boolean[jointEmbeddingsFromSamples.size()];
+ kdTreeJointAtSamples.findPointsWithinR(radiusJointSpikes + eps,
+ new double[][] { actualJointEmbeddingsFromSpikes.elementAt(i) },
+ true,
+ isWithinR,
+ indicesWithinR);
+ distanceAndNumPoints temp = findMaxDistanceAndNumPointsFromIndices(actualJointEmbeddingsFromSpikes.elementAt(i), indicesWithinR,
+ jointEmbeddingsFromSamples);
+ kJointSamples = temp.numPoints;
+ radiusJointSamples = temp.distance;
+ } else {
+ /*
+ The maximum was the radius in the set of embeddings at randomly sampled points, so redo search in the set of embeddings
+ at spikes, using this larger radius.
+ */
+ kJointSamples = Knns;
+ int[] indicesWithinR = new int[jointEmbeddingsFromSamples.size()];
+ boolean[] isWithinR = new boolean[jointEmbeddingsFromSamples.size()];
+ actualKdTreeJointAtSpikes.findPointsWithinR(radiusJointSamples + eps,
+ new double[][] { actualJointEmbeddingsFromSpikes.elementAt(i) },
+ true,
+ isWithinR,
+ indicesWithinR);
+ distanceAndNumPoints temp = findMaxDistanceAndNumPointsFromIndices(actualJointEmbeddingsFromSpikes.elementAt(i), indicesWithinR,
+ actualJointEmbeddingsFromSpikes);
+ // -1 due to the point itself being in the set
+ kJointSpikes = temp.numPoints - 1;
+ radiusJointSpikes = temp.distance;
+ }
+
+ // Repeat the above steps, but in the conditioning (rather than joint) space.
+ double radiusConditioningSpikes = kdTreeConditioningAtSpikes.findKNearestNeighbours(Knns, i).poll().norms[0];
+ double radiusConditioningSamples = kdTreeConditioningAtSamples.findKNearestNeighbours(Knns,
+ new double[][] { conditioningEmbeddingsFromSpikes.elementAt(i) }).poll().norms[0];
+ int kConditioningSpikes = 0;
+ int kConditioningSamples = 0;
+ if (radiusConditioningSpikes >= radiusConditioningSamples) {
+ kConditioningSpikes = Knns;
+ int[] indicesWithinR = new int[conditioningEmbeddingsFromSamples.size()];
+ boolean[] isWithinR = new boolean[conditioningEmbeddingsFromSamples.size()];
+ kdTreeConditioningAtSamples.findPointsWithinR(radiusConditioningSpikes + eps,
+ new double[][] { conditioningEmbeddingsFromSpikes.elementAt(i) },
+ true,
+ isWithinR,
+ indicesWithinR);
+ distanceAndNumPoints temp = findMaxDistanceAndNumPointsFromIndices(conditioningEmbeddingsFromSpikes.elementAt(i), indicesWithinR,
+ conditioningEmbeddingsFromSamples);
+ kConditioningSamples = temp.numPoints;
+ radiusConditioningSamples = temp.distance;
+ } else {
+ kConditioningSamples = Knns;
+ int[] indicesWithinR = new int[conditioningEmbeddingsFromSamples.size()];
+ boolean[] isWithinR = new boolean[conditioningEmbeddingsFromSamples.size()];
+ kdTreeConditioningAtSpikes.findPointsWithinR(radiusConditioningSamples + eps,
+ new double[][] { conditioningEmbeddingsFromSpikes.elementAt(i) },
+ true,
+ isWithinR,
+ indicesWithinR);
+ distanceAndNumPoints temp = findMaxDistanceAndNumPointsFromIndices(conditioningEmbeddingsFromSpikes.elementAt(i), indicesWithinR,
+ conditioningEmbeddingsFromSpikes);
+ // -1 due to the point itself being in the set
+ kConditioningSpikes = temp.numPoints - 1;
+ radiusConditioningSpikes = temp.distance;
+ }
+
+ /*
+ * TODO
+ * The KdTree class defaults to the squared euclidean distance when the euclidean norm is specified. This is fine for Kraskov estimators
+ * (as the radii are never used, just the numbers of points within radii). It causes problems here though, as we do use the radii and the
+ * squared euclidean distance is not a distance metric. We get around this by just taking the square root here, but it might be better to
+ * fix this in the KdTree class.
+ */
+ if (normType == EuclideanUtils.NORM_EUCLIDEAN) {
+ radiusJointSpikes = Math.sqrt(radiusJointSpikes);
+ radiusJointSamples = Math.sqrt(radiusJointSamples);
+ radiusConditioningSpikes = Math.sqrt(radiusConditioningSpikes);
+ radiusConditioningSamples = Math.sqrt(radiusConditioningSamples);
+ }
+
+ currentSum += (MathsUtils.digamma(kJointSpikes) - MathsUtils.digamma(kJointSamples) +
+ ((numDestPastIntervals + numCondPastIntervals + numSourcePastIntervals) * (-Math.log(radiusJointSpikes) + Math.log(radiusJointSamples))) -
+ MathsUtils.digamma(kConditioningSpikes) + MathsUtils.digamma(kConditioningSamples) +
+ + ((numDestPastIntervals + numCondPastIntervals) * (Math.log(radiusConditioningSpikes) - Math.log(radiusConditioningSamples))));
+ if (Double.isNaN(currentSum)) {
+ throw new Exception("NaNs in TE clac");
+ }
+ }
+ // Normalise by time
+ double timeSum = 0;
+ for (Double time : processTimeLengths) {
+ timeSum += time;
+ }
+ currentSum /= timeSum;
+ return currentSum;
+ }
+
+
+ @Override
+ public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck, double estimatedValue) throws Exception {
+ return computeSignificance(numPermutationsToCheck,
+ estimatedValue,
+ System.currentTimeMillis());
+ }
+
+ @Override
+ public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck,
+ double estimatedValue, long randomSeed) throws Exception{
+
+ Random random = new Random(randomSeed);
+ double[] surrogateTEValues = new double[numPermutationsToCheck];
+
+ for (int permutationNumber = 0; permutationNumber < numPermutationsToCheck; permutationNumber++) {
+ Vector resampledConditioningEmbeddingsFromSamples = new Vector();
+ Vector resampledJointEmbeddingsFromSamples = new Vector();
+
+ // Send all of the observations through:
+ Iterator sourceIterator = vectorOfSourceSpikeTimes.iterator();
+ Iterator conditionalIterator = null;
+ if (vectorOfConditionalSpikeTimes.size() > 0) {
+ conditionalIterator = vectorOfConditionalSpikeTimes.iterator();
+ }
+ int timeSeriesIndex = 0;
+ for (double[] destSpikeTimes : vectorOfDestinationSpikeTimes) {
+ double[] sourceSpikeTimes = sourceIterator.next();
+ double[][] conditionalSpikeTimes = null;
+ if (vectorOfConditionalSpikeTimes.size() > 0) {
+ conditionalSpikeTimes = conditionalIterator.next();
+ } else {
+ conditionalSpikeTimes = new double[][] {};
+ }
+ processEventsFromSpikingTimeSeries(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes,
+ resampledConditioningEmbeddingsFromSamples, resampledJointEmbeddingsFromSamples,
+ surrogateNumSamplesMultiplier);
+ }
+ // Convert the vectors to arrays so that they can be put in the trees
+ double[][] arrayedResampledConditioningEmbeddingsFromSamples = new double[resampledConditioningEmbeddingsFromSamples.size()][numDestPastIntervals];
+ for (int i = 0; i < resampledConditioningEmbeddingsFromSamples.size(); i++) {
+ arrayedResampledConditioningEmbeddingsFromSamples[i] = resampledConditioningEmbeddingsFromSamples.elementAt(i);
+ }
+
+ KdTree resampledKdTreeConditioningAtSamples = new KdTree(arrayedResampledConditioningEmbeddingsFromSamples);
+ resampledKdTreeConditioningAtSamples.setNormType(normType);
+
+ Vector conditionallyPermutedJointEmbeddingsFromSpikes = new Vector(jointEmbeddingsFromSpikes);
+
+ Vector usedIndices = new Vector();
+ for (int i = 0; i < conditionallyPermutedJointEmbeddingsFromSpikes.size(); i++) {
+ PriorityQueue neighbours =
+ resampledKdTreeConditioningAtSamples.findKNearestNeighbours(kPerm,
+ new double[][] {conditioningEmbeddingsFromSpikes.elementAt(i)});
+ ArrayList foundIndices = new ArrayList();
+ for (int j = 0; j < kPerm; j++) {
+ foundIndices.add(neighbours.poll().sampleIndex);
+ }
+
+ ArrayList prunedIndices = new ArrayList(foundIndices);
+ prunedIndices.removeAll(usedIndices);
+ int chosenIndex = 0;
+ if (prunedIndices.size() > 0) {
+ chosenIndex = prunedIndices.get(random.nextInt(prunedIndices.size()));
+ } else {
+ chosenIndex = foundIndices.get(random.nextInt(foundIndices.size()));
+ }
+ usedIndices.add(chosenIndex);
+ int embeddingLength = conditionallyPermutedJointEmbeddingsFromSpikes.elementAt(i).length;
+ for(int j = 0; j < numSourcePastIntervals; j++) {
+ conditionallyPermutedJointEmbeddingsFromSpikes.elementAt(i)[embeddingLength - numSourcePastIntervals + j] =
+ resampledJointEmbeddingsFromSamples.elementAt(chosenIndex)[embeddingLength - numSourcePastIntervals + j];
+ }
+
+ }
+
+ double[][] arrayedConditionallyPermutedJointEmbeddingsFromSpikes = new double[conditionallyPermutedJointEmbeddingsFromSpikes.size()][];
+ for (int i = 0; i < conditionallyPermutedJointEmbeddingsFromSpikes.size(); i++) {
+ arrayedConditionallyPermutedJointEmbeddingsFromSpikes[i] = conditionallyPermutedJointEmbeddingsFromSpikes.elementAt(i);
+ }
+ KdTree conditionallyPermutedKdTreeJointFromSpikes = new KdTree(arrayedConditionallyPermutedJointEmbeddingsFromSpikes);
+ conditionallyPermutedKdTreeJointFromSpikes.setNormType(normType);
+ surrogateTEValues[permutationNumber] = computeAverageLocalOfObservations(conditionallyPermutedKdTreeJointFromSpikes,
+ conditionallyPermutedJointEmbeddingsFromSpikes);
+ }
+ return new EmpiricalMeasurementDistribution(surrogateTEValues, estimatedValue);
+ }
+
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#
+ * computeLocalOfPreviousObservations()
*/
@Override
- public SpikingLocalInformationValues computeLocalOfPreviousObservations()
- throws Exception {
+ public SpikingLocalInformationValues computeLocalOfPreviousObservations() throws Exception {
// TODO Auto-generated method stub
return null;
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeSignificance(int)
- */
- @Override
- public EmpiricalMeasurementDistribution computeSignificance(
- int numPermutationsToCheck) throws Exception {
- // TODO Auto-generated method stub
- return null;
- }
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeSignificance(int[][])
- */
- @Override
- public EmpiricalMeasurementDistribution computeSignificance(
- int[][] newOrderings) throws Exception {
- // TODO Auto-generated method stub
- return null;
- }
-
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#setDebug(boolean)
+ /*
+ * (non-Javadoc)
+ *
+ * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#setDebug(
+ * boolean)
*/
@Override
public void setDebug(boolean debug) {
this.debug = debug;
}
- /* (non-Javadoc)
- * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getLastAverage()
+ /*
+ * (non-Javadoc)
+ *
+ * @see
+ * infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getLastAverage
+ * ()
*/
@Override
public double getLastAverage() {
// TODO Auto-generated method stub
return 0;
}
-
}
diff --git a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegrationWindows2.java b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegrationWindows2.java
index 2349b8a..64f25d6 100644
--- a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegrationWindows2.java
+++ b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegrationWindows2.java
@@ -723,33 +723,33 @@ public class TransferEntropyCalculatorSpikingIntegrationWindows2 implements
//kdTreeDestNext.setNormType(3);
- nnPQ =
- kdTreeDestNext.findKNearestNeighbours(
- Knns, eventIndexWithinType);
+ // nnPQ =
+ // kdTreeDestNext.findKNearestNeighbours(
+ // Knns, eventIndexWithinType);
- // Find eps_{x,y,z} as the maximum x, y and z norms amongst this set:
- radius_sourcePast = 0.0;
- radius_destPast = 0.0;
- radius_destNext = 0.0;
- radius_destNext_sampleIndex = -1;
- radius_max = 0;
- for (int j = 0; j < Knns; j++) {
- // Take the furthest remaining of the nearest neighbours from the PQ:
- NeighbourNodeData nnData = nnPQ.poll();
- /*if (nnData.norms[0] > radius_sourcePast) {
- radius_sourcePast = nnData.norms[0];
- }*/
- if (nnData.norms[0] > radius_destPast) {
- radius_destPast = nnData.norms[0];
- }
- if (nnData.norms[1] > radius_destNext) {
- radius_destNext = nnData.norms[1];
- radius_destNext_sampleIndex = nnData.sampleIndex;
- }
- if(nnData.distance > radius_max) {
- radius_max = nnData.distance;
- }
- }
+ // // Find eps_{x,y,z} as the maximum x, y and z norms amongst this set:
+ // radius_sourcePast = 0.0;
+ // radius_destPast = 0.0;
+ // radius_destNext = 0.0;
+ // radius_destNext_sampleIndex = -1;
+ // radius_max = 0.0;
+ // for (int j = 0; j < Knns; j++) {
+ // // Take the furthest remaining of the nearest neighbours from the PQ:
+ // NeighbourNodeData nnData = nnPQ.poll();
+ // /*if (nnData.norms[0] > radius_sourcePast) {
+ // radius_sourcePast = nnData.norms[0];
+ // }*/
+ // if (nnData.norms[0] > radius_destPast) {
+ // radius_destPast = nnData.norms[0];
+ // }
+ // if (nnData.norms[1] > radius_destNext) {
+ // radius_destNext = nnData.norms[1];
+ // radius_destNext_sampleIndex = nnData.sampleIndex;
+ // }
+ // if(nnData.distance > radius_max) {
+ // radius_max = nnData.distance;
+ // }
+ // }
//radius_max = Math.max(radius_destPast, radius_destNext);
@@ -834,8 +834,8 @@ public class TransferEntropyCalculatorSpikingIntegrationWindows2 implements
double rawRateGivenSourceAndDest = (double) Knns / timeInWindowWithMatchingJointHistories;
double rawRateGivenDest = (double) countOfDestNextMatched / timeInWindowWithMatchingDestHistory;
// Attempt at bias correction:
- double logRateGivenSourceAndDestCorrected = -Math.log(timeInWindowWithMatchingJointHistories);
- double logRateGivenDestCorrected = -Math.log(timeInWindowWithMatchingDestHistory);
+ double logRateGivenSourceAndDestCorrected = MathsUtils.digamma(Knns) - Math.log(timeInWindowWithMatchingJointHistories);
+ double logRateGivenDestCorrected = MathsUtils.digamma(countOfDestNextMatched) - Math.log(timeInWindowWithMatchingDestHistory);
//System.out.println(logRateGivenSourceAndDestCorrected + " " + logRateGivenDestCorrected);
//============================
diff --git a/java/source/infodynamics/utils/KdTree.java b/java/source/infodynamics/utils/KdTree.java
old mode 100755
new mode 100644
index 1076af7..a8c8f09
--- a/java/source/infodynamics/utils/KdTree.java
+++ b/java/source/infodynamics/utils/KdTree.java
@@ -275,10 +275,10 @@ public class KdTree extends NearestNeighbourSearcher {
// Could remove these since the code is now functional,
// but may be better to leave them in just in case the code breaks:
if (leftIndex > leftStart + leftNumPoints) {
- throw new RuntimeException("Exceeded expected number of points on left");
+ throw new RuntimeException("Exceeded expected number of points on left - likely had an NaN in the data");
}
if (rightIndex > rightStart + rightNumPoints) {
- throw new RuntimeException("Exceeded expected number of points on right");
+ throw new RuntimeException("Exceeded expected number of points on right - likely had an NaN in the data");
}
// Update the pointer for the sorted indices for this dimension,
// and keep the new temporary array
@@ -1384,19 +1384,19 @@ public class KdTree extends NearestNeighbourSearcher {
}
@Override
- public void findPointsWithinR(int sampleIndex,
+ public int findPointsWithinR(int sampleIndex,
double r, boolean allowEqualToR, boolean[] isWithinR,
int[] indicesWithinR) {
if (rootNode == null) {
indicesWithinR[0] = -1;
- return;
+ return 0;
}
int currentIndexInIndicesWithinR =
findPointsWithinR(sampleIndex,
rootNode, 0, r, allowEqualToR, isWithinR,
indicesWithinR, 0);
indicesWithinR[currentIndexInIndicesWithinR] = -1;
- return;
+ return currentIndexInIndicesWithinR;
}
/**
@@ -2510,20 +2510,21 @@ public class KdTree extends NearestNeighbourSearcher {
* @param allowEqualToR if true, then count points at radius r also, otherwise only those strictly within r
* @param isWithinR the array MUST be passed in with all points set to false initially, and is returned indicating whether each sample was found to be within r of that at sampleIndex.
* @param indicesWithinR a list of array indices for points marked as true in isWithinR, terminated with a -1 value.
+ * @return the number of points found to match
*/
- public void findPointsWithinRs(int sampleIndex,
+ public int findPointsWithinRs(int sampleIndex,
double[] rs, int dynCorrExclTime, boolean allowEqualToR,
boolean[] isWithinR, int[] indicesWithinR) {
if (rootNode == null) {
indicesWithinR[0] = -1;
- return;
+ return 0;
}
int currentIndexInIndicesWithinR =
findPointsWithinRs(sampleIndex,
rootNode, 0, rs, dynCorrExclTime, allowEqualToR, isWithinR,
indicesWithinR, 0);
indicesWithinR[currentIndexInIndicesWithinR] = -1;
- return;
+ return currentIndexInIndicesWithinR;
}
/**
diff --git a/java/source/infodynamics/utils/NearestNeighbourSearcher.java b/java/source/infodynamics/utils/NearestNeighbourSearcher.java
old mode 100755
new mode 100644
index 5f66dcb..3f2c90a
--- a/java/source/infodynamics/utils/NearestNeighbourSearcher.java
+++ b/java/source/infodynamics/utils/NearestNeighbourSearcher.java
@@ -326,8 +326,9 @@ public abstract class NearestNeighbourSearcher {
* found to be within r of that at sampleIndex.
* @param indicesWithinR a list of array indices
* for points marked as true in isWithinR, terminated with a -1 value.
+ * @return the number of matching points found
*/
- public abstract void findPointsWithinR(
+ public abstract int findPointsWithinR(
int sampleIndex, double r,
boolean allowEqualToR, boolean[] isWithinR, int[] indicesWithinR);
diff --git a/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java b/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java
index 38b0731..a06c456 100644
--- a/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java
+++ b/java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java
@@ -685,7 +685,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
*
*
*/
- public void findPointsWithinR(
+ public int findPointsWithinR(
int sampleIndex, double r, boolean allowEqualToR,
boolean[] isWithinR, int[] indicesWithinR) {
int indexInIndicesWithinR = 0;
@@ -718,7 +718,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
break;
}
// Write the terminating integer into the indicesWithinR array:
- indicesWithinR[indexInIndicesWithinR++] = -1;
+ indicesWithinR[indexInIndicesWithinR] = -1;
+ return indexInIndicesWithinR;
}
/**
diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/GPUPerformanceTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/GPUPerformanceTester.java
index 6be9705..8436052 100644
--- a/java/unittests/infodynamics/measures/continuous/kraskov/GPUPerformanceTester.java
+++ b/java/unittests/infodynamics/measures/continuous/kraskov/GPUPerformanceTester.java
@@ -86,14 +86,14 @@ public class GPUPerformanceTester extends TestCase {
startTime = Calendar.getInstance().getTimeInMillis();
miCalc.initialise(1,1);
miCalc.setObservations(source, dest);
- cpu_val = miCalc.computeSignificance(nb_surrogates).getMeanOfDistribution();
+ cpu_val = miCalc.computeSignificance(nb_surrogates).actualValue;
cpu_duration = Calendar.getInstance().getTimeInMillis() - startTime;
miCalc.setProperty("USE_GPU", "true");
miCalc.initialise(1,1);
startTime = Calendar.getInstance().getTimeInMillis();
miCalc.setObservations(source, dest);
- gpu_val = miCalc.computeSignificance(nb_surrogates).getMeanOfDistribution();
+ gpu_val = miCalc.computeSignificance(nb_surrogates).actualValue;
gpu_duration = Calendar.getInstance().getTimeInMillis() - startTime;
assertEquals(cpu_val, gpu_val, 0.0001);
@@ -124,12 +124,12 @@ public class GPUPerformanceTester extends TestCase {
return;
}
- int timeSteps = 10000;
+ int timeSteps = 1000;
RandomGenerator rg = new RandomGenerator();
double[] source = rg.generateNormalData(timeSteps, 0, 1);
double[] dest = rg.generateNormalData(timeSteps, 0, 1);
- compareGPUPerformance(source, dest, "Random low-dimensional data");
+ compareGPUPerformanceSurrogates(source, dest, 150, "Random low-dimensional data");
return;
}
diff --git a/readme-template.txt b/readme-template.txt
index 44351a6..1a11192 100755
--- a/readme-template.txt
+++ b/readme-template.txt
@@ -134,6 +134,18 @@ Notices for this software are found in the notices/JAMA directory.
Release notes
===============
+v1.5 26/11/2018
+---------------
+(after 753 commits recorded by github, repository as at https://github.com/jlizier/jidt/tree/603445651cc0bf155a42c9ba336141bc7f29bccd save for this file update)
+Added GPU (cuda) capability for KSG Conditional Mutual Information calculator (proper documentation to come), including unit tests and brief wiki page;
+Added auto-embedding for TE/AIS with multivariate KSG, and univariate and multivariate Gaussian estimator (plus unit tests), for Ragwitz criteria and Maximum bias-corrected AIS, and also added Maximum bias corrected AIS and TE to handle source embedding as well;
+Kozachenko entropy estimator adds noise to data by default;
+Added bias-correction property to Gaussian and Kernel estimators for MI and conditional MI, including with surrogates (only option for kernel);
+Enabled use of different bases for different variables in MI discrete estimator;
+All new above features enabled in AutoAnalyser;
+Added drop-down menus for parameters in AutoAnalyser;
+Included long-form lecture slides in course folder;
+
v1.4 26/11/2017
---------------
(after 638 commits recorded by github, repository as at https://github.com/jlizier/jidt/tree/589d51674e6a9cfb569432679e515bea17092876 save for this file update)
diff --git a/tester.py b/tester.py
new file mode 100755
index 0000000..bd8265b
--- /dev/null
+++ b/tester.py
@@ -0,0 +1,186 @@
+##
+## 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 .
+##
+
+# = Example 4 - Transfer entropy on continuous data using Kraskov estimators =
+
+# Simple transfer entropy (TE) calculation on continuous-valued data using the Kraskov-estimator TE calculator.
+
+from jpype import *
+import random
+import math
+import os
+import numpy as np
+
+
+NUM_REPS = 2
+NUM_SPIKES = int(2e3)
+NUM_OBSERVATIONS = 2
+NUM_SURROGATES = 10
+
+# Params for canonical example generation
+RATE_Y = 1.0
+RATE_X_MAX = 10
+
+def generate_canonical_example_processes(num_y_events):
+ event_train_x = []
+ event_train_x.append(0)
+
+ event_train_y = np.random.uniform(0, int(num_y_events / RATE_Y), int(num_y_events))
+ event_train_y.sort()
+
+ most_recent_y_index = 0
+ previous_x_candidate = 0
+ while most_recent_y_index < (len(event_train_y) - 1):
+
+ this_x_candidate = previous_x_candidate + random.expovariate(RATE_X_MAX)
+
+ while most_recent_y_index < (len(event_train_y) - 1) and this_x_candidate > event_train_y[most_recent_y_index + 1]:
+ most_recent_y_index += 1
+
+ delta_t = this_x_candidate - event_train_y[most_recent_y_index]
+
+ rate = 0
+
+ if delta_t > 1:
+ rate = 0.5
+ else:
+ rate = 0.5 + 5.0 * math.exp(-50 * (delta_t - 0.5)**2) - 5.0 * math.exp(-50 * (0.5)**2)
+ if random.random() < rate/float(RATE_X_MAX):
+ event_train_x.append(this_x_candidate)
+ previous_x_candidate = this_x_candidate
+
+ event_train_x.sort()
+ event_train_y.sort()
+
+ return event_train_x, event_train_y
+
+# Change location of jar to match yours (we assume script is called from demos/python):
+jarLocation = os.path.join(os.getcwd(), "infodynamics.jar");
+if (not(os.path.isfile(jarLocation))):
+ exit("infodynamics.jar not found (expected at " + os.path.abspath(jarLocation) + ") - are you running from demos/python?")
+# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
+startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
+teCalcClass = JPackage("infodynamics.measures.spiking.integration").TransferEntropyCalculatorSpikingIntegration
+
+
+teCalc = teCalcClass()
+teCalc.setProperty("knns", "4")
+print("Independent Poisson Processes")
+teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
+teCalc.setProperty("SOURCE_PAST_INTERVALS", "1,2")
+teCalc.appendConditionalIntervals(JArray(JInt, 1)([1, 2]))
+teCalc.appendConditionalIntervals(JArray(JInt, 1)([1, 2]))
+teCalc.setProperty("NORM_TYPE", "MAX_NORM")
+
+results_poisson = np.zeros(NUM_REPS)
+for i in range(NUM_REPS):
+ teCalc.startAddObservations()
+ for j in range(NUM_OBSERVATIONS):
+ sourceArray = NUM_SPIKES*np.random.random(NUM_SPIKES)
+ sourceArray.sort()
+ destArray = NUM_SPIKES*np.random.random(NUM_SPIKES)
+ destArray.sort()
+ condArray = NUM_SPIKES*np.random.random((2, NUM_SPIKES))
+ condArray.sort(axis = 1)
+ teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray), JArray(JDouble, 2)(condArray))
+ teCalc.finaliseAddObservations();
+ result = teCalc.computeAverageLocalOfObservations()
+ print("TE result %.4f nats" % (result,))
+ sig = teCalc.computeSignificance(NUM_SURROGATES, result)
+ print(sig.pValue)
+ results_poisson[i] = result
+print("Summary: mean ", np.mean(results_poisson), " std dev ", np.std(results_poisson))
+
+teCalc = teCalcClass()
+teCalc.setProperty("knns", "4")
+print("Noisy copy zero TE")
+#teCalc.appendConditionalIntervals(JArray(JInt, 1)([1]))
+teCalc.setProperty("DEST_PAST_INTERVALS", "1")
+teCalc.setProperty("SOURCE_PAST_INTERVALS", "1")
+#teCalc.setProperty("NORM_TYPE", "MAX_NORM")
+
+results_noisy_zero = np.zeros(NUM_REPS)
+for i in range(NUM_REPS):
+ teCalc.startAddObservations()
+ for j in range(NUM_OBSERVATIONS):
+ condArray = np.ones((1, NUM_SPIKES)) + 0.05 * np.random.random((1, NUM_SPIKES))
+ condArray = np.cumsum(condArray, axis = 1)
+ condArray.sort(axis = 1)
+ sourceArray = condArray[0, :] + 0.25 + 0.05 * np.random.normal(size = condArray.shape[1])
+ sourceArray.sort()
+ destArray = condArray[0, :] + 0.5 + 0.05 * np.random.normal(size = condArray.shape[1])
+ destArray.sort()
+ #teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray), JArray(JDouble, 2)(condArray))
+ teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray))
+ teCalc.finaliseAddObservations();
+ result = teCalc.computeAverageLocalOfObservations()
+ print("TE result %.4f nats" % (result,))
+ sig = teCalc.computeSignificance(NUM_SURROGATES, result)
+ print(sig.pValue)
+ results_noisy_zero[i] = result
+print("Summary: mean ", np.mean(results_noisy_zero), " std dev ", np.std(results_noisy_zero))
+
+
+
+teCalc = teCalcClass()
+teCalc.setProperty("knns", "4")
+print("Noisy copy non-zero TE")
+teCalc.appendConditionalIntervals(JArray(JInt, 1)([1]))
+teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
+teCalc.setProperty("SOURCE_PAST_INTERVALS", "1")
+#teCalc.setProperty("NORM_TYPE", "MAX_NORM")
+
+results_noisy_non_zero = np.zeros(NUM_REPS)
+for i in range(NUM_REPS):
+ teCalc.startAddObservations()
+ for j in range(NUM_OBSERVATIONS):
+ sourceArray = np.ones((1, NUM_SPIKES)) + 0.05 * np.random.random((1, NUM_SPIKES))
+ sourceArray = np.cumsum(sourceArray)
+ sourceArray.sort()
+ condArray = sourceArray + 0.25 + 0.05 * np.random.normal(size = sourceArray.shape)
+ condArray.sort()
+ condArray = np.expand_dims(condArray, 0)
+ destArray = sourceArray + 0.5 + 0.05 * np.random.normal(size = sourceArray.shape)
+ destArray.sort()
+ teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray), JArray(JDouble, 2)(condArray))
+ teCalc.finaliseAddObservations();
+ result = teCalc.computeAverageLocalOfObservations()
+ print("TE result %.4f nats" % (result,))
+ sig = teCalc.computeSignificance(NUM_SURROGATES, result)
+ print(sig.pValue)
+ results_noisy_non_zero[i] = result
+print("Summary: mean ", np.mean(results_noisy_non_zero), " std dev ", np.std(results_noisy_zero))
+
+print("Canonical example")
+teCalc = teCalcClass()
+teCalc.setProperty("knns", "4")
+teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
+teCalc.setProperty("SOURCE_PAST_INTERVALS", "1")
+#teCalc.setProperty("NUM_SAMPLES_MULTIPLIER", "1")
+#teCalc.setProperty("NORM_TYPE", "MAX_NORM")
+
+results_canonical = np.zeros(NUM_REPS)
+for i in range(NUM_REPS):
+ event_train_x, event_train_y = generate_canonical_example_processes(NUM_SPIKES)
+ teCalc.setObservations(JArray(JDouble, 1)(event_train_y), JArray(JDouble, 1)(event_train_x))
+ result = teCalc.computeAverageLocalOfObservations()
+ results_canonical[i] = result
+ print("TE result %.4f nats" % (result,))
+ sig = teCalc.computeSignificance(NUM_SURROGATES, result)
+ print(sig.pValue)
+print("Summary: mean ", np.mean(results_canonical), " std dev ", np.std(results_canonical))
diff --git a/tutorial/README-TutorialAndExercise.pdf b/tutorial/README-TutorialAndExercise.pdf
index 7e487b8..596baa3 100644
Binary files a/tutorial/README-TutorialAndExercise.pdf and b/tutorial/README-TutorialAndExercise.pdf differ