From dc85001ff717396f294cf4dd77f9cd4e1aa4b82c Mon Sep 17 00:00:00 2001 From: David Shorten Date: Mon, 2 Aug 2021 22:06:07 +1000 Subject: [PATCH] refactoring + conditional processes + radius sharing + euclidean norm working --- ...erEntropyCalculatorSpikingIntegration.java | 663 +++++++++++++----- java/source/infodynamics/utils/KdTree.java | 16 +- tester.py | 98 ++- 3 files changed, 564 insertions(+), 213 deletions(-) diff --git a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java index 67ccdf9..103cd27 100644 --- a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java +++ b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java @@ -7,6 +7,7 @@ import java.util.PriorityQueue; import java.util.Random; import java.util.Vector; +//import infodynamics.measures.continuous.kraskov.EuclideanUtils; import infodynamics.measures.spiking.TransferEntropyCalculatorSpiking; import infodynamics.utils.EmpiricalMeasurementDistribution; import infodynamics.utils.KdTree; @@ -15,61 +16,85 @@ import infodynamics.utils.MatrixUtils; 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 { +public class TransferEntropyCalculatorSpikingIntegration implements TransferEntropyCalculatorSpiking { /** - * Number of past destination spikes to consider (akin to embedding length) + * Number of past destination interspike intervals to consider (akin to embedding length) */ protected int k = 1; /** - * Number of past source spikes to consider (akin to embedding length) + * Number of past source interspike intervals to consider (akin to embedding length) */ protected int l = 1; + /** + * Property name for number of interspike intervals for the conditional variables + */ + public static final String COND_EMBED_LENGTHS_PROP_NAME = "COND_EMBED_LENGTHS"; + /** + * Array of history interspike interval embedding lengths for the conditional variables. + * Can be an empty array or null if there are no conditional variables. + */ + protected int[] condEmbedDims = new int[] {}; + /** * 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; - - Vector targetEmbeddingsFromSpikes = null; + + /** + * Storage for conditional observations supplied via + * {@link #addObservations(double[], double[])} etc. + */ + protected Vector vectorOfConditionalSpikeTimes = null; + + Vector ConditioningEmbeddingsFromSpikes = null; Vector jointEmbeddingsFromSpikes = null; - Vector targetEmbeddingsFromSamples = null; + Vector ConditioningEmbeddingsFromSamples = null; Vector jointEmbeddingsFromSamples = null; + Vector processTimeLengths = null; protected KdTree kdTreeJointAtSpikes = null; protected KdTree kdTreeJointAtSamples = null; protected KdTree kdTreeConditioningAtSpikes = null; protected KdTree kdTreeConditioningAtSamples = null; - + public static final String KNNS_PROP_NAME = "Knns"; /** - * Property name for an amount of random Gaussian noise to be - * added to the data (default is 1e-8, matching the MILCA toolkit). + * 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 PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD"; + public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD"; + /** * Whether to add an amount of random noise to the incoming data */ @@ -79,35 +104,61 @@ public class TransferEntropyCalculatorSpikingIntegration implements */ protected double noiseLevel = (double) 1e-8; - protected boolean trimToPosNextSpikeTimes = false; - /** * 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 num_samples_multiplier = 1.0; + /** + * 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(k, l); } - - /* (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(k, this.l); } - /* (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 { @@ -120,40 +171,74 @@ public class TransferEntropyCalculatorSpikingIntegration implements 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); + int k_temp = Integer.parseInt(propertyValue); + if (k_temp < 1) { + throw new Exception ("Invalid k value less than 1."); + } else { + k = k_temp; + } } else if (propertyName.equalsIgnoreCase(L_PROP_NAME)) { - l = Integer.parseInt(propertyValue); + int l_temp = Integer.parseInt(propertyValue); + if (l_temp < 1) { + throw new Exception ("Invalid l value less than 1."); + } else { + l = l_temp; + } + } else if (propertyName.equalsIgnoreCase(COND_EMBED_LENGTHS_PROP_NAME)) { + int[] condEmbedDims_temp = ParsedProperties.parseStringArrayOfInts(propertyValue); + for (int dim : condEmbedDims_temp) { + if (dim < 1) { + throw new Exception ("Invalid conditional embedding value less than 1."); + } + } + condEmbedDims = condEmbedDims_temp; } else if (propertyName.equalsIgnoreCase(KNNS_PROP_NAME)) { Knns = 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_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 temp_num_samples_multiplier = Double.parseDouble(propertyValue); + if (temp_num_samples_multiplier <= 0) { + throw new Exception ("Num samples multiplier must be greater than 0."); + } else { + num_samples_multiplier = temp_num_samples_multiplier; + } + } else if (propertyName.equalsIgnoreCase(PROP_NORM_TYPE)) { + normType = KdTree.validateNormType(propertyValue); } 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 { @@ -163,121 +248,152 @@ public class TransferEntropyCalculatorSpikingIntegration implements return Integer.toString(l); } else if (propertyName.equalsIgnoreCase(KNNS_PROP_NAME)) { return Integer.toString(Knns); - } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) { - return Double.toString(noiseLevel); + } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) { + return Double.toString(noiseLevel); + } else if (propertyName.equalsIgnoreCase(PROP_SAMPLE_MULTIPLIER)) { + return Double.toString(num_samples_multiplier); } else { // No property matches for this class return null; } } - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#setObservations(double[], double[]) + /* + * (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 { - targetEmbeddingsFromSpikes = new Vector(); + ConditioningEmbeddingsFromSpikes = new Vector(); jointEmbeddingsFromSpikes = new Vector(); - targetEmbeddingsFromSamples = new Vector(); + ConditioningEmbeddingsFromSamples = new Vector(); jointEmbeddingsFromSamples = new Vector(); - + processTimeLengths = new Vector(); + // Send all of the observations through: Iterator sourceIterator = vectorOfSourceSpikeTimes.iterator(); int timeSeriesIndex = 0; - for (double[] destSpikeTimes : vectorOfDestinationSpikeTimes) { - double[] sourceSpikeTimes = sourceIterator.next(); - processEventsFromSpikingTimeSeries(sourceSpikeTimes, destSpikeTimes, - targetEmbeddingsFromSpikes, jointEmbeddingsFromSpikes, - targetEmbeddingsFromSamples, jointEmbeddingsFromSamples); + if (vectorOfConditionalSpikeTimes.size() > 0) { + Iterator conditionalIterator = vectorOfConditionalSpikeTimes.iterator(); + for (double[] destSpikeTimes : vectorOfDestinationSpikeTimes) { + double[] sourceSpikeTimes = sourceIterator.next(); + double[][] conditionalSpikeTimes = conditionalIterator.next(); + processEventsFromSpikingTimeSeries(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes, ConditioningEmbeddingsFromSpikes, + jointEmbeddingsFromSpikes, ConditioningEmbeddingsFromSamples, jointEmbeddingsFromSamples, + processTimeLengths); + } + } else { + for (double[] destSpikeTimes : vectorOfDestinationSpikeTimes) { + double[] sourceSpikeTimes = sourceIterator.next(); + double[][] conditionalSpikeTimes = new double[][] {}; + processEventsFromSpikingTimeSeries(sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes, ConditioningEmbeddingsFromSpikes, + jointEmbeddingsFromSpikes, ConditioningEmbeddingsFromSamples, jointEmbeddingsFromSamples, + processTimeLengths); + } } // Convert the vectors to arrays so that they can be put in the trees - double[][] arrayedTargetEmbeddingsFromSpikes = new double[targetEmbeddingsFromSpikes.size()][k]; - double[][] arrayedJointEmbeddingsFromSpikes = new double[targetEmbeddingsFromSpikes.size()][k + l]; - for (int i = 0; i < targetEmbeddingsFromSpikes.size(); i++) { - arrayedTargetEmbeddingsFromSpikes[i] = targetEmbeddingsFromSpikes.elementAt(i); + double[][] arrayedTargetEmbeddingsFromSpikes = new double[ConditioningEmbeddingsFromSpikes.size()][k]; + double[][] arrayedJointEmbeddingsFromSpikes = new double[ConditioningEmbeddingsFromSpikes.size()][k + l]; + for (int i = 0; i < ConditioningEmbeddingsFromSpikes.size(); i++) { + arrayedTargetEmbeddingsFromSpikes[i] = ConditioningEmbeddingsFromSpikes.elementAt(i); arrayedJointEmbeddingsFromSpikes[i] = jointEmbeddingsFromSpikes.elementAt(i); } - double[][] arrayedTargetEmbeddingsFromSamples = new double[targetEmbeddingsFromSamples.size()][k]; - double[][] arrayedJointEmbeddingsFromSamples = new double[targetEmbeddingsFromSamples.size()][k + l]; - for (int i = 0; i < targetEmbeddingsFromSamples.size(); i++) { - arrayedTargetEmbeddingsFromSamples[i] = targetEmbeddingsFromSamples.elementAt(i); + double[][] arrayedTargetEmbeddingsFromSamples = new double[ConditioningEmbeddingsFromSamples.size()][k]; + double[][] arrayedJointEmbeddingsFromSamples = new double[ConditioningEmbeddingsFromSamples.size()][k + l]; + for (int i = 0; i < ConditioningEmbeddingsFromSamples.size(); i++) { + arrayedTargetEmbeddingsFromSamples[i] = ConditioningEmbeddingsFromSamples.elementAt(i); arrayedJointEmbeddingsFromSamples[i] = jointEmbeddingsFromSamples.elementAt(i); } - kdTreeJointAtSpikes = new KdTree( - new int[] {k + l}, - new double[][][] {arrayedJointEmbeddingsFromSpikes}); - kdTreeJointAtSamples = new KdTree( - new int[] {k + l}, - new double[][][] {arrayedJointEmbeddingsFromSamples}); - kdTreeConditioningAtSpikes = new KdTree( - new int[] {k}, - new double[][][] {arrayedTargetEmbeddingsFromSpikes}); - kdTreeConditioningAtSamples = new KdTree( - new int[] {k}, - new double[][][] {arrayedTargetEmbeddingsFromSamples}); + kdTreeJointAtSpikes = new KdTree(arrayedJointEmbeddingsFromSpikes); + kdTreeJointAtSamples = new KdTree(arrayedJointEmbeddingsFromSamples); + kdTreeConditioningAtSpikes = new KdTree(arrayedTargetEmbeddingsFromSpikes); + kdTreeConditioningAtSamples = new KdTree(arrayedTargetEmbeddingsFromSamples); - /*kdTreeJointAtSpikes.setNormType("EUCLIDEAN"); - kdTreeJointAtSamples.setNormType("EUCLIDEAN"); - kdTreeConditioningAtSpikes.setNormType("EUCLIDEAN"); - kdTreeConditioningAtSamples.setNormType("EUCLIDEAN");*/ + kdTreeJointAtSpikes.setNormType(normType); + kdTreeJointAtSamples.setNormType(normType); + kdTreeConditioningAtSpikes.setNormType(normType); + kdTreeConditioningAtSamples.setNormType(normType); } - protected void makeEmbeddingsAtPoints(double[] pointsAtWhichToMakeEmbeddings, double[] sourceSpikeTimes, double[] destSpikeTimes, - Vector targetEmbeddings, Vector jointEmbeddings) { - //System.out.println("foo"); + protected void makeEmbeddingsAtPoints(double[] pointsAtWhichToMakeEmbeddings, int index_of_first_point_to_use, + double[] sourceSpikeTimes, double[] destSpikeTimes, + double[][] conditionalSpikeTimes, + Vector ConditioningEmbeddings, + Vector jointEmbeddings) { Random random = new Random(); - int embedding_point_index = 0; + int embedding_point_index = index_of_first_point_to_use; int most_recent_dest_index = k; int most_recent_source_index = l; - - // Make sure that the first point at which an embedding is made has enough preceding spikes in both source and - // target for embeddings to be made. - while (pointsAtWhichToMakeEmbeddings[embedding_point_index] <= destSpikeTimes[most_recent_dest_index] | - pointsAtWhichToMakeEmbeddings[embedding_point_index] <= sourceSpikeTimes[most_recent_source_index]) { - embedding_point_index++; + int[] most_recent_conditioning_indices = Arrays.copyOf(condEmbedDims, condEmbedDims.length); + int total_length_of_conditioning_embeddings = 0; + for (int i = 0; i < condEmbedDims.length; i++) { + total_length_of_conditioning_embeddings += condEmbedDims[i]; } // Loop through the points at which embeddings need to be made - for (;embedding_point_index < pointsAtWhichToMakeEmbeddings.length; embedding_point_index++) { + for (; embedding_point_index < pointsAtWhichToMakeEmbeddings.length; embedding_point_index++) { // Advance the tracker of the most recent dest index - while (most_recent_dest_index < (destSpikeTimes.length - 1)) { + while (most_recent_dest_index < (destSpikeTimes.length - 1)) { if (destSpikeTimes[most_recent_dest_index + 1] < pointsAtWhichToMakeEmbeddings[embedding_point_index]) { most_recent_dest_index++; } else { @@ -285,163 +401,328 @@ public class TransferEntropyCalculatorSpikingIntegration implements } } // Do the same for the most recent source index - while (most_recent_source_index < (sourceSpikeTimes.length - 1)) { - if (sourceSpikeTimes[most_recent_source_index + 1] < pointsAtWhichToMakeEmbeddings[embedding_point_index]) { + while (most_recent_source_index < (sourceSpikeTimes.length - 1)) { + if (sourceSpikeTimes[most_recent_source_index + + 1] < pointsAtWhichToMakeEmbeddings[embedding_point_index]) { most_recent_source_index++; } else { break; } } + // Now advance the trackers for the most recent conditioning indices + for (int j = 0; j < most_recent_conditioning_indices.length; j++) { + while (most_recent_conditioning_indices[j] < (conditionalSpikeTimes[j].length - 1)) { + if (conditionalSpikeTimes[j][most_recent_conditioning_indices[j] + 1] < pointsAtWhichToMakeEmbeddings[embedding_point_index]) { + most_recent_conditioning_indices[j]++; + } else { + break; + } + } + } - double[] destPast = new double[k]; - double[] jointPast = new double[k + l]; - destPast[0] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - - destSpikeTimes[most_recent_dest_index]; - jointPast[0] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - - destSpikeTimes[most_recent_dest_index]; - jointPast[k] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - - sourceSpikeTimes[most_recent_source_index]; + double[] conditioningPast = new double[k + total_length_of_conditioning_embeddings]; + double[] jointPast = new double[k + total_length_of_conditioning_embeddings + l]; + + // Add the embedding intervals from the target process + conditioningPast[0] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - destSpikeTimes[most_recent_dest_index]; + jointPast[0] = pointsAtWhichToMakeEmbeddings[embedding_point_index] + - destSpikeTimes[most_recent_dest_index]; for (int i = 1; i < k; i++) { - destPast[i] = destSpikeTimes[most_recent_dest_index - i + 1] - - destSpikeTimes[most_recent_dest_index - i]; - jointPast[i] = destSpikeTimes[most_recent_dest_index - i + 1] - - destSpikeTimes[most_recent_dest_index - i]; - } - for (int i = 1; i < l; i++) { - jointPast[k + i] = sourceSpikeTimes[most_recent_source_index - i + 1] - - sourceSpikeTimes[most_recent_source_index - i]; + conditioningPast[i] = destSpikeTimes[most_recent_dest_index - i + 1] + - destSpikeTimes[most_recent_dest_index - i]; + jointPast[i] = destSpikeTimes[most_recent_dest_index - i + 1] + - destSpikeTimes[most_recent_dest_index - i]; } - if (addNoise) { - for (int i = 0; i < k; i++) { - destPast[i] += random.nextGaussian()*noiseLevel; - } - for (int i = 0; i < l; i++) { - jointPast[i] += random.nextGaussian()*noiseLevel; + // Add the embeding intervals from the conditional processes + int index_of_next_embedding_interval = k; + for (int i = 0; i < condEmbedDims.length; i++) { + conditioningPast[index_of_next_embedding_interval] = + pointsAtWhichToMakeEmbeddings[embedding_point_index] - conditionalSpikeTimes[i][most_recent_conditioning_indices[i]]; + jointPast[index_of_next_embedding_interval] = + pointsAtWhichToMakeEmbeddings[embedding_point_index] - conditionalSpikeTimes[i][most_recent_conditioning_indices[i]]; + index_of_next_embedding_interval += 1; + for (int j = 1; j < condEmbedDims[i]; j++) { + conditioningPast[index_of_next_embedding_interval] = + conditionalSpikeTimes[i][most_recent_conditioning_indices[i] - j + 1] - + conditionalSpikeTimes[i][most_recent_conditioning_indices[i] - j]; + jointPast[index_of_next_embedding_interval] = + conditionalSpikeTimes[i][most_recent_conditioning_indices[i] - j + 1] - + conditionalSpikeTimes[i][most_recent_conditioning_indices[i] - j]; + index_of_next_embedding_interval += 1; } } - - targetEmbeddings.add(destPast); + + // Add the embedding intervals from the source process (this only gets added to the joint embeddings) + jointPast[k + total_length_of_conditioning_embeddings] = pointsAtWhichToMakeEmbeddings[embedding_point_index] + - sourceSpikeTimes[most_recent_source_index]; + for (int i = 1; i < l; i++) { + jointPast[k + total_length_of_conditioning_embeddings + i] = sourceSpikeTimes[most_recent_source_index - i + 1] + - sourceSpikeTimes[most_recent_source_index - i]; + } + + // Add Gaussian noise, if necessary + if (addNoise) { + for (int i = 0; i < conditioningPast.length; i++) { + conditioningPast[i] += random.nextGaussian() * noiseLevel; + } + for (int i = 0; i < jointPast.length; i++) { + jointPast[i] += random.nextGaussian() * noiseLevel; + } + } + + ConditioningEmbeddings.add(conditioningPast); jointEmbeddings.add(jointPast); } } - protected void processEventsFromSpikingTimeSeries(double[] sourceSpikeTimes, double[] destSpikeTimes, - Vector targetEmbeddingsFromSpikes, Vector jointEmbeddingsFromSpikes, - Vector targetEmbeddingsFromSamples, Vector jointEmbeddingsFromSamples) - throws Exception { - // addObservationsAfterParamsDetermined(sourceSpikeTimes, destSpikeTimes); - + protected void processEventsFromSpikingTimeSeries(double[] sourceSpikeTimes, double[] destSpikeTimes, double[][] conditionalSpikeTimes, + Vector ConditioningEmbeddingsFromSpikes, Vector jointEmbeddingsFromSpikes, + Vector ConditioningEmbeddingsFromSamples, Vector jointEmbeddingsFromSamples, + Vector processTimeLengths) + throws Exception { + // First sort the spike times in case they were not properly in ascending order: Arrays.sort(sourceSpikeTimes); Arrays.sort(destSpikeTimes); + + int first_target_index_of_embedding = k; + while (destSpikeTimes[first_target_index_of_embedding] <= sourceSpikeTimes[l - 1]) { + first_target_index_of_embedding++; + } + if (conditionalSpikeTimes.length != condEmbedDims.length) { + throw new Exception("Number of conditional embedding lengths does not match the number of conditional processes"); + } + for (int i = 0; i < conditionalSpikeTimes.length; i++) { + while (destSpikeTimes[first_target_index_of_embedding] <= conditionalSpikeTimes[i][condEmbedDims[i]]) { + first_target_index_of_embedding++; + } + } - double sample_lower_bound = Arrays.stream(sourceSpikeTimes).min().getAsDouble(); - double sample_upper_bound = Arrays.stream(sourceSpikeTimes).max().getAsDouble(); - double[] randomSampleTimes = new double[sourceSpikeTimes.length]; + //processTimeLengths.add(destSpikeTimes[sourceSpikeTimes.length - 1] - destSpikeTimes[first_target_index_of_embedding]); + processTimeLengths.add(destSpikeTimes[destSpikeTimes.length - 1] - destSpikeTimes[first_target_index_of_embedding]); + + double sample_lower_bound = destSpikeTimes[first_target_index_of_embedding]; + double sample_upper_bound = destSpikeTimes[destSpikeTimes.length - 1]; + int num_samples = (int) Math.round(num_samples_multiplier * (destSpikeTimes.length - first_target_index_of_embedding + 1)); + double[] randomSampleTimes = new double[num_samples]; Random rand = new Random(); for (int i = 0; i < randomSampleTimes.length; i++) { randomSampleTimes[i] = sample_lower_bound + rand.nextDouble() * (sample_upper_bound - sample_lower_bound); } Arrays.sort(randomSampleTimes); - makeEmbeddingsAtPoints(destSpikeTimes, sourceSpikeTimes, destSpikeTimes, targetEmbeddingsFromSpikes, jointEmbeddingsFromSpikes); - makeEmbeddingsAtPoints(randomSampleTimes, sourceSpikeTimes, destSpikeTimes, targetEmbeddingsFromSamples, jointEmbeddingsFromSamples); + makeEmbeddingsAtPoints(destSpikeTimes, first_target_index_of_embedding, sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes, + ConditioningEmbeddingsFromSpikes, jointEmbeddingsFromSpikes); + makeEmbeddingsAtPoints(randomSampleTimes, 0, sourceSpikeTimes, destSpikeTimes, conditionalSpikeTimes, + ConditioningEmbeddingsFromSamples, jointEmbeddingsFromSamples); } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getAddedMoreThanOneObservationSet() + + /* + * (non-Javadoc) + * + * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking# + * getAddedMoreThanOneObservationSet() */ @Override public boolean getAddedMoreThanOneObservationSet() { - return (vectorOfDestinationSpikeTimes != null) && - (vectorOfDestinationSpikeTimes.size() > 1); + return (vectorOfDestinationSpikeTimes != null) && (vectorOfDestinationSpikeTimes.size() > 1); } - private double max_neighbour_distance(PriorityQueue nnPQ) { - double max_val = -1e9; - while (nnPQ.peek() != null) { - NeighbourNodeData nnData = nnPQ.poll(); - if (nnData.norms[0] > max_val) { - max_val = nnData.norms[0]; + // Class to allow returning two values in the subsequent method + private static class distanceAndNumPoints { + public double distance; + public int numPoints; + + public distanceAndNumPoints(double distance, int numPoints) { + this.distance = distance; + this.numPoints = numPoints; + } + } + + 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 max_val; + return new distanceAndNumPoints(maxDistance, i); } - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeAverageLocalOfObservations() + /* + * (non-Javadoc) + * + * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking# + * computeAverageLocalOfObservations() */ @Override public double computeAverageLocalOfObservations() throws Exception { double currentSum = 0; - for (int i = 0; i < targetEmbeddingsFromSpikes.size(); i++) { + for (int i = 0; i < ConditioningEmbeddingsFromSpikes.size(); i++) { - PriorityQueue nnPQJointSpikes = - kdTreeJointAtSpikes.findKNearestNeighbours(Knns + 1, new double[][] {jointEmbeddingsFromSpikes.elementAt(i)}); - PriorityQueue nnPQJointSamples = - kdTreeJointAtSamples.findKNearestNeighbours(Knns, new double[][] {jointEmbeddingsFromSpikes.elementAt(i)}); - PriorityQueue nnPQConditioningSpikes = - kdTreeConditioningAtSpikes.findKNearestNeighbours(Knns + 1, new double[][] {targetEmbeddingsFromSpikes.elementAt(i)}); - PriorityQueue nnPQConditioningSamples = - kdTreeConditioningAtSamples.findKNearestNeighbours(Knns, new double[][] {targetEmbeddingsFromSpikes.elementAt(i)}); + double radiusJointSpikes = kdTreeJointAtSpikes.findKNearestNeighbours(Knns, i).poll().norms[0]; + double radiusJointSamples = kdTreeJointAtSamples.findKNearestNeighbours(Knns, + new double[][] { jointEmbeddingsFromSpikes.elementAt(i) }).poll().norms[0]; - double radiusJointSpikes = max_neighbour_distance(nnPQJointSpikes); - double radiusJointSamples = max_neighbour_distance(nnPQJointSamples); - double radiusConditioningSpikes = max_neighbour_distance(nnPQConditioningSpikes); - double radiusConditioningSamples = max_neighbour_distance(nnPQConditioningSamples); + /* + 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.01; + // 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[][] { jointEmbeddingsFromSpikes.elementAt(i) }, + true, + isWithinR, + indicesWithinR); + distanceAndNumPoints temp = findMaxDistanceAndNumPointsFromIndices(jointEmbeddingsFromSpikes.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()]; + kdTreeJointAtSpikes.findPointsWithinR(radiusJointSamples + eps, + new double[][] { jointEmbeddingsFromSpikes.elementAt(i) }, + true, + isWithinR, + indicesWithinR); + distanceAndNumPoints temp = findMaxDistanceAndNumPointsFromIndices(jointEmbeddingsFromSpikes.elementAt(i), indicesWithinR, + jointEmbeddingsFromSpikes); + // -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; + } - currentSum += ((k + l) * (- Math.log(radiusJointSpikes) + Math.log(radiusJointSamples)) - + k * (Math.log(radiusConditioningSpikes) - Math.log(radiusConditioningSamples))); + + + currentSum += (MathsUtils.digamma(kJointSpikes) - MathsUtils.digamma(kJointSamples) + + ((k + l) * (-Math.log(radiusJointSpikes) + Math.log(radiusJointSamples))) - + MathsUtils.digamma(kConditioningSpikes) + MathsUtils.digamma(kConditioningSamples) + + + (k * (Math.log(radiusConditioningSpikes) - Math.log(radiusConditioningSamples)))); + if (Double.isNaN(currentSum)) { + throw new Exception(kJointSpikes + " " + kJointSamples + " " + kConditioningSpikes + " " + kConditioningSamples + "\n" + + radiusJointSpikes + " " + radiusJointSamples + " " + radiusConditioningSpikes + " " + radiusConditioningSamples); + } } // Normalise by time - currentSum /= (vectorOfDestinationSpikeTimes.elementAt(0)[vectorOfDestinationSpikeTimes.elementAt(0).length - 1] - - vectorOfDestinationSpikeTimes.elementAt(0)[0]); + double time_sum = 0; + for (Double time : processTimeLengths) { + time_sum += time; + } + currentSum /= time_sum; return currentSum; } - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeLocalOfPreviousObservations() + /* + * (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) + /* + * (non-Javadoc) + * + * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking# + * computeSignificance(int) */ @Override - public EmpiricalMeasurementDistribution computeSignificance( - int numPermutationsToCheck) throws Exception { + public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck) throws Exception { // TODO Auto-generated method stub return null; } - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeSignificance(int[][]) + /* + * (non-Javadoc) + * + * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking# + * computeSignificance(int[][]) */ @Override - public EmpiricalMeasurementDistribution computeSignificance( - int[][] newOrderings) throws Exception { + 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() { diff --git a/java/source/infodynamics/utils/KdTree.java b/java/source/infodynamics/utils/KdTree.java index 8e5b6f6..0f58755 100644 --- a/java/source/infodynamics/utils/KdTree.java +++ b/java/source/infodynamics/utils/KdTree.java @@ -337,7 +337,7 @@ public class KdTree extends NearestNeighbourSearcher { double difference = x1[d] - x2[d]; distance += difference * difference; } - return distance; + return Math.sqrt(distance); } } @@ -359,6 +359,7 @@ public class KdTree extends NearestNeighbourSearcher { */ public final static double normWithAbort(double[] x1, double[] x2, double limit, int normToUse) { + double distance = 0.0; switch (normToUse) { case EuclideanUtils.NORM_MAX_NORM: @@ -379,15 +380,17 @@ public class KdTree extends NearestNeighbourSearcher { return distance; // case EuclideanUtils.NORM_EUCLIDEAN_SQUARED: default: + // Limit is often r, so must square + limit = limit * limit; // Inlined from {@link EuclideanUtils}: for (int d = 0; d < x1.length; d++) { double difference = x1[d] - x2[d]; distance += difference * difference; if (distance > limit) { - return Double.POSITIVE_INFINITY; + return Double.POSITIVE_INFINITY; } } - return distance; + return Math.sqrt(distance); } } @@ -1305,9 +1308,10 @@ public class KdTree extends NearestNeighbourSearcher { if (normTypeToUse == EuclideanUtils.NORM_MAX_NORM) { absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; } else { + absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; // norm type is EuclideanUtils#NORM_EUCLIDEAN_SQUARED // Track the square distance - absDistOnThisDim = distOnThisDim * distOnThisDim; + //absDistOnThisDim = distOnThisDim; } if ((node.indexOfThisPoint != sampleIndex) && @@ -2053,7 +2057,7 @@ public class KdTree extends NearestNeighbourSearcher { KdTreeNode node, int level, double r, boolean allowEqualToR, boolean[] isWithinR, int[] indicesWithinR, int nextIndexInIndicesWithinR) { - + //System.out.println("foo"); // Point to the correct array for the data at this level int currentDim = level % totalDimensions; double[][] data = dimensionToArray[currentDim]; @@ -2068,9 +2072,9 @@ public class KdTree extends NearestNeighbourSearcher { if (normTypeToUse == EuclideanUtils.NORM_MAX_NORM) { absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; } else { + absDistOnThisDim = (distOnThisDim > 0) ? distOnThisDim : - distOnThisDim; // norm type is EuclideanUtils#NORM_EUCLIDEAN_SQUARED // Track the square distance - absDistOnThisDim = distOnThisDim * distOnThisDim; } if ((absDistOnThisDim < r) || diff --git a/tester.py b/tester.py index a010ee4..09817c4 100755 --- a/tester.py +++ b/tester.py @@ -27,15 +27,14 @@ import os import numpy as np -NUM_REPS = 20 -NUM_SPIKES = int(1e4) +NUM_REPS = 5 +NUM_SPIKES = int(5e3) +NUM_OBSERVATIONS = 2 # 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) @@ -76,30 +75,97 @@ if (not(os.path.isfile(jarLocation))): # 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") + +teCalc = teCalcClass() +teCalc.setProperty("knns", "4") print("Independent Poisson Processes") -teCalc.setProperty("k_HISTORY", "1") -teCalc.setProperty("l_HISTORY", "1") +teCalc.setProperty("COND_EMBED_LENGTHS", "2,2") +teCalc.setProperty("k_HISTORY", "2") +teCalc.setProperty("l_HISTORY", "2") +teCalc.setProperty("NORM_TYPE", "MAX_NORM") results_poisson = np.zeros(NUM_REPS) for i in range(NUM_REPS): - sourceArray = NUM_SPIKES*np.random.random(NUM_SPIKES) - sourceArray.sort() - destArray = NUM_SPIKES*np.random.random(NUM_SPIKES) - destArray.sort() - - teCalc.setObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray)) + 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,)) results_poisson[i] = result print("Summary: mean ", np.mean(results_poisson), " std dev ", np.std(results_poisson)) -print("Canonical example") + + + + + +teCalc = teCalcClass() +teCalc.setProperty("knns", "4") +print("Noisy copy zero TE") +teCalc.setProperty("COND_EMBED_LENGTHS", "2") teCalc.setProperty("k_HISTORY", "2") -teCalc.setProperty("l_HISTORY", "1") +teCalc.setProperty("l_HISTORY", "2") + +results_noisy_zero = np.zeros(NUM_REPS) +for i in range(NUM_REPS): + teCalc.startAddObservations() + for j in range(NUM_OBSERVATIONS): + condArray = NUM_SPIKES*np.random.random((1, NUM_SPIKES)) + condArray.sort(axis = 1) + sourceArray = condArray[0, :] + 0.25 + 0.1 * np.random.normal(size = condArray.shape[1]) + sourceArray.sort() + destArray = condArray[0, :] + 0.5 + 0.1 * np.random.normal(size = condArray.shape[1]) + 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,)) + 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 non-zero TE") +teCalc.setProperty("COND_EMBED_LENGTHS", "2") +teCalc.setProperty("k_HISTORY", "2") +teCalc.setProperty("l_HISTORY", "2") + +results_noisy_zero = 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() + condArray = sourceArray + 0.25 + 0.1 * np.random.normal(size = sourceArray.shape) + condArray.sort() + condArray = np.expand_dims(condArray, 0) + destArray = sourceArray + 0.5 + 0.1 * 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,)) + results_poisson[i] = result +print("Summary: mean ", np.mean(results_poisson), " std dev ", np.std(results_poisson)) + +print("Canonical example") +teCalc = teCalcClass() +teCalc.setProperty("knns", "4") +teCalc.setProperty("k_HISTORY", "2") +teCalc.setProperty("l_HISTORY", "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):