Targeting Java version 8 in our builds now, in part to be compatible with latest Jpype in python

sensible parameter choices for SpikeTrainTETesting.py

README.md

@@ -19,8 +19,8 @@ JIDT is distributed under the [GNU GPL v3 license](http://www.gnu.org/licenses/g
 # Getting started
 
 1. [Download](../../wiki/Downloads) and [Installation](../../wiki/Installation) is very easy!
-    1. _Quick start_: take a `git clone` (then build via [AntScripts](../../wiki/AntScripts)) OR download the latest [v1.5 full distribution](http://lizier.me/joseph/software/jidt/download.php?file=infodynamics-dist-1.5.zip) (suitable for all platforms) and see the readme.txt file therein.
-1. [Documentation](../../wiki/Documentation) including: the paper describing JIDT at [arXiv:1408.3270](http://arxiv.org/abs/1408.3270) (distributed with the toolkit), a (beta) [Course](../../wiki/Course) including short video lectures and a shorter [Tutorial](../../wiki/Tutorial), and [Javadocs (v1.5 here)](http://lizier.me/joseph/software/jidt/javadocs/v1.5/);
+    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)
 
@@ -45,6 +45,10 @@ 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;

build.xml

@@ -5,7 +5,7 @@
 	</description>
 	
 	<!-- set global properties for this build -->
-	<property name="version" value="1.6"/>
+	<property name="version" value="1.6.1"/>
 	<property name="mainfilename" value="infodynamics"/>
 	<property name="jarplainname" value="${mainfilename}.jar" />
 	<property name="jarversiondistnamezip" value="${mainfilename}-jar-${version}.zip" />
@@ -43,17 +43,17 @@
 
 	<!-- Compile the java toolkit -->
 	<target name="compile" depends="init" description="compile the source">
-		<!-- Compile to Java 7 to provide compatibility for users with older JREs.
+		<!-- Compile to Java 8 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 7.
-		     Indeed, one gets the warning: "bootstrap class path not set in conjunction with -source 1.7"
+		     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"
 		     To fix this, one would use the bootstrap classpath to point our JDK to an rt.jar
-		     for Java 7. 
-		     At this stage, I'm sure I'm not using new library calls from Java 8+, so we can
+		     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 8 libraries so as not to cause
+		     like this. I'll endeavour not to use JDK 9+ libraries so as not to cause
 		     any issues here ... -->
-		<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.7" source="1.7" encoding="UTF8">
+		<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.8" source="1.8" encoding="UTF8">
 			<classpath refid="apache-classpath"/>
 		</javac>
 

demos/AutoAnalyser/launchAutoAnalyser.py

@@ -0,0 +1,11 @@
+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)
+

demos/octave/FlockingAnalysis/NetLogoExample/loadProperties.m

@@ -43,6 +43,8 @@ properties.files = 'positions%s.txt';
 % 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;
 
@@ -102,7 +104,10 @@ 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 Matlan in order to call autocorr. Set to false if you don't have that to avoid this being called.
+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;

demos/octave/FlockingAnalysis/generateObservations.m

@@ -145,7 +145,7 @@ for fileIndex = 1:length(files)
 	if (acfTime > maxAcfTime)
 		maxAcfTime = acfTime;
 	end
-	fprintf('ACF time of %d from heading %d and speed %d\n', acfTime, maxAcfHeadingsDecayTime, maxAcfSpeedDecayTime);
+	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.

demos/octave/FlockingAnalysis/loadBasic2d.m

@@ -1,9 +1,11 @@
 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 timestamp in column 1,
-%  then fish 1's X and Y coordinates in columns 2 and 3, 
-%  then fish 2's X and Y coordinates in columns 4 and 5, 
+% 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.
 
 %%
@@ -40,7 +42,14 @@ data = load(dataFileName);
 
 %%% FORMAT DATA IN A MORE CONVENIENT WAY %%%
 
-numFish = (size(data,2)-1) ./ 2; % number of fish (we know it from the raw data)
+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);
 
@@ -52,8 +61,8 @@ posY = nan(numCycles,numFish);
 % fill the position tables
 for f = 1 : numFish % cycle over all fish
     % copy into new variables
-    posX(:,f) = data(:,2+(f-1)*2);
-    posY(:,f) = data(:,3+(f-1)*2);
+    posX(:,f) = data(:,startCol+(f-1)*2);
+    posY(:,f) = data(:,startCol+1+(f-1)*2);
 end
 
 % display to check

demos/octave/FlockingAnalysis/loadxls2d.m

@@ -1,4 +1,4 @@
-function [x,y] = loadxls3d(filename, properties)
+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.:

demos/octave/FlockingAnalysis/runAnalysis.m

@@ -81,11 +81,11 @@ if (isfield(properties, 'kRange') || isfield(properties, 'tauRange'))
 	maxAIStau = properties.tauRange(1);
 	aisForKAndTau = zeros(length(properties.kRange), length(properties.tauRange));
 	kIndex = 0;
-	tauIndex = 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))
@@ -118,7 +118,7 @@ if (isfield(properties, 'kRange') || isfield(properties, 'tauRange'))
 	properties.tau = maxAIStau;
 	properties.ais = maxAIS;
 	properties.destSamplesOnly = false;
-	fprintf('*** Optmised k=%d and tau=%d (giving AIS=%.4f)\n', properties.k, ...
+	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 :
@@ -140,6 +140,15 @@ else
 	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

demos/octave/FlockingAnalysis/watchRawMovement.m

@@ -31,6 +31,10 @@ 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);

demos/python/EffectiveNetworkInference/spiking/net_inf.py

@@ -49,23 +49,23 @@ LOG_FILE_NAME = "logs/" + net_type_name + "_" + num_spikes_string + "_" + repeat
 log = open(LOG_FILE_NAME, "w")
 sys.stdout = log
 
-def prepare_conditional_trains(calc_object, cond_set):
+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])
-                        teCalc.appendConditionalIntervals(JArray(JInt, 1)(cond_set[key]))
+                        calc_object.appendConditionalIntervals(JArray(JInt, 1)(cond_set[key]))
         return cond_trains
 
-def set_target_embeddings(embedding_list):
+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])
-                teCalc.setProperty("DEST_PAST_INTERVALS", embedding_string)
+                calc_object.setProperty("DEST_PAST_INTERVALS", embedding_string)
         else:
-                teCalc.setProperty("DEST_PAST_INTERVALS", "")
+                calc_object.setProperty("DEST_PAST_INTERVALS", "")
 
 
 target_index = int(target_index_string)                
@@ -100,7 +100,7 @@ next_target_interval = 2
 still_significant = True
 print("**** Determining target embedding set ****\n")
 while still_significant:
-        set_target_embeddings(target_embedding_set)
+        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]))
@@ -136,7 +136,7 @@ while still_significant:
         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)
+        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))
@@ -218,7 +218,7 @@ while not everything_significant:
                 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)
+                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))

demos/python/SpikingTE/SpikeTrainTETesting.py

@@ -81,8 +81,10 @@ teCalc.setProperty("knns", "4")
 print("Independent Poisson Processes")
 teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
 teCalc.setProperty("SOURCE_PAST_INTERVALS", "1,2")
-teCalc.setProperty("DO_JITTERED_SAMPLING", "true")
-teCalc.setProperty("JITTERED_SAMPLING_NOISE_LEVEL", "0")
+# 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") 
@@ -112,8 +114,7 @@ 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", "true")
-teCalc.setProperty("JITTERED_SAMPLING_NOISE_LEVEL", "0")
+teCalc.setProperty("DO_JITTERED_SAMPLING", "false")
 #teCalc.setProperty("NORM_TYPE", "MAX_NORM")
 
 results_noisy_zero = np.zeros(NUM_REPS)
@@ -145,8 +146,7 @@ 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", "true")
-teCalc.setProperty("JITTERED_SAMPLING_NOISE_LEVEL", "0")
+teCalc.setProperty("DO_JITTERED_SAMPLING", "false")
 #teCalc.setProperty("NORM_TYPE", "MAX_NORM") 
 
 results_noisy_non_zero = np.zeros(NUM_REPS)
@@ -175,8 +175,7 @@ teCalc = teCalcClass()
 teCalc.setProperty("knns", "4")
 teCalc.setProperty("DEST_PAST_INTERVALS", "1,2")
 teCalc.setProperty("SOURCE_PAST_INTERVALS", "1")
-teCalc.setProperty("DO_JITTERED_SAMPLING", "true")
-teCalc.setProperty("JITTERED_SAMPLING_NOISE_LEVEL", "0")
+teCalc.setProperty("DO_JITTERED_SAMPLING", "false")
 #teCalc.setProperty("NUM_SAMPLES_MULTIPLIER", "1")
 #teCalc.setProperty("NORM_TYPE", "MAX_NORM") 
 

demos/python/create_env.bat

@@ -0,0 +1,87 @@
+@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
+

demos/python/create_env.sh

@@ -0,0 +1,84 @@
+#!/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
+

demos/python/example1TeBinaryData.py

@@ -52,8 +52,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=numpy.int)
-destNumpy = numpy.array(destArray, dtype=numpy.int)
+sourceNumpy = numpy.array(sourceArray, dtype=int)
+destNumpy = numpy.array(destArray, dtype=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:

demos/python/start_env.sh

@@ -0,0 +1,42 @@
+#!/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
+

demos/python/start_jupyter_env.sh

@@ -0,0 +1,41 @@
+#!/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
+

java/source/infodynamics/demos/autoanalysis/AutoAnalyser.java

@@ -808,10 +808,11 @@ public abstract class AutoAnalyser extends JFrame
 		} else {
 			pythonCode.append("import readFloatsFile\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");
+		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");
 		// 3. Matlab:
 		StringBuffer matlabCode = new StringBuffer();
 		matlabCode.append("% Add JIDT jar library to the path, and disable warnings that it's already there:\n");
@@ -1240,7 +1241,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 +

java/source/infodynamics/measures/continuous/ActiveInfoStorageCalculatorViaMutualInfo.java

@@ -96,6 +96,12 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	 */
 	protected Vector<boolean[]> vectorOfValidityOfObservations;
 
+	/**
+	 * Store the time index at which we were asked to start taking time-series
+	 *  observations for each  observation set.
+	 */
+	protected Vector<Integer> vectorOfOffsetsInTimeSeries;
+
 	/**
 	 * Property name for the auto-embedding method. Defaults to {@link #AUTO_EMBED_METHOD_NONE}.
 	 * Other valid values are {@link #AUTO_EMBED_METHOD_RAGWITZ} or
@@ -242,6 +248,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 		this.tau = tau;
 		vectorOfObservationTimeSeries = null;
 		vectorOfValidityOfObservations = null;
+		vectorOfOffsetsInTimeSeries = null;
 	}
 
 	/**
@@ -258,7 +265,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	 * 	   and embedding delay ({@link #TAU_PROP_NAME}). Default is {@link #AUTO_EMBED_METHOD_NONE} meaning
 	 * 	   values are set manually; other accepted values include: {@link #AUTO_EMBED_METHOD_RAGWITZ} for use
 	 * 		of the Ragwitz criteria and {@link #AUTO_EMBED_METHOD_MAX_CORR_AIS} for using
-	 * 		the maz bias-corrected AIS criteria (both searching up to {@link #PROP_K_SEARCH_MAX} and 
+	 * 		the max bias-corrected AIS criteria (both searching up to {@link #PROP_K_SEARCH_MAX} and 
 	 * 		{@link #PROP_TAU_SEARCH_MAX}, as outlined by Garland et al. in the references list above).
 	 * 		Use of any value other than {@link #AUTO_EMBED_METHOD_NONE}
 	 * 		will lead to any previous settings for k and tau (via e.g. {@link #initialise(int, int)} or
@@ -366,6 +373,18 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 		finaliseAddObservations();
 	}
 
+	/* (non-Javadoc)
+	 * @see infodynamics.measures.continuous.ActiveInfoStorageCalculator#setObservations(double[], boolean[])
+	 */
+	@Override
+	public void setObservations(double[] observations, boolean[] valid)
+			throws Exception {
+		startAddObservations();
+		// Add these observations and the indication of their validity
+		addObservations(observations, valid);
+		finaliseAddObservations();
+	}
+
 	/* (non-Javadoc)
 	 * @see infodynamics.measures.continuous.ActiveInfoStorageCalculator#startAddObservations()
 	 */
@@ -373,6 +392,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	public void startAddObservations() {
 		vectorOfObservationTimeSeries = new Vector<double[]>();
 		vectorOfValidityOfObservations = new Vector<boolean[]>();
+		vectorOfOffsetsInTimeSeries = new Vector<Integer>();
 	}
 
 	/* (non-Javadoc)
@@ -380,9 +400,15 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	 */
 	@Override
 	public void addObservations(double[] observations) throws Exception {
+		addObservationsParsed(observations, 0);
+	}
+	
+	protected void addObservationsParsed(double[] observations, int startTimeStep)
+			throws Exception {
 		// Store these observations in our vector for now
 		vectorOfObservationTimeSeries.add(observations);
 		vectorOfValidityOfObservations.add(null); // All observations were valid
+		vectorOfOffsetsInTimeSeries.add(startTimeStep);
 	}
 	
 	/* (non-Javadoc)
@@ -393,6 +419,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 		// Store these observations in our vector for now
 		vectorOfObservationTimeSeries.add(observations);
 		vectorOfValidityOfObservations.add(valid); // All observations were valid
+		vectorOfOffsetsInTimeSeries.add(0); // no offset here
 	}
 
 	/**
@@ -403,10 +430,13 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	 * @param k_in_use k embedding dimension to use
 	 * @param tau_in_use tau embedding delay to use
 	 * @param observations time series of observations
+	 * @param observationSetIndex which observation set these samples came from
+	 * @param offsetInOriginalTimeSeries offset of the samples in their original time series
 	 * @throws Exception
 	 */
 	protected void addObservationsWithGivenParams(MutualInfoCalculatorMultiVariate miCalc_in_use,
-			int k_in_use, int tau_in_use, double[] observations) throws Exception {
+			int k_in_use, int tau_in_use, double[] observations,
+			int observationSetIndex, int offsetInOriginalTimeSeries) throws Exception {
 		if (observations.length - (k_in_use-1)*tau_in_use - 1 <= 0) {
 			// There are no observations to add here
 			// Don't throw an exception, do nothing since more observations
@@ -419,7 +449,8 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 		double[][] currentDestNextVectors =
 				MatrixUtils.makeDelayEmbeddingVector(observations, 1, (k_in_use-1)*tau_in_use + 1,
 						observations.length - (k_in_use-1)*tau_in_use - 1);
-		miCalc_in_use.addObservations(currentDestPastVectors, currentDestNextVectors);
+		miCalc_in_use.addObservationsTrackObservationIDs(currentDestPastVectors, currentDestNextVectors,
+				observationSetIndex, offsetInOriginalTimeSeries + (k_in_use-1)*tau_in_use + 1);
 	}
 	
 	/**
@@ -436,10 +467,12 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	 *  whether the entry in observations at that index is valid; we only take vectors
 	 *  as samples to add to the observation set where all points in the time series
 	 *  (even between points in the embedded k-vector with embedding delays) are valid.
+	 * @param observationSetIndex which observation set these samples came from
 	 * @throws Exception
 	 */
 	protected void addObservationsWithGivenParams(MutualInfoCalculatorMultiVariate miCalc_in_use,
-			int k_in_use, int tau_in_use, double[] observations, boolean[] valid) throws Exception {
+			int k_in_use, int tau_in_use, double[] observations, boolean[] valid,
+			int observationSetIndex) throws Exception {
 		
 		// compute the start and end times using our determined embedding parameters:
 		Vector<int[]> startAndEndTimePairs = computeStartAndEndTimePairs(k_in_use, tau_in_use, valid);
@@ -448,7 +481,8 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 			int startTime = timePair[0];
 			int endTime = timePair[1];
 			addObservationsWithGivenParams(miCalc_in_use, k_in_use, tau_in_use,
-					MatrixUtils.select(observations, startTime, endTime - startTime + 1));
+					MatrixUtils.select(observations, startTime, endTime - startTime + 1),
+					observationSetIndex, startTime);
 		}
 	}
 	
@@ -458,7 +492,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	@Override
 	public void addObservations(double[] observations, int startTime,
 			int numTimeSteps) throws Exception {
-		addObservations(MatrixUtils.select(observations, startTime, numTimeSteps));
+		addObservationsParsed(MatrixUtils.select(observations, startTime, numTimeSteps), startTime);
 	}
 
 	/**
@@ -616,6 +650,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 		
 		vectorOfObservationTimeSeries = null; // No longer required
 		vectorOfValidityOfObservations = null;
+		vectorOfOffsetsInTimeSeries = null;
 	}
 
 	/**
@@ -639,34 +674,25 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 		miCalc_in_use.startAddObservations();
 		// Send all of the observations through:
 		Iterator<boolean[]> validityIterator = vectorOfValidityOfObservations.iterator();
+		Iterator<Integer> offsetsInTimeSeriesIterator = vectorOfOffsetsInTimeSeries.iterator();
+		int setNum = 0;
 		for (double[] observations : vectorOfObservationTimeSeries) {
 			boolean[] validity = validityIterator.next();
 			if (validity == null) {
 				// Add the whole time-series
 				addObservationsWithGivenParams(miCalc_in_use, k_in_use,
-						tau_in_use, observations);
+						tau_in_use, observations,
+						setNum, offsetsInTimeSeriesIterator.next());
 			} else {
 				addObservationsWithGivenParams(miCalc_in_use, k_in_use,
-						tau_in_use, observations, validity);
+						tau_in_use, observations, validity, setNum);
 			}
+			setNum++;
 		}
 		// TODO do we need to throw an exception if there are no observations to add?
 		miCalc_in_use.finaliseAddObservations();
 	}
 	
-	/* (non-Javadoc)
-	 * @see infodynamics.measures.continuous.ActiveInfoStorageCalculator#setObservations(double[], boolean[])
-	 */
-	@Override
-	public void setObservations(double[] observations, boolean[] valid)
-			throws Exception {
-		startAddObservations();
-		// Add these observations and the indication of their validity
-		vectorOfObservationTimeSeries.add(observations);
-		vectorOfValidityOfObservations.add(valid);
-		finaliseAddObservations();
-	}
-
 	/**
 	 * Compute a vector of start and end pairs of time points, between which we have
 	 *  valid series of observations.
@@ -855,4 +881,20 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
 	public int getNumObservations() throws Exception {
 		return miCalc.getNumObservations();
 	}
+
+	/**
+	 * Retrieve an array indicating which observation set each sample came from
+	 * @return
+	 */
+	public int[] getObservationSetIndices() {
+		return miCalc.getObservationSetIndices();
+	}
+
+	/**
+	 * Retrieve an array indicating which time index within its observation set that sample came from
+	 * @return
+	 */
+	public int[] getObservationTimePoints() {
+		return miCalc.getObservationTimePoints();
+	}
 }

java/source/infodynamics/measures/continuous/ConditionalMutualInfoCalculatorMultiVariate.java

@@ -92,6 +92,16 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 	 *  where it is 1e-8, matching the MILCA toolkit)
 	 */
 	public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD";
+	/**
+	 * Property name for the seed for the random number generator for noise to be
+	 *  added to the data (default is no seed)
+	 */
+	public static final String PROP_NOISE_SEED = "NOISE_SEED";
+	/**
+	 * Property value to indicate no seed for the random number generator for noise to be
+	 *  added to the data
+	 */
+	public static final String NOISE_NO_SEED_VALUE = "NONE";
 
 	/**
 	 * Initialise the calculator for (re-)use, clearing PDFs,
@@ -129,17 +139,8 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 			double[][] cond) throws Exception;
 
 	/**
-	 * Sets a single series from which to compute the PDF.
-	 * Cannot be called in conjunction with
-	 * {@link #startAddObservations()} / {@link #addObservations(double[][], double[][], double[][])} or
-	 * {@link #addObservations(double[][], double[][], double[][], int, int)} /
-	 * {@link #finaliseAddObservations()}.</p>
-	 * 
-	 * <p>The supplied series may be (multivariate) time-series or 
-	 * simply a set of separate observations without a time interpretation.
-	 * 
-	 * <p>This method can only be used where dimensions of all
-	 * variables have been set to 1.</p>
+	 * As per {@link #setObservations(double[][], double[][], double[][])}
+	 * but where dimensions of all variables have been set to 1 
 	 * 
 	 * @param var1 univariate observations for variable 1
 	 * @param var2 univariate observations for variable 2
@@ -151,6 +152,90 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 	public void setObservations(double[] var1, double[] var2,
 			double[] cond) throws Exception;
 
+	/**
+	 * As per {@link #setObservations(double[][], double[][], double[][])}
+	 * but where dimensions of var2 and cond are 1 
+	 * 
+	 * @param var1 observations for variable 1
+	 * @param var2 univariate observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond univariate observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void setObservations(double[][] var1, double[] var2,
+			double[] cond) throws Exception;
+
+	/**
+	 * As per {@link #setObservations(double[][], double[][], double[][])}
+	 * but where dimensions of var1 and cond are 1 
+	 * 
+	 * @param var1 univariate observations for variable 1
+	 * @param var2 observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond univariate observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void setObservations(double[] var1, double[][] var2,
+			double[] cond) throws Exception;
+
+	/**
+	 * As per {@link #setObservations(double[][], double[][], double[][])}
+	 * but where dimensions of var1 and var2 are 1 
+	 * 
+	 * @param var1 univariate observations for variable 1
+	 * @param var2 univariate observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void setObservations(double[] var1, double[] var2,
+			double[][] cond) throws Exception;
+
+	/**
+	 * As per {@link #setObservations(double[][], double[][], double[][])}
+	 * but where dimension of cond is 1 
+	 * 
+	 * @param var1 observations for variable 1
+	 * @param var2 observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond univariate observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void setObservations(double[][] var1, double[][] var2,
+			double[] cond) throws Exception;
+
+	/**
+	 * As per {@link #setObservations(double[][], double[][], double[][])}
+	 * but where dimension of var2 is 1 
+	 * 
+	 * @param var1 observations for variable 1
+	 * @param var2 univariate observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void setObservations(double[][] var1, double[] var2,
+			double[][] cond) throws Exception;
+
+	/**
+	 * As per {@link #setObservations(double[][], double[][], double[][])}
+	 * but where dimension of var1 is 1 
+	 * 
+	 * @param var1 univariate observations for variable 1
+	 * @param var2 observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void setObservations(double[] var1, double[][] var2,
+			double[][] cond) throws Exception;
+
 	/**
 	 * Sets a single series from which to compute the PDF,
 	 *  where all the various observations are valid.
@@ -248,18 +333,8 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 			double[][] cond) throws Exception;
 	
 	/**
-	 * <p>Adds a new set of observations to update the PDFs with - is
-	 * intended to be called multiple times.
-	 * Must be called after {@link #startAddObservations()}; call
-	 * {@link #finaliseAddObservations()} once all observations have
-	 * been supplied.</p>
-	 * 
-	 * <p>Note that the arrays must not be over-written by the user
-	 *  until after finaliseAddObservations() has been called
-	 *  (they are not copied by this method necessarily, but the method
-	 *  may simply hold a pointer to them).</p>
-	 * 
-	 * <p>This method can only be used where dimensions of all
+	 * <p>As per {@link #addObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of all
 	 * variables have been set to 1.</p>
 	 * 
 	 * @param var1 univariate observations for variable 1
@@ -272,6 +347,96 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 	public void addObservations(double[] var1, double[] var2,
 			double[] cond) throws Exception;
 	
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of
+	 * var2 and cond have been set to 1.</p>
+	 * 
+	 * @param var1 observations for variable 1
+	 * @param var2 univariate observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond univariate observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void addObservations(double[][] var1, double[] var2,
+			double[] cond) throws Exception;
+
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of
+	 * var1 and cond have been set to 1.</p>
+	 * 
+	 * @param var1 univariate observations for variable 1
+	 * @param var2 observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond univariate observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void addObservations(double[] var1, double[][] var2,
+			double[] cond) throws Exception;
+
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of
+	 * var1 and var2 have been set to 1.</p>
+	 * 
+	 * @param var1 univariate observations for variable 1
+	 * @param var2 univariate observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void addObservations(double[] var1, double[] var2,
+			double[][] cond) throws Exception;
+
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of
+	 * var1 have been set to 1.</p>
+	 * 
+	 * @param var1 univariate observations for variable 1
+	 * @param var2 observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void addObservations(double[] var1, double[][] var2,
+			double[][] cond) throws Exception;
+
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of
+	 * var2 have been set to 1.</p>
+	 * 
+	 * @param var1 observations for variable 1
+	 * @param var2 univariate observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void addObservations(double[][] var1, double[] var2,
+			double[][] cond) throws Exception;
+
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of
+	 * cond have been set to 1.</p>
+	 * 
+	 * @param var1 observations for variable 1
+	 * @param var2 observations for variable 2
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond univariate observations for the conditional
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @throws Exception
+	 */
+	public void addObservations(double[][] var1, double[][] var2,
+			double[] cond) throws Exception;
+
 	/**
 	 * <p>Adds a new sub-series of observations to update the PDFs with - is
 	 * intended to be called multiple times.
@@ -300,6 +465,28 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 			double[][] cond,
 			int startTime, int numTimeSteps) throws Exception;
 
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][], duoble[][])};
+     *  but also includes parameters to track which observation set
+     *  the samples came from. Intended to only be used by other 
+     *  estimator classes here and not by users directly.</p>
+	 * 
+	 * @param var1 multivariate observations for variable 1
+	 *  (first index is time or observation index, second is variable number)
+	 * @param var2 multivariate observations for variable 2
+	 *  (first index is time or observation index, second is variable number)
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param cond multivariate observations for the conditional
+	 *  (first index is time or observation index, second is variable number)
+	 *  Length must match <code>var1</code>, and their indices must correspond.
+	 * @param observationSetIndexToUse which set of observations these came fmor
+	 * @param startTimeIndex which was the first time index of these
+	 *  samples within that observation set.
+	 * @throws Exception
+	 */
+	public void addObservationsTrackObservationIDs(double[][] var1, double[][] var2, double[][] cond,
+			int observationSetIndexToUse, int startTimeIndex) throws Exception;
+
 	/**
 	 * Signal that the observations are now all added, PDFs can now be constructed.
 	 * 
@@ -436,6 +623,38 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 	public double[] computeLocalUsingPreviousObservations(double states1[][], double states2[][], double[][] condStates)
 		throws Exception;
 
+	/**
+	 * <p>As per {@link #computeLocalUsingPreviousObservations(double[][], double[][], double[][])}
+	 * but can only be used where dimensions of
+	 * states1 and states2 have been set to 1.</p>
+	 * 
+	 * @param states1 series of univariate observations for variable 1
+	 * @param states2 series of univariate observations for variable 2
+	 *  Length must match <code>states1</code>, and their indices must correspond.
+	 * @param condStates series of multivariate observations for the conditional
+	 *  (first index is time or observation index, second is variable number).
+	 *  Length must match <code>states1</code>, and their indices must correspond.
+	 * @return the series of local conditional MI values.
+	 * @throws Exception
+	 */
+	public double[] computeLocalUsingPreviousObservations(double states1[], double states2[], double[][] condStates)
+		throws Exception;
+
+	/**
+	 * <p>As per {@link #computeLocalUsingPreviousObservations(double[][], double[][], double[][])}
+	 * but can only be used all dimensions have been set to 1.</p>
+	 * 
+	 * @param states1 series of univariate observations for variable 1
+	 * @param states2 series of univariate observations for variable 2.
+	 *  Length must match <code>states1</code>, and their indices must correspond.
+	 * @param condStates series of univariate observations for the conditional.
+	 *  Length must match <code>states1</code>, and their indices must correspond.
+	 * @return the series of local conditional MI values.
+	 * @throws Exception
+	 */
+	public double[] computeLocalUsingPreviousObservations(double states1[], double states2[], double[] condStates)
+		throws Exception;
+
 	/**
 	 * @throws Exception if the implementing class computes MI without
 	 * explicit observations (e.g. see
@@ -454,4 +673,17 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
 	 */
 	public boolean getAddedMoreThanOneObservationSet();
 
+	/**
+	 * Retrieve an array indicating which observation set each sample came from
+	 * 
+	 * @return array of integers
+	 */
+	public int[] getObservationSetIndices();
+
+	/**
+	 * Retrieve an array indicating which time index within its observation set that sample came from
+	 * 
+	 * @return array of integers
+	 */
+	public int[] getObservationTimePoints();
 }

java/source/infodynamics/measures/continuous/ConditionalMutualInfoMultiVariateCommon.java

@@ -22,6 +22,7 @@ import infodynamics.utils.EmpiricalMeasurementDistribution;
 import infodynamics.utils.MatrixUtils;
 import infodynamics.utils.RandomGenerator;
 
+import java.util.Arrays;
 import java.util.Iterator;
 import java.util.Random;
 import java.util.Vector;
@@ -82,6 +83,16 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	 */
 	protected double[][] condObservations;
 
+	/**
+	 * Track which observation set each sample came from
+	 */
+	protected int[] observationSetIndices;
+
+	/**
+	 * Track which sample index within an observation set that each sample came from
+	 */
+	protected int[] observationTimePoints;
+	
 	/**
 	 * Total number of observations supplied.
 	 * Only valid after {@link #finaliseAddObservations()} is called.
@@ -119,6 +130,21 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	 * {@link #addObservations(double[][], double[][], double[][])} etc
 	 */
 	protected Vector<double[][]> vectorOfCondObservations;
+	/**
+	 * Tracks separate (time-series) observation sets
+	 *  we are taking samples from
+	 */
+	protected int observationSetIndex = 0;
+	/**
+	 * Storage for which observation set each
+	 *  block of samples comes from
+	 */
+	protected Vector<Integer> vectorOfObservationSetIndices;
+	/**
+	 * Storage for start time point for the observation
+	 *  set within its block of samples
+	 */
+	protected Vector<Integer> vectorOfObservationStartTimePoints;
 	
 	/**
 	 * Whether the user has added more than one disjoint observation set
@@ -140,6 +166,14 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	 *  and 1e-8 is used to match MILCA toolkit)
 	 */
 	protected double noiseLevel = (double) 0;
+	/**
+	 * Has the user set a seed for the random noise
+	 */
+	protected boolean noiseSeedSet = false;
+	/**
+	 * Seed that the user set for the random noise
+	 */
+	protected long noiseSeed = 0;
 
 	/**
 	 * Cache for the means of each dimension in variable 1, in case we need to normalise 
@@ -188,6 +222,11 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		var1Observations = null;
 		var2Observations = null;
 		condObservations = null;
+		observationSetIndices = null;
+		observationTimePoints = null;
+		observationSetIndex = 0;
+		vectorOfObservationSetIndices = null;
+		vectorOfObservationStartTimePoints = null;
 		addedMoreThanOneObservationSet = false;
 		var1Means = null;
 		var1Stds = null;
@@ -219,6 +258,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	 *      (Default is 0, except for KSG estimators where it is recommended by Kraskov
 	 *      and so they use 1e-8 to match the MILCA toolkit, although that adds in
 	 *      a random amount of noise in [0,noiseLevel) ).</li>
+	 *  <li>{@link #PROP_NOISE_SEED} -- a long value seed for the random noise generator or
+	 *      the string {@link ConditionalMutualInfoCalculatorMultiVariate#NOISE_NO_SEED_VALUE} for no seed (default)</li>
 	 * </ul>
 	 * 
 	 * <p>Unknown property values are ignored.</p>
@@ -229,6 +270,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	 */
 	@Override
 	public void setProperty(String propertyName, String propertyValue) {
+
+		boolean propertySet = true;
 		if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) {
 			normalise = Boolean.parseBoolean(propertyValue);
 		} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
@@ -240,6 +283,20 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 				addNoise = true;
 				noiseLevel = Double.parseDouble(propertyValue);
 			}
+	    } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) {
+	    	if (propertyValue.equals(NOISE_NO_SEED_VALUE)) {
+	    		noiseSeedSet = false;
+	    	} else {
+	    		noiseSeedSet = true;
+	    		noiseSeed = Long.parseLong(propertyValue);
+	    	}
+		} else {
+			// No property was set here
+			propertySet = false;
+		}
+		if (debug && propertySet) {
+			System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
+					" to " + propertyValue);
 		}
 	}
 
@@ -249,6 +306,12 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 			return Boolean.toString(normalise);
 		} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
 			return Double.toString(noiseLevel);
+	    } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) {
+	    	if (noiseSeedSet) {
+	    		return Long.toString(noiseSeed);
+	    	} else {
+	    		return NOISE_NO_SEED_VALUE;
+	    	}
 		} else {
 			// No property matches for this class
 			return null;
@@ -256,8 +319,7 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	}
 
 	@Override
-	public void setObservations(double[][] var1, double[][] var2,
-			double[][] cond) throws Exception {
+	public void setObservations(double[][] var1, double[][] var2, double[][] cond) throws Exception {
 		startAddObservations();
 		addObservations(var1, var2, cond);
 		finaliseAddObservations();
@@ -266,7 +328,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 
 	/**
 	 * A non-overloaded method signature for setObservations with 2D arguments, as there have been
-	 *  some problems calling overloaded versions of setObservations from jpype. 
+	 *  some problems calling overloaded versions of setObservations from python jpype.
+	 *  Resolved if one follows the AutoAnalyser generated code, but left for back compatibility.
 	 * 
 	 * @param var1
 	 * @param var2
@@ -278,8 +341,55 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	}
 
 	@Override
-	public void setObservations(double[] var1, double[] var2,
-			double[] cond) throws Exception {
+	public void setObservations(double[] var1, double[] var2, double[] cond) throws Exception {
+		startAddObservations();
+		addObservations(var1, var2, cond);
+		finaliseAddObservations();
+		addedMoreThanOneObservationSet = false;
+	}
+
+	@Override
+	public void setObservations(double[][] var1, double[] var2, double[] cond) throws Exception {
+		startAddObservations();
+		addObservations(var1, var2, cond);
+		finaliseAddObservations();
+		addedMoreThanOneObservationSet = false;
+	}
+
+	@Override
+	public void setObservations(double[] var1, double[][] var2, double[] cond) throws Exception {
+		startAddObservations();
+		addObservations(var1, var2, cond);
+		finaliseAddObservations();
+		addedMoreThanOneObservationSet = false;
+	}
+
+	@Override
+	public void setObservations(double[] var1, double[] var2, double[][] cond) throws Exception {
+		startAddObservations();
+		addObservations(var1, var2, cond);
+		finaliseAddObservations();
+		addedMoreThanOneObservationSet = false;
+	}
+
+	@Override
+	public void setObservations(double[] var1, double[][] var2, double[][] cond) throws Exception {
+		startAddObservations();
+		addObservations(var1, var2, cond);
+		finaliseAddObservations();
+		addedMoreThanOneObservationSet = false;
+	}
+
+	@Override
+	public void setObservations(double[][] var1, double[] var2, double[][] cond) throws Exception {
+		startAddObservations();
+		addObservations(var1, var2, cond);
+		finaliseAddObservations();
+		addedMoreThanOneObservationSet = false;
+	}
+
+	@Override
+	public void setObservations(double[][] var1, double[][] var2, double[] cond) throws Exception {
 		startAddObservations();
 		addObservations(var1, var2, cond);
 		finaliseAddObservations();
@@ -288,7 +398,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 
 	/**
 	 * A non-overloaded method signature for setObservations with 1D arguments, as there have been
-	 *  some problems calling overloaded versions of setObservations from jpype. 
+	 *  some problems calling overloaded versions of setObservations from python jpype.
+	 *  Resolved if one follows the AutoAnalyser generated code, but left for back compatibility.
 	 * 
 	 * @param var1
 	 * @param var2
@@ -299,16 +410,47 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		setObservations(var1, var, cond);
 	}
 	
+	@Override
+	public void setObservations(double[][] var1, double[][] var2,
+			double[][] cond,
+			boolean[] var1Valid, boolean[] var2Valid,
+			boolean[] condValid) throws Exception {
+		
+		startAddObservations();
+		addObservations(var1, var2, cond, var1Valid, var2Valid, condValid);
+		finaliseAddObservations();
+	}
+
+	@Override
+	public void setObservations(double[][] var1, double[][] var2,
+			double[][] cond,
+			boolean[][] var1Valid, boolean[][] var2Valid,
+			boolean[][] condValid) throws Exception {
+
+		startAddObservations();
+		addObservations(var1, var2, cond, var1Valid, var2Valid, condValid);
+		finaliseAddObservations();
+	}
+
 	@Override
 	public void startAddObservations() {
 		vectorOfVar1Observations = new Vector<double[][]>();
 		vectorOfVar2Observations = new Vector<double[][]>();
 		vectorOfCondObservations = new Vector<double[][]>();
+		vectorOfObservationSetIndices = new Vector<Integer>();
+		vectorOfObservationStartTimePoints = new Vector<Integer>();
 	}
 	
 	@Override
 	public void addObservations(double[][] var1, double[][] var2,
 			double[][] cond) throws Exception {
+		// Use the current observationSetIndex and increment for next use:
+		addObservationsTrackObservationIDs(var1, var2, cond, observationSetIndex++, 0);
+	}
+	
+	@Override
+	public void addObservationsTrackObservationIDs(double[][] var1, double[][] var2,
+			double[][] cond, int observationSetIndexToUse, int startTimeIndex) throws Exception {
 		if (vectorOfVar1Observations == null) {
 			// startAddObservations was not called first
 			throw new RuntimeException("User did not call startAddObservations before addObservations");
@@ -333,6 +475,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		vectorOfVar1Observations.add(var1);
 		vectorOfVar2Observations.add(var2);
 		vectorOfCondObservations.add(cond);
+		vectorOfObservationSetIndices.add(observationSetIndexToUse);
+		vectorOfObservationStartTimePoints.add(startTimeIndex);
 		if (vectorOfVar1Observations.size() > 1) {
 			addedMoreThanOneObservationSet = true;
 		}
@@ -364,6 +508,7 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		}
 		double[][] reshapedConditional = null;
 		if (dimensionsCond == 1) {
+			// This won't execute if dimensionsCond == 0
 			reshapedConditional = MatrixUtils.reshape(cond, cond.length, 1);
 		}
 		addObservations(MatrixUtils.reshape(var1, var1.length, 1),
@@ -372,6 +517,110 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 
 	}
 	
+	@Override
+	public void addObservations(double[][] var1, double[] var2,
+			double[] cond) throws Exception {
+		if ((dimensionsVar2 != 1) ||
+				((dimensionsCond != 1) && (dimensionsCond != 0))) {
+			throw new Exception("The number of dimensions for variables var2 and cond (having been initialised to " +
+					dimensionsVar2 + " & " +
+					dimensionsCond + ") can only be 1 (or 0 for conditional) when " +
+					"the addObservations(double[][],double[],double[]) and " + 
+					"setObservations(double[][],double[],double[]) methods are called");
+		}
+		double[][] reshapedConditional = null;
+		if (dimensionsCond == 1) {
+			// This won't execute if dimensionsCond == 0
+			reshapedConditional = MatrixUtils.reshape(cond, cond.length, 1);
+		}
+		addObservations(var1,
+						MatrixUtils.reshape(var2, var2.length, 1),
+						reshapedConditional);
+	}
+
+	@Override
+	public void addObservations(double[] var1, double[][] var2,
+			double[] cond) throws Exception {
+		if ((dimensionsVar1 != 1) ||
+				((dimensionsCond != 1) && (dimensionsCond != 0))) {
+			throw new Exception("The number of dimensions for variables var1 and cond (having been initialised to " +
+					dimensionsVar1 + " & " +
+					dimensionsCond + ") can only be 1 (or 0 for conditional) when " +
+					"the addObservations(double[],double[][],double[]) and " + 
+					"setObservations(double[],double[][],double[]) methods are called");
+		}
+		double[][] reshapedConditional = null;
+		if (dimensionsCond == 1) {
+			// This won't execute if dimensionsCond == 0
+			reshapedConditional = MatrixUtils.reshape(cond, cond.length, 1);
+		}
+		addObservations(MatrixUtils.reshape(var1, var1.length, 1),
+						var2,
+						reshapedConditional);
+	}
+
+	@Override
+	public void addObservations(double[] var1, double[] var2,
+			double[][] cond) throws Exception {
+		if ((dimensionsVar1 != 1) || (dimensionsVar2 != 1)) {
+			throw new Exception("The number of dimensions for variables var1 and var2 (having been initialised to " +
+					dimensionsVar1 + " & " +
+					dimensionsVar2 + ") can only be 1 when " +
+					"the addObservations(double[],double[],double[][]) and " + 
+					"setObservations(double[],double[],double[][]) methods are called");
+		}
+		addObservations(MatrixUtils.reshape(var1, var1.length, 1),
+						MatrixUtils.reshape(var2, var2.length, 1),
+						cond);
+	}
+
+	@Override
+	public void addObservations(double[] var1, double[][] var2,
+			double[][] cond) throws Exception {
+		if (dimensionsVar1 != 1) {
+			throw new Exception("The number of dimensions for variable var1 (having been initialised to " +
+					dimensionsVar1 + ") can only be 1 when " +
+					"the addObservations(double[],double[][],double[][]) and " + 
+					"setObservations(double[],double[][],double[][]) methods are called");
+		}
+		addObservations(MatrixUtils.reshape(var1, var1.length, 1),
+						var2,
+						cond);
+	}
+
+	@Override
+	public void addObservations(double[][] var1, double[] var2,
+			double[][] cond) throws Exception {
+		if (dimensionsVar2 != 1) {
+			throw new Exception("The number of dimensions for variable var2 (having been initialised to " +
+					dimensionsVar2 + ") can only be 1 when " +
+					"the addObservations(double[][],double[],double[][]) and " + 
+					"setObservations(double[][],double[],double[][]) methods are called");
+		}
+		addObservations(var1,
+						MatrixUtils.reshape(var2, var2.length, 1),
+						cond);
+	}
+
+	@Override
+	public void addObservations(double[][] var1, double[][] var2,
+			double[] cond) throws Exception {
+		if ((dimensionsCond != 1) && (dimensionsCond != 0)) {
+			throw new Exception("The number of dimensions for variable cond (having been initialised to " +
+					dimensionsCond + ") can only be 1 or 0 when " +
+					"the addObservations(double[][],double[][],double[]) and " + 
+					"setObservations(double[][],double[][],double[]) methods are called");
+		}
+		double[][] reshapedConditional = null;
+		if (dimensionsCond == 1) {
+			// This won't execute if dimensionsCond == 0
+			reshapedConditional = MatrixUtils.reshape(cond, cond.length, 1);
+		}
+		addObservations(var1,
+						var2,
+						reshapedConditional);
+	}
+
 	/**
 	 * A non-overloaded method signature for addObservations with 1D arguments, as there have been
 	 *  some problems calling overloaded versions of setObservations from jpype. 
@@ -389,6 +638,12 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 	public void addObservations(double[][] var1, double[][] var2,
 			double[][] cond,
 			int startTime, int numTimeSteps) throws Exception {
+		addObservations(var1, var2, cond, startTime, numTimeSteps, observationSetIndex++);
+	}
+	
+	protected void addObservations(double[][] var1, double[][] var2,
+			double[][] cond,
+			int startTime, int numTimeSteps, int observationSetIndexToUse) throws Exception {
 		if (vectorOfVar1Observations == null) {
 			// startAddObservations was not called first
 			throw new RuntimeException("User did not call startAddObservations before addObservations");
@@ -402,15 +657,14 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 			condToAdd = new double[numTimeSteps][];
 			System.arraycopy(cond, startTime, condToAdd, 0, numTimeSteps);
 		}
-		addObservations(var1ToAdd, var2ToAdd, condToAdd);
+		addObservationsTrackObservationIDs(var1ToAdd, var2ToAdd, condToAdd, observationSetIndexToUse, startTime);
 	}
 
-	@Override
-	public void setObservations(double[][] var1, double[][] var2,
+	public void addObservations(double[][] var1, double[][] var2,
 			double[][] cond,
 			boolean[] var1Valid, boolean[] var2Valid,
 			boolean[] condValid) throws Exception {
-		
+
 		Vector<int[]> startAndEndTimePairs =
 				computeStartAndEndTimePairs(var1Valid, var2Valid, condValid);
 		
@@ -419,23 +673,22 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		for (int[] timePair : startAndEndTimePairs) {
 			int startTime = timePair[0];
 			int endTime = timePair[1];
-			addObservations(var1, var2, cond, startTime, endTime - startTime + 1);
+			addObservations(var1, var2, cond, startTime, endTime - startTime + 1, observationSetIndex);
 		}
+		observationSetIndex++;
 		finaliseAddObservations();
 	}
-
-	@Override
-	public void setObservations(double[][] var1, double[][] var2,
+	
+	public void addObservations(double[][] var1, double[][] var2,
 			double[][] cond,
 			boolean[][] var1Valid, boolean[][] var2Valid,
 			boolean[][] condValid) throws Exception {
-
 		boolean[] allVar1Valid = MatrixUtils.andRows(var1Valid);
 		boolean[] allVar2Valid = MatrixUtils.andRows(var2Valid);
 		boolean[] allCondValid = MatrixUtils.andRows(condValid);
-		setObservations(var1, var2, cond, allVar1Valid, allVar2Valid, allCondValid);
+		addObservations(var1, var2, cond, allVar1Valid, allVar2Valid, allCondValid);
 	}
-
+	
 	/**
 	 * Signal that the observations are now all added, PDFs can now be constructed.
 	 * 
@@ -458,10 +711,14 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		var1Observations = new double[totalObservations][dimensionsVar1];
 		var2Observations = new double[totalObservations][dimensionsVar2];
 		condObservations = new double[totalObservations][dimensionsCond];
+		observationSetIndices = new int[totalObservations];
+		observationTimePoints = new int[totalObservations];
 		
 		int startObservation = 0;
 		Iterator<double[][]> iteratorVar2 = vectorOfVar2Observations.iterator();
 		Iterator<double[][]> iteratorCond = vectorOfCondObservations.iterator();
+		Iterator<Integer> iteratorObsSetIndices = vectorOfObservationSetIndices.iterator();
+		Iterator<Integer> iteratorObsStartTimePoints = vectorOfObservationStartTimePoints.iterator();
 		for (double[][] var1 : vectorOfVar1Observations) {
 			double[][] var2 = iteratorVar2.next();
 			double[][] cond = iteratorCond.next();
@@ -478,6 +735,12 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 					condObservations, startObservation, 0,
 					cond.length, dimensionsCond);
 			} // else we can do nothing there
+			// And update which observation set and time index each sample came from:
+			Arrays.fill(observationSetIndices, startObservation, startObservation+var1.length, iteratorObsSetIndices.next());
+			int firstTimeSampleId = iteratorObsStartTimePoints.next();
+			for (int i = 0; i < var1.length; i++) {
+				observationTimePoints[startObservation + i] = firstTimeSampleId + i;
+			}
 			startObservation += var2.length;
 		}
 		
@@ -502,6 +765,9 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		// Add Gaussian noise of std dev noiseLevel to the data if required
 		if (addNoise) {
 			Random random = new Random();
+	    	if (noiseSeedSet) {
+	    		random.setSeed(noiseSeed);
+	    	}
 			for (int r = 0; r < var1Observations.length; r++) {
 				for (int c = 0; c < dimensionsVar1; c++) {
 					var1Observations[r][c] +=
@@ -668,6 +934,34 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		// Compute the MI
 		return miSurrogateCalculator.computeAverageLocalOfObservations();
 	}
+	
+	@Override
+	public double[] computeLocalUsingPreviousObservations(double[] states1, double[] states2, double[][] condStates)
+			throws Exception {
+		if ((dimensionsVar1 != 1) || (dimensionsVar2 != 1)) {
+			throw new Exception("The number of source and dest dimensions (having been initialised to " +
+					dimensionsVar1 + " and " + dimensionsVar2 + ") can only be 1 when " +
+					"the univariate computeLocalUsingPreviousObservations(double[],double[],double[][]) " + 
+					"method is called");
+		}
+		return computeLocalUsingPreviousObservations(MatrixUtils.reshape(states1, states1.length, 1),
+						MatrixUtils.reshape(states2, states2.length, 1),
+						condStates);
+	}
+
+	@Override
+	public double[] computeLocalUsingPreviousObservations(double[] states1, double[] states2, double[] condStates)
+			throws Exception {
+		if ((dimensionsVar1 != 1) || (dimensionsVar2 != 1) || (dimensionsCond != 1)) {
+			throw new Exception("The number of source, dest and conditional dimensions (having been initialised to " +
+					dimensionsVar1 + ", " + dimensionsVar2 + " and " + dimensionsCond + ") can only be 1 when " +
+					"the univariate computeLocalUsingPreviousObservations(double[],double[],double[]) " + 
+					"method is called");
+		}
+		return computeLocalUsingPreviousObservations(MatrixUtils.reshape(states1, states1.length, 1),
+						MatrixUtils.reshape(states2, states2.length, 1),
+						MatrixUtils.reshape(condStates, condStates.length, 1));
+	}
 
 	@Override
 	public void setDebug(boolean debug) {
@@ -763,4 +1057,13 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
 		return addedMoreThanOneObservationSet;
 	}	
 
+	@Override
+	public int[] getObservationSetIndices() {
+		return observationSetIndices;
+	}
+
+	@Override
+	public int[] getObservationTimePoints() {
+		return observationTimePoints;
+	}
 }

java/source/infodynamics/measures/continuous/EntropyCalculator.java

@@ -48,8 +48,11 @@ Theory' (John Wiley & Sons, New York, 1991).</li>
  */
 public interface EntropyCalculator extends InfoMeasureCalculatorContinuous {
 
-	// TODO Add addObservations() methods for entropy calculator
-	
+	// TODO Add addObservations() methods for entropy calculator.
+	//  It's currently in EntropyCalculatorMultivariate; bring in
+	//  here when we make all univariate entropy calculators use their 
+	//  underlying multivariate form.
+
 	/**
 	 * Sets the samples from which to compute the PDF for the entropy.
 	 * Should only be called once, the last call contains the

java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariate.java

@@ -61,6 +61,26 @@ public interface EntropyCalculatorMultiVariate
 	 * Property name for the number of dimensions
 	 */
 	public static final String NUM_DIMENSIONS_PROP_NAME = "NUM_DIMENSIONS";
+	/**
+	 * Property for whether we normalise the incoming observations to mean 0,
+	 * standard deviation 1.
+	 */
+	public static final String NORMALISE_PROP_NAME = "NORMALISE";
+	/**
+	 * Property name for an amount of random Gaussian noise to be
+	 *  added to the data (default is 1e-8, matching the MILCA toolkit).
+	 */
+	public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD";
+	/**
+	 * Property name for the seed for the random number generator for noise to be
+	 *  added to the data (default is no seed)
+	 */
+	public static final String PROP_NOISE_SEED = "NOISE_SEED";
+	/**
+	 * Property value to indicate no seed for the random number generator for noise to be
+	 *  added to the data
+	 */
+	public static final String NOISE_NO_SEED_VALUE = "NONE";
 	
 	/**
 	 * Set properties for the underlying calculator implementation.
@@ -72,6 +92,20 @@ public interface EntropyCalculatorMultiVariate
 	 * <ul>
 	 *  <li>{@link #NUM_DIMENSIONS_PROP_NAME} -- number of dimensions in the joint
 	 * 		variable that we are computing the entropy of.</li>
+	 * 	<li>{@link #NORMALISE_PROP_NAME} -- whether to normalise the incoming variable values
+	 * 			to mean 0, standard deviation 1, or not (default false). Sets {@link #normalise}.</li>
+	 *  <li>{@link #PROP_ADD_NOISE} -- a standard deviation for an amount of
+	 *      random Gaussian noise to add to
+	 *      each variable, to avoid having neighbourhoods with artificially
+	 *      large counts. (We also accept "false" to indicate "0".)
+	 *      The amount is added in after any normalisation,
+	 *      so can be considered as a number of standard deviations of the data.
+	 *      Default is 0 for most estimators; this is strongly recommended by
+	 *      by Kraskov for the KSG method though, so for the Kozachenko estimator we 
+	 *      use 1e-8 to match the MILCA toolkit (though note it adds in
+	 *      a random amount of noise in [0,noiseLevel) ).</li>
+	 *  <li>{@link #PROP_NOISE_SEED} -- a long value seed for the random noise generator or
+	 *      the string {@link MutualInfoCalculatorMultiVariate#NOISE_NO_SEED_VALUE} for no seed (default)</li>
 	 * </ul>
 	 *  
 	 * <p>Unknown property values are ignored.</p>
@@ -94,6 +128,46 @@ public interface EntropyCalculatorMultiVariate
 	 */
 	public void initialise(int dimensions);
 	
+	/**
+	 * Signal that we will add in the samples for computing the PDF 
+	 * from several disjoint time-series or trials via calls to
+	 * "addObservations" rather than "setObservations" type methods
+	 * (defined by the child interfaces and classes).
+	 */
+	public void startAddObservations();
+	
+	/**
+	 * Add more observations for which to compute the PDFs for the entropy.
+	 * May be called multiple times between {@link #startAddObservations()} and
+	 * {@link #finaliseAddObservations()}.
+	 * 
+	 * @param observations multivariate time series of observations; first index
+	 *  is time step, second index is variable number (total should match dimensions
+	 *  supplied to {@link #initialise(int)}
+	 * @throws Exception if the dimensions of the observations do not match 
+	 *  the expected value supplied in {@link #initialise(int)}; implementations
+	 *  may throw other more specific exceptions also.
+	 */
+	public void addObservations(double[][] observations) throws Exception;
+
+	/**
+	 * Add more samples from which to compute the PDF for the entropy.
+	 * Only allowed to be called when set up for dimension == 1.
+	 * May be called multiple times between {@link #startAddObservations()} and
+	 * {@link #finaliseAddObservations()}.
+	 * 
+	 * @param observations array of (univariate) samples
+	 * @throws Exception if the expected dimensions were not 1.
+	 */
+	public void addObservations(double[] observations) throws Exception;
+
+	/**
+	 * Signal that the observations are now all added, PDFs can now be constructed.
+	 * 
+	 * @throws Exception 
+	 */
+	public void finaliseAddObservations() throws Exception;
+	
 	/**
 	 * Set the observations for which to compute the PDFs for the entropy
 	 * Should only be called once, the last call contains the

java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariateCommon.java

@@ -0,0 +1,342 @@
+/*
+ *  Java Information Dynamics Toolkit (JIDT)
+ *  Copyright (C) 2024, 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/>.
+ */
+package infodynamics.measures.continuous;
+
+import java.util.Random;
+import java.util.Vector;
+
+import infodynamics.utils.MatrixUtils;
+
+/**
+ * Implements {@link EntropyCalculatorMultiVariate} to provide a base
+ * class with common functionality for child class implementations of
+ * {@link EntropyCalculatorMultiVariate}
+ * via various estimators. 
+ * 
+ * <p>These various estimators include: e.g. box-kernel estimation, Kozachenko, etc
+ * (see the child classes linked above).
+ * </p>
+ * 
+ * <p>Usage is as outlined in {@link EntropyCalculatorMultiVariate}.</p>
+ * 
+  * @author Joseph Lizier (<a href="joseph.lizier at gmail.com">email</a>,
+ * <a href="http://lizier.me/joseph/">www</a>)
+ */
+public abstract class EntropyCalculatorMultiVariateCommon implements EntropyCalculatorMultiVariate {
+	/**
+	 * Whether we're in debug mode
+	 */
+	protected boolean debug = false;
+	/**
+	 * Number of observations supplied
+	 */
+	protected int totalObservations = 0;
+	/**
+	 * Number of joint variables/dimensions
+	 */
+	protected int dimensions = 1;
+	/**
+	 * Track whether we've computed the average for the supplied
+	 *  observations yet
+	 */
+	protected boolean isComputed = false;
+	/**
+	 * Last computed average entropy
+	 */
+	protected double lastAverage = 0.0;
+	/**
+	 * Whether we normalise the incoming observations to mean 0,
+	 * standard deviation 1.
+	 */
+	protected boolean normalise = true;
+	/**
+	 * Whether to add an amount of random noise to the incoming data
+	 */
+	protected boolean addNoise = false;
+	/**
+	 * Amount of random Gaussian noise to add to the incoming data
+	 */
+	protected double noiseLevel = (double) 0.0;
+	/**
+	 * Has the user set a seed for the random noise
+	 */
+	protected boolean noiseSeedSet = false;
+	/**
+	 * Seed that the user set for the random noise
+	 */
+	protected long noiseSeed = 0;
+	/**
+	 * Storage for observations supplied via {@link #addObservations(double[][])}
+	 * type calls
+	 */
+	protected Vector<double[][]> vectorOfObservations;
+	/**
+	 * The set of observations, retained in case the user wants to retrieve the local
+	 *  entropy values of these.
+	 */
+	protected double[][] observations;
+
+	/**
+	 * Default constructor
+	 */
+	public EntropyCalculatorMultiVariateCommon() {
+		// Nothing to do
+	}
+	
+	@Override
+	public void initialise() throws Exception {
+		initialise(dimensions);
+	}
+
+	@Override
+	public void initialise(int dimensions) {
+		this.dimensions = dimensions;
+		observations = null;
+		totalObservations = 0;
+		isComputed = false;
+		lastAverage = 0;
+	}
+
+	@Override
+	public void setProperty(String propertyName, String propertyValue) throws Exception {
+
+		boolean propertySet = true;
+
+		if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
+			dimensions = Integer.parseInt(propertyValue);
+		} else if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
+			normalise = Boolean.parseBoolean(propertyValue);
+		} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
+			if (propertyValue.equals("0") ||
+					propertyValue.equalsIgnoreCase("false")) {
+				addNoise = false;
+				noiseLevel = 0;
+			} else {
+				addNoise = true;
+				noiseLevel = Double.parseDouble(propertyValue);
+			}
+	    } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) {
+	    	if (propertyValue.equals(NOISE_NO_SEED_VALUE)) {
+	    		noiseSeedSet = false;
+	    	} else {
+	    		noiseSeedSet = true;
+	    		noiseSeed = Long.parseLong(propertyValue);
+	    	}
+		} else {
+			// No property was set
+			propertySet = false;
+		}
+		if (debug && propertySet) {
+			System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
+					" to " + propertyValue);
+		}
+	}
+	
+	@Override
+	public String getProperty(String propertyName) throws Exception {
+		if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
+			return Integer.toString(dimensions);
+		} else if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
+			return Boolean.toString(normalise);
+	    } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
+	        return Double.toString(noiseLevel);
+	    } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) {
+	    	if (noiseSeedSet) {
+	    		return Long.toString(noiseSeed);
+	    	} else {
+	    		return NOISE_NO_SEED_VALUE;
+	    	}
+		} else {
+			// No property was set, and no superclass to call:
+			return null;
+		}
+	}
+
+	@Override
+	public void startAddObservations() {
+		isComputed = false;
+		totalObservations = 0;
+		observations = null;
+		vectorOfObservations = new Vector<double[][]>();
+	}
+
+	@Override
+	public void finaliseAddObservations() throws Exception {
+
+		observations = new double[totalObservations][dimensions];
+		
+		// Construct the joint vectors from the given observations
+		int startObservation = 0;
+		for (double[][] obs : vectorOfObservations) {
+			if ((obs == null) || (obs.length < 1)) {
+				continue;
+			}
+			// Dimension was checked earlier by addObservations.
+			// Copy the data from these given observations into our master array
+			MatrixUtils.arrayCopy(obs, 0, 0,
+					observations, startObservation, 0,
+					obs.length, dimensions);
+			startObservation += obs.length;
+		}
+
+		// We don't need to keep the vector of observation sets anymore:
+		vectorOfObservations = null;
+
+		if (addNoise) {
+			Random random = new Random();
+			if (noiseSeedSet) {
+				random.setSeed(noiseSeed);
+			}
+			// Add Gaussian noise of std dev noiseLevel to the data
+			for (int r = 0; r < totalObservations; r++) {
+				for (int c = 0; c < dimensions; c++) {
+					observations[r][c] += random.nextGaussian()*noiseLevel;
+				}
+			}
+		}
+	}
+
+	@Override
+	public void setObservations(double[][] observations) throws Exception {
+		startAddObservations();
+		addObservations(observations);
+		finaliseAddObservations();
+	}
+
+	/*
+	 * (non-Javadoc)
+	 * @see infodynamics.measures.continuous.EntropyCalculator#setObservations(double[])
+	 * 
+	 * This method here to ensure we make compatibility with the 
+	 *  EntropyCalculator interface; only valid if dimension > 1
+	 */
+	@Override
+	public void setObservations(double[] observations) throws Exception {
+		startAddObservations();
+		addObservations(observations);
+		finaliseAddObservations();
+	}
+
+	/**
+	 * Shortcut method to join two sets of samples into a joint multivariate 
+	 * 
+	 * Each row of the data is an observation; each column of
+	 *  the row is a new variable in the multivariate observation.
+	 * This method signature allows the user to call setObservations for
+	 *  joint time series without combining them into a single joint time
+	 *  series (we do the combining for them).
+	 * 
+	 * @param data1 observations1 few variables in the joint data
+	 * @param observations2 the other variables in the joint data
+	 * @throws Exception When the length of the two arrays of observations do not match.
+	 * @see #setObservations(double[][])
+	 */
+	@Deprecated
+	public void setObservations(double[][] observations1, double[][] observations2)
+			throws Exception {
+		startAddObservations();
+		addObservations(observations1, observations2);
+		finaliseAddObservations();
+	}
+
+	@Override
+	public void addObservations(double[][] observations) throws Exception {
+		if (vectorOfObservations == null) {
+		// startAddObservations was not called first
+			throw new RuntimeException("User did not call startAddObservations before addObservations");
+		}
+		if ((observations != null) && (observations.length > 0)) {
+			// Check the dimension:
+			if (observations[0].length != dimensions) {
+				throw new Exception(String.format("Cannot supply observations with dimension %d when expected dimension = %d",
+						observations[0].length, dimensions));
+			}
+			totalObservations += observations.length;
+		}
+		// Add the observations whether empty or else valid:
+		vectorOfObservations.add(observations);
+	}
+
+	/*
+	 * (non-Javadoc)
+	 * @see infodynamics.measures.continuous.EntropyCalculator#addObservations(double[])
+	 * 
+	 * This method here to ensure we make compatibility with the 
+	 *  EntropyCalculator interface; only valid if dimension > 1
+	 */
+	@Override
+	public void addObservations(double[] observations) throws Exception {
+		if (dimensions != 1) {
+			throw new Exception(String.format("Cannot set univariate observations when expected dimension = %d", dimensions));
+		}
+		double[][] observations2D = MatrixUtils.reshape(observations, observations.length, 1);
+		addObservations(observations2D);
+	}
+	
+	/**
+	 * Shortcut method to join two sets of samples into a joint multivariate. 
+	 * 
+	 * Each row of the data is an observation; each column of
+	 *  the row is a new variable in the multivariate observation.
+	 * This method signature allows the user to call addObservations for
+	 *  joint time series without combining them into a single joint time
+	 *  series (we do the combining for them).
+	 * 
+	 * @param data1 first few variables in the joint data
+	 * @param data2 the other variables in the joint data
+	 * @throws Exception When the length of the two arrays of observations do not match.
+	 * @see #addObservations(double[][])
+	 */
+	@Deprecated
+	public void addObservations(double[][] data1,
+			double[][] data2) throws Exception {
+		int timeSteps = data1.length;
+		if ((data1 == null) || (data2 == null)) {
+			throw new Exception("Cannot have null data arguments");
+		}
+		if (data1.length != data2.length) {
+			throw new Exception("Length of data1 (" + data1.length + ") is not equal to the length of data2 (" +
+					data2.length + ")");
+		}
+		int data1Variables = data1[0].length;
+		int data2Variables = data2[0].length;
+		double[][] data = new double[timeSteps][data1Variables + data2Variables];
+		for (int t = 0; t < timeSteps; t++) {
+			System.arraycopy(data1[t], 0, data[t], 0, data1Variables);
+			System.arraycopy(data2[t], 0, data[t], data1Variables, data2Variables);
+		}
+		// Now defer to the normal setObservations method
+		addObservations(data);
+	}
+
+	@Override
+	public int getNumObservations() throws Exception {
+		return totalObservations;
+	}
+
+	@Override
+	public double getLastAverage() {
+		return lastAverage;
+	}
+	
+	@Override
+	public void setDebug(boolean debug) {
+		this.debug = debug;
+	}
+
+}

java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariate.java

@@ -91,6 +91,35 @@ public interface MutualInfoCalculatorMultiVariate
 	 *  if the data is to be normalised, that will be done before adding this noise.
 	 */
 	public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD";
+	/**
+	 * Property name for the seed for the random number generator for noise to be
+	 *  added to the data (default is no seed)
+	 */
+	public static final String PROP_NOISE_SEED = "NOISE_SEED";
+	/**
+	 * Property value to indicate no seed for the random number generator for noise to be
+	 *  added to the data
+	 */
+	public static final String NOISE_NO_SEED_VALUE = "NONE";
+
+	/**
+	 * <p>As per {@link #addObservations(double[][], double[][])};
+     *  but also includes parameters to track which observation set
+     *  the samples came from. Intended to only be used by other 
+     *  estimator classes here and not by users directly.</p>
+	 * 
+	 * @param source multivariate observations for variable 1
+	 *  (first index is time or observation index, second is variable number)
+	 * @param destination multivariate observations for variable 2
+	 *  (first index is time or observation index, second is variable number)
+	 *  Length must match <code>source</code>, and their indices must correspond.
+	 * @param observationSetIndexToUse which set of observations these came fmor
+	 * @param startTimeIndex which was the first time index of these
+	 *  samples within that observation set.
+	 * @throws Exception
+	 */
+	public void addObservationsTrackObservationIDs(double[][] source, double[][] destination,
+			int observationSetIndexToUse, int startTimeIndex) throws Exception;
 
 	/**
 	 * Compute the mutual information if the observations of the
@@ -130,4 +159,18 @@ public interface MutualInfoCalculatorMultiVariate
 	 * @throws Exception for invalid property values
 	 */
 	public String getProperty(String propertyName) throws Exception;
+
+	/**
+	 * Retrieve an array indicating which observation set each sample came from
+	 * 
+	 * @return array of integers
+	 */
+	public int[] getObservationSetIndices();
+
+	/**
+	 * Retrieve an array indicating which time index within its observation set that sample came from
+	 * 
+	 * @return array of integers
+	 */
+	public int[] getObservationTimePoints();
 }

java/source/infodynamics/measures/continuous/MutualInfoMultiVariateCommon.java

@@ -23,6 +23,7 @@ import infodynamics.utils.EmpiricalMeasurementDistribution;
 import infodynamics.utils.MatrixUtils;
 import infodynamics.utils.RandomGenerator;
 
+import java.util.Arrays;
 import java.util.Iterator;
 import java.util.Random;
 import java.util.Vector;
@@ -86,6 +87,16 @@ public abstract class MutualInfoMultiVariateCommon implements
 	 */
 	protected double[] destStdsBeforeNorm;
 
+	/**
+	 * Track which observation set each sample came from
+	 */
+	protected int[] observationSetIndices;
+
+	/**
+	 * Track which sample index within an observation set that each sample came from
+	 */
+	protected int[] observationTimePoints;
+
 	/**
 	 * Total number of observations supplied.
 	 * Only valid after {@link #finaliseAddObservations()} is called.
@@ -127,6 +138,22 @@ public abstract class MutualInfoMultiVariateCommon implements
 	 * type calls
 	 */
 	protected Vector<double[][]> vectorOfDestinationObservations;
+	/**
+	 * Tracks separate (time-series) observation sets
+	 *  we are taking samples from
+	 */
+	protected int observationSetIndex = 0;
+	/**
+	 * Storage for which observation set each
+	 *  block of samples comes from
+	 */
+	protected Vector<Integer> vectorOfObservationSetIndices;
+	/**
+	 * Storage for start time point for the observation
+	 *  set within its block of samples
+	 */
+	protected Vector<Integer> vectorOfObservationStartTimePoints;
+
 	/**
 	 * Whether the user has supplied more than one (disjoint) set of samples
 	 */
@@ -140,9 +167,19 @@ public abstract class MutualInfoMultiVariateCommon implements
 	 */
 	protected boolean addNoise = false;
 	/**
-	 * Amount of random Gaussian noise to add to the incoming data
+	 * Amount of random Gaussian noise to add to the incoming data.
+	 * 0 by default except for KSG estimators (where it is recommended
+	 *  and 1e-8 is used to match MILCA toolkit)
 	 */
 	protected double noiseLevel = (double) 0.0;
+	/**
+	 * Has the user set a seed for the random noise
+	 */
+	protected boolean noiseSeedSet = false;
+	/**
+	 * Seed that the user set for the random noise
+	 */
+	protected long noiseSeed = 0;
 
 	/* (non-Javadoc)
 	 * @see infodynamics.measures.continuous.ChannelCalculatorCommon#initialise()
@@ -163,6 +200,11 @@ public abstract class MutualInfoMultiVariateCommon implements
 		sourceStdsBeforeNorm = null;
 		destMeansBeforeNorm = null;
 		destStdsBeforeNorm = null;
+		observationSetIndices = null;
+		observationTimePoints = null;
+		observationSetIndex = 0;
+		vectorOfObservationSetIndices = null;
+		vectorOfObservationStartTimePoints = null;
 		addedMoreThanOneObservationSet = false;
 	}
 
@@ -190,6 +232,8 @@ public abstract class MutualInfoMultiVariateCommon implements
 	 *      by Kraskov for the KSG method though, so for that estimator we 
 	 *      use 1e-8 to match the MILCA toolkit (though note it adds in
 	 *      a random amount of noise in [0,noiseLevel) ).</li>
+	 *  <li>{@link #PROP_NOISE_SEED} -- a long value seed for the random noise generator or
+	 *      the string {@link MutualInfoCalculatorMultiVariate#NOISE_NO_SEED_VALUE} for no seed (default)</li>
 	 * </ul>
 	 * 
 	 * <p>Unknown property values are ignored.</p>
@@ -218,6 +262,13 @@ public abstract class MutualInfoMultiVariateCommon implements
 	          addNoise = true;
 	          noiseLevel = Double.parseDouble(propertyValue);
 	        }
+	    } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) {
+	    	if (propertyValue.equals(NOISE_NO_SEED_VALUE)) {
+	    		noiseSeedSet = false;
+	    	} else {
+	    		noiseSeedSet = true;
+	    		noiseSeed = Long.parseLong(propertyValue);
+	    	}
 		} else {
 			// No property was set here
 			propertySet = false;
@@ -238,6 +289,12 @@ public abstract class MutualInfoMultiVariateCommon implements
 	        return Boolean.toString(normalise);
 	    } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
 	        return Double.toString(noiseLevel);
+	    } else if (propertyName.equalsIgnoreCase(PROP_NOISE_SEED)) {
+	    	if (noiseSeedSet) {
+	    		return Long.toString(noiseSeed);
+	    	} else {
+	    		return NOISE_NO_SEED_VALUE;
+	    	}
 		} else {
 			// No property was recognised here
 			return null;
@@ -296,14 +353,54 @@ public abstract class MutualInfoMultiVariateCommon implements
 		finaliseAddObservations();
 	}
 
+	@Override
+	public void setObservations(double[][] source, double[][] destination,
+			boolean[] sourceValid, boolean[] destValid) throws Exception {
+		startAddObservations();
+		addObservations(source, destination, sourceValid, destValid);
+		finaliseAddObservations();
+	}
+	
+	@Override
+	public void setObservations(double[] source, double[] destination,
+			boolean[] sourceValid, boolean[] destValid) throws Exception {
+
+		if ((dimensionsDest != 1) || (dimensionsSource != 1)) {
+			throw new Exception("The number of source and dest dimensions (having been initialised to " +
+					dimensionsSource + " and " + dimensionsDest + ") can only be 1 when " +
+					"the univariate addObservations(double[],double[]) and " + 
+					"setObservations(double[],double[]) methods are called");
+		}
+		setObservations(MatrixUtils.reshape(source, source.length, 1),
+				MatrixUtils.reshape(destination, destination.length, 1),
+				sourceValid, destValid);
+	}
+	
+	@Override
+	public void setObservations(double[][] source, double[][] destination,
+			boolean[][] sourceValid, boolean[][] destValid) throws Exception {
+
+		boolean[] allSourceValid = MatrixUtils.andRows(sourceValid);
+		boolean[] allDestValid = MatrixUtils.andRows(destValid);
+		setObservations(source, destination, allSourceValid, allDestValid);
+	}
+
 	@Override
 	public void startAddObservations() {
 		vectorOfSourceObservations = new Vector<double[][]>();
 		vectorOfDestinationObservations = new Vector<double[][]>();
+		vectorOfObservationSetIndices = new Vector<Integer>();
+		vectorOfObservationStartTimePoints = new Vector<Integer>();
 	}
 	
 	@Override
 	public void addObservations(double[][] source, double[][] destination) throws Exception {
+		// Use the current observationSetIndex and increment for next use:
+		addObservationsTrackObservationIDs(source, destination, observationSetIndex++, 0);
+	}
+	
+	public void addObservationsTrackObservationIDs(double[][] source, double[][] destination,
+			int observationSetIndexToUse, int startTimeIndex) throws Exception {
 		if (vectorOfSourceObservations == null) {
 			// startAddObservations was not called first
 			throw new RuntimeException("User did not call startAddObservations before addObservations");
@@ -326,6 +423,8 @@ public abstract class MutualInfoMultiVariateCommon implements
 		}
 		vectorOfSourceObservations.add(source);
 		vectorOfDestinationObservations.add(destination);
+		vectorOfObservationSetIndices.add(observationSetIndexToUse);
+		vectorOfObservationStartTimePoints.add(startTimeIndex);
 	}
 
 	/**
@@ -391,8 +490,14 @@ public abstract class MutualInfoMultiVariateCommon implements
 				MatrixUtils.reshape(destination, destination.length, 1));
 	}
 	
+	@Override
 	public void addObservations(double[][] source, double[][] destination,
 			int startTime, int numTimeSteps) throws Exception {
+		addObservations(source, destination, startTime, numTimeSteps, observationSetIndex++);
+	}
+	
+	protected void addObservations(double[][] source, double[][] destination,
+			int startTime, int numTimeSteps, int observationSetIndexToUse) throws Exception {
 		if (vectorOfSourceObservations == null) {
 			// startAddObservations was not called first
 			throw new RuntimeException("User did not call startAddObservations before addObservations");
@@ -403,12 +508,12 @@ public abstract class MutualInfoMultiVariateCommon implements
 		}
 		double[][] sourceToAdd = new double[numTimeSteps][];
 		System.arraycopy(source, startTime, sourceToAdd, 0, numTimeSteps);
-		vectorOfSourceObservations.add(sourceToAdd);
 		double[][] destToAdd = new double[numTimeSteps][];
 		System.arraycopy(destination, startTime, destToAdd, 0, numTimeSteps);
-		vectorOfDestinationObservations.add(destToAdd);
+		addObservationsTrackObservationIDs(sourceToAdd, destToAdd, observationSetIndexToUse, startTime);
 	}
 
+	@Override
 	public void addObservations(double[] source, double[] destination,
 			int startTime, int numTimeSteps) throws Exception {
 		
@@ -423,8 +528,14 @@ public abstract class MutualInfoMultiVariateCommon implements
 						startTime, numTimeSteps);
 	}
 
-
-	public void setObservations(double[][] source, double[][] destination,
+	public void addObservations(double[] source, double[] destination,
+			boolean[] sourceValid, boolean[] destValid) throws Exception {
+		addObservations(MatrixUtils.reshape(source, source.length, 1),
+				MatrixUtils.reshape(destination, destination.length, 1),
+				sourceValid, destValid);
+	}
+	
+	public void addObservations(double[][] source, double[][] destination,
 			boolean[] sourceValid, boolean[] destValid) throws Exception {
 		
 		Vector<int[]> startAndEndTimePairs = computeStartAndEndTimePairs(sourceValid, destValid);
@@ -434,32 +545,18 @@ public abstract class MutualInfoMultiVariateCommon implements
 		for (int[] timePair : startAndEndTimePairs) {
 			int startTime = timePair[0];
 			int endTime = timePair[1];
-			addObservations(source, destination, startTime, endTime - startTime + 1);
+			addObservations(source, destination, startTime, endTime - startTime + 1, observationSetIndex);
 		}
+		observationSetIndex++;
 		finaliseAddObservations();
 	}
 
-	@Override
-	public void setObservations(double[] source, double[] destination,
-			boolean[] sourceValid, boolean[] destValid) throws Exception {
-
-		if ((dimensionsDest != 1) || (dimensionsSource != 1)) {
-			throw new Exception("The number of source and dest dimensions (having been initialised to " +
-					dimensionsSource + " and " + dimensionsDest + ") can only be 1 when " +
-					"the univariate addObservations(double[],double[]) and " + 
-					"setObservations(double[],double[]) methods are called");
-		}
-		setObservations(MatrixUtils.reshape(source, source.length, 1),
-				MatrixUtils.reshape(destination, destination.length, 1),
-				sourceValid, destValid);
-	}
-	
-	public void setObservations(double[][] source, double[][] destination,
+	public void addObservations(double[][] source, double[][] destination,
 			boolean[][] sourceValid, boolean[][] destValid) throws Exception {
-
+		
 		boolean[] allSourceValid = MatrixUtils.andRows(sourceValid);
 		boolean[] allDestValid = MatrixUtils.andRows(destValid);
-		setObservations(source, destination, allSourceValid, allDestValid);
+		addObservations(source, destination, allSourceValid, allDestValid);
 	}
 
 	/**
@@ -482,11 +579,15 @@ public abstract class MutualInfoMultiVariateCommon implements
 		}
 		destObservations = new double[totalObservations][dimensionsDest];
 		sourceObservations = new double[totalObservations][dimensionsSource];
+		observationSetIndices = new int[totalObservations];
+		observationTimePoints = new int[totalObservations];
 		
 		// Construct the joint vectors from the given observations
 		//  (removing redundant data which is outside any timeDiff)
 		int startObservation = 0;
 		Iterator<double[][]> iterator = vectorOfDestinationObservations.iterator();
+		Iterator<Integer> iteratorObsSetIndices = vectorOfObservationSetIndices.iterator();
+		Iterator<Integer> iteratorObsStartTimePoints = vectorOfObservationStartTimePoints.iterator();
 		for (double[][] source : vectorOfSourceObservations) {
 			double[][] destination = iterator.next();
 			// Copy the data from these given observations into our master 
@@ -497,7 +598,14 @@ public abstract class MutualInfoMultiVariateCommon implements
 			MatrixUtils.arrayCopy(destination, timeDiff, 0,
 					destObservations, startObservation, 0,
 					destination.length - timeDiff, dimensionsDest);
-			startObservation += destination.length - timeDiff;
+			int numNewObservations = destination.length - timeDiff;
+			// And update which observation set and time index each sample came from:
+			Arrays.fill(observationSetIndices, startObservation, startObservation + numNewObservations, iteratorObsSetIndices.next());
+			int firstTimeSampleId = iteratorObsStartTimePoints.next(); // This is the first sample of the source; storing the time index for destination below.
+			for (int i = 0; i < numNewObservations; i++) {
+				observationTimePoints[startObservation + i] = firstTimeSampleId + timeDiff + i;
+			}
+			startObservation += numNewObservations;
 		}
 		if (vectorOfSourceObservations.size() > 1) {
 			addedMoreThanOneObservationSet = true;
@@ -518,6 +626,9 @@ public abstract class MutualInfoMultiVariateCommon implements
 		// Add Gaussian noise of std dev noiseLevel to the data if required
 		if (addNoise) {
 			Random random = new Random();
+	    	if (noiseSeedSet) {
+	    		random.setSeed(noiseSeed);
+	    	}
 			for (int r = 0; r < sourceObservations.length; r++) {
 				for (int c = 0; c < dimensionsSource; c++) {
 					sourceObservations[r][c] +=
@@ -712,6 +823,7 @@ public abstract class MutualInfoMultiVariateCommon implements
 		return miSurrogateCalculator.computeAverageLocalOfObservations();
 	}
 
+	@Override
 	public double[] computeLocalUsingPreviousObservations(
 			double[] newSourceObservations, double[] newDestObservations)
 					throws Exception {
@@ -834,4 +946,14 @@ public abstract class MutualInfoMultiVariateCommon implements
 		}
 		return startAndEndTimePairs;
 	}	
+
+	@Override
+	public int[] getObservationSetIndices() {
+		return observationSetIndices;
+	}
+
+	@Override
+	public int[] getObservationTimePoints() {
+		return observationTimePoints;
+	}
 }

java/source/infodynamics/measures/continuous/TransferEntropyCalculatorViaCondMutualInfo.java

@@ -148,6 +148,12 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 	 */
 	protected Vector<boolean[]> vectorOfValidityOfDestination;
 
+	/**
+	 * Store the time index at which we were asked to start taking time-series
+	 *  observations for each  observation set.
+	 */
+	protected Vector<Integer> vectorOfOffsetsInTimeSeries;
+	
 	/**
 	 * Array of the number of observations added by each separate call to
 	 * {@link #addObservations(double[], double[])}, in order of which those calls
@@ -352,6 +358,7 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 		vectorOfDestinationTimeSeries = null;
 		vectorOfValidityOfSource = null;
 		vectorOfValidityOfDestination = null;
+		vectorOfOffsetsInTimeSeries = null;
 		separateNumObservations = new int[] {};
 	}
 
@@ -501,22 +508,38 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 		finaliseAddObservations();
 	}
 
+	@Override
+	public void setObservations(double[] source, double[] destination,
+			boolean[] sourceValid, boolean[] destValid) throws Exception {
+		
+		startAddObservations();
+		addObservations(source, destination, sourceValid, destValid);
+		finaliseAddObservations();
+	}
+
 	@Override
 	public void startAddObservations() {
 		vectorOfSourceTimeSeries = new Vector<double[]>();
 		vectorOfDestinationTimeSeries = new Vector<double[]>();
 		vectorOfValidityOfSource = new Vector<boolean[]>();
 		vectorOfValidityOfDestination = new Vector<boolean[]>();
+		vectorOfOffsetsInTimeSeries = new Vector<Integer>();
 	}
 	
 	@Override
 	public void addObservations(double[] source, double[] destination)
 			throws Exception {
+		addObservationsParsed(source, destination, 0); // No offset here
+	}
+
+	protected void addObservationsParsed(double[] source, double[] destination, int startTimeStep)
+			throws Exception {
 		// Store these observations in our vectors for now
 		vectorOfSourceTimeSeries.add(source);
 		vectorOfDestinationTimeSeries.add(destination);
 		vectorOfValidityOfSource.add(null); // All observations were valid
 		vectorOfValidityOfDestination.add(null); // All observations were valid
+		vectorOfOffsetsInTimeSeries.add(startTimeStep);
 	}
 
 	/**
@@ -536,6 +559,7 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 		vectorOfDestinationTimeSeries.add(destination);
 		vectorOfValidityOfSource.add(sourceValid); // All observations were valid
 		vectorOfValidityOfDestination.add(destValid); // All observations were valid
+		vectorOfOffsetsInTimeSeries.add(0); // no offset here
 	}
 
 	/**
@@ -554,6 +578,8 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 	 * @param delay_in_use source-target delay to use
 	 * @param source time series of source observations
 	 * @param destination time series of destination observations
+	 * @param observationSetIndex which observation set these samples came from
+	 * @param offsetInOriginalTimeSeries offset of the samples in their original time series
 	 * @return the number of observations added
 	 * @throws Exception
 	 */
@@ -561,7 +587,8 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 			ConditionalMutualInfoCalculatorMultiVariate condMiCalc_in_use,
 			int k_in_use, int k_tau_in_use, int l_in_use, int l_tau_in_use,
 			int delay_in_use,
-			double[] source, double[] destination) throws Exception {
+			double[] source, double[] destination,
+			int observationSetIndex, int offsetInOriginalTimeSeries) throws Exception {
 		if (source.length != destination.length) {
 			throw new Exception(String.format("Source and destination lengths (%d and %d) must match!",
 					source.length, destination.length));
@@ -587,7 +614,8 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 				MatrixUtils.makeDelayEmbeddingVector(source, l_in_use, l_tau_in_use,
 						startTimeForFirstDestEmbedding_in_use + 1 - delay_in_use,
 						source.length - startTimeForFirstDestEmbedding_in_use - 1);
-		condMiCalc_in_use.addObservations(currentSourcePastVectors, currentDestNextVectors, currentDestPastVectors);
+		condMiCalc_in_use.addObservationsTrackObservationIDs(currentSourcePastVectors, currentDestNextVectors, currentDestPastVectors,
+				observationSetIndex, offsetInOriginalTimeSeries + startTimeForFirstDestEmbedding_in_use + 1);
 		return destination.length - startTimeForFirstDestEmbedding_in_use - 1;
 	}
 
@@ -615,6 +643,7 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 	 * 	the source at that index is valid.
 	 * @param destValid array (with indices the same as destination) indicating whether
 	 * 	the destination at that index is valid.
+	 * @param observationSetIndex which observation set these samples came from
 	 * @return total number of observations added 
 	 * @throws Exception
 	 */
@@ -623,7 +652,8 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 			int k_in_use, int k_tau_in_use, int l_in_use, int l_tau_in_use,
 			int delay_in_use,
 			double[] source, double[] destination,
-			boolean[] sourceValid, boolean[] destValid) throws Exception {
+			boolean[] sourceValid, boolean[] destValid,
+			int observationSetIndex) throws Exception {
 		
 		// Compute the start and end time pairs using our embedding parameters:
 		Vector<int[]> startAndEndTimePairs =
@@ -638,7 +668,8 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 					condMiCalc_in_use, k_in_use, k_tau_in_use, l_in_use,
 					l_tau_in_use, delay_in_use,
 					MatrixUtils.select(source, startTime, endTime - startTime + 1),
-					MatrixUtils.select(destination, startTime, endTime - startTime + 1));
+					MatrixUtils.select(destination, startTime, endTime - startTime + 1),
+					observationSetIndex, startTime);
 		}
 		return totalObservationsAdded;
 	}
@@ -654,8 +685,8 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 			// There are not enough observations given the arguments here
 			throw new Exception("Not enough observations to set here given startTime and numTimeSteps parameters");
 		}
-		addObservations(MatrixUtils.select(source, startTime, numTimeSteps),
-					    MatrixUtils.select(destination, startTime, numTimeSteps));
+		addObservationsParsed(MatrixUtils.select(source, startTime, numTimeSteps),
+					    MatrixUtils.select(destination, startTime, numTimeSteps), startTime);
 	}
 
 	/**
@@ -841,6 +872,7 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 		vectorOfDestinationTimeSeries = null; // No longer required
 		vectorOfValidityOfSource = null;
 		vectorOfValidityOfDestination = null;
+		vectorOfOffsetsInTimeSeries = null;
 	}
 
 	/**
@@ -900,6 +932,7 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 		Iterator<double[]> destIterator = vectorOfDestinationTimeSeries.iterator();
 		Iterator<boolean[]> sourceValidityIterator = vectorOfValidityOfSource.iterator();
 		Iterator<boolean[]> destValidityIterator = vectorOfValidityOfDestination.iterator();
+		Iterator<Integer> offsetsInTimeSeriesIterator = vectorOfOffsetsInTimeSeries.iterator();
 		int[] separateNumObservationsArray = new int[vectorOfDestinationTimeSeries.size()];
 		int setNum = 0;
 		for (double[] source : vectorOfSourceTimeSeries) {
@@ -911,12 +944,12 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 				// Add the whole time-series
 				observationsAddedThisTime = addObservationsWithGivenParams(
 						condMiCalc_in_use, k_in_use, k_tau_in_use, l_in_use,
-						l_tau_in_use, delay_in_use, source, destination);
+						l_tau_in_use, delay_in_use, source, destination, setNum, offsetsInTimeSeriesIterator.next());
 			} else {
 				observationsAddedThisTime = addObservationsWithGivenParams(
 						condMiCalc_in_use, k_in_use, k_tau_in_use, l_in_use,
 						l_tau_in_use, delay_in_use, source, destination,
-						sourceValidity, destValidity);
+						sourceValidity, destValidity, setNum);
 			}
 			separateNumObservationsArray[setNum++] = observationsAddedThisTime;
 		}
@@ -927,20 +960,6 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 		return separateNumObservationsArray;
 	}
 	
-	@Override
-	public void setObservations(double[] source, double[] destination,
-			boolean[] sourceValid, boolean[] destValid) throws Exception {
-		
-		startAddObservations();
-		// Add these observations and the indication of their validity
-		//  for later analysis:
-		vectorOfSourceTimeSeries.add(source);
-		vectorOfDestinationTimeSeries.add(destination);
-		vectorOfValidityOfSource.add(sourceValid);
-		vectorOfValidityOfDestination.add(destValid);
-		finaliseAddObservations();
-	}
-
 	/**
 	 * Compute a vector of start and end pairs of time points, between which we have
 	 *  valid series of both source and destinations. (I.e. all points within the
@@ -1131,6 +1150,22 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
 		return separateNumObservations;
 	}
 	
+	/**
+	 * Retrieve an array indicating which observation set each sample came from
+	 * @return
+	 */
+	public int[] getObservationSetIndices() {
+		return condMiCalc.getObservationSetIndices();
+	}
+
+	/**
+	 * Retrieve an array indicating which time index within its observation set that sample came from
+	 * @return
+	 */
+	public int[] getObservationTimePoints() {
+		return condMiCalc.getObservationTimePoints();
+	}
+
 	@Override
 	public boolean getAddedMoreThanOneObservationSet() {
 		return condMiCalc.getAddedMoreThanOneObservationSet();

java/source/infodynamics/measures/continuous/gaussian/EntropyCalculatorMultiVariateGaussian.java

@@ -19,6 +19,7 @@
 package infodynamics.measures.continuous.gaussian;
 
 import infodynamics.measures.continuous.EntropyCalculatorMultiVariate;
+import infodynamics.measures.continuous.EntropyCalculatorMultiVariateCommon;
 import infodynamics.utils.MatrixUtils;
 
 /**
@@ -50,6 +51,7 @@ Theory' (John Wiley & Sons, New York, 1991).</li>
  * <a href="http://lizier.me/joseph/">www</a>)
  */
 public class EntropyCalculatorMultiVariateGaussian
+	extends EntropyCalculatorMultiVariateCommon 
 	implements EntropyCalculatorMultiVariate, Cloneable {
 
 	/**
@@ -63,84 +65,34 @@ public class EntropyCalculatorMultiVariateGaussian
 	 */
 	protected double[] means;
 	
-	/**
-	 * The set of observations, retained in case the user wants to retrieve the local
-	 *  entropy values of these
-	 */
-	protected double[][] observations;
-	
-	/**
-	 * Number of dimensions for our multivariate data
-	 */
-	protected int dimensions = 1;
-	
 	/**
 	 * Determinant of the covariance matrix; stored to save computation time
 	 */
 	protected double detCovariance = 0.0;
-
-	/**
-	 * Last average entropy we computed
-	 */
-	protected double lastAverage = 0;
-	
-	/**
-	 * Whether we are in debug mode
-	 */
-	protected boolean debug = false;
 	
 	/**
 	 * Construct an instance
 	 */
 	public EntropyCalculatorMultiVariateGaussian() {
-		// Nothing to do
+		super();
 	}
 	
-	@Override
-	public void initialise() throws Exception {
-		initialise(dimensions);
-	}
-
 	public void initialise(int dimensions) {
+		super.initialise(dimensions);
 		means = null;
 		L = null;
-		observations = null;
-		this.dimensions = dimensions;
 		detCovariance = 0;
-		lastAverage = 0.0;
 	}
 
-	/**
-	 * @throws Exception where the observations do not match the expected number of 
-	 *  dimensions, or covariance matrix is not positive definite (reflecting 
-	 *  redundant variables in the observations)
-	 */
 	@Override
-	public void setObservations(double[][] observations) throws Exception {
-		// Check that the observations was of the correct number of dimensions:
-		if (observations[0].length != dimensions) {
-			means = null;
-			L = null;
-			throw new Exception("Supplied observations does not match initialised number of dimensions");
-		}
+	public void finaliseAddObservations() throws Exception {
+		super.finaliseAddObservations();
 		means = MatrixUtils.means(observations);
+		double[][] originalObservations = observations; // Keep a reference as below
 		setCovariance(MatrixUtils.covarianceMatrix(observations, means));
 		// And keep a reference to the observations used here (must set this
 		//  *after* setCovariance, since setCovariance sets the observations to null
-		this.observations = observations;
-	}
-
-	/**
-	 * @throws Exception where the observations do not match the expected number of 
-	 *  dimensions, or covariance matrix is not positive definite (reflecting 
-	 *  redundant variables in the observations)
-	 */
-	@Override
-	public void setObservations(double[] observations) throws Exception {
-		if (dimensions != 1) {
-			throw new Exception(String.format("Cannot set univariate observations when expected dimension = %d", dimensions));
-		}
-		setObservations(MatrixUtils.reshape(observations, observations.length, 1));
+		observations = originalObservations;
 	}
 
 	/**
@@ -209,49 +161,45 @@ public class EntropyCalculatorMultiVariateGaussian
 	 */
 	@Override
 	public double computeAverageLocalOfObservations() {
+		if (isComputed) {
+			return lastAverage;
+		}
 		// Simple way:
 		// detCovariance = MatrixUtils.determinantSymmPosDefMatrix(covariance);
 		// Using cached Cholesky decomposition:
 		detCovariance = MatrixUtils.determinantViaCholeskyResult(L);
 		lastAverage =  0.5 * (dimensions* (1 + Math.log(2.0*Math.PI)) +
 			Math.log(detCovariance));
+		isComputed = true;
 		return lastAverage;
 	}
 
-	@Override
-	public void setDebug(boolean debug) {
-		this.debug = debug;
-	}
-
+	/**
+	 * <p>Set properties for this calculator.
+	 *  New property values are not guaranteed to take effect until the next call
+	 *  to an initialise method.
+	 * 
+	 * <p>Valid property names, and what their
+	 * values should represent, include:</p>
+	 * <ul>
+	 * 		<li>{@link #NORMALISE_PROP_NAME} as per {@link EntropyCalculatorMultiVariate#setProperty()}
+	 *         however note that for this Gaussian estimator setting this to true would fix the entropy
+	 *         values.</li>
+	 *  	<li>any valid properties for {@link EntropyCalculatorMultiVariate#setProperty(String, String)}.</li>
+	 * </ul> 
+	 * 
+	 * <p>Unknown property values are ignored.</p>
+	 * 
+	 * @param propertyName name of the property
+	 * @param propertyValue value of the property
+	 * @throws Exception for invalid property values
+	 */
 	@Override
 	public void setProperty(String propertyName, String propertyValue)
 			throws Exception {
-		boolean propertySet = true;
-		if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
-			dimensions = Integer.parseInt(propertyValue);
-		} else {
-			// No property was set
-			propertySet = false;
-		}
-		if (debug && propertySet) {
-			System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
-					" to " + propertyValue);
-		}
-	}
-
-	@Override
-	public String getProperty(String propertyName) throws Exception {
-		if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
-			return Integer.toString(dimensions);
-		} else {
-			// No property was set, and no superclass to call:
-			return null;
-		}
-	}
-	
-	@Override
-	public double getLastAverage() {
-		return lastAverage;
+		// Actually don't need to implement this method, but want 
+		//  the javadocs to get generated to warn user about the normalise property.
+		super.setProperty(propertyName, propertyValue);
 	}
 
 	/**
@@ -316,9 +264,13 @@ public class EntropyCalculatorMultiVariateGaussian
 		if (observations == null) {
 			throw new Exception("Cannot compute local values since no observations were supplied");
 		}
-		return computeLocalUsingPreviousObservations(observations);
+		double[] localValues = computeLocalUsingPreviousObservations(observations);
+		lastAverage = MatrixUtils.mean(localValues);
+		isComputed = true;
+		return localValues;
 	}
 
+	@Override
 	public int getNumObservations() throws Exception {
 		if (observations == null) {
 			throw new Exception("Cannot return number of observations because either " +

java/source/infodynamics/measures/continuous/kernel/EntropyCalculatorMultiVariateKernel.java

@@ -19,7 +19,7 @@
 package infodynamics.measures.continuous.kernel;
 
 import infodynamics.measures.continuous.EntropyCalculatorMultiVariate;
-import infodynamics.utils.MatrixUtils;
+import infodynamics.measures.continuous.EntropyCalculatorMultiVariateCommon;
 
 /**
  * <p>Computes the differential entropy of a given set of observations
@@ -44,42 +44,13 @@ import infodynamics.utils.MatrixUtils;
  * @author Joseph Lizier (<a href="joseph.lizier at gmail.com">email</a>,
  * <a href="http://lizier.me/joseph/">www</a>)
  */
-public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMultiVariate {
+public class EntropyCalculatorMultiVariateKernel
+	extends EntropyCalculatorMultiVariateCommon 
+	implements EntropyCalculatorMultiVariate
+	{
 
 	private KernelEstimatorMultiVariate mvke = null;
-	/**
-	 * Number of observations supplied
-	 */
-	private int totalObservations = 0;
-	// private int dimensions = 0;
-	/**
-	 * Whether we're in debug mode
-	 */
-	private boolean debug = false;
-	/**
-	 * The supplied observations
-	 */
-	private double[][] observations = null;
-	/**
-	 * Last computed average entropy
-	 */
-	private double lastEntropy;
-	/**
-	 * Whether we normalise the incoming observations to mean 0,
-	 * standard deviation 1.
-	 */
-	private boolean normalise = true;
-	/**
-	 * Property for whether we normalise the incoming observations to mean 0,
-	 * standard deviation 1.
-	 */
-	public static final String NORMALISE_PROP_NAME = "NORMALISE";
 	
-	/**
-	 * Number of joint variables/dimensions
-	 */
-	protected int dimensions = 1;
-
 	/**
 	 * Default value for kernel width
 	 */
@@ -101,15 +72,10 @@ public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMul
 	 * Construct an instance
 	 */
 	public EntropyCalculatorMultiVariateKernel() {
+		super();
 		mvke = new KernelEstimatorMultiVariate();
 		mvke.setDebug(debug);
 		mvke.setNormalise(normalise);
-		lastEntropy = 0.0;
-	}
-
-	@Override
-	public void initialise() throws Exception {
-		initialise(dimensions);
 	}
 
 	public void initialise(int dimensions) {
@@ -128,34 +94,22 @@ public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMul
 	 *  standard deviations from the mean (otherwise it is an absolute value)
 	 */
 	public void initialise(int dimensions, double kernelWidth) {
+		super.initialise(dimensions);
 		this.kernelWidth = kernelWidth;
-		this.dimensions = dimensions;
 		mvke.initialise(dimensions, kernelWidth);
-		// this.dimensions = dimensions;
-		lastEntropy = 0.0;
 	}
 
 	@Override
-	public void setObservations(double observations[][]) {
+	public void finaliseAddObservations() throws Exception {
+		super.finaliseAddObservations();
 		mvke.setObservations(observations);
-		totalObservations = observations.length;
-		this.observations = observations;
-	}
-
-	/**
-	 * @throws Exception where the observations do not match the expected number of 
-	 *  dimensions
-	 */
-	@Override
-	public void setObservations(double[] observations) throws Exception {
-		if (dimensions != 1) {
-			throw new Exception(String.format("Cannot set univariate observations when expected dimension = %d", dimensions));
-		}
-		setObservations(MatrixUtils.reshape(observations, observations.length, 1));
 	}
 
 	@Override
 	public double computeAverageLocalOfObservations() {
+		if (isComputed) {
+			return lastAverage;
+		}
 		double entropy = 0.0;
 		for (int b = 0; b < totalObservations; b++) {
 			double prob = mvke.getProbability(observations[b]);
@@ -165,8 +119,9 @@ public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMul
 				System.out.println(b + ": " + prob + " -> " + (-cont/Math.log(2.0)) + " -> sum: " + (entropy/Math.log(2.0)));
 			}
 		}
-		lastEntropy = entropy / (double) totalObservations / Math.log(2.0);
-		return lastEntropy;
+		lastAverage = entropy / (double) totalObservations / Math.log(2.0);
+		isComputed = true;
+		return lastAverage;
 	}
 
 	@Override
@@ -204,22 +159,18 @@ public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMul
 		}
 		entropy = entropy / (double) totalObservations / Math.log(2.0); // Don't use /= as I'm not sure of the order of operations it applies.
 		if (isPreviousObservations) {
-			lastEntropy = entropy; 
+			lastAverage = entropy;
+			isComputed = true;
 		}
 		return localEntropy;
 	}
 
 	@Override
 	public void setDebug(boolean debug) {
-		this.debug = debug;
+		super.setDebug(debug);
 		mvke.setDebug(debug);
 	}
 
-	@Override
-	public double getLastAverage() {
-		return lastEntropy;
-	}
-	
 	/**
 	 * <p>Set properties for the kernel entropy calculator.
 	 *  New property values are not guaranteed to take effect until the next call
@@ -232,8 +183,6 @@ public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMul
 	 * 			kernel width to be used in the calculation. If {@link #normalise} is set,
 	 * 		    then this is a number of standard deviations; otherwise it
 	 * 			is an absolute value. Default is {@link #DEFAULT_KERNEL_WIDTH}.</li>
-	 * 		<li>{@link #NORMALISE_PROP_NAME} -- whether to normalise the incoming variable values
-	 * 			to mean 0, standard deviation 1, or not (default false). Sets {@link #normalise}.</li>
 	 *  	<li>any valid properties for {@link EntropyCalculatorMultiVariate#setProperty(String, String)}.</li>
 	 * </ul> 
 	 * 
@@ -254,14 +203,15 @@ public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMul
 		if (propertyName.equalsIgnoreCase(KERNEL_WIDTH_PROP_NAME) ||
 				propertyName.equalsIgnoreCase(EPSILON_PROP_NAME)) {
 			kernelWidth = Double.parseDouble(propertyValue);
-		} else if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
-			normalise = Boolean.parseBoolean(propertyValue);
-			mvke.setNormalise(normalise);
-		} else if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
-			dimensions = Integer.parseInt(propertyValue);
 		} else {
 			// No property was set
 			propertySet = false;
+			// try the superclass:
+			super.setProperty(propertyName, propertyValue);
+			if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
+				// Handle an additional step for this one:
+				mvke.setNormalise(normalise);
+			}
 		}
 		if (debug && propertySet) {
 			System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
@@ -274,19 +224,10 @@ public class EntropyCalculatorMultiVariateKernel implements EntropyCalculatorMul
 		if (propertyName.equalsIgnoreCase(KERNEL_WIDTH_PROP_NAME) ||
 				propertyName.equalsIgnoreCase(EPSILON_PROP_NAME)) {
 			return Double.toString(kernelWidth);
-		} else if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
-			return Boolean.toString(normalise);
-		} else if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
-			return Integer.toString(dimensions);
 		} else {
-			// No property was set, and no superclass to call:
-			return null;
+		      // try the superclass:
+		      return super.getProperty(propertyName);
 		}
 	}
 
-	@Override
-	public int getNumObservations() throws Exception {
-		return totalObservations;
-	}
-
 }

java/source/infodynamics/measures/continuous/kozachenko/EntropyCalculatorMultiVariateKozachenko.java

@@ -18,13 +18,10 @@
 
 package infodynamics.measures.continuous.kozachenko;
 
-import java.util.Random;
-import java.util.Vector;
-
 import infodynamics.measures.continuous.EntropyCalculatorMultiVariate;
+import infodynamics.measures.continuous.EntropyCalculatorMultiVariateCommon;
 import infodynamics.utils.EuclideanUtils;
 import infodynamics.utils.MathsUtils;
-import infodynamics.utils.MatrixUtils;
 
 /**
  * <p>Computes the differential entropy of a given set of observations
@@ -59,59 +56,11 @@ import infodynamics.utils.MatrixUtils;
  * @author Joseph Lizier (<a href="joseph.lizier at gmail.com">email</a>,
  * <a href="http://lizier.me/joseph/">www</a>)
  */
-public class EntropyCalculatorMultiVariateKozachenko  
+public class EntropyCalculatorMultiVariateKozachenko
+	extends EntropyCalculatorMultiVariateCommon 
 	implements EntropyCalculatorMultiVariate {
 
-	/**
-	 * Total number of observations supplied.
-	 */
-	private int totalObservations = 0;
-	/**
-	 * Number of dimensions of our multivariate data set
-	 */
-	private int dimensions = 1;
-	/**
-	 * The set of observations, retained in case the user wants to retrieve the local
-	 *  entropy values of these.
-	 */
-	protected double[][] rawData;
-	/**
-	 * Store the last computed average H
-	 */
-	private double lastAverage = 0.0;
-	/**
-	 * Store the last computed local H
-	 */
-	private double[] lastLocalEntropy;
-	/**
-	 * Track whether we've computed the average for the supplied
-	 *  observations yet
-	 */
-	private boolean isComputed;
-  /**
-	 * Storage for observations supplied via {@link #addObservations(double[][])}
-	 * type calls
-	 */
-	protected Vector<double[][]> vectorOfObservations;
-	/**
-	 * Whether to report debug messages or not
-	 */
-	protected boolean debug = false;
 
-  /**
-   * Property name for an amount of random Gaussian noise to be
-   *  added to the data (default is 1e-8, matching the MILCA toolkit).
-   */
-  public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD";
-
-  /**
-   * Whether to add an amount of random noise to the incoming data
-   */
-  protected boolean addNoise = true;
-  /**
-   * Amount of random Gaussian noise to add to the incoming data
-   */
-  protected double noiseLevel = (double) 1e-8;
   /**
    * Stored pre-computed value of the Euler-Mascheroni constant
    */
@@ -121,171 +70,9 @@ public class EntropyCalculatorMultiVariateKozachenko
 	 * Construct an instance
 	 */
 	public EntropyCalculatorMultiVariateKozachenko() {
-		totalObservations = 0;
-		isComputed = false;
-		lastLocalEntropy = null;
-	}
-
-	@Override
-	public void initialise() throws Exception {
-		initialise(dimensions);
-	}
-
-	public void initialise(int dimensions) {
-		this.dimensions = dimensions;
-		rawData = null;
-		totalObservations = 0;
-		isComputed = false;
-		lastLocalEntropy = null;
-	}
-
-	@Override
-	public void setProperty(String propertyName, String propertyValue)
-			throws Exception {
-		boolean propertySet = true;
-		if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
-			dimensions = Integer.parseInt(propertyValue);
-    } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
-      if (propertyValue.equals("0") ||
-          propertyValue.equalsIgnoreCase("false")) {
-        addNoise = false;
-        noiseLevel = 0;
-      } else {
-        addNoise = true;
-        noiseLevel = Double.parseDouble(propertyValue);
-      }
-		} else {
-			// No property was set, and no superclass to call.
-			propertySet = false;
-		}
-		if (debug && propertySet) {
-			System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
-					" to " + propertyValue);
-		}
-	}
-
-	@Override
-	public String getProperty(String propertyName) throws Exception {
-		if (propertyName.equalsIgnoreCase(NUM_DIMENSIONS_PROP_NAME)) {
-			return Integer.toString(dimensions);
-    } else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
-      return Double.toString(noiseLevel);
-		} else {
-			// No property was set, and no superclass to call.
-			return null;
-		}
-	}
-
-	public void startAddObservations() {
-		isComputed = false;
-    totalObservations = 0;
-		lastLocalEntropy = null;
-    rawData = null;
-		vectorOfObservations = new Vector<double[][]>();
-	}
-
-	public void finaliseAddObservations() {
-
-		rawData = new double[totalObservations][dimensions];
-		
-		// Construct the joint vectors from the given observations
-		//  (removing redundant data which is outside any timeDiff)
-		int startObservation = 0;
-		for (double[][] obs : vectorOfObservations) {
-			// Copy the data from these given observations into our master array
-			MatrixUtils.arrayCopy(obs, 0, 0,
-					rawData, startObservation, 0,
-					obs.length, dimensions);
-			startObservation += obs.length;
-		}
-
-		// We don't need to keep the vector of observation sets anymore:
-		vectorOfObservations = null;
-
-    if (addNoise) {
-      Random random = new Random();
-      // Add Gaussian noise of std dev noiseLevel to the data
-      for (int r = 0; r < totalObservations; r++) {
-        for (int c = 0; c < dimensions; c++) {
-          rawData[r][c] += random.nextGaussian()*noiseLevel;
-        }
-      }
-    }
-	}
-
-	@Override
-  public void setObservations(double[][] observations) {
-    startAddObservations();
-    addObservations(observations);
-    finaliseAddObservations();
-  }
-
-  public void setObservations(double[][] observations1, double[][] observations2)
-      throws Exception {
-    startAddObservations();
-    addObservations(observations1, observations2);
-    finaliseAddObservations();
-  }
-
-	/*
-	 * (non-Javadoc)
-	 * @see infodynamics.measures.continuous.EntropyCalculator#setObservations(double[])
-	 * 
-	 * This method here to ensure we make compatibility with the 
-	 *  EntropyCalculator interface.
-	 */
-	@Override
-  public void setObservations(double[] observations) {
-    startAddObservations();
-    addObservations(observations);
-    finaliseAddObservations();
-  }
-
-  public void addObservations(double[][] observations) {
-		if (vectorOfObservations == null) {
-			// startAddObservations was not called first
-			throw new RuntimeException("User did not call startAddObservations before addObservations");
-		}
-    vectorOfObservations.add(observations);
-    totalObservations += observations.length;
-  }
-
-	public void addObservations(double[] observations) {
-		rawData = MatrixUtils.reshape(observations, observations.length, 1);
-		addObservations(rawData);
-	}
-	
-	/**
-	 * Each row of the data is an observation; each column of
-	 *  the row is a new variable in the multivariate observation.
-	 * This method signature allows the user to call setObservations for
-	 *  joint time series without combining them into a single joint time
-	 *  series (we do the combining for them).
-	 * 
-	 * @param data1 first few variables in the joint data
-	 * @param data2 the other variables in the joint data
-	 * @throws Exception When the length of the two arrays of observations do not match.
-	 * @see #addObservations(double[][])
-	 */
-	public void addObservations(double[][] data1,
-			double[][] data2) throws Exception {
-		int timeSteps = data1.length;
-		if ((data1 == null) || (data2 == null)) {
-			throw new Exception("Cannot have null data arguments");
-		}
-		if (data1.length != data2.length) {
-			throw new Exception("Length of data1 (" + data1.length + ") is not equal to the length of data2 (" +
-					data2.length + ")");
-		}
-		int data1Variables = data1[0].length;
-		int data2Variables = data2[0].length;
-		double[][] data = new double[timeSteps][data1Variables + data2Variables];
-		for (int t = 0; t < timeSteps; t++) {
-			System.arraycopy(data1[t], 0, data[t], 0, data1Variables);
-			System.arraycopy(data2[t], 0, data[t], data1Variables, data2Variables);
-		}
-		// Now defer to the normal setObservations method
-		addObservations(data);
+		super();
+		noiseLevel = (double) 1e-8; // Default to align with KSG estimators
+	    addNoise = true;
 	}
 
 	/**
@@ -298,12 +85,12 @@ public class EntropyCalculatorMultiVariateKozachenko
 		}
 		double sdTermHere = sdTerm(totalObservations, dimensions);
 		double emConstHere = eulerMascheroniTerm(totalObservations);
-		double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(rawData);
+		double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(observations);
 		double entropy = 0.0;
 		if (debug) {
 			System.out.println("t,\tminDist,\tlogMinDist,\tsum");
 		}
-		for (int t = 0; t < rawData.length; t++) {
+		for (int t = 0; t < observations.length; t++) {
 			entropy += Math.log(2.0 * minDistance[t]);
 			if (debug) {
 				System.out.println(t + ",\t" + 
@@ -332,21 +119,18 @@ public class EntropyCalculatorMultiVariateKozachenko
 	 */
 	@Override
 	public double[] computeLocalOfPreviousObservations() {
-		if (lastLocalEntropy != null) {
-			return lastLocalEntropy;
-		}
 
 		double sdTermHere = sdTerm(totalObservations, dimensions);
 		double emConstHere = eulerMascheroniTerm(totalObservations);
 		double constantToAddIn = sdTermHere + emConstHere;
 		
-		double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(rawData);
+		double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(observations);
 		double entropy = 0.0;
-		double[] localEntropy = new double[rawData.length];
+		double[] localEntropy = new double[observations.length];
 		if (debug) {
 			System.out.println("t,\tminDist,\tlogMinDist,\tlocal,\tsum");
 		}
-		for (int t = 0; t < rawData.length; t++) {
+		for (int t = 0; t < observations.length; t++) {
 			localEntropy[t] = Math.log(2.0 * minDistance[t]) * (double) dimensions;
 			// using natural units
 			// localEntropy[t] /= Math.log(2);
@@ -362,7 +146,7 @@ public class EntropyCalculatorMultiVariateKozachenko
 		}
 		entropy /= (double) totalObservations;
 		lastAverage = entropy;
-		lastLocalEntropy = localEntropy;
+		isComputed = true;
 		return localEntropy;
 	}
 
@@ -375,13 +159,13 @@ public class EntropyCalculatorMultiVariateKozachenko
 		double emConstHere = eulerMascheroniTerm(totalObservations);
 		double constantToAddIn = sdTermHere + emConstHere;
 		
-		double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(rawData);
+		double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(observations);
 		double entropy = 0.0;
-		double[] localEntropy = new double[rawData.length];
+		double[] localEntropy = new double[observations.length];
 		if (debug) {
 			System.out.println("t,\tminDist,\tlogMinDist,\tlocal,\tsum");
 		}
-		for (int t = 0; t < rawData.length; t++) {
+		for (int t = 0; t < observations.length; t++) {
 			localEntropy[t] = Math.log(2.0 * minDistance[t]) * (double) dimensions;
 			// using natural units
 			// localEntropy[t] /= Math.log(2);
@@ -395,9 +179,6 @@ public class EntropyCalculatorMultiVariateKozachenko
 						entropy);
 			}
 		}
-		entropy /= (double) totalObservations;
-		lastAverage = entropy;
-		lastLocalEntropy = localEntropy;
 		return localEntropy;
 	}
 
@@ -473,18 +254,4 @@ public class EntropyCalculatorMultiVariateKozachenko
 		return result;
 	}
 	
-	@Override
-	public void setDebug(boolean debug) {
-		this.debug = debug;
-	}
-
-	@Override
-	public double getLastAverage() {
-		return lastAverage;
-	}
-
-	@Override
-	public int getNumObservations() {
-		return totalObservations;
-	}
 }

java/source/infodynamics/measures/continuous/kraskov/ActiveInfoStorageCalculatorKraskov.java

@@ -23,6 +23,7 @@ import java.util.Hashtable;
 import infodynamics.measures.continuous.ActiveInfoStorageCalculator;
 import infodynamics.measures.continuous.ActiveInfoStorageCalculatorViaMutualInfo;
 import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
+import infodynamics.utils.MatrixUtils;
 
 /**
  * An Active Information Storage (AIS) calculator (implementing {@link ActiveInfoStorageCalculator})
@@ -230,4 +231,48 @@ public class ActiveInfoStorageCalculatorKraskov
 		}
 	}
 
+	/**
+     * Debug method to return the k nearest neighbour distances that 
+     *  would be utilised for each sample point here.
+     * Note that this is specifically the max-norm across the target-targetPast variables, which is used for each
+     *  range search in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+     *  
+     * @param startTimePoint
+     * @param numTimePoints
+     * @return
+     * @throws Exception
+     */
+	public double[] kNNDistances(int startTimePoint, int numTimePoints) throws Exception {
+		// Defer the call to the underlying KSG CMI estimator
+		return ((MutualInfoCalculatorMultiVariateKraskov) miCalc).kNNDistances(startTimePoint, numTimePoints);
+	}
+
+	/**
+	 * Debug method to return the k nearest neighbour distances that
+	 *  would be utilised in {@link #computeLocalUsingPreviousObservations(double[])}
+	 *  for a cross AIS.
+     * Note that this is specifically the max-norm across the target-targetPast variables, which is used for each
+     *  range search in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+	 *  
+	 * @param startTimePoint
+	 * @param numTimePoints
+	 * @param newObservations
+	 * @return
+	 * @throws Exception
+	 */
+	public double[] kNNDistancesForNewSamples(int startTimePoint, int numTimePoints,
+			double[] newObservations) throws Exception {
+		if (newObservations.length - (k-1)*tau - 1 <= 0) {
+			// There are no observations to compute for here
+			throw new Exception("Not enough samples to embed");
+		}
+		double[][] newDestPastVectors = 
+				MatrixUtils.makeDelayEmbeddingVector(newObservations, k, tau, (k-1)*tau, newObservations.length - (k-1)*tau - 1);
+		double[][] newDestNextVectors =
+				MatrixUtils.makeDelayEmbeddingVector(newObservations, 1, (k-1)*tau + 1, newObservations.length - (k-1)*tau - 1);
+		return ((MutualInfoCalculatorMultiVariateKraskov) miCalc).
+				kNNDistancesForNewSamples(startTimePoint, numTimePoints, newDestPastVectors, newDestNextVectors);
+	}
 }

java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java

@@ -19,6 +19,7 @@
 package infodynamics.measures.continuous.kraskov;
 
 import java.util.Calendar;
+import java.util.PriorityQueue;
 
 import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariate;
 import infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateCommon;
@@ -27,6 +28,7 @@ import infodynamics.utils.KdTree;
 import infodynamics.utils.MathsUtils;
 import infodynamics.utils.MatrixUtils;
 import infodynamics.utils.NearestNeighbourSearcher;
+import infodynamics.utils.NeighbourNodeData;
 import infodynamics.utils.UnivariateNearestNeighbourSearcher;
 import infodynamics.utils.EmpiricalMeasurementDistribution;
 import infodynamics.utils.NativeUtils;
@@ -313,13 +315,6 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
 		// Allow the parent to generate the data for us first
 		super.finaliseAddObservations();
 		
-		if (dynCorrExcl && addedMoreThanOneObservationSet) {
-			// We have not properly implemented dynamic correlation exclusion for
-			//  multiple observation sets, so throw an error
-			throw new RuntimeException("Addition of multiple observation sets is not currently " +
-					"supported with property " + PROP_DYN_CORR_EXCL_TIME + " set");
-		}
-		
 		if (totalObservations <= k + 2*dynCorrExclTime) {
 			throw new Exception("There are less observations provided (" +
 					totalObservations +
@@ -725,35 +720,41 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
     if (kdTreeJoint == null) {
       kdTreeJoint = new KdTree(
           new int[] {dimensionsVar1, dimensionsVar2, dimensionsCond},
-          new double[][][] {var1Observations, var2Observations, condObservations});
+          new double[][][] {var1Observations, var2Observations, condObservations},
+          observationSetIndices, observationTimePoints);
       kdTreeJoint.setNormType(normType);
     }
     if (dimensionsVar1 > 1) {
       if (kdTreeVar1Conditional == null) {
         kdTreeVar1Conditional = new KdTree(
             new int[] {dimensionsVar1, dimensionsCond},
-            new double[][][] {var1Observations, condObservations});
+            new double[][][] {var1Observations, condObservations},
+            observationSetIndices, observationTimePoints);
         kdTreeVar1Conditional.setNormType(normType);
       } 
     } else { // Univariate variable 1, so we'll search its space alone as this is faster
       if (uniNNSearcherVar1 == null) {
-        uniNNSearcherVar1 = new UnivariateNearestNeighbourSearcher(var1Observations);
+        uniNNSearcherVar1 = new UnivariateNearestNeighbourSearcher(var1Observations,
+                observationSetIndices, observationTimePoints);
       }
     }
     if (dimensionsVar2 > 1) {
       if (kdTreeVar2Conditional == null) {
         kdTreeVar2Conditional = new KdTree(
             new int[] {dimensionsVar2, dimensionsCond},
-            new double[][][] {var2Observations, condObservations});
+            new double[][][] {var2Observations, condObservations},
+            observationSetIndices, observationTimePoints);
         kdTreeVar2Conditional.setNormType(normType);
       }
     } else { // Univariate variable 2, so we'll search its space alone as this is faster
       if (uniNNSearcherVar2 == null) {
-        uniNNSearcherVar2 = new UnivariateNearestNeighbourSearcher(var2Observations);
+        uniNNSearcherVar2 = new UnivariateNearestNeighbourSearcher(var2Observations,
+                observationSetIndices, observationTimePoints);
       }
     }
     if ((nnSearcherConditional == null) && (dimensionsCond > 0)) {
-      nnSearcherConditional = NearestNeighbourSearcher.create(condObservations);
+      nnSearcherConditional = NearestNeighbourSearcher.create(condObservations,
+              observationSetIndices, observationTimePoints);
       nnSearcherConditional.setNormType(normType);
     }
   }
@@ -911,4 +912,136 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
 	// public ConditionalMutualInfoCalculatorMultiVariateKraskov clone() {
 	//	return this;
 	// }
+	
+    /**
+     * Debug method to return the k nearest neighbour distances that 
+     *  would be utilised for each sample point here.
+     * Note that this is specifically the max-norm across the three variables, which is used for each
+     *  range search in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+     *  
+     * @param startTimePoint
+     * @param numTimePoints
+     * @return
+     * @throws Exception
+     */
+	public double[] kNNDistances(int startTimePoint, int numTimePoints) throws Exception {
+		double[] kNNdistances = new double[numTimePoints];
+		
+		for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
+			// Compute eps for this time step by
+			//  finding the kth closest neighbour for point t:
+			PriorityQueue<NeighbourNodeData> nnPQ =
+					kdTreeJoint.findKNearestNeighbours(k, t, dynCorrExclTime);
+			// First element in the PQ is the kth NN,
+			//  and epsilon = kthNnData.distance
+			NeighbourNodeData kthNnData = nnPQ.poll();
+
+			kNNdistances[t - startTimePoint] = kthNnData.distance;
+		}
+		return kNNdistances;
+	}
+
+	/**
+	 * Debug method to return the k nearest neighbour distances that
+	 *  would be utilised in {@link #computeLocalUsingPreviousObservations(double[][], double[][], double[][])}
+	 *  for a cross conditional MI.
+     * Note that this is specifically the max-norm across the three variables, which is used for each
+     *  range search in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+	 *  
+	 * @param startTimePoint
+	 * @param numTimePoints
+	 * @param newVar1Observations
+	 * @param newVar2Observations
+	 * @return
+	 * @throws Exception
+	 */
+	public double[] kNNDistancesForNewSamples(int startTimePoint, int numTimePoints,
+			double[][] newStates1, double[][] newStates2, double[][] newCondStates) throws Exception {
+
+		// Do normalisation of the incoming data if required:
+		double[][] states1ToUse, states2ToUse, condStatesToUse;
+		if (normalise) {
+			states1ToUse = MatrixUtils.normaliseIntoNewArray(newStates1, var1Means, var1Stds);
+			states2ToUse = MatrixUtils.normaliseIntoNewArray(newStates2, var2Means, var2Stds);
+			if (dimensionsCond != 0) {
+				condStatesToUse = MatrixUtils.normaliseIntoNewArray(newCondStates, condMeans, condStds);
+			} else {
+				condStatesToUse = null;
+			}
+		} else {
+			states1ToUse = newStates1;
+			states2ToUse = newStates2;
+			condStatesToUse = newCondStates;
+		}
+		
+		double[] kNNdistances = new double[numTimePoints];
+		
+		for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
+			// Compute eps for this time step by
+			//  finding the kth closest neighbour for the new sample:
+			PriorityQueue<NeighbourNodeData> nnPQ =
+					kdTreeJoint.findKNearestNeighbours(k,
+							new double[][] {states1ToUse[t], states2ToUse[t], condStatesToUse[t]});
+			// First element in the PQ is the kth NN,
+			//  and epsilon = kthNnData.distance
+			NeighbourNodeData kthNnData = nnPQ.poll();
+			
+			kNNdistances[t - startTimePoint] = kthNnData.distance;
+		}
+		return kNNdistances;
+	}
+
+	/**
+	 * As per {@link #kNNDistancesForNewSamples(int, int, double[][], double[][], double[][])}
+	 * but for univariates
+	 *  
+	 * @param startTimePoint
+	 * @param numTimePoints
+	 * @param newVar1Observations
+	 * @param newVar2Observations
+	 * @return
+	 * @throws Exception
+	 */
+	public double[] kNNDistancesForNewSamples(int startTimePoint, int numTimePoints,
+			double[] newStates1, double[] newStates2, double[][] newCondStates) throws Exception {
+
+		if ((dimensionsVar1 != 1) || (dimensionsVar2 != 1)) {
+			throw new Exception("The number of source and dest dimensions (having been initialised to " +
+					dimensionsVar1 + " and " + dimensionsVar2 + ") can only be 1 when " +
+					"the univariate kNNDistancesForNewSamples(int, int, double[],double[],double[][]) " + 
+					"method is called");
+		}
+		return kNNDistancesForNewSamples(startTimePoint, numTimePoints,
+				MatrixUtils.reshape(newStates1, newStates1.length, 1),
+				MatrixUtils.reshape(newStates2, newStates2.length, 1),
+				newCondStates);
+	}
+
+	/**
+	 * As per {@link #kNNDistancesForNewSamples(int, int, double[][], double[][], double[][])}
+	 * but for univariates
+	 *  
+	 * @param startTimePoint
+	 * @param numTimePoints
+	 * @param newVar1Observations
+	 * @param newVar2Observations
+	 * @return
+	 * @throws Exception
+	 */
+	public double[] kNNDistancesForNewSamples(int startTimePoint, int numTimePoints,
+			double[] newStates1, double[] newStates2, double[] newCondStates) throws Exception {
+
+		if ((dimensionsVar1 != 1) || (dimensionsVar2 != 1) || (dimensionsCond != 1)) {
+			throw new Exception("The number of source, dest and conditional dimensions (having been initialised to " +
+					dimensionsVar1 + ", " + dimensionsVar2 + " and " + dimensionsCond + ") can only be 1 when " +
+					"the univariate kNNDistancesForNewSamples(int, int, double[],double[],double[]) " + 
+					"method is called");
+		}
+		return kNNDistancesForNewSamples(startTimePoint, numTimePoints,
+				MatrixUtils.reshape(newStates1, newStates1.length, 1),
+				MatrixUtils.reshape(newStates2, newStates2.length, 1),
+				MatrixUtils.reshape(newCondStates, newCondStates.length, 1));
+	}
 }

java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java

@@ -18,19 +18,13 @@
 
 package infodynamics.measures.continuous.kraskov;
 
-import java.util.Arrays;
 import java.util.Calendar;
 import java.util.PriorityQueue;
 
 import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariate;
-import infodynamics.utils.EmpiricalMeasurementDistribution;
-import infodynamics.utils.FirstIndexComparatorDouble;
-import infodynamics.utils.KdTree;
 import infodynamics.utils.MathsUtils;
 import infodynamics.utils.MatrixUtils;
-import infodynamics.utils.NearestNeighbourSearcher;
 import infodynamics.utils.NeighbourNodeData;
-import infodynamics.utils.UnivariateNearestNeighbourSearcher;
 
 /**
  * <p>Computes the differential conditional mutual information of two multivariate
@@ -340,7 +334,7 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1
 
 		for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
 			// Compute eps for this time step by
-			//  finding the kth closest neighbour for point t:
+			//  finding the kth closest neighbour for the new sample:
 			long methodStartTime = Calendar.getInstance().getTimeInMillis();
 			PriorityQueue<NeighbourNodeData> nnPQ =
 					kdTreeJoint.findKNearestNeighbours(k,

java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java

@@ -295,13 +295,6 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
     // Allow the parent to generate the data for us first
     super.finaliseAddObservations();
     
-    if (dynCorrExcl && addedMoreThanOneObservationSet) {
-      // We have not properly implemented dynamic correlation exclusion for
-      //  multiple observation sets, so throw an error
-      throw new RuntimeException("Addition of multiple observation sets is not currently " +
-          "supported with property " + PROP_DYN_CORR_EXCL_TIME + " set");
-    }
-    
     if (totalObservations <= k + 2*dynCorrExclTime) {
       throw new Exception("There are less observations provided (" +
           totalObservations +
@@ -729,15 +722,18 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
     //  the tree construction itself afterwards can't really be well parallelised.
     if (kdTreeJoint == null) {
       kdTreeJoint = new KdTree(new int[] {dimensionsSource, dimensionsDest},
-            new double[][][] {sourceObservations, destObservations});
+            new double[][][] {sourceObservations, destObservations},
+            observationSetIndices, observationTimePoints);
       kdTreeJoint.setNormType(normType);
     }
     if (nnSearcherSource == null) {
-      nnSearcherSource = NearestNeighbourSearcher.create(sourceObservations);
+      nnSearcherSource = NearestNeighbourSearcher.create(sourceObservations,
+    		  observationSetIndices, observationTimePoints);
       nnSearcherSource.setNormType(normType);
     }
     if (nnSearcherDest == null) {
-      nnSearcherDest = NearestNeighbourSearcher.create(destObservations);
+      nnSearcherDest = NearestNeighbourSearcher.create(destObservations,
+    		  observationSetIndices, observationTimePoints);
       nnSearcherDest.setNormType(normType);
     }
   }
@@ -913,7 +909,8 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
     public static final int INDEX_SUM_NY = 2;
     public static final int INDEX_SUM_2X_NEIGH_DIST = 3; // Not used by GPU, and only read for conditional entropy
     public static final int RETURN_ARRAY_LENGTH_MI = 3;
-    public static final int RETURN_ARRAY_LENGTH = 4;
+    @SuppressWarnings("unused") // Might be used by caller later
+	public static final int RETURN_ARRAY_LENGTH = 4;
 
     public MiKraskovThreadRunner(
         MutualInfoCalculatorMultiVariateKraskov miCalc,
@@ -1110,5 +1107,110 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
     }
     
     return totalErrors;
-  } 
+  }
+
+    /**
+     * Debug method to return the k nearest neighbour distances that 
+     *  would be utilised for each sample point here.
+     * Note that this is specifically the max-norm across the two variables, which is used for both
+     *  variables in the range searches in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+     *  
+     * @param startTimePoint
+     * @param numTimePoints
+     * @return
+     * @throws Exception
+     */
+	public double[] kNNDistances(int startTimePoint, int numTimePoints) throws Exception {
+		
+		double[] kNNdistances = new double[numTimePoints];
+		
+		for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
+			// Compute eps for this time step by
+			//  finding the kth closest neighbour for point t:
+			PriorityQueue<NeighbourNodeData> nnPQ =
+					kdTreeJoint.findKNearestNeighbours(k, t, dynCorrExclTime);
+			// First element in the PQ is the kth NN,
+			//  and epsilon = kthNnData.distance
+			NeighbourNodeData kthNnData = nnPQ.poll();
+
+			kNNdistances[t - startTimePoint] = kthNnData.distance;
+		}
+		return kNNdistances;
+	}
+	
+
+	/**
+	 * Debug method to return the k nearest neighbour distances that
+	 *  would be utilised in {@link #computeLocalUsingPreviousObservations(double[][], double[][])}
+	 *  for a cross MI.
+     * Note that this is specifically the max-norm across the two variables, which is used for both
+     *  variables in the range searches in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+	 *  
+	 * @param startTimePoint
+	 * @param numTimePoints
+	 * @param newVar1Observations
+	 * @param newVar2Observations
+	 * @return
+	 * @throws Exception
+	 */
+	public double[] kNNDistancesForNewSamples(int startTimePoint, int numTimePoints,
+			double[][] newVar1Observations, double[][] newVar2Observations) throws Exception {
+
+		double[][] states1ToUse, states2ToUse;
+		if (normalise) {
+			states1ToUse = MatrixUtils.normaliseIntoNewArray(newVar1Observations, sourceMeansBeforeNorm, sourceStdsBeforeNorm, 0, newVar1Observations.length-timeDiff);
+			states2ToUse = MatrixUtils.normaliseIntoNewArray(newVar2Observations, destMeansBeforeNorm, destStdsBeforeNorm, timeDiff, newVar2Observations.length-timeDiff);
+		} else {
+			if (timeDiff > 0) {
+				states1ToUse = MatrixUtils.selectRows(newVar1Observations, 0, newVar1Observations.length-timeDiff);
+				states2ToUse = MatrixUtils.selectRows(newVar2Observations, 0, newVar2Observations.length-timeDiff);
+			} else {
+				states1ToUse = newVar1Observations;
+				states2ToUse = newVar2Observations;
+			}
+		}
+
+		double[] kNNdistances = new double[numTimePoints];
+		
+		for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
+			// Compute eps for this time step by
+			//  finding the kth closest neighbour for the new sample:
+			PriorityQueue<NeighbourNodeData> nnPQ =
+					kdTreeJoint.findKNearestNeighbours(k,
+							new double[][] {states1ToUse[t], states2ToUse[t]});
+			// First element in the PQ is the kth NN,
+			//  and epsilon = kthNnData.distance
+			NeighbourNodeData kthNnData = nnPQ.poll();
+			
+			kNNdistances[t - startTimePoint] = kthNnData.distance;
+		}
+		return kNNdistances;
+	}
+
+	/**
+	 * As per {@link #kNNDistancesForNewSamples(int, int, double[][], double[][])}
+	 * but for univariates
+	 *  
+	 * @param startTimePoint
+	 * @param numTimePoints
+	 * @param newVar1Observations
+	 * @param newVar2Observations
+	 * @return
+	 * @throws Exception
+	 */
+	public double[] kNNDistancesForNewSamples(int startTimePoint, int numTimePoints,
+			double[] newVar1Observations, double[] newVar2Observations) throws Exception {
+
+		if ((dimensionsSource != 1) || (dimensionsDest != 1)) {
+			throw new Exception("The number of source and dest dimensions (having been initialised to " +
+					dimensionsSource + " and " + dimensionsDest + ") can only be 1 when " +
+					"the univariate kNNDistancesForNewSamples(int, int, double[],double[]) " + 
+					"method is called");
+		}
+		return kNNDistancesForNewSamples(startTimePoint, numTimePoints,
+				MatrixUtils.reshape(newVar1Observations, newVar1Observations.length, 1),
+				MatrixUtils.reshape(newVar2Observations, newVar2Observations.length, 1));
+	}
 }

java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java

@@ -72,7 +72,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1
 		double sumNx = 0;
 		double sumNy = 0;
 		double sum2xkNNDist = 0;
-				
+		
 		for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
 			// Compute eps for this time step by
 			//  finding the kth closest neighbour for point t:
@@ -196,7 +196,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1
 		
 		for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
 			// Compute eps for this time step by
-			//  finding the kth closest neighbour for point t:
+			//  finding the kth closest neighbour for the new sample:
 			PriorityQueue<NeighbourNodeData> nnPQ =
 					kdTreeJoint.findKNearestNeighbours(k,
 							new double[][] {newVar1Observations[t], newVar2Observations[t]});

java/source/infodynamics/measures/continuous/kraskov/TransferEntropyCalculatorKraskov.java

@@ -21,6 +21,7 @@ package infodynamics.measures.continuous.kraskov;
 import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariate;
 import infodynamics.measures.continuous.TransferEntropyCalculator;
 import infodynamics.measures.continuous.TransferEntropyCalculatorViaCondMutualInfo;
+import infodynamics.utils.MatrixUtils;
 
 /**
  * <p>Computes the differential transfer entropy (TE) between two univariate
@@ -250,4 +251,64 @@ public class TransferEntropyCalculatorKraskov
 			return super.getProperty(propertyName);
 		}
 	}
+	
+	/**
+     * Debug method to return the k nearest neighbour distances that 
+     *  would be utilised for each sample point here.
+     * Note that this is specifically the max-norm across the source-target-targetPast variables, which is used for each
+     *  range search in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+     *  
+     * @param startTimePoint
+     * @param numTimePoints
+     * @return
+     * @throws Exception
+     */
+	public double[] kNNDistances(int startTimePoint, int numTimePoints) throws Exception {
+		// Defer the call to the underlying KSG CMI estimator
+		return ((ConditionalMutualInfoCalculatorMultiVariateKraskov) condMiCalc).kNNDistances(startTimePoint, numTimePoints);
+	}
+	
+	/**
+	 * Debug method to return the k nearest neighbour distances that
+	 *  would be utilised in {@link #computeLocalUsingPreviousObservations(double[], double[])}
+	 *  for a cross TE.
+     * Note that this is specifically the max-norm across the source-target-targetPast variables, which is used for each
+     *  range search in algorithm 1 (although algorithm 2 would use the max distance
+     *  for each variable within the kNNs in their separate range searches). 
+	 *  
+	 * @param startTimePoint
+	 * @param numTimePoints
+	 * @param newSourceObservations
+	 * @param newDestObservations
+	 * @return
+	 * @throws Exception
+	 */
+	public double[] kNNDistancesForNewSamples(int startTimePoint, int numTimePoints,
+			double[] newSourceObservations, double[] newDestObservations) throws Exception {
+		if (newSourceObservations.length != newDestObservations.length) {
+			throw new Exception(String.format("Source and destination lengths (%d and %d) must match!",
+					newSourceObservations.length, newDestObservations.length));
+		}
+		if (newDestObservations.length < startTimeForFirstDestEmbedding + 2) {
+			// There are no observations to compute for here
+			throw new Exception("Not enough samples to embed");
+		}
+		// Now embed as per computeLocalUsingPreviousObservations() in super:
+		double[][] newDestPastVectors = 
+				MatrixUtils.makeDelayEmbeddingVector(newDestObservations, k, k_tau,
+						startTimeForFirstDestEmbedding,
+						newDestObservations.length - startTimeForFirstDestEmbedding - 1);
+		double[][] newDestNextVectors =
+				MatrixUtils.makeDelayEmbeddingVector(newDestObservations, 1,
+						startTimeForFirstDestEmbedding + 1,
+						newDestObservations.length - startTimeForFirstDestEmbedding - 1);
+		double[][] newSourcePastVectors = 
+				MatrixUtils.makeDelayEmbeddingVector(newSourceObservations, l, l_tau,
+						startTimeForFirstDestEmbedding + 1 - delay,
+						newSourceObservations.length - startTimeForFirstDestEmbedding - 1);
+		return ((ConditionalMutualInfoCalculatorMultiVariateKraskov) condMiCalc).
+				kNNDistancesForNewSamples(startTimePoint, numTimePoints, newSourcePastVectors, newDestNextVectors, newDestPastVectors);
+	}
+
 }

java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java

@@ -1,7 +1,6 @@
 package infodynamics.measures.spiking.integration;
 
 import java.util.Arrays;
-import java.util.Collections;
 import java.util.Iterator;
 import java.util.PriorityQueue;
 import java.util.Random;
@@ -16,8 +15,6 @@ import infodynamics.utils.KdTree;
 import infodynamics.utils.MathsUtils;
 import infodynamics.utils.MatrixUtils;
 import infodynamics.utils.NeighbourNodeData;
-import infodynamics.utils.FirstIndexComparatorDouble;
-import infodynamics.utils.UnivariateNearestNeighbourSearcher;
 import infodynamics.utils.EuclideanUtils;
 import infodynamics.utils.ParsedProperties;
 
@@ -103,6 +100,9 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 	protected KdTree kdTreeConditioningAtSpikes = null;
 	protected KdTree kdTreeConditioningAtSamples = null;
 
+	// Cache these to return with the local values:
+	protected double[][] arrayedJointEmbeddingsFromSpikes;
+	
 	public static final String KNNS_PROP_NAME = "Knns";
 
 	/**
@@ -400,7 +400,6 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 		if (vectorOfConditionalSpikeTimes.size() > 0) {
 			conditionalIterator = vectorOfConditionalSpikeTimes.iterator();
 		}
-		int timeSeriesIndex = 0;
 		for (double[] destSpikeTimes : vectorOfDestinationSpikeTimes) {
 			double[] sourceSpikeTimes = sourceIterator.next();
 			double[][] conditionalSpikeTimes = null;
@@ -417,13 +416,15 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 
 
 		// Convert the vectors to arrays so that they can be put in the trees
+		// Embeddings from target spikes:
 		double[][] arrayedTargetEmbeddingsFromSpikes = new double[conditioningEmbeddingsFromSpikes.size()][numDestPastIntervals + numCondPastIntervals];
-		double[][] arrayedJointEmbeddingsFromSpikes = new double[conditioningEmbeddingsFromSpikes.size()][numDestPastIntervals +
+		arrayedJointEmbeddingsFromSpikes = new double[conditioningEmbeddingsFromSpikes.size()][numDestPastIntervals +
 														  numCondPastIntervals + numSourcePastIntervals];
 		for (int i = 0; i < conditioningEmbeddingsFromSpikes.size(); i++) {
 			arrayedTargetEmbeddingsFromSpikes[i] = conditioningEmbeddingsFromSpikes.elementAt(i);
 			arrayedJointEmbeddingsFromSpikes[i] = jointEmbeddingsFromSpikes.elementAt(i);
 		}
+		// Sample points:
 		double[][] arrayedTargetEmbeddingsFromSamples = new double[conditioningEmbeddingsFromSamples.size()][numDestPastIntervals + numCondPastIntervals];
 		double[][] arrayedJointEmbeddingsFromSamples = new double[conditioningEmbeddingsFromSamples.size()][numDestPastIntervals +
 														  numCondPastIntervals + numSourcePastIntervals];
@@ -697,19 +698,51 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 		return new distanceAndNumPoints(maxDistance, i);
 	}
 
-	@Override
-	public double computeAverageLocalOfObservations() throws Exception {
-		return computeAverageLocalOfObservations(kdTreeJointAtSpikes, jointEmbeddingsFromSpikes);
+	/**
+	 * Data structure to store the interim results from a TE calculation 
+	 */
+	protected class TERateResults {
+		protected double[] localContributions;
+		protected double meanRate;
+		protected double totalTime;
+		protected int numTargetSpikes;
+		
+		protected TERateResults(double meanRate, double[] localValues, int numTargetSpikes, double totalTime) {
+			this.meanRate = meanRate;
+			this.localContributions = localValues;
+			this.numTargetSpikes = numTargetSpikes;
+			this.totalTime = totalTime;
+		}
+	}
+	
+	@Override
+	public double computeAverageLocalOfObservations() throws Exception {
+		TERateResults teRateResults = computeAverageLocalOfObservations(kdTreeJointAtSpikes, jointEmbeddingsFromSpikes);
+		return teRateResults.meanRate;
 	}
 	
-
 	/*
 	 * We take the actual joint tree at spikes (along with the associated embeddings) as an argument, as we will need to swap these out when
 	 * computing surrogates.
 	 */
-	public double computeAverageLocalOfObservations(KdTree actualKdTreeJointAtSpikes, Vector<double[]> actualJointEmbeddingsFromSpikes) throws Exception {
+	protected TERateResults computeAverageLocalOfObservations(KdTree actualKdTreeJointAtSpikes, Vector<double[]> actualJointEmbeddingsFromSpikes) throws Exception {
 
 		double currentSum = 0;
+		int numTargetSpikesToSkip = MatrixUtils.max(destPastIntervals);
+		
+		double[] localValues;
+		int indexInLocals = 0;
+		if (processTimeLengths.size() == 1) {
+			// Pad the local values for the first spikes (which we can't compute a TE contribution for) to zero
+			localValues = new double[conditioningEmbeddingsFromSpikes.size() + numTargetSpikesToSkip];
+			indexInLocals = numTargetSpikesToSkip;
+		} else {
+			// Don't pad with zeros, would need to properly return locals for each separate spike train
+			// TODO Do this separately for each separate spike train when we implement in a similar fashion
+			//  for regular time series estimators
+			localValues = new double[conditioningEmbeddingsFromSpikes.size()];
+		}
+		
 		for (int i = 0; i < conditioningEmbeddingsFromSpikes.size(); i++) {
 
 			double radiusJointSpikes = actualKdTreeJointAtSpikes.findKNearestNeighbours(Knns, i).poll().norms[0];
@@ -812,16 +845,18 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 				radiusConditioningSamples = Math.sqrt(radiusConditioningSamples);
 			}
 			
-			currentSum += (MathsUtils.digamma(kJointSpikes) - MathsUtils.digamma(kJointSamples) + 
+			localValues[indexInLocals] = (MathsUtils.digamma(kJointSpikes) - MathsUtils.digamma(kJointSamples) + 
 				       ((numDestPastIntervals + numCondPastIntervals + numSourcePastIntervals) * (-Math.log(radiusJointSpikes) + Math.log(radiusJointSamples))) -
 				       MathsUtils.digamma(kConditioningSpikes) + MathsUtils.digamma(kConditioningSamples) + 
-				       + ((numDestPastIntervals + numCondPastIntervals) * (Math.log(radiusConditioningSpikes) - Math.log(radiusConditioningSamples))));
-			if (Double.isNaN(currentSum)) {
+				       + ((numDestPastIntervals + numCondPastIntervals) * (Math.log(radiusConditioningSpikes) - Math.log(radiusConditioningSamples)))); 
+			currentSum += localValues[indexInLocals];
+			if (Double.isNaN(localValues[indexInLocals])) {
 				for (double[] embed : jointEmbeddingsFromSamples) {
 					System.out.println(Arrays.toString(embed));
 				}
 				throw new Exception("NaNs in TE clac " + radiusJointSpikes + " " + radiusJointSamples + " " + tempRadiusJointSamples);
 			}
+			indexInLocals++; // Move to record the next local value
 		}
 		// Normalise by time
 		double timeSum = 0;
@@ -829,7 +864,9 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 			timeSum += time;
 		}
 		currentSum /= timeSum;
-		return currentSum;
+		
+		TERateResults teRateResults = new TERateResults(currentSum, localValues, conditioningEmbeddingsFromSpikes.size(), timeSum);
+		return teRateResults;
 	}
 
 
@@ -857,7 +894,6 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 			if (vectorOfConditionalSpikeTimes.size() > 0) {
 				conditionalIterator = vectorOfConditionalSpikeTimes.iterator();
 			}
-			int timeSeriesIndex = 0;
 			for (double[] destSpikeTimes : vectorOfDestinationSpikeTimes) {
 				double[] sourceSpikeTimes = sourceIterator.next();
 				double[][] conditionalSpikeTimes = null;
@@ -885,7 +921,7 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 			KdTree resampledKdTreeConditioningAtSamples = new KdTree(arrayedResampledConditioningEmbeddingsFromSamples);
 			resampledKdTreeConditioningAtSamples.setNormType(normType);
 
-			Vector<double[]> conditionallyPermutedJointEmbeddingsFromSpikes = new Vector(jointEmbeddingsFromSpikes);
+			Vector<double[]> conditionallyPermutedJointEmbeddingsFromSpikes = new Vector<double[]>(jointEmbeddingsFromSpikes);
 
 			Vector<Integer> usedIndices = new Vector<Integer>();
 			for (int i = 0; i < conditionallyPermutedJointEmbeddingsFromSpikes.size(); i++) {
@@ -920,13 +956,28 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 			}
 			KdTree conditionallyPermutedKdTreeJointFromSpikes = new KdTree(arrayedConditionallyPermutedJointEmbeddingsFromSpikes);
 			conditionallyPermutedKdTreeJointFromSpikes.setNormType(normType);
-			surrogateTEValues[permutationNumber] = computeAverageLocalOfObservations(conditionallyPermutedKdTreeJointFromSpikes,
-												 conditionallyPermutedJointEmbeddingsFromSpikes);
+			TERateResults teRateResults = computeAverageLocalOfObservations(conditionallyPermutedKdTreeJointFromSpikes,
+					 conditionallyPermutedJointEmbeddingsFromSpikes);
+			surrogateTEValues[permutationNumber] = teRateResults.meanRate;
 		}
 		return new EmpiricalMeasurementDistribution(surrogateTEValues, estimatedValue);
 	}
 	
-
+	// Data structure to return local values from this estimator
+	public class SpikingTELocalValues implements SpikingLocalInformationValues {
+		public double[] contributionsAtEachSpike; // Will have zeros for the spikes which are part of the first embedding
+		public double[][] jointEmbeddingsAtSpikes; // Rows here correspond to each target spike we computed a local TE for.
+					// First numDestPastIntervals columns are for the target embedded ISIs,
+					// next numCondPastIntervals columns are for any ISIs relevant to conditional variables,
+					// and last numSourcePastIntervals columns are for the ISIs relevant to the source spikes
+		
+		public SpikingTELocalValues(double[] contributionsAtEachSpike) {
+			this.contributionsAtEachSpike = contributionsAtEachSpike;
+			// int numTargetSpikesToSkip = MatrixUtils.max(destPastIntervals);
+			this.jointEmbeddingsAtSpikes = arrayedJointEmbeddingsFromSpikes;
+		}
+	}
+	
 	/*
 	 * (non-Javadoc)
 	 * 
@@ -935,8 +986,8 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
 	 */
 	@Override
 	public SpikingLocalInformationValues computeLocalOfPreviousObservations() throws Exception {
-		// TODO Auto-generated method stub
-		return null;
+		TERateResults teRateResults = computeAverageLocalOfObservations(kdTreeJointAtSpikes, jointEmbeddingsFromSpikes);
+		return new SpikingTELocalValues(teRateResults.localContributions);
 	}
 
 

java/source/infodynamics/utils/KdTree.java

@@ -113,6 +113,20 @@ public class KdTree extends NearestNeighbourSearcher {
 		this(new int[] {data[0].length}, new double[][][] {data});
 	}
 	
+	/**
+	 * Construct the k-d tree from a set of double[][] data.
+	 * 
+	 * @param data a double[][] 2D data set, first indexed
+	 *  by time, second index by variable number.
+	 * @param observationSetIndices array indicating for each sample which 
+	 *  observation set is came from (only used for dynamic correlation exclusion)
+	 * @param observationTimePoints array indicating for each sample which
+	 *  time index it had in the observation set it came from
+	 */
+	public KdTree(double[][] data, int[] observationSetIndices, int[] observationTimePoints) {
+		this(new int[] {data[0].length}, new double[][][] {data}, observationSetIndices, observationTimePoints);
+	}
+
 	/**
 	 * Construct the k-d tree from a <b>set</b> of double[][] data,
 	 * considered jointly.
@@ -126,7 +140,29 @@ public class KdTree extends NearestNeighbourSearcher {
 	 *   within this data set)
 	 */
 	public KdTree(int[] dimensions, double[][][] data) {
-				
+		this(dimensions, data, null, null);
+	}
+	
+	/**
+	 * Construct the k-d tree from a <b>set</b> of double[][] data,
+	 * considered jointly.
+	 * 
+	 * @param dimensions an array of dimensions for each
+	 *  of the 2D data sets.
+	 * @param data an array of double[][] 2D data sets
+	 *  for each data[i]
+	 *  (where i is the main variable number within data, then
+	 *   after that the first index is sample number, second is dimension
+	 *   within this data set)
+	 * @param observationSetIndices array indicating for each sample which 
+	 *  observation set is came from (only used for dynamic correlation exclusion).
+	 *  null means only a single observation set used
+	 * @param observationTimePoints array indicating for each sample which
+	 *  time index it had in the observation set it came from
+	 *  null means only a single observation set used
+	 */
+	public KdTree(int[] dimensions, double[][][] data, int[] observationSetIndices, int[] observationTimePoints) {
+			
 		this.originalDataSets = data;
 		int numObservations = data[0].length;
 				
@@ -183,6 +219,18 @@ public class KdTree extends NearestNeighbourSearcher {
 		rootNode = constructKdTree(0, 0, numObservations, masterSortedArrayIndices);
 		// And destroy the temporary storage of sorted array indices:
 		masterSortedArrayIndices = null;
+		
+		if (observationSetIndices == null) {
+			// observationSetIndices and observationTimePoints are
+			//  not supplied, so by default we assume only the one observation set:
+			observationSetIndices = new int[numObservations];
+			observationTimePoints = new int[numObservations];
+			for (int n = 0; n < numObservations; n++) {
+				observationTimePoints[n] = n;
+			}
+		}
+		this.observationSetIndices = observationSetIndices;
+		this.observationTimePoints = observationTimePoints;
 	}
 
 	/**
@@ -667,7 +715,7 @@ public class KdTree extends NearestNeighbourSearcher {
 	 * @param sampleIndex sample index in the data to find a nearest neighbour
 	 *  for
 	 * @param dynCorrExclTime size of dynamic correlation exclusion time window
-	 *  on either side of sampleIndex. 0 means exclude only sampleIndex itself.
+	 *  on either side of sampleIndex in the same observation set. 0 means exclude only sampleIndex itself.
 	 * @param node node to start searching from in the kd-tree. Cannot be null
 	 * @param level which level we're currently at in the tree
 	 * @param currentKBest a PriorityQueue of NeighbourNodeData objects
@@ -700,8 +748,8 @@ public class KdTree extends NearestNeighbourSearcher {
 		//  (will not throw an Exception if the PQ is empty)
 		NeighbourNodeData furthestCached = currentKBest.peek();
 		
-		if (((node.indexOfThisPoint - sampleIndex > dynCorrExclTime)
-			|| (node.indexOfThisPoint - sampleIndex < -dynCorrExclTime)) &&
+		if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex])
+			|| (Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
 			((currentKBest.size() < K) || (absDistOnThisDim < furthestCached.distance))) {
 			// Preliminary check says we need to compute the full distance
 			//  to use or at least to check if it should be 
@@ -1073,7 +1121,8 @@ public class KdTree extends NearestNeighbourSearcher {
 			absDistOnThisDim = distOnThisDim * distOnThisDim;
 		}
 		
-		if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+		if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex])
+				|| (Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
 			((absDistOnThisDim <  r) ||
 			 ( allowEqualToR && (absDistOnThisDim == r)))) {
 			// Preliminary check says we need to compute the full distance
@@ -1591,7 +1640,8 @@ public class KdTree extends NearestNeighbourSearcher {
 			absDistOnThisDim = distOnThisDim * distOnThisDim;
 		}
 		
-		if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+		if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) ||
+				(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
 			((absDistOnThisDim <  r) ||
 			 ( allowEqualToR && (absDistOnThisDim == r)))) {
 			// Preliminary check says we need to compute the full distance
@@ -1738,7 +1788,8 @@ public class KdTree extends NearestNeighbourSearcher {
 			absDistOnThisDim = distOnThisDim * distOnThisDim;
 		}
 		
-		if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+		if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) ||
+				(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
 			((absDistOnThisDim <  r) ||
 			 ( allowEqualToR && (absDistOnThisDim == r)))) {
 			// Preliminary check says we need to compute the full distance
@@ -1912,7 +1963,8 @@ public class KdTree extends NearestNeighbourSearcher {
 		}
 		
 		if (testResultsForGivenVariable[node.indexOfThisPoint] &&
-			(Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+			((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) ||
+					(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
 			((absDistOnThisDim <  r) ||
 			 ( allowEqualToR && (absDistOnThisDim == r)))) {
 			// Preliminary check says we need to compute the full distance
@@ -2574,7 +2626,8 @@ public class KdTree extends NearestNeighbourSearcher {
 			absDistOnThisDim = distOnThisDim * distOnThisDim;
 		}
 		
-		if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
+		if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) || 
+				(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
 			((absDistOnThisDim <  rs[variableNumber]) ||
 			 ( allowEqualToR && (absDistOnThisDim == rs[variableNumber])))) {
 			// Preliminary check says we need to compute the full distance

java/source/infodynamics/utils/NearestNeighbourSearcher.java

@@ -40,6 +40,19 @@ public abstract class NearestNeighbourSearcher {
 	 */
 	protected int normTypeToUse = EuclideanUtils.NORM_MAX_NORM;
 
+	/**
+	 * observationSetIndices is an array indicating for each sample which 
+	 *  observation set is came from (only used for dynamic correlation exclusion).
+	 *  null means only a single observation set used
+	 */
+	protected int[] observationSetIndices;
+	/**
+	 * observationTimePoints is an array indicating for each sample which
+	 *  time index it had in the observation set it came from
+	 *  null means only a single observation set used
+	 */
+	protected int[] observationTimePoints;
+
 	/**
 	 * Factory method to construct the searcher from a set of double[][] data.
 	 * This will return a {@link KdTree} or if the data is univaraite
@@ -51,17 +64,37 @@ public abstract class NearestNeighbourSearcher {
 	public static NearestNeighbourSearcher create(double[][] data)
 			throws Exception {
 		
+		return create(data, null, null);
+	}
+	
+	/**
+	 * Factory method to construct the searcher from a set of double[][] data.
+	 * This will return a {@link KdTree} or if the data is univaraite
+	 * (i.e. only one column) a {@link UnivariateNearestNeighbourSearcher}
+	 * 
+	 * @param data a double[][] 2D data set, first indexed
+	 *  by time, second index by variable number.
+	 * @param observationSetIndices array indicating for each sample which 
+	 *  observation set is came from (only used for dynamic correlation exclusion)
+	 * @param observationTimePoints array indicating for each sample which
+	 *  time index it had in the observation set it came from
+	 */
+	public static NearestNeighbourSearcher create(double[][] data,
+			int[] observationSetIndices, int[] observationTimePoints)
+			throws Exception {
+		
 		if ((data == null) || (data[0].length == 0)) {
 			// We have null data:
 			return null;
 		} else if (data[0].length == 1) {
 			// We have univariate data:
-			return new UnivariateNearestNeighbourSearcher(MatrixUtils.selectColumn(data, 0));
+			return new UnivariateNearestNeighbourSearcher(MatrixUtils.selectColumn(data, 0),
+					observationSetIndices, observationTimePoints);
 		} else {
-			return new KdTree(data);
+			return new KdTree(data, observationSetIndices, observationTimePoints);
 		}
 	}
-	
+
 	/**
 	 * Factory method to construct the searcher from a set of double[][][] data.
 	 * 
@@ -71,11 +104,29 @@ public abstract class NearestNeighbourSearcher {
 	public static NearestNeighbourSearcher create(int[] dimensions, double[][][] data)
 			throws Exception {
 		
+		return create(dimensions, data, null, null);
+	}
+	
+	/**
+	 * Factory method to construct the searcher from a set of double[][][] data.
+	 * 
+	 * @param data an array of double[][] 2D data sets, first indexed
+	 *  by time, second index by variable number.
+	 * @param observationSetIndices array indicating for each sample which 
+	 *  observation set is came from (only used for dynamic correlation exclusion)
+	 * @param observationTimePoints array indicating for each sample which
+	 *  time index it had in the observation set it came from
+	 */
+	public static NearestNeighbourSearcher create(int[] dimensions, double[][][] data,
+			int[] observationSetIndices, int[] observationTimePoints)
+			throws Exception {
+		
 		if ((dimensions.length == 1) && (dimensions[0] == 1)) {
 			// We have univariate data:
-			return new UnivariateNearestNeighbourSearcher(MatrixUtils.selectColumn(data[0], 0));
+			return new UnivariateNearestNeighbourSearcher(MatrixUtils.selectColumn(data[0], 0),
+					observationSetIndices, observationTimePoints);
 		} else {
-			return new KdTree(dimensions, data);
+			return new KdTree(dimensions, data, observationSetIndices, observationTimePoints);
 		}
 	}
 	

java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java

@@ -63,15 +63,25 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 	protected double[] sortedValues = null;
 	
 	public UnivariateNearestNeighbourSearcher(double[][] data) throws Exception {
+		this(data, null, null);
+	}
+	
+	public UnivariateNearestNeighbourSearcher(double[][] data,
+			int[] observationSetIndices, int[] observationTimePoints) throws Exception {
 		// Ideally we would not call the constructor until after the following check,
 		//  but the constructor must come first in Java.
-		this(MatrixUtils.selectColumn(data, 0));
+		this(MatrixUtils.selectColumn(data, 0), observationSetIndices, observationTimePoints);
 		if (data[0].length != 1) {
 			throw new Exception("Cannot define UnivariateNearestNeighbourSearcher for multivariate data");
 		}
 	}
-	
+
 	public UnivariateNearestNeighbourSearcher(double[] data) throws Exception {
+		this(data, null, null);
+	}
+	
+	public UnivariateNearestNeighbourSearcher(double[] data,
+			int[] observationSetIndices, int[] observationTimePoints) throws Exception {
 		this.originalDataSet = data;
 		numObservations = data.length;
 		if (numObservations <= 1) {
@@ -102,6 +112,18 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		for (int i = 0; i < numObservations; i++) {
 			sortedValues[i] = originalDataSet[sortedArrayIndices[i]];
 		}
+
+		if (observationSetIndices == null) {
+			// observationSetIndices and observationTimePoints are
+			//  not supplied, so by default we assume only the one observation set:
+			observationSetIndices = new int[numObservations];
+			observationTimePoints = new int[numObservations];
+			for (int n = 0; n < numObservations; n++) {
+				observationTimePoints[n] = n;
+			}
+		}
+		this.observationSetIndices = observationSetIndices;
+		this.observationTimePoints = observationTimePoints;
 	}
 
 	/**
@@ -252,7 +274,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 				lowerCandidate >= 0;
 				lowerCandidate--) {
 			indexOfLowerCandidate = sortedArrayIndices[lowerCandidate];
-			if (Math.abs(sampleIndex - indexOfLowerCandidate) > dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfLowerCandidate]) || 
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfLowerCandidate]) > dynCorrExclTime)) {
 				// This sample is outside the dynamic correlation exclusion window
 				break;
 			}
@@ -266,7 +289,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 				upperCandidate <= numObservations - 1;
 				upperCandidate++) {
 			indexOfUpperCandidate = sortedArrayIndices[upperCandidate];
-			if (Math.abs(sampleIndex - indexOfUpperCandidate) > dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfUpperCandidate]) ||
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfUpperCandidate]) > dynCorrExclTime)) {
 				// This sample is outside the dynamic correlation exclusion window
 				break;
 			}
@@ -301,7 +325,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 						upperCandidate <= numObservations - 1;
 						upperCandidate++) {
 					indexOfUpperCandidate = sortedArrayIndices[upperCandidate];
-					if (Math.abs(sampleIndex - indexOfUpperCandidate) > dynCorrExclTime) {
+					if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfUpperCandidate]) ||
+							(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfUpperCandidate]) > dynCorrExclTime)) {
 						// This sample is outside the dynamic correlation exclusion window
 						break;
 					}
@@ -317,7 +342,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 						lowerCandidate >= 0;
 						lowerCandidate--) {
 					indexOfLowerCandidate = sortedArrayIndices[lowerCandidate];
-					if (Math.abs(sampleIndex - indexOfLowerCandidate) > dynCorrExclTime) {
+					if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfLowerCandidate]) || 
+							(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfLowerCandidate]) > dynCorrExclTime)) {
 						// This sample is outside the dynamic correlation exclusion window
 						break;
 					}
@@ -474,7 +500,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		int indexInSortedArray = indicesInSortedArray[sampleIndex];
 		// Check the points with smaller data values first:
 		for (int i = indexInSortedArray - 1; i >= 0; i--) {
-			if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
 				// Can't count this point, but keep checking:
 				continue;
 			}
@@ -490,7 +517,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		}
 		// Next check the points with larger data values:
 		for (int i = indexInSortedArray + 1; i < numObservations; i++) {
-			if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
 				// Can't count this point, but keep checking:
 				continue;
 			}
@@ -747,7 +775,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		int indexInSortedArray = indicesInSortedArray[sampleIndex];
 		// Check the points with smaller data values first:
 		for (int i = indexInSortedArray - 1; i >= 0; i--) {
-			if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
 				// Can't count this point, but keep checking:
 				continue;
 			}
@@ -764,7 +793,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		}
 		// Next check the points with larger data values:
 		for (int i = indexInSortedArray + 1; i < numObservations; i++) {
-			if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
 				// Can't count this point, but keep checking:
 				continue;
 			}
@@ -793,7 +823,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		int indexInSortedArray = indicesInSortedArray[sampleIndex];
 		// Check the points with smaller data values first:
 		for (int i = indexInSortedArray - 1; i >= 0; i--) {
-			if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
 				// Can't count this point, but keep checking:
 				continue;
 			}
@@ -812,7 +843,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		}
 		// Next check the points with larger data values:
 		for (int i = indexInSortedArray + 1; i < numObservations; i++) {
-			if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
+			if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+					(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
 				// Can't count this point, but keep checking:
 				continue;
 			}
@@ -868,7 +900,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		// Check the points with smaller data values first:
 		for (int i = indexInSortedArray - 1; i >= 0; i--) {
 			if (!additionalCriteria[sortedArrayIndices[i]] || 
-					(Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime)) {
+					((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+							(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime))) {
 				// Can't count this point, but keep checking:
 				continue;
 			}
@@ -889,7 +922,8 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
 		// Next check the points with larger data values:
 		for (int i = indexInSortedArray + 1; i < numObservations; i++) {
 			if (!additionalCriteria[sortedArrayIndices[i]] || 
-					(Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime)) {
+					((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
+							(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime))) {
 				// Can't count this point, but keep checking:
 				continue;
 			}

java/unittests/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoMultiVariateTester.java

@@ -159,6 +159,40 @@ public class ConditionalMutualInfoMultiVariateTester extends
 				condMi, 0.0000000001);
 	}
 	
+	public void testNoConditional() throws Exception {
+		ArrayFileReader afr = new ArrayFileReader("demos/data/4ColsPairedOneStepNoisyDependence-1.txt");
+		double[][] data = afr.getDouble2DMatrix();
+		double[] source = MatrixUtils.selectColumn(data, 1);
+		double[] dest = MatrixUtils.selectColumn(data, 2);
+		
+		// Set up the value we expect from MI:
+		MutualInfoCalculatorMultiVariateGaussian miCalc =
+				new MutualInfoCalculatorMultiVariateGaussian();
+		miCalc.initialise(1, 1);
+		miCalc.setObservations(source, dest);	
+		double mi = miCalc.computeAverageLocalOfObservations();
+		
+		// Now compute via CMI calculator with null passed:
+		ConditionalMutualInfoCalculatorMultiVariateGaussian condMiCalc =
+				new ConditionalMutualInfoCalculatorMultiVariateGaussian();
+		condMiCalc.initialise(1, 1, 0);
+		condMiCalc.setObservations(source, dest, (double[][]) null);	
+		double condMi = condMiCalc.computeAverageLocalOfObservations();
+		assertEquals(mi, condMi, 0.000001);
+
+		// Now compute via CMI calculator with dummy column passed:
+		condMiCalc.initialise(1, 1, 0);
+		condMiCalc.setObservations(source, dest, MatrixUtils.selectColumn(data, 3));	
+		condMi = condMiCalc.computeAverageLocalOfObservations();
+		assertEquals(mi, condMi, 0.000001);
+
+		// Now compute via CMI calculator with empty column passed (all as 2D):
+		condMiCalc.initialise(1, 1, 0);
+		condMiCalc.setObservations(MatrixUtils.selectColumns(data, 1, 1), MatrixUtils.selectColumns(data, 2, 1), new double[1000][0]);	
+		condMi = condMiCalc.computeAverageLocalOfObservations();
+		assertEquals(mi, condMi, 0.000001);
+	}
+	
 	public void testBiasCorrectionDoesNotChangeAnalyticPValue() throws Exception {
 		ConditionalMutualInfoCalculatorMultiVariateGaussian cmiCalc =
 				new ConditionalMutualInfoCalculatorMultiVariateGaussian();
@@ -189,7 +223,7 @@ public class ConditionalMutualInfoMultiVariateTester extends
 		ChiSquareMeasurementDistribution distroBiasCorrected = cmiCalc.computeSignificance();
 		assertEquals(avBiasCorrected, distroBiasCorrected.actualValue, 0.0000001);
 		// And now check that the pValues are unchanged whether we bias correct or not:
-		assertEquals(distroNotBiasCorrected.pValue, distroBiasCorrected.pValue);
+		assertEquals(distroNotBiasCorrected.pValue, distroBiasCorrected.pValue, 0.0000001);
 	}
 
 	protected int timeStepsDepCheck = 100;
@@ -395,7 +429,7 @@ public class ConditionalMutualInfoMultiVariateTester extends
 		//  - both dimensions are copied
 		double[][] sourceData = rg.generateNormalData(timeStepsDepCheck, dimensions,
 				0, 1);
-		double[][] condData = rg.generateNormalData(timeStepsDepCheck, dimensions,
+		double[][] condData = rg.generateNormalData(timeStepsDepCheck, conditionalDims,
 				0, 1);
 		double[][] destData = MatrixUtils.arrayCopy(condData);
 
@@ -447,4 +481,5 @@ public class ConditionalMutualInfoMultiVariateTester extends
 		condMiCalc.setObservations(MatrixUtils.selectColumns(destData, 0, 1), sourceData, condData);
 		assertTrue(Double.isInfinite(condMiCalc.computeAverageLocalOfObservations()));
 	}
+	
 }

java/unittests/infodynamics/measures/continuous/kraskov/ActiveInfoStorageTester.java

@@ -5,6 +5,7 @@ import java.util.Arrays;
 import infodynamics.measures.continuous.ActiveInfoStorageCalculator;
 import infodynamics.utils.ArrayFileReader;
 import infodynamics.utils.MatrixUtils;
+import infodynamics.utils.RandomGenerator;
 import junit.framework.TestCase;
 
 public class ActiveInfoStorageTester extends TestCase {
@@ -116,4 +117,74 @@ public class ActiveInfoStorageTester extends TestCase {
 				ais.getProperty(ActiveInfoStorageCalculator.TAU_PROP_NAME), differentNNResult);
 	}
 	
-}
+	/**
+	 * Test that observationSetIndices and observationStartTimePoints 
+	 * for the underlying MI estimator are written properly
+	 * 
+	 * @throws Exception 
+	 */
+	@SuppressWarnings("static-access")
+	public void testObservationSetIndices() throws Exception {
+		
+		int timeSteps = 1000;
+		
+		ActiveInfoStorageCalculatorKraskov aisCalc = 
+				new ActiveInfoStorageCalculatorKraskov();
+		RandomGenerator rg = new RandomGenerator();
+		double[] data = rg.generateNormalData(timeSteps, 0, 1);
+		aisCalc.initialise();
+		
+		// First check that for a simple single observation set everything works:
+		aisCalc.setObservations(data);
+
+		int[] observationSetIds = aisCalc.getObservationSetIndices();
+		int[] timeSeriesIndices = aisCalc.getObservationTimePoints();
+		assert(observationSetIds.length == timeSteps - 1);
+		for (int t = 0; t < timeSteps - 1; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], 1 + t); // Should be offset by 1 for k history
+		}
+		
+		// Now add the same one twice:
+		aisCalc.initialise();
+		aisCalc.startAddObservations();
+		aisCalc.addObservations(data);
+		aisCalc.addObservations(data);
+		aisCalc.finaliseAddObservations();
+		observationSetIds = aisCalc.getObservationSetIndices();
+		timeSeriesIndices = aisCalc.getObservationTimePoints();
+		assert(observationSetIds.length == 2*timeSteps - 2);
+		for (int t = 0; t < timeSteps - 1; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], 1 + t);
+		}
+		for (int t = 0; t < timeSteps - 1; t++) {
+			assertEquals(observationSetIds[timeSteps - 1 + t], 1);
+			assertEquals(timeSeriesIndices[timeSteps - 1 + t], 1 + t);
+		}
+		
+		// Now add NUM_SEGMENTS segments:
+		int NUM_SEGMENTS = 10;
+		aisCalc.setProperty(aisCalc.K_PROP_NAME, "3");
+		aisCalc.setProperty(aisCalc.TAU_PROP_NAME, "3");
+		aisCalc.initialise();
+		aisCalc.startAddObservations();
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			aisCalc.addObservations(rg.generateNormalData(timeSteps, 0, 1));
+		}
+		aisCalc.finaliseAddObservations();
+		observationSetIds = aisCalc.getObservationSetIndices();
+		timeSeriesIndices = aisCalc.getObservationTimePoints();
+		int k = Integer.valueOf(aisCalc.getProperty(aisCalc.K_PROP_NAME)); // In case we change to autoembedding above
+		int tau = Integer.valueOf(aisCalc.getProperty(aisCalc.TAU_PROP_NAME));
+		System.out.printf("Embedding dimension %d and delay %d\n", k, tau);
+		assert(observationSetIds.length == (timeSteps - (k - 1)*tau - 1) * NUM_SEGMENTS);
+		int t = 0;
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			for (int i = 0; i < timeSteps - (k-1)*tau - 1; i++) {
+				assertEquals(observationSetIds[t], r);
+				assertEquals(timeSeriesIndices[t], (k-1)*tau + 1 + i);
+				t++;
+			}
+		}
+	}}

java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java

@@ -21,6 +21,7 @@ package infodynamics.measures.continuous.kraskov;
 import infodynamics.utils.ArrayFileReader;
 import infodynamics.utils.MathsUtils;
 import infodynamics.utils.MatrixUtils;
+import infodynamics.utils.RandomGenerator;
 
 public class ConditionalMutualInfoMultiVariateTester
 	extends infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateAbstractTester {
@@ -97,10 +98,11 @@ public class ConditionalMutualInfoMultiVariateTester
 	 * @param var2 dest multivariate data set
 	 * @param kNNs array of Kraskov k nearest neighbours parameter to check
 	 * @param expectedResults array of expected results for each k
+	 * @return errors of the computed values against expectedResults
 	 */
-	protected void checkTEForGivenData(double[][] var1, double[][] var2,
+	protected double[] checkTEForGivenData(double[][] var1, double[][] var2,
 			int[] kNNs, double[] expectedResults) throws Exception {
-		checkTEForGivenData(var1, var2, 1, 1, kNNs, expectedResults);
+		return checkTEForGivenData(var1, var2, 1, 1, kNNs, expectedResults);
 	}
 	
 	/**
@@ -114,12 +116,38 @@ public class ConditionalMutualInfoMultiVariateTester
 	 * @param historyL history length l of source
 	 * @param kNNs array of Kraskov k nearest neighbours parameter to check
 	 * @param expectedResults array of expected results for each k
+	 * @return errors of the computed values against expectedResults
 	 */
-	protected void checkTEForGivenData(double[][] var1, double[][] var2,
+	protected double[] checkTEForGivenData(double[][] var1, double[][] var2,
 			int historyK, int historyL, int[] kNNs, double[] expectedResults) throws Exception {
+		return checkTEForGivenData(var1, var2, historyK, historyL, kNNs, expectedResults,
+				0, "NONE", 0.000001);
+	}
+	
+	/**
+	 * Utility function to run Kraskov conditional MI algorithm 1
+	 *  as transfer entropy for data with known results
+	 *  from TRENTOOL.
+	 * 
+	 * @param var1 source multivariate data set
+	 * @param var2 dest multivariate data set
+	 * @param historyK history length k of destination
+	 * @param historyL history length l of source
+	 * @param kNNs array of Kraskov k nearest neighbours parameter to check
+	 * @param expectedResults array of expected results for each k
+	 * @param noiseLevel noise to add to the data - set to 0 if we 
+	 *   need to exactly reproduce calculations (most cases in unit tests for consistency)
+	 * @param noiseSeed seed for the random noise generator (either "NONE" or Long string)
+	 * @param tolerance tolerance to accept the calculation
+	 * @return errors of the computed values against expectedResults
+	 */
+	protected double[] checkTEForGivenData(double[][] var1, double[][] var2,
+			int historyK, int historyL, int[] kNNs, double[] expectedResults,
+			double noiseLevel, String noiseSeed, double tolerance) throws Exception {
 				
 		ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = getNewCalc(1);
-		
+		double[] errors = new double[expectedResults.length]; 
+
 		// Which is the first time index for the dest next state?
 		//  It depends on the values of k and l for embedding the past state
 		//  of destination and source.
@@ -147,7 +175,9 @@ public class ConditionalMutualInfoMultiVariateTester
 			condMiCalc.setProperty(
 					ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
 					NUM_THREADS_TO_USE);
-			condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency of results for unit test
+			condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE,
+					Double.toString(noiseLevel));
+			condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NOISE_SEED, noiseSeed);
 			condMiCalc.initialise(var1[0].length * historyL,
 					var2[0].length, var2[0].length * historyK);
 			// Construct the joint vectors of the source states
@@ -192,8 +222,10 @@ public class ConditionalMutualInfoMultiVariateTester
 			System.out.printf("k=%d: Average MI %.8f (expected %.8f)\n",
 					k, condMi, expectedResults[kIndex]);
 			// 6 decimal places is Matlab accuracy
-			assertEquals(expectedResults[kIndex], condMi, 0.000001);			
+			assertEquals(expectedResults[kIndex], condMi, tolerance);
+			errors[kIndex] = condMi - expectedResults[kIndex];
 		}
+		return errors;
 	}
 
 	/**
@@ -699,7 +731,7 @@ public class ConditionalMutualInfoMultiVariateTester
 			//		EuclideanUtils.NORM_MAX_NORM_STRING);
 			condMiCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests
 			condMiCalc.initialise(var1[0].length, var2[0].length, 0);
-			condMiCalc.setObservations(var1, var2, null);
+			condMiCalc.setObservations(var1, var2, (double[][]) null);
 			// condMiCalc.setDebug(true);
 			double mi = condMiCalc.computeAverageLocalOfObservations();
 			// condMiCalc.setDebug(false);
@@ -736,4 +768,143 @@ public class ConditionalMutualInfoMultiVariateTester
 			assertEquals(expectedMi, miFromEmptyVectors, 0.0000001);
 		}
 	}
+	
+	/**
+	 * Test that observationSetIndices and observationStartTimePoints are written properly
+	 * 
+	 * @throws Exception 
+	 */
+	public void testObservationSetIndices() throws Exception {
+		
+		int dimensions = 1;
+		int timeSteps = 100;
+		
+		ConditionalMutualInfoCalculatorMultiVariateKraskov condMiCalc = getNewCalc(1);
+		condMiCalc.initialise(dimensions, dimensions, dimensions);
+		
+		// generate some random data
+		RandomGenerator rg = new RandomGenerator();
+		double[][] sourceData = rg.generateNormalData(timeSteps, dimensions,
+				0, 1);
+		double[][] destData = rg.generateNormalData(timeSteps, dimensions,
+				0, 1);
+		double[][] condData = rg.generateNormalData(timeSteps, dimensions,
+				0, 1);
+		
+		// First check that for a simple single observation set everything works:
+		condMiCalc.setObservations(sourceData, destData, condData);
+
+		int[] observationSetIds = condMiCalc.getObservationSetIndices();
+		int[] timeSeriesIndices = condMiCalc.getObservationTimePoints();
+		assert(observationSetIds.length == timeSteps);
+		for (int t = 0; t < timeSteps; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], t);
+		}
+		
+		// Now add the same one twice:
+		condMiCalc.initialise(dimensions, dimensions, dimensions);
+		condMiCalc.startAddObservations();
+		condMiCalc.addObservations(sourceData, destData, condData);
+		condMiCalc.addObservations(sourceData, destData, condData);
+		condMiCalc.finaliseAddObservations();
+		observationSetIds = condMiCalc.getObservationSetIndices();
+		timeSeriesIndices = condMiCalc.getObservationTimePoints();
+		assert(observationSetIds.length == 2*timeSteps);
+		for (int t = 0; t < timeSteps; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], t);
+		}
+		for (int t = 0; t < timeSteps; t++) {
+			assertEquals(observationSetIds[timeSteps + t], 1);
+			assertEquals(timeSeriesIndices[timeSteps + t], t);
+		}
+		
+		// Now add NUM_SEGMENTS randomly chosen segments:
+		int NUM_SEGMENTS = 10;
+		int maxLength = 10;
+		int[] startPoints = rg.generateRandomInts(NUM_SEGMENTS, timeSteps - maxLength);
+		int[] lengthsMinus1 = rg.generateRandomInts(NUM_SEGMENTS, maxLength - 1); // ensures we don't add segments of length 0
+		condMiCalc.initialise(dimensions, dimensions, dimensions);
+		condMiCalc.startAddObservations();
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			condMiCalc.addObservations(sourceData, destData, condData, startPoints[r], lengthsMinus1[r]+1);
+		}
+		condMiCalc.finaliseAddObservations();
+		observationSetIds = condMiCalc.getObservationSetIndices();
+		timeSeriesIndices = condMiCalc.getObservationTimePoints();
+		assert(observationSetIds.length == MatrixUtils.sum(lengthsMinus1) + NUM_SEGMENTS);
+		int t = 0;
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			for (int i = 0; i < lengthsMinus1[r]+1; i++) {
+				assertEquals(observationSetIds[t], r);
+				assertEquals(timeSeriesIndices[t], startPoints[r] + i);
+				t++;
+			}
+		}
+	}	
+
+	/**
+	 * Extends testUnivariateTEforCoupledVariablesFromFile to test
+	 * using seed for random number generator
+	 * 
+	 * @throws Exception if file not found 
+	 * 
+	 */
+	public void testWithSeed() throws Exception {
+		
+		// Test set 1:
+		
+		ArrayFileReader afr = new ArrayFileReader("demos/data/2coupledRandomCols-1.txt");
+		double[][] data = afr.getDouble2DMatrix();
+		
+		// Use various Kraskov k nearest neighbours parameter
+		int[] kNNs = {4};
+		// Expected values from TRENTOOL:
+		double[] expectedFromTRENTOOL = {0.3058006};
+		
+		System.out.println("Kraskov Cond MI as TE comparison 1 - univariate coupled data 1");
+		double[] noNoiseError = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+				MatrixUtils.selectColumns(data, new int[] {1}),
+				kNNs, expectedFromTRENTOOL);
+		double noNoiseResult = expectedFromTRENTOOL[0] + noNoiseError[0];
+		
+		// And now in the reverse direction:
+		double[] expectedFromTRENTOOLRev = new double[] {-0.0029744};
+		
+		System.out.println("  reverse direction:");
+		double[] noNoiseErrorRev = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
+				MatrixUtils.selectColumns(data, new int[] {0}),
+				kNNs, expectedFromTRENTOOLRev);
+		double noNoiseResultRev = expectedFromTRENTOOLRev[0] + noNoiseErrorRev[0];
+
+		// Check that changing the random number generator still returns close to those with no noise 
+		//  results, but with larger tolerance
+		System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1");
+		double[] withSeed1Error = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+				MatrixUtils.selectColumns(data, new int[] {1}),
+				1, 1, kNNs, expectedFromTRENTOOL,
+				1e-8, "1", 0.01);
+		double[] withSeed1Results = new double[] {expectedFromTRENTOOL[0] + withSeed1Error[0]};
+		// Now check that the results are exact when we repeat with the same seed:
+		System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1 repeat want exact");
+		checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+				MatrixUtils.selectColumns(data, new int[] {1}),
+				1, 1, kNNs, withSeed1Results,
+				1e-8, "1", 1e-10);
+
+		// And in reverse direction:
+		System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1 reverse");
+		double[] withSeed1ErrorRev = checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
+				MatrixUtils.selectColumns(data, new int[] {0}),
+				1, 1, kNNs, expectedFromTRENTOOLRev,
+				1e-8, "1", 0.01);
+		double[] withSeed1ResultsRev = new double[] {expectedFromTRENTOOLRev[0] + withSeed1ErrorRev[0]};
+		// Now check that the results are exact when we repeat with the same seed:
+		System.out.println("\n Kraskov Cond MI as TE comparison 1 - univariate coupled data 1 - seed 1 reverse repeat want exact");
+		checkTEForGivenData(MatrixUtils.selectColumns(data, new int[] {1}),
+				MatrixUtils.selectColumns(data, new int[] {0}),
+				1, 1, kNNs, withSeed1ResultsRev,
+				1e-8, "1", 1e-10);
+}
 }

java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java

@@ -21,6 +21,7 @@ package infodynamics.measures.continuous.kraskov;
 import infodynamics.utils.ArrayFileReader;
 import infodynamics.utils.MathsUtils;
 import infodynamics.utils.MatrixUtils;
+import infodynamics.utils.RandomGenerator;
 
 public class MutualInfoMultiVariateTester
 	extends infodynamics.measures.continuous.MutualInfoMultiVariateAbstractTester {
@@ -95,6 +96,23 @@ public class MutualInfoMultiVariateTester
 		checkComputeSignificanceDoesntAlterAverage(2);
 	}
 	
+	/**
+	 * Utility function to run Kraskov MI for data with known results.
+	 * Sets to use no noise in the calculation and a tolerance of 0.0000001
+	 * 
+	 * @param var1
+	 * @param var2
+	 * @param kNNs array of Kraskov k nearest neighbours parameter to check
+	 * @param expectedResults array of expected results for each k
+	 * @return errors of the computed values against expectedResults
+	 */
+	protected double[] checkMIForGivenData(double[][] var1, double[][] var2,
+			int[] kNNs, double[] expectedResults) throws Exception {
+		// Dropping required accuracy by one order of magnitude, due
+		//  to faster but slightly less accurate digamma estimator change
+		return checkMIForGivenData(var1, var2, kNNs, expectedResults, 0, "NONE", 0.0000001);
+	}
+	
 	/**
 	 * Utility function to run Kraskov MI for data with known results
 	 * 
@@ -102,13 +120,20 @@ public class MutualInfoMultiVariateTester
 	 * @param var2
 	 * @param kNNs array of Kraskov k nearest neighbours parameter to check
 	 * @param expectedResults array of expected results for each k
+	 * @param noiseLevel noise to add to the data - set to 0 if we 
+	 *   need to exactly reproduce calculations (most cases in unit tests for consistency)
+	 * @param noiseSeed seed for the random noise generator (either "NONE" or Long string)
+	 * @param tolerance tolerance to accept the calculation
+	 * @return errors of the computed values against expectedResults
 	 */
-	protected void checkMIForGivenData(double[][] var1, double[][] var2,
-			int[] kNNs, double[] expectedResults) throws Exception {
-				
+	protected double[] checkMIForGivenData(double[][] var1, double[][] var2,
+			int[] kNNs, double[] expectedResults,
+			double noiseLevel, String noiseSeed, double tolerance) throws Exception {
+			
 		// The Kraskov MILCA toolkit MIhigherdim executable 
 		//  uses algorithm 2 by default (this is what it means by rectangular):
 		MutualInfoCalculatorMultiVariateKraskov miCalc = getNewCalc(2);
+		double[] errors = new double[expectedResults.length]; 
 		
 		for (int kIndex = 0; kIndex < kNNs.length; kIndex++) {
 			int k = kNNs[kIndex];
@@ -121,7 +146,9 @@ public class MutualInfoMultiVariateTester
 			// No longer need to set this property as it's set by default:
 			//miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE,
 			//		EuclideanUtils.NORM_MAX_NORM_STRING);
-			miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests
+			miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE,
+					Double.toString(noiseLevel));
+			miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_NOISE_SEED, noiseSeed);
 			miCalc.initialise(var1[0].length, var2[0].length);
 			miCalc.setObservations(var1, var2);
 			miCalc.setDebug(true);
@@ -130,10 +157,10 @@ public class MutualInfoMultiVariateTester
 			
 			System.out.printf("k=%d: Average MI %.8f (expected %.8f)\n",
 					k, mi, expectedResults[kIndex]);
-			// Dropping required accuracy by one order of magnitude, due
-			//  to faster but slightly less accurate digamma estimator change
-			assertEquals(expectedResults[kIndex], mi, 0.0000001);			
+			assertEquals(expectedResults[kIndex], mi, tolerance);
+			errors[kIndex] = mi - expectedResults[kIndex];
 		}
+		return errors;
 	}
 	
 	/**
@@ -652,5 +679,195 @@ public class MutualInfoMultiVariateTester
 		assertEquals(expected_H_X_given_Y, conditionalEnt, 0.02);
 
 	}
+	
+	/**
+	 * Test that observationSetIndices and observationStartTimePoints are written properly
+	 * 
+	 * @throws Exception 
+	 */
+	public void testObservationSetIndices() throws Exception {
+		
+		int dimensions = 1;
+		int timeSteps = 100;
+		
+		MutualInfoCalculatorMultiVariateKraskov miCalc = getNewCalc(1);
+		miCalc.initialise(dimensions, dimensions);
+		
+		// generate some random data
+		RandomGenerator rg = new RandomGenerator();
+		double[][] sourceData = rg.generateNormalData(timeSteps, dimensions,
+				0, 1);
+		double[][] destData = rg.generateNormalData(timeSteps, dimensions,
+				0, 1);
+		
+		// First check that for a simple single observation set everything works:
+		miCalc.setObservations(sourceData, destData);
+
+		int[] observationSetIds = miCalc.getObservationSetIndices();
+		int[] timeSeriesIndices = miCalc.getObservationTimePoints();
+		assert(observationSetIds.length == timeSteps);
+		for (int t = 0; t < timeSteps; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], t);
+		}
+		
+		// Now add the same one twice:
+		miCalc.initialise(dimensions, dimensions);
+		miCalc.startAddObservations();
+		miCalc.addObservations(sourceData, destData);
+		miCalc.addObservations(sourceData, destData);
+		miCalc.finaliseAddObservations();
+		observationSetIds = miCalc.getObservationSetIndices();
+		timeSeriesIndices = miCalc.getObservationTimePoints();
+		assert(observationSetIds.length == 2*timeSteps);
+		for (int t = 0; t < timeSteps; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], t);
+		}
+		for (int t = 0; t < timeSteps; t++) {
+			assertEquals(observationSetIds[timeSteps + t], 1);
+			assertEquals(timeSeriesIndices[timeSteps + t], t);
+		}
+		
+		// Now add NUM_SEGMENTS randomly chosen segments:
+		int NUM_SEGMENTS = 10;
+		int maxLength = 10;
+		int[] startPoints = rg.generateRandomInts(NUM_SEGMENTS, timeSteps - maxLength);
+		int[] lengthsMinus1 = rg.generateRandomInts(NUM_SEGMENTS, maxLength - 1); // ensures we don't add segments of length 0
+		miCalc.initialise(dimensions, dimensions);
+		miCalc.startAddObservations();
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			miCalc.addObservations(sourceData, destData, startPoints[r], lengthsMinus1[r]+1);
+		}
+		miCalc.finaliseAddObservations();
+		observationSetIds = miCalc.getObservationSetIndices();
+		timeSeriesIndices = miCalc.getObservationTimePoints();
+		assert(observationSetIds.length == MatrixUtils.sum(lengthsMinus1) + NUM_SEGMENTS);
+		int t = 0;
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			for (int i = 0; i < lengthsMinus1[r]+1; i++) {
+				assertEquals(observationSetIds[t], r);
+				assertEquals(timeSeriesIndices[t], startPoints[r] + i);
+				t++;
+			}
+		}
+	}
+
+	/**
+	 * Test the MI values returned, and that we can basically turn dyn corr excl
+	 *  off by adding each data point in separately
+	 * @throws Exception 
+	 */
+	public void testDynCorrExclAndSetIndices2() throws Exception {
+		int nForEach = 100;
+		double[] x = new double[nForEach * 4];
+		double[] y = new double[nForEach * 4];
+		
+		int k = 4; // Make sure this is even
+		int exclWindow = 2; // Make sure this is capped by k/2. Leave this at 2 - these of the code assumes this
+		
+		int segment = 0;
+		for (int rx = 0; rx < 2; rx++) {
+			for (int ry = 0; ry < 2; ry++) {
+				for (int t = 0; t < nForEach; t++) {
+					// x and y coordinate will be == t % nForEach when rx==0 or ry==0,
+					//  else we will push them up by nForEach*2..
+					// So we will basically have diagonal lines in four quadrants.
+					x[segment*nForEach + t] = t + ((rx == 1) ? nForEach*2 : 0);
+					y[segment*nForEach + t] = t + ((ry == 1) ? nForEach*2 : 0);
+				}
+				segment++;
+			}
+		}
+		// We're going to set kNNs to be even, so that the nearest neighbours must be 
+		//  the kNN closest points on the same line segement. n_x and n_y will be 2*(kNN) + 1 respectively with algorithm 2,
+		//  since they will include the kNN in the same line segment (>= than the kNNth), plus the 
+		//  corresponding point and it's kNNs in the segment with the same x or y range.
+		// So we have:
+		double expectedMIWithoutDynCorrExcl = MathsUtils.digamma(k) - 1.0 / (double)k - 2*MathsUtils.digamma(2*k+1) + MathsUtils.digamma(nForEach*4);
+		// But if we turn on dynamic correlation exclusion, say for 2 points, then the kNNs will now step outside the previous k
+		//  and the kernel widths for n_x and n_y will be wider by approx 2x the exclusion window size.
+		// For points not on the edges, n_x and n_y will increase by 2x the exclusion window size up to k/2 since they
+		//  will include this many more points in the corresponding segments.
+		// For points near the edges, n_x and n_y will increase only by 1x the exclusion window on one side (where it doesn't abut the edge), and the amount of the 
+		//  exclusion window that includes points on the other side.
+		// The maths we're using here just assumes that the window is size 2 to make coding it fast (we're not going to bother changing it)
+		double expectedMIWithDynCorrExcl =
+				(double) ((nForEach-4) * 4) / (double) (nForEach * 4) * (MathsUtils.digamma(k) - 1.0 / (double)k - 2*MathsUtils.digamma(2*k+1+2*exclWindow) + MathsUtils.digamma(nForEach*4)) +
+				(double) ((2) * 4) / (double) (nForEach * 4) * (MathsUtils.digamma(k) - 1.0 / (double)k - 2*MathsUtils.digamma(2*k+1+exclWindow + exclWindow-1) + MathsUtils.digamma(nForEach*4)) +
+				(double) ((2) * 4) / (double) (nForEach * 4) * (MathsUtils.digamma(k) - 1.0 / (double)k - 2*MathsUtils.digamma(2*k+1+exclWindow + exclWindow-2) + MathsUtils.digamma(nForEach*4));
+		
+		MutualInfoCalculatorMultiVariateKraskov miCalc = getNewCalc(2);
+		miCalc.setProperty(miCalc.PROP_K, Integer.toString(k));
+		miCalc.setProperty(miCalc.PROP_ADD_NOISE, Integer.toString(0)); // Need to noise so that our analytic results hold
+
+		// First check that for a simple single observation set everything works:
+		miCalc.initialise(1, 1);
+		miCalc.setObservations(x, y);
+		double miWithoutDynCorrExcl = miCalc.computeAverageLocalOfObservations();
+		assertEquals(expectedMIWithoutDynCorrExcl, miWithoutDynCorrExcl, 0.00001);
+
+		// Now turn dynamic correlation exclusion on:
+		miCalc.setProperty(miCalc.PROP_DYN_CORR_EXCL_TIME, Integer.toString(exclWindow));
+		// miCalc.setProperty(miCalc.PROP_NUM_THREADS, Integer.toString(1)); // To simplify debugging
+		miCalc.initialise(1, 1);
+		miCalc.setObservations(x, y);
+		double miWithDynCorrExcl = miCalc.computeAverageLocalOfObservations();
+		assertEquals(expectedMIWithDynCorrExcl, miWithDynCorrExcl, 0.00001);
+
+		// This time, we leave dynamic correlation exclusion turned on, but add each sample as
+		// a separate data set -- this will serve to effectively turn dynamic correlation exclusion off again!
+		// Gives a good test that exclusion only considers points within the same data set.
+		miCalc.initialise(1, 1);
+		miCalc.startAddObservations();
+		for (int t = 0; t < x.length; t++) {
+			miCalc.addObservations(new double[] {x[t]}, new double[] {y[t]});
+		}
+		miCalc.finaliseAddObservations();
+		double miWithDynCorrExclAndSeperateSets = miCalc.computeAverageLocalOfObservations();
+		assertEquals(expectedMIWithoutDynCorrExcl, miWithDynCorrExclAndSeperateSets, 0.00001);
+	}
+
+	/**
+	 * Extends the tests from testUnivariateMIforRandomVariablesFromFile to use
+	 *  a random seed for noise and check that results are repeatable.
+	 * 
+	 * @throws Exception if file not found 
+	 * 
+	 */
+	public void testUnivariateMIWithSeed() throws Exception {
+		
+		// Test set 1:
+		
+		ArrayFileReader afr = new ArrayFileReader("demos/data/2randomCols-1.txt");
+		double[][] data = afr.getDouble2DMatrix();
+		
+		// Use various Kraskov k nearest neighbours parameter
+		int[] kNNs = {1, 2, 3, 4, 5, 6, 10, 15};
+		// Expected values from Kraskov's MILCA toolkit:
+		double[] expectedFromMILCA = {-0.05294175, -0.03944338, -0.02190217,
+				0.00120807, -0.00924771, -0.00316402, -0.00778205, -0.00565778};
+		
+		System.out.println("Kraskov comparison 1 - univariate random data 1 - no seed");
+		double[] noNoiseErrors = checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+				MatrixUtils.selectColumns(data, new int[] {1}),
+				kNNs, expectedFromMILCA);
+		double[] noNoiseResults = MatrixUtils.add(expectedFromMILCA, noNoiseErrors);
+		
+		// Check that changing the random number generator still returns close to those with no noise 
+		//  results, but with larger tolerance
+		System.out.println("\n Kraskov comparison 1 - univariate random data 1 - seed 1");
+		double[] withSeed1Errors = checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+				MatrixUtils.selectColumns(data, new int[] {1}),
+				kNNs, noNoiseResults,
+				1e-8, "1", 0.01);
+		double[] withSeed1Results = MatrixUtils.add(noNoiseResults, withSeed1Errors);
+		// Now check that the results are exact when we repeat with the same seed:
+		System.out.println("\n Kraskov comparison 1 - univariate random data 1 - seed 1 repeat");
+		checkMIForGivenData(MatrixUtils.selectColumns(data, new int[] {0}),
+				MatrixUtils.selectColumns(data, new int[] {1}),
+				kNNs, withSeed1Results,
+				1e-8, "1", 1e-10);
+	}
 }
 

java/unittests/infodynamics/measures/continuous/kraskov/TransferEntropyTester.java

@@ -737,4 +737,76 @@ public class TransferEntropyTester
 			assertEquals(correctL, optimisedL);
 	  }
 
+	/**
+	 * Test that observationSetIndices and observationStartTimePoints 
+	 * for the underlying CMI estimator are written properly
+	 * 
+	 * @throws Exception 
+	 */
+	@SuppressWarnings("static-access")
+	public void testObservationSetIndices() throws Exception {
+		
+		int timeSteps = 1000;
+		
+		TransferEntropyCalculatorKraskov teCalc = 
+				new TransferEntropyCalculatorKraskov();
+		RandomGenerator rg = new RandomGenerator();
+		double[] source = rg.generateNormalData(timeSteps, 0, 1);
+		double[] target = rg.generateNormalData(timeSteps, 0, 1);
+		teCalc.initialise();
+		
+		// First check that for a simple single observation set everything works:
+		teCalc.setObservations(source, target);
+
+		int[] observationSetIds = teCalc.getObservationSetIndices();
+		int[] timeSeriesIndices = teCalc.getObservationTimePoints();
+		assert(observationSetIds.length == timeSteps - 1);
+		for (int t = 0; t < timeSteps - 1; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], 1 + t); // Should be offset by 1 for k history
+		}
+		
+		// Now add the same one twice:
+		teCalc.initialise();
+		teCalc.startAddObservations();
+		teCalc.addObservations(source, target);
+		teCalc.addObservations(source, target);
+		teCalc.finaliseAddObservations();
+		observationSetIds = teCalc.getObservationSetIndices();
+		timeSeriesIndices = teCalc.getObservationTimePoints();
+		assert(observationSetIds.length == 2*timeSteps - 2);
+		for (int t = 0; t < timeSteps - 1; t++) {
+			assertEquals(observationSetIds[t], 0);
+			assertEquals(timeSeriesIndices[t], 1 + t);
+		}
+		for (int t = 0; t < timeSteps - 1; t++) {
+			assertEquals(observationSetIds[timeSteps - 1 + t], 1);
+			assertEquals(timeSeriesIndices[timeSteps - 1 + t], 1 + t);
+		}
+		
+		// Now add NUM_SEGMENTS segments:
+		int NUM_SEGMENTS = 10;
+		teCalc.setProperty(teCalc.K_PROP_NAME, "3");
+		teCalc.setProperty(teCalc.K_TAU_PROP_NAME, "3");
+		teCalc.initialise();
+		teCalc.startAddObservations();
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			teCalc.addObservations(rg.generateNormalData(timeSteps, 0, 1), rg.generateNormalData(timeSteps, 0, 1));
+		}
+		teCalc.finaliseAddObservations();
+		observationSetIds = teCalc.getObservationSetIndices();
+		timeSeriesIndices = teCalc.getObservationTimePoints();
+		int k = Integer.valueOf(teCalc.getProperty(teCalc.K_PROP_NAME)); // In case we change to autoembedding above
+		int tau = Integer.valueOf(teCalc.getProperty(teCalc.K_TAU_PROP_NAME));
+		System.out.printf("Embedding dimension %d and delay %d\n", k, tau);
+		assert(observationSetIds.length == (timeSteps - (k - 1)*tau - 1) * NUM_SEGMENTS);
+		int t = 0;
+		for (int r = 0; r < NUM_SEGMENTS; r++) {
+			for (int i = 0; i < timeSteps - (k-1)*tau - 1; i++) {
+				assertEquals(observationSetIds[t], r);
+				assertEquals(timeSeriesIndices[t], (k-1)*tau + 1 + i);
+				t++;
+			}
+		}
+	}
 }

java/unittests/infodynamics/utils/KdTreeTest.java

@@ -422,6 +422,72 @@ public class KdTreeTest extends TestCase {
 		}
 	}
 
+	public void testFindKNearestNeighboursWithExclusionWindowAndMultiDataSets() throws Exception {
+		int dimension = 4;
+		int numSamplesPerSet = 400;
+		int numSets = 5;
+		int exclusionWindow = 50;
+		
+		for (int K = 1; K < 5; K++) {
+			double[][] data = rg.generateNormalData(numSamplesPerSet*numSets, dimension, 0, 1);
+			int[] obsSetIds = new int[numSamplesPerSet*numSets];
+			int[] timeIndicesInSets = new int[numSamplesPerSet*numSets];
+			int ti = 0;
+			for (int s = 0; s < numSets; s++) {
+				for (int s2 = 0; s2 < numSamplesPerSet; s2++) {
+					obsSetIds[ti] = s;
+					timeIndicesInSets[ti] = s2;
+					ti++;
+				}
+			}
+			
+			long startTime = Calendar.getInstance().getTimeInMillis();
+			KdTree kdTree = new KdTree(data, obsSetIds, timeIndicesInSets);
+			long endTimeTree = Calendar.getInstance().getTimeInMillis();
+			System.out.printf("Tree of %d points for %d NNs constructed in: %.3f sec\n",
+					data.length, K, ((double) (endTimeTree - startTime)/1000.0));
+			
+			EuclideanUtils normCalculator = new EuclideanUtils(EuclideanUtils.NORM_MAX_NORM);
+			startTime = Calendar.getInstance().getTimeInMillis();
+			for (int t = 0; t < data.length; t++) {
+				PriorityQueue<NeighbourNodeData> nnPQ =
+						kdTree.findKNearestNeighbours(K, t, exclusionWindow);
+				assertTrue(nnPQ.size() == K);
+				// Now find the K nearest neighbours with a naive all-pairs comparison
+				double[][] distancesAndIndices = new double[data.length][2];
+				for (int t2 = 0; t2 < data.length; t2++) {
+					boolean inDifferentSets = ((t / numSamplesPerSet) != (t2 / numSamplesPerSet));
+					if (inDifferentSets || (Math.abs(t2 - t) > exclusionWindow)) {
+						distancesAndIndices[t2][0] = normCalculator.norm(data[t], data[t2]);
+					} else {
+						distancesAndIndices[t2][0] = Double.POSITIVE_INFINITY;
+					}
+					distancesAndIndices[t2][1] = t2;
+				}
+				int[] timeStepsOfKthMins =
+						MatrixUtils.kMinIndices(distancesAndIndices, 0, K);
+				for (int i = 0; i < K; i++) {
+					// Check that the ith nearest neighbour matches for each method.
+					// Note that these two method provide a different sorting order
+					NeighbourNodeData nnData = nnPQ.poll();
+					if (timeStepsOfKthMins[K - 1 - i] != nnData.sampleIndex) {
+						// We have an error:
+						System.out.printf("Erroneous match between indices %d (expected) " +
+						 " and %d\n", timeStepsOfKthMins[K - 1 - i], nnData.sampleIndex);
+					}
+					assertEquals(timeStepsOfKthMins[K - 1 - i], nnData.sampleIndex);
+					// And check that none of the nearest neighbours were within the window
+					//  and from the same data set
+					boolean inDifferentSets = ((t / numSamplesPerSet) != (nnData.sampleIndex / numSamplesPerSet));
+					assertTrue(inDifferentSets || (Math.abs(nnData.sampleIndex - t) > exclusionWindow));
+				}
+			}		
+			long endTimeValidate = Calendar.getInstance().getTimeInMillis();
+			System.out.printf("All %d nearest neighbours found in: %.3f sec\n",
+					K, ((double) (endTimeValidate - startTime)/1000.0));
+		}
+	}
+
 	public void testFindKNearestNeighboursForSeparateArrays() throws Exception {
 		int variables = 3;
 		int dimensionsPerVariable = 3;

java/unittests/infodynamics/utils/UnivariateNearestNeighbourTest.java

@@ -625,4 +625,66 @@ public class UnivariateNearestNeighbourTest extends TestCase {
 					K, ((double) (endTimeValidate - nnEndTime)/1000.0));
 		}
 	}
+
+	public void testFindKNearestNeighboursWithExclusionWindowAndMultiDataSets() throws Exception {
+		int numSamplesPerSet = 400;
+		int numSets = 5;
+		int exclusionWindow = 50;
+		
+		for (int K = 1; K < 5; K++) {
+			double[] data = rg.generateNormalData(numSamplesPerSet*numSets, 0, 1);
+			int[] obsSetIds = new int[numSamplesPerSet*numSets];
+			int[] timeIndicesInSets = new int[numSamplesPerSet*numSets];
+			int ti = 0;
+			for (int s = 0; s < numSets; s++) {
+				for (int s2 = 0; s2 < numSamplesPerSet; s2++) {
+					obsSetIds[ti] = s;
+					timeIndicesInSets[ti] = s2;
+					ti++;
+				}
+			}
+			
+			long startTime = Calendar.getInstance().getTimeInMillis();
+			UnivariateNearestNeighbourSearcher searcher = new UnivariateNearestNeighbourSearcher(data, obsSetIds, timeIndicesInSets);
+			long endTimeTree = Calendar.getInstance().getTimeInMillis();
+			System.out.printf("Searcher of %d points for %d NNs constructed in: %.3f sec\n",
+					data.length, K, ((double) (endTimeTree - startTime)/1000.0));
+			
+			startTime = Calendar.getInstance().getTimeInMillis();
+			for (int t = 0; t < data.length; t++) {
+				PriorityQueue<NeighbourNodeData> nnPQ =
+						searcher.findKNearestNeighbours(K, t, exclusionWindow);
+				assertTrue(nnPQ.size() == K);
+				// Now find the K nearest neighbours with a naive all-pairs comparison
+				double[][] distancesAndIndices = new double[data.length][2];
+				for (int t2 = 0; t2 < data.length; t2++) {
+					boolean inDifferentSets = ((t / numSamplesPerSet) != (t2 / numSamplesPerSet));
+					// If we weren't catering for different sample sets, it would run like this (so if you run this it will lead to an error):
+					// if ((Math.abs(t2 - t) > exclusionWindow)) {
+					if (inDifferentSets || (Math.abs(t2 - t) > exclusionWindow)) {
+						distancesAndIndices[t2][0] = Math.abs(data[t] - data[t2]);
+					} else {
+						distancesAndIndices[t2][0] = Double.POSITIVE_INFINITY;
+					}
+					distancesAndIndices[t2][1] = t2;
+				}
+				int[] timeStepsOfKthMins =
+						MatrixUtils.kMinIndices(distancesAndIndices, 0, K);
+				for (int i = 0; i < K; i++) {
+					// Check that the ith nearest neighbour matches for each method.
+					// Note that these two method provide a different sorting order
+					NeighbourNodeData nnData = nnPQ.poll();
+					if (timeStepsOfKthMins[K - 1 - i] != nnData.sampleIndex) {
+						// We have an error:
+						System.out.printf("Erroneous match between indices %d (expected) " +
+						 " and %d\n", timeStepsOfKthMins[K - 1 - i], nnData.sampleIndex);
+					}
+					assertEquals(timeStepsOfKthMins[K - 1 - i], nnData.sampleIndex);
+				}
+			}		
+			long endTimeValidate = Calendar.getInstance().getTimeInMillis();
+			System.out.printf("All %d nearest neighbours found in: %.3f sec\n",
+					K, ((double) (endTimeValidate - startTime)/1000.0));
+		}
+	}
 }

readme-template.txt

@@ -136,6 +136,12 @@ Notices for this software are found in the notices/JAMA directory.
  Release notes
 ===============
 
+v1.6.1 22/8/2023
+-------------
+(after 909 commits recorded by github, repository as at https://github.com/jlizier/jidt/tree/90baf68ee7332e15030447b44d262a0fc54773f6 save for this file update)
+Minor updates to supporting use in Python, including virtual environments;
+Minor tweaks to fish schooling examples (mostly comments)
+
 v1.6 5/9/2022
 -------------
 (after 889 commits recorded by github, repository as at https://github.com/jlizier/jidt/tree/d750a737bea2a8b1f33b7cd0ad167ec999d907ef save for this file update)