diff --git a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegrationWindows2.java b/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegrationWindows2.java deleted file mode 100644 index 64f25d6..0000000 --- a/java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegrationWindows2.java +++ /dev/null @@ -1,1282 +0,0 @@ -package infodynamics.measures.spiking.integration; - -import java.util.Arrays; -import java.util.Iterator; -import java.util.PriorityQueue; -import java.util.Random; -import java.util.Vector; - -import infodynamics.measures.spiking.TransferEntropyCalculatorSpiking; -import infodynamics.utils.EmpiricalMeasurementDistribution; -import infodynamics.utils.KdTree; -import infodynamics.utils.MathsUtils; -import infodynamics.utils.MatrixUtils; -import infodynamics.utils.NeighbourNodeData; -import infodynamics.utils.FirstIndexComparatorDouble; -import infodynamics.utils.UnivariateNearestNeighbourSearcher; - -/** - * 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}

- * - * @author Joseph Lizier (email, - * www) - */ -public class TransferEntropyCalculatorSpikingIntegrationWindows2 implements - TransferEntropyCalculatorSpiking { - - /** - * Number of past destination spikes to consider (akin to embedding length) - */ - protected int k = 1; - /** - * Number of past source spikes to consider (akin to embedding length) - */ - protected int l = 1; - - /** - * 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. - */ - protected Vector vectorOfSourceSpikeTimes = null; - - /** - * 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 - */ - 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; - - /** - * KdTrees for searching the joint past spaces and time to next spike, - * for each possibility of which spiked next - */ - protected KdTree[] kdTreesJoint = null; - - /** - * KdTrees for searching the joint past spaces, - * for each possibility of which spiked next - */ - protected KdTree[] kdTreesSourceDestHistories = 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) - */ - public static final String TRIM_TO_POS_PROP_NAME = "TRIM_RANGE_TO_POS_TIMES"; - /** - * 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"; - /** - * 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; - - protected boolean trimToPosNextSpikeTimes = false; - - /** - * Stores whether we are in debug mode - */ - protected boolean debug = false; - - public TransferEntropyCalculatorSpikingIntegrationWindows2() { - super(); - } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(int) - */ - @Override - public void initialise() throws Exception { - initialise(k,l); - } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#initialise(int) - */ - @Override - public void initialise(int k) throws Exception { - initialise(k,this.l); - } - - /* (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) - */ - @Override - 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); - } 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_ADD_NOISE)) { - if (propertyValue.equals("0") || - propertyValue.equalsIgnoreCase("false")) { - addNoise = false; - noiseLevel = 0; - } else { - addNoise = true; - noiseLevel = Double.parseDouble(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); - } - } - - /* (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)) { - 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 { - // No property matches for this class - return null; - } - } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#setObservations(double[], double[]) - */ - @Override - public void setObservations(double[] source, double[] destination) - throws Exception { - startAddObservations(); - addObservations(source, destination); - finaliseAddObservations(); - } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#startAddObservations() - */ - @Override - public void startAddObservations() { - vectorOfSourceSpikeTimes = new Vector(); - vectorOfDestinationSpikeTimes = new Vector(); - } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#addObservations(double[], double[]) - */ - @Override - 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() - */ - @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(); - } - destPastAndNextTimings = new Vector(); - eventTypeLocator = new Vector(); - eventIndexLocator = new Vector(); - numEventsPerObservationSet = new Vector(); - - // Send all of the observations through: - Iterator sourceIterator = vectorOfSourceSpikeTimes.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++; - } - // 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}); - } - // 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++; - } - 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); - } - } - - protected void processEventsFromSpikingTimeSeries(double[] sourceSpikeTimes, double[] destSpikeTimes, - int timeSeriesIndex, Vector[] eventTimings, - Vector destPastAndNextTimings, Vector eventTypeLocator, - Vector eventIndexLocator, Vector numEventsPerObservationSet) throws Exception { - // addObservationsAfterParamsDetermined(sourceSpikeTimes, destSpikeTimes); - - // 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); - } - } - // Post-condition: dest_index and source_index are set correctly for the first set of pasts - - double timeToNextSpike; - boolean nextIsDest = false; - double[] spikeTimesForNextSpiker; - double timeOfPrevDestSpike = destSpikeTimes[dest_index]; - int numEvents = 0; - Random random = null; - if (addNoise) { - random = new Random(); - } - 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; - if (addNoise) { - timeToNextSpike += random.nextGaussian()*noiseLevel; - } - // 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]; - if (addNoise) { - sourcePast[0] += random.nextGaussian()*noiseLevel; - } - for (int i = 1; i < k; i++) { - destPast[i - 1] = destSpikeTimes[dest_index - i + 1] - - destSpikeTimes[dest_index - i]; - if (addNoise) { - destPast[i - 1] += random.nextGaussian()*noiseLevel; - } - } - for (int i = 1; i < l; i++) { - sourcePast[i] = sourceSpikeTimes[source_index - i + 1] - - sourceSpikeTimes[source_index - i]; - if (addNoise) { - sourcePast[i] += random.nextGaussian()*noiseLevel; - } - } - // 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); - } - } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getAddedMoreThanOneObservationSet() - */ - @Override - public boolean getAddedMoreThanOneObservationSet() { - return (vectorOfDestinationSpikeTimes != null) && - (vectorOfDestinationSpikeTimes.size() > 1); - } - - /* (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")); - } - - // Select only events where the destination spiked next: - if (eventType != NEXT_DEST) { - // Pre-condition: next event is a source spike so we'll continue to check next event - if (debug && (eventIndex < 10000)) { - System.out.println(); - } - continue; - } - // Post-condition: the next event is a destination spike: - - // Find the Knns nearest neighbour matches to this event, - // with the same previous spiker and the next. - // TODO Add dynamic exclusion time later - //kdTreesJoint[eventType].setNormType(3); - 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; - double 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[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(nnData.distance > radius_max) { - radius_max = nnData.distance; - } - } - - //System.out.println(radius_max); - - // Postcondition: radius_* variables hold the search radius for each sourcePast, destPast and destNext matches. - - 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].setNormType(3); - kdTreesSourceDestHistories[NEXT_DEST]. - findPointsWithinRs(eventIndexWithinType, - new double[] {radius_max, radius_max}, 0, - true, isWithinR, indicesWithinR); - //Vector founds = findPointsWithinR(eventIndexWithinType, - // double r, boolean allowEqualToR) { - // And check which of these samples had spike time in dest after ours: - int countOfDestNextAndGreater = 0; - double timeInWindowWithMatchingJointHistories = 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]); - - // We need to check how long we spent in the window matching the next spike - // with a matching history - if (matchedHistoryEventTimings[2][0] >= timeToNextSpikeSincePreviousDestSpike - radius_max - && - (-matchedHistoryEventTimings[0][0]) < timeToNextSpikeSincePreviousDestSpike + radius_max) { - // This sample had a matched history and the next spike was a - // destination spike in the window or after it. - - // Real start of window cannot be before previous destination spike: - double realStartOfWindow = Math.max(timeToNextSpikeSincePreviousDestSpike - radius_max, 0); - // Also, real start cannot be before previous source spike: - // (previous source spike occurs at -matchedHistoryEventTimings[0][0] relative - // to previous destination spike) - realStartOfWindow = Math.max(realStartOfWindow, -matchedHistoryEventTimings[0][0]); - - // Real end of window happened either when the spike occurred or at - // the end of the window: - double realEndOfWindow = Math.min(matchedHistoryEventTimings[2][0], - timeToNextSpikeSincePreviousDestSpike + radius_max); - - // Add in how much time with a matching history we spent in this window: - timeInWindowWithMatchingJointHistories += - realEndOfWindow - realStartOfWindow; - - 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].setNormType(3); - kdTreesSourceDestHistories[NEXT_SOURCE]. - findPointsWithinRs( - new double[] {radius_max, radius_max}, 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] >= timeToNextSpikeSincePreviousDestSpike - radius_max - && - (-matchedHistoryEventTimings[0][0]) < timeToNextSpikeSincePreviousDestSpike + radius_max) { - // 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 - // (i.e. in the window or after it). - // (The "equal to" is why we look for matches within radius_destNext here as well.) - - // Real start of window cannot be before previous destination spike: - // 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.) - - double realStartOfWindow = Math.max(timeToNextSpikeSincePreviousDestSpike - radius_max, 0); - // Also, real start cannot be before previous source spike: - // (previous source spike occurs at -matchedHistoryEventTimings[0][0] relative - // to previous destination spike) - realStartOfWindow = Math.max(realStartOfWindow, -matchedHistoryEventTimings[0][0]); - - // Real end of window happened either when the spike occurred or at - // the end of the window: - double realEndOfWindow = Math.min(matchedHistoryEventTimings[2][0], - timeToNextSpikeSincePreviousDestSpike + radius_max); - - // Add in how much time with a matching history we spent in this window: - timeInWindowWithMatchingJointHistories += - realEndOfWindow - realStartOfWindow; - - countOfSourceNextAndGreater++; - } - // Reset the isWithinR array while we're here - isWithinR[indicesWithinR[nIndex]] = false; - } - - // We need to count spike rate for all the times we're actually within the matching window for the next spike - // This is kind of inspired by the Greg Ver Steeg et al. approach in - // http://www.jmlr.org/proceedings/papers/v38/gao15.pdf - // which is thinking about where the space is actually being explored. - // This is where the window correction code was placed, which we're now replacing - // with computing the length of actual time we spend in the next window. - - - 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; - double timeInWindowWithMatchingDestHistory = 0; - // Real start of window cannot be before previous destination spike: - double realStartOfWindow = 0; - if (k > 1) { - - //kdTreeDestNext.setNormType(3); - - // 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.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); - - //System.out.println(radius_destPast + " " + radius_destNext); - - // 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.setNormType(3); - kdTreeDestHistory.findPointsWithinR(indexForNextIsDest, radius_max, - true, isWithinR, indicesWithinR); - // And check which of these samples had next spike time after our window starts: - 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_max - && - (-matchedHistoryEventTimings[0][0]) < timeToNextSpikeSincePreviousDestSpike + radius_max) { - // 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.) - - realStartOfWindow = Math.max(timeToNextSpikeSincePreviousDestSpike - radius_max, 0); - // Also, real start cannot be before previous source spike: - // (previous source spike occurs at -matchedHistoryEventTimings[0][0] relative - // to previous destination spike) - realStartOfWindow = Math.max(realStartOfWindow, -matchedHistoryEventTimings[0][0]); - - // Real end of window happened either when the spike occurred or at - // the end of the window: - double realEndOfWindow = Math.min(matchedHistoryEventTimings[1][0], - timeToNextSpikeSincePreviousDestSpike + radius_max); - - // Add in how much time with a matching history we spent in this window: - timeInWindowWithMatchingDestHistory += - realEndOfWindow - realStartOfWindow; - - countOfDestNextAndGreaterMatchedDest++; - if (matchedHistoryEventTimings[1][0] <= timeToNextSpikeSincePreviousDestSpike + radius_max) { - // Then we also have a match on the next spike itself - countOfDestNextMatched++; - } - } - // Reset the isWithinR array while we're here - isWithinR[indicesWithinR[nIndex]] = false; - } - //System.out.println(Knns + " " + countOfDestNextMatched + " " + timeInWindowWithMatchingDestHistory); - } else { - - // For k = 1, we only care about time since last spike. - // So count how many of the next spikes were within the window first: - // -- 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, Math.min(radius_destNext, timeToNextSpikeSincePreviousDestSpike), true); - // And also check how long each of these spent in the window: - timeInWindowWithMatchingDestHistory = - nnSearcherDestTimeToNextSpike.sumDistanceAboveThresholdForPointsWithinRs(indexForNextIsDest, - radius_destNext, Math.min(radius_destNext, timeToNextSpikeSincePreviousDestSpike), true); - // Now 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, - Math.min(radius_destNext, timeToNextSpikeSincePreviousDestSpike), true); - // And we need to add time in for all of the (countOfDestNextAndGreaterMatchedDest - countOfDestNextMatched) - // points which didn't spike in the window - timeInWindowWithMatchingDestHistory += (double) (countOfDestNextAndGreaterMatchedDest - countOfDestNextMatched) * - (timeToNextSpikeSincePreviousDestSpike + radius_destNext - realStartOfWindow); - } - - if (debug && (eventIndex < 10000)) { - System.out.printf(", and %d of %d points for D history only; ", - countOfDestNextMatched, countOfDestNextAndGreaterMatchedDest); - } - - //============================ - // This code section takes the counts of spikes and total intervals, and - // estimates the log rates. - // Inferred rates raw: - double rawRateGivenSourceAndDest = (double) Knns / timeInWindowWithMatchingJointHistories; - double rawRateGivenDest = (double) countOfDestNextMatched / timeInWindowWithMatchingDestHistory; - // Attempt at bias correction: - double logRateGivenSourceAndDestCorrected = MathsUtils.digamma(Knns) - Math.log(timeInWindowWithMatchingJointHistories); - double logRateGivenDestCorrected = MathsUtils.digamma(countOfDestNextMatched) - Math.log(timeInWindowWithMatchingDestHistory); - //System.out.println(logRateGivenSourceAndDestCorrected + " " + logRateGivenDestCorrected); - //============================ - - if (debug && (eventIndex < 10000)) { - System.out.printf(" te ~~ %.4f - %.4f = %.4f, log (%.4f)/(%.4f) = %.4f (counts %d/%d = %.4f, %d/%d = %.4f -> te %.4f)\n", - logRateGivenSourceAndDestCorrected, - logRateGivenDestCorrected, - logRateGivenSourceAndDestCorrected - logRateGivenDestCorrected, - rawRateGivenSourceAndDest, - rawRateGivenDest, - Math.log(rawRateGivenSourceAndDest / rawRateGivenDest), - Knns, - Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater, - (double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater), - countOfDestNextMatched, countOfDestNextAndGreaterMatchedDest, - (double) countOfDestNextMatched / (double) countOfDestNextAndGreaterMatchedDest, - Math.log(((double) Knns / (double) (Knns + countOfSourceNextAndGreater + countOfDestNextAndGreater)) / - ((double) (countOfDestNextMatched) / (double) (countOfDestNextAndGreaterMatchedDest)))); - } - - //====================== - // Add the contribution in: - // a.If we were using digamma logs: - contributionFromSpikes += logRateGivenSourceAndDestCorrected - logRateGivenDestCorrected; - // contributionFromSpikes += Math.log(rawRateGivenSourceAndDest / rawRateGivenDest); - - - // 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)); - } - System.out.println("All done!"); - contributionFromSpikes /= totalTimeLength; - return contributionFromSpikes; - } - - /* - * 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); - } - - // 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); - - } - // 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; - } - 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; - } - */ - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#computeLocalOfPreviousObservations() - */ - @Override - 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) - */ - @Override - public void setDebug(boolean debug) { - this.debug = debug; - } - - /* (non-Javadoc) - * @see infodynamics.measures.spiking.TransferEntropyCalculatorSpiking#getLastAverage() - */ - @Override - public double getLastAverage() { - // TODO Auto-generated method stub - return 0; - } - -}