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

build.xml

@@ -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/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;
 

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
+

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));
+		}
+	}
 }