natas_hw
/
jidt
mirror of https://github.com/jlizier/jidt
0
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Compare commits

base: natas_hw:master
Branches Tags
natas_hw:master
natas_hw:course
natas_hw:autoanalyser
natas_hw:gh-pages
natas_hw:working
natas_hw:v1.6.1
natas_hw:v1.6
natas_hw:v1.5
natas_hw:v1.4
natas_hw:v1.3.1
natas_hw:v1.3
..
compare: natas_hw:v1.4
Branches Tags
natas_hw:course
natas_hw:master
natas_hw:autoanalyser
natas_hw:gh-pages
natas_hw:working
natas_hw:v1.6.1
natas_hw:v1.6
natas_hw:v1.5
natas_hw:v1.4
natas_hw:v1.3.1
natas_hw:v1.3

No commits in common. "master" and "v1.4" have entirely different histories.
master ... v1.4

228 changed files with 2918 additions and 33547 deletions
Show Stats Expand all files Collapse all files

52
README.md
View File

@ -1,6 +1,6 @@
# Java Information Dynamics Toolkit (JIDT)
Copyright (C) 2012- [Joseph T. Lizier](http://lizier.me/joseph/); 2014- Ipek Özdemir; 2017- [Pedro Mediano](https://www.doc.ic.ac.uk/~pam213/); 2019- Emanuele Crosato, Sooraj Sekhar, Oscar Huaigu Xu; 2022- [David Shorten](https://scholar.google.com/citations?user=ggF3Lt4AAAAJ&hl=en)
Copyright (C) 2012-2014 [Joseph T. Lizier](http://lizier.me/joseph/); 2014-2016 [Joseph T. Lizier](http://lizier.me/joseph/) and Ipek Özdemir; 2017- [Joseph T. Lizier](http://lizier.me/joseph/), Ipek Özdemir and [Pedro Mediano](https://www.doc.ic.ac.uk/~pam213/)
*JIDT* provides a stand-alone, open-source code Java implementation (also usable in [Matlab, Octave](../../wiki/UseInOctaveMatlab), [Python](../../wiki/UseInPython), [R](../../wiki/UseInR), [Julia](../../wiki/UseInJulia) and [Clojure](../../wiki/UseInClojure)) of information-theoretic measures of distributed computation in complex systems: i.e. information storage, transfer and modification.
@ -8,25 +8,19 @@ JIDT includes implementations:
* principally for the measures **transfer entropy**, **mutual information**, and their conditional variants, as well as **active information storage**, entropy, etc;
* for both _discrete_ and _continuous_-valued data;
* using various types of estimators (e.g. _Kraskov-Stögbauer-Grassberger estimators_, _box-kernel estimation_, _linear-Gaussian_),
as described in full at [ImplementedMeasures](../../wiki/ImplementedMeasures).
JIDT is easy to use:
* It ships with a **GUI application** -- the [AutoAnalyser](../../wiki/AutoAnalyser), see picture below -- to facilitate point-and-click analysis, as well as code template generation for more complex analysis.
* We provide **short video lectures** and corresponding slides in a (beta) [Course](../../wiki/Course) on how to understand using information-theoretic tools to analyse complex systems, and to implement such analysis with JIDT.
as described in full at ImplementedMeasures.
JIDT is distributed under the [GNU GPL v3 license](http://www.gnu.org/licenses/gpl.html) (or later).
# Getting started
1. [Download](../../wiki/Downloads) and [Installation](../../wiki/Installation) is very easy!
1. _Quick start_: take a `git clone` (then build via [AntScripts](../../wiki/AntScripts)) OR download the latest [v1.6.1 full distribution](https://lizier.me/joseph/software/jidt/download.php?file=infodynamics-dist-1.6.1.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 (beta) [Course](../../wiki/Course) including short video lectures and a shorter [Tutorial](../../wiki/Tutorial), and [Javadocs (v1.6.1 here)](https://lizier.me/joseph/software/jidt/javadocs/v1.6.1/);
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_: 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)
[![Computing in the GUI app image](https://raw.githubusercontent.com/jlizier/jidt/master/web/AutoAnalyserGUI-2-Compute.png)](../../wiki/AutoAnalyser)
[![Course and video lectures](https://raw.githubusercontent.com/jlizier/jidt/master/web/YouTubePlaylist.png)](../../wiki/Course)
![Computing in the GUI app image](https://raw.githubusercontent.com/jlizier/jidt/master/web/AutoAnalyserGUI-2-Compute.png)
For further information or announcements:
* Join our discussion group: http://groups.google.com/d/forum/jidt-discuss
@ -45,28 +39,6 @@ See other [PublicationsUsingThisToolkit](../../wiki/PublicationsUsingThisToolkit
# News
_22/08/2023_ - New full distribution files available for **release v1.6.1**; Changes for v1.6.1 include:
Minor updates to supporting use in Python, including virtual environments;
Minor tweaks to fish schooling examples (mostly comments).
_5/09/2022_ - New full distribution files available for **release v1.6**; Changes for v1.6 include:
Adding Flocking/Schooling/Swarming demo;
Included Pedro's code on IIT and O-/S-Information measures;
Spiking TE estimator added from David;
Fixed up AutoAnalyser to work well for Python3 and numpy;
Links to lecture videos included in the beta wiki for the course;
Added rudimentary effective network inference (simplified version of the IDTxl full algorithm) in demos/octave/EffectiveNetworkInference;
_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;
@ -150,11 +122,3 @@ _19/11/2012_ - New jar and full distribution files available for v0.1.2, includi
_31/10/2012_ - Jar and full distribution files available for v0.1.1 (first distribution)
_7/5/2012_ - JIDT project created and code uploaded
# Acknowledgements
This project has been supported by funding through:
* Australian Research Council Discovery Early Career Researcher Award (DECRA) "Relating function of complex networks to structure using information theory", J.T. Lizier, 2016-19 DE160100630
* Universities Australia - Deutscher Akademischer Austauschdienst (German Academic Exchange Service) UA-DAAD Australia-Germany Joint Research Co-operation grant "Measuring neural information synthesis and its impairment", Wibral, Lizier, Priesemann, Wollstadt, Finn, 2016-17
* University of Sydney Research Accelerator (SOAR) Fellowship 2019 Scheme, J.T. Lizier (CI), 2019-2020
* Australian Research Council Discovery Project "Large-scale computational modelling of epidemics in Australia: analysis, prediction and mitigation", M. Prokopenko, P. Pattison, M. Gambhir, J.T. Lizier, M. Piraveenan, 2016-19 DP160102742

49
build.xml
View File

@ -5,13 +5,12 @@
</description>
<!-- set global properties for this build -->
<property name="version" value="1.6.1"/>
<property name="version" value="1.4"/>
<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"/>
@ -43,17 +42,17 @@
<!-- Compile the java toolkit -->
<target name="compile" depends="init" description="compile the source">
<!-- Compile to Java 8 to provide compatibility for users with older JREs.
<!-- Compile to Java 6 to provide compatibility for users with older JREs.
Caveat: The flags here only check the language compatibility, but
may still use newer libraries which may cause issues for users with JDK 8.
Indeed, one gets the warning: "bootstrap class path not set in conjunction with -source 1.8"
may still use newer libraries which may cause issues for users with JDK 6.
Indeed, one gets the warning: "bootstrap class path not set in conjunction with -source 1.6"
To fix this, one would use the bootstrap classpath to point our JDK to an rt.jar
for Java 8.
At this stage, I'm sure I'm not using new library calls from Java 9+, so we can
ignore the warning, and I don't want to bother installing newer Java just to compile
like this. I'll endeavour not to use JDK 9+ libraries so as not to cause
for Java 6.
At this stage, I'm sure I'm not using new library calls from Java 7/8, so we can
ignore the warning, and I don't want to bother installing Java 7 just to compile
like this. I'll endeavour not to use JDK 8 libraries so as not to cause
any issues here ... -->
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.8" source="1.8" encoding="UTF8">
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.6" source="1.6">
<classpath refid="apache-classpath"/>
</javac>
@ -116,12 +115,12 @@
<target name="javadocs" depends="compile" description="Make the javadocs for the toolkit">
<delete dir="${javadocsdir}"/>
<javadoc destdir="${javadocsdir}">
<packageset dir="${src}">
<include name="**"/>
<exclude name="infodynamics/demos/**"/>
<exclude name="**/commonsmath3/*"/>
<exclude name="**/commonsmath3/**"/>
</packageset>
<fileset dir="${src}">
<include name="**/*.java"/>
<exclude name="infodynamics/demos/**/*.java"/>
<exclude name="**/commonsmath3/*.java"/>
<exclude name="**/commonsmath3/**/*.java"/>
</fileset>
</javadoc>
<!-- Change some of the style in the javadocs css for our lists: -->
<concat destfile="${javadocsdir}/stylesheet.css" append="true">
@ -149,7 +148,6 @@
<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.*"/>
@ -177,7 +175,7 @@
<!-- Compile and jar the toolkit with debug symbols -->
<target name="debug" depends="init" description="compile and jar with debug symbols">
<echo message="Compiling for debug"/>
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.7" source="1.7" debug="true">
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.6" source="1.6" debug="true">
<classpath refid="apache-classpath"/>
</javac>
@ -295,30 +293,23 @@
</zip>
</target>
<!-- Developer build - builds everything and makes the full distribution file in zip and tar.gz -->
<!-- Developer build - builds everything and makes the full distribution file -->
<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"/>
<zipfileset file="${jarplainname}" filemode="755"/>
<fileset file="${jarplainname}"/>
<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,**/*.sh,**/*.bat" prefix="demos"/>
<zipfileset dir="demos" includes="**/*.sh,**/*.bat" prefix="demos" filemode="755"/> <!-- Do these separately to get executable permissions -->
<zipfileset dir="demos" includes="**/*.*,**/*" excludes="clojure/deploy,clojure/deploy/*.*,python/*.pyc" prefix="demos"/>
<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="tutorial" prefix="tutorial"/>
<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
course/Module0-Overview/Lec0a-IntroToCourse.pdf
View File

Binary file not shown.

BIN
course/Module1-IntroductionToInformationTheory/Lec1-IntroductionToInformationTheory.pdf
View File

Binary file not shown.

BIN
course/Module10-EffectiveNetworkInference/Lec10-EffectiveNetworkInference.pdf
View File

Binary file not shown.

BIN
course/Module11-WrapUp/Lec11-InfoTheoryWrapUp.pdf
View File

Binary file not shown.

BIN
course/Module2-WhatIsInformation/Lec2-WhatIsInformation.pdf
View File

Binary file not shown.

BIN
course/Module3-JIDT/Lec3-IntroToJIDT.pdf
View File

Binary file not shown.

BIN
course/Module4-EstimatorsAndJIDT/Lec4-EstimatorsAndJIDT.pdf
View File

Binary file not shown.

BIN
course/Module5-StatisticalSignificance/Lec5-StatSignificanceAndUndersampling.pdf
View File

Binary file not shown.

BIN
course/Module6-InfoTheoryAndSelfOrganisation/Lec6-SelfOrganisation.pdf
View File

Binary file not shown.

BIN
course/Module7-InfoProcessingInComplexSystems/Lec7-InfoProcessing.pdf
View File

Binary file not shown.

BIN
course/Module8-InformationStorage/Lec8-InfoStorage.pdf
View File

Binary file not shown.

BIN
course/Module9-InformationTransfer/Lec9-InfoTransfer.pdf
View File

Binary file not shown.

BIN
course/README-Course.pdf
View File

Binary file not shown.

3
course/README.md
View File

@ -1,3 +0,0 @@
# Course
Please see the [Course page](../../../wiki/Course) on our wiki for full information about the course, including lecture slides and videos.

6
cuda/Makefile
View File

@ -74,7 +74,7 @@ endif
# Common includes and paths for CUDA. This assumes the CUDA toolkit is in PATH
INCLUDES := -I. -I./cub -I${JNI_INCLUDE} -I${JNI_PLATFORM_INCLUDE}
NVCCLDFLAGS := -L. -L$(BIN) -lcuda -lcudart
LDFLAGS += -L. -L$(BIN) -lcuda -lcudart
.PHONY: all clean test
@ -93,7 +93,7 @@ $(BIN)/libgpuKnnLibrary.a: $(BIN)/gpuKnnLibrary.o
## Compile host code
#-------------------
c_objects = $(addprefix $(BIN)/,digamma.o gpuMILibrary.o gpuCMILibrary.o kraskovCuda.o)
c_objects = $(addprefix $(BIN)/,digamma.o gpuMILibrary.o kraskovCuda.o)
$(BIN)/%.o: %.c
${GCC} ${INCLUDES} ${CCFLAGS} -x c -std=c99 -fPIC -c $< -o $@
@ -102,7 +102,7 @@ $(BIN)/%.o: %.c
## Final shared library linking
#------------------------------
$(BIN)/libKraskov.so: $(BIN)/libgpuKnnLibrary.a $(c_objects)
${NVCC} ${NVCCFLAGS} ${INCLUDES} -Xcompiler -fPIC -shared -o $@ $(c_objects) ${NVCCLDFLAGS} -lgpuKnnLibrary
${NVCC} ${NVCCFLAGS} ${INCLUDES} -Xcompiler -fPIC -shared -o $@ $(c_objects) ${LDFLAGS} -lgpuKnnLibrary
## Test binary targets

BIN
cuda/README-GPU.pdf
View File

Binary file not shown.

199
cuda/gpuCMILibrary.c
View File

@ -1,199 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "gpuCMILibrary.h"
#include "gpuKnnLibrary.h"
#include "digamma.h"
#include "ctimer.h"
jidt_error_t CMIKraskov_C(int N, float *source, int dimx, float *dest, int dimy,
float *cond, int dimz, int k, int thelier, int nb_surrogates,
int returnLocals, int useMaxNorm, int isAlgorithm1, float *result,
int variableToReorder) {
return CMIKraskovWithReorderings(N, source, dimx, dest, dimy, cond, dimz,
k, thelier, nb_surrogates, returnLocals, useMaxNorm, isAlgorithm1, result,
0, NULL, variableToReorder);
}
/**
* Calculate Mutual Information using the KSG algorithm.
*/
jidt_error_t CMIKraskovWithReorderings(int N, float *source, int dimx,
float *dest, int dimy, float *cond, int dimz, int k, int thelier,
int nb_surrogates, int returnLocals, int useMaxNorm,
int isAlgorithm1, float *result, int reorderingsGiven, int **reorderings,
int variableToReorder) {
CPerfTimer pt = startTimer("Rearranging pointset");
// Allocate more space if surrogates are requested
int nchunks = nb_surrogates + 1;
int dims = dimx + dimy + dimz;
float *pointset = (float *) malloc(N * dims * nchunks * sizeof(float));
if (nb_surrogates == 0) {
memcpy( pointset, source, N*dimx*sizeof(float));
memcpy( pointset + N*dimx, cond, N*dimz*sizeof(float));
memcpy(pointset + N*(dimx+dimz), dest, N*dimy*sizeof(float));
}
if (nb_surrogates > 0) {
for (int i = 0; i < N; i++) {
for (int j = 0; j < dimx; j++) {
pointset[j*N*nchunks+i] = source[N*j+i];
}
for (int j = 0; j < dimz; j++) {
pointset[nchunks*N*dimx + j*N*nchunks + i] = cond[N*j+i];
}
for (int j = 0; j < dimy; j++) {
pointset[nchunks*N*(dimx+dimz) + j*N*nchunks + i] = dest[N*j+i];
}
}
// If surrogates requested, copy permutations as well
int *order;
int perm[N];
if (!reorderingsGiven) {
for (int i = 0; i < N; i++) {
perm[i] = i;
}
}
for (int s = 0; s < nb_surrogates; s++) {
if (reorderingsGiven) {
order = reorderings[s];
} else {
randperm(perm, N);
order = perm;
}
for (int i = 0; i < N; i++) {
if (variableToReorder == 1) {
for (int j = 0; j < dimx; j++) {
pointset[(s+1)*N + N*j*nchunks + i] = source[N*j + order[i]];
}
for (int j = 0; j < dimz; j++) {
pointset[nchunks*N*dimx + (s+1)*N + N*j*nchunks + i] = cond[N*j + i];
}
} else {
for (int j = 0; j < dimx; j++) {
pointset[(s+1)*N + N*j*nchunks + i] = source[N*j + i];
}
for (int j = 0; j < dimz; j++) {
pointset[nchunks*N*dimx + (s+1)*N + N*j*nchunks + i] = cond[N*j + order[i]];
}
}
for (int j = 0; j < dimy; j++) {
pointset[nchunks*N*(dimx+dimz) + (s+1)*N + N*j*nchunks + i] = dest[N*j + i];
}
}
}
}
stopTimer(pt);
jidt_error_t err = CMIKraskovByPointsetChunks(N*nchunks, source, dimx,
dest, dimy, cond, dimz, k, thelier,
nchunks, returnLocals, useMaxNorm,
isAlgorithm1, result, pointset);
FREE(pointset);
return err;
}
jidt_error_t CMIKraskovByPointsetChunks(int signalLength, float *source, int dimx,
float *dest, int dimy, float *cond, int dimz, int k, int thelier, int nchunks,
int returnLocals, int useMaxNorm, int isAlgorithm1, float *result,
float *pointset) {
int dims = dimx + dimy + dimz;
int err;
int trialLength = signalLength/((float) nchunks);
float *d_source, *d_dest, *d_cond, *d_distances, *d_radii, *d_digammas;
int *d_nx, *d_ny, *d_nz, *d_indexes;
{
CPerfTimer pt = startTimer("GPU_warmup");
gpuWarmUp();
stopTimer(pt);
}
// 1. Allocate space in GPU and transfer memory
// ======================
allocateDeviceMemoryCMI(signalLength, k, dimx, dimy, dimz, &d_source, &d_dest, &d_cond,
&d_distances, &d_indexes, &d_radii, &d_nx, &d_ny, &d_nz, &d_digammas, pointset);
// 2. Find nearest neighbours in joint space
// ======================
{
CPerfTimer pt = startTimer("kNN_full");
d_cudaFindKnn(d_indexes, d_distances, d_source, d_source, k,
thelier, nchunks, dims, signalLength, useMaxNorm);
stopTimer(pt);
}
// 4. Count points strictly within R in the XZ-, YZ- and Z-spaces
// ======================
{
CPerfTimer pt = startTimer("RS_full");
d_cudaFindRSAll(d_nx, d_source, d_source, d_radii, thelier, nchunks, dimx + dimz, signalLength, useMaxNorm);
d_cudaFindRSAll(d_ny, d_cond, d_cond, d_radii, thelier, nchunks, dimy + dimz, signalLength, useMaxNorm);
d_cudaFindRSAll(d_nz, d_cond, d_cond, d_radii, thelier, nchunks, dimz, signalLength, useMaxNorm);
stopTimer(pt);
}
// 6. Set locals, surrogates or digammas for return
// ======================
{
CPerfTimer pt = startTimer("Digammas_full");
if (returnLocals) {
float digammaK = cpuDigamma(k);
float digammas[trialLength];
d_parallelDigammasCMI(digammas, d_digammas, d_nx, d_ny, d_nz, signalLength);
for (int i = 0; i < trialLength; i++) {
result[i] = digammaK - digammas[i];
}
} else {
float digammaK = cpuDigamma(k);
float sumDigammas[nchunks];
d_cudaSumDigammasCMI(sumDigammas, d_nx, d_ny, d_nz, d_digammas, trialLength, nchunks);
if (nchunks > 1) {
for (int ii = 0; ii < nchunks; ii++) {
result[ii] = digammaK - sumDigammas[ii]/((float) trialLength);
}
} else {
// Sign changed to comply with the returnValues processing in the Java
// KSG CMI calc, which is different from the one in the MI calc.
result[0] = -1 * sumDigammas[0];
result[1] = -1;
result[2] = -1;
result[3] = -1;
result[4] = -1;
result[5] = -1;
}
}
stopTimer(pt);
}
err = JIDT_SUCCESS;
freeDeviceMemory(d_source);
return err;
}

30
cuda/gpuCMILibrary.h
View File

@ -1,30 +0,0 @@
#ifndef GPUCMILIBRARY_H
#define GPUCMILIBRARY_H
#include "gpuKnnLibrary.h"
#define FREE(x) { if (x) free(x); x = NULL; }
#ifdef __cplusplus
extern "C" {
#endif
jidt_error_t CMIKraskovWithReorderings(int N, float *source, int dimx,
float *dest, int dimy, float *cond, int dimz,
int k, int thelier, int nb_surrogates, int returnLocals, int useMaxNorm,
int isAlgorithm1, float *result, int reorderingsGiven, int **reorderings,
int variableToReorder);
jidt_error_t CMIKraskov_C(int N, float *source, int dimx, float *dest, int dimy,
float *cond, int dimz, int k, int thelier, int nb_surrogates,
int returnLocals, int useMaxNorm, int isAlgorithm1, float *result,
int variableToReorder);
jidt_error_t CMIKraskovByPointsetChunks(int N, float *source, int dimx,
float *dest, int dimy, float *cond, int dimz, int k, int thelier, int nb_surrogates,
int returnLocals, int useMaxNorm, int isAlgorithm1, float *result, float *pointset);
#ifdef __cplusplus
}
#endif
#endif

20
cuda/gpuKnnBF_kernel.cu
View File

@ -469,26 +469,6 @@ __global__ void gpuDigammas(float *g_digammas, int *g_nx, int *g_ny, int signall
}
__global__ void gpuDigammasCMI(float *g_digammas, int *g_nx, int *g_ny, int *g_nz, int signallength) {
const unsigned int i = threadIdx.x + blockDim.x*blockIdx.x;
if(i < signallength){
// Fetch n and put it in thread memory
double dgX = (double) g_nx[i];
double dgY = (double) g_ny[i];
double dgZ = (double) g_nz[i];
// In-place digamma calculation
digammaXp1(&dgX);
digammaXp1(&dgY);
digammaXp1(&dgZ);
// Copy back to global memory
g_digammas[i] = (float) (dgX + dgY - dgZ);
}
return;
}

142
cuda/gpuKnnLibrary.c
View File

@ -57,57 +57,6 @@ int allocateDeviceMemory(int signalLength, int k, int dimx, int dimy,
return 1;
}
/**
* Allocate all necessary memory for the whole CMI calculation in a single call
* to cudaMalloc, and point the pointers to the right place.
*
* @param signalLength total number of samples given, including surrogates
* @param k nunmber of neighbours to find
* @param dimx dimension of source points
* @param dimy dimension of dest points
* @param dimz dimension of cond points
* @param source,dest,cond,distances,indexes,radii,nx,ny,nz,digammas device pointers
* @param pointset pointer to the data array in host memory
*
* @return error code
*/
int allocateDeviceMemoryCMI(int signalLength, int k, int dimx, int dimy, int dimz,
float **source, float **dest, float **cond, float **distances, int **indexes,
float **radii, int **nx, int **ny, int **nz, float **digammas, float *pointset) {
float *d_pointset;
int dims = dimx + dimy + dimz;
size_t mempointset = signalLength * dims * sizeof(float);
size_t memdistances = signalLength * k * sizeof(float);
size_t memindexes = signalLength * k * sizeof(int);
size_t memcounts = 3 * signalLength * sizeof(int);
size_t memdigammas = signalLength * sizeof(float);
size_t memtotal = mempointset + memdistances + memindexes + memcounts + memdigammas;
checkCudaErrors( cudaMalloc((void **) &d_pointset, memtotal) );
cudaError_t error = cudaGetLastError();
if(error!=cudaSuccess){
fprintf(stderr,"%s",cudaGetErrorString(error));
return 0;
}
checkCudaErrors( cudaMemcpy(d_pointset, pointset, mempointset, cudaMemcpyHostToDevice) );
*source = d_pointset;
*cond = *source + signalLength*dimx;
*dest = *cond + signalLength*dimz;
*distances = *dest + signalLength*dimy;
*radii = *distances + (k-1)*signalLength;
*indexes = (int *) (*distances + k*signalLength);
*nx = *indexes + signalLength;
*ny = *nx + signalLength;
*nz = *ny + signalLength;
*digammas = (float *) (*nz + signalLength);
return 1;
}
/**
* Free all the memory used in GPU (if allocated using allocateDeviceMemory.
*
@ -531,26 +480,6 @@ int d_parallelDigammas(float *digammas, float *d_digammas, int *d_nx,
}
int d_parallelDigammasCMI(float *digammas, float *d_digammas, int *d_nx,
int *d_ny, int *d_nz, int signalLength) {
// Kernel parameters
dim3 threads(1,1,1);
dim3 grid(1,1,1);
threads.x = 512;
grid.x = (signalLength-1)/threads.x + 1;
// Launch kernel
gpuDigammasCMI<<<grid.x, threads.x>>>(d_digammas, d_nx, d_ny, d_nz, signalLength);
checkCudaErrors( cudaDeviceSynchronize() );
checkCudaErrors( cudaMemcpy(digammas, d_digammas, signalLength * sizeof(float), cudaMemcpyDeviceToHost) );
checkCudaErrors( cudaDeviceSynchronize() );
return 1;
}
int parallelDigammas(float *digammas, int *nx, int *ny, int signalLength) {
int *d_nx, *d_ny;
@ -584,12 +513,24 @@ int parallelDigammas(float *digammas, int *nx, int *ny, int signalLength) {
return 1;
}
int cudaBlockReduce(float *sumDigammas, float *d_digammas, int trialLength, int nchunks) {
int d_cudaSumDigammas(float *sumDigammas, int *d_nx, int *d_ny,
float *d_digammas, int trialLength, int nchunks) {
float *d_sumDigammas;
int signalLength = trialLength * nchunks;
// Kernel parameters
dim3 threads(1,1,1);
dim3 grid(1,1,1);
threads.x = 512;
grid.x = (signalLength-1)/threads.x + 1;
checkCudaErrors( cudaMalloc((void **) &d_sumDigammas, nchunks * sizeof(int)) );
// Launch kernel to calculate (digamma(nx+1) + digamma(ny+1)), and leave
// results in GPU
gpuDigammas<<<grid.x, threads.x>>>(d_digammas, d_nx, d_ny, signalLength);
checkCudaErrors( cudaDeviceSynchronize() );
int offset_size = nchunks + 1;
int offsets[offset_size];
for (int i = 0; i < (nchunks+1); i++) { offsets[i] = i*trialLength; }
@ -617,46 +558,6 @@ int cudaBlockReduce(float *sumDigammas, float *d_digammas, int trialLength, int
checkCudaErrors( cudaFree(d_sumDigammas) );
return 1;
}
int d_cudaSumDigammas(float *sumDigammas, int *d_nx, int *d_ny,
float *d_digammas, int trialLength, int nchunks) {
int signalLength = trialLength * nchunks;
// Kernel parameters
dim3 threads(1,1,1);
dim3 grid(1,1,1);
threads.x = 512;
grid.x = (signalLength-1)/threads.x + 1;
// Launch kernel to calculate (digamma(nx+1) + digamma(ny+1)), and leave
// results in GPU
gpuDigammas<<<grid.x, threads.x>>>(d_digammas, d_nx, d_ny, signalLength);
checkCudaErrors( cudaDeviceSynchronize() );
return cudaBlockReduce(sumDigammas, d_digammas, trialLength, nchunks);
}
int d_cudaSumDigammasCMI(float *sumDigammas, int *d_nx, int *d_ny, int *d_nz,
float *d_digammas, int trialLength, int nchunks) {
int signalLength = trialLength * nchunks;
// Kernel parameters
dim3 threads(1,1,1);
dim3 grid(1,1,1);
threads.x = 512;
grid.x = (signalLength-1)/threads.x + 1;
// Launch kernel to calculate (digamma(nx+1) + digamma(ny+1)), and leave
// results in GPU
gpuDigammasCMI<<<grid.x, threads.x>>>(d_digammas, d_nx, d_ny, d_nz, signalLength);
checkCudaErrors( cudaDeviceSynchronize() );
return cudaBlockReduce(sumDigammas, d_digammas, trialLength, nchunks);
}
/**
@ -693,21 +594,6 @@ void device_reset(void) {
void gpuWarmUp(void) {
cudaSetDevice(0);
}
/**
* Make random permutation of perm[].
*
* @param perm preallocated and prefilled integer array to be shuffled
* @param n number of elements in perm
*/
void randperm(int perm[], int n) {
// Random permutation the order
for (int i = 0; i < n; i++) {
int j, t;
j = rand() % (n-i) + i;
t = perm[j]; perm[j] = perm[i]; perm[i] = t; // Swap i and j
}
}
#ifdef __cplusplus
}
#endif

16
cuda/gpuKnnLibrary.h
View File

@ -4,16 +4,10 @@
#ifdef __cplusplus
extern "C" {
#endif
typedef enum { JIDT_SUCCESS, JIDT_ERROR } jidt_error_t;
int allocateDeviceMemory(int signalLength, int kth, int dimx, int dimy,
float **source, float **dest, float **distances, int **indexes,
float **radii, int **nx, int **ny, float **digammas, float *pointset);
int allocateDeviceMemoryCMI(int signalLength, int k, int dimx, int dimy, int dimz,
float **source, float **dest, float **cond, float **distances, int **indexes,
float **radii, int **nx, int **ny, int **nz, float **digammas, float *pointset);
int freeDeviceMemory(float *d_pointset);
int cudaFindKnn(int* h_bf_indexes, float* h_bf_distances, float* h_pointset,
@ -53,22 +47,12 @@ int d_cudaFindRSAll(int* d_bf_npointsrange, float* d_bf_pointset, float* d_bf_qu
int d_parallelDigammas(float *digammas, float *d_digammas, int *d_nx,
int *d_ny, int signalLength);
int d_parallelDigammasCMI(float *digammas, float *d_digammas, int *d_nx,
int *d_ny, int *d_nz, int signalLength);
int cudaBlockReduce(float *sumDigammas, float *d_digammas, int trialLength, int nchunks);
int d_cudaSumDigammas(float *sumDigammas, int *d_nx, int *d_ny,
float *d_digammas, int trialLength, int nchunks);
int d_cudaSumDigammasCMI(float *sumDigammas, int *d_nx, int *d_ny, int *d_nz,
float *d_digammas, int trialLength, int nchunks);
void device_reset(void);
void gpuWarmUp(void);
void randperm(int perm[], int n);
#ifdef __cplusplus
}
#endif

16
cuda/gpuMILibrary.c
View File

@ -8,6 +8,22 @@
#include "ctimer.h"
/**
* Make random permutation of perm[].
*
* @param perm preallocated and prefilled integer array to be shuffled
* @param n number of elements in perm
*/
void randperm(int perm[], int n) {
// Random permutation the order
for (int i = 0; i < n; i++) {
int j, t;
j = rand() % (n-i) + i;
t = perm[j]; perm[j] = perm[i]; perm[i] = t; // Swap i and j
}
}
jidt_error_t MIKraskov_C(int N, float *source, int dimx, float *dest, int dimy,
int k, int thelier, int nb_surrogates, int returnLocals, int useMaxNorm,
int isAlgorithm1, float *result) {

6
cuda/gpuMILibrary.h
View File

@ -1,13 +1,13 @@
#ifndef GPUMILIBRARY_H
#define GPUMILIBRARY_H
#include "gpuKnnLibrary.h"
#define FREE(x) { if (x) free(x); x = NULL; }
#ifdef __cplusplus
extern "C" {
#endif
typedef enum { JIDT_SUCCESS, JIDT_ERROR } jidt_error_t;
jidt_error_t MIKraskovWithReorderings(int N, float *source, int dimx, float *dest, int dimy,
int k, int thelier, int nb_surrogates, int returnLocals, int useMaxNorm,
int isAlgorithm1, float *result, int reorderingsGiven, int **reorderings);
@ -20,6 +20,8 @@ jidt_error_t MIKraskovByPointsetChunks(int N, float *source, int dimx,
float *dest, int dimy, int k, int thelier, int nb_surrogates,
int returnLocals, int useMaxNorm, int isAlgorithm1, float *result,
float *pointset);
void randperm(int perm[], int n);
#ifdef __cplusplus
}
#endif

188
cuda/kraskovCuda.c
View File

@ -5,7 +5,6 @@
#define check(ans) { _check((ans), __FILE__, __LINE__); }
#include "gpuMILibrary.h"
#include "gpuCMILibrary.h"
#include "ctimer.h"
#ifdef __cplusplus
@ -197,193 +196,6 @@ JNIEXPORT jdoubleArray JNICALL
} // End of function MIKraskov
/*
* Class: infodynamics_measures_continuous_kraskov_ConditionalMutualInfoCalculatorMultiVariateKraskov
* Method: CMIKraskov
* Signature: (I[DI[DI[DIIIZZZIZ[II)[D
*/
JNIEXPORT jdoubleArray JNICALL
Java_infodynamics_measures_continuous_kraskov_ConditionalMutualInfoCalculatorMultiVariateKraskov_CMIKraskov(
JNIEnv *env, jobject thisObj, jint j_N,
jobjectArray j_sourceArray, jint j_dimx,
jobjectArray j_destArray, jint j_dimy,
jobjectArray j_condArray, jint j_dimz,
jint j_k, jint j_theiler, jboolean j_returnLocals,
jboolean j_useMaxNorm, jboolean j_isAlgorithm1, jint j_nbSurrogates,
jboolean j_reorderingsGiven, jobjectArray j_orderings,
jint j_variableToReorder) {
// Check that incoming data has correct size
// =====================
jsize sourceLength = (*env)->GetArrayLength(env, j_sourceArray);
jsize destLength = (*env)->GetArrayLength(env, j_destArray);
jsize condLength = (*env)->GetArrayLength(env, j_condArray);
// if (sourceLength != j_N || destLength != j_N || (j_N%(j_nbSurrogates+1) != 0)) {
if (sourceLength != j_N || destLength != j_N || condLength != j_N) {
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
(*env)->ThrowNew(env, Exception, "Data has wrong length.");
}
if (!j_isAlgorithm1) {
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
(*env)->ThrowNew(env, Exception, "Only algorithm 1 is supported.");
}
if ((j_returnLocals || !j_isAlgorithm1) && (j_nbSurrogates > 0)) {
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
(*env)->ThrowNew(env, Exception, "Surrogates only supported for average MI with KSG1.");
}
// Copy variables from Java
// =====================
int N = j_N;
int k = j_k;
int dimx = j_dimx;
int dimy = j_dimy;
int dimz = j_dimz;
int theiler = j_theiler;
int returnLocals = j_returnLocals ? 1 : 0;
int useMaxNorm = j_useMaxNorm ? 1 : 0;
int isAlgorithm1 = j_isAlgorithm1 ? 1 : 0;
int nb_surrogates = j_nbSurrogates;
int reorderingsGiven = j_reorderingsGiven ? 1 : 0;
int variableToReorder = j_variableToReorder;
CPerfTimer pt = startTimer("Java array copy");
float *source = (float *) malloc(N * dimx * sizeof(float));
float *dest = (float *) malloc(N * dimy * sizeof(float));
float *cond = (float *) malloc(N * dimz * sizeof(float));
if (NULL == source || NULL == dest || NULL == cond) {
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
(*env)->ThrowNew(env, Exception, "Error allocating data.");
}
for (int i = 0; i < N; i++) {
jdoubleArray j_sourceRow = (jdoubleArray) (*env)->GetObjectArrayElement(env, j_sourceArray, i);
jdoubleArray j_destRow = (jdoubleArray) (*env)->GetObjectArrayElement(env, j_destArray, i);
jdoubleArray j_condRow = (jdoubleArray) (*env)->GetObjectArrayElement(env, j_condArray, i);
jdouble *sourceRow = (*env)->GetDoubleArrayElements(env, j_sourceRow, NULL);
jdouble *destRow = (*env)->GetDoubleArrayElements(env, j_destRow, NULL);
jdouble *condRow = (*env)->GetDoubleArrayElements(env, j_condRow, NULL);
// Data in java are doubles, but GPUs need floats.
// We have to cast them manually, so we can't memcopy
// The following for-loops get two matrices in T-by-D indexing (i.e.
// first dimension is time, second is variable) and return the data in
// column-major form
for (int j = 0; j < dimx; j++) {
source[N*j + i] = (float) sourceRow[j];
}
for (int j = 0; j < dimy; j++) {
dest[N*j + i] = (float) destRow[j];
}
for (int j = 0; j < dimz; j++) {
cond[N*j + i] = (float) condRow[j];
}
(*env)->ReleaseDoubleArrayElements(env, j_sourceRow, sourceRow, 0);
(*env)->ReleaseDoubleArrayElements(env, j_destRow, destRow, 0);
(*env)->ReleaseDoubleArrayElements(env, j_condRow, condRow, 0);
(*env)->DeleteLocalRef(env, j_sourceRow);
(*env)->DeleteLocalRef(env, j_destRow);
(*env)->DeleteLocalRef(env, j_condRow);
// FIXME: I'm not entirely sure I'm freeing all the memory here. I should
// check for memory leaks more carefully.
}
int **reorderings = NULL;
if (reorderingsGiven) {
reorderings = (int **) malloc(nb_surrogates * sizeof(int *));
for (int i = 0; i < nb_surrogates; i++) {
jintArray j_order = (jdoubleArray) (*env)->GetObjectArrayElement(env, j_orderings, i);
jint *order = (*env)->GetIntArrayElements(env, j_order, NULL);
reorderings[i] = (int *) malloc(N * sizeof(int));
for (int j = 0; j < N; j++) {
reorderings[i][j] = order[j];
}
(*env)->ReleaseIntArrayElements(env, j_order, order, 0);
(*env)->DeleteLocalRef(env, j_order);
}
}
stopTimer(pt);
// Call C function
// =========================
int resultSize;
if (returnLocals) {
resultSize = N;
} else if (nb_surrogates > 0) {
resultSize = nb_surrogates + 1;
} else {
resultSize = 6;
}
float *result = (float *) malloc(resultSize * sizeof(float));
jidt_error_t ret;
if (!reorderingsGiven) {
ret = CMIKraskov_C(N, source, dimx, dest, dimy, cond, dimz,
k, theiler, nb_surrogates, returnLocals, useMaxNorm,
isAlgorithm1, result, variableToReorder);
} else {
ret = CMIKraskovWithReorderings(N, source, dimx, dest, dimy, cond, dimz, k, theiler,
nb_surrogates, returnLocals, useMaxNorm,
isAlgorithm1, result, reorderingsGiven,
reorderings, variableToReorder);
}
if (JIDT_ERROR == ret) {
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
(*env)->ThrowNew(env, Exception, "Error in GPU execution.");
}
// Free memory and return
// =========================
if (source) free(source);
if (dest) free(dest);
if (cond) free(cond);
if (reorderingsGiven) {
for (int i = 0; i < nb_surrogates; i++) {
if (reorderings[i]) free(reorderings[i]);
}
if (reorderings) free (reorderings);
}
jdouble outCArray[resultSize];
for (int i = 0; i < resultSize; i++) {
outCArray[i] = result[i];
}
// Set Java array for return
jdoubleArray outJNIArray = (*env)->NewDoubleArray(env, resultSize); // allocate
if (NULL == outJNIArray) {
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");
(*env)->ThrowNew(env, Exception, "Error creating return array.");
}
(*env)->SetDoubleArrayRegion(env, outJNIArray, 0 , resultSize, outCArray); // copy
if (result) { free(result); }
return outJNIArray;
} // End of function CMIKraskov
#ifdef __cplusplus
}
#endif

246
cuda/unittest.cpp
View File

@ -5,7 +5,6 @@
#include "gpuKnnLibrary.h"
#include "gpuMILibrary.h"
#include "gpuCMILibrary.h"
#include "digamma.h"
using lest::approx;
@ -607,46 +606,6 @@ CASE("Test correct pointset arrangement without reorderings")
EXPECT(result1[2] == approx(result2[2]));
},
CASE("Test correct pointset arrangement without reorderings in CMI")
{
int N = 10;
int dimx = 1;
int dimy = 1;
int dimz = 1;
int k = 2;
int thelier = 0;
int returnLocals = 0;
int useMaxNorm = 1;
int isAlgorithm1 = 1;
float source[10] = {0.4, 1, -4, 1, 1, 0.2, 98, 12, 1.2, 1.3};
float dest[10] = { -3, 1, 3, -2, 2.1, 8.5, 4.2, 100, 12, 0};
float cond[10] = { -1, 4, 3, -8, 0.3, 2.1, 3.2, 111, 32, 7};
float pointset[30] = {0.4, 1, -4, 1, 1, 0.2, 98, 12, 1.2, 1.3,
-1, 4, 3, -8, 0.3, 2.1, 3.2, 111, 32, 7,
-3, 1, 3, -2, 2.1, 8.5, 4.2, 100, 12, 0};
float result1[6];
float result2[6];
jidt_error_t err;
err = CMIKraskov_C(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
0, returnLocals, useMaxNorm, isAlgorithm1, result1, 1);
EXPECT(err == JIDT_SUCCESS);
err = CMIKraskovByPointsetChunks(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result2, pointset);
EXPECT(err == JIDT_SUCCESS);
EXPECT(result1[0] == approx(result2[0]));
EXPECT(result1[1] == approx(result2[1]));
EXPECT(result1[2] == approx(result2[2]));
EXPECT(result1[3] == approx(result2[3]));
EXPECT(result1[4] == approx(result2[4]));
EXPECT(result1[5] == approx(result2[5]));
},
CASE("Test correct pointset arrangement in more than one dimension")
{
int N = 5;
@ -695,67 +654,6 @@ CASE("Test correct pointset arrangement in more than one dimension")
EXPECT(result1[2] == approx(result2[2]));
},
CASE("Test correct pointset arrangement in more than one dimension for CMI")
{
int N = 5;
int dimx = 2;
int dimy = 2;
int dimz = 2;
int k = 2;
int thelier = 0;
int returnLocals = 0;
int useMaxNorm = 1;
int isAlgorithm1 = 1;
// Source points: X Y
// 0.4 0.2
// 1 98
// -4 12
// 1 1.2
// 1 1.3
//
// Dest points: X Y
// -3 8.5
// 1 4.2
// 3 100
// -2 12
// 2.1 0
//
// Cond points: X Y
// -1 2.1
// 4 3.2
// 3 111
// -8 32
// 0.3 7
float source[10] = {0.4, 1, -4, 1, 1, 0.2, 98, 12, 1.2, 1.3};
float dest[10] = { -3, 1, 3, -2, 2.1, 8.5, 4.2, 100, 12, 0};
float cond[10] = { -1, 4, 3, -8, 0.3, 2.1, 3.2, 111, 32, 7};
float pointset[30] = {0.4, 1, -4, 1, 1, 0.2, 98, 12, 1.2, 1.3,
-1, 4, 3, -8, 0.3, 2.1, 3.2, 111, 32, 7,
-3, 1, 3, -2, 2.1, 8.5, 4.2, 100, 12, 0};
float result1[6];
float result2[6];
jidt_error_t err;
err = CMIKraskov_C(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
0, returnLocals, useMaxNorm, isAlgorithm1, result1, 1);
EXPECT(err == JIDT_SUCCESS);
err = CMIKraskovByPointsetChunks(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result2, pointset);
EXPECT(err == JIDT_SUCCESS);
EXPECT(result1[0] == approx(result2[0]));
EXPECT(result1[1] == approx(result2[1]));
EXPECT(result1[2] == approx(result2[2]));
EXPECT(result1[3] == approx(result2[3]));
EXPECT(result1[4] == approx(result2[4]));
EXPECT(result1[5] == approx(result2[5]));
},
CASE("Test that same sample in repeated chunks gives same result")
{
int N = 5;
@ -793,46 +691,6 @@ CASE("Test that same sample in repeated chunks gives same result")
},
CASE("Test that same sample in repeated chunks gives same result in CMI")
{
int N = 5;
int dimx = 1;
int dimy = 1;
int dimz = 1;
// Sample source and dest data
float source[5] = {0.4, 1, -4, 1, 1};
float dest[5] = { -3, 1, 3, -2, 2.1};
float cond[5] = { -1, 4, 3, -8, 0.3};
// Pointset with source and dest repeated twice
float double_pointset[30] = {0.4, 1, -4, 1, 1, 0.4, 1, -4, 1, 1,
-1, 4, 3, -8, 0.3, -1, 4, 3, -8, 0.3,
-3, 1, 3, -2, 2.1, -3, 1, 3, -2, 2.1};
int k = 2;
int thelier = 0;
int returnLocals = 0;
int useMaxNorm = 1;
int isAlgorithm1 = 1;
float result1[3];
float result2[2];
jidt_error_t err;
err = CMIKraskov_C(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
0, returnLocals, useMaxNorm, isAlgorithm1, result1, 1);
err = CMIKraskovByPointsetChunks(N*2, source, dimx, dest, dimy, cond, dimz, k, thelier,
2, returnLocals, useMaxNorm, isAlgorithm1, result2, double_pointset);
float CMI1 = cpuDigamma(k) + result1[0]/((double) N);
EXPECT(err == JIDT_SUCCESS);
EXPECT(result2[0] == approx(CMI1));
EXPECT(result2[0] == approx(result2[1]));
},
CASE("Test that same sample of 2D data in repeated chunks gives same result")
{
int N = 5;
@ -905,39 +763,13 @@ CASE("Test that sample and identity reordering have same MI")
int **reorderings = &order_p;
jidt_error_t err;
printf("====================================\n");
err = MIKraskovWithReorderings(N, source, dimx, dest, dimy, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result, reorderingsGiven, reorderings);
EXPECT(err == JIDT_SUCCESS);
EXPECT(result[0] == approx(result[1]));
},
CASE("Test that sample and identity reordering have same CMI")
{
int N = 5;
int dimx = 1;
int dimy = 1;
int dimz = 1;
float source[5] = {0.4, 1, -4, 1, 1};
float dest[5] = { -3, 1, 3, -2, 2.1};
float cond[5] = { -1, 4, 3, -8, 0.3};
int k = 2;
int thelier = 0;
int returnLocals = 0;
int useMaxNorm = 1;
int isAlgorithm1 = 1;
float result[2];
int reorderingsGiven = 1;
int order[5] = {0, 1, 2, 3, 4};
int *order_p = order;
int **reorderings = &order_p;
jidt_error_t err;
err = CMIKraskovWithReorderings(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result, reorderingsGiven, reorderings, 1);
EXPECT(err == JIDT_SUCCESS);
printf("Test results: %f\t%f\n", result[0], result[1]);
EXPECT(result[0] == approx(result[1]));
},
@ -976,6 +808,7 @@ CASE("Test identity reordering with more than one dimension")
int **reorderings = &order_p;
jidt_error_t err;
printf("====================================\n");
err = MIKraskovWithReorderings(N, source, dimx, dest, dimy, k, thelier,
0, returnLocals, useMaxNorm, isAlgorithm1, result1, reorderingsGiven, reorderings);
@ -985,6 +818,8 @@ CASE("Test identity reordering with more than one dimension")
1, returnLocals, useMaxNorm, isAlgorithm1, result2, reorderingsGiven, reorderings);
EXPECT(err == JIDT_SUCCESS);
printf("Test results: %f\t%f\n", result2[0], result2[1]);
EXPECT(result2[0] == approx(result2[1]));
},
@ -1014,12 +849,15 @@ CASE("Test identity reordering with larger dataset")
float result[2];
jidt_error_t err;
printf("====================================\n");
err = MIKraskovWithReorderings(N, source, dimx, dest, dimy, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result, reorderingsGiven, reorderings);
free(source); free(dest); free(order);
EXPECT(err == JIDT_SUCCESS);
printf("Test results: %f\t%f\n", result[0], result[1]);
EXPECT(result[0] == approx(result[1]));
},
@ -1061,10 +899,13 @@ CASE("Test non-identity reordering in 2D")
float result[2];
jidt_error_t err;
printf("====================================\n");
err = MIKraskovWithReorderings(N, source, dimx, dest, dimy, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result, reorderingsGiven, reorderings);
EXPECT(err == JIDT_SUCCESS);
printf("Test results: %f\t%f\n", result[0], result[1]);
EXPECT(result[0] != approx(result[1]));
},
@ -1090,12 +931,15 @@ CASE("Test random surrogates in 2D")
float result[3];
jidt_error_t err;
printf("====================================\n");
err = MIKraskov_C(N, source, dimx, dest, dimy, k, thelier,
2, returnLocals, useMaxNorm, isAlgorithm1, result);
free(source); free(dest);
EXPECT(err == JIDT_SUCCESS);
printf("Test results: %f\t%f\n", result[0], result[1]);
EXPECT(result[0] != result[1]);
EXPECT(result[0] != result[2]);
EXPECT(result[1] != result[2]);
@ -1134,77 +978,23 @@ CASE("Test that the first result of calculation with surrogates is the same as w
float result2[2];
jidt_error_t err;
printf("====================================\n");
err = MIKraskov_C(N, source, dimx, dest, dimy, k, thelier,
0, returnLocals, useMaxNorm, isAlgorithm1, result1);
EXPECT(err == JIDT_SUCCESS);
float MI1 = cpuDigamma(k) + cpuDigamma(N) - result1[0]/((double) N);
printf("====================================\n");
err = MIKraskov_C(N, source, dimx, dest, dimy, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result2);
EXPECT(err == JIDT_SUCCESS);
printf("No surrogates: %f, with 1 surrogate: %f, %f\n", MI1, result2[0], result2[1]);
EXPECT(MI1 == approx(result2[0]));
EXPECT(result2[0] != result2[1]);
},
CASE("Test that the first result of calculation with surrogates is the same as without in CMI")
{
int thelier = 0;
int useMaxNorm = 1;
int N = 10;
int dimx = 1;
int dimy = 1;
int dimz = 1;
int k = 2;
int isAlgorithm1 = 1;
int returnLocals = 0;
// Source points: X Y Z
// 0.4 0.2 0
// 1 98 13
// -4 12 7
// 1 1.2 -1
// 1 1.3 0
//
// Dest points: X Y
// -3 8.5
// 1 4.2
// 3 100
// -2 13
// 2.1 0
//
// Cond points: X Y
// -1 2.1
// 4 3.2
// 3 111
// -8 32
// 0.3 7
float source[15] = {0.4, 1, -4, 1, 1, 0.2, 98, 12, 1.2, 1.3, 0, 13, 7, -1, 0};
float dest[10] = { -3, 1, 3, -2, 2.1, 8.5, 4.2, 100, 13, 0};
float cond[10] = { -1, 4, 3, -8, 0.3, 2.1, 3.2, 111, 32, 7};
float result1[6];
float result2[2];
jidt_error_t err;
err = CMIKraskov_C(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
0, returnLocals, useMaxNorm, isAlgorithm1, result1, 1);
EXPECT(err == JIDT_SUCCESS);
float CMI1 = cpuDigamma(k) + result1[0]/((double) N);
err = CMIKraskov_C(N, source, dimx, dest, dimy, cond, dimz, k, thelier,
1, returnLocals, useMaxNorm, isAlgorithm1, result2, 1);
EXPECT(err == JIDT_SUCCESS);
EXPECT(CMI1 == approx(result2[0]));
EXPECT(result2[0] != result2[1]);
},
CASE("Basic digamma sum test")
{
int N = 2;

BIN
demos/AutoAnalyser/README-AutoAnalyserDemo.pdf
View File

Binary file not shown.

8
demos/AutoAnalyser/launchAutoAnalyser.m
View File

@ -1,8 +0,0 @@
% Launch the AutoAnalyser from within Matlab. You will need to have Matlab open at the demos/AutoAnalyser folder.
% This is useful where you don't have a separate Java Runtime Environment installed, and so utilise Matlab's
warning('off','MATLAB:Java:DuplicateClass');
javaaddpath('../../infodynamics.jar');
autoAnalyser = javaObject('infodynamics.demos.autoanalysis.AutoAnalyserLauncher', false);
fprintf('Be warned - closing the AutoAnalyser applet seems to close Matlab completely, so only do so when ready to exit Matlab!\n');

11
demos/AutoAnalyser/launchAutoAnalyser.py
View File

@ -1,11 +0,0 @@
from jpype import *
# Add JIDT jar library to the path
jarLocation = "../../infodynamics.jar"
# 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)
# Construct an AutoAnalyser
calcClass = JPackage("infodynamics.demos.autoanalysis").AutoAnalyserLauncher
calc = calcClass(False)

5
demos/AutoAnalyser/runAISAutoAnalyser.bat Executable file
View File

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

5
demos/AutoAnalyser/runAISAutoAnalyser.sh Executable file
View File

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

5
demos/AutoAnalyser/runCMIAutoAnalyser.bat Executable file
View File

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

5
demos/AutoAnalyser/runCMIAutoAnalyser.sh Executable file
View File

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

5
demos/AutoAnalyser/runCTEAutoAnalyser.bat Executable file
View File

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

5
demos/AutoAnalyser/runCTEAutoAnalyser.sh Executable file
View File

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

5
demos/AutoAnalyser/runHAutoAnalyser.bat Executable file
View File

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

5
demos/AutoAnalyser/runHAutoAnalyser.sh Executable file
View File

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

5
demos/AutoAnalyser/runMIAutoAnalyser.bat Executable file
View File

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

5
demos/AutoAnalyser/runMIAutoAnalyser.sh Executable file
View File

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

5
demos/AutoAnalyser/runTEAutoAnalyser.bat Executable file
View File

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

5
demos/AutoAnalyser/runTEAutoAnalyser.sh Executable file
View File

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

BIN
demos/README-Demos.pdf
View File

Binary file not shown.

10000
demos/data/10ColsGaussianSelfDependence2Steps.txt
View File

File diff suppressed because it is too large Load Diff

1210
demos/data/SFI-heartRate_breathVol_bloodOx-extract.txt
View File

File diff suppressed because it is too large Load Diff

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

3
demos/java/infodynamics/demos/Example10GPUBenchmark.java
View File

@ -67,10 +67,9 @@ 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.computeSignificance(nb_surrogates).actualValue;
timeAndValue[1] = miCalc.computeAverageLocalOfObservations();
timeAndValue[0] = (System.nanoTime() - startTime)/1000000.0;
return timeAndValue;
}

7
demos/octave/CellularAutomata/plotLocalInfoMeasureForCA.m
View File

@ -45,7 +45,6 @@
% - measureParams.j - we measure information transfer across j cells to the right per time step
% - options - a stucture containing a range of other options, i.e.:
% - plotOptions - structure as defined for the plotLocalInfoValues function
% - plotOptions.figNum - figure number to plot the info value profile in
% - seed - state for the random number generator used to set the initial condition of the CA (use this
% for reproducibility of plots, or to produce profiles for several different measures of the same CA raw states).
% We set rand('state', options.seed) if options.seed is supplied, and restore the previous seed afterwards.
@ -108,11 +107,7 @@ function [caStates, localValues] = plotLocalInfoMeasureForCA(neighbourhood, base
printDriver = options.saveImagesFormat;
fontSize = 13;
end
if (not(isfield(options.plotOptions, 'figNum')))
figNum = 2;
else
figNum = options.plotOptions.figNum;
end
figNum = 2;
toc
% The offsets of the parents (see runCA for how this is computed, especially for even neighbourhood):
fullSetOfParents = ceil(-neighbourhood / 2) : ceil(-neighbourhood / 2) + (neighbourhood-1);

6
demos/octave/CellularAutomata/plotLocalInfoValues.m
View File

@ -176,11 +176,7 @@ function h = plotLocalInfoValues(localResults, plotOptions)
bluemap = prepareColourmap(mainSignVectorLength, true, scalingMainComponent, scalingScdryComponent, gammaPower);
colormap(bluemap);
% Now, plot the local values with the pre-prepared colormap
if (maxLocal == 0)
h = imagesc(localResultsToPlot, [0, 1]);
else
h = imagesc(localResultsToPlot, [0, maxLocal]);
end
h = imagesc(localResultsToPlot, [0, maxLocal]);
end
axis([0.5 (plotOptions.plotCols+0.5) 0.5 (plotOptions.plotRows+0.5)]);
colorbar

30
demos/octave/CellularAutomata/runCA.m
View File

@ -46,24 +46,17 @@
% - debug - turn on various debug messages
% - seedOrState (optional) - if a scalar, it is the state input for the random number generator (so one can repeat CA investigations for the same initial state).
% - if a vector, it is the initial state for the CA (must be of length cells)
% - alterProbability (optional, default 0) - defines stochastic probability of altering the value in the update of any given cell
%
% Outputs:
% - caStates - a run, from random initial conditions, of a CA of the given parameters.
% - ruleTable - the lookup table for each neighbourhood configuration, constructed from the rule number
% - executedRules - which CA rule was executed for every cell update that occurred for the CA.
function [caStates, ruleTable, executedRules] = runCA(neighbourhood, base, rule, cells, steps, debug, seedOrState, alterProbability)
function [caStates, ruleTable, executedRules] = runCA(neighbourhood, base, rule, cells, steps, debug, seedOrState)
% Check arguments:
ca = [];
if (nargin < 8) || isempty(alterProbability)
alterProbability = 0;
end
if (nargin < 7)
seedOrState = [];
end
if (~isempty(seedOrState))
if (nargin >= 7)
if (isscalar(seedOrState))
% User has supplied seed for the random number generator:
fprintf('Generating initial random CA state from seed %d\n', seedOrState);
@ -79,13 +72,13 @@ function [caStates, ruleTable, executedRules] = runCA(neighbourhood, base, rule,
else
fprintf('Generating initial random CA state\n');
end
if (nargin < 6) || isempty(debug)
if (nargin < 6)
debug = false;
end
if (nargin < 5) || isempty(steps)
if (nargin < 5)
steps = 100;
end
if (nargin < 4) || isempty(cells)
if (nargin < 4)
cells = 100;
end
if (nargin < 3)
@ -212,19 +205,6 @@ function [caStates, ruleTable, executedRules] = runCA(neighbourhood, base, rule,
% Need to add 1 to the ruleToRun because of the indexing starting from 1 not 0.
ca = ruleTable(ruleToRun + 1)';
% Check whether any of the CA states should be altered
if (alterProbability > 0)
% Could code this in a faster way, but won't matter too much
for c = 1 : cells
if (rand() < alterProbability)
% Alter this cell
chooseFrom = 0:(base-1);
chooseFrom(ca(c)+1) = []; % Remove the current value as an option so that we change this one
ca(c) = chooseFrom(randsample(length(chooseFrom), 1));
end
end
end
if (debug)
ca
end

357
demos/octave/EffectiveNetworkInference/greedyInferParents.m
View File

@ -1,357 +0,0 @@
%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2021, Joseph T. Lizier
% Distributed under GNU General Public License v3
%
% Infer the parent source variables to a given target, using the greedy/iterative/multivariate algorithm with TE.
% This a simplistic implementation of the full algorithm implemented in IDTxl - https://github.com/pwollstadt/idtxl -
% you are referred to IDTxl for an implementation with full features available.
%
% Inputs:
% - data: multivariate data, indexed by time,variableNumber, and possibly trialNumber.
% - parameters: an object containing the expected properties, or a string
% describing the filename to run load this object in. Can include:
% - parameters.calcType: which estimator type to use, from options
% 'discrete', 'gaussian', 'ksg', 'granger' --
% granger is equivalent to gaussian, but performed with Oliver Cliff's
% toolkit, to include proper autocorrelation correction.
% - parameters.timePointsToSkipAtStart: number of time points we'll skip
% at the start (default 0)
% - parameters.timePointsToSkipAtEnd: number of time points we'll skip
% at the end (default 0)
% - parameters.numDiscreteBins: alphabet size for the discrete variables when used. (default 2)
% - parameters.k: target embedding length, can be 'auto' to indicate auto-embedding (default 1)
% - parameters.k_max: max target embedding length to use when parameters.k == 'auto' (default 10)
% - parameters.numSurrogates: number of surrogates to run, or 0 for analytic surrogates (default 1000)
% - parameters.maxDynCorrExclLags: maximum length of dynamic correlation exclusion, which will be auto-fitted (default 50)
% - parameters.jidtLocation: path to the JIDT folder
% - parameters.gcToolkitLocation: Oliver's toolkit location for Granger
% - parameters.debug: whether to print debugging results as parents are inferred (default true)
% - targetIndex: which column number to run the inference for
% - uncorrectedThresholdForOneTarget: p-value threshold (where 0 is most significant) to select sources. We will Bonferroni correct this over sources here (but user should correct over targets if they wish)
function [parentSet, results, pValues, otherResults] = greedyInferParents(data, parameters, targetIndex, uncorrectedThresholdForOneTarget)
tic
if ischar(parameters)
% Assume that this string contains a filename which when run will load
% a properties object for this run
eval(['run ', parameters]);
end
% Add JIDT jar library to the path, and disable warnings that it's already there:
warning('off','MATLAB:Java:DuplicateClass');
%javaaddpath('/home/joseph/temp/jidt-master/jidt/infodynamics.jar');
javaaddpath([parameters.jidtLocation, 'infodynamics.jar']);
% Add utilities to the path
% addpath('/home/joseph/temp/jidt-master/jidt/demos/octave');
addpath([parameters.jidtLocation, 'demos/octave']);
% Add Oliver Cliff's toolkit to path
addpath(genpath(parameters.gcToolkitLocation));
% Set parameter defaults:
if ~isfield(parameters, 'numDiscreteBins')
parameters.numDiscreteBins = 2;
end
if ~isfield(parameters, 'k')
parameters.k = 1;
end
if ~isfield(parameters, 'k_max')
parameters.k_max = 10;
end
if ~isfield(parameters, 'numSurrogates')
parameters.numSurrogates = 1000;
end
if ~isfield(parameters, 'verbose')
parameters.verbose = true;
end
if ~isfield(parameters, 'timePointsToSkipAtStart')
parameters.timePointsToSkipAtStart = 0;
end
if ~isfield(parameters, 'timePointsToSkipAtEnd')
parameters.timePointsToSkipAtEnd = 0;
end
if ~isfield(parameters, 'maxDynCorrExclLags')
parameters.maxDynCorrExclLags = 50;
end
T = size(data,1); % TIMEPOINTS
N = size(data,2); % NODES
hasMultipleTrials = (length(size(data)) > 2);
if (hasMultipleTrials)
R = size(data,3); % TRIALS
% Analyse data(time, variables), but skip the first
% and last few steps
data = data(parameters.timePointsToSkipAtStart+1:end-parameters.timePointsToSkipAtEnd,:,:);
else
% Analyse data(time, variables), but skip the first
% and last few steps
data = data(parameters.timePointsToSkipAtStart+1:end-parameters.timePointsToSkipAtEnd,:);
end
threshold = uncorrectedThresholdForOneTarget / (N-1); % Bonferroni correcting the threshold
parentSet = [];
pValues = [];
results = [];
otherResults.k_history = 1;
candidates = [1:targetIndex-1, targetIndex+1:N]; % All sources except the targets are candidates as parents
% Set boolean flags for which calculator we are doing
is_jidt = true;
is_discrete = false;
is_ksg= false;
javaConverterSingleArray = 'octaveToJavaDoubleArray';
javaConverterMatrix = 'octaveToJavaDoubleMatrix';
if (strcmp(parameters.calcType, 'granger'))
is_jidt = false;
if (~ischar(parameters.k))
% if (ischar(parameters.k) && ~strcmp(parameters.k, 'auto')) % I don't think this logic was correct
parameters.k = char(string(parameters.k));
end
if (hasMultipleTrials)
error('Granger calculator does not currently support multiple trials');
end
elseif (strcmp(parameters.calcType, 'discrete'))
is_discrete = true;
javaConverterSingleArray = 'octaveToJavaIntArray';
javaConverterMatrix = 'octaveToJavaIntMatrix';
elseif (strcmp(parameters.calcType, 'ksg'))
is_ksg = true;
end
acfDecayTimes = -1 * ones(N,1);
% Grab the ACF decay time for the target
if (is_jidt)
if (~hasMultipleTrials)
acfDecayTimes(targetIndex) = computeAcfDecayTime(data(:, targetIndex), parameters);
else
acfDecayTimes(targetIndex) = computeAcfDecayTime(squeeze(data(:, targetIndex, :)), parameters);
end
k_history = parameters.k;
end
if (parameters.verbose)
fprintf('Beginning greedy selection of parents for %d with threshold %.6f\n', targetIndex, threshold);
end
% LOOP 1 -- iterating over rounds of source selection
while ~isempty(candidates)
% Whilst there are other candidates left (and we haven't quit)
% Construct the calculator and set properties:
if (is_jidt)
if (isempty(parentSet))
% Just doing pairwise TEs this round
if (is_discrete)
if (ischar(k_history)) % assume is 'auto'
error('Autoembedding not supported for discrete calculator at the moment');
end
calc = javaObject('infodynamics.measures.discrete.TransferEntropyCalculatorDiscrete', parameters.numDiscreteBins, k_history);
else
if (is_ksg)
calc = javaObject('infodynamics.measures.continuous.kraskov.TransferEntropyCalculatorKraskov');
else
calc = javaObject('infodynamics.measures.continuous.gaussian.TransferEntropyCalculatorGaussian');
end
if (ischar(k_history)) % assume is 'auto'
calc.setProperty('AUTO_EMBED_METHOD', 'MAX_CORR_AIS_DEST_ONLY');
calc.setProperty('AUTO_EMBED_K_SEARCH_MAX', string(parameters.k_max));
else
calc.setProperty('k_HISTORY', string(k_history));
end
end
else
% We're doing conditional TEs this round, conditioned on
% length(parentSet) other sources
if (is_discrete)
calc = javaObject('infodynamics.measures.discrete.ConditionalTransferEntropyCalculatorDiscrete', parameters.numDiscreteBins, k_history, length(parentSet));
else
if (is_ksg)
calc = javaObject('infodynamics.measures.continuous.kraskov.ConditionalTransferEntropyCalculatorKraskov');
else
calc = javaObject('infodynamics.measures.continuous.gaussian.ConditionalTransferEntropyCalculatorGaussian');
end
% Assume we have saved the relevant k after the first pairwise calculation
calc.setProperty('k_HISTORY', string(k_history));
% Set up the correct number of conditionals
calc.setProperty(calc.COND_EMBED_LENGTHS_PROP_NAME, strjoin(string(ones(length(parentSet), 1)), ','));
calc.setProperty(calc.COND_EMBED_DELAYS_PROP_NAME, strjoin(string(ones(length(parentSet), 1)), ','));
calc.setProperty(calc.COND_DELAYS_PROP_NAME, strjoin(string(ones(length(parentSet), 1)), ','));
end
end
if (~hasMultipleTrials)
destination = feval(javaConverterSingleArray, data(:, targetIndex));
conditionals = feval(javaConverterMatrix, data(:, parentSet));
end
else
destination = data(:, targetIndex);
conditionals = data(:, parentSet);
end
thisRoundTEResults = zeros(1, length(candidates));
thisRoundpValResults = zeros(1, length(candidates)); % Only used for Granger
% LOOP 2 -- checking (conditional) TE from all current candidates,
% given current parent set
for sIndexInCandidates = 1:length(candidates)
% For each candidate:
if is_jidt
if (acfDecayTimes(candidates(sIndexInCandidates)) < 0)
% We haven't computed the ACF decay time for this source yet
if (~hasMultipleTrials)
acfDecayTimes(candidates(sIndexInCandidates)) = ...
computeAcfDecayTime(data(:,candidates(sIndexInCandidates)), parameters);
else
acfDecayTimes(candidates(sIndexInCandidates)) = ...
computeAcfDecayTime(squeeze(data(:,candidates(sIndexInCandidates),:)), parameters);
end
end
if is_ksg
calc.setProperty('DYN_CORR_EXCL', num2str(max(acfDecayTimes([targetIndex,candidates(sIndexInCandidates),parentSet]))));
end
% 3. Initialise the calculator for (re-)use:
calc.initialise();
% 4. Supply the sample data:
calc.setDebug(true);
if (~hasMultipleTrials)
sourceTimeSeries = feval(javaConverterSingleArray, data(:, candidates(sIndexInCandidates)));
if (isempty(parentSet))
calc.setObservations(sourceTimeSeries, destination);
else
calc.setObservations(sourceTimeSeries, destination, conditionals);
end
else
calc.startAddObservations();
for numTrial = 1 : R
destination = feval(javaConverterSingleArray, squeeze(data(:, targetIndex, numTrial)));
conditionals = feval(javaConverterMatrix, squeeze(data(:, parentSet, numTrial)));
sourceTimeSeries = feval(javaConverterSingleArray, squeeze(data(:, candidates(sIndexInCandidates),numTrial)));
if (isempty(parentSet))
calc.addObservations(sourceTimeSeries, destination);
else
calc.addObservations(sourceTimeSeries, destination, conditionals);
end
end
calc.finaliseAddObservations();
end
calc.setDebug(false);
% 5. Compute the estimate:
result = calc.computeAverageLocalOfObservations();
thisRoundTEResults(sIndexInCandidates) = result;
if ischar(k_history)
k_history = calc.getProperty('k_HISTORY');
% We autoembedded the target history if we were going to - now grab the
% determined value to use next time
calc.setProperty('AUTO_EMBED_METHOD', 'NONE');
fprintf('Target history embedding set to %s\n', k_history);
end
otherResults.k_history = k_history;
else
% Compute Granger via Oliver's toolkit:
sourceTimeSeries = data(:, candidates(sIndexInCandidates));
[result,pval] = mvgc(destination,sourceTimeSeries,conditionals, ...
'p',parameters.k,'q','1','test','modified','surrogates',parameters.numSurrogates);
thisRoundTEResults(sIndexInCandidates) = result;
thisRoundpValResults(sIndexInCandidates) = pval;
% TODO need to readout the k history here
end
end
% Find the strongest source out of these candidates:
if is_jidt
[maxTE, maxIndex] = max(thisRoundTEResults);
else
% pval is higher for more significant. We will use this to
% determine the source selection, since it corrects the raw
% measure values for autocorrelation in this toolkit
[pValue, maxIndex] = max(thisRoundpValResults);
maxTE = thisRoundTEResults(maxIndex);
end
strongestSource = candidates(maxIndex);
% fprintf('Strongest source is %d with TE %.5f\n', strongestSource, maxTE);
% Check significance of this source (first need to set up its calculator again):
if is_jidt
calc.initialise();
if (~hasMultipleTrials)
sourceTimeSeries = feval(javaConverterSingleArray, data(:, strongestSource));
if (isempty(parentSet))
calc.setObservations(sourceTimeSeries, destination);
else
calc.setObservations(sourceTimeSeries, destination, conditionals);
end
else
calc.startAddObservations();
for numTrial = 1 : R
destination = feval(javaConverterSingleArray, squeeze(data(:, targetIndex, numTrial)));
conditionals = feval(javaConverterMatrix, squeeze(data(:, parentSet, numTrial)));
sourceTimeSeries = feval(javaConverterSingleArray, squeeze(data(:, strongestSource, numTrial)));
if (isempty(parentSet))
calc.addObservations(sourceTimeSeries, destination);
else
calc.addObservations(sourceTimeSeries, destination, conditionals);
end
end
calc.finaliseAddObservations();
end
maxTE = calc.computeAverageLocalOfObservations();
if (parameters.numSurrogates == 0)
measDist = calc.computeSignificance();
else
measDist = calc.computeSignificance(parameters.numSurrogates);
end
pValue = 1 - measDist.pValue; % Complementing the p value so it's the proportion of null the measure is greater than
end
if (pValue > 1 - threshold)
% We add this source to the parent set
if (parameters.verbose)
fprintf('Selected source %d, with TE(%d->%d | %s)=%.5f, p-value=%.5f (conditioning on %d parents)\n', ...
strongestSource, strongestSource, targetIndex, strjoin(string(parentSet)), maxTE, pValue, length(parentSet));
end
candidates(maxIndex) = []; % Remove this source from the candidates
parentSet = [parentSet, strongestSource];
results = [results, maxTE];
pValues = [pValues, pValue];
else
% Source was not significant, so we quit
if (parameters.verbose)
fprintf('-- Max TE source %d was not significant (TE(%d->%d | %s)=%.5f, p-value=%.6f (threshold %.6f)), quitting\n', ...
strongestSource, strongestSource, targetIndex, strjoin(string(parentSet)), maxTE, pValue, 1-threshold);
end
break;
end
toc
end
if (parameters.verbose)
fprintf('\nFinal selected parents: %s -> %d\n', strjoin(string(parentSet)), targetIndex);
end
end
% Returns the first time the ACF dips below 1/e for the series x,
% or if there are multiple series for x we take the mean across all of them
function acfDecayTime = computeAcfDecayTime(x, parameters)
numTrials = size(x,2);
acfDecayTimes = zeros(1,numTrials);
for trial = 1 : numTrials
[acfValues, ~] = autocorr(x(:, trial), 'NumLags', parameters.maxDynCorrExclLags);
acfDecayTimes(trial) = parameters.maxDynCorrExclLags; % Default is max value
for t = 1 : parameters.maxDynCorrExclLags
if (acfValues(t) < exp(-1))
acfDecayTimes(trial) = t;
break;
end
end
end
acfDecayTime = round(mean(acfDecayTimes));
end

740
demos/octave/FlockingAnalysis/NetLogoExample/Flocking.nlogo
View File

@ -1,740 +0,0 @@
;; Adapted by J.T. Lizier from the original CC BY NC SA 3.0 licensed Flocking model from NetLogo
;; and released under the same license.
;; See license info on the Info tab (and below if viewing text file)
globals [
headings-file
positionsx-file
positionsy-file
sorted-agents
]
turtles-own [
flockmates ;; agentset of nearby turtles
nearest-neighbor ;; closest one of our flockmates
]
to setup
clear-all
create-turtles population
[ set color yellow - 2 + random 7 ;; random shades look nice
set size 1.5 ;; easier to see
setxy random-xcor random-ycor
set flockmates no-turtles ]
reset-ticks
;; Set up for file output of data:
set headings-file "headings.txt"
set positionsx-file "positionsx.txt"
set positionsy-file "positionsy.txt"
if file-exists? headings-file [
file-delete headings-file ;; clear file
]
if file-exists? positionsx-file [
file-delete positionsx-file ;; clear file
]
if file-exists? positionsy-file [
file-delete positionsy-file ;; clear file
]
set sorted-agents sort turtles
end
to go
ask turtles [ flock ]
;; the following line is used to make the turtles
;; animate more smoothly.
repeat 5 [ ask turtles [ fd 0.2 ] display ]
;; for greater efficiency, at the expense of smooth
;; animation, substitute the following line instead:
;; ask turtles [ fd 1 ]
capture-data ;; write the current position+heading data of turtles to file
tick
end
to flock ;; turtle procedure
find-flockmates
if any? flockmates
[ find-nearest-neighbor
ifelse distance nearest-neighbor < minimum-separation
[ separate ]
[ align
cohere ] ]
end
to find-flockmates ;; turtle procedure
set flockmates other turtles in-radius vision
end
to find-nearest-neighbor ;; turtle procedure
set nearest-neighbor min-one-of flockmates [distance myself]
end
;;; SEPARATE
to separate ;; turtle procedure
turn-away ([heading] of nearest-neighbor) max-separate-turn
end
;;; ALIGN
to align ;; turtle procedure
turn-towards average-flockmate-heading max-align-turn
end
to-report average-flockmate-heading ;; turtle procedure
;; We can't just average the heading variables here.
;; For example, the average of 1 and 359 should be 0,
;; not 180. So we have to use trigonometry.
let x-component sum [dx] of flockmates
let y-component sum [dy] of flockmates
ifelse x-component = 0 and y-component = 0
[ report heading ]
[ report atan x-component y-component ]
end
;;; COHERE
to cohere ;; turtle procedure
turn-towards average-heading-towards-flockmates max-cohere-turn
end
to-report average-heading-towards-flockmates ;; turtle procedure
;; "towards myself" gives us the heading from the other turtle
;; to me, but we want the heading from me to the other turtle,
;; so we add 180
let x-component mean [sin (towards myself + 180)] of flockmates
let y-component mean [cos (towards myself + 180)] of flockmates
ifelse x-component = 0 and y-component = 0
[ report heading ]
[ report atan x-component y-component ]
end
;;; HELPER PROCEDURES
to turn-towards [new-heading max-turn] ;; turtle procedure
turn-at-most (subtract-headings new-heading heading) max-turn
end
to turn-away [new-heading max-turn] ;; turtle procedure
turn-at-most (subtract-headings heading new-heading) max-turn
end
;; turn right by "turn" degrees (or left if "turn" is negative),
;; but never turn more than "max-turn" degrees
to turn-at-most [turn max-turn] ;; turtle procedure
ifelse abs turn > max-turn
[ ifelse turn > 0
[ rt max-turn ]
[ lt max-turn ] ]
[ rt turn ]
end
to capture-data
;; To do a one-off:
;; file-open "headings.txt"
;; ;; Does not select turtles in order:
;; ask turtles [ file-print heading ]
;; file-close
if ticks > 4000
[ stop ]
;; Write the current headings
file-open headings-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write heading
]
]
file-print " " ;; Terminate the line
file-close
;; Write the x positions
file-open positionsx-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write xcor
]
]
file-print " " ;; Terminate the line
file-close
;; Write the y positions
file-open positionsy-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write ycor
]
]
file-print " " ;; Terminate the line
file-close
end
; Copyright 1998 Uri Wilensky.
; See Info tab for full copyright and license.
@#$#@#$#@
GRAPHICS-WINDOW
250
10
755
516
-1
-1
7.0
1
10
1
1
1
0
1
1
1
-35
35
-35
35
1
1
1
ticks
30.0
BUTTON
39
93
116
126
NIL
setup
NIL
1
T
OBSERVER
NIL
NIL
NIL
NIL
1
BUTTON
122
93
203
126
NIL
go
T
1
T
OBSERVER
NIL
NIL
NIL
NIL
0
SLIDER
9
51
232
84
population
population
1.0
1000.0
300.0
1.0
1
NIL
HORIZONTAL
SLIDER
4
217
237
250
max-align-turn
max-align-turn
0.0
20.0
5.0
0.25
1
degrees
HORIZONTAL
SLIDER
4
251
237
284
max-cohere-turn
max-cohere-turn
0.0
20.0
3.0
0.25
1
degrees
HORIZONTAL
SLIDER
4
285
237
318
max-separate-turn
max-separate-turn
0.0
20.0
1.5
0.25
1
degrees
HORIZONTAL
SLIDER
9
135
232
168
vision
vision
0.0
10.0
3.0
0.5
1
patches
HORIZONTAL
SLIDER
9
169
232
202
minimum-separation
minimum-separation
0.0
5.0
1.0
0.25
1
patches
HORIZONTAL
@#$#@#$#@
## WHAT IS IT?
This model is an attempt to mimic the flocking of birds. (The resulting motion also resembles schools of fish.) The flocks that appear in this model are not created or led in any way by special leader birds. Rather, each bird is following exactly the same set of rules, from which flocks emerge.
## HOW IT WORKS
The birds follow three rules: "alignment", "separation", and "cohesion".
"Alignment" means that a bird tends to turn so that it is moving in the same direction that nearby birds are moving.
"Separation" means that a bird will turn to avoid another bird which gets too close.
"Cohesion" means that a bird will move towards other nearby birds (unless another bird is too close).
When two birds are too close, the "separation" rule overrides the other two, which are deactivated until the minimum separation is achieved.
The three rules affect only the bird's heading. Each bird always moves forward at the same constant speed.
## HOW TO USE IT
First, determine the number of birds you want in the simulation and set the POPULATION slider to that value. Press SETUP to create the birds, and press GO to have them start flying around.
The default settings for the sliders will produce reasonably good flocking behavior. However, you can play with them to get variations:
Three TURN-ANGLE sliders control the maximum angle a bird can turn as a result of each rule.
VISION is the distance that each bird can see 360 degrees around it.
## THINGS TO NOTICE
Central to the model is the observation that flocks form without a leader.
There are no random numbers used in this model, except to position the birds initially. The fluid, lifelike behavior of the birds is produced entirely by deterministic rules.
Also, notice that each flock is dynamic. A flock, once together, is not guaranteed to keep all of its members. Why do you think this is?
After running the model for a while, all of the birds have approximately the same heading. Why?
Sometimes a bird breaks away from its flock. How does this happen? You may need to slow down the model or run it step by step in order to observe this phenomenon.
## THINGS TO TRY
Play with the sliders to see if you can get tighter flocks, looser flocks, fewer flocks, more flocks, more or less splitting and joining of flocks, more or less rearranging of birds within flocks, etc.
You can turn off a rule entirely by setting that rule's angle slider to zero. Is one rule by itself enough to produce at least some flocking? What about two rules? What's missing from the resulting behavior when you leave out each rule?
Will running the model for a long time produce a static flock? Or will the birds never settle down to an unchanging formation? Remember, there are no random numbers used in this model.
## EXTENDING THE MODEL
Currently the birds can "see" all around them. What happens if birds can only see in front of them? The `in-cone` primitive can be used for this.
Is there some way to get V-shaped flocks, like migrating geese?
What happens if you put walls around the edges of the world that the birds can't fly into?
Can you get the birds to fly around obstacles in the middle of the world?
What would happen if you gave the birds different velocities? For example, you could make birds that are not near other birds fly faster to catch up to the flock. Or, you could simulate the diminished air resistance that birds experience when flying together by making them fly faster when in a group.
Are there other interesting ways you can make the birds different from each other? There could be random variation in the population, or you could have distinct "species" of bird.
## NETLOGO FEATURES
Notice the need for the `subtract-headings` primitive and special procedure for averaging groups of headings. Just subtracting the numbers, or averaging the numbers, doesn't give you the results you'd expect, because of the discontinuity where headings wrap back to 0 once they reach 360.
## RELATED MODELS
* Moths
* Flocking Vee Formation
* Flocking - Alternative Visualizations
## CREDITS AND REFERENCES
(Note: This is an adaptation by J.T. Lizier of the original Flocking model distributed in the NetLogo Models Library, under CC BY NC SA license (see below))
This model is inspired by the Boids simulation invented by Craig Reynolds. The algorithm we use here is roughly similar to the original Boids algorithm, but it is not the same. The exact details of the algorithm tend not to matter very much -- as long as you have alignment, separation, and cohesion, you will usually get flocking behavior resembling that produced by Reynolds' original model. Information on Boids is available at http://www.red3d.com/cwr/boids/.
## HOW TO CITE
If you mention this model or the NetLogo software in a publication, we ask that you include the citations below.
For the model itself:
* Wilensky, U. (1998). NetLogo Flocking model. http://ccl.northwestern.edu/netlogo/models/Flocking. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.
Please cite the NetLogo software as:
* Wilensky, U. (1999). NetLogo. http://ccl.northwestern.edu/netlogo/. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.
## COPYRIGHT AND LICENSE
Copyright 1998 Uri Wilensky.
![CC BY-NC-SA 3.0](http://ccl.northwestern.edu/images/creativecommons/byncsa.png)
This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License. To view a copy of this license, visit https://creativecommons.org/licenses/by-nc-sa/3.0/ or send a letter to Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.
Commercial licenses are also available. To inquire about commercial licenses, please contact Uri Wilensky at uri@northwestern.edu.
This model was created as part of the project: CONNECTED MATHEMATICS: MAKING SENSE OF COMPLEX PHENOMENA THROUGH BUILDING OBJECT-BASED PARALLEL MODELS (OBPML). The project gratefully acknowledges the support of the National Science Foundation (Applications of Advanced Technologies Program) -- grant numbers RED #9552950 and REC #9632612.
This model was converted to NetLogo as part of the projects: PARTICIPATORY SIMULATIONS: NETWORK-BASED DESIGN FOR SYSTEMS LEARNING IN CLASSROOMS and/or INTEGRATED SIMULATION AND MODELING ENVIRONMENT. The project gratefully acknowledges the support of the National Science Foundation (REPP & ROLE programs) -- grant numbers REC #9814682 and REC-0126227. Converted from StarLogoT to NetLogo, 2002.
<!-- 1998 2002 -->
@#$#@#$#@
default
true
0
Polygon -7500403 true true 150 5 40 250 150 205 260 250
airplane
true
0
Polygon -7500403 true true 150 0 135 15 120 60 120 105 15 165 15 195 120 180 135 240 105 270 120 285 150 270 180 285 210 270 165 240 180 180 285 195 285 165 180 105 180 60 165 15
arrow
true
0
Polygon -7500403 true true 150 0 0 150 105 150 105 293 195 293 195 150 300 150
box
false
0
Polygon -7500403 true true 150 285 285 225 285 75 150 135
Polygon -7500403 true true 150 135 15 75 150 15 285 75
Polygon -7500403 true true 15 75 15 225 150 285 150 135
Line -16777216 false 150 285 150 135
Line -16777216 false 150 135 15 75
Line -16777216 false 150 135 285 75
bug
true
0
Circle -7500403 true true 96 182 108
Circle -7500403 true true 110 127 80
Circle -7500403 true true 110 75 80
Line -7500403 true 150 100 80 30
Line -7500403 true 150 100 220 30
butterfly
true
0
Polygon -7500403 true true 150 165 209 199 225 225 225 255 195 270 165 255 150 240
Polygon -7500403 true true 150 165 89 198 75 225 75 255 105 270 135 255 150 240
Polygon -7500403 true true 139 148 100 105 55 90 25 90 10 105 10 135 25 180 40 195 85 194 139 163
Polygon -7500403 true true 162 150 200 105 245 90 275 90 290 105 290 135 275 180 260 195 215 195 162 165
Polygon -16777216 true false 150 255 135 225 120 150 135 120 150 105 165 120 180 150 165 225
Circle -16777216 true false 135 90 30
Line -16777216 false 150 105 195 60
Line -16777216 false 150 105 105 60
car
false
0
Polygon -7500403 true true 300 180 279 164 261 144 240 135 226 132 213 106 203 84 185 63 159 50 135 50 75 60 0 150 0 165 0 225 300 225 300 180
Circle -16777216 true false 180 180 90
Circle -16777216 true false 30 180 90
Polygon -16777216 true false 162 80 132 78 134 135 209 135 194 105 189 96 180 89
Circle -7500403 true true 47 195 58
Circle -7500403 true true 195 195 58
circle
false
0
Circle -7500403 true true 0 0 300
circle 2
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240
cow
false
0
Polygon -7500403 true true 200 193 197 249 179 249 177 196 166 187 140 189 93 191 78 179 72 211 49 209 48 181 37 149 25 120 25 89 45 72 103 84 179 75 198 76 252 64 272 81 293 103 285 121 255 121 242 118 224 167
Polygon -7500403 true true 73 210 86 251 62 249 48 208
Polygon -7500403 true true 25 114 16 195 9 204 23 213 25 200 39 123
cylinder
false
0
Circle -7500403 true true 0 0 300
dot
false
0
Circle -7500403 true true 90 90 120
face happy
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 255 90 239 62 213 47 191 67 179 90 203 109 218 150 225 192 218 210 203 227 181 251 194 236 217 212 240
face neutral
false
0
Circle -7500403 true true 8 7 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Rectangle -16777216 true false 60 195 240 225
face sad
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 168 90 184 62 210 47 232 67 244 90 220 109 205 150 198 192 205 210 220 227 242 251 229 236 206 212 183
fish
false
0
Polygon -1 true false 44 131 21 87 15 86 0 120 15 150 0 180 13 214 20 212 45 166
Polygon -1 true false 135 195 119 235 95 218 76 210 46 204 60 165
Polygon -1 true false 75 45 83 77 71 103 86 114 166 78 135 60
Polygon -7500403 true true 30 136 151 77 226 81 280 119 292 146 292 160 287 170 270 195 195 210 151 212 30 166
Circle -16777216 true false 215 106 30
flag
false
0
Rectangle -7500403 true true 60 15 75 300
Polygon -7500403 true true 90 150 270 90 90 30
Line -7500403 true 75 135 90 135
Line -7500403 true 75 45 90 45
flower
false
0
Polygon -10899396 true false 135 120 165 165 180 210 180 240 150 300 165 300 195 240 195 195 165 135
Circle -7500403 true true 85 132 38
Circle -7500403 true true 130 147 38
Circle -7500403 true true 192 85 38
Circle -7500403 true true 85 40 38
Circle -7500403 true true 177 40 38
Circle -7500403 true true 177 132 38
Circle -7500403 true true 70 85 38
Circle -7500403 true true 130 25 38
Circle -7500403 true true 96 51 108
Circle -16777216 true false 113 68 74
Polygon -10899396 true false 189 233 219 188 249 173 279 188 234 218
Polygon -10899396 true false 180 255 150 210 105 210 75 240 135 240
house
false
0
Rectangle -7500403 true true 45 120 255 285
Rectangle -16777216 true false 120 210 180 285
Polygon -7500403 true true 15 120 150 15 285 120
Line -16777216 false 30 120 270 120
leaf
false
0
Polygon -7500403 true true 150 210 135 195 120 210 60 210 30 195 60 180 60 165 15 135 30 120 15 105 40 104 45 90 60 90 90 105 105 120 120 120 105 60 120 60 135 30 150 15 165 30 180 60 195 60 180 120 195 120 210 105 240 90 255 90 263 104 285 105 270 120 285 135 240 165 240 180 270 195 240 210 180 210 165 195
Polygon -7500403 true true 135 195 135 240 120 255 105 255 105 285 135 285 165 240 165 195
line
true
0
Line -7500403 true 150 0 150 300
line half
true
0
Line -7500403 true 150 0 150 150
pentagon
false
0
Polygon -7500403 true true 150 15 15 120 60 285 240 285 285 120
person
false
0
Circle -7500403 true true 110 5 80
Polygon -7500403 true true 105 90 120 195 90 285 105 300 135 300 150 225 165 300 195 300 210 285 180 195 195 90
Rectangle -7500403 true true 127 79 172 94
Polygon -7500403 true true 195 90 240 150 225 180 165 105
Polygon -7500403 true true 105 90 60 150 75 180 135 105
plant
false
0
Rectangle -7500403 true true 135 90 165 300
Polygon -7500403 true true 135 255 90 210 45 195 75 255 135 285
Polygon -7500403 true true 165 255 210 210 255 195 225 255 165 285
Polygon -7500403 true true 135 180 90 135 45 120 75 180 135 210
Polygon -7500403 true true 165 180 165 210 225 180 255 120 210 135
Polygon -7500403 true true 135 105 90 60 45 45 75 105 135 135
Polygon -7500403 true true 165 105 165 135 225 105 255 45 210 60
Polygon -7500403 true true 135 90 120 45 150 15 180 45 165 90
square
false
0
Rectangle -7500403 true true 30 30 270 270
square 2
false
0
Rectangle -7500403 true true 30 30 270 270
Rectangle -16777216 true false 60 60 240 240
star
false
0
Polygon -7500403 true true 151 1 185 108 298 108 207 175 242 282 151 216 59 282 94 175 3 108 116 108
target
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240
Circle -7500403 true true 60 60 180
Circle -16777216 true false 90 90 120
Circle -7500403 true true 120 120 60
tree
false
0
Circle -7500403 true true 118 3 94
Rectangle -6459832 true false 120 195 180 300
Circle -7500403 true true 65 21 108
Circle -7500403 true true 116 41 127
Circle -7500403 true true 45 90 120
Circle -7500403 true true 104 74 152
triangle
false
0
Polygon -7500403 true true 150 30 15 255 285 255
triangle 2
false
0
Polygon -7500403 true true 150 30 15 255 285 255
Polygon -16777216 true false 151 99 225 223 75 224
truck
false
0
Rectangle -7500403 true true 4 45 195 187
Polygon -7500403 true true 296 193 296 150 259 134 244 104 208 104 207 194
Rectangle -1 true false 195 60 195 105
Polygon -16777216 true false 238 112 252 141 219 141 218 112
Circle -16777216 true false 234 174 42
Rectangle -7500403 true true 181 185 214 194
Circle -16777216 true false 144 174 42
Circle -16777216 true false 24 174 42
Circle -7500403 false true 24 174 42
Circle -7500403 false true 144 174 42
Circle -7500403 false true 234 174 42
turtle
true
0
Polygon -10899396 true false 215 204 240 233 246 254 228 266 215 252 193 210
Polygon -10899396 true false 195 90 225 75 245 75 260 89 269 108 261 124 240 105 225 105 210 105
Polygon -10899396 true false 105 90 75 75 55 75 40 89 31 108 39 124 60 105 75 105 90 105
Polygon -10899396 true false 132 85 134 64 107 51 108 17 150 2 192 18 192 52 169 65 172 87
Polygon -10899396 true false 85 204 60 233 54 254 72 266 85 252 107 210
Polygon -7500403 true true 119 75 179 75 209 101 224 135 220 225 175 261 128 261 81 224 74 135 88 99
wheel
false
0
Circle -7500403 true true 3 3 294
Circle -16777216 true false 30 30 240
Line -7500403 true 150 285 150 15
Line -7500403 true 15 150 285 150
Circle -7500403 true true 120 120 60
Line -7500403 true 216 40 79 269
Line -7500403 true 40 84 269 221
Line -7500403 true 40 216 269 79
Line -7500403 true 84 40 221 269
x
false
0
Polygon -7500403 true true 270 75 225 30 30 225 75 270
Polygon -7500403 true true 30 75 75 30 270 225 225 270
@#$#@#$#@
NetLogo 6.0
@#$#@#$#@
set population 200
setup
repeat 200 [ go ]
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
default
0.0
-0.2 0 0.0 1.0
0.0 1 1.0 0.0
0.2 0 0.0 1.0
link direction
true
0
Line -7500403 true 150 150 90 180
Line -7500403 true 150 150 210 180
@#$#@#$#@
0
@#$#@#$#@

810
demos/octave/FlockingAnalysis/NetLogoExample/Flocking_withSPEED.nlogo
View File

@ -1,810 +0,0 @@
;; Adapted by J.T. Lizier and S. Sekhar from the original CC BY NC SA 3.0 licensed Flocking model from NetLogo
;; and released under the same license.
;; See license info on the Info tab (and below if viewing text file)
globals [
headings-file
positionsx-file
positionsy-file
sorted-agents
]
turtles-own [
flockmates ;; agentset of nearby turtles
nearest-neighbor ;; closest one of our flockmates
speed
]
to setup
clear-all
create-turtles population
[ set color yellow - 2 + random 7 ;; random shades look nice
set size 1.5 ;; easier to see
setxy random-xcor random-ycor
set flockmates no-turtles
set speed ((random-float 0.5) + 0.25)
]
reset-ticks
;; Set up for file output of data:
set headings-file "headings.txt"
set positionsx-file "positionsx.txt"
set positionsy-file "positionsy.txt"
if file-exists? headings-file [
file-delete headings-file ;; clear file
]
if file-exists? positionsx-file [
file-delete positionsx-file ;; clear file
]
if file-exists? positionsy-file [
file-delete positionsy-file ;; clear file
]
set sorted-agents sort turtles
end
to go
ask turtles [
fd speed
flock ]
if ticks > 0 [ ;; You can take data after a larger number of ticks to remove a transient
capture-data ];; write the current position+heading data of turtles to file
tick
end
to flock ;; turtle procedure
find-flockmates
if any? flockmates
[ find-nearest-neighbor
ifelse distance nearest-neighbor < minimum-separation
[ separate ]
[ align
cohere
speedmatch
] ]
end
to find-flockmates ;; turtle procedure
set flockmates other turtles in-radius vision
end
to find-nearest-neighbor ;; turtle procedure
set nearest-neighbor min-one-of flockmates [distance myself]
end
;;; SEPARATE
to separate ;; turtle procedure
turn-away ([heading] of nearest-neighbor) max-separate-turn
end
;;;SPEEDMATCH
to speedmatch
let speed-change (abs(speed - average-flockmate-speed) * 0.80 )
( ifelse average-flockmate-speed > speed
[ set speed (speed + speed-change) ]
[ set speed (speed - speed-change) ] )
set speed (speed + speed-change * ((random-float 0.2) - 0.1))
end
to-report average-flockmate-speed
let speed-component mean [speed] of flockmates
report speed-component
end
;;; ALIGN
to align ;; turtle procedure
turn-towards average-flockmate-heading max-align-turn
end
to-report average-flockmate-heading ;; turtle procedure
;; We can't just average the heading variables here.
;; For example, the average of 1 and 359 should be 0,
;; not 180. So we have to use trigonometry.
let x-component sum [dx] of flockmates
let y-component sum [dy] of flockmates
ifelse x-component = 0 and y-component = 0
[ report heading ]
[ report atan x-component y-component ]
end
;;; COHERE
to cohere ;; turtle procedure
turn-towards average-heading-towards-flockmates max-cohere-turn
end
to-report average-heading-towards-flockmates ;; turtle procedure
;; "towards myself" gives us the heading from the other turtle
;; to me, but we want the heading from me to the other turtle,
;; so we add 180
let x-component mean [sin (towards myself + 180)] of flockmates
let y-component mean [cos (towards myself + 180)] of flockmates
ifelse x-component = 0 and y-component = 0
[ report heading ]
[ report atan x-component y-component ]
end
;;; HELPER PROCEDURES
to turn-towards [new-heading max-turn] ;; turtle procedure
turn-at-most (subtract-headings new-heading heading) max-turn
end
to turn-away [new-heading max-turn] ;; turtle procedure
turn-at-most (subtract-headings heading new-heading) max-turn
end
;; turn right by "turn" degrees (or left if "turn" is negative),
;; but never turn more than "max-turn" degrees
to turn-at-most [turn max-turn] ;; turtle procedure
ifelse abs turn > max-turn
[ ifelse turn > 0
[ rt max-turn ]
[ lt max-turn ] ]
[ rt turn ]
end
to capture-data
;; To do a one-off:
;; file-open "headings.txt"
;; ;; Does not select turtles in order:
;; ask turtles [ file-print heading ]
;; file-close
if ticks > 4000
[ stop ]
;; Write the current headings
file-open headings-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write heading
]
]
file-print " " ;; Terminate the line
file-close
;; Write the x positions
file-open positionsx-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write xcor
]
]
file-print " " ;; Terminate the line
file-close
;; Write the y positions
file-open positionsy-file ;; Opening file for writing
foreach sorted-agents [
[the-turtle] ->
ask the-turtle [
file-write ycor
]
]
file-print " " ;; Terminate the line
file-close
end
; Copyright 1998 Uri Wilensky.
; See Info tab for full copyright and license.
@#$#@#$#@
GRAPHICS-WINDOW
250
10
755
516
-1
-1
7.0
1
10
1
1
1
0
1
1
1
-35
35
-35
35
1
1
1
ticks
30.0
BUTTON
39
93
116
126
NIL
setup
NIL
1
T
OBSERVER
NIL
NIL
NIL
NIL
1
BUTTON
122
93
203
126
NIL
go
T
1
T
OBSERVER
NIL
NIL
NIL
NIL
0
SLIDER
9
51
232
84
population
population
1.0
1000.0
300.0
1.0
1
NIL
HORIZONTAL
SLIDER
4
217
237
250
max-align-turn
max-align-turn
0.0
20.0
5.0
0.25
1
degrees
HORIZONTAL
SLIDER
4
251
237
284
max-cohere-turn
max-cohere-turn
0.0
20.0
3.0
0.25
1
degrees
HORIZONTAL
SLIDER
4
285
237
318
max-separate-turn
max-separate-turn
0.0
20.0
1.5
0.25
1
degrees
HORIZONTAL
SLIDER
9
135
232
168
vision
vision
0.0
10.0
3.0
0.5
1
patches
HORIZONTAL
SLIDER
9
169
232
202
minimum-separation
minimum-separation
0.0
5.0
1.0
0.25
1
patches
HORIZONTAL
PLOT
757
42
1051
447
plot 1
Speed
Turtles
0.2
0.8
0.0
100.0
true
false
"\n\nset-histogram-num-bars 20" ""
PENS
"default" 1.0 1 -16777216 true "" "histogram [speed] of turtles\n"
MONITOR
30
434
173
479
NIL
mean [speed] of turtles
10
1
11
MONITOR
95
376
216
421
NIL
[speed] of turtle 25
10
1
11
MONITOR
23
329
151
374
NIL
[speed] of turtle 213
10
1
11
@#$#@#$#@
## WHAT IS IT?
This model is an attempt to mimic the flocking of birds. (The resulting motion also resembles schools of fish.) The flocks that appear in this model are not created or led in any way by special leader birds. Rather, each bird is following exactly the same set of rules, from which flocks emerge.
## HOW IT WORKS
The birds follow three rules: "alignment", "separation", and "cohesion".
"Alignment" means that a bird tends to turn so that it is moving in the same direction that nearby birds are moving.
"Separation" means that a bird will turn to avoid another bird which gets too close.
"Cohesion" means that a bird will move towards other nearby birds (unless another bird is too close).
When two birds are too close, the "separation" rule overrides the other two, which are deactivated until the minimum separation is achieved.
The three rules affect only the bird's heading. Each bird always moves forward at the same constant speed.
## HOW TO USE IT
First, determine the number of birds you want in the simulation and set the POPULATION slider to that value. Press SETUP to create the birds, and press GO to have them start flying around.
The default settings for the sliders will produce reasonably good flocking behavior. However, you can play with them to get variations:
Three TURN-ANGLE sliders control the maximum angle a bird can turn as a result of each rule.
VISION is the distance that each bird can see 360 degrees around it.
## THINGS TO NOTICE
Central to the model is the observation that flocks form without a leader.
There are no random numbers used in this model, except to position the birds initially. The fluid, lifelike behavior of the birds is produced entirely by deterministic rules.
Also, notice that each flock is dynamic. A flock, once together, is not guaranteed to keep all of its members. Why do you think this is?
After running the model for a while, all of the birds have approximately the same heading. Why?
Sometimes a bird breaks away from its flock. How does this happen? You may need to slow down the model or run it step by step in order to observe this phenomenon.
## THINGS TO TRY
Play with the sliders to see if you can get tighter flocks, looser flocks, fewer flocks, more flocks, more or less splitting and joining of flocks, more or less rearranging of birds within flocks, etc.
You can turn off a rule entirely by setting that rule's angle slider to zero. Is one rule by itself enough to produce at least some flocking? What about two rules? What's missing from the resulting behavior when you leave out each rule?
Will running the model for a long time produce a static flock? Or will the birds never settle down to an unchanging formation? Remember, there are no random numbers used in this model.
## EXTENDING THE MODEL
Currently the birds can "see" all around them. What happens if birds can only see in front of them? The `in-cone` primitive can be used for this.
Is there some way to get V-shaped flocks, like migrating geese?
What happens if you put walls around the edges of the world that the birds can't fly into?
Can you get the birds to fly around obstacles in the middle of the world?
What would happen if you gave the birds different velocities? For example, you could make birds that are not near other birds fly faster to catch up to the flock. Or, you could simulate the diminished air resistance that birds experience when flying together by making them fly faster when in a group.
Are there other interesting ways you can make the birds different from each other? There could be random variation in the population, or you could have distinct "species" of bird.
## NETLOGO FEATURES
Notice the need for the `subtract-headings` primitive and special procedure for averaging groups of headings. Just subtracting the numbers, or averaging the numbers, doesn't give you the results you'd expect, because of the discontinuity where headings wrap back to 0 once they reach 360.
## RELATED MODELS
* Moths
* Flocking Vee Formation
* Flocking - Alternative Visualizations
## CREDITS AND REFERENCES
(Note: This is an adaptation by J.T. Lizier of the original Flocking model distributed in the NetLogo Models Library, under CC BY NC SA license (see below))
This model is inspired by the Boids simulation invented by Craig Reynolds. The algorithm we use here is roughly similar to the original Boids algorithm, but it is not the same. The exact details of the algorithm tend not to matter very much -- as long as you have alignment, separation, and cohesion, you will usually get flocking behavior resembling that produced by Reynolds' original model. Information on Boids is available at http://www.red3d.com/cwr/boids/.
## HOW TO CITE
If you mention this model or the NetLogo software in a publication, we ask that you include the citations below.
For the model itself:
* Wilensky, U. (1998). NetLogo Flocking model. http://ccl.northwestern.edu/netlogo/models/Flocking. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.
Please cite the NetLogo software as:
* Wilensky, U. (1999). NetLogo. http://ccl.northwestern.edu/netlogo/. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL.
## COPYRIGHT AND LICENSE
Copyright 1998 Uri Wilensky.
![CC BY-NC-SA 3.0](http://ccl.northwestern.edu/images/creativecommons/byncsa.png)
This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License. To view a copy of this license, visit https://creativecommons.org/licenses/by-nc-sa/3.0/ or send a letter to Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.
Commercial licenses are also available. To inquire about commercial licenses, please contact Uri Wilensky at uri@northwestern.edu.
This model was created as part of the project: CONNECTED MATHEMATICS: MAKING SENSE OF COMPLEX PHENOMENA THROUGH BUILDING OBJECT-BASED PARALLEL MODELS (OBPML). The project gratefully acknowledges the support of the National Science Foundation (Applications of Advanced Technologies Program) -- grant numbers RED #9552950 and REC #9632612.
This model was converted to NetLogo as part of the projects: PARTICIPATORY SIMULATIONS: NETWORK-BASED DESIGN FOR SYSTEMS LEARNING IN CLASSROOMS and/or INTEGRATED SIMULATION AND MODELING ENVIRONMENT. The project gratefully acknowledges the support of the National Science Foundation (REPP & ROLE programs) -- grant numbers REC #9814682 and REC-0126227. Converted from StarLogoT to NetLogo, 2002.
<!-- 1998 2002 -->
@#$#@#$#@
default
true
0
Polygon -7500403 true true 150 5 40 250 150 205 260 250
airplane
true
0
Polygon -7500403 true true 150 0 135 15 120 60 120 105 15 165 15 195 120 180 135 240 105 270 120 285 150 270 180 285 210 270 165 240 180 180 285 195 285 165 180 105 180 60 165 15
arrow
true
0
Polygon -7500403 true true 150 0 0 150 105 150 105 293 195 293 195 150 300 150
box
false
0
Polygon -7500403 true true 150 285 285 225 285 75 150 135
Polygon -7500403 true true 150 135 15 75 150 15 285 75
Polygon -7500403 true true 15 75 15 225 150 285 150 135
Line -16777216 false 150 285 150 135
Line -16777216 false 150 135 15 75
Line -16777216 false 150 135 285 75
bug
true
0
Circle -7500403 true true 96 182 108
Circle -7500403 true true 110 127 80
Circle -7500403 true true 110 75 80
Line -7500403 true 150 100 80 30
Line -7500403 true 150 100 220 30
butterfly
true
0
Polygon -7500403 true true 150 165 209 199 225 225 225 255 195 270 165 255 150 240
Polygon -7500403 true true 150 165 89 198 75 225 75 255 105 270 135 255 150 240
Polygon -7500403 true true 139 148 100 105 55 90 25 90 10 105 10 135 25 180 40 195 85 194 139 163
Polygon -7500403 true true 162 150 200 105 245 90 275 90 290 105 290 135 275 180 260 195 215 195 162 165
Polygon -16777216 true false 150 255 135 225 120 150 135 120 150 105 165 120 180 150 165 225
Circle -16777216 true false 135 90 30
Line -16777216 false 150 105 195 60
Line -16777216 false 150 105 105 60
car
false
0
Polygon -7500403 true true 300 180 279 164 261 144 240 135 226 132 213 106 203 84 185 63 159 50 135 50 75 60 0 150 0 165 0 225 300 225 300 180
Circle -16777216 true false 180 180 90
Circle -16777216 true false 30 180 90
Polygon -16777216 true false 162 80 132 78 134 135 209 135 194 105 189 96 180 89
Circle -7500403 true true 47 195 58
Circle -7500403 true true 195 195 58
circle
false
0
Circle -7500403 true true 0 0 300
circle 2
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240
cow
false
0
Polygon -7500403 true true 200 193 197 249 179 249 177 196 166 187 140 189 93 191 78 179 72 211 49 209 48 181 37 149 25 120 25 89 45 72 103 84 179 75 198 76 252 64 272 81 293 103 285 121 255 121 242 118 224 167
Polygon -7500403 true true 73 210 86 251 62 249 48 208
Polygon -7500403 true true 25 114 16 195 9 204 23 213 25 200 39 123
cylinder
false
0
Circle -7500403 true true 0 0 300
dot
false
0
Circle -7500403 true true 90 90 120
face happy
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 255 90 239 62 213 47 191 67 179 90 203 109 218 150 225 192 218 210 203 227 181 251 194 236 217 212 240
face neutral
false
0
Circle -7500403 true true 8 7 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Rectangle -16777216 true false 60 195 240 225
face sad
false
0
Circle -7500403 true true 8 8 285
Circle -16777216 true false 60 75 60
Circle -16777216 true false 180 75 60
Polygon -16777216 true false 150 168 90 184 62 210 47 232 67 244 90 220 109 205 150 198 192 205 210 220 227 242 251 229 236 206 212 183
fish
false
0
Polygon -1 true false 44 131 21 87 15 86 0 120 15 150 0 180 13 214 20 212 45 166
Polygon -1 true false 135 195 119 235 95 218 76 210 46 204 60 165
Polygon -1 true false 75 45 83 77 71 103 86 114 166 78 135 60
Polygon -7500403 true true 30 136 151 77 226 81 280 119 292 146 292 160 287 170 270 195 195 210 151 212 30 166
Circle -16777216 true false 215 106 30
flag
false
0
Rectangle -7500403 true true 60 15 75 300
Polygon -7500403 true true 90 150 270 90 90 30
Line -7500403 true 75 135 90 135
Line -7500403 true 75 45 90 45
flower
false
0
Polygon -10899396 true false 135 120 165 165 180 210 180 240 150 300 165 300 195 240 195 195 165 135
Circle -7500403 true true 85 132 38
Circle -7500403 true true 130 147 38
Circle -7500403 true true 192 85 38
Circle -7500403 true true 85 40 38
Circle -7500403 true true 177 40 38
Circle -7500403 true true 177 132 38
Circle -7500403 true true 70 85 38
Circle -7500403 true true 130 25 38
Circle -7500403 true true 96 51 108
Circle -16777216 true false 113 68 74
Polygon -10899396 true false 189 233 219 188 249 173 279 188 234 218
Polygon -10899396 true false 180 255 150 210 105 210 75 240 135 240
house
false
0
Rectangle -7500403 true true 45 120 255 285
Rectangle -16777216 true false 120 210 180 285
Polygon -7500403 true true 15 120 150 15 285 120
Line -16777216 false 30 120 270 120
leaf
false
0
Polygon -7500403 true true 150 210 135 195 120 210 60 210 30 195 60 180 60 165 15 135 30 120 15 105 40 104 45 90 60 90 90 105 105 120 120 120 105 60 120 60 135 30 150 15 165 30 180 60 195 60 180 120 195 120 210 105 240 90 255 90 263 104 285 105 270 120 285 135 240 165 240 180 270 195 240 210 180 210 165 195
Polygon -7500403 true true 135 195 135 240 120 255 105 255 105 285 135 285 165 240 165 195
line
true
0
Line -7500403 true 150 0 150 300
line half
true
0
Line -7500403 true 150 0 150 150
pentagon
false
0
Polygon -7500403 true true 150 15 15 120 60 285 240 285 285 120
person
false
0
Circle -7500403 true true 110 5 80
Polygon -7500403 true true 105 90 120 195 90 285 105 300 135 300 150 225 165 300 195 300 210 285 180 195 195 90
Rectangle -7500403 true true 127 79 172 94
Polygon -7500403 true true 195 90 240 150 225 180 165 105
Polygon -7500403 true true 105 90 60 150 75 180 135 105
plant
false
0
Rectangle -7500403 true true 135 90 165 300
Polygon -7500403 true true 135 255 90 210 45 195 75 255 135 285
Polygon -7500403 true true 165 255 210 210 255 195 225 255 165 285
Polygon -7500403 true true 135 180 90 135 45 120 75 180 135 210
Polygon -7500403 true true 165 180 165 210 225 180 255 120 210 135
Polygon -7500403 true true 135 105 90 60 45 45 75 105 135 135
Polygon -7500403 true true 165 105 165 135 225 105 255 45 210 60
Polygon -7500403 true true 135 90 120 45 150 15 180 45 165 90
square
false
0
Rectangle -7500403 true true 30 30 270 270
square 2
false
0
Rectangle -7500403 true true 30 30 270 270
Rectangle -16777216 true false 60 60 240 240
star
false
0
Polygon -7500403 true true 151 1 185 108 298 108 207 175 242 282 151 216 59 282 94 175 3 108 116 108
target
false
0
Circle -7500403 true true 0 0 300
Circle -16777216 true false 30 30 240
Circle -7500403 true true 60 60 180
Circle -16777216 true false 90 90 120
Circle -7500403 true true 120 120 60
tree
false
0
Circle -7500403 true true 118 3 94
Rectangle -6459832 true false 120 195 180 300
Circle -7500403 true true 65 21 108
Circle -7500403 true true 116 41 127
Circle -7500403 true true 45 90 120
Circle -7500403 true true 104 74 152
triangle
false
0
Polygon -7500403 true true 150 30 15 255 285 255
triangle 2
false
0
Polygon -7500403 true true 150 30 15 255 285 255
Polygon -16777216 true false 151 99 225 223 75 224
truck
false
0
Rectangle -7500403 true true 4 45 195 187
Polygon -7500403 true true 296 193 296 150 259 134 244 104 208 104 207 194
Rectangle -1 true false 195 60 195 105
Polygon -16777216 true false 238 112 252 141 219 141 218 112
Circle -16777216 true false 234 174 42
Rectangle -7500403 true true 181 185 214 194
Circle -16777216 true false 144 174 42
Circle -16777216 true false 24 174 42
Circle -7500403 false true 24 174 42
Circle -7500403 false true 144 174 42
Circle -7500403 false true 234 174 42
turtle
true
0
Polygon -10899396 true false 215 204 240 233 246 254 228 266 215 252 193 210
Polygon -10899396 true false 195 90 225 75 245 75 260 89 269 108 261 124 240 105 225 105 210 105
Polygon -10899396 true false 105 90 75 75 55 75 40 89 31 108 39 124 60 105 75 105 90 105
Polygon -10899396 true false 132 85 134 64 107 51 108 17 150 2 192 18 192 52 169 65 172 87
Polygon -10899396 true false 85 204 60 233 54 254 72 266 85 252 107 210
Polygon -7500403 true true 119 75 179 75 209 101 224 135 220 225 175 261 128 261 81 224 74 135 88 99
wheel
false
0
Circle -7500403 true true 3 3 294
Circle -16777216 true false 30 30 240
Line -7500403 true 150 285 150 15
Line -7500403 true 15 150 285 150
Circle -7500403 true true 120 120 60
Line -7500403 true 216 40 79 269
Line -7500403 true 40 84 269 221
Line -7500403 true 40 216 269 79
Line -7500403 true 84 40 221 269
x
false
0
Polygon -7500403 true true 270 75 225 30 30 225 75 270
Polygon -7500403 true true 30 75 75 30 270 225 225 270
@#$#@#$#@
NetLogo 6.0
@#$#@#$#@
set population 200
setup
repeat 200 [ go ]
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
default
0.0
-0.2 0 0.0 1.0
0.0 1 1.0 0.0
0.2 0 0.0 1.0
link direction
true
0
Line -7500403 true 150 150 90 180
Line -7500403 true 150 150 210 180
@#$#@#$#@
0
@#$#@#$#@

113
demos/octave/FlockingAnalysis/NetLogoExample/loadProperties.m
View File

@ -1,113 +0,0 @@
%
% This script loads properties for transfer entropy analysis of the data from
% the NetLogo Flocking model.
%
% Author: Joseph T. Lizier, Sooraj Sekhar, 2019
%
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, 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/>.
%%
% This script loads the default properties for the transfer entropy processing:
clear('properties');
%%%%%%%%%
% FILENAMES
%%%%%%%%%
% Input data files:
% properties.files can be:
% a. a cell array of file names, e.g.: {'file1.xlsx', 'file2.xlsx'}
% b. a call to ls or ls with an argument, e.g. ls('*.xlsx')
% c. a space or tab separated character row vector of file names
% d. a character matrix of filenames (each filename on a separate row)
properties.files = 'positions%s.txt';
% Function to read in the data files:
% loadScript must point to a function .m file that accepts two arguments
% (the name of a file, and properties object) and returns [x,y,z] (z optional, only when 3D)
% data where each is an array, e.g. x(time, fishIndex) indexed first by time and second by fish index.
% Where an individual is not present at any given time step, set each x,y,z position
% to nan for these time steps. The subsequent scripts will then ignore these samples.
% Use the name of the .m file after an "@" character:
properties.loadScript = @loadseparatexy;
% Is the data returned by the loadScript 3D (true) or 2D (false)?
properties.data3d = false;
% Results file - will hold the parsed velocities / relative positions, plus the
% local transfer entropy results
properties.resultsFile = 'results.mat';
%%%%%%%%%
% PARAMETERS
%%%%%%%%%
% Distance within which to consider a pair for the info theoretic analysis (units are as per what is used in the data files)
properties.pairRange = 4; % These ones have a causal range of 3
%%%%%%%%%
% INFORMATION THEORETIC Parameters
% Only lag is used for computing lagged mutual information
% All lag, k and tau are used for transfer entropy
% k - embedding dimension of the past of the destination array.
% tau - embedding delay: time cycles separating each element in the past of the destination.
% lag - time delay between the source and target in cycles
% You can set kRange, tauRange and lagRange to ask that these are optimised by runAnalysis:
% properties.kRange = 1:10;
% properties.tauRange = 1:4;
% properties.lagRange = 1:10;
% You can also set k, tau and lag to values that generateObservations should use
% (although note that if this is called via runAnalysis then it will overwrite them):
properties.k = 1;
properties.tau = 1;
properties.lag = 1;
% Do we compute TE based on headings or speeds or both for the info calculations? (need at least one of these set to true)
properties.speedcalc = false;
properties.headingcalc = true;
% Additionally: Do we include the relative source position in the transfer entropy calculation (true), or
% only the relative source heading (false)
properties.includeSourcePositionInTransfer = false;
% Do we take relative source heading and position with respect to dest heading at that same
% time point (true, this is what we did for Crosato paper) or with respect
% to dest heading just previous to state update (false)?
% UPDATE THIS COMMENT TO MAKE MORE SENSE. Ash and Liss say false is more biologically realistic
properties.sourceWrtSameDestTime = true;
% JIDT location:
properties.jidtJarLocation = '../../../../infodynamics.jar';
% Which estimator to use.
% Valid values are 'gaussian' (linear) or 'kraskov' (non-linear)
% properties.estimator = 'gaussian';
properties.estimator = 'kraskov';
% Properties for JIDT estimators:
properties.jidt.kNNs = 4; % Number of nearest neighbours for Kraskov algorithm: just use 4 (default)
properties.jidt.autoDynamicCorrelationExclusion = true; % Exclude nearest neighbours from at least the same target transition from being included in counts for TE. This requires the Econometrics Toolbox in Matlab in order to call autocorr. Set to false if you don't have that to avoid this being called.
properties.aisNumSurrogates = 0; % Number of surrogate calculations to run for AIS (just to see the noise floor. 0 means skip)
properties.teNumSurrogates = 0; % Number of surrogate calculations to run for TE (just to see the noise floor. 0 means skip)
% Exit the analysis after calculating the AIS (i.e. don't compute TE)
properties.computeAISOnly = false;

BIN
demos/octave/FlockingAnalysis/README-Flocking.pdf
View File

Binary file not shown.

73
demos/octave/FlockingAnalysis/computeAIS.m
View File

@ -1,73 +0,0 @@
function ais = computeAIS(D, Dpast, properties)
% Computes active information storage from the pre-processed velocity data
%
% Author: Joseph T. Lizier, 2019
%
% Inputs:
% - D - target samples (may be multivariate as per generateObservations)
% - Dpast - target past samples (multivariate, and embedded up to k previous samples)
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are loaded from loadProperties.m
%
% Outputs:
% - ais - active information storage value (MI between D and Dpast)
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
% -- STEP 0 : create java object
% check out http://lizier.me/joseph/software/jidt/javadocs/v1.3/
% for the description of all classes and methods
javaaddpath(properties.jidtJarLocation); % add JIDT path
% Select correct MI class for this estimator type
if (strcmp('kraskov', properties.estimator))
MI_CLASS = 'infodynamics.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov1';
else
MI_CLASS = 'infodynamics.measures.continuous.gaussian.MutualInfoCalculatorMultiVariateGaussian';
end
% Compute the information stored in the past of the target.
% Creates java object of the given class
AIScalculator = javaObject(MI_CLASS);
% -- STEP 1 : set properties
AIScalculator.setProperty('k', num2str(properties.jidt.kNNs));
AIScalculator.setProperty('BIAS_CORRECTION', 'true'); % Used for Gaussian only
% -- STEP 2 : initialise
% here the parameters are the dimensionality of the series
% in this case taken directly from the number of columns in each variable
AIScalculator.initialise(size(D,2), size(Dpast,2));
% -- STEP 3 : add in observations
AIScalculator.setObservations(D, Dpast);
% -- STEP 4 : compute the local AIS (will be bias-corrected now for either Gaussian or KSG)
ais = AIScalculator.computeAverageLocalOfObservations(); % global (average) value
fprintf('Mean AIS_%s (k=%d,tau=%d) = %.3f\n', properties.estimator, properties.k, properties.tau, ais);
if (properties.aisNumSurrogates > 0)
% Compute the (statistical significance via) null distribution empirically (e.g. with 100 permutations),
% and use this for empirical bias correction (otherwise we're relying on analytic)
aisMeasDist = AIScalculator.computeSignificance(properties.aisNumSurrogates);
fprintf('Null distribution: %.4f +/- %.4f std dev.; p(surrogate > measured)=%.3f from %d surrogates)\n', ...
aisMeasDist.getMeanOfDistribution(), aisMeasDist.getStdOfDistribution(), ...
aisMeasDist.pValue, properties.aisNumSurrogates);
ais = ais - aisMeasDist.getMeanOfDistribution();
fprintf('Bias corrected Mean AIS_%s (k=%d,tau=%d) = %.3f\n', properties.estimator, properties.k, properties.tau, ais);
end
end

101
demos/octave/FlockingAnalysis/computeTE.m
View File

@ -1,101 +0,0 @@
function tranEntropy = computeTE(S, D, Dpast, properties)
% Computes transfer entropy from the pre-processed velocity data
%
% Author: Emanuele Crosato, Joseph T. Lizier, 2019
%
% Inputs:
% - S - source samples (may be multivariate as per generate3DObservations)
% - D - target samples (may be multivariate as per generate3DObservations)
% - Dpast - target past samples (multivariate, and embedded up to k previous samples)
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are loaded from loadProperties.m
%
% Outputs:
% - te - transfer entropy value (conditional MI from S (maybe plus RelSourcePos) to D given Dpast). If not requested, then
% the te (and an array of local values) is saved to properties.resultsFile
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
% -- STEP 0 : create java object
% check out http://lizier.me/joseph/software/jidt/javadocs/v1.3/
% for the description of all classes and methods
javaaddpath(properties.jidtJarLocation); % add JIDT path
% the following is the java class for conditional mutual information in JIDT
% transfer entropy is infact mutual information conditioning on the past
% of the destination
if (strcmp('kraskov', properties.estimator))
CMI_CLASS = 'infodynamics.measures.continuous.kraskov.ConditionalMutualInfoCalculatorMultiVariateKraskov1';
else
CMI_CLASS = 'infodynamics.measures.continuous.gaussian.ConditionalMutualInfoCalculatorMultiVariateGaussian';
end
TEcalculator = javaObject(CMI_CLASS); % creates a java object of the given class
% -- STEP 1 : set properties
TEcalculator.setProperty('k', num2str(properties.jidt.kNNs));
TEcalculator.setProperty('BIAS_CORRECTION', 'true'); % Used for Gaussian only
if (isfield(properties.jidt, 'dynamicCorrelationExclusion'))
% We'll ensure samples from the same target transition aren't included in nearest neighbour counts
% (it will exclude some others as well, but this only adds some small noise to the calculation)
TEcalculator.setProperty('DYN_CORR_EXCL', num2str(properties.jidt.dynamicCorrelationExclusion));
end
% -- STEP 2 : initialise
% here the parameters are the dimensionality of the series
% in this case taken directly from the number of columns in each variable
TEcalculator.initialise(size(S,2), size(D,2), size(Dpast,2));
% -- STEP 3 : add in observations
TEcalculator.setObservations(S, D, Dpast);
% -- STEP 4 : compute the local entropies
tranEntropy = TEcalculator.computeAverageLocalOfObservations(); % global (average) value
fprintf('Mean TE_%s (k=%d,tau=%d,lag=%d) = %.4f\n', ...
properties.estimator, properties.k, properties.tau, properties.lag, tranEntropy);
if (properties.teNumSurrogates > 0)
% Compute the (statistical significance via) null distribution empirically (e.g. with 100 permutations):
measDist = TEcalculator.computeSignificance(properties.teNumSurrogates);
fprintf('Null distribution: %.4f +/- %.4f std dev.; p(surrogate > measured)=%.5f from %d surrogates)\n', ...
measDist.getMeanOfDistribution(), measDist.getStdOfDistribution(), ...
measDist.pValue, properties.teNumSurrogates);
pValue = measDist.pValue;
meanOfSurrogates = measDist.getMeanOfDistribution();
stdOfSurrogates = measDist.getStdOfDistribution();
else
pValue = 1;
meanOfSurrogates = 0;
stdOfSurrogates = 0;
end
if (nargout >= 1)
% Supply the samples back to the caller
% (the caller is probably trying to optimise parameters at the moment)
% Nothing to do then actually...
else
% We're going to save the results instead
% First generate the local values to save as well:
localTranEntropy = TEcalculator.computeLocalOfPreviousObservations(); % local values
% save results
save(properties.resultsFile, 'tranEntropy', 'localTranEntropy', 'pValue', 'meanOfSurrogates', 'stdOfSurrogates', '-append');
fprintf('Transfer entropy saved in %s\n', properties.resultsFile);
end
end

457
demos/octave/FlockingAnalysis/generateObservations.m
View File

@ -1,457 +0,0 @@
function [D, Dpast, S, RelSourcePos, safeDynamicCorrelationExclusionSamples] = generateObservations(properties)
% This function generates the observations from which
% we can then compute information dynamics with JIDT.
% This will work for either 2D or 3D samples (as specified by the properties)
%
% Author: Emanuele Crosato, Joseph T. Lizier, Sooraj Sekhar, 2019
%
% Inputs:
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file.
%
% Outputs:
% - D - target samples (may be multivariate as per below)
% - Dpast - target past samples (multivariate, and embedded up to k previous samples)
% - S - source relative headings samples (may be multivariate as per below)
% - RelSourcePos - relative source position (may be multivariate as below)
% - safeDynamicCorrelationExclusionSamples - maximum number of samples with potential dynamic correlation, used for dynamic correlation exclusion externally for the calculations.
% If no outputs are requested, these are saved to properties.resultsFile
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
% Call utility to put filename lists in a common format
files = processFilenames(properties.files);
% initialize series for storing the samples:
S = []; % initialise source observations
D = []; % initialise destination observations
Dpast = []; % initialise destination past observations
% only S and D are necessary for lagged mutual information
% all S, D and Dpast are necessary for transfer entropy
fileTimeAndPair = []; % to save time and pair in the format [file_index, time_index, target_index, source_index]
RelSourcePos = []; % initialise relative source positions
sample = 1; % initialise the current sample number
maxSourceSamplesForATarget = 0; % Track the maximum number of in range sources for a given target sample
maxAcfTime = 1; % Track the maximum ACF time for targets
% Need to loop over fileIndex rather than for file = files (which doesn't work properly for cell array of length 1)
for fileIndex = 1:length(files)
dataFileName = files{fileIndex};
% load preprocessed data using the function specified in properties.loadScript
if (properties.data3d)
[x,y,z] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y,z])));
else
[x,y] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y])));
end
fprintf('Loading data in %s (%d missing values)\n', dataFileName, numMissing);
samplesBeforeThisFile = sample;
% Translate the raw positions into delta Positions
velX = x(2:end,:) - x(1:end-1,:);
velY = y(2:end,:) - y(1:end-1,:);
velXY = sqrt((velX .* velX) + (velY .* velY));
velALL = velXY;
if (properties.data3d)
velZ = z(2:end,:) - z(1:end-1,:);
velALL = sqrt((velXY .* velXY) + (velZ .* velZ));
end
% Translate velocities into headings:
if (properties.data3d)
% Spherical polars:
[headingXY,headingZ,speed] = cart2sph(velX,velY,velZ);
else
% Polar coordinates:
[headingXY,speed] = cart2pol(velX, velY);
end
% Manual way is (can verify these are same except for Nans on x,y and z):
% headingXY = atan(velY ./ velX) + (((velX < 0).*(velY>0)) .* pi) + ...
% (((velX < 0).*(velY<=0)) .* (- pi));
% But the above fails if x and y are *both* zero leaving Nan:
% velYOnVelX = velY ./ velX; % Don't replace all nans in headingXY as some mean missing values
% headingXY(isnan(velYOnVelX(:))) = 0;
% xyMagnitude = sqrt(velX.^2 + velY.^2);
% Since xyMagnitude can only be positive, we can take the straight atan to
% compute headingZ
% headingZ = atan(velZ ./ xyMagnitude);
% get number of fish and update time cycles
numCycles = size(velX,1);
numFish = size(velX,2);
% Calculate the autocorrelation decay time in the changes in headings and the speed:
acfHeadingsDecayTime = zeros(numFish, 1);
acfSpeedDecayTime = zeros(numFish, 1);
for f = 1:numFish
numLags = min(200, size(headingXY,1)-1);
try
% The following sometimes breaks I think if too many Nans
[acfHeading, lags] = autocorr(headingXY(2:end,f) - headingXY(1:end-1,f), 'NumLags', numLags);
catch ME
fprintf('ACF Heading broke for file %d fish %d: either you need to install the Econometrics toolbox (if properties.jidt.autoDynamicCorrelationExclusion = true) or there are too many NaNs in the file\n', fileIndex, f);
continue;
end
acfHeadingsDecayTime(f) = numLags;
for t = 1 : numLags
if (acfHeading(t) < exp(-1))
acfHeadingsDecayTime(f) = t;
break;
end
end
try
[acfSpeed, lags] = autocorr(speed(:,f), 'NumLags', numLags);
catch ME
fprintf('ACF Speed broke for file %d fish %d: either you need to install the Econometrics toolbox (if properties.jidt.autoDynamicCorrelationExclusion = true) or there are too many NaNs in the file\n', fileIndex, f);
continue;
end
acfSpeedDecayTime(f) = numLags;
for t = 1 : numLags
if (acfSpeed(t) < exp(-1))
acfSpeedDecayTime(f) = t;
break;
end
end
end
maxAcfHeadingsDecayTime = mean(acfHeadingsDecayTime); % max(acfHeadingsDecayTime);
maxAcfSpeedDecayTime = mean(acfSpeedDecayTime); % max(acfSpeedDecayTime);
if properties.headingcalc == true
acfTime = maxAcfHeadingsDecayTime;
end
if properties.speedcalc == true
acfTime = maxAcfSpeedDecayTime;
end
if properties.speedcalc == true && properties.headingcalc == true
acfTime = max(maxAcfHeadingsDecayTime, maxAcfSpeedDecayTime);
end
if (acfTime > maxAcfTime)
maxAcfTime = acfTime;
end
fprintf('ACF time of %d from heading %.2f and speed %.2f\n', acfTime, maxAcfHeadingsDecayTime, maxAcfSpeedDecayTime);
% Initialising destPastSample is only important in terms of ensuring it is a row vector.
% If we have 3D data, the vector will get padded out to the appropriate length with the first sample below.
destPastSampleHead = zeros(1, properties.k);
destPastSampleSpeed = zeros(1, properties.k);
startTime = max(1+properties.lag, (properties.k-1)*properties.tau + 3); % Adding 3: one for target, one for first target past, one for taking differences
for i = startTime : numCycles % cycle over time
timePointForSourceHeading = i-properties.lag; % This is indexed into velX and headingXY, hence no extra +1 !
timePointForSourcePosition = i+1-properties.lag; % This is indexed into x not velX, hence the extra +1 !
% Compute relative position of source (at time timePointForSourcePosition) to
% target either at this same time step or the current time at which it is updating.
% Note this position difference is relative to absolute Cartesian coordinates
% (we'll convert to relative to source heading later):
if (properties.sourceWrtSameDestTime)
destPositionTimePointRef = timePointForSourcePosition;
else
destPositionTimePointRef = i;
end
% cycle over fish pairs
for idxFD = 1 : numFish % Target/Destination
% check destination variable
if isnan(headingXY(i,idxFD)) || isnan(headingXY(i-1,idxFD))
continue;
end
if properties.data3d && (isnan(headingZ(i,idxFD)) || isnan(headingZ(i-1,idxFD)))
continue;
end
% check destination past vector
missingFound = false;
for h = 1 : properties.k
idx = i-1-(h-1)*properties.tau;
if isnan(headingXY(idx,idxFD)) || isnan(headingXY(idx-1,idxFD))
missingFound = true;
break;
end
if properties.data3d && (isnan(headingZ(idx,idxFD)) || isnan(headingZ(idx-1,idxFD)))
missingFound = true;
break;
end
end
if missingFound
continue;
end
% Postcondition: All destination variables are ok
% We will create **source** observation as source headings relative to target headings at appropriate time point:
if (properties.sourceWrtSameDestTime)
% Take reference dest heading at same time as source:
theta_FDXY_ref = headingXY(timePointForSourceHeading,idxFD);
else
% Take reference dest heading at prev time step:
theta_FDXY_ref = headingXY(i-1,idxFD);
end
if properties.data3d
if (properties.sourceWrtSameDestTime)
theta_FDZ_ref = headingZ(timePointForSourceHeading,idxFD);
else
theta_FDZ_ref = headingZ(i-1,idxFD);
end
end
% create **destination** observation as change in headings:
if properties.headingcalc == true
theta_FDXY_curr = headingXY(i,idxFD);
theta_FDXY_prev = headingXY(i-1,idxFD);
if properties.data3d
theta_FDZ_curr = headingZ(i,idxFD);
theta_FDZ_prev = headingZ(i-1,idxFD);
destSampleHead = [angleDifference(theta_FDXY_curr, theta_FDXY_prev), ...
angleDifference(theta_FDZ_curr, theta_FDZ_prev)];
else
destSampleHead = angleDifference(theta_FDXY_curr, theta_FDXY_prev);
end
destSample = destSampleHead;
end
if properties.speedcalc == true
destSampleSpeed = velALL(i,idxFD);
destSample = destSampleSpeed;
end
if properties.speedcalc == true && properties.headingcalc == true
destSample = [destSampleHead, destSampleSpeed];
end
% create **destination past** observation as changes in headings at each step:
% TODO: we could take differences to previous sample amongst the k rather than only
% one back from each sample: I'm not sure if this would be a more wholistic embedding or not
% (only makes a difference if tau>1)
if (properties.headingcalc == true)
DpastColIndexHead = 1;
for h = 1 : properties.k
idx = i-1-(h-1)*properties.tau;
theta_FDXY_curr = headingXY(idx,idxFD);
theta_FDXY_prev = headingXY(idx-1,idxFD);
destPastSampleHead(DpastColIndexHead) = angleDifference(theta_FDXY_curr, theta_FDXY_prev);
DpastColIndexHead = DpastColIndexHead + 1;
if properties.data3d
theta_FDZ_curr = headingZ(idx,idxFD);
theta_FDZ_prev = headingZ(idx-1,idxFD);
destPastSampleHead(DpastColIndexHead) = angleDifference(theta_FDZ_curr, theta_FDZ_prev);
DpastColIndexHead = DpastColIndexHead + 1;
end
end
destPastSample = destPastSampleHead;
end
if (properties.speedcalc == true)
DpastColIndexSpeed = 1;
for h = 1 : properties.k
idx = i - 1 - (h - 1) * properties.tau;
destPastSampleSpeed (DpastColIndexSpeed) = velALL(idx,idxFD);
DpastColIndexSpeed = DpastColIndexSpeed + 1;
end
destPastSample = destPastSampleSpeed;
end
if (properties.speedcalc == true && properties.headingcalc == true)
destPastSample = [destPastSampleHead, destPastSampleSpeed];
end
if (isfield(properties, 'destSamplesOnly'))
if (properties.destSamplesOnly)
% User has asked for [D,Dpast] samples only to be returned,
% so we can do these now (without looping over sources):
% Fill in the destination and destination next samples now from above:
D(sample, :) = destSample;
Dpast(sample, :) = destPastSample;
fileTimeAndPair(sample,:) = [fileIndex i idxFD nan];
% increment sample number
sample = sample + 1;
continue; % skip looping over the sources
end
end
numSourceSamplesForThisTarget = 0;
for idxFS = 1 : numFish % Source
% check not the same fish
if (idxFD == idxFS)
continue;
end
relXOfSource = x(timePointForSourcePosition,idxFS) - x(destPositionTimePointRef,idxFD);
relYOfSource = y(timePointForSourcePosition,idxFS) - y(destPositionTimePointRef,idxFD);
if (properties.data3d)
relZOfSource = z(timePointForSourcePosition,idxFS) - z(destPositionTimePointRef,idxFD);
[xyAbsoluteAngleOfSource,zAbsoluteAngleOfSource,distanceBetween] = ...
cart2sph(relXOfSource,relYOfSource,relZOfSource);
else
[xyAbsoluteAngleOfSource,distanceBetween] = ...
cart2pol(relXOfSource,relYOfSource);
end
% Manually: (verified this matches cart2sph):
% xyAbsoluteAngleOfSource = atan(relYOfSource ./ relXOfSource) + ...
% (((relXOfSource < 0).*(relYOfSource>0)) .* pi) + ...
% (((relXOfSource < 0).*(relYOfSource<=0)) .* (- pi));
% xyRelMagnitude = sqrt(relXOfSource.^2 + relYOfSource.^2);
% zAbsoluteAngleOfSource = atan(relZOfSource ./ xyRelMagnitude);
% distanceBetween = sqrt(relXOfSource.^2 + relYOfSource.^2 + ...
% relZOfSource.^2);
% check in range
if (distanceBetween > properties.pairRange)
continue;
end
% check source variable
if isnan(headingXY(timePointForSourceHeading,idxFD)) || isnan(headingXY(timePointForSourceHeading,idxFS))
continue;
end
if properties.data3d && (isnan(headingZ(timePointForSourceHeading,idxFD)) || isnan(headingZ(timePointForSourceHeading,idxFS)))
continue;
end
% Postcondition: There are no missing headings so we can generate an observation.
numSourceSamplesForThisTarget = numSourceSamplesForThisTarget + 1;
% Now compose the data that will be saved for this sample:
% fileTimeAndPair is [file_index, time_index, target_index, source_index]
fileTimeAndPair(sample,:) = [fileIndex i idxFD idxFS];
% create **source** observation as source headings relative to target headings at appropriate time point:
theta_FSXY_lag = headingXY(timePointForSourceHeading,idxFS);
% And convert the absolute angular positions into
% relative angular positions compared to the target's
% heading.
xyRelativeAngleOfSource = angleDifference(...
xyAbsoluteAngleOfSource, ...
theta_FDXY_ref);
if properties.data3d
theta_FSZ_lag = headingZ(timePointForSourceHeading,idxFS);
sourceSampleHeading = [angleDifference(theta_FDXY_ref, theta_FSXY_lag), ...
angleDifference(theta_FDZ_ref, theta_FSZ_lag)];
sourceSample = sourceSampleHeading;
% Elevation angle differences need to be in -pi/2,pi/2 range
zRelativeAngleOfSource = angleDifferencePiOn2(...
zAbsoluteAngleOfSource, ...
theta_FDZ_ref);
% Store these relative polar coordinates of source at timePointForSourcePosition
RelSourcePos(sample,:) = [distanceBetween, xyRelativeAngleOfSource, zRelativeAngleOfSource];
else
sourceSampleHeading = angleDifference(theta_FDXY_ref, theta_FSXY_lag);
sourceSample = sourceSampleHeading;
% Store these relative polar coordinates of source at timePointForSourcePosition
RelSourcePos(sample,:) = [distanceBetween, xyRelativeAngleOfSource];
end
if properties.speedcalc == true
sourceSampleSpeed = velALL(timePointForSourceHeading,idxFS);
sourceSample = sourceSampleSpeed;
end
if properties.speedcalc == true && properties.headingcalc == true
sourceSample = [sourceSampleHeading, sourceSampleSpeed];
end
if properties.includeSourcePositionInTransfer
S(sample,:) = [sourceSample, RelSourcePos(sample,:)];
else
S(sample,:) = sourceSample;
end
% Fill in the destination and destination next samples now from above:
D(sample, :) = destSample;
Dpast(sample, :) = destPastSample;
% increment sample number
sample = sample + 1;
end
if (numSourceSamplesForThisTarget > maxSourceSamplesForATarget)
maxSourceSamplesForATarget = numSourceSamplesForThisTarget;
end
end
% fprintf('Run time step %d\n', i);
end
fprintf(' added %d samples\n', sample - samplesBeforeThisFile);
end
if (sample - 1 == 0)
% We've added no samples
warning('No samples added for the given parameters!');
end
if (isfield(properties, 'destSamplesOnly'))
if (properties.destSamplesOnly)
safeDynamicCorrelationExclusionSamples = maxAcfTime * numFish;
else
safeDynamicCorrelationExclusionSamples = maxAcfTime * numFish * maxSourceSamplesForATarget;
end
end
if (nargout > 1)
% Supply the samples back to the caller
% (the caller is probably trying to optimise parameters at the moment)
% Nothing to do then actually...
else
% We're going to save the samples instead
% display to check
fprintf('Displaying first 5 samples for source, target, target past and fileTimeAndPair as a check:');
disp(S(1:5,:));
disp(D(1:5,:));
disp(Dpast(1:5,:));
disp(fileTimeAndPair(1:5,:));
% input(prompt);
% save series and properties
save(properties.resultsFile, 'S', 'D', 'Dpast', 'files', 'fileTimeAndPair', 'RelSourcePos', ...
'maxSourceSamplesForATarget', 'safeDynamicCorrelationExclusionSamples', 'properties');
fprintf('Series saved in %s (%d samples in total)\n', properties.resultsFile, sample - 1);
end
end
% End function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for computing the difference between two angles.
% the differrence must be between pi and -pi
function [diff] = angleDifference(angleA, angleB)
diff = angleA - angleB; % subtract angles
if abs(diff) > pi % if absolute value is larger than pi
% replace with the complementary angle and switch sign
diff = (2*pi - abs(diff) ) * (-sign(diff));
end
end
% function for computing the difference between two angles in [-pi/2,pi/2].
% the differrence must be returned between pi/2 and -pi/2.
% This is used for differences in elevation angles
function [diff] = angleDifferencePiOn2(angleA, angleB)
diff = angleA - angleB; % subtract angles
% Pre-condition: differences between angles which were in range of
% [-pi/2,pi/2] can only be in range [-pi,pi]
if (diff > pi/2)
diff = pi/2 - (diff - pi/2);
elseif (diff < -pi/2)
diff = -pi/2 + (-pi/2 - diff);
end
end

75
demos/octave/FlockingAnalysis/loadBasic2d.m
View File

@ -1,75 +0,0 @@
function [posX,posY] = loadBasic2d(dataFileName, properties)
% This script loads the raw data from a .txt file,
% preprocesses the data and save it as a .mat file.
% The txt data is assumed to have early columns with other data
% (e.g. a timestamp in column 1), with the fish coordinates
% starting from properties.loadBasic2d.startColumn (defaults to 2)
% with fish 1's X and Y coordinates in the first of those columns (defaults to 2 and 3),
% then fish 2's X and Y coordinates in the next of those columns (defaults to 4 and 5),
% and so on.
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
% prompt = 'Press a key to continue';
%%% LOAD RAW DATA %%%
%fprintf('Loading data in %s\n', dataFileName); % print to check
% input(prompt);
% read the .txt as a Matlab table (assuming .txt is tab-separated)
data = load(dataFileName);
%fprintf('Size of data is %d - %d\n', ... % print table's size - %d is for int
% size(data,1), size(data,2)); % see also %f (real) and %s (string)
% disp(data(1:5,:)); % display first 10 rows to check
% input(prompt);
%%% FORMAT DATA IN A MORE CONVENIENT WAY %%%
if (nargin == 1) || (~isfield(properties, 'loadBasic2d'))
% Assume the first column is a datestamp
properties.loadBasic2d.startColumn = 2;
end
startCol = properties.loadBasic2d.startColumn;
colsToSkip = startCol - 1;
numFish = (size(data,2)-colsToSkip) ./ 2; % number of fish (we know it from the raw data)
numCycles = size(data,1); % number of time steps (we also know it)
fprintf('Number of fish %d and cycles %d\n', numFish, numCycles);
% prepare variables for x and y position
% as a table [numFish x numCycles]
posX = nan(numCycles,numFish);
posY = nan(numCycles,numFish);
% fill the position tables
for f = 1 : numFish % cycle over all fish
% copy into new variables
posX(:,f) = data(:,startCol+(f-1)*2);
posY(:,f) = data(:,startCol+1+(f-1)*2);
end
% display to check
% disp(size(posX));
% disp(posX(1:5,:));
% disp(size(posY));
% disp(posY(1:5,:));
% input(prompt);

40
demos/octave/FlockingAnalysis/loadseparatexy.m
View File

@ -1,40 +0,0 @@
function [x,y] = loadseparatexy(filename, properties)
% loadseparatexy loads 2D fish data from 2 separate txt files (one for x, one for y)
% where in each file time increases down the
% rows and then across the columns we have position columns for each
% fish in turn, i.e. in position x file:
% <fish1x>, <fish2x>, <fish3x>, etc
%
% Inputs:
% - filename - the name template of the file to load, with %s where 'x' and 'y' should be filled in
% - properties (not required) - properties object (may be required for other file loaders)
% Outputs:
% - x - 2D array, each row contains x position for each fish (in columns)
% - y - as per x
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
xfilename = sprintf(filename, 'x');
x = load(xfilename);
yfilename = sprintf(filename, 'y');
y = load(yfilename);
end

43
demos/octave/FlockingAnalysis/loadseparatexyz.m
View File

@ -1,43 +0,0 @@
function [x,y,z] = loadseparatexy(filename, properties)
% loadseparatexy loads 3D fish data from 3 separate txt files (one for x, one for y, one for z)
% where in each file time increases down the
% rows and then across the columns we have position columns for each
% fish in turn, i.e. in position x file:
% <fish1x>, <fish2x>, <fish3x>, etc
%
% Inputs:
% - filename - the name template of the file to load, with %s where 'x', 'y' and 'z' should be filled in
% - properties (not required) - properties object (may be required for other file loaders)
% Outputs:
% - x - 2D array, each row contains x position for each fish (in columns)
% - y - as per x
% - z - as per x
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
xfilename = sprintf(filename, 'x');
x = load(xfilename);
yfilename = sprintf(filename, 'y');
y = load(yfilename);
zfilename = sprintf(filename, 'z');
z = load(zfilename);
end

43
demos/octave/FlockingAnalysis/loadxls2d.m
View File

@ -1,43 +0,0 @@
function [x,y] = loadxls2d(filename, properties)
% loadxls2d loads 2D fish data from an xls file where time increases down the
% rows and then across the columns we have 2 x,y position columns for each
% fish in turn, i.e.:
% <timestamp1>, <fish1x>, <fish1y>, <fish2x>, <fish2y>, etc
% <timestamp2>, <fish1x>, <fish1y>, <fish2x>, <fish2y>, etc
%
% Inputs:
% - filename - the name of the file to load
% - properties (not required) - properties object (may be required for other file loaders)
% Outputs:
% - x - 2D array, each row contains x position for each fish (in columns)
% - y - as per x
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
[data,txt,raw] = xlsread(filename);
% Make sure we preprocess data to have x,y, as 2D arrays of x(timeStep, fishID):
% timeSteps = size(data,1); % number of rows; not required
xFishIndex = 2 : 2 : size(data,2); % which indices are the x values for different fish
x = data(:,xFishIndex);
y = data(:,xFishIndex+1);
end

45
demos/octave/FlockingAnalysis/loadxls3d.m
View File

@ -1,45 +0,0 @@
function [x,y,z] = loadxls3d(filename, properties)
% loadxls3d loads 3D fish data from an xls file where time increases down the
% rows and then across the columns we have 3 x,y,z position columns for each
% fish in turn, i.e.:
% <timestamp1>, <fish1x>, <fish1y>, <fish1z>, <fish2x>, <fish2y>, <fish2z>, etc
% <timestamp2>, <fish1x>, <fish1y>, <fish1z>, <fish2x>, <fish2y>, <fish2z>, etc
%
% Inputs:
% - filename - the name of the file to load
% - properties (not required) - properties object (may be required for other file loaders)
% Outputs:
% - x - 2D array, each row contains x position for each fish (in columns)
% - y - as per x
% - z - as per x
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
[data,txt,raw] = xlsread(filename);
% Make sure we preprocess data to have x,y, and z as 2D arrays of x(timeStep, fishID):
% timeSteps = size(data,1); % number of rows; not required
xFishIndex = 2 : 3 : size(data,2); % which indices are the x values for different fish
x = data(:,xFishIndex);
y = data(:,xFishIndex+1);
z = data(:,xFishIndex+2);
end

213
demos/octave/FlockingAnalysis/plotLocalTEs.m
View File

@ -1,213 +0,0 @@
function plotLocalTEs(properties)
% Plot the local TEs to show where the information transfer hotspots are from target fish relative to each source
%
% Author: Joseph T. Lizier, 2019
%
% Inputs:
% - properties - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are loaded from loadProperties.m
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
if (nargin < 1)
fprintf('No properties object supplied, attempting to load properties via a loadProperties script ...');
% By default, just try to load properties locally
if (exist('loadProperties') == 2)
% there is a loadProperties script
loadProperties;
else
% there is not a loadProperties script
error('No properties object supplied, and no loadProperties script found.');
end
end
load(properties.resultsFile);
% Loads:
% S -- source samples
% D -- target samples
% Dpast -- target past samples (embedded up to k previous samples)
% files -- cell array of file names that we took samples from
% fileTimeAndPair -- each row holds file index, time index, target index, source index
% RelSourcePos -- each row holds distance between the pair for this sample,
% their xyRelativeAngleOfSource, and zRelativeAngleOfSource
% lag -- source-target lag that is in use
% k -- embedding length for target that is in use
% tau -- embedding delay for target that is in use
% pairRange -- range within which we've pulled source-target interactions
% tranEntropy -- average transfer entropy
% localTranEntropy -- local TE for each sample
fprintf('%d samples in total for %d fish\n', length(S), length(unique(fileTimeAndPair(:,3))));
relDistance = RelSourcePos(:,1);
relTheta = RelSourcePos(:,2);
% I think this is giving us the right conversions:
if (properties.data3d)
relPhi = RelSourcePos(:,3);
distXY = relDistance .* cos(relPhi);
distZ = relDistance .* sin(relPhi);
else
distXY = relDistance;
end
distInFront = distXY .* cos(relTheta); % X coordinate
distToLeft = distXY .* sin(relTheta); % Y coordinate
% Plot where all the raw positions are:
% Will need to turn this off when we have too many
figure()
% polar(relTheta, distXY, '.r'); % This is equivalent to below:
scatter(distInFront, distToLeft, 2, localTranEntropy);
title('Relative position of source in XY plane for target heading, coloured for TE');
colorbar;
% Plot the density of samples
makePolarBinnedPlot(relTheta, distXY, ones(length(relTheta), 1), 12, 10, true, false);
title('Density of samples in each bin (r_{XY},\theta)');
% Plot the TE in XY plane
makePolarBinnedPlot(relTheta, distXY, localTranEntropy, 12, 10, true, true);
title('Average TE in each bin (r_{XY},\theta)');
% Plot raw positions in phi-z:
% figure()
% polar(relPhi, relDistance, '.r');
% title('Relative position of source in Z-phi plane for target heading');
if (properties.data3d)
% Plot the density of samples in distance-phi plane
makePolarBinnedPlot(relPhi, relDistance, ones(length(relTheta), 1), 12, 10, true, false);
title('Density of samples in each bin (r, \phi)');
% Plot the TE in distance-phi plane
makePolarBinnedPlot(relPhi, relDistance, localTranEntropy, 12, 10, true, true);
title('Average TE in each bin (r, \phi)');
xlabel('r_{XY} [mm]');
ylabel('z [mm]');
end
end
% Inputs:
% - thetas - angles for each sample
% - radii - radius for each sample
% - numAngleBins - how many bins to make across 2*pi
% - numRadialBins - how many bins to make up to the maximum radii
% - useMaxEntBinning - whether to make bins with approx same numbers of points (true)
% or same size (false)
% - plotMean - if true (default) plot the mean within each bin, else plot the total (divded by area)
% The latter is used for densities for example
function makePolarBinnedPlot(thetas, radii, valuesToPlot, numAngleBins, numRadialBins, useMaxEntBinning, plotMean)
if (nargin < 6)
useMaxEntBinning = false;
end
if (nargin < 7)
plotMean = true;
end
if ((min(thetas) < -pi/2) || (max(thetas) > pi/2))
% We're using full angular range -pi : pi
minAngle = -pi;
maxAngle = pi;
extraBinForPlotWrap = true;
else
% We're only using -pi/2:pi/2
minAngle = -pi/2;
maxAngle = pi/2;
extraBinForPlotWrap = false;
end
angleStep = (maxAngle - minAngle) / numAngleBins;
radiusStep = max(radii) / numRadialBins;
% Simple way to do the binning for even bins:
% binnedAngles = floor(thetas ./ angleStep); % Gives the discrete bin for the angle
% binnedRadii = floor(radii ./ radiusStep); % Gives the discrete bin for the radius
% More general, and allowing bins to spread with points:
if (useMaxEntBinning)
% Space the bins for roughly same
% numbers of points (when examined marginally):
sortedAngles = sort(thetas);
binAngleEdges = [minAngle; sortedAngles(floor((1:(numAngleBins-1)).*length(sortedAngles)./numAngleBins)); maxAngle]';
sortedRadii = sort(radii);
binRadiusEdges = [0; sortedRadii(floor((1:(numRadialBins-1)).*length(sortedRadii)./numRadialBins)); max(radii)]';
else
% Space the bins equally
binAngleEdges = minAngle:angleStep:maxAngle;
binRadiusEdges = 0:radiusStep:max(radii);
end
[angleHistCounts,binnedAngles] = histc(thetas, binAngleEdges);
[radiiHistCounts,binnedRadii] = histc(radii, binRadiusEdges);
angleBinValues = 1:numAngleBins; % unique(binnedAngles);
radiusBinValues = 1:numRadialBins; % unique(binnedRadii);
if (extraBinForPlotWrap)
valuesForEachBin = zeros(length(angleBinValues) + 1, length(radiusBinValues));
else
valuesForEachBin = zeros(length(angleBinValues), length(radiusBinValues));
end
numberOfSamples = 0;
for aIndex = 1 : length(angleBinValues)
for rIndex = 1 : length(radiusBinValues)
indicesForThisBin = find((binnedAngles == angleBinValues(aIndex)) & (binnedRadii == radiusBinValues(rIndex)));
if (plotMean)
valueForThisBin = mean(valuesToPlot(indicesForThisBin));
else
% Plot a density: compute total then divide by area.
valueForThisBin = sum(valuesToPlot(indicesForThisBin));
areaOfBin = pi .* (binRadiusEdges(rIndex+1).^2 - binRadiusEdges(rIndex).^2) .* ...
mod(abs(binAngleEdges(aIndex+1) - binAngleEdges(aIndex)), 2.*pi) ./ (2.*pi);
valueForThisBin = valueForThisBin ./ areaOfBin;
end
if (length(indicesForThisBin) == 0)
valueForThisBin = 0;
end
% fprintf('Mean value for r=%.1f+,theta=%.3f+ is %.3f (from %d samples)\n', binRadiusEdges(rIndex), ...
% binAngleEdges(aIndex), valueForThisBin, length(indicesForThisBin));
numberOfSamples = numberOfSamples + length(indicesForThisBin);
valuesForEachBin(aIndex, rIndex) = valueForThisBin;
end
end
% Now convert these so we can plot them:
fprintf('Found %d in total in the bins\n', numberOfSamples);
if (extraBinForPlotWrap)
% And add for first angle again to complete the plot
valuesForEachBin(end,:) = valuesForEachBin(1,:);
[THETA,RR] = meshgrid([(binAngleEdges(1:end-1)+binAngleEdges(2:end))./2, (binAngleEdges(1)+binAngleEdges(2))./2], ...
(binRadiusEdges(1:end-1)+binRadiusEdges(2:end))./2);
else
[THETA,RR] = meshgrid([(binAngleEdges(1:end-1)+binAngleEdges(2:end))./2], ...
(binRadiusEdges(1:end-1)+binRadiusEdges(2:end))./2);
end
[A,B] = pol2cart(THETA,RR);
figure();
% Old way which pinned TE values on the vertices of polygons (looks yuck)
% surf(A,B,valuesForEachBin','edgecolor','none')
% New way, smoothed visualisation:
plot = pcolor(A,B,valuesForEachBin');
plot.FaceColor = 'interp';
set(plot, 'EdgeColor', 'none');
xlabel('x [mm]');
ylabel('y [mm]');
colorbar;
view(0,90)
end

41
demos/octave/FlockingAnalysis/processFilenames.m
View File

@ -1,41 +0,0 @@
function fileCellArray = processFilenames(fileList)
%
% Turns the fileList from the properties file (usually properties.files) into a cell array of file names. The fileList can be either:
% a. a cell array of file names, e.g.: {'file1.xlsx', 'file2.xlsx'}
% b. a call to ls or ls with an argument, e.g. ls('*.xlsx')
% c. a space or tab separated character row vector of file names
% d. a character matrix of filenames (each filename on a separate row)
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2022, Joseph T. Lizier et al.
%%
%% 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/>.
%%
if (iscell(fileList))
% We're done already:
fileCellArray = fileList;
elseif (isvector(fileList))
% We have a row vector of space/tab separate filenames:
fileCellArray = strsplit(strtrim(fileList)); % extra strtrim to remove trailing \n's
elseif (ismatrix(fileList))
fileCellArray = {};
for r = 1 : size(fileList, 1)
fileCellArray{r} = strtrim(fileList(r,:));
end
else
error('fileList appears to be of an incorrect format\n');
end
end

220
demos/octave/FlockingAnalysis/runAnalysis.m
View File

@ -1,220 +0,0 @@
function runAnalysis(properties)
% This high-level function generates the local transfer entropy results, first optimising
% parameters (i.e. embedding length and delay, and source-target lag), then
% storing local transfer entropy values for the optimised parameters.
%
% Author: Joseph T. Lizier, 2019
%
% Inputs:
% - properties (required) - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are attempted to be loaded from loadProperties.m
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
if (nargin < 1)
fprintf('No properties object supplied, attempting to load properties via a loadProperties script ...');
% By default, just try to load properties locally
if (exist('loadProperties') == 2)
% there is a loadProperties script so attempt to run it to load a properties object
loadProperties;
else
% there is not a loadProperties script
error('No properties object supplied, and no loadProperties script found.');
end
else
% A properties argument was supplied
if (ischar(properties))
% We're assuming it was the name of a properties file
if ((length(properties) > 2) && (strcmp(properties(end-1:end), '.m')))
% Remove the '.m':
properties(end-1:end) = [];
end
if (exist(properties) == 2)
% attempt to run the properties .m file: (after making sure the properties variable is cleared; not necessary but is clean)
propertiesFile = properties;
clear properties;
eval(propertiesFile);
else
error('%s is not an .m file we can find that can be used to load a properties object', properties);
end
% else
% We assume it was the properties object.
end
end
% Step 1: Auto-embed if required:
% Set the lag to 1 as a dummy if we are optimising over that later as well:
if (~isfield(properties, 'lag'))
properties.lag = 1;
end
if (isfield(properties, 'kRange') || isfield(properties, 'tauRange'))
% If any one of these two ranges weren't supplied, set the range variables
% to the value of the corresponding non-range variable:
if (~isfield(properties, 'kRange'))
properties.kRange = properties.k;
end
if (~isfield(properties, 'tauRange'))
properties.tauRange = properties.tau;
end
% Ask generateObservations to only return the target samples for the AIS calculation
properties.destSamplesOnly = true;
% Optimise k and tau
maxAIS = -inf;
maxAISk = properties.kRange(1);
maxAIStau = properties.tauRange(1);
aisForKAndTau = zeros(length(properties.kRange), length(properties.tauRange));
kIndex = 0;
for k = properties.kRange
kIndex = kIndex + 1;
properties.k = k;
minTau = min(properties.tauRange);
tauIndex = 0;
for tau = properties.tauRange
tauIndex = tauIndex + 1;
if ((k == 1) && (tau > minTau))
% We only need compute k=1 for a single tau
continue;
end
properties.tau = tau;
% Generate the observations for k,tau:
[D, Dpast, ~, ~, safeDynamicCorrelationExclusionSamples] = generateObservations(properties);
if (isempty(D))
% There were no samples found for the given parameters, presumably k etc are too long
continue;
end
% Check if we're turning on dynamic correlation exclusion:
if (isfield(properties.jidt, 'autoDynamicCorrelationExclusion'))
properties.jidt.dynamicCorrelationExclusion = safeDynamicCorrelationExclusionSamples;
end
% Compute the AIS:
ais = computeAIS(D, Dpast, properties);
if (ais > maxAIS)
maxAIS = ais;
maxAISk = k;
maxAIStau = tau;
end
aisForKAndTau(kIndex, tauIndex) = ais;
end
end
% Optimisation is complete:
properties.k = maxAISk;
properties.tau = maxAIStau;
properties.ais = maxAIS;
properties.destSamplesOnly = false;
fprintf('*** Optmised k=%d and tau=%d (giving AIS=%.4f - see above for null distribution for these parameters)\n', properties.k, ...
properties.tau, properties.ais);
else
% Hard coded embedding parameters: compute the AIS to be saved anyway :
% Ask generateObservations to only return the target samples for the AIS calculation
properties.destSamplesOnly = true;
[D, Dpast, ~, ~, safeDynamicCorrelationExclusionSamples] = generateObservations(properties);
if (isempty(D))
% There were no samples found for the given parameters, presumably k etc are too long
error('No embeddings possible for the given hard coded k and tau\n');
end
% Check if we're turning on dynamic correlation exclusion:
if (isfield(properties.jidt, 'autoDynamicCorrelationExclusion'))
properties.jidt.dynamicCorrelationExclusion = safeDynamicCorrelationExclusionSamples;
end
% Compute the AIS:
ais = computeAIS(D, Dpast, properties);
properties.ais = ais;
properties.destSamplesOnly = false;
fprintf('*** Hard-coded values for k=%d and tau=%d to be used (giving AIS=%.4f)\n', properties.k, ...
properties.tau, properties.ais);
end
% Check whether to continue to TE calculations:
if (~isfield(properties, 'computeAISOnly'))
properties.computeAISOnly = false;
end
if (properties.computeAISOnly)
fprintf('Finishing after AIS calculation only, as requested.\n');
return;
end
% Step 2: automatically select the correct lag if required:
teNumSurrogates = properties.teNumSurrogates; % Store this for later, turn it off now
properties.teNumSurrogates = 0; % No need to run any surrogates during parameter fitting
if (isfield(properties, 'lagRange'))
% Caller asks us to maximise the TE over a given range:
maxTE = -inf;
maxTElag = properties.lagRange(1);
teForLag = zeros(length(properties.lagRange), 1);
lagIndex = 0;
for lag = properties.lagRange
lagIndex = lagIndex + 1;
properties.lag = lag;
% Generate the observations for k,tau,lag:
[D, Dpast, S, ~, safeDynamicCorrelationExclusionSamples] = generateObservations(properties);
if (isempty(S))
% There were no samples found for the given parameters, presumably k etc are too long
continue;
end
% Check if we're turning on dynamic correlation exclusion:
if (isfield(properties.jidt, 'autoDynamicCorrelationExclusion'))
properties.jidt.dynamicCorrelationExclusion = safeDynamicCorrelationExclusionSamples;
end
% Compute the TE:
te = computeTE(S, D, Dpast, properties);
if (te > maxTE)
maxTE = te;
maxTElag = lag;
end
teForLag(lagIndex) = te;
end
% Optimisation is complete:
properties.lag = maxTElag;
properties.tranEntropy = maxTE;
fprintf('*** Optmised lag=%d (giving TE=%.4f)\n', properties.lag, ...
properties.tranEntropy);
else
fprintf('*** Hard-coded value for lag=%d to be used\n', properties.lag);
end
% 3. Compute TE with the correct parameters
% Now, once again pre-process the positional data into velocities, this time
% saving them into the results file (by not requesting [S,D,Dpast] outputs):
generateObservations(properties);
% And load these samples (S, D, Dpast, properties, maxSourceSamplesForATarget, safeDynamicCorrelationExclusionSamples, etc) in from the saved file:
load(properties.resultsFile);
properties.teNumSurrogates = teNumSurrogates; % Allow surrogates to be computed for this final run with correct parameters
% And compute the TE again for the optimal parameters, this time
% saving the files:
% Turn on dynamic correlation exclusion if required:
if (isfield(properties.jidt, 'autoDynamicCorrelationExclusion'))
if (properties.jidt.autoDynamicCorrelationExclusion)
properties.jidt.dynamicCorrelationExclusion = safeDynamicCorrelationExclusionSamples; % safeDynamicCorrelationExclusionSamples was loaded from the results file
else
properties.jidt.dynamicCorrelationExclusion = 0; % no dynamic correlation exclusion
end
end
% Compute TE with no output arguments so that results are saved
computeTE(S, D, Dpast, properties);
save(properties.resultsFile, 'ais', '-append'); % Add the AIS into the results file as well
if (isfield(properties, 'kRange') || isfield(properties, 'tauRange'))
save(properties.resultsFile, 'aisForKAndTau', '-append'); % Add the AISs computed in auto-embedding
end
if (isfield(properties, 'lagRange'))
save(properties.resultsFile, 'teForLag', '-append'); % Add the TEs computed in optimising the source-target lag
end
end

151
demos/octave/FlockingAnalysis/watchLocalTEs.m
View File

@ -1,151 +0,0 @@
function watchLocalTEs(dataFileIndex, properties, refreshRate, zoomIn, plotTEIn)
% Plots a movie of local TEs in the swarm from the given input file
%
% Author: Joseph T. Lizier, 2020
%
% Inputs:
% - dataFileIndex - the index of the file to load, from the list of files listed in the properties file. Default 1.
% - properties - object with properties for the calculations,
% with sub-members as specificied in the loadProperties.m file. If not supplied
% the properties are loaded from loadProperties.m
% - refreshRate - how often to change the plotted data (default 0.1 sec)
% - zoomIn - whether to zoom into the individuals (default), or have plot take in whole field
% - plotTEIn - whether to plot the average TE into a target (true, default) or average TE out from a source (false)
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
if (nargin < 1)
dataFileIndex = 1;
end
if (nargin < 2)
fprintf('No properties object supplied, attempting to load properties via a loadProperties script ...\n');
% By default, just try to load properties locally
if (exist('loadProperties') == 2)
% there is a loadProperties script
loadProperties;
else
% there is not a loadProperties script
error('No properties object supplied, and no loadProperties script found.');
end
end
load(properties.resultsFile);
% Loads:
% S -- source samples
% D -- target samples
% Dpast -- target past samples (embedded up to k previous samples)
% files -- cell array of file names that we took samples from
% fileTimeAndPair -- each row holds file index, time index, target index, source index
% RelSourcePos -- each row holds distance between the pair for this sample,
% their xyRelativeAngleOfSource, and zRelativeAngleOfSource
% lag -- source-target lag that is in use
% k -- embedding length for target that is in use
% tau -- embedding delay for target that is in use
% pairRange -- range within which we've pulled source-target interactions
% tranEntropy -- average transfer entropy
% localTranEntropy -- local TE for each sample
fprintf('%d samples in total for %d fish\n', length(S), length(unique(fileTimeAndPair(:,3))));
if (nargin < 3)
refreshRate = 0.1;
end
% Call utility to put filename lists in a common format
files = processFilenames(properties.files);
dataFileName = files{dataFileIndex};
% load preprocessed data using the function specified in properties.loadScript
if (properties.data3d)
[x,y,z] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y,z])));
maxZ = max(z(:));
minZ = min(z(:));
else
[x,y] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y])));
end
maxX = max(x(:));
minX = min(x(:));
maxY = max(y(:));
minY = min(y(:));
if (nargin < 4)
zoomIn = true;
end
figure()
if (~zoomIn)
fprintf('Focussing...\n');
if (properties.data3d)
axis([minX, maxX, minY, maxY, minZ, maxZ]);
else
axis([minX, maxX, minY, maxY]);
end
% Need to set this so that the axes don't keep updating
end
fprintf('%d missing values\n', numMissing);
if (nargin < 5)
plotTEIn = true;
end
% Work out the range of TEs for this data file:
teMin = min(localTranEntropy(find(fileTimeAndPair(:,1) == dataFileIndex)));
teMax = max(localTranEntropy(find(fileTimeAndPair(:,1) == dataFileIndex)));
% Loop over all the time steps in this data file
numFish = size(x,2);
cb = colorbar;
xlabel('x');
ylabel('y');
cb.Label.String = 'Av local TE';
% caxis([teMin teMax]); % These are likely too extreme for the averages
for t = 1:size(x, 1)
set(gca,'NextPlot','replacechildren') ;
% Now loop over all fish as either source or target:
averageTEs = zeros(1,numFish);
for f = 1:numFish
% Find which interactions involve TE into or out from this fish at this time step
if (plotTEIn)
% For TE in, match the target:
rowIDs = find((fileTimeAndPair(:,1) == dataFileIndex) & (fileTimeAndPair(:,2) == t) & ...
(fileTimeAndPair(:,3) == f));
else
% For TE out, match the source:
rowIDs = find((fileTimeAndPair(:,1) == dataFileIndex) & (fileTimeAndPair(:,2) == t) & ...
(fileTimeAndPair(:,4) == f));
end
% Now average the TE into or out of this fish:
averageTEs(f) = mean(localTranEntropy(rowIDs));
end
if (properties.data3d)
scatter3(x(t,:), y(t,:), z(t,:), 5, averageTEs);
else
scatter(x(t,:), y(t,:), 5, averageTEs)
end
if (plotTEIn)
title(sprintf('Average TEs into each individual at time %d, coloured for TE', t));
else
title(sprintf('Average TEs out from each individual at time %d, coloured for TE', t));
end
pause(refreshRate)
end
end

74
demos/octave/FlockingAnalysis/watchRawMovement.m
View File

@ -1,74 +0,0 @@
function watchRawMovement(dataFileName, properties, refreshRate, zoomIn)
% Plots a movie of raw movement data from the given input file
%
% Author: Joseph T. Lizier, 2020
%
% Inputs:
% - dataFileName - the name of the file to load
% - properties - properties object
% - refreshRate - how often to change the plotted data (default 0.1 sec)
% - zoomIn - whether to zoom into the individuals (default), or have plot take in whole field
%%
%% Java Information Dynamics Toolkit (JIDT)
%% Copyright (C) 2019, Joseph T. Lizier et al.
%%
%% 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/>.
%%
if (nargin < 3)
refreshRate = 0.1;
end
if (nargin < 4)
zoomIn = true;
end
% load preprocessed data using the function specified in properties.loadScript
if (properties.data3d)
[x,y,z] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y,z])));
maxZ = max(z(:));
minZ = min(z(:));
else
[x,y] = feval(properties.loadScript, dataFileName, properties);
numMissing = sum(sum(isnan([x,y])));
end
maxX = max(x(:));
minX = min(x(:));
maxY = max(y(:));
minY = min(y(:));
if (~zoomIn)
fprintf('Focussing...\n');
if (properties.data3d)
axis([minX, maxX, minY, maxY, minZ, maxZ]);
else
axis([minX, maxX, minY, maxY]);
end
% Need to set this so that the axes don't keep updating
end
fprintf('Loading data in %s (%d missing values)\n', dataFileName, numMissing);
for t = 1:size(x, 1)
set(gca,'NextPlot','replacechildren') ;
if (properties.data3d)
plot3(x(t,:), y(t,:), z(t,:), 'x')
else
plot(x(t,:), y(t,:), 'x')
end
pause(refreshRate)
end
end

15
demos/octave/javaMatrixToOctave.m Normal file → Executable file
View File

@ -47,12 +47,7 @@ function octaveMatrix = javaMatrixToOctave(javaMatrix, startRow, startCol, numRo
% Convert whole matrix first:
tmp = javaObject('org.octave.Matrix', javaMatrix);
% Make sure tmp.ident() is converted to native octave:
if (exist('java_matrix_autoconversion') > 0)
oldFlag = java_matrix_autoconversion (1);
else
% Must be old octave version:
oldFlag = java_convert_matrix (1);
end
oldFlag = java_convert_matrix (1);
converted = false;
unwind_protect
octaveMatrix = tmp.ident(tmp);
@ -60,12 +55,7 @@ function octaveMatrix = javaMatrixToOctave(javaMatrix, startRow, startCol, numRo
unwind_protect_cleanup
% restore to non-default conversion, otherwise we get
% bad errors on other calls
if (exist('java_matrix_autoconversion') > 0)
java_matrix_autoconversion(oldFlag);
else
% Must be old octave version:
java_convert_matrix(oldFlag);
end
java_convert_matrix(oldFlag);
end_unwind_protect
if (converted)
if (nargin >= 2)
@ -77,7 +67,6 @@ function octaveMatrix = javaMatrixToOctave(javaMatrix, startRow, startCol, numRo
else
% Else we're in matlab, in which case the native java type can be handled, so return it directly:
octaveMatrix = javaMatrix;
return;
end
% Else, we encountered an error in the octave resizing, so fall through to element by element conversion:

7
demos/octave/octaveToJavaIntArray.m
View File

@ -32,12 +32,7 @@ function jIntArray = octaveToJavaIntArray(octaveArray)
if (length(octaveArray) > 1)
% Do this the normal way
tmp = javaObject('infodynamics.utils.OctaveMatrix');
try
tmp.loadIntData(octaveArray,[1, length(octaveArray)]);
catch
% Most likely error here is that octaveArray is interpreted as booleans, so try loading as booleans:
tmp.loadBooleanAsIntData(octaveArray,[1, length(octaveArray)]);
end
tmp.loadIntData(octaveArray,[1, length(octaveArray)]);
jIntArray = tmp.asIntVector();
else
% For length 1 arrays, we need to perform a hack here or else

7
demos/octave/octaveToJavaIntMatrix.m
View File

@ -32,12 +32,7 @@ function jIntMatrix = octaveToJavaIntMatrix(octaveMatrix)
if ((rows(octaveMatrix)*columns(octaveMatrix)) > 1)
% Do this the normal way
tmp = javaObject('infodynamics.utils.OctaveMatrix');
try
tmp.loadIntData(reshape(octaveMatrix,1,rows(octaveMatrix)*columns(octaveMatrix)),[rows(octaveMatrix), columns(octaveMatrix)]);
catch
% Most likely error here is that octaveMatrix is interpreted as booleans, so try loading as booleans:
tmp.loadBooleanAsIntData(reshape(octaveMatrix,1,rows(octaveMatrix)*columns(octaveMatrix)),[rows(octaveMatrix), columns(octaveMatrix)]);
end
tmp.loadIntData(reshape(octaveMatrix,1,rows(octaveMatrix)*columns(octaveMatrix)),[rows(octaveMatrix), columns(octaveMatrix)]);
jIntMatrix = tmp.asIntMatrix();
else
% For length 1 arrays, we need to perform a hack here or else

259
demos/python/EffectiveNetworkInference/spiking/net_inf.py
View File

@ -1,259 +0,0 @@
# Argument order: network_type_name num_spikes sim_number target_index
from jpype import *
import random
import math
import os
import numpy as np
import pickle
import copy
import sys
# net_type_name is useful if you are iterating over multiple files with different network types.
# Looking at the definition of SPIKES_FILE_NAME and OUTPUT_FILE_PREFIX will imply what the purpose of
# these command line arguments is.
net_type_name = sys.argv[1]
num_spikes_string = sys.argv[2]
repeat_num_string = sys.argv[3]
target_index_string = sys.argv[4]
# The number of surrogates to create for each significance test of a TE value
NUM_SURROGATES_PER_TE_VAL = 100
# The p level below which the null hypothesis will be rejected.
P_LEVEL = 0.05
# The number of nearest neighbours to consider in the TE estimation.
KNNS = 10
# The number of random sample points laid down will be NUM_SAMPLES_MULTIPLIER * length_of_target_train
NUM_SAMPLES_MULTIPLIER = 5.0
#SURROGATE_NUM_SAMPLES_MULTIPLIER = 5.0
# As above, but for the creation of surrogates
SURROGATE_NUM_SAMPLES_MULTIPLIER = 5.0
# The number of nearest neighbours to consider when using the local permutation method to create surrogates
K_PERM = 20
# The level of the noise to add to the random sample points used in creating surrogates
JITTERING_LEVEL = 2000
# When MAX_NUM_SECOND_INTERVALS sources have 2 or more history intervals added into the conditioning set, the inference stops
MAX_NUM_SECOND_INTERVALS = 2
# Exclude target spikes beyond this number
MAX_NUM_TARGET_SPIKES = int(num_spikes_string)
# The spikes file with the below name is expected to contain a single pickled Python list. This list contains numpy arrays. Each
# numpy array contains the spike times of each candidate target.
SPIKES_FILE_NAME = "spikes_LIF_" + net_type_name + "_" + repeat_num_string + ".pk"
# The ground truth file of the below name is expected to contain a single pickled Python list. This list contains tuples of the format(source, target).
# source and target are integers of the indices of true connections.
GROUND_TRUTH_FILE_NAME = "connections_LIF_"+ net_type_name + "_" + repeat_num_string + ".pk"
OUTPUT_FILE_PREFIX = "results/inferred_sources_target_2_" + net_type_name + "_" + num_spikes_string + "_" + repeat_num_string + "_" + target_index_string
LOG_FILE_NAME = "logs/" + net_type_name + "_" + num_spikes_string + "_" + repeat_num_string + "_" + target_index_string + ".log"
log = open(LOG_FILE_NAME, "w")
sys.stdout = log
def prepare_conditional_trains(calc_object, cond_set, spikes):
cond_trains = []
calc_object.clearConditionalIntervals()
if len(cond_set) > 0:
for key in cond_set.keys():
cond_trains.append(spikes[key])
calc_object.appendConditionalIntervals(JArray(JInt, 1)(cond_set[key]))
return cond_trains
def set_target_embeddings(embedding_list, calc_object):
if len(embedding_list) > 0:
embedding_string = str(embedding_list[0])
for i in range(2, len(embedding_list)):
embedding_string += "," + str(embedding_list[i])
calc_object.setProperty("DEST_PAST_INTERVALS", embedding_string)
else:
calc_object.setProperty("DEST_PAST_INTERVALS", "")
target_index = int(target_index_string)
print("\n****** Network inference for target neuron", target_index, "******\n\n")
# Setup JIDT
jarLocation = os.path.join(os.getcwd(), "../jidt/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?")
startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
teCalcClass = JPackage("infodynamics.measures.spiking.integration").TransferEntropyCalculatorSpikingIntegration
teCalc = teCalcClass()
teCalc.setProperty("knns", str(KNNS))
teCalc.setProperty("NUM_SAMPLES_MULTIPLIER", str(NUM_SAMPLES_MULTIPLIER))
teCalc.setProperty("SURROGATE_NUM_SAMPLES_MULTIPLIER", str(SURROGATE_NUM_SAMPLES_MULTIPLIER))
teCalc.setProperty("K_PERM", str(K_PERM))
teCalc.setProperty("DO_JITTERED_SAMPLING", "true")
teCalc.setProperty("JITTERED_SAMPLING_NOISE_LEVEL", str(JITTERING_LEVEL))
# Load spikes and ground truth connectivity
spikes = pickle.load(open(SPIKES_FILE_NAME, 'rb'))
cons = pickle.load(open(GROUND_TRUTH_FILE_NAME, 'rb'))
if MAX_NUM_TARGET_SPIKES < len(spikes[target_index]):
spikes[target_index] = spikes[target_index][:MAX_NUM_TARGET_SPIKES]
print("Number of target spikes: ", len(spikes[target_index]), "\n\n")
# First determine the correct target embedding
target_embedding_set = [1]
next_target_interval = 2
still_significant = True
print("**** Determining target embedding set ****\n")
while still_significant:
set_target_embeddings(target_embedding_set, teCalc)
teCalc.setProperty("SOURCE_PAST_INTERVALS", str(next_target_interval))
teCalc.startAddObservations()
teCalc.addObservations(JArray(JDouble, 1)(spikes[target_index]), JArray(JDouble, 1)(spikes[target_index]))
teCalc.finaliseAddObservations();
TE = teCalc.computeAverageLocalOfObservations()
sig = teCalc.computeSignificance(NUM_SURROGATES_PER_TE_VAL, TE)
print("candidate interval:", next_target_interval, " TE:", TE, " p val:", sig.pValue)
if sig.pValue > P_LEVEL:
print("Lost significance, end of target embedding determination")
still_significant = False
else:
target_embedding_set.append(next_target_interval)
next_target_interval += 1
print("target embedding set:", target_embedding_set, "\n\n")
# Now add the sources
# cond_set is a dictionary where keys are added sources and values are lists of included intervals for the
# source key.
cond_set = dict()
# next_interval_for_each_candidate will be a matrix with two columns
# first column has the source indices, second has the next interval that will be considered
next_interval_for_each_candidate = np.arange(0, len(spikes), dtype = np.intc)
next_interval_for_each_candidate = next_interval_for_each_candidate[next_interval_for_each_candidate != target_index]
next_interval_for_each_candidate = np.column_stack((next_interval_for_each_candidate, np.ones(len(next_interval_for_each_candidate), dtype = np.intc)))
still_significant = True
TE_vals_at_each_round = []
surrogate_vals_at_each_round = []
print("**** Adding Sources ****\n")
num_twos = 0
while still_significant:
print("Current conditioning set:")
for key in cond_set.keys():
print("source", key, "intervals", cond_set[key])
print("\nEstimating TE on candidate sources")
cond_trains = prepare_conditional_trains(teCalc, cond_set, spikes)
TE_vals = np.zeros(next_interval_for_each_candidate.shape[0])
debiased_TE_vals = -1 * np.ones(next_interval_for_each_candidate.shape[0])
surrogate_vals = -1 * np.ones((next_interval_for_each_candidate.shape[0], NUM_SURROGATES_PER_TE_VAL))
debiased_surrogate_vals = 1 - np.ones((next_interval_for_each_candidate.shape[0], NUM_SURROGATES_PER_TE_VAL))
is_con = np.zeros(next_interval_for_each_candidate.shape[0])
for i in range(next_interval_for_each_candidate.shape[0]):
if len(spikes[next_interval_for_each_candidate[i, 0]]) < 10:
continue
teCalc.startAddObservations()
teCalc.setProperty("SOURCE_PAST_INTERVALS", str(next_interval_for_each_candidate[i, 1]))
if len(cond_set) > 0:
teCalc.addObservations(JArray(JDouble, 1)(spikes[next_interval_for_each_candidate[i, 0]]),
JArray(JDouble, 1)(spikes[target_index]), JArray(JDouble, 2)(cond_trains))
else:
teCalc.addObservations(JArray(JDouble, 1)(spikes[next_interval_for_each_candidate[i, 0]]),
JArray(JDouble, 1)(spikes[target_index]))
teCalc.finaliseAddObservations();
TE_vals[i] = teCalc.computeAverageLocalOfObservations()
is_con[i] = ([next_interval_for_each_candidate[i, 0], target_index] in cons)
sig = teCalc.computeSignificance(NUM_SURROGATES_PER_TE_VAL, TE_vals[i])
surrogate_vals[i] = sig.distribution
debiased_TE_vals[i] = TE_vals[i] - np.mean(surrogate_vals[i])
debiased_surrogate_vals[i] = sig.distribution - np.mean(surrogate_vals[i])
print("Source", next_interval_for_each_candidate[i, 0], "Interval", next_interval_for_each_candidate[i, 1],
" TE:", str(debiased_TE_vals[i]))
log.flush()
TE_vals_at_each_round.append(TE_vals)
surrogate_vals_at_each_round.append(surrogate_vals)
sorted_TE_indices = np.argsort(debiased_TE_vals)
print("\nSorted order of sources:\n", next_interval_for_each_candidate[:, 0][sorted_TE_indices[:]])
print("Ground truth for sorted order:\n", is_con[sorted_TE_indices[:]])
index_of_max_candidate = sorted_TE_indices[-1]
samples_from_max_dist = np.max(debiased_surrogate_vals, axis = 0)
np.sort(samples_from_max_dist)
index_of_first_greater_than_estimate = np.searchsorted(samples_from_max_dist > debiased_TE_vals[index_of_max_candidate], 1)
p_val = (NUM_SURROGATES_PER_TE_VAL - index_of_first_greater_than_estimate)/float(NUM_SURROGATES_PER_TE_VAL)
print("\nMaximum candidate is source", next_interval_for_each_candidate[index_of_max_candidate, 0],
"interval", next_interval_for_each_candidate[index_of_max_candidate, 1])
print("p: ", p_val)
if p_val <= P_LEVEL:
if (next_interval_for_each_candidate[index_of_max_candidate, 0]) in cond_set:
cond_set[next_interval_for_each_candidate[index_of_max_candidate, 0]].append(next_interval_for_each_candidate[index_of_max_candidate, 1])
else:
cond_set[next_interval_for_each_candidate[index_of_max_candidate, 0]] = [next_interval_for_each_candidate[index_of_max_candidate, 1]]
if next_interval_for_each_candidate[index_of_max_candidate, 1] == 2:
num_twos += 1
if num_twos >= MAX_NUM_SECOND_INTERVALS:
print("\nMaximum number of second intervals reached\n\n")
still_significant = False
next_interval_for_each_candidate[index_of_max_candidate, 1] += 1
print("\nCandidate added\n\n")
else:
still_significant = False
print("\nLost Significance\n\n")
print("**** Pruning Sources ****\n")
# Repeatedly removes the connection that has the lowest TE out of all insignificant connections.
# Only considers the furthest intervals as candidates in each round.
everything_significant = False
while not everything_significant:
print("Current conditioning set:")
for key in cond_set.keys():
print("source", key, "intervals", cond_set[key])
print("\nEstimating TE on candidate sources")
everything_significant = True
insignificant_sources = []
insignificant_sources_TE = []
for candidate_source in cond_set:
cond_set_minus_candidate = copy.deepcopy(cond_set)
# If more than one interval, remove the last
if len(cond_set_minus_candidate[candidate_source]) > 1:
cond_set_minus_candidate[candidate_source] = cond_set_minus_candidate[candidate_source][:-1]
# Otherwise, remove source from dict
else:
cond_set_minus_candidate.pop(candidate_source)
teCalc.setProperty("SOURCE_PAST_INTERVALS", str(cond_set[candidate_source][-1]))
cond_trains = prepare_conditional_trains(teCalc, cond_set_minus_candidate, spikes)
teCalc.startAddObservations()
if len(cond_set_minus_candidate) > 0:
teCalc.addObservations(JArray(JDouble, 1)(spikes[candidate_source]), JArray(JDouble, 1)(spikes[target_index]), JArray(JDouble, 2)(cond_trains))
else:
teCalc.addObservations(JArray(JDouble, 1)(spikes[candidate_source]), JArray(JDouble, 1)(spikes[target_index]))
teCalc.finaliseAddObservations();
TE = teCalc.computeAverageLocalOfObservations()
sig = teCalc.computeSignificance(NUM_SURROGATES_PER_TE_VAL, TE)
print("Source", candidate_source, "Interval", cond_set[candidate_source][-1],
" TE:", str(round(TE, 2)), " p val:", sig.pValue)
if sig.pValue > P_LEVEL:
everything_significant = False
insignificant_sources.append(candidate_source)
insignificant_sources_TE.append(TE)
if not everything_significant:
min_TE_source = insignificant_sources[np.argmin(insignificant_sources_TE)]
print("removing source", min_TE_source, "interval", cond_set[min_TE_source][-1])
if len(cond_set[min_TE_source]) > 1:
cond_set[min_TE_source] = cond_set[min_TE_source][:-1]
else:
cond_set.pop(min_TE_source)
print("\n\n****** Final Inferred Source Set ******\n")
for key in cond_set.keys():
print("source", key, "intervals", cond_set[key])
print("\nTrue Sources:")
for con in cons:
if con[1] == target_index:
print(con[0], " ",)
output_file = open(OUTPUT_FILE_PREFIX + ".pk", 'wb')
pickle.dump(cond_set, output_file)
#pickle.dump(surrogate_vals_at_each_round, output_file)
#pickle.dump(TE_vals_at_each_round, output_file)
output_file.close()

46
demos/python/EffectiveNetworkInference/spiking/spk_to_pk.py
View File

@ -1,46 +0,0 @@
# This script converts CSV files of spike times (e.g. from the Wagenaar data set) into
# pickle files of spike times in the format that the net_inf.py script expects
import numpy as np
import pickle
import sys
import ast
import matplotlib.pyplot as plt
RUN = "1-1-20.2"
spk_file = open('extracted_data_wagenaar/1-1/' + RUN + '.spk', 'r')
time_upper = 8 * 60 * 60 * 2.5e4
spikes = []
for line in spk_file:
line = line.strip()
line = line.split(",")
line = [float(time) for time in line if time != ""]
spikes.append(np.array(line))
start_times = [train[0] for train in spikes if len(train) > 0]
lowest_start_time = min(start_times)
cutoff_time = lowest_start_time + time_upper
for i in range(len(spikes)):
spikes[i] = spikes[i][spikes[i] < cutoff_time]
spikes[i] = spikes[i] - lowest_start_time
spikes[i] = spikes[i] + np.random.uniform(size = spikes[i].shape) - 0.5
spikes[i] = np.sort(spikes[i])
print(len(spikes))
for i in range(len(spikes)):
print(spikes[i].shape)
print(spikes[i][:10])
#plt.eventplot(spikes, linewidth = 0.5)
#plt.show()
spikes_file = open("spikes_LIF_" + RUN + "_" + sys.argv[1] + ".pk", "wb")
pickle.dump(spikes, spikes_file)
cons = [[0, 0]]
connections_file = open("connections_LIF_" + RUN + "_" + sys.argv[1] + ".pk", "wb")
pickle.dump(cons, connections_file)
spikes_file.close()
connections_file.close()

BIN
demos/python/README-PythonDemos.pdf Normal file → Executable file
View File

Binary file not shown.

BIN
demos/python/README-UseInPython.pdf
View File

Binary file not shown.

191
demos/python/SpikingTE/SpikeTrainTETesting.py
View File

@ -1,191 +0,0 @@
##
## Java Information Dynamics Toolkit (JIDT)
## Copyright (C) 2022, David P. Shorten, 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/>.
##
# Transfer entropy (TE) calculation on generated spike train data using the continuous-time TE estimator.
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")
# It is recommended that this is never set to 'true', apart from cases of extremely bursty spiking (that is, long periods
# of no activity and short periods of intense spiking). In such cases, care must also be taken in the setting of the
# parameter JITTERED_SAMPLING_NOISE_LEVEL.
teCalc.setProperty("DO_JITTERED_SAMPLING", "false")
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("DO_JITTERED_SAMPLING", "false")
#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("DO_JITTERED_SAMPLING", "false")
#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("DO_JITTERED_SAMPLING", "false")
#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))

87
demos/python/create_env.bat
View File

@ -1,87 +0,0 @@
@ECHO OFF
REM
REM Java Information Dynamics Toolkit (JIDT)
REM Copyright (C) 2022, Joseph T. Lizier
REM
REM This program is free software: you can redistribute it and/or modify
REM it under the terms of the GNU General Public License as published by
REM the Free Software Foundation, either version 3 of the License, or
REM (at your option) any later version.
REM
REM This program is distributed in the hope that it will be useful,
REM but WITHOUT ANY WARRANTY; without even the implied warranty of
REM MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
REM GNU General Public License for more details.
REM
REM You should have received a copy of the GNU General Public License
REM along with this program. If not, see <http://www.gnu.org/licenses/>.
REM
REM Create a python environment (stored in folder %folder%) with jpype1, numpy and scipy installed
REM Name of folder to use and python commands -- change if required:
set folder=jpype_env
set pythonCmd=python
set pipCmd=pip
REM First make sure that the virtualenv package is installed.
%pipCmd% show virtualenv >nul 2>&1
if %errorlevel% == 0 (
echo virtualenv already installed, proceeding
) else (
echo installing virtualenv with %pipCmd% ...
%pythonCmd% -m pip install --user virtualenv
REM %errorlevel% doesnt seem to return as expect from the above, so checking success via pip:
%pipCmd% show virtualenv >nul 2>&1
if %errorlevel% neq 0 (
echo pip install of virtualenv failed
exit /b 1
) else (
echo pip install of virtualenv succeeded
)
)
REM Create a python environment (stored in folder %folder%)
%pythonCmd% -m venv %folder%
if %errorlevel% neq 0 (
REM Virtual environment creation did not work:
echo Virtual environment creation did not work. Do you need to pip install virtualenv? >&2
exit /b 2
) else (
echo Virtual environment created in %folder%
)
REM enter the environment
call %folder%\Scripts\activate.bat
if %errorlevel% neq 0 (
echo Virtual environment unable to be activated
exit /b 3
) else (
echo Python environment started and activated.
echo Beginning pip installations for the environment
)
REM install jpype1 and numpy (does not matter if they are already installed)
%pipCmd% install jpype1
%pipCmd% install numpy
echo.
echo jpype1 and numpy installed - you have a functional installation.
echo.
echo Now trying scipy, matplotlib and jupyter, but they are optional...
echo.
%pipCmd% install scipy
%pipCmd% install matplotlib
%pipCmd% install jupyter
echo.
echo scipy, matplotlib and jupyter installed
echo.
echo.
echo In Powershell activate the environment via calling: %folder%\Scripts\Activate.ps1
echo Otherwise activate the environment via calling: %folder%\Scripts\activate.bat
deactivate

84
demos/python/create_env.sh
View File

@ -1,84 +0,0 @@
#!/bin/bash
##
## Java Information Dynamics Toolkit (JIDT)
## Copyright (C) 2022, 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/>.
##
# Create a python environment (stored in folder $folder) with jpype1, numpy and scipy installed
# Name of folder to use and python commands
folder=jpype_env
pythonCmd=python3
pipCmd=pip3
# First make sure that the virtualenv package is installed.
$pipCmd show virtualenv > /dev/null 2>&1
if [ $? -eq 0 ]; then
echo "virtualenv already installed, proceeding"
else
echo "installing virtualenv with $pipCmd ...".
# On ubuntu, one could also install via the main package manager, e.g. sudo apt-get install python3-venv (I think this takes care of the followng anyway, but am unsure)
$pythonCmd -m pip install --user virtualenv
if [ $? -ne 0 ]; then
echo "pip install of virtualenv failed"
exit 1
else
echo "pip install of virtualenv succeeded"
fi
fi
# Create a python environment (stored in folder $folder)
$pythonCmd -m venv $folder
if [ $? -ne 0 ]; then
# Virtual environment creation did not work:
echo "Virtual environment creation did not work." >&2
echo "If you are on ubuntu you should now run: sudo apt-get install python3-venv" >&2
echo "Then run this script again" >&2
exit 2
else
echo "Virtual environment created in $folder"
fi
# enter the environment
source $folder/bin/activate
if [ $? -ne 0 ]; then
echo "Virtual environment unable to be activated" >&2
exit 3
else
echo "Python environment started and activated."
echo "Beginning pip installations for the environment"
fi
# install jpype1 and numpy (does not matter if they are already installed)
$pipCmd install jpype1
$pipCmd install numpy
echo
echo "Jpype1 and numpy installed - you have a functional installation."
echo
echo "Now trying scipy, matplotlib and jupyter, but they are optional..."
echo
$pipCmd install scipy
$pipCmd install matplotlib
$pipCmd install jupyter
echo
echo "Scipy, matplotlib and jupyter installed"
echo
deactivate

11
demos/python/example1TeBinaryData.py Normal file → Executable file
View File

@ -23,12 +23,9 @@
import jpype
import random
import numpy
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
jpype.startJVM(jpype.getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
@ -52,8 +49,8 @@ print("For random source, result should be close to 0 bits: %.4f" % teCalc.compu
# Next, demonstrate how to do this with a numpy array
teCalc.initialise()
# Create the numpy arrays:
sourceNumpy = numpy.array(sourceArray, dtype=int)
destNumpy = numpy.array(destArray, dtype=int)
sourceNumpy = numpy.array(sourceArray, dtype=numpy.int)
destNumpy = numpy.array(destArray, dtype=numpy.int)
# The above can be passed straight through to JIDT in python 2:
# teCalc.addObservations(sourceNumpy, destNumpy)
# But you need to do this in python 3:

7
demos/python/example2TeMultidimBinaryData.py
View File

@ -24,12 +24,9 @@
from jpype import *
import random
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# 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)

7
demos/python/example3TeContinuousDataKernel.py
View File

@ -23,12 +23,9 @@
from jpype import *
import random
import math
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# 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)

7
demos/python/example4TeContinuousDataKraskov.py
View File

@ -23,12 +23,9 @@
from jpype import *
import random
import math
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# 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)

7
demos/python/example5TeBinaryMultivarTransfer.py
View File

@ -23,12 +23,9 @@
from jpype import *
import random
from operator import xor
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# 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)

11
demos/python/example6DynamicCallingMutualInfo.py
View File

@ -31,12 +31,9 @@ import random
import string
import numpy
import readFloatsFile
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# 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)
@ -92,7 +89,7 @@ miCalc = miCalcClass()
# a. Initialise the calculator for a univariate calculation:
miCalc.initialise(1, 1)
# b. Supply the observations to compute the PDFs from:
miCalc.setObservations(JArray(JDouble, 1)(univariateSeries1.tolist()), JArray(JDouble, 1)(univariateSeries2.tolist()))
miCalc.setObservations(univariateSeries1, univariateSeries2)
# c. Make the MI calculation:
miUnivariateValue = miCalc.computeAverageLocalOfObservations()
@ -103,7 +100,7 @@ miUnivariateValue = miCalc.computeAverageLocalOfObservations()
# to use the required number of dimensions for each variable:
miCalc.initialise(len(jointVariable1Columns), len(jointVariable2Columns))
# b. Supply the observations to compute the PDFs from:
miCalc.setObservations(JArray(JDouble, 2)(jointVariable1.tolist()), JArray(JDouble, 2)(jointVariable2.tolist()))
miCalc.setObservations(jointVariable1, jointVariable2)
# c. Make the MI calculation:
miJointValue = miCalc.computeAverageLocalOfObservations()

7
demos/python/example7EnsembleMethodTeContinuousDataKraskov.py
View File

@ -24,12 +24,9 @@
from jpype import *
import random
import math
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# 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)

8
demos/python/example9TeKraskovAutoEmbedding.py
View File

@ -21,17 +21,15 @@
# Transfer entropy (TE) calculation on continuous-valued data using the Kraskov-estimator TE calculator,
# with automatic selection of embedding parameters
from jpype import *
import random
import math
import numpy
import readFloatsFile
import os
# 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?")
# Change location of jar to match yours:
jarLocation = "../../infodynamics.jar"
# 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)

42
demos/python/start_env.sh
View File

@ -1,42 +0,0 @@
#!/bin/bash
##
## Java Information Dynamics Toolkit (JIDT)
## Copyright (C) 2022, 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/>.
##
# Start the python environment (stored in folder $folder) with jpype1, numpy and others installed
# Name of folder to use and python commands
folder=jpype_env
# enter the environment
source $folder/bin/activate
if [ $? -ne 0 ]; then
echo "Virtual environment unable to be activated" >&2
# Try return first in case this script was sourced.
return 3 2>/dev/null
exit 3
else
echo "Python environment from $folder started and activated."
fi
echo
echo "Make sure you called this script as: source start_env.sh"
echo
echo "If you called it like that, you will have your python environment activated."
echo "If you just called ./start_env.sh go back and run again as above"
echo

41
demos/python/start_jupyter_env.sh
View File

@ -1,41 +0,0 @@
#!/bin/bash
##
## Java Information Dynamics Toolkit (JIDT)
## Copyright (C) 2022, 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/>.
##
# Start the python environment (stored in folder $folder) with jpype1, numpy and others installed
# and then launch jupyter
# Name of folder to use and python commands
folder=jpype_env
# enter the environment
source $folder/bin/activate
if [[ "$VIRTUAL_ENV" != "" ]]
then
echo "Started virtual environment, now starting jupyter ..."
else
echo "Virtual environment not started, exiting"
# Try return first in case this script was sourced.
return 1 2>/dev/null
exit 1
fi
cd ../..
jupyter notebook

45
demos/python/writeFloatsFile.py
View File

@ -1,45 +0,0 @@
##
## Java Information Dynamics Toolkit (JIDT)
## Copyright (C) 2020, 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/>.
##
def writeFloatsFile(filename, array):
"Write a 2D array of floats to a given file"
with open(filename, "w") as f:
# Space separate numbers, one time step per line, each column is a variable
for item in array:
# write all items
if iterable(item):
# Assume this item is a row with several columns of data
first = True;
for subitem in item:
if (not(first)):
f.write(" ");
f.write("%f" % subitem)
first = False;
else:
f.write("%f" % item)
f.write("\n")
f.close()
def iterable(a):
try:
iter(a)
except Exception:
return False
else:
return True

45
demos/python/writeIntsFile.py
View File

@ -1,45 +0,0 @@
##
## Java Information Dynamics Toolkit (JIDT)
## Copyright (C) 2020, 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/>.
##
def writeIntsFile(filename, array):
"Write a 2D array of ints to a given file"
with open(filename, "w") as f:
# Space separate numbers, one time step per line, each column is a variable
for item in array:
# write all items
if iterable(item):
# Assume this item is a row with several columns of data
first = True;
for subitem in item:
if (not(first)):
f.write(" ");
f.write("%d" % subitem)
first = False;
else:
f.write("%d" % item)
f.write("\n")
f.close()
def iterable(a):
try:
iter(a)
except Exception:
return False
else:
return True

BIN
demos/r/README-R_Demos.pdf Normal file → Executable file
View File

Binary file not shown.

BIN
demos/r/README-UseInR.pdf
View File

Binary file not shown.

3
demos/r/example3TeContinuousDataKernel.r
View File

@ -60,6 +60,5 @@ nullDist <- .jcall(teCalc,"Linfodynamics/utils/EmpiricalMeasurementDistribution;
"computeSignificance", 100L)
cat("Null distribution for unrelated source and destination",
"(i.e. the bias) has mean", .jcall(nullDist, "D", "getMeanOfDistribution"),
"bits and standard deviation", .jcall(nullDist, "D", "getStdOfDistribution"),
", while the above measurement is beaten by a proportion of", nullDist$pValue, "of the null distribution\n")
"bits and standard deviation", .jcall(nullDist, "D", "getStdOfDistribution"), "\n")

203
java/source/infodynamics/demos/autoanalysis/AutoAnalyser.java
View File

@ -28,7 +28,6 @@ import infodynamics.utils.EmpiricalNullDistributionComputer;
import infodynamics.utils.MatrixUtils;
import javax.swing.BorderFactory;
import javax.swing.DefaultCellEditor;
import javax.swing.ImageIcon;
import javax.swing.JComponent;
import javax.swing.JFrame;
@ -50,7 +49,6 @@ import javax.swing.event.ChangeListener;
import javax.swing.event.DocumentEvent;
import javax.swing.event.DocumentListener;
import javax.swing.table.AbstractTableModel;
import javax.swing.table.TableCellEditor;
import javax.swing.table.TableCellRenderer;
import javax.swing.table.TableColumn;
@ -110,13 +108,11 @@ public abstract class AutoAnalyser extends JFrame
protected String[] discreteProperties; // Children to initialise
protected String[] discretePropertyDefaultValues; // Children to initialise
protected String[] discretePropertyDescriptions; // Children to initialise
protected String[][] discretePropertyValueChoices; // Children to initialise
// Common property names for all continuous calculators:
protected String[] commonContPropertyNames;
protected String[] commonContPropertiesFieldNames;
protected String[] commonContPropertyDescriptions;
protected String[][] commonContPropertyValueChoices;
// Children can define properties for specific continuous
// calculators
@ -183,8 +179,6 @@ public abstract class AutoAnalyser extends JFrame
protected int numPermutationsToCheck = 100;
// Table for the properties
protected JTable propertiesTable;
// Default editor for the property values
protected TableCellEditor propertiesDefaultEditor;
// Table model (local class) for the table
protected PropertiesTableModel propertiesTableModel;
// Names of the properties
@ -193,8 +187,6 @@ public abstract class AutoAnalyser extends JFrame
protected Vector<String> propertyFieldNames;
// Descriptions of the fields for the properties
protected Vector<String> propertyDescriptions;
// Lists of drop-down options for the properties
protected Vector<String[]> propertyValueChoices;
// Values of the properties
protected HashMap<String,String> propertyValues;
// CheckBox for "Compute result?"
@ -221,7 +213,7 @@ public abstract class AutoAnalyser extends JFrame
protected String pathToAutoAnalyserDir = "";
// Main JIDT git/distribution folder, inferred from pathToAutoAnalyserDir
protected String jidtFolder = "";
public class TextAreaWithImage extends JTextArea {
/**
@ -302,7 +294,7 @@ public abstract class AutoAnalyser extends JFrame
Image watermarkImage = (new ImageIcon(pathToAutoAnalyserDir + "JIDT-logo-watermark.png")).getImage();
setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
setSize(1100,680);
setSize(1100,670);
setTitle(appletTitle);
// Centre in the middle of the screen
setLocationRelativeTo(null);
@ -374,17 +366,13 @@ public abstract class AutoAnalyser extends JFrame
putCalcPropertiesInTable();
propertiesTableModel = new PropertiesTableModel();
propertiesTable = new TableWithToolTip(propertiesTableModel);
// Get the default editor for the properties values:
propertiesDefaultEditor = propertiesTable.getDefaultEditor(
propertiesTable.getColumnClass(1));
System.out.println("Default properties editor is " + propertiesDefaultEditor.getClass().getName());
// Make sure any properties are saved when the compute button is clicked
propertiesTable.putClientProperty("terminateEditOnFocusLost", Boolean.TRUE);
Font headerFont = propertiesTable.getTableHeader().getFont();
propertiesTable.getTableHeader().setFont(headerFont.deriveFont(Font.BOLD));
TableColumn valueColumn = propertiesTable.getColumn("Property value");
valueColumn.setMinWidth(170);
valueColumn.setMaxWidth(170);
valueColumn.setMinWidth(130);
valueColumn.setMaxWidth(130);
JScrollPane propsTableScrollPane = new JScrollPane(propertiesTable);
// Set up for ~18 rows maximum (the +6 is exact to fit all props
// for Kraskov TE in without scrollbar)
@ -395,7 +383,7 @@ public abstract class AutoAnalyser extends JFrame
propsTableScrollPane.setMinimumSize(
new Dimension(d.width,rowHeight*17+6));
System.out.println("Row height was " + rowHeight);
// Checkbox for compute result
computeResultCheckBox = new JCheckBox("Compute result?");
computeResultCheckBox.setToolTipText("Compute result or only generate code?");
@ -430,7 +418,7 @@ public abstract class AutoAnalyser extends JFrame
javaAreaScrollPane.setHorizontalScrollBarPolicy(
JScrollPane.HORIZONTAL_SCROLLBAR_AS_NEEDED);
int codeTextAreaWidth = 560;
int codeTextAreaHeight = 530;
int codeTextAreaHeight = 480;
Dimension codeTextAreaDimension =
new Dimension(codeTextAreaWidth, codeTextAreaHeight);
javaAreaScrollPane.setPreferredSize(codeTextAreaDimension);
@ -665,8 +653,8 @@ public abstract class AutoAnalyser extends JFrame
}
protected void loadData(boolean isInts) {
ArrayFileReader afr = new ArrayFileReader(dataFile);
try {
ArrayFileReader afr = new ArrayFileReader(dataFile);
if (isInts) {
dataDiscrete = afr.getInt2DMatrix();
dataRows = dataDiscrete.length;
@ -728,35 +716,33 @@ public abstract class AutoAnalyser extends JFrame
int[] singleCalcColumns = new int[numVariables];
Vector<int[]> variableCombinations = new Vector<int[]>();
try {
if (!allCombosCheckBox.isSelected()) {
// we're doing a single combination
for (int i = 0; i < numVariables; i++) {
singleCalcColumns[i] = Integer.parseInt(variableColTextFields[i].getText());
if ((singleCalcColumns[i] < 0) || (singleCalcColumns[i] >= dataColumns)) {
JOptionPane.showMessageDialog(this,
String.format("%s column must be between 0 and %d for this data set",
variableColNumLabels[i], dataColumns-1));
resultsLabel.setText(" ");
return;
}
}
}
if (computeResultCheckBox.isSelected()) {
// Only need variableCombinations filled out if we're computing
if (computeResultCheckBox.isSelected()) {
// Only need variableCombinations filled out if we're computing
try {
if (allCombosCheckBox.isSelected()) {
// We're doing all combinations
fillOutAllCombinations(variableCombinations);
} else {
// we're doing a single combination
for (int i = 0; i < numVariables; i++) {
singleCalcColumns[i] = Integer.parseInt(variableColTextFields[i].getText());
if ((singleCalcColumns[i] < 0) || (singleCalcColumns[i] >= dataColumns)) {
JOptionPane.showMessageDialog(this,
String.format("%s column must be between 0 and %d for this data set",
variableColNumLabels[i], dataColumns-1));
resultsLabel.setText(" ");
return;
}
}
variableCombinations.add(singleCalcColumns);
}
} catch (Exception e) {
// Catches number format exception, and column number out of bounds
JOptionPane.showMessageDialog(this,
e.getMessage());
resultsLabel.setText("Cannot parse a column number from input: " + e.getMessage());
return;
}
} catch (Exception e) {
// Catches number format exception, and column number out of bounds
JOptionPane.showMessageDialog(this,
e.getMessage());
resultsLabel.setText("Cannot parse a column number from input: " + e.getMessage());
return;
}
// Generate headers:
@ -785,11 +771,11 @@ public abstract class AutoAnalyser extends JFrame
String jarLocation, pythonDemosLocation, matlabDemosLocation;
try {
File jarLocationFile = new File(jidtFolder + "infodynamics.jar");
jarLocation = jarLocationFile.getCanonicalPath().replace("\\", "\\\\");
jarLocation = jarLocationFile.getCanonicalPath();
File pythonDemosLocationFile = new File(pathToAutoAnalyserDir + "../python");
pythonDemosLocation = pythonDemosLocationFile.getCanonicalPath().replace("\\", "\\\\");
pythonDemosLocation = pythonDemosLocationFile.getCanonicalPath();
File matlabDemosLocationFile = new File(pathToAutoAnalyserDir + "../octave");
matlabDemosLocation = matlabDemosLocationFile.getCanonicalPath().replace("\\", "\\\\");
matlabDemosLocation = matlabDemosLocationFile.getCanonicalPath();
} catch (IOException ioex) {
JOptionPane.showMessageDialog(this,
ioex.getMessage());
@ -800,7 +786,6 @@ public abstract class AutoAnalyser extends JFrame
StringBuffer pythonCode = new StringBuffer();
pythonCode.append("from jpype import *\n");
pythonCode.append("import numpy\n");
pythonCode.append("import sys\n");
pythonCode.append("# Our python data file readers are a bit of a hack, python users will do better on this:\n");
pythonCode.append("sys.path.append(\"" + pythonDemosLocation + "\")\n");
if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE)) {
@ -808,11 +793,10 @@ public abstract class AutoAnalyser extends JFrame
} else {
pythonCode.append("import readFloatsFile\n\n");
}
pythonCode.append("if (not isJVMStarted()):\n");
pythonCode.append(" # Add JIDT jar library to the path\n");
pythonCode.append(" jarLocation = \"" + jarLocation + "\"\n");
pythonCode.append(" # Start the JVM (add the \"-Xmx\" option with say 1024M if you get crashes due to not enough memory space)\n");
pythonCode.append(" startJVM(getDefaultJVMPath(), \"-ea\", \"-Djava.class.path=\" + jarLocation, convertStrings=True)\n\n");
pythonCode.append("# Add JIDT jar library to the path\n");
pythonCode.append("jarLocation = \"" + jarLocation + "\"\n");
pythonCode.append("# Start the JVM (add the \"-Xmx\" option with say 1024M if you get crashes due to not enough memory space)\n");
pythonCode.append("startJVM(getDefaultJVMPath(), \"-ea\", \"-Djava.class.path=\" + jarLocation)\n\n");
// 3. Matlab:
StringBuffer matlabCode = new StringBuffer();
matlabCode.append("% Add JIDT jar library to the path, and disable warnings that it's already there:\n");
@ -918,9 +902,9 @@ public abstract class AutoAnalyser extends JFrame
// 2. Python
pythonCode.append("# " + loadDataComment);
if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE)) {
pythonCode.append("dataRaw = readIntsFile.readIntsFile(\"" + filenameAsEscapedString + "\")\n");
pythonCode.append("dataRaw = readIntsFile.readIntsFile(\"" + dataFile.getAbsolutePath() + "\")\n");
} else {
pythonCode.append("dataRaw = readFloatsFile.readFloatsFile(\"" + filenameAsEscapedString + "\")\n");
pythonCode.append("dataRaw = readFloatsFile.readFloatsFile(\"" + dataFile.getAbsolutePath() + "\")\n");
}
pythonCode.append("# As numpy array:\n");
pythonCode.append("data = numpy.array(dataRaw)\n");
@ -938,15 +922,15 @@ public abstract class AutoAnalyser extends JFrame
}
} else {
for (int i=0; i < numVariables; i++) {
pythonCode.append(variableColNumLabels[i].toLowerCase() + " = JArray(JDouble, 1)(data[:," +
singleCalcColumns[i] + "].tolist())\n");
pythonCode.append(variableColNumLabels[i].toLowerCase() + " = data[:," +
singleCalcColumns[i] + "]\n");
}
}
pythonCode.append("\n");
}
// 3. Matlab
matlabCode.append("% " + loadDataComment);
matlabCode.append("data = load('" + filenameAsEscapedString + "');\n");
matlabCode.append("data = load('" + dataFile.getAbsolutePath() + "');\n");
if (! allCombosCheckBox.isSelected()) {
matlabCode.append("% Column indices start from 1 in Matlab:\n");
if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE)) {
@ -1128,14 +1112,13 @@ public abstract class AutoAnalyser extends JFrame
columnVariables[i] + "].tolist()), " + propertyValues.get(DISCRETE_PROPNAME_BASE) + ")\n");
}
} else {
// Continuous data
for (int i = 0; i < numVariables; i++) {
javaCode.append(javaPrefix + "double[] " + variableColNumLabels[i].toLowerCase() +
" = MatrixUtils.selectColumn(data, " + columnVariables[i] + ");\n");
matlabCode.append(matlabPrefix + variableColNumLabels[i].toLowerCase() +
" = octaveToJavaDoubleArray(data(:, " + columnVariables[i] + "));\n");
pythonCode.append(pythonPrefix + variableColNumLabels[i].toLowerCase() +
" = JArray(JDouble, 1)(data[:, " + columnVariables[i] + "].tolist())\n");
" = data[:, " + columnVariables[i] + "]\n");
}
}
javaCode.append("\n");
@ -1177,7 +1160,7 @@ public abstract class AutoAnalyser extends JFrame
javaCode.append(javaPrefix + "calc." + setObservationsMethod + "(" + methodArguments + ");\n");
// 2. Python
pythonCode.append(pythonPrefix + "# " + supplyDataComment);
pythonCode.append(pythonPrefix + "calc." + setObservationsMethod + pythonSetObsSuffix() + "(" + methodArguments + ")\n");
pythonCode.append(pythonPrefix + "calc." + setObservationsMethod + "(" + methodArguments + ")\n");
// 3. Matlab
matlabCode.append(matlabPrefix + "% " + supplyDataComment);
matlabCode.append(matlabPrefix + "calc." + setObservationsMethod + "(" + methodArguments + ");\n");
@ -1241,7 +1224,7 @@ public abstract class AutoAnalyser extends JFrame
extraFormatTerms + "result" + statSigFormatTerms + ");\n");
// 2. Python
pythonCode.append("\n" + pythonPrefix + "print(\"" + resultsPrefixString +
"%.4f " + units + resultsSuffixString + "\" %\\\n " + pythonPrefix + "(" +
"%.4f " + units + resultsSuffixString + "\" %\n " + pythonPrefix + "(" +
extraFormatTerms + "result" + statSigFormatTerms + "))\n");
// 3. Matlab
matlabCode.append("\n" + matlabPrefix + "fprintf('" + resultsPrefixString +
@ -1392,7 +1375,7 @@ public abstract class AutoAnalyser extends JFrame
System.out.println(resultsText);
}
if ((!allCombosCheckBox.isSelected()) &&
if ((allCombosCheckBox.isSelected()) &&
!(selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE) ||
(selectedCalcType.equalsIgnoreCase(CALC_TYPE_BINNED)))) {
// Read the current property values back out (in case of
@ -1520,15 +1503,7 @@ public abstract class AutoAnalyser extends JFrame
InfoMeasureCalculatorContinuous calcContinuous,
int[] columnCombo) throws Exception;
/**
* Method to return any suffix if required on the set/addObservations method for python calls
*
* @return
*/
protected String pythonSetObsSuffix() {
return "";
}
/**
* Extends JTable to add ToolTipText to the property names
*
@ -1553,22 +1528,9 @@ public abstract class AutoAnalyser extends JFrame
if (vColIndex == 0) {
JComponent jc = (JComponent)c;
try {
String toolTipText;
String selectedCalcType = (String)
calcTypeComboBox.getSelectedItem();
if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE) ||
selectedCalcType.equalsIgnoreCase(CALC_TYPE_BINNED)) {
toolTipText = "<html>" + propertyNames.get(rowIndex) + ": " +
propertyDescriptions.get(rowIndex) + "</html>";
} else {
toolTipText = "<html>" + propertyFieldNames.get(rowIndex) + ": " +
propertyDescriptions.get(rowIndex) + "</html>";
}
jc.setToolTipText(toolTipText);
jc.setToolTipText("<html>" + propertyFieldNames.get(rowIndex) + ": " + propertyDescriptions.get(rowIndex) + "</html>");
} catch (ArrayIndexOutOfBoundsException aioobe) {
// Catch if the row number was outside our array of descriptions (e.g. empty row)
System.out.println("prepareRenderer: Row number " + rowIndex +
" was outside our array of names/fieldnames/descriptions");
}
}
}
@ -1594,70 +1556,6 @@ public abstract class AutoAnalyser extends JFrame
}
}
*/
/**
* Use this method to set combo box options for editing cells
*/
@Override
public TableCellEditor getCellEditor(int row, int column) {
if ((propertyValueChoices.get(row) != null) && (column == 1)) {
// We need to construct a combo box for the selection for this property:
try {
JComboBox<String> paramChoiceComboBox = new JComboBox<String>();
// Set font to not bold and one size less than the default (to fit better)
Font font = paramChoiceComboBox.getFont();
paramChoiceComboBox.setFont(font.deriveFont(Font.PLAIN, font.getSize()-1));
String[] choices = propertyValueChoices.get(row);
for (int c = 0; c < choices.length; c++) {
paramChoiceComboBox.addItem(choices[c]);
}
return new DefaultCellEditor(paramChoiceComboBox);
} catch (Exception e) {
e.printStackTrace();
// But now allow this to be handled by the default cell editor
}
}
// I think the following would do the default behaviour:
// return this.getDefaultEditor(this.getColumnClass(column));
// however it should be safer to just allow the parent to handle:
return super.getCellEditor(row, column);
}
/**
* This method allows us to keep a combo box visible when
* the property value is no longer selected.
* Adapted from answer at https://stackoverflow.com/questions/30744524/how-to-make-the-jcombobox-dropdown-always-visible-in-a-jtable
*/
@Override
public TableCellRenderer getCellRenderer(int row, int column) {
if ((propertyValueChoices.get(row) != null) && (column == 1)) {
try {
return new TableCellRenderer() {
JComboBox<String> box = new JComboBox<String>();
int defaultFontSize = box.getFont().getSize();
@Override
public Component getTableCellRendererComponent(JTable table,
Object value, boolean isSelected, boolean hasFocus, int row,
int column) {
// Set font to not bold and one size less than the default (to fit better)
Font font = box.getFont();
box.setFont(font.deriveFont(Font.PLAIN, defaultFontSize-1));
// Now empty all items out and just put the value we currently have.
// (This is only for displaying, the editor will override with available values if
// user wants to edit).
box.removeAllItems();
box.addItem(value.toString());
return box;
}
};
} catch (Exception e) {
e.printStackTrace();
// But now allow this to be handled by the default cell renderer
}
}
return super.getCellRenderer(row, column);
}
}
protected class PropertiesTableModel extends AbstractTableModel {
@ -1745,8 +1643,6 @@ public abstract class AutoAnalyser extends JFrame
calcProperties.classSpecificPropertiesFieldNames;
String[] classSpecificPropertyDescriptions =
calcProperties.classSpecificPropertyDescriptions;
String[][] classSpecificPropertyValueChoices =
calcProperties.classSpecificPropertyValueChoices;
calcClass = calcProperties.calcClass;
Object calc = calcProperties.calc;
@ -1754,7 +1650,6 @@ public abstract class AutoAnalyser extends JFrame
propertyNames = new Vector<String>();
propertyFieldNames = new Vector<String>();
propertyDescriptions = new Vector<String>();
propertyValueChoices = new Vector<String[]>();
if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE) ||
selectedCalcType.equalsIgnoreCase(CALC_TYPE_BINNED)) {
@ -1765,12 +1660,10 @@ public abstract class AutoAnalyser extends JFrame
for (String propName : discreteProperties) {
String propertyDescription = discretePropertyDescriptions[i];
String defaultPropertyValue = discretePropertyDefaultValues[i];
String[] propertyValueChoiceSet = discretePropertyValueChoices[i];
i++;
propertyNames.add(propName);
propertyDescriptions.add(propertyDescription);
propertyValues.put(propName, defaultPropertyValue);
propertyValueChoices.add(propertyValueChoiceSet);
System.out.println("Adding property name " + propName);
}
} else {
@ -1779,7 +1672,6 @@ public abstract class AutoAnalyser extends JFrame
for (String fieldName : commonContPropertiesFieldNames) {
String propName = commonContPropertyNames[i];
String propertyDescription = commonContPropertyDescriptions[i];
String[] propertyValueChoiceSet = commonContPropertyValueChoices[i];
i++;
System.out.println("Adding property name " +
abstractContinuousClass.getSimpleName() + "." + fieldName +
@ -1787,7 +1679,6 @@ public abstract class AutoAnalyser extends JFrame
propertyFieldNames.add(abstractContinuousClass.getSimpleName() + "." + fieldName);
propertyNames.add(propName);
propertyDescriptions.add(propertyDescription);
propertyValueChoices.add(propertyValueChoiceSet);
}
// Then for the specific estimator types
@ -1795,11 +1686,9 @@ public abstract class AutoAnalyser extends JFrame
for (String fieldName : classSpecificPropertiesFieldNames) {
String propName = classSpecificPropertyNames[i];
String propertyDescription = classSpecificPropertyDescriptions[i];
String[] propertyValueChoiceSet = classSpecificPropertyValueChoices[i];
i++;
propertyNames.add(propName);
propertyDescriptions.add(propertyDescription);
propertyValueChoices.add(propertyValueChoiceSet);
if (fieldName.contains(".")) {
System.out.println("Adding property name " + fieldName +
" = \"" + propName + "\"");
@ -1845,7 +1734,6 @@ public abstract class AutoAnalyser extends JFrame
String[] classSpecificPropertyNames;
String[] classSpecificPropertiesFieldNames;
String[] classSpecificPropertyDescriptions;
String[][] classSpecificPropertyValueChoices;
}
/**
@ -1868,7 +1756,6 @@ public abstract class AutoAnalyser extends JFrame
calcProperties.classSpecificPropertyNames = discreteProperties;
calcProperties.classSpecificPropertiesFieldNames = null; // Not used
calcProperties.classSpecificPropertyDescriptions = discretePropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = discretePropertyValueChoices;
// TODO Later can add binning method to the properties for
// binned calculator here.
return calcProperties;

114
java/source/infodynamics/demos/autoanalysis/AutoAnalyserAIS.java
View File

@ -19,12 +19,9 @@
package infodynamics.demos.autoanalysis;
import infodynamics.measures.continuous.ActiveInfoStorageCalculator;
import infodynamics.measures.continuous.ActiveInfoStorageCalculatorViaMutualInfo;
import infodynamics.measures.continuous.InfoMeasureCalculatorContinuous;
import infodynamics.measures.continuous.gaussian.ActiveInfoStorageCalculatorGaussian;
import infodynamics.measures.continuous.gaussian.MutualInfoCalculatorMultiVariateGaussian;
import infodynamics.measures.continuous.kernel.ActiveInfoStorageCalculatorKernel;
import infodynamics.measures.continuous.kernel.ActiveInfoStorageCalculatorMultiVariateKernel;
import infodynamics.measures.continuous.kernel.MutualInfoCalculatorMultiVariateKernel;
import infodynamics.measures.continuous.kraskov.ActiveInfoStorageCalculatorKraskov;
import infodynamics.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov;
@ -58,15 +55,12 @@ public class AutoAnalyserAIS extends AutoAnalyser {
protected String[] gaussianProperties;
protected String[] gaussianPropertiesFieldNames;
protected String[] gaussianPropertyDescriptions;
protected String[][] gaussianPropertyValueChoices;
protected String[] kernelProperties;
protected String[] kernelPropertiesFieldNames;
protected String[] kernelPropertyDescriptions;
protected String[][] kernelPropertyValueChoices;
protected String[] kraskovProperties;
protected String[] kraskovPropertiesFieldNames;
protected String[] kraskovPropertyDescriptions;
protected String[][] kraskovPropertyValueChoices;
public AutoAnalyserAIS() {
super();
@ -112,84 +106,39 @@ public class AutoAnalyserAIS extends AutoAnalyser {
"Number of discrete states available for each variable (i.e. 2 for binary)",
"History embedding length (k_HISTORY)"
};
discretePropertyValueChoices = new String[][] {
null,
null
};
// Continuous:
abstractContinuousClass = ActiveInfoStorageCalculator.class;
// Common properties for all continuous calcs:
commonContPropertyNames = new String[] {
ActiveInfoStorageCalculator.K_PROP_NAME,
ActiveInfoStorageCalculator.TAU_PROP_NAME,
ActiveInfoStorageCalculatorViaMutualInfo.PROP_AUTO_EMBED_METHOD,
ActiveInfoStorageCalculatorViaMutualInfo.PROP_K_SEARCH_MAX,
ActiveInfoStorageCalculatorViaMutualInfo.PROP_TAU_SEARCH_MAX,
ActiveInfoStorageCalculator.TAU_PROP_NAME
};
commonContPropertiesFieldNames = new String[] {
"K_PROP_NAME",
"TAU_PROP_NAME",
"ActiveInfoStorageCalculatorViaMutualInfo.PROP_AUTO_EMBED_METHOD",
"ActiveInfoStorageCalculatorViaMutualInfo.PROP_K_SEARCH_MAX",
"ActiveInfoStorageCalculatorViaMutualInfo.PROP_TAU_SEARCH_MAX",
"TAU_PROP_NAME"
};
commonContPropertyDescriptions = new String[] {
"History embedding length (k_HISTORY)",
"History embedding delay (k_TAU)",
"Method to automatically determine embedding length (k_HISTORY)<br/> and delay (k_TAU) for " +
"the samples. Default is \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_NONE +
"\" meaning values are set manually; other values include: <br/> -- \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_RAGWITZ +
"\" for use of the Ragwitz criteria for both source and destination (searching up to \"" + ActiveInfoStorageCalculatorKraskov.PROP_K_SEARCH_MAX +
"\" and \"" + ActiveInfoStorageCalculatorKraskov.PROP_TAU_SEARCH_MAX + "\"); <br/> -- \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_MAX_CORR_AIS +
"\" for maximising the (bias corrected) Active Info Storage (searching up to \"" + ActiveInfoStorageCalculatorKraskov.PROP_K_SEARCH_MAX +
"\" and \"" + ActiveInfoStorageCalculatorKraskov.PROP_TAU_SEARCH_MAX + "\"); <br/>Use of values other than \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_NONE +
"\" leads to any previous settings for embedding lengths and delays to be overwritten after observations are supplied",
"Max. embedding length to search to <br/>if auto embedding (as determined by " + ActiveInfoStorageCalculatorKraskov.PROP_AUTO_EMBED_METHOD + ")",
"Max. embedding delay to search to <br/>if auto embedding (as determined by " + ActiveInfoStorageCalculatorKraskov.PROP_AUTO_EMBED_METHOD + ")",
};
commonContPropertyValueChoices = new String[][] {
null,
null,
{ActiveInfoStorageCalculatorViaMutualInfo.AUTO_EMBED_METHOD_NONE,
ActiveInfoStorageCalculatorViaMutualInfo.AUTO_EMBED_METHOD_RAGWITZ,
ActiveInfoStorageCalculatorViaMutualInfo.AUTO_EMBED_METHOD_MAX_CORR_AIS},
null,
null,
"History embedding delay (k_TAU)"
};
// Gaussian properties:
gaussianProperties = new String[] {
MutualInfoCalculatorMultiVariateGaussian.PROP_BIAS_CORRECTION,
ActiveInfoStorageCalculatorGaussian.PROP_MAX_CORR_AIS_NUM_SURROGATES
};
gaussianPropertiesFieldNames = new String[] {
"MutualInfoCalculatorMultiVariateGaussian.PROP_BIAS_CORRECTION",
"ActiveInfoStorageCalculatorGaussian.PROP_MAX_CORR_AIS_NUM_SURROGATES"
};
gaussianPropertyDescriptions = new String[] {
"Whether the analytically determined bias (as the mean of the<br/>" +
"surrogate distribution) will be subtracted from all" +
"calculated values. Default is false.",
"Number of surrogates to use in computing the bias correction<br/>if required for " +
ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_MAX_CORR_AIS + " auto-embedding method.<br/>" +
"(default is 0, meaning to use analytic method -- recommended)"
};
gaussianPropertyValueChoices = new String[][] {
{"true", "false"},
null
};
// Kernel:
kernelProperties = new String[] {
MutualInfoCalculatorMultiVariateKernel.KERNEL_WIDTH_PROP_NAME,
MutualInfoCalculatorMultiVariateKernel.DYN_CORR_EXCL_TIME_NAME,
MutualInfoCalculatorMultiVariateKernel.NORMALISE_PROP_NAME,
ActiveInfoStorageCalculatorMultiVariateKernel.PROP_MAX_CORR_AIS_NUM_SURROGATES
MutualInfoCalculatorMultiVariateKernel.NORMALISE_PROP_NAME,
};
kernelPropertiesFieldNames = new String[] {
"MutualInfoCalculatorMultiVariateKernel.KERNEL_WIDTH_PROP_NAME",
"MutualInfoCalculatorMultiVariateKernel.DYN_CORR_EXCL_TIME_NAME",
"MutualInfoCalculatorMultiVariateKernel.NORMALISE_PROP_NAME",
"ActiveInfoStorageCalculatorMultiVariateKernel.PROP_MAX_CORR_AIS_NUM_SURROGATES"
"KERNEL_WIDTH_PROP_NAME",
"DYN_CORR_EXCL_TIME_NAME",
"NORMALISE_PROP_NAME"
};
kernelPropertyDescriptions = new String[] {
"Kernel width to be used in the calculation. <br/>If the property " +
@ -198,16 +147,7 @@ public class AutoAnalyserAIS extends AutoAnalyser {
"otherwise it is an absolute value.",
"Dynamic correlation exclusion time or <br/>Theiler window (see Kantz and Schreiber); " +
"0 (default) means no dynamic exclusion window",
"(boolean) whether to normalise <br/>each incoming time-series to mean 0, standard deviation 1, or not (default true, recommended)",
"Number of surrogates to use in computing the bias correction<br/>if required for " +
ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_MAX_CORR_AIS + " auto-embedding method.<br/>" +
"(default is 20)"
};
kernelPropertyValueChoices = new String[][] {
null,
null,
{"true", "false"},
null
"(boolean) whether to normalise <br/>each incoming time-series to mean 0, standard deviation 1, or not (recommended)",
};
// KSG (Kraskov):
kraskovProperties = new String[] {
@ -217,8 +157,9 @@ public class AutoAnalyserAIS extends AutoAnalyser {
MutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME,
MutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE,
MutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
MutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU,
ActiveInfoStorageCalculatorKraskov.PROP_KRASKOV_ALG_NUM,
ActiveInfoStorageCalculatorKraskov.PROP_AUTO_EMBED_METHOD,
ActiveInfoStorageCalculatorKraskov.PROP_K_SEARCH_MAX,
ActiveInfoStorageCalculatorKraskov.PROP_TAU_SEARCH_MAX,
ActiveInfoStorageCalculatorKraskov.PROP_RAGWITZ_NUM_NNS,
};
kraskovPropertiesFieldNames = new String[] {
@ -228,9 +169,10 @@ public class AutoAnalyserAIS extends AutoAnalyser {
"MutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME",
"MutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE",
"MutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS",
"MutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU",
"PROP_KRASKOV_ALG_NUM",
"PROP_RAGWITZ_NUM_NNS"
"PROP_AUTO_EMBED_METHOD",
"PROP_K_SEARCH_MAX",
"PROP_TAU_SEARCH_MAX",
"PROP_RAGWITZ_NUM_NNS"
};
kraskovPropertyDescriptions = new String[] {
"(boolean) whether to normalise <br/>each incoming time-series to mean 0, standard deviation 1, or not (recommended)",
@ -244,21 +186,18 @@ public class AutoAnalyserAIS extends AutoAnalyser {
"\"MAX_NORM\" (default), otherwise \"EUCLIDEAN\" or \"EUCLIDEAN_SQUARED\" (both equivalent here)",
"Number of parallel threads to use <br/>in computation: an integer > 0 or \"USE_ALL\" " +
"(default, to indicate to use all available processors)",
"Whether to enable the GPU module (number of threads then has no bearing); boolean, default false",
"Which KSG algorithm to use (1 or 2)",
"Method to automatically determine embedding length (k_HISTORY)<br/> and delay (k_TAU) for " +
"the samples. Default is \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_NONE +
"\" meaning values are set manually; other values include: <br/> -- \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_RAGWITZ +
"\" for use of the Ragwitz criteria for both source and destination (searching up to \"" + ActiveInfoStorageCalculatorKraskov.PROP_K_SEARCH_MAX +
"\" and \"" + ActiveInfoStorageCalculatorKraskov.PROP_TAU_SEARCH_MAX + "\"); <br/> -- \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_MAX_CORR_AIS +
"\" for maximising the (bias corrected) Active Info Storage (searching up to \"" + ActiveInfoStorageCalculatorKraskov.PROP_K_SEARCH_MAX +
"\" and \"" + ActiveInfoStorageCalculatorKraskov.PROP_TAU_SEARCH_MAX + "\"); <br/>Use of values other than \"" + ActiveInfoStorageCalculatorKraskov.AUTO_EMBED_METHOD_NONE +
"\" leads to any previous settings for embedding lengths and delays to be overwritten after observations are supplied",
"Max. embedding length to search to <br/>if auto embedding (as determined by " + ActiveInfoStorageCalculatorKraskov.PROP_AUTO_EMBED_METHOD + ")",
"Max. embedding delay to search to <br/>if auto embedding (as determined by " + ActiveInfoStorageCalculatorKraskov.PROP_AUTO_EMBED_METHOD + ")",
"Number of k nearest neighbours for <br/>Ragwitz auto embedding (if used; defaults to match property \"k\")"
};
kraskovPropertyValueChoices = new String[][] {
{"true", "false"},
null,
null,
null,
{"MAX_NORM", "EUCLIDEAN", "EUCLIDEAN_SQUARED"},
null,
{"true", "false"},
{"1", "2"},
null,
};
}
@Override
@ -365,17 +304,14 @@ public class AutoAnalyserAIS extends AutoAnalyser {
calcProperties.classSpecificPropertyNames = gaussianProperties;
calcProperties.classSpecificPropertiesFieldNames = gaussianPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = gaussianPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = gaussianPropertyValueChoices;
} else if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_KRASKOV)) {
calcProperties.classSpecificPropertyNames = kraskovProperties;
calcProperties.classSpecificPropertiesFieldNames = kraskovPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = kraskovPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = kraskovPropertyValueChoices;
} else if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_KERNEL)) {
calcProperties.classSpecificPropertyNames = kernelProperties;
calcProperties.classSpecificPropertiesFieldNames = kernelPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = kernelPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = kernelPropertyValueChoices;
} else {
calcProperties = null;
throw new Exception("No recognised calculator selected: " +

46
java/source/infodynamics/demos/autoanalysis/AutoAnalyserCMI.java
View File

@ -56,11 +56,9 @@ public class AutoAnalyserCMI extends AutoAnalyser
protected String[] gaussianProperties;
protected String[] gaussianPropertiesFieldNames;
protected String[] gaussianPropertyDescriptions;
protected String[][] gaussianPropertyValueChoices;
protected String[] kraskovProperties;
protected String[] kraskovPropertiesFieldNames;
protected String[] kraskovPropertyDescriptions;
protected String[][] kraskovPropertyValueChoices;
protected static final String CALC_TYPE_KRASKOV_ALG1 = CALC_TYPE_KRASKOV + " alg. 1";
protected static final String CALC_TYPE_KRASKOV_ALG2 = CALC_TYPE_KRASKOV + " alg. 2";
@ -109,9 +107,6 @@ public class AutoAnalyserCMI extends AutoAnalyser
"Number of discrete states available for each variable (i.e. 2 for binary).<br/>" +
"Can be set individually for each variable -- see code."
};
discretePropertyValueChoices = new String[][] {
null
};
// Continuous:
abstractContinuousClass = ConditionalMutualInfoCalculatorMultiVariate.class;
@ -125,23 +120,12 @@ public class AutoAnalyserCMI extends AutoAnalyser
commonContPropertyDescriptions = new String[] {
// None
};
commonContPropertyValueChoices = new String[][] {
// None
};
// Gaussian properties:
gaussianProperties = new String[] {
ConditionalMutualInfoCalculatorMultiVariateGaussian.PROP_BIAS_CORRECTION,
};
gaussianPropertiesFieldNames = new String[] {
"PROP_BIAS_CORRECTION"
};
gaussianPropertyDescriptions = new String[] {
"Whether the analytically determined bias (as the mean of the<br/>" +
"surrogate distribution) will be subtracted from all" +
"calculated values. Default is false."
};
gaussianPropertyValueChoices = new String[][] {
{"true", "false"}
};
// KSG (Kraskov):
kraskovProperties = new String[] {
@ -151,7 +135,6 @@ public class AutoAnalyserCMI extends AutoAnalyser
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME,
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE,
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU,
};
kraskovPropertiesFieldNames = new String[] {
"ConditionalMutualInfoMultiVariateCommon.PROP_NORMALISE",
@ -160,7 +143,6 @@ public class AutoAnalyserCMI extends AutoAnalyser
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME",
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE",
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS",
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU"
};
kraskovPropertyDescriptions = new String[] {
"(boolean) whether to normalise <br/>each incoming time-series to mean 0, standard deviation 1, or not (recommended)",
@ -174,17 +156,8 @@ public class AutoAnalyserCMI extends AutoAnalyser
"\"MAX_NORM\" (default), otherwise \"EUCLIDEAN\" or \"EUCLIDEAN_SQUARED\" (both equivalent here)",
"Number of parallel threads to use <br/>in computation: an integer > 0 or \"USE_ALL\" " +
"(default, to indicate to use all available processors)",
"Whether to enable the GPU module (number of threads then has no bearing); boolean, default false"
};
kraskovPropertyValueChoices = new String[][] {
{"true", "false"},
null,
null,
null,
{"MAX_NORM", "EUCLIDEAN", "EUCLIDEAN_SQUARED"},
null,
{"true", "false"}
};
}
@Override
@ -337,13 +310,11 @@ public class AutoAnalyserCMI extends AutoAnalyser
calcProperties.classSpecificPropertyNames = gaussianProperties;
calcProperties.classSpecificPropertiesFieldNames = gaussianPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = gaussianPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = gaussianPropertyValueChoices;
} else if (selectedCalcType.startsWith(CALC_TYPE_KRASKOV)) {
// The if statement will work for both MI Kraskov calculators
calcProperties.classSpecificPropertyNames = kraskovProperties;
calcProperties.classSpecificPropertiesFieldNames = kraskovPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = kraskovPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = kraskovPropertyValueChoices;
} else {
calcProperties = null;
throw new Exception("No recognised calculator selected: " +
@ -392,21 +363,6 @@ public class AutoAnalyserCMI extends AutoAnalyser
base + ", " + base + ", " + base);
}
@Override
protected String pythonSetObsSuffix() {
String selectedCalcType = (String)
calcTypeComboBox.getSelectedItem();
if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE) ||
selectedCalcType.equalsIgnoreCase(CALC_TYPE_BINNED)) {
return "";
} else {
// For the moment we could direct all calls to the 1D arrays version,
// but it is working fine with JPype 0.7; later
// when we have 2D inputs we should dynamically detect that and return "2D"
return "";
}
}
/**
* @param args
*/

41
java/source/infodynamics/demos/autoanalysis/AutoAnalyserCTE.java
View File

@ -21,7 +21,6 @@ package infodynamics.demos.autoanalysis;
import infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateCommon;
import infodynamics.measures.continuous.ConditionalTransferEntropyCalculator;
import infodynamics.measures.continuous.InfoMeasureCalculatorContinuous;
import infodynamics.measures.continuous.gaussian.ConditionalMutualInfoCalculatorMultiVariateGaussian;
import infodynamics.measures.continuous.gaussian.ConditionalTransferEntropyCalculatorGaussian;
import infodynamics.measures.continuous.kraskov.ConditionalMutualInfoCalculatorMultiVariateKraskov;
import infodynamics.measures.continuous.kraskov.ConditionalTransferEntropyCalculatorKraskov;
@ -58,11 +57,9 @@ public class AutoAnalyserCTE extends AutoAnalyser
protected String[] gaussianProperties;
protected String[] gaussianPropertiesFieldNames;
protected String[] gaussianPropertyDescriptions;
protected String[][] gaussianPropertyValueChoices;
protected String[] kraskovProperties;
protected String[] kraskovPropertiesFieldNames;
protected String[] kraskovPropertyDescriptions;
protected String[][] kraskovPropertyValueChoices;
public AutoAnalyserCTE() {
super();
@ -111,10 +108,6 @@ public class AutoAnalyserCTE extends AutoAnalyser
"Can be set individually for each variable -- see code.",
"Destination history embedding length (k_HISTORY)",
};
discretePropertyValueChoices = new String[][] {
null,
null
};
// Continuous:
abstractContinuousClass = ConditionalTransferEntropyCalculator.class;
@ -149,30 +142,12 @@ public class AutoAnalyserCTE extends AutoAnalyser
"Conditional history embeding delay",
"Delay from conditional to destination (in time steps)"
};
commonContPropertyValueChoices = new String[][] {
null,
null,
null,
null,
null,
null,
null,
null
};
// Gaussian properties:
gaussianProperties = new String[] {
ConditionalMutualInfoCalculatorMultiVariateGaussian.PROP_BIAS_CORRECTION,
};
gaussianPropertiesFieldNames = new String[] {
"ConditionalMutualInfoCalculatorMultiVariateGaussian.PROP_BIAS_CORRECTION",
};
gaussianPropertyDescriptions = new String[] {
"Whether the analytically determined bias (as the mean of the<br/>" +
"surrogate distribution) will be subtracted from all" +
"calculated values. Default is false."
};
gaussianPropertyValueChoices = new String[][] {
{"true", "false"},
};
// KSG (Kraskov):
kraskovProperties = new String[] {
@ -182,7 +157,6 @@ public class AutoAnalyserCTE extends AutoAnalyser
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME,
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE,
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU,
ConditionalTransferEntropyCalculatorKraskov.PROP_KRASKOV_ALG_NUM
};
kraskovPropertiesFieldNames = new String[] {
@ -192,7 +166,6 @@ public class AutoAnalyserCTE extends AutoAnalyser
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME",
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE",
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS",
"ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU",
"PROP_KRASKOV_ALG_NUM"
};
kraskovPropertyDescriptions = new String[] {
@ -207,19 +180,9 @@ public class AutoAnalyserCTE extends AutoAnalyser
"\"MAX_NORM\" (default), otherwise \"EUCLIDEAN\" or \"EUCLIDEAN_SQUARED\" (both equivalent here)",
"Number of parallel threads to use <br/>in computation: an integer > 0 or \"USE_ALL\" " +
"(default, to indicate to use all available processors)",
"Whether to enable the GPU module (number of threads then has no bearing); boolean, default false",
"Which KSG algorithm to use (1 or 2)",
};
kraskovPropertyValueChoices = new String[][] {
{"true", "false"},
null,
null,
null,
{"MAX_NORM", "EUCLIDEAN", "EUCLIDEAN_SQUARED"},
null,
{"true", "false"},
{"1", "2"}
};
}
@Override
@ -372,13 +335,11 @@ public class AutoAnalyserCTE extends AutoAnalyser
calcProperties.classSpecificPropertyNames = gaussianProperties;
calcProperties.classSpecificPropertiesFieldNames = gaussianPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = gaussianPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = gaussianPropertyValueChoices;
} else if (selectedCalcType.startsWith(CALC_TYPE_KRASKOV)) {
// The if statement will work for both MI Kraskov calculators
calcProperties.classSpecificPropertyNames = kraskovProperties;
calcProperties.classSpecificPropertiesFieldNames = kraskovPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = kraskovPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = kraskovPropertyValueChoices;
} else {
calcProperties = null;
throw new Exception("No recognised calculator selected: " +

8
java/source/infodynamics/demos/autoanalysis/AutoAnalyserChannelCalculator.java
View File

@ -47,16 +47,13 @@ public abstract class AutoAnalyserChannelCalculator extends AutoAnalyser {
protected String[] gaussianProperties;
protected String[] gaussianPropertiesFieldNames;
protected String[] gaussianPropertyDescriptions;
protected String[][] gaussianPropertyValueChoices;
protected String[] kernelProperties;
protected String[] kernelPropertiesFieldNames;
protected String[] kernelPropertyDescriptions;
protected String[][] kernelPropertyValueChoices;
protected String[] kraskovProperties;
protected String[] kraskovPropertiesFieldNames;
protected String[] kraskovPropertyDescriptions;
protected String[][] kraskovPropertyValueChoices;
public AutoAnalyserChannelCalculator() {
super();
@ -210,18 +207,15 @@ public abstract class AutoAnalyserChannelCalculator extends AutoAnalyser {
calcProperties.classSpecificPropertyNames = gaussianProperties;
calcProperties.classSpecificPropertiesFieldNames = gaussianPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = gaussianPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = gaussianPropertyValueChoices;
} else if (selectedCalcType.startsWith(CALC_TYPE_KRASKOV)) {
// The if statement will work for both MI Kraskov calculators
calcProperties.classSpecificPropertyNames = kraskovProperties;
calcProperties.classSpecificPropertiesFieldNames = kraskovPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = kraskovPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = kraskovPropertyValueChoices;
} else if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_KERNEL)) {
calcProperties.classSpecificPropertyNames = kernelProperties;
calcProperties.classSpecificPropertiesFieldNames = kernelPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = kernelPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = kernelPropertyValueChoices;
} else {
calcProperties = null;
throw new Exception("No recognised calculator selected: " +

22
java/source/infodynamics/demos/autoanalysis/AutoAnalyserEntropy.java
View File

@ -52,15 +52,12 @@ public class AutoAnalyserEntropy extends AutoAnalyser {
protected String[] gaussianProperties;
protected String[] gaussianPropertiesFieldNames;
protected String[] gaussianPropertyDescriptions;
protected String[][] gaussianPropertyValueChoices;
protected String[] kernelProperties;
protected String[] kernelPropertiesFieldNames;
protected String[] kernelPropertyDescriptions;
protected String[][] kernelPropertyValueChoices;
protected String[] klProperties;
protected String[] klPropertiesFieldNames;
protected String[] klPropertyDescriptions;
protected String[][] klPropertyValueChoices;
public AutoAnalyserEntropy() {
super();
@ -103,9 +100,6 @@ public class AutoAnalyserEntropy extends AutoAnalyser {
discretePropertyDescriptions = new String[] {
"Number of discrete states available for each variable (i.e. 2 for binary)"
};
discretePropertyValueChoices = new String[][] {
null
};
// Continuous:
abstractContinuousClass = EntropyCalculator.class;
@ -119,9 +113,6 @@ public class AutoAnalyserEntropy extends AutoAnalyser {
commonContPropertyDescriptions = new String[] {
// None
};
commonContPropertyValueChoices = new String[][] {
// None
};
// Gaussian properties:
gaussianProperties = new String[] {
};
@ -129,8 +120,6 @@ public class AutoAnalyserEntropy extends AutoAnalyser {
};
gaussianPropertyDescriptions = new String[] {
};
gaussianPropertyValueChoices = new String[][] {
};
// Kernel:
kernelProperties = new String[] {
EntropyCalculatorKernel.KERNEL_WIDTH_PROP_NAME,
@ -145,11 +134,7 @@ public class AutoAnalyserEntropy extends AutoAnalyser {
EntropyCalculatorKernel.NORMALISE_PROP_NAME +
" is set, then this is a number of standard deviations; " +
"otherwise it is an absolute value.",
"(boolean) whether to normalise <br/>the incoming time-series to mean 0, standard deviation 1, or not (default true, recommended)",
};
kernelPropertyValueChoices = new String[][] {
null,
{"true", "false"}
"(boolean) whether to normalise <br/>the incoming time-series to mean 0, standard deviation 1, or not (recommended)",
};
// KSG (Kraskov):
klProperties = new String[] {
@ -160,8 +145,6 @@ public class AutoAnalyserEntropy extends AutoAnalyser {
};
klPropertyDescriptions = new String[] {
};
klPropertyValueChoices = new String[][] {
};
}
@Override
@ -268,17 +251,14 @@ public class AutoAnalyserEntropy extends AutoAnalyser {
calcProperties.classSpecificPropertyNames = gaussianProperties;
calcProperties.classSpecificPropertiesFieldNames = gaussianPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = gaussianPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = gaussianPropertyValueChoices;
} else if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_KOZ_LEO)) {
calcProperties.classSpecificPropertyNames = klProperties;
calcProperties.classSpecificPropertiesFieldNames = klPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = klPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = klPropertyValueChoices;
} else if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_KERNEL)) {
calcProperties.classSpecificPropertyNames = kernelProperties;
calcProperties.classSpecificPropertiesFieldNames = kernelPropertiesFieldNames;
calcProperties.classSpecificPropertyDescriptions = kernelPropertyDescriptions;
calcProperties.classSpecificPropertyValueChoices = kernelPropertyValueChoices;
} else {
calcProperties = null;
throw new Exception("No recognised calculator selected: " +

49
java/source/infodynamics/demos/autoanalysis/AutoAnalyserMI.java
View File

@ -92,10 +92,6 @@ public class AutoAnalyserMI extends AutoAnalyserChannelCalculator
"Number of discrete states available for each variable (i.e. 2 for binary)",
"Time-lag from source to dest to consider MI across; must be >= 0 (0 for standard MI)",
};
discretePropertyValueChoices = new String[][] {
null,
null
};
// Continuous:
abstractContinuousClass = MutualInfoCalculatorMultiVariate.class;
@ -109,23 +105,12 @@ public class AutoAnalyserMI extends AutoAnalyserChannelCalculator
commonContPropertyDescriptions = new String[] {
"Time-lag from source to dest to consider MI across; must be >= 0 (0 for standard MI)"
};
commonContPropertyValueChoices = new String[][] {
null
};
// Gaussian properties:
gaussianProperties = new String[] {
MutualInfoCalculatorMultiVariateGaussian.PROP_BIAS_CORRECTION
};
gaussianPropertiesFieldNames = new String[] {
"PROP_BIAS_CORRECTION"
};
gaussianPropertyDescriptions = new String[] {
"Whether the analytically determined bias (as the mean of the<br/>" +
"surrogate distribution) will be subtracted from all" +
"calculated values. Default is false."
};
gaussianPropertyValueChoices = new String[][] {
{"true", "false"}
};
// Kernel:
kernelProperties = new String[] {
@ -147,11 +132,6 @@ public class AutoAnalyserMI extends AutoAnalyserChannelCalculator
"0 (default) means no dynamic exclusion window",
"(boolean) whether to normalise <br/>each incoming time-series to mean 0, standard deviation 1, or not (recommended)",
};
kernelPropertyValueChoices = new String[][] {
null,
null,
{"true", "false"}
};
// KSG (Kraskov):
kraskovProperties = new String[] {
MutualInfoCalculatorMultiVariateKraskov.PROP_NORMALISE,
@ -185,15 +165,7 @@ public class AutoAnalyserMI extends AutoAnalyserChannelCalculator
"(default, to indicate to use all available processors)",
"Whether to enable the GPU module (number of threads then has no bearing); boolean, default false"
};
kraskovPropertyValueChoices = new String[][] {
{"true", "false"},
null,
null,
null,
{"MAX_NORM", "EUCLIDEAN", "EUCLIDEAN_SQUARED"},
null,
{"true", "false"},
};
}
/**
@ -241,26 +213,11 @@ public class AutoAnalyserMI extends AutoAnalyserChannelCalculator
}
return new DiscreteCalcAndArguments(
new MutualInformationCalculatorDiscrete(base, base, timeDiff),
new MutualInformationCalculatorDiscrete(base, timeDiff),
base,
base + ", " + base + ", " + timeDiff);
base + ", " + timeDiff);
}
@Override
protected String pythonSetObsSuffix() {
String selectedCalcType = (String)
calcTypeComboBox.getSelectedItem();
if (selectedCalcType.equalsIgnoreCase(CALC_TYPE_DISCRETE) ||
selectedCalcType.equalsIgnoreCase(CALC_TYPE_BINNED)) {
return "";
} else {
// For the moment we could direct all calls to the 1D arrays version,
// but it is working fine with JPype 0.7; later
// when we have 2D inputs we should dynamically detect that and return "2D"
return "";
}
}
/**
* @param args
*/

Some files were not shown because too many files have changed in this diff Show More