This commit is contained in:
David Shorten 2022-02-22 11:55:09 +11:00
commit b4e32727c9
31 changed files with 1032 additions and 1232 deletions

View File

@ -14,11 +14,11 @@ JIDT is distributed under the [GNU GPL v3 license](http://www.gnu.org/licenses/g
# Getting started
1. [Download](../../wiki/Downloads) and [Installation](../../wiki/Installation) is very easy!
1. _Quick start_: download the latest [v1.4 full distribution](http://lizier.me/joseph/software/jidt/download.php?file=infodynamics-dist-1.4.zip) (suitable for all platforms) and see the readme.txt file therein.
1. [Documentation](../../wiki/Documentation) including: the paper describing JIDT at [arXiv:1408.3270](http://arxiv.org/abs/1408.3270) (distributed with the toolkit), a [Tutorial](../../wiki/Tutorial), and [Javadocs (v1.4 here)](http://lizier.me/joseph/software/jidt/javadocs/v1.4/);
1. [Demos](../../wiki/Demos) are included with the full distribution, including a [GUI app](../../wiki/AutoAnalyser) for automatic analysis and code generation (see picture below), [simple java demos](../../wiki/SimpleJavaExamples) and [cellular automata (CA) demos](../../wiki/CellularAutomataDemos).
1. These Java tools can easily be used in [Matlab/Octave](../../wiki/OctaveMatlabExamples), [Python](../../wiki/PythonExamples), [R](../../wiki/R_Examples), [Julia](../../wiki/JuliaExamples) and [Clojure](../../wiki/Clojure_Examples)! (click on each language here for examples)
1. [Download](../../wiki/Downloads) and [Installation](../../wiki/Installation) is very easy!
1. _Quick start_: take a `git clone` or download the latest [v1.5 full distribution](http://lizier.me/joseph/software/jidt/download.php?file=infodynamics-dist-1.5.zip) (suitable for all platforms) and see the readme.txt file therein.
1. [Documentation](../../wiki/Documentation) including: the paper describing JIDT at [arXiv:1408.3270](http://arxiv.org/abs/1408.3270) (distributed with the toolkit), a [Tutorial](../../wiki/Tutorial), and [Javadocs (v1.5 here)](http://lizier.me/joseph/software/jidt/javadocs/v1.5/);
1. [Demos](../../wiki/Demos) are included with the full distribution, including a [GUI app](../../wiki/AutoAnalyser) for automatic analysis and code generation (see picture below), [simple java demos](../../wiki/SimpleJavaExamples) and [cellular automata (CA) demos](../../wiki/CellularAutomataDemos).
1. These Java tools can easily be used in [Matlab/Octave](../../wiki/OctaveMatlabExamples), [Python](../../wiki/PythonExamples), [R](../../wiki/R_Examples), [Julia](../../wiki/JuliaExamples) and [Clojure](../../wiki/Clojure_Examples)! (click on each language here for examples)
![Computing in the GUI app image](https://raw.githubusercontent.com/jlizier/jidt/master/web/AutoAnalyserGUI-2-Compute.png)
@ -39,6 +39,16 @@ See other [PublicationsUsingThisToolkit](../../wiki/PublicationsUsingThisToolkit
# News
_26/11/2018_ - New jar and full distribution files available for **release v1.5**; Changes for v1.5 include:
Added GPU (cuda) capability for KSG Conditional Mutual Information calculator (proper documentation to come), brief [wiki page](../../wiki/GPU) and unit tests included;
Added auto-embedding for TE/AIS with multivariate KSG, and univariate and multivariate Gaussian estimator (plus unit tests), for Ragwitz criteria and Maximum bias-corrected AIS, and also added Maximum bias corrected AIS and TE to handle source embedding as well;
Kozachenko entropy estimator adds noise to data by default;
Added bias-correction property to Gaussian and Kernel estimators for MI and conditional MI, including with surrogates (only option for kernel);
Enabled use of different bases for different variables in MI discrete estimator;
All new above features enabled in AutoAnalyser;
Added drop-down menus for parameters in AutoAnalyser;
Included long-form lecture slides in course folder;
_26/11/2017_ - New jar and full distribution files available for **release v1.4**; Changes for v1.4 include:
Major expansion of functionality for AutoAnalysers: adding Launcher applet and capability to double click jar to launch, added Entropy, CMI, CTE and AIS AutoAnalysers, also added binned estimator type, added all variables/pairs analysis, added statistical significance analysis, and ensured functionality of generated Python code with Python3;
Added GPU (cuda) capability for KSG Mutual Information calculator (proper documentation and wiki page to come), including unit tests;

View File

@ -5,12 +5,13 @@
</description>
<!-- set global properties for this build -->
<property name="version" value="1.4"/>
<property name="version" value="1.5"/>
<property name="mainfilename" value="infodynamics"/>
<property name="jarplainname" value="${mainfilename}.jar" />
<property name="jarversiondistnamezip" value="${mainfilename}-jar-${version}.zip" />
<property name="distname" value="${mainfilename}-dist-${version}" />
<property name="distnamezip" value="${distname}.zip" />
<property name="distnametargz" value="${distname}.tar.gz" />
<property name="src" location="java/source"/>
<property name="cudasrc" location="cuda"/>
<property name="bin" location="bin"/>
@ -148,6 +149,7 @@
<delete file="demos/clojure/project.clj"/>
<delete file="${jarversiondistnamezip}"/>
<delete file="${distnamezip}"/>
<delete file="${distnametargz}"/>
<delete file="${jarplainname}"/>
<delete>
<fileset dir="demos/AutoAnalyser" includes="GeneratedCalculator.*"/>
@ -293,24 +295,30 @@
</zip>
</target>
<!-- Developer build - builds everything and makes the full distribution file -->
<!-- Developer build - builds everything and makes the full distribution file in zip and tar.gz -->
<target name="dist" depends="jar,junit,javadocs,readmefiles" description="developer: generate the full distribution">
<echo message="${ant.project.name}: ${ant.file}"/>
<zip destfile="${distnamezip}">
<fileset file="build.xml"/>
<fileset file="${jarplainname}"/>
<zipfileset file="${jarplainname}" filemode="755"/>
<fileset file="license-gplv3.txt"/>
<fileset file="readme.txt"/>
<fileset file="InfoDynamicsToolkit.pdf"/>
<fileset file="JIDT-logo.png" erroronmissingdir="false"/> <!-- This file is missing in full repository versions -->
<fileset file="${versionfile}"/>
<zipfileset dir="java" includes="**/*.java" prefix="java"/>
<zipfileset dir="demos" includes="**/*.*,**/*" excludes="clojure/deploy,clojure/deploy/*.*,python/*.pyc" prefix="demos"/>
<zipfileset dir="demos" includes="**/*.*,**/*" excludes="clojure/deploy,clojure/deploy/*.*,python/*.pyc,**/*.sh,**/*.bat" prefix="demos"/>
<zipfileset dir="demos" includes="**/*.sh,**/*.bat" prefix="demos" filemode="755"/> <!-- Do these separately to get executable permissions -->
<zipfileset dir="javadocs" includes="**/*.*,**/*" prefix="javadocs"/>
<zipfileset dir="notices" includes="**/*.*,**/*" prefix="notices"/>
<zipfileset dir="cuda" prefix="cuda" excludes="benchmark.sh"/>
<zipfileset dir="cuda" prefix="cuda" includes="benchmark.sh" filemode="755"/> <!-- Do this separately to get executable permissions -->
<zipfileset dir="course" prefix="course"/>
<zipfileset dir="tutorial" prefix="tutorial"/> <!-- Get rid of this when tutorial is subsumed in course... -->
<zipfileset dir="web" includes="JIDT-logo.png" prefix="" erroronmissingdir="false"/> <!-- This file is missing in zip dist versions -->
</zip>
<tar destfile="${distnametargz}" compression="gzip" longfile="posix"> <!-- for longfiles could also use "gnu" but apparently is slightly less widely supported -->
<zipfileset src="${distnamezip}"/>
</tar>
</target>
</project>

BIN
cuda/README-GPU.pdf Normal file

Binary file not shown.

View File

@ -1,5 +0,0 @@
@ECHO OFF
REM Run the example:
java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserAIS

View File

@ -1,5 +0,0 @@
#!/bin/bash
# Run the example:
java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserAIS

View File

@ -1,5 +0,0 @@
@ECHO OFF
REM Run the example:
java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCMI

View File

@ -1,5 +0,0 @@
#!/bin/bash
# Run the example:
java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCMI

View File

@ -1,5 +0,0 @@
@ECHO OFF
REM Run the example:
java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCTE

View File

@ -1,5 +0,0 @@
#!/bin/bash
# Run the example:
java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserCTE

View File

@ -1,5 +0,0 @@
@ECHO OFF
REM Run the example:
java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserEntropy

View File

@ -1,5 +0,0 @@
#!/bin/bash
# Run the example:
java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserEntropy

View File

@ -1,5 +0,0 @@
@ECHO OFF
REM Run the example:
java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserMI

View File

@ -1,5 +0,0 @@
#!/bin/bash
# Run the example:
java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserMI

View File

@ -1,5 +0,0 @@
@ECHO OFF
REM Run the example:
java -classpath "..\..\infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserTE

View File

@ -1,5 +0,0 @@
#!/bin/bash
# Run the example:
java -classpath "../../infodynamics.jar" infodynamics.demos.autoanalysis.AutoAnalyserTE

Binary file not shown.

0
demos/java/example10GPUBenchmark.sh Normal file → Executable file
View File

View File

@ -67,9 +67,10 @@ public class Example10GPUBenchmark {
miCalc.setProperty("k", "4");
miCalc.initialise(src[0].length, tgt[0].length);
miCalc.setObservations(src, tgt);
int nb_surrogates = 150;
double[] timeAndValue = new double[2];
long startTime = System.nanoTime();
timeAndValue[1] = miCalc.computeAverageLocalOfObservations();
timeAndValue[1] = miCalc.computeSignificance(nb_surrogates).actualValue;
timeAndValue[0] = (System.nanoTime() - startTime)/1000000.0;
return timeAndValue;
}

Binary file not shown.

View File

@ -717,9 +717,10 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
System.load(gpuLibraryPath);
}
} catch (Throwable e) {
String errmsg = "GPU library not found. To compile GPU code set the enablegpu flag to true in build.xml";
String errmsg = "GPU library not found. To compile GPU code set the enablegpu flag to true in build.xml, or run `ant gpu jar`.";
errmsg += "\nFor more information see the JIDT GPU wiki page: https://github.com/jlizier/jidt/wiki/GPU";
if (gpuLibraryPath.length() > 0) {
errmsg += "\nGPU library was not found in the path provided. Provide full path including library file name.";
errmsg += "\n\nGPU library was not found in the path provided. Provide full path including library file name.";
errmsg += "\nExample: /home/johndoe/myfolder/libKraskov.so";
}
throw new Exception(errmsg);

View File

@ -282,68 +282,11 @@ public interface TransferEntropyCalculatorSpiking {
* will vary depending on the underlying implementation
*/
public SpikingLocalInformationValues computeLocalOfPreviousObservations() throws Exception;
/**
* Generate a bootstrapped distribution of what the TE would look like,
* under a null hypothesis that the source values of our
* samples had no temporal relation to the destination value.
*
* <p>See Section II.E "Statistical significance testing" of
* the JIDT paper below for a description of how this is done for MI and TE in general.
* </p>
*
* <p>Note that if several disjoint time-series have been added
* as observations using {@link #addObservations(double[])} etc.,
* then these separate "trials" will be mixed up in the generation
* of surrogates here.</p>
*
* <p>This method (in contrast to {@link #computeSignificance(int[][])})
* creates <i>random</i> shufflings of the source embedding vectors for the surrogate
* calculations.</p>
*
* @param numPermutationsToCheck number of surrogate samples to bootstrap
* to generate the distribution.
* @return the distribution of TE scores under this null hypothesis.
* @see "J.T. Lizier, 'JIDT: An information-theoretic
* toolkit for studying the dynamics of complex systems', 2014."
*/
public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck) throws Exception;
/**
* Generate a bootstrapped distribution of what the TE would look like,
* under a null hypothesis that the source values of our
* samples had no relation to the destination value.
*
* <p>See Section II.E "Statistical significance testing" of
* the JIDT paper below for a description of how this is done for MI and TE.
* </p>
*
* <p>Note that if several disjoint time-series have been added
* as observations using {@link #addObservations(double[])} etc.,
* then these separate "trials" will be mixed up in the generation
* of surrogates here.</p>
*
* <p>This method (in contrast to {@link #computeSignificance(int)})
* allows the user to specify how to construct the surrogates,
* such that repeatable results may be obtained.</p>
*
* @param newOrderings a specification of how to shuffle the source embedding vectors
* between all of our samples to create the surrogates to generate the distribution with. The first
* index is the permutation number (i.e. newOrderings.length is the number
* of surrogate samples we use to bootstrap to generate the distribution here.)
* Each array newOrderings[i] should be an array of length N (where
* would be the value returned by {@link #getNumObservations()}),
* containing a permutation of the values in 0..(N-1).
* TODO Need to think this through a little more before implementing.
* @return the distribution of channel measure scores under this null hypothesis.
* @see "J.T. Lizier, 'JIDT: An information-theoretic
* toolkit for studying the dynamics of complex systems', 2014."
* @throws Exception where the length of each permutation in newOrderings
* is not equal to the number N samples that were previously supplied.
*/
public EmpiricalMeasurementDistribution computeSignificance(
int[][] newOrderings) throws Exception;
public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck, double estimatedValue) throws Exception;
public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck,
double estimatedValue, long randomSeed) throws Exception;
/**
* Set or clear debug mode for extra debug printing to stdout
*

View File

@ -723,33 +723,33 @@ public class TransferEntropyCalculatorSpikingIntegrationWindows2 implements
//kdTreeDestNext.setNormType(3);
nnPQ =
kdTreeDestNext.findKNearestNeighbours(
Knns, eventIndexWithinType);
// nnPQ =
// kdTreeDestNext.findKNearestNeighbours(
// Knns, eventIndexWithinType);
// Find eps_{x,y,z} as the maximum x, y and z norms amongst this set:
radius_sourcePast = 0.0;
radius_destPast = 0.0;
radius_destNext = 0.0;
radius_destNext_sampleIndex = -1;
radius_max = 0;
for (int j = 0; j < Knns; j++) {
// Take the furthest remaining of the nearest neighbours from the PQ:
NeighbourNodeData nnData = nnPQ.poll();
/*if (nnData.norms[0] > radius_sourcePast) {
radius_sourcePast = nnData.norms[0];
}*/
if (nnData.norms[0] > radius_destPast) {
radius_destPast = nnData.norms[0];
}
if (nnData.norms[1] > radius_destNext) {
radius_destNext = nnData.norms[1];
radius_destNext_sampleIndex = nnData.sampleIndex;
}
if(nnData.distance > radius_max) {
radius_max = nnData.distance;
}
}
// // Find eps_{x,y,z} as the maximum x, y and z norms amongst this set:
// radius_sourcePast = 0.0;
// radius_destPast = 0.0;
// radius_destNext = 0.0;
// radius_destNext_sampleIndex = -1;
// radius_max = 0.0;
// for (int j = 0; j < Knns; j++) {
// // Take the furthest remaining of the nearest neighbours from the PQ:
// NeighbourNodeData nnData = nnPQ.poll();
// /*if (nnData.norms[0] > radius_sourcePast) {
// radius_sourcePast = nnData.norms[0];
// }*/
// if (nnData.norms[0] > radius_destPast) {
// radius_destPast = nnData.norms[0];
// }
// if (nnData.norms[1] > radius_destNext) {
// radius_destNext = nnData.norms[1];
// radius_destNext_sampleIndex = nnData.sampleIndex;
// }
// if(nnData.distance > radius_max) {
// radius_max = nnData.distance;
// }
// }
//radius_max = Math.max(radius_destPast, radius_destNext);
@ -834,8 +834,8 @@ public class TransferEntropyCalculatorSpikingIntegrationWindows2 implements
double rawRateGivenSourceAndDest = (double) Knns / timeInWindowWithMatchingJointHistories;
double rawRateGivenDest = (double) countOfDestNextMatched / timeInWindowWithMatchingDestHistory;
// Attempt at bias correction:
double logRateGivenSourceAndDestCorrected = -Math.log(timeInWindowWithMatchingJointHistories);
double logRateGivenDestCorrected = -Math.log(timeInWindowWithMatchingDestHistory);
double logRateGivenSourceAndDestCorrected = MathsUtils.digamma(Knns) - Math.log(timeInWindowWithMatchingJointHistories);
double logRateGivenDestCorrected = MathsUtils.digamma(countOfDestNextMatched) - Math.log(timeInWindowWithMatchingDestHistory);
//System.out.println(logRateGivenSourceAndDestCorrected + " " + logRateGivenDestCorrected);
//============================

17
java/source/infodynamics/utils/KdTree.java Executable file → Normal file
View File

@ -275,10 +275,10 @@ public class KdTree extends NearestNeighbourSearcher {
// Could remove these since the code is now functional,
// but may be better to leave them in just in case the code breaks:
if (leftIndex > leftStart + leftNumPoints) {
throw new RuntimeException("Exceeded expected number of points on left");
throw new RuntimeException("Exceeded expected number of points on left - likely had an NaN in the data");
}
if (rightIndex > rightStart + rightNumPoints) {
throw new RuntimeException("Exceeded expected number of points on right");
throw new RuntimeException("Exceeded expected number of points on right - likely had an NaN in the data");
}
// Update the pointer for the sorted indices for this dimension,
// and keep the new temporary array
@ -1384,19 +1384,19 @@ public class KdTree extends NearestNeighbourSearcher {
}
@Override
public void findPointsWithinR(int sampleIndex,
public int findPointsWithinR(int sampleIndex,
double r, boolean allowEqualToR, boolean[] isWithinR,
int[] indicesWithinR) {
if (rootNode == null) {
indicesWithinR[0] = -1;
return;
return 0;
}
int currentIndexInIndicesWithinR =
findPointsWithinR(sampleIndex,
rootNode, 0, r, allowEqualToR, isWithinR,
indicesWithinR, 0);
indicesWithinR[currentIndexInIndicesWithinR] = -1;
return;
return currentIndexInIndicesWithinR;
}
/**
@ -2510,20 +2510,21 @@ public class KdTree extends NearestNeighbourSearcher {
* @param allowEqualToR if true, then count points at radius r also, otherwise only those strictly within r
* @param isWithinR the array MUST be passed in with all points set to false initially, and is returned indicating whether each sample was found to be within r of that at sampleIndex.
* @param indicesWithinR a list of array indices for points marked as true in isWithinR, terminated with a -1 value.
* @return the number of points found to match
*/
public void findPointsWithinRs(int sampleIndex,
public int findPointsWithinRs(int sampleIndex,
double[] rs, int dynCorrExclTime, boolean allowEqualToR,
boolean[] isWithinR, int[] indicesWithinR) {
if (rootNode == null) {
indicesWithinR[0] = -1;
return;
return 0;
}
int currentIndexInIndicesWithinR =
findPointsWithinRs(sampleIndex,
rootNode, 0, rs, dynCorrExclTime, allowEqualToR, isWithinR,
indicesWithinR, 0);
indicesWithinR[currentIndexInIndicesWithinR] = -1;
return;
return currentIndexInIndicesWithinR;
}
/**

View File

@ -326,8 +326,9 @@ public abstract class NearestNeighbourSearcher {
* found to be within r of that at sampleIndex.
* @param indicesWithinR a list of array indices
* for points marked as true in isWithinR, terminated with a -1 value.
* @return the number of matching points found
*/
public abstract void findPointsWithinR(
public abstract int findPointsWithinR(
int sampleIndex, double r,
boolean allowEqualToR, boolean[] isWithinR, int[] indicesWithinR);

View File

@ -685,7 +685,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
* </p>
*
*/
public void findPointsWithinR(
public int findPointsWithinR(
int sampleIndex, double r, boolean allowEqualToR,
boolean[] isWithinR, int[] indicesWithinR) {
int indexInIndicesWithinR = 0;
@ -718,7 +718,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
break;
}
// Write the terminating integer into the indicesWithinR array:
indicesWithinR[indexInIndicesWithinR++] = -1;
indicesWithinR[indexInIndicesWithinR] = -1;
return indexInIndicesWithinR;
}
/**

View File

@ -86,14 +86,14 @@ public class GPUPerformanceTester extends TestCase {
startTime = Calendar.getInstance().getTimeInMillis();
miCalc.initialise(1,1);
miCalc.setObservations(source, dest);
cpu_val = miCalc.computeSignificance(nb_surrogates).getMeanOfDistribution();
cpu_val = miCalc.computeSignificance(nb_surrogates).actualValue;
cpu_duration = Calendar.getInstance().getTimeInMillis() - startTime;
miCalc.setProperty("USE_GPU", "true");
miCalc.initialise(1,1);
startTime = Calendar.getInstance().getTimeInMillis();
miCalc.setObservations(source, dest);
gpu_val = miCalc.computeSignificance(nb_surrogates).getMeanOfDistribution();
gpu_val = miCalc.computeSignificance(nb_surrogates).actualValue;
gpu_duration = Calendar.getInstance().getTimeInMillis() - startTime;
assertEquals(cpu_val, gpu_val, 0.0001);
@ -124,12 +124,12 @@ public class GPUPerformanceTester extends TestCase {
return;
}
int timeSteps = 10000;
int timeSteps = 1000;
RandomGenerator rg = new RandomGenerator();
double[] source = rg.generateNormalData(timeSteps, 0, 1);
double[] dest = rg.generateNormalData(timeSteps, 0, 1);
compareGPUPerformance(source, dest, "Random low-dimensional data");
compareGPUPerformanceSurrogates(source, dest, 150, "Random low-dimensional data");
return;
}

View File

@ -134,6 +134,18 @@ Notices for this software are found in the notices/JAMA directory.
Release notes
===============
v1.5 26/11/2018
---------------
(after 753 commits recorded by github, repository as at https://github.com/jlizier/jidt/tree/603445651cc0bf155a42c9ba336141bc7f29bccd save for this file update)
Added GPU (cuda) capability for KSG Conditional Mutual Information calculator (proper documentation to come), including unit tests and brief wiki page;
Added auto-embedding for TE/AIS with multivariate KSG, and univariate and multivariate Gaussian estimator (plus unit tests), for Ragwitz criteria and Maximum bias-corrected AIS, and also added Maximum bias corrected AIS and TE to handle source embedding as well;
Kozachenko entropy estimator adds noise to data by default;
Added bias-correction property to Gaussian and Kernel estimators for MI and conditional MI, including with surrogates (only option for kernel);
Enabled use of different bases for different variables in MI discrete estimator;
All new above features enabled in AutoAnalyser;
Added drop-down menus for parameters in AutoAnalyser;
Included long-form lecture slides in course folder;
v1.4 26/11/2017
---------------
(after 638 commits recorded by github, repository as at https://github.com/jlizier/jidt/tree/589d51674e6a9cfb569432679e515bea17092876 save for this file update)

186
tester.py Executable file
View File

@ -0,0 +1,186 @@
##
## Java Information Dynamics Toolkit (JIDT)
## Copyright (C) 2012, Joseph T. Lizier
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
##
# = Example 4 - Transfer entropy on continuous data using Kraskov estimators =
# Simple transfer entropy (TE) calculation on continuous-valued data using the Kraskov-estimator TE calculator.
from jpype import *
import random
import math
import os
import numpy as np
NUM_REPS = 2
NUM_SPIKES = int(2e3)
NUM_OBSERVATIONS = 2
NUM_SURROGATES = 10
# Params for canonical example generation
RATE_Y = 1.0
RATE_X_MAX = 10
def generate_canonical_example_processes(num_y_events):
event_train_x = []
event_train_x.append(0)
event_train_y = np.random.uniform(0, int(num_y_events / RATE_Y), int(num_y_events))
event_train_y.sort()
most_recent_y_index = 0
previous_x_candidate = 0
while most_recent_y_index < (len(event_train_y) - 1):
this_x_candidate = previous_x_candidate + random.expovariate(RATE_X_MAX)
while most_recent_y_index < (len(event_train_y) - 1) and this_x_candidate > event_train_y[most_recent_y_index + 1]:
most_recent_y_index += 1
delta_t = this_x_candidate - event_train_y[most_recent_y_index]
rate = 0
if delta_t > 1:
rate = 0.5
else:
rate = 0.5 + 5.0 * math.exp(-50 * (delta_t - 0.5)**2) - 5.0 * math.exp(-50 * (0.5)**2)
if random.random() < rate/float(RATE_X_MAX):
event_train_x.append(this_x_candidate)
previous_x_candidate = this_x_candidate
event_train_x.sort()
event_train_y.sort()
return event_train_x, event_train_y
# Change location of jar to match yours (we assume script is called from demos/python):
jarLocation = os.path.join(os.getcwd(), "infodynamics.jar");
if (not(os.path.isfile(jarLocation))):
exit("infodynamics.jar not found (expected at " + os.path.abspath(jarLocation) + ") - are you running from demos/python?")
# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
teCalcClass = JPackage("infodynamics.measures.spiking.integration").TransferEntropyCalculatorSpikingIntegration
teCalc = teCalcClass()
teCalc.setProperty("knns", "4")
print("Independent Poisson Processes")
teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
teCalc.setProperty("SOURCE_PAST_INTERVALS", "1,2")
teCalc.appendConditionalIntervals(JArray(JInt, 1)([1, 2]))
teCalc.appendConditionalIntervals(JArray(JInt, 1)([1, 2]))
teCalc.setProperty("NORM_TYPE", "MAX_NORM")
results_poisson = np.zeros(NUM_REPS)
for i in range(NUM_REPS):
teCalc.startAddObservations()
for j in range(NUM_OBSERVATIONS):
sourceArray = NUM_SPIKES*np.random.random(NUM_SPIKES)
sourceArray.sort()
destArray = NUM_SPIKES*np.random.random(NUM_SPIKES)
destArray.sort()
condArray = NUM_SPIKES*np.random.random((2, NUM_SPIKES))
condArray.sort(axis = 1)
teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray), JArray(JDouble, 2)(condArray))
teCalc.finaliseAddObservations();
result = teCalc.computeAverageLocalOfObservations()
print("TE result %.4f nats" % (result,))
sig = teCalc.computeSignificance(NUM_SURROGATES, result)
print(sig.pValue)
results_poisson[i] = result
print("Summary: mean ", np.mean(results_poisson), " std dev ", np.std(results_poisson))
teCalc = teCalcClass()
teCalc.setProperty("knns", "4")
print("Noisy copy zero TE")
#teCalc.appendConditionalIntervals(JArray(JInt, 1)([1]))
teCalc.setProperty("DEST_PAST_INTERVALS", "1")
teCalc.setProperty("SOURCE_PAST_INTERVALS", "1")
#teCalc.setProperty("NORM_TYPE", "MAX_NORM")
results_noisy_zero = np.zeros(NUM_REPS)
for i in range(NUM_REPS):
teCalc.startAddObservations()
for j in range(NUM_OBSERVATIONS):
condArray = np.ones((1, NUM_SPIKES)) + 0.05 * np.random.random((1, NUM_SPIKES))
condArray = np.cumsum(condArray, axis = 1)
condArray.sort(axis = 1)
sourceArray = condArray[0, :] + 0.25 + 0.05 * np.random.normal(size = condArray.shape[1])
sourceArray.sort()
destArray = condArray[0, :] + 0.5 + 0.05 * np.random.normal(size = condArray.shape[1])
destArray.sort()
#teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray), JArray(JDouble, 2)(condArray))
teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray))
teCalc.finaliseAddObservations();
result = teCalc.computeAverageLocalOfObservations()
print("TE result %.4f nats" % (result,))
sig = teCalc.computeSignificance(NUM_SURROGATES, result)
print(sig.pValue)
results_noisy_zero[i] = result
print("Summary: mean ", np.mean(results_noisy_zero), " std dev ", np.std(results_noisy_zero))
teCalc = teCalcClass()
teCalc.setProperty("knns", "4")
print("Noisy copy non-zero TE")
teCalc.appendConditionalIntervals(JArray(JInt, 1)([1]))
teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
teCalc.setProperty("SOURCE_PAST_INTERVALS", "1")
#teCalc.setProperty("NORM_TYPE", "MAX_NORM")
results_noisy_non_zero = np.zeros(NUM_REPS)
for i in range(NUM_REPS):
teCalc.startAddObservations()
for j in range(NUM_OBSERVATIONS):
sourceArray = np.ones((1, NUM_SPIKES)) + 0.05 * np.random.random((1, NUM_SPIKES))
sourceArray = np.cumsum(sourceArray)
sourceArray.sort()
condArray = sourceArray + 0.25 + 0.05 * np.random.normal(size = sourceArray.shape)
condArray.sort()
condArray = np.expand_dims(condArray, 0)
destArray = sourceArray + 0.5 + 0.05 * np.random.normal(size = sourceArray.shape)
destArray.sort()
teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray), JArray(JDouble, 2)(condArray))
teCalc.finaliseAddObservations();
result = teCalc.computeAverageLocalOfObservations()
print("TE result %.4f nats" % (result,))
sig = teCalc.computeSignificance(NUM_SURROGATES, result)
print(sig.pValue)
results_noisy_non_zero[i] = result
print("Summary: mean ", np.mean(results_noisy_non_zero), " std dev ", np.std(results_noisy_zero))
print("Canonical example")
teCalc = teCalcClass()
teCalc.setProperty("knns", "4")
teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
teCalc.setProperty("SOURCE_PAST_INTERVALS", "1")
#teCalc.setProperty("NUM_SAMPLES_MULTIPLIER", "1")
#teCalc.setProperty("NORM_TYPE", "MAX_NORM")
results_canonical = np.zeros(NUM_REPS)
for i in range(NUM_REPS):
event_train_x, event_train_y = generate_canonical_example_processes(NUM_SPIKES)
teCalc.setObservations(JArray(JDouble, 1)(event_train_y), JArray(JDouble, 1)(event_train_x))
result = teCalc.computeAverageLocalOfObservations()
results_canonical[i] = result
print("TE result %.4f nats" % (result,))
sig = teCalc.computeSignificance(NUM_SURROGATES, result)
print(sig.pValue)
print("Summary: mean ", np.mean(results_canonical), " std dev ", np.std(results_canonical))