From 379f4496f0e6705f0c318cbd2fee908b6f99b676 Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Tue, 3 Feb 2015 11:13:56 +0000 Subject: [PATCH] Adding Schreiber TE example for Python, with Kraskov estimator, mirroring the main demo here in Matlab/Octave --- .../heartRateBreathRateKraskov.sh | 2 + .../runHeartBreathRateKraskov.py | 131 ++++++++++++++++++ demos/python/readFloatsFile.py | 6 + demos/python/readIntsFile.py | 6 + 4 files changed, 145 insertions(+) create mode 100644 demos/python/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.py diff --git a/demos/java/SchreiberTransferEntropyExamples/heartRateBreathRateKraskov.sh b/demos/java/SchreiberTransferEntropyExamples/heartRateBreathRateKraskov.sh index ca5ef07..2d5be28 100755 --- a/demos/java/SchreiberTransferEntropyExamples/heartRateBreathRateKraskov.sh +++ b/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 + diff --git a/demos/python/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.py b/demos/python/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.py new file mode 100644 index 0000000..be582f4 --- /dev/null +++ b/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 . +## + +# 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 + diff --git a/demos/python/readFloatsFile.py b/demos/python/readFloatsFile.py index 3c0df03..64d2b9b 100755 --- a/demos/python/readFloatsFile.py +++ b/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 diff --git a/demos/python/readIntsFile.py b/demos/python/readIntsFile.py index 3eaf1a4..8b1d99a 100755 --- a/demos/python/readIntsFile.py +++ b/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