Adding Schreiber TE example for Python, with Kraskov estimator, mirroring the main demo here in Matlab/Octave

demos/java/SchreiberTransferEntropyExamples/heartRateBreathRateKraskov.sh

@@ -17,3 +17,5 @@ javac -classpath "../../infodynamics.jar" "infodynamics/demos/schreiberTransferE
 # Run the example:
 java -classpath ".:../../infodynamics.jar" infodynamics.demos.schreiberTransferEntropyExamples.HeartBreathRateKraskovRunner $1 $2 $3 $4
 
+cd SchreiberTransferEntropyExamples
+

demos/python/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.py

@@ -0,0 +1,131 @@
+##
+##  Java Information Dynamics Toolkit (JIDT)
+##  Copyright (C) 2015, 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/>.
+##
+
+# runHeartBreathRateKraskov.py kHistory lHistory knns numSurrogates
+#
+# runHeartBreathRateKraskov.py
+# Version 1.0
+# Joseph Lizier
+# 3/2/2015
+#
+# Used to explore information transfer in the heart rate / breath rate example of Schreiber --
+#  but estimates TE using Kraskov-Stoegbauer-Grassberger estimation.
+#
+#
+# Inputs
+# - kHistory - destination embedding length
+# - lHistory - source embedding length
+# - knns - a scalar specifying a single, or vector specifying a comma separated list, for values of K nearest neighbours to evaluate TE (Kraskov) with.
+# - numSurrogates - a scalar specifying the number of surrogates to evaluate TE from null distribution
+#
+# Run e.g. python runHeartBreathRateKraskov.py 2 2 1,2,3,4,5,6,7,8,9,10
+
+from jpype import *
+import sys
+import os
+import random
+import math
+import string
+import numpy
+# Import our readFloatsFile utility in the above directory:
+sys.path.append(os.path.relpath(".."))
+import readFloatsFile
+
+# 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)
+
+# Read in the command line arguments and assign default if required.
+# first argument in argv is the filename, so program arguments start from index 1.
+if (len(sys.argv) < 2):
+	kHistory = 1;
+else:
+	kHistory = int(sys.argv[1]);
+if (len(sys.argv) < 3):
+	lHistory = 1;
+else:
+	lHistory = int(sys.argv[2]);
+if (len(sys.argv) < 4):
+	knns = [4];
+else:
+	knnsStrings = sys.argv[3].split(",");
+	knns = [int(i) for i in knnsStrings]
+if (len(sys.argv) < 5):
+	numSurrogates = 0;
+else:
+	numSurrogates = int(sys.argv[4]);
+
+# Read in the data
+datafile = '../../data/SFI-heartRate_breathVol_bloodOx.txt'
+rawData = readFloatsFile.readFloatsFile(datafile)
+# As numpy array:
+data = numpy.array(rawData)
+# Heart rate is first column, and we restrict to the samples that Schreiber mentions (2350:3550)
+heart = data[2349:3550,0]; # Extracts what Matlab does with 2350:3550 argument there.
+# Chest vol is second column
+chestVol = data[2349:3550,1];
+# bloodOx = data[2349:3550,2];
+
+timeSteps = len(heart);
+	
+print("TE for heart rate <-> breath rate for Kraskov estimation with %d samples:" % timeSteps);
+
+# Using a KSG estimator for TE is the least biased way to run this:
+teCalcClass = JPackage("infodynamics.measures.continuous.kraskov").TransferEntropyCalculatorKraskov
+teCalc = teCalcClass();
+
+teHeartToBreath = [];
+teBreathToHeart = [];
+
+for knnIndex in range(len(knns)):
+	knn = knns[knnIndex];
+	# Compute a TE value for knn nearest neighbours
+				
+	# Perform calculation for heart -> breath (lag 1)
+	teCalc.initialise(kHistory,1,lHistory,1,1);
+	teCalc.setProperty("k", str(knn));
+	teCalc.setObservations(JArray(JDouble, 1)(heart),
+				JArray(JDouble, 1)(chestVol));
+	teHeartToBreath.append( teCalc.computeAverageLocalOfObservations() );
+	if (numSurrogates > 0):
+		teHeartToBreathNullDist = teCalc.computeSignificance(numSurrogates);
+		teHeartToBreathNullMean = teHeartToBreathNullDist.getMeanOfDistribution();
+		teHeartToBreathNullStd = teHeartToBreathNullDist.getStdOfDistribution();
+	
+	# Perform calculation for breath -> heart (lag 1)
+	teCalc.initialise(kHistory,1,lHistory,1,1);
+	teCalc.setProperty("k", str(knn));
+	teCalc.setObservations(JArray(JDouble, 1)(chestVol),
+				JArray(JDouble, 1)(heart));
+	teBreathToHeart.append( teCalc.computeAverageLocalOfObservations() );
+	if (numSurrogates > 0):
+		teBreathToHeartNullDist = teCalc.computeSignificance(numSurrogates);
+		teBreathToHeartNullMean = teBreathToHeartNullDist.getMeanOfDistribution();
+		teBreathToHeartNullStd = teBreathToHeartNullDist.getStdOfDistribution();
+	
+	print("TE(k=%d,l=%d,knn=%d): h->b = %.3f" % (kHistory, lHistory, knn, teHeartToBreath[knnIndex])), # , for no newline
+	if (numSurrogates > 0):
+		print(" (null = %.3f +/- %.3f)" % (teHeartToBreathNullMean, teHeartToBreathNullStd)),
+	print(", b->h = %.3f nats" % teBreathToHeart[knnIndex]),
+	if (numSurrogates > 0):
+		print("(null = %.3f +/- %.3f)" % (teBreathToHeartNullMean, teBreathToHeartNullStd)),
+	print
+	
+# Exercise: plot the results
+

demos/python/readFloatsFile.py

@@ -22,6 +22,12 @@ def readFloatsFile(filename):
  	    # Space separate numbers, one time step per line, each column is a variable
  	    array = []
 	    for line in f: # read all lines
+	    	if (line.startswith("%") or line.startswith("#")):
+	    		# Assume this is a comment line
+	    		continue;
+	    	if (len(line.split()) == 0):
+	    		# Line is empty
+	    		continue;
 		array.append([float(x) for x in line.split()])
 	    return array
     

demos/python/readIntsFile.py

@@ -22,6 +22,12 @@ def readIntsFile(filename):
  	    # Space separate numbers, one time step per line, each column is a variable
  	    array = []
 	    for line in f: # read all lines
+	    	if (line.startswith("%") or line.startswith("#")):
+	    		# Assume this is a comment line
+	    		continue;
+	    	if (len(line.split()) == 0):
+	    		# Line is empty
+	    		continue;
 		array.append([int(x) for x in line.split()])
 	    return array