From 4ece80b5efde9e9985bbd84d183c76db6f81685a Mon Sep 17 00:00:00 2001 From: David Shorten Date: Sat, 24 Jul 2021 23:32:57 +1000 Subject: [PATCH] appear to have 4kl estimator working --- ...erEntropyCalculatorSpikingIntegration.java | 105 ++++++++++++++---- tester.py | 8 +- 2 files changed, 88 insertions(+), 25 deletions(-) diff --git a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java index fbe2e63..5dbc9d9 100644 --- a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java +++ b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java @@ -73,9 +73,9 @@ public class TransferEntropyCalculatorSpikingIntegration implements */ Vector targetEmbeddingsFromSpikes = null; - Vector sourceEmbeddingsFromSpikes = null; + Vector jointEmbeddingsFromSpikes = null; Vector targetEmbeddingsFromSamples = null; - Vector sourceEmbeddingsFromSamples = null; + Vector jointEmbeddingsFromSamples = null; protected KdTree kdTreeJointAtSpikes = null; protected KdTree kdTreeJointAtSamples = null; @@ -145,7 +145,7 @@ public class TransferEntropyCalculatorSpikingIntegration implements /** * Whether to add an amount of random noise to the incoming data */ - protected boolean addNoise = false; + protected boolean addNoise = true; /** * Amount of random Gaussian noise to add to the incoming data */ @@ -303,9 +303,9 @@ public class TransferEntropyCalculatorSpikingIntegration implements // New targetEmbeddingsFromSpikes = new Vector(); - sourceEmbeddingsFromSpikes = new Vector(); + jointEmbeddingsFromSpikes = new Vector(); targetEmbeddingsFromSamples = new Vector(); - sourceEmbeddingsFromSamples = new Vector(); + jointEmbeddingsFromSamples = new Vector(); // Send all of the observations through: Iterator sourceIterator = vectorOfSourceSpikeTimes.iterator(); @@ -317,18 +317,39 @@ public class TransferEntropyCalculatorSpikingIntegration implements processEventsFromSpikingTimeSeries(sourceSpikeTimes, destSpikeTimes, timeSeriesIndex, eventTimings, destPastAndNextTimings, eventTypeLocator, eventIndexLocator, numEventsPerObservationSet, - targetEmbeddingsFromSpikes, sourceEmbeddingsFromSpikes, - targetEmbeddingsFromSamples, sourceEmbeddingsFromSamples); + targetEmbeddingsFromSpikes, jointEmbeddingsFromSpikes, + targetEmbeddingsFromSamples, jointEmbeddingsFromSamples); } double[][] arrayedTargetEmbeddingsFromSpikes = new double[targetEmbeddingsFromSpikes.size()][k]; + double[][] arrayedJointEmbeddingsFromSpikes = new double[targetEmbeddingsFromSpikes.size()][k]; for (int i = 0; i < targetEmbeddingsFromSpikes.size(); i++) { arrayedTargetEmbeddingsFromSpikes[i] = targetEmbeddingsFromSpikes.elementAt(i); + arrayedJointEmbeddingsFromSpikes[i] = jointEmbeddingsFromSpikes.elementAt(i); + } + double[][] arrayedTargetEmbeddingsFromSamples = new double[targetEmbeddingsFromSamples.size()][k]; + double[][] arrayedJointEmbeddingsFromSamples = new double[targetEmbeddingsFromSamples.size()][k]; + for (int i = 0; i < targetEmbeddingsFromSamples.size(); i++) { + arrayedTargetEmbeddingsFromSamples[i] = targetEmbeddingsFromSamples.elementAt(i); + arrayedJointEmbeddingsFromSamples[i] = jointEmbeddingsFromSamples.elementAt(i); } for (int i = 0; i < 5; i++) { - System.out.println(arrayedTargetEmbeddingsFromSpikes[i][0] + " " + arrayedTargetEmbeddingsFromSpikes[i][1]); + System.out.println(arrayedTargetEmbeddingsFromSamples[i][0] + " " + arrayedTargetEmbeddingsFromSamples[i][1]); } + + 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}); // Now we have collected all the events. // Load up the search structures: @@ -383,7 +404,7 @@ public class TransferEntropyCalculatorSpikingIntegration implements } protected void makeEmbeddingsAtPoints(double[] pointsAtWhichToMakeEmbeddings, double[] sourceSpikeTimes, double[] destSpikeTimes, - Vector targetEmbeddings, Vector sourceEmbeddings) { + Vector targetEmbeddings, Vector jointEmbeddings) { //System.out.println("foo"); Random random = new Random(); @@ -425,36 +446,42 @@ public class TransferEntropyCalculatorSpikingIntegration implements //System.out.println(k + " " + l); double[] destPast = new double[k]; - double[] sourcePast = new double[l]; + double[] jointPast = new double[k + l]; destPast[0] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - destSpikeTimes[most_recent_dest_index]; - sourcePast[0] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - + jointPast[0] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - + destSpikeTimes[most_recent_dest_index]; + jointPast[k] = pointsAtWhichToMakeEmbeddings[embedding_point_index] - sourceSpikeTimes[most_recent_source_index]; if (addNoise) { destPast[0] += random.nextGaussian()*noiseLevel; - sourcePast[0] += random.nextGaussian()*noiseLevel; + jointPast[0] += random.nextGaussian()*noiseLevel; + jointPast[k] += random.nextGaussian()*noiseLevel; } 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]; if (addNoise) { - destPast[i - 1] += random.nextGaussian()*noiseLevel; + destPast[i] += random.nextGaussian()*noiseLevel; + jointPast[i] += random.nextGaussian()*noiseLevel; } } for (int i = 1; i < l; i++) { - sourcePast[i] = sourceSpikeTimes[most_recent_source_index - i + 1] - + jointPast[k + i] = sourceSpikeTimes[most_recent_source_index - i + 1] - sourceSpikeTimes[most_recent_source_index - i]; if (addNoise) { - sourcePast[i - 1] += random.nextGaussian()*noiseLevel; + jointPast[k + i] += random.nextGaussian()*noiseLevel; } } //System.out.println(Arrays.toString(destPast)); - //System.out.println(Arrays.toString(sourcePast)); + //System.out.println(Arrays.toString(jointPast)); //System.out.println(); targetEmbeddings.add(destPast); - sourceEmbeddings.add(sourcePast); + jointEmbeddings.add(jointPast); } } @@ -462,8 +489,8 @@ public class TransferEntropyCalculatorSpikingIntegration implements int timeSeriesIndex, Vector[] eventTimings, Vector destPastAndNextTimings, Vector eventTypeLocator, Vector eventIndexLocator, Vector numEventsPerObservationSet, - Vector targetEmbeddingsFromSpikes, Vector sourceEmbeddingsFromSpikes, - Vector targetEmbeddingsFromSamples, Vector sourceEmbeddingsFromSamples) + Vector targetEmbeddingsFromSpikes, Vector jointEmbeddingsFromSpikes, + Vector targetEmbeddingsFromSamples, Vector jointEmbeddingsFromSamples) throws Exception { // addObservationsAfterParamsDetermined(sourceSpikeTimes, destSpikeTimes); @@ -473,7 +500,7 @@ public class TransferEntropyCalculatorSpikingIntegration implements // New - int NUM_SAMPLES = 1000; + int NUM_SAMPLES = sourceSpikeTimes.length; double sample_lower_bound = Arrays.stream(sourceSpikeTimes).min().getAsDouble(); double sample_upper_bound = Arrays.stream(sourceSpikeTimes).max().getAsDouble(); double[] randomSampleTimes = new double[NUM_SAMPLES]; @@ -485,8 +512,8 @@ public class TransferEntropyCalculatorSpikingIntegration implements // End New - makeEmbeddingsAtPoints(destSpikeTimes, sourceSpikeTimes, destSpikeTimes, targetEmbeddingsFromSpikes, sourceEmbeddingsFromSpikes); - makeEmbeddingsAtPoints(randomSampleTimes, sourceSpikeTimes, destSpikeTimes, targetEmbeddingsFromSamples, sourceEmbeddingsFromSamples); + makeEmbeddingsAtPoints(destSpikeTimes, sourceSpikeTimes, destSpikeTimes, targetEmbeddingsFromSpikes, jointEmbeddingsFromSpikes); + makeEmbeddingsAtPoints(randomSampleTimes, sourceSpikeTimes, destSpikeTimes, targetEmbeddingsFromSamples, jointEmbeddingsFromSamples); @@ -786,11 +813,45 @@ public class TransferEntropyCalculatorSpikingIntegration implements (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]; + } + } + return max_val; + } + /* (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++) { + + 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 = max_neighbour_distance(nnPQJointSpikes); + double radiusJointSamples = max_neighbour_distance(nnPQJointSamples); + double radiusConditioningSpikes = max_neighbour_distance(nnPQConditioningSpikes); + double radiusConditioningSamples = max_neighbour_distance(nnPQConditioningSamples); + + currentSum += (- Math.log(radiusJointSpikes) + Math.log(radiusJointSamples) + + Math.log(radiusConditioningSpikes) - Math.log(radiusConditioningSamples)); + } + currentSum /= targetEmbeddingsFromSpikes.size(); + System.out.println("New estimate " + currentSum); int numberOfEvents = eventTypeLocator.size(); diff --git a/tester.py b/tester.py index 6c9f084..ab1824b 100755 --- a/tester.py +++ b/tester.py @@ -42,10 +42,12 @@ sourceArray = [random.normalvariate(0,1) for r in range(numObservations)] destArray = [0] + [sum(pair) for pair in zip([covariance*y for y in sourceArray[0:numObservations-1]], \ [(1-covariance)*y for y in [random.normalvariate(0,1) for r in range(numObservations-1)]] ) ] -sourceArray = 1e3*np.random.random(int(1e3)) +sourceArray = 1e5*np.random.random(int(1e5)) sourceArray.sort() -destArray = 1e3*np.random.random(int(1e3)) -destArray.sort() +#destArray = 1e5*np.random.random(int(1e5)) +#destArray.sort() +destArray = sourceArray + 1 +destArray += np.random.normal(scale = 0.01, size = destArray.shape) # Uncorrelated source array: sourceArray2 = [random.normalvariate(0,1) for r in range(numObservations)]