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

Compare commits

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

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

32 changed files with 640 additions and 2803 deletions
Show Stats Expand all files Collapse all files

14
build.xml
View File

@ -43,17 +43,17 @@
<!-- Compile the java toolkit -->
<target name="compile" depends="init" description="compile the source">
<!-- Compile to Java 8 to provide compatibility for users with older JREs.
<!-- Compile to Java 7 to provide compatibility for users with older JREs.
Caveat: The flags here only check the language compatibility, but
may still use newer libraries which may cause issues for users with JDK 8.
Indeed, one gets the warning: "bootstrap class path not set in conjunction with -source 1.8"
may still use newer libraries which may cause issues for users with JDK 7.
Indeed, one gets the warning: "bootstrap class path not set in conjunction with -source 1.7"
To fix this, one would use the bootstrap classpath to point our JDK to an rt.jar
for Java 8.
At this stage, I'm sure I'm not using new library calls from Java 9+, so we can
for Java 7.
At this stage, I'm sure I'm not using new library calls from Java 8+, so we can
ignore the warning, and I don't want to bother installing newer Java just to compile
like this. I'll endeavour not to use JDK 9+ libraries so as not to cause
like this. I'll endeavour not to use JDK 8 libraries so as not to cause
any issues here ... -->
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.8" source="1.8" encoding="UTF8">
<javac srcdir="${src}" destdir="${bin}" includeAntRuntime="false" target="1.7" source="1.7" encoding="UTF8">
<classpath refid="apache-classpath"/>
</javac>

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

@ -43,8 +43,6 @@ 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;

87
demos/python/create_env.bat
View File

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

84
java/source/infodynamics/measures/continuous/ActiveInfoStorageCalculatorViaMutualInfo.java
View File

@ -96,12 +96,6 @@ 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
@ -248,7 +242,6 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
this.tau = tau;
vectorOfObservationTimeSeries = null;
vectorOfValidityOfObservations = null;
vectorOfOffsetsInTimeSeries = null;
}
/**
@ -265,7 +258,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 max bias-corrected AIS criteria (both searching up to {@link #PROP_K_SEARCH_MAX} and
* the maz 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
@ -373,18 +366,6 @@ 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()
*/
@ -392,7 +373,6 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
public void startAddObservations() {
vectorOfObservationTimeSeries = new Vector<double[]>();
vectorOfValidityOfObservations = new Vector<boolean[]>();
vectorOfOffsetsInTimeSeries = new Vector<Integer>();
}
/* (non-Javadoc)
@ -400,15 +380,9 @@ 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)
@ -419,7 +393,6 @@ 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
}
/**
@ -430,13 +403,10 @@ 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,
int observationSetIndex, int offsetInOriginalTimeSeries) throws Exception {
int k_in_use, int tau_in_use, double[] observations) 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
@ -449,8 +419,7 @@ 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.addObservationsTrackObservationIDs(currentDestPastVectors, currentDestNextVectors,
observationSetIndex, offsetInOriginalTimeSeries + (k_in_use-1)*tau_in_use + 1);
miCalc_in_use.addObservations(currentDestPastVectors, currentDestNextVectors);
}
/**
@ -467,12 +436,10 @@ 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,
int observationSetIndex) throws Exception {
int k_in_use, int tau_in_use, double[] observations, boolean[] valid) throws Exception {
// compute the start and end times using our determined embedding parameters:
Vector<int[]> startAndEndTimePairs = computeStartAndEndTimePairs(k_in_use, tau_in_use, valid);
@ -481,8 +448,7 @@ 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),
observationSetIndex, startTime);
MatrixUtils.select(observations, startTime, endTime - startTime + 1));
}
}
@ -492,7 +458,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
@Override
public void addObservations(double[] observations, int startTime,
int numTimeSteps) throws Exception {
addObservationsParsed(MatrixUtils.select(observations, startTime, numTimeSteps), startTime);
addObservations(MatrixUtils.select(observations, startTime, numTimeSteps));
}
/**
@ -650,7 +616,6 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
vectorOfObservationTimeSeries = null; // No longer required
vectorOfValidityOfObservations = null;
vectorOfOffsetsInTimeSeries = null;
}
/**
@ -674,25 +639,34 @@ 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,
setNum, offsetsInTimeSeriesIterator.next());
tau_in_use, observations);
} else {
addObservationsWithGivenParams(miCalc_in_use, k_in_use,
tau_in_use, observations, validity, setNum);
tau_in_use, observations, validity);
}
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.
@ -881,20 +855,4 @@ 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();
}
}

278
java/source/infodynamics/measures/continuous/ConditionalMutualInfoCalculatorMultiVariate.java
View File

@ -92,16 +92,6 @@ 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,
@ -139,8 +129,17 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
double[][] cond) throws Exception;
/**
* As per {@link #setObservations(double[][], double[][], double[][])}
* but where dimensions of all variables have been set to 1
* 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>
*
* @param var1 univariate observations for variable 1
* @param var2 univariate observations for variable 2
@ -152,90 +151,6 @@ 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.
@ -333,8 +248,18 @@ public interface ConditionalMutualInfoCalculatorMultiVariate
double[][] cond) throws Exception;
/**
* <p>As per {@link #addObservations(double[][], double[][], double[][])}
* but can only be used where dimensions of all
* <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
* variables have been set to 1.</p>
*
* @param var1 univariate observations for variable 1
@ -347,96 +272,6 @@ 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.
@ -465,28 +300,6 @@ 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.
*
@ -623,38 +436,6 @@ 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
@ -673,17 +454,4 @@ 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();
}

337
java/source/infodynamics/measures/continuous/ConditionalMutualInfoMultiVariateCommon.java
View File

@ -22,7 +22,6 @@ 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;
@ -83,16 +82,6 @@ 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.
@ -130,21 +119,6 @@ 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
@ -166,14 +140,6 @@ 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
@ -222,11 +188,6 @@ 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;
@ -258,8 +219,6 @@ 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>
@ -270,8 +229,6 @@ 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)) {
@ -283,20 +240,6 @@ 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);
}
}
@ -306,12 +249,6 @@ 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;
@ -319,7 +256,8 @@ 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();
@ -328,8 +266,7 @@ 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 python jpype.
* Resolved if one follows the AutoAnalyser generated code, but left for back compatibility.
* some problems calling overloaded versions of setObservations from jpype.
*
* @param var1
* @param var2
@ -341,55 +278,8 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
}
@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();
addedMoreThanOneObservationSet = false;
}
@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();
@ -398,8 +288,7 @@ 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 python jpype.
* Resolved if one follows the AutoAnalyser generated code, but left for back compatibility.
* some problems calling overloaded versions of setObservations from jpype.
*
* @param var1
* @param var2
@ -410,47 +299,16 @@ 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");
@ -475,8 +333,6 @@ 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;
}
@ -508,7 +364,6 @@ 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),
@ -517,110 +372,6 @@ 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.
@ -638,12 +389,6 @@ 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");
@ -657,14 +402,15 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
condToAdd = new double[numTimeSteps][];
System.arraycopy(cond, startTime, condToAdd, 0, numTimeSteps);
}
addObservationsTrackObservationIDs(var1ToAdd, var2ToAdd, condToAdd, observationSetIndexToUse, startTime);
addObservations(var1ToAdd, var2ToAdd, condToAdd);
}
public void addObservations(double[][] var1, double[][] var2,
@Override
public void setObservations(double[][] var1, double[][] var2,
double[][] cond,
boolean[] var1Valid, boolean[] var2Valid,
boolean[] condValid) throws Exception {
Vector<int[]> startAndEndTimePairs =
computeStartAndEndTimePairs(var1Valid, var2Valid, condValid);
@ -673,22 +419,23 @@ 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, observationSetIndex);
addObservations(var1, var2, cond, startTime, endTime - startTime + 1);
}
observationSetIndex++;
finaliseAddObservations();
}
public void addObservations(double[][] var1, double[][] var2,
@Override
public void setObservations(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);
addObservations(var1, var2, cond, allVar1Valid, allVar2Valid, allCondValid);
setObservations(var1, var2, cond, allVar1Valid, allVar2Valid, allCondValid);
}
/**
* Signal that the observations are now all added, PDFs can now be constructed.
*
@ -711,14 +458,10 @@ 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();
@ -735,12 +478,6 @@ 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;
}
@ -765,9 +502,6 @@ 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] +=
@ -934,34 +668,6 @@ 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) {
@ -1057,13 +763,4 @@ public abstract class ConditionalMutualInfoMultiVariateCommon implements
return addedMoreThanOneObservationSet;
}
@Override
public int[] getObservationSetIndices() {
return observationSetIndices;
}
@Override
public int[] getObservationTimePoints() {
return observationTimePoints;
}
}

7
java/source/infodynamics/measures/continuous/EntropyCalculator.java
View File

@ -48,11 +48,8 @@ Theory' (John Wiley & Sons, New York, 1991).</li>
*/
public interface EntropyCalculator extends InfoMeasureCalculatorContinuous {
// 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.
// TODO Add addObservations() methods for entropy calculator
/**
* Sets the samples from which to compute the PDF for the entropy.
* Should only be called once, the last call contains the

74
java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariate.java
View File

@ -61,26 +61,6 @@ 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.
@ -92,20 +72,6 @@ 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>
@ -128,46 +94,6 @@ 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

342
java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariateCommon.java
View File

@ -1,342 +0,0 @@
/*
* 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;
}
}

43
java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariate.java
View File

@ -91,35 +91,6 @@ 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
@ -159,18 +130,4 @@ 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();
}

172
java/source/infodynamics/measures/continuous/MutualInfoMultiVariateCommon.java
View File

@ -23,7 +23,6 @@ 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;
@ -87,16 +86,6 @@ 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.
@ -138,22 +127,6 @@ 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
*/
@ -167,19 +140,9 @@ public abstract class MutualInfoMultiVariateCommon implements
*/
protected boolean addNoise = false;
/**
* 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)
* 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;
/* (non-Javadoc)
* @see infodynamics.measures.continuous.ChannelCalculatorCommon#initialise()
@ -200,11 +163,6 @@ public abstract class MutualInfoMultiVariateCommon implements
sourceStdsBeforeNorm = null;
destMeansBeforeNorm = null;
destStdsBeforeNorm = null;
observationSetIndices = null;
observationTimePoints = null;
observationSetIndex = 0;
vectorOfObservationSetIndices = null;
vectorOfObservationStartTimePoints = null;
addedMoreThanOneObservationSet = false;
}
@ -232,8 +190,6 @@ 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>
@ -262,13 +218,6 @@ 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;
@ -289,12 +238,6 @@ 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;
@ -353,54 +296,14 @@ 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");
@ -423,8 +326,6 @@ public abstract class MutualInfoMultiVariateCommon implements
}
vectorOfSourceObservations.add(source);
vectorOfDestinationObservations.add(destination);
vectorOfObservationSetIndices.add(observationSetIndexToUse);
vectorOfObservationStartTimePoints.add(startTimeIndex);
}
/**
@ -490,14 +391,8 @@ 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");
@ -508,12 +403,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);
addObservationsTrackObservationIDs(sourceToAdd, destToAdd, observationSetIndexToUse, startTime);
vectorOfDestinationObservations.add(destToAdd);
}
@Override
public void addObservations(double[] source, double[] destination,
int startTime, int numTimeSteps) throws Exception {
@ -528,14 +423,8 @@ public abstract class MutualInfoMultiVariateCommon implements
startTime, numTimeSteps);
}
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,
public void setObservations(double[][] source, double[][] destination,
boolean[] sourceValid, boolean[] destValid) throws Exception {
Vector<int[]> startAndEndTimePairs = computeStartAndEndTimePairs(sourceValid, destValid);
@ -545,18 +434,32 @@ public abstract class MutualInfoMultiVariateCommon implements
for (int[] timePair : startAndEndTimePairs) {
int startTime = timePair[0];
int endTime = timePair[1];
addObservations(source, destination, startTime, endTime - startTime + 1, observationSetIndex);
addObservations(source, destination, startTime, endTime - startTime + 1);
}
observationSetIndex++;
finaliseAddObservations();
}
public void addObservations(double[][] source, double[][] destination,
@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,
boolean[][] sourceValid, boolean[][] destValid) throws Exception {
boolean[] allSourceValid = MatrixUtils.andRows(sourceValid);
boolean[] allDestValid = MatrixUtils.andRows(destValid);
addObservations(source, destination, allSourceValid, allDestValid);
setObservations(source, destination, allSourceValid, allDestValid);
}
/**
@ -579,15 +482,11 @@ 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
@ -598,14 +497,7 @@ public abstract class MutualInfoMultiVariateCommon implements
MatrixUtils.arrayCopy(destination, timeDiff, 0,
destObservations, startObservation, 0,
destination.length - timeDiff, dimensionsDest);
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;
startObservation += destination.length - timeDiff;
}
if (vectorOfSourceObservations.size() > 1) {
addedMoreThanOneObservationSet = true;
@ -626,9 +518,6 @@ 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] +=
@ -823,7 +712,6 @@ public abstract class MutualInfoMultiVariateCommon implements
return miSurrogateCalculator.computeAverageLocalOfObservations();
}
@Override
public double[] computeLocalUsingPreviousObservations(
double[] newSourceObservations, double[] newDestObservations)
throws Exception {
@ -946,14 +834,4 @@ public abstract class MutualInfoMultiVariateCommon implements
}
return startAndEndTimePairs;
}
@Override
public int[] getObservationSetIndices() {
return observationSetIndices;
}
@Override
public int[] getObservationTimePoints() {
return observationTimePoints;
}
}

79
java/source/infodynamics/measures/continuous/TransferEntropyCalculatorViaCondMutualInfo.java
View File

@ -148,12 +148,6 @@ 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
@ -358,7 +352,6 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
vectorOfDestinationTimeSeries = null;
vectorOfValidityOfSource = null;
vectorOfValidityOfDestination = null;
vectorOfOffsetsInTimeSeries = null;
separateNumObservations = new int[] {};
}
@ -508,38 +501,22 @@ 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);
}
/**
@ -559,7 +536,6 @@ 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
}
/**
@ -578,8 +554,6 @@ 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
*/
@ -587,8 +561,7 @@ 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,
int observationSetIndex, int offsetInOriginalTimeSeries) throws Exception {
double[] source, double[] destination) 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));
@ -614,8 +587,7 @@ 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.addObservationsTrackObservationIDs(currentSourcePastVectors, currentDestNextVectors, currentDestPastVectors,
observationSetIndex, offsetInOriginalTimeSeries + startTimeForFirstDestEmbedding_in_use + 1);
condMiCalc_in_use.addObservations(currentSourcePastVectors, currentDestNextVectors, currentDestPastVectors);
return destination.length - startTimeForFirstDestEmbedding_in_use - 1;
}
@ -643,7 +615,6 @@ 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
*/
@ -652,8 +623,7 @@ 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,
int observationSetIndex) throws Exception {
boolean[] sourceValid, boolean[] destValid) throws Exception {
// Compute the start and end time pairs using our embedding parameters:
Vector<int[]> startAndEndTimePairs =
@ -668,8 +638,7 @@ 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),
observationSetIndex, startTime);
MatrixUtils.select(destination, startTime, endTime - startTime + 1));
}
return totalObservationsAdded;
}
@ -685,8 +654,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");
}
addObservationsParsed(MatrixUtils.select(source, startTime, numTimeSteps),
MatrixUtils.select(destination, startTime, numTimeSteps), startTime);
addObservations(MatrixUtils.select(source, startTime, numTimeSteps),
MatrixUtils.select(destination, startTime, numTimeSteps));
}
/**
@ -872,7 +841,6 @@ public class TransferEntropyCalculatorViaCondMutualInfo implements
vectorOfDestinationTimeSeries = null; // No longer required
vectorOfValidityOfSource = null;
vectorOfValidityOfDestination = null;
vectorOfOffsetsInTimeSeries = null;
}
/**
@ -932,7 +900,6 @@ 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) {
@ -944,12 +911,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, setNum, offsetsInTimeSeriesIterator.next());
l_tau_in_use, delay_in_use, source, destination);
} 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, setNum);
sourceValidity, destValidity);
}
separateNumObservationsArray[setNum++] = observationsAddedThisTime;
}
@ -960,6 +927,20 @@ 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
@ -1150,22 +1131,6 @@ 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();

136
java/source/infodynamics/measures/continuous/gaussian/EntropyCalculatorMultiVariateGaussian.java
View File

@ -19,7 +19,6 @@
package infodynamics.measures.continuous.gaussian;
import infodynamics.measures.continuous.EntropyCalculatorMultiVariate;
import infodynamics.measures.continuous.EntropyCalculatorMultiVariateCommon;
import infodynamics.utils.MatrixUtils;
/**
@ -51,7 +50,6 @@ Theory' (John Wiley & Sons, New York, 1991).</li>
* <a href="http://lizier.me/joseph/">www</a>)
*/
public class EntropyCalculatorMultiVariateGaussian
extends EntropyCalculatorMultiVariateCommon
implements EntropyCalculatorMultiVariate, Cloneable {
/**
@ -65,34 +63,84 @@ 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() {
super();
// Nothing to do
}
public void initialise(int dimensions) {
super.initialise(dimensions);
means = null;
L = null;
detCovariance = 0;
@Override
public void initialise() throws Exception {
initialise(dimensions);
}
public void initialise(int 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 finaliseAddObservations() throws Exception {
super.finaliseAddObservations();
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");
}
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
observations = originalObservations;
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));
}
/**
@ -161,45 +209,49 @@ 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;
}
/**
* <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 setDebug(boolean debug) {
this.debug = debug;
}
@Override
public void setProperty(String propertyName, String propertyValue)
throws Exception {
// 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);
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;
}
/**
@ -264,13 +316,9 @@ public class EntropyCalculatorMultiVariateGaussian
if (observations == null) {
throw new Exception("Cannot compute local values since no observations were supplied");
}
double[] localValues = computeLocalUsingPreviousObservations(observations);
lastAverage = MatrixUtils.mean(localValues);
isComputed = true;
return localValues;
return computeLocalUsingPreviousObservations(observations);
}
@Override
public int getNumObservations() throws Exception {
if (observations == null) {
throw new Exception("Cannot return number of observations because either " +

111
java/source/infodynamics/measures/continuous/kernel/EntropyCalculatorMultiVariateKernel.java
View File

@ -19,7 +19,7 @@
package infodynamics.measures.continuous.kernel;
import infodynamics.measures.continuous.EntropyCalculatorMultiVariate;
import infodynamics.measures.continuous.EntropyCalculatorMultiVariateCommon;
import infodynamics.utils.MatrixUtils;
/**
* <p>Computes the differential entropy of a given set of observations
@ -44,13 +44,42 @@ import infodynamics.measures.continuous.EntropyCalculatorMultiVariateCommon;
* @author Joseph Lizier (<a href="joseph.lizier at gmail.com">email</a>,
* <a href="http://lizier.me/joseph/">www</a>)
*/
public class EntropyCalculatorMultiVariateKernel
extends EntropyCalculatorMultiVariateCommon
implements EntropyCalculatorMultiVariate
{
public class EntropyCalculatorMultiVariateKernel 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
*/
@ -72,10 +101,15 @@ public class EntropyCalculatorMultiVariateKernel
* 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) {
@ -94,22 +128,34 @@ public class EntropyCalculatorMultiVariateKernel
* 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 finaliseAddObservations() throws Exception {
super.finaliseAddObservations();
public void setObservations(double observations[][]) {
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]);
@ -119,9 +165,8 @@ public class EntropyCalculatorMultiVariateKernel
System.out.println(b + ": " + prob + " -> " + (-cont/Math.log(2.0)) + " -> sum: " + (entropy/Math.log(2.0)));
}
}
lastAverage = entropy / (double) totalObservations / Math.log(2.0);
isComputed = true;
return lastAverage;
lastEntropy = entropy / (double) totalObservations / Math.log(2.0);
return lastEntropy;
}
@Override
@ -159,18 +204,22 @@ public class EntropyCalculatorMultiVariateKernel
}
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) {
lastAverage = entropy;
isComputed = true;
lastEntropy = entropy;
}
return localEntropy;
}
@Override
public void setDebug(boolean debug) {
super.setDebug(debug);
this.debug = 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
@ -183,6 +232,8 @@ public class EntropyCalculatorMultiVariateKernel
* 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>
*
@ -203,15 +254,14 @@ public class EntropyCalculatorMultiVariateKernel
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 +
@ -224,10 +274,19 @@ public class EntropyCalculatorMultiVariateKernel
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 {
// try the superclass:
return super.getProperty(propertyName);
// No property was set, and no superclass to call:
return null;
}
}
@Override
public int getNumObservations() throws Exception {
return totalObservations;
}
}

263
java/source/infodynamics/measures/continuous/kozachenko/EntropyCalculatorMultiVariateKozachenko.java
View File

@ -18,10 +18,13 @@
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
@ -56,11 +59,59 @@ import infodynamics.utils.MathsUtils;
* @author Joseph Lizier (<a href="joseph.lizier at gmail.com">email</a>,
* <a href="http://lizier.me/joseph/">www</a>)
*/
public class EntropyCalculatorMultiVariateKozachenko
extends EntropyCalculatorMultiVariateCommon
public class EntropyCalculatorMultiVariateKozachenko
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
*/
@ -70,9 +121,171 @@ public class EntropyCalculatorMultiVariateKozachenko
* Construct an instance
*/
public EntropyCalculatorMultiVariateKozachenko() {
super();
noiseLevel = (double) 1e-8; // Default to align with KSG estimators
addNoise = true;
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);
}
/**
@ -85,12 +298,12 @@ public class EntropyCalculatorMultiVariateKozachenko
}
double sdTermHere = sdTerm(totalObservations, dimensions);
double emConstHere = eulerMascheroniTerm(totalObservations);
double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(observations);
double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(rawData);
double entropy = 0.0;
if (debug) {
System.out.println("t,\tminDist,\tlogMinDist,\tsum");
}
for (int t = 0; t < observations.length; t++) {
for (int t = 0; t < rawData.length; t++) {
entropy += Math.log(2.0 * minDistance[t]);
if (debug) {
System.out.println(t + ",\t" +
@ -119,18 +332,21 @@ 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(observations);
double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(rawData);
double entropy = 0.0;
double[] localEntropy = new double[observations.length];
double[] localEntropy = new double[rawData.length];
if (debug) {
System.out.println("t,\tminDist,\tlogMinDist,\tlocal,\tsum");
}
for (int t = 0; t < observations.length; t++) {
for (int t = 0; t < rawData.length; t++) {
localEntropy[t] = Math.log(2.0 * minDistance[t]) * (double) dimensions;
// using natural units
// localEntropy[t] /= Math.log(2);
@ -146,7 +362,7 @@ public class EntropyCalculatorMultiVariateKozachenko
}
entropy /= (double) totalObservations;
lastAverage = entropy;
isComputed = true;
lastLocalEntropy = localEntropy;
return localEntropy;
}
@ -159,13 +375,13 @@ public class EntropyCalculatorMultiVariateKozachenko
double emConstHere = eulerMascheroniTerm(totalObservations);
double constantToAddIn = sdTermHere + emConstHere;
double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(observations);
double[] minDistance = EuclideanUtils.computeMinEuclideanDistances(rawData);
double entropy = 0.0;
double[] localEntropy = new double[observations.length];
double[] localEntropy = new double[rawData.length];
if (debug) {
System.out.println("t,\tminDist,\tlogMinDist,\tlocal,\tsum");
}
for (int t = 0; t < observations.length; t++) {
for (int t = 0; t < rawData.length; t++) {
localEntropy[t] = Math.log(2.0 * minDistance[t]) * (double) dimensions;
// using natural units
// localEntropy[t] /= Math.log(2);
@ -179,6 +395,9 @@ public class EntropyCalculatorMultiVariateKozachenko
entropy);
}
}
entropy /= (double) totalObservations;
lastAverage = entropy;
lastLocalEntropy = localEntropy;
return localEntropy;
}
@ -254,4 +473,18 @@ 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;
}
}

45
java/source/infodynamics/measures/continuous/kraskov/ActiveInfoStorageCalculatorKraskov.java
View File

@ -23,7 +23,6 @@ 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})
@ -231,48 +230,4 @@ 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);
}
}

159
java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java
View File

@ -19,7 +19,6 @@
package infodynamics.measures.continuous.kraskov;
import java.util.Calendar;
import java.util.PriorityQueue;
import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariate;
import infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateCommon;
@ -28,7 +27,6 @@ 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;
@ -315,6 +313,13 @@ 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 +
@ -720,41 +725,35 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
if (kdTreeJoint == null) {
kdTreeJoint = new KdTree(
new int[] {dimensionsVar1, dimensionsVar2, dimensionsCond},
new double[][][] {var1Observations, var2Observations, condObservations},
observationSetIndices, observationTimePoints);
new double[][][] {var1Observations, var2Observations, condObservations});
kdTreeJoint.setNormType(normType);
}
if (dimensionsVar1 > 1) {
if (kdTreeVar1Conditional == null) {
kdTreeVar1Conditional = new KdTree(
new int[] {dimensionsVar1, dimensionsCond},
new double[][][] {var1Observations, condObservations},
observationSetIndices, observationTimePoints);
new double[][][] {var1Observations, condObservations});
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,
observationSetIndices, observationTimePoints);
uniNNSearcherVar1 = new UnivariateNearestNeighbourSearcher(var1Observations);
}
}
if (dimensionsVar2 > 1) {
if (kdTreeVar2Conditional == null) {
kdTreeVar2Conditional = new KdTree(
new int[] {dimensionsVar2, dimensionsCond},
new double[][][] {var2Observations, condObservations},
observationSetIndices, observationTimePoints);
new double[][][] {var2Observations, condObservations});
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,
observationSetIndices, observationTimePoints);
uniNNSearcherVar2 = new UnivariateNearestNeighbourSearcher(var2Observations);
}
}
if ((nnSearcherConditional == null) && (dimensionsCond > 0)) {
nnSearcherConditional = NearestNeighbourSearcher.create(condObservations,
observationSetIndices, observationTimePoints);
nnSearcherConditional = NearestNeighbourSearcher.create(condObservations);
nnSearcherConditional.setNormType(normType);
}
}
@ -912,136 +911,4 @@ 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));
}
}

8
java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java
View File

@ -18,13 +18,19 @@
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
@ -334,7 +340,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 the new sample:
// finding the kth closest neighbour for point t:
long methodStartTime = Calendar.getInstance().getTimeInMillis();
PriorityQueue<NeighbourNodeData> nnPQ =
kdTreeJoint.findKNearestNeighbours(k,

126
java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java
View File

@ -295,6 +295,13 @@ 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 +
@ -722,18 +729,15 @@ 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},
observationSetIndices, observationTimePoints);
new double[][][] {sourceObservations, destObservations});
kdTreeJoint.setNormType(normType);
}
if (nnSearcherSource == null) {
nnSearcherSource = NearestNeighbourSearcher.create(sourceObservations,
observationSetIndices, observationTimePoints);
nnSearcherSource = NearestNeighbourSearcher.create(sourceObservations);
nnSearcherSource.setNormType(normType);
}
if (nnSearcherDest == null) {
nnSearcherDest = NearestNeighbourSearcher.create(destObservations,
observationSetIndices, observationTimePoints);
nnSearcherDest = NearestNeighbourSearcher.create(destObservations);
nnSearcherDest.setNormType(normType);
}
}
@ -909,8 +913,7 @@ 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;
@SuppressWarnings("unused") // Might be used by caller later
public static final int RETURN_ARRAY_LENGTH = 4;
public static final int RETURN_ARRAY_LENGTH = 4;
public MiKraskovThreadRunner(
MutualInfoCalculatorMultiVariateKraskov miCalc,
@ -1107,110 +1110,5 @@ 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));
}
}
}

4
java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java
View File

@ -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 the new sample:
// finding the kth closest neighbour for point t:
PriorityQueue<NeighbourNodeData> nnPQ =
kdTreeJoint.findKNearestNeighbours(k,
new double[][] {newVar1Observations[t], newVar2Observations[t]});

61
java/source/infodynamics/measures/continuous/kraskov/TransferEntropyCalculatorKraskov.java
View File

@ -21,7 +21,6 @@ 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
@ -251,64 +250,4 @@ 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);
}
}

89
java/source/infodynamics/measures/spiking/integration/TransferEntropyCalculatorSpikingIntegration.java
View File

@ -1,6 +1,7 @@
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;
@ -15,6 +16,8 @@ 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;
@ -100,9 +103,6 @@ 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,6 +400,7 @@ 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;
@ -416,15 +417,13 @@ 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];
arrayedJointEmbeddingsFromSpikes = new double[conditioningEmbeddingsFromSpikes.size()][numDestPastIntervals +
double[][] 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];
@ -698,51 +697,19 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
return new distanceAndNumPoints(maxDistance, i);
}
/**
* 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;
return computeAverageLocalOfObservations(kdTreeJointAtSpikes, jointEmbeddingsFromSpikes);
}
/*
* 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.
*/
protected TERateResults computeAverageLocalOfObservations(KdTree actualKdTreeJointAtSpikes, Vector<double[]> actualJointEmbeddingsFromSpikes) throws Exception {
public double 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];
@ -845,18 +812,16 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
radiusConditioningSamples = Math.sqrt(radiusConditioningSamples);
}
localValues[indexInLocals] = (MathsUtils.digamma(kJointSpikes) - MathsUtils.digamma(kJointSamples) +
currentSum += (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))));
currentSum += localValues[indexInLocals];
if (Double.isNaN(localValues[indexInLocals])) {
+ ((numDestPastIntervals + numCondPastIntervals) * (Math.log(radiusConditioningSpikes) - Math.log(radiusConditioningSamples))));
if (Double.isNaN(currentSum)) {
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;
@ -864,9 +829,7 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
timeSum += time;
}
currentSum /= timeSum;
TERateResults teRateResults = new TERateResults(currentSum, localValues, conditioningEmbeddingsFromSpikes.size(), timeSum);
return teRateResults;
return currentSum;
}
@ -894,6 +857,7 @@ 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;
@ -921,7 +885,7 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
KdTree resampledKdTreeConditioningAtSamples = new KdTree(arrayedResampledConditioningEmbeddingsFromSamples);
resampledKdTreeConditioningAtSamples.setNormType(normType);
Vector<double[]> conditionallyPermutedJointEmbeddingsFromSpikes = new Vector<double[]>(jointEmbeddingsFromSpikes);
Vector<double[]> conditionallyPermutedJointEmbeddingsFromSpikes = new Vector(jointEmbeddingsFromSpikes);
Vector<Integer> usedIndices = new Vector<Integer>();
for (int i = 0; i < conditionallyPermutedJointEmbeddingsFromSpikes.size(); i++) {
@ -956,28 +920,13 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
}
KdTree conditionallyPermutedKdTreeJointFromSpikes = new KdTree(arrayedConditionallyPermutedJointEmbeddingsFromSpikes);
conditionallyPermutedKdTreeJointFromSpikes.setNormType(normType);
TERateResults teRateResults = computeAverageLocalOfObservations(conditionallyPermutedKdTreeJointFromSpikes,
conditionallyPermutedJointEmbeddingsFromSpikes);
surrogateTEValues[permutationNumber] = teRateResults.meanRate;
surrogateTEValues[permutationNumber] = computeAverageLocalOfObservations(conditionallyPermutedKdTreeJointFromSpikes,
conditionallyPermutedJointEmbeddingsFromSpikes);
}
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)
*
@ -986,8 +935,8 @@ public class TransferEntropyCalculatorSpikingIntegration implements TransferEntr
*/
@Override
public SpikingLocalInformationValues computeLocalOfPreviousObservations() throws Exception {
TERateResults teRateResults = computeAverageLocalOfObservations(kdTreeJointAtSpikes, jointEmbeddingsFromSpikes);
return new SpikingTELocalValues(teRateResults.localContributions);
// TODO Auto-generated method stub
return null;
}

71
java/source/infodynamics/utils/KdTree.java
View File

@ -113,20 +113,6 @@ 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.
@ -140,29 +126,7 @@ 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;
@ -219,18 +183,6 @@ 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;
}
/**
@ -715,7 +667,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 in the same observation set. 0 means exclude only sampleIndex itself.
* on either side of sampleIndex. 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
@ -748,8 +700,8 @@ public class KdTree extends NearestNeighbourSearcher {
// (will not throw an Exception if the PQ is empty)
NeighbourNodeData furthestCached = currentKBest.peek();
if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex])
|| (Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
if (((node.indexOfThisPoint - sampleIndex > dynCorrExclTime)
|| (node.indexOfThisPoint - 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
@ -1121,8 +1073,7 @@ public class KdTree extends NearestNeighbourSearcher {
absDistOnThisDim = distOnThisDim * distOnThisDim;
}
if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex])
|| (Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
((absDistOnThisDim < r) ||
( allowEqualToR && (absDistOnThisDim == r)))) {
// Preliminary check says we need to compute the full distance
@ -1640,8 +1591,7 @@ public class KdTree extends NearestNeighbourSearcher {
absDistOnThisDim = distOnThisDim * distOnThisDim;
}
if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) ||
(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
((absDistOnThisDim < r) ||
( allowEqualToR && (absDistOnThisDim == r)))) {
// Preliminary check says we need to compute the full distance
@ -1788,8 +1738,7 @@ public class KdTree extends NearestNeighbourSearcher {
absDistOnThisDim = distOnThisDim * distOnThisDim;
}
if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) ||
(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
((absDistOnThisDim < r) ||
( allowEqualToR && (absDistOnThisDim == r)))) {
// Preliminary check says we need to compute the full distance
@ -1963,8 +1912,7 @@ public class KdTree extends NearestNeighbourSearcher {
}
if (testResultsForGivenVariable[node.indexOfThisPoint] &&
((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) ||
(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
(Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
((absDistOnThisDim < r) ||
( allowEqualToR && (absDistOnThisDim == r)))) {
// Preliminary check says we need to compute the full distance
@ -2626,8 +2574,7 @@ public class KdTree extends NearestNeighbourSearcher {
absDistOnThisDim = distOnThisDim * distOnThisDim;
}
if (((observationSetIndices[node.indexOfThisPoint] != observationSetIndices[sampleIndex]) ||
(Math.abs(observationTimePoints[node.indexOfThisPoint] - observationTimePoints[sampleIndex]) > dynCorrExclTime)) &&
if ((Math.abs(node.indexOfThisPoint - sampleIndex) > dynCorrExclTime) &&
((absDistOnThisDim < rs[variableNumber]) ||
( allowEqualToR && (absDistOnThisDim == rs[variableNumber])))) {
// Preliminary check says we need to compute the full distance

61
java/source/infodynamics/utils/NearestNeighbourSearcher.java
View File

@ -40,19 +40,6 @@ 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
@ -64,37 +51,17 @@ 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),
observationSetIndices, observationTimePoints);
return new UnivariateNearestNeighbourSearcher(MatrixUtils.selectColumn(data, 0));
} else {
return new KdTree(data, observationSetIndices, observationTimePoints);
return new KdTree(data);
}
}
/**
* Factory method to construct the searcher from a set of double[][][] data.
*
@ -104,29 +71,11 @@ 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),
observationSetIndices, observationTimePoints);
return new UnivariateNearestNeighbourSearcher(MatrixUtils.selectColumn(data[0], 0));
} else {
return new KdTree(dimensions, data, observationSetIndices, observationTimePoints);
return new KdTree(dimensions, data);
}
}

62
java/source/infodynamics/utils/UnivariateNearestNeighbourSearcher.java
View File

@ -63,25 +63,15 @@ 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), observationSetIndices, observationTimePoints);
this(MatrixUtils.selectColumn(data, 0));
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 {
public UnivariateNearestNeighbourSearcher(double[] data) throws Exception {
this.originalDataSet = data;
numObservations = data.length;
if (numObservations <= 1) {
@ -112,18 +102,6 @@ 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;
}
/**
@ -274,8 +252,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
lowerCandidate >= 0;
lowerCandidate--) {
indexOfLowerCandidate = sortedArrayIndices[lowerCandidate];
if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfLowerCandidate]) ||
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfLowerCandidate]) > dynCorrExclTime)) {
if (Math.abs(sampleIndex - indexOfLowerCandidate) > dynCorrExclTime) {
// This sample is outside the dynamic correlation exclusion window
break;
}
@ -289,8 +266,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
upperCandidate <= numObservations - 1;
upperCandidate++) {
indexOfUpperCandidate = sortedArrayIndices[upperCandidate];
if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfUpperCandidate]) ||
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfUpperCandidate]) > dynCorrExclTime)) {
if (Math.abs(sampleIndex - indexOfUpperCandidate) > dynCorrExclTime) {
// This sample is outside the dynamic correlation exclusion window
break;
}
@ -325,8 +301,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
upperCandidate <= numObservations - 1;
upperCandidate++) {
indexOfUpperCandidate = sortedArrayIndices[upperCandidate];
if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfUpperCandidate]) ||
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfUpperCandidate]) > dynCorrExclTime)) {
if (Math.abs(sampleIndex - indexOfUpperCandidate) > dynCorrExclTime) {
// This sample is outside the dynamic correlation exclusion window
break;
}
@ -342,8 +317,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
lowerCandidate >= 0;
lowerCandidate--) {
indexOfLowerCandidate = sortedArrayIndices[lowerCandidate];
if ((observationSetIndices[sampleIndex] != observationSetIndices[indexOfLowerCandidate]) ||
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[indexOfLowerCandidate]) > dynCorrExclTime)) {
if (Math.abs(sampleIndex - indexOfLowerCandidate) > dynCorrExclTime) {
// This sample is outside the dynamic correlation exclusion window
break;
}
@ -500,8 +474,7 @@ 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 ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
// Can't count this point, but keep checking:
continue;
}
@ -517,8 +490,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
}
// Next check the points with larger data values:
for (int i = indexInSortedArray + 1; i < numObservations; i++) {
if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
// Can't count this point, but keep checking:
continue;
}
@ -775,8 +747,7 @@ 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 ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
// Can't count this point, but keep checking:
continue;
}
@ -793,8 +764,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
}
// Next check the points with larger data values:
for (int i = indexInSortedArray + 1; i < numObservations; i++) {
if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
// Can't count this point, but keep checking:
continue;
}
@ -823,8 +793,7 @@ 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 ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
// Can't count this point, but keep checking:
continue;
}
@ -843,8 +812,7 @@ public class UnivariateNearestNeighbourSearcher extends NearestNeighbourSearcher
}
// Next check the points with larger data values:
for (int i = indexInSortedArray + 1; i < numObservations; i++) {
if ((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime)) {
if (Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime) {
// Can't count this point, but keep checking:
continue;
}
@ -900,8 +868,7 @@ 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]] ||
((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime))) {
(Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime)) {
// Can't count this point, but keep checking:
continue;
}
@ -922,8 +889,7 @@ 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]] ||
((observationSetIndices[sampleIndex] == observationSetIndices[sortedArrayIndices[i]]) &&
(Math.abs(observationTimePoints[sampleIndex] - observationTimePoints[sortedArrayIndices[i]]) <= dynCorrExclTime))) {
(Math.abs(sampleIndex - sortedArrayIndices[i]) <= dynCorrExclTime)) {
// Can't count this point, but keep checking:
continue;
}

39
java/unittests/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoMultiVariateTester.java
View File

@ -159,40 +159,6 @@ 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();
@ -223,7 +189,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, 0.0000001);
assertEquals(distroNotBiasCorrected.pValue, distroBiasCorrected.pValue);
}
protected int timeStepsDepCheck = 100;
@ -429,7 +395,7 @@ public class ConditionalMutualInfoMultiVariateTester extends
// - both dimensions are copied
double[][] sourceData = rg.generateNormalData(timeStepsDepCheck, dimensions,
0, 1);
double[][] condData = rg.generateNormalData(timeStepsDepCheck, conditionalDims,
double[][] condData = rg.generateNormalData(timeStepsDepCheck, dimensions,
0, 1);
double[][] destData = MatrixUtils.arrayCopy(condData);
@ -481,5 +447,4 @@ public class ConditionalMutualInfoMultiVariateTester extends
condMiCalc.setObservations(MatrixUtils.selectColumns(destData, 0, 1), sourceData, condData);
assertTrue(Double.isInfinite(condMiCalc.computeAverageLocalOfObservations()));
}
}

73
java/unittests/infodynamics/measures/continuous/kraskov/ActiveInfoStorageTester.java
View File

@ -5,7 +5,6 @@ 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 {
@ -117,74 +116,4 @@ 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++;
}
}
}}
}

185
java/unittests/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoMultiVariateTester.java
View File

@ -21,7 +21,6 @@ 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 {
@ -98,11 +97,10 @@ 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 double[] checkTEForGivenData(double[][] var1, double[][] var2,
protected void checkTEForGivenData(double[][] var1, double[][] var2,
int[] kNNs, double[] expectedResults) throws Exception {
return checkTEForGivenData(var1, var2, 1, 1, kNNs, expectedResults);
checkTEForGivenData(var1, var2, 1, 1, kNNs, expectedResults);
}
/**
@ -116,38 +114,12 @@ 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 double[] checkTEForGivenData(double[][] var1, double[][] var2,
protected void 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.
@ -175,9 +147,7 @@ public class ConditionalMutualInfoMultiVariateTester
condMiCalc.setProperty(
ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
NUM_THREADS_TO_USE);
condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE,
Double.toString(noiseLevel));
condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_NOISE_SEED, noiseSeed);
condMiCalc.setProperty(ConditionalMutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency of results for unit test
condMiCalc.initialise(var1[0].length * historyL,
var2[0].length, var2[0].length * historyK);
// Construct the joint vectors of the source states
@ -222,10 +192,8 @@ 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, tolerance);
errors[kIndex] = condMi - expectedResults[kIndex];
assertEquals(expectedResults[kIndex], condMi, 0.000001);
}
return errors;
}
/**
@ -731,7 +699,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, (double[][]) null);
condMiCalc.setObservations(var1, var2, null);
// condMiCalc.setDebug(true);
double mi = condMiCalc.computeAverageLocalOfObservations();
// condMiCalc.setDebug(false);
@ -768,143 +736,4 @@ 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);
}
}

231
java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java
View File

@ -21,7 +21,6 @@ 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 {
@ -96,23 +95,6 @@ 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
*
@ -120,20 +102,13 @@ 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 double[] checkMIForGivenData(double[][] var1, double[][] var2,
int[] kNNs, double[] expectedResults,
double noiseLevel, String noiseSeed, double tolerance) throws Exception {
protected void checkMIForGivenData(double[][] var1, double[][] var2,
int[] kNNs, double[] expectedResults) 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];
@ -146,9 +121,7 @@ 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,
Double.toString(noiseLevel));
miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_NOISE_SEED, noiseSeed);
miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests
miCalc.initialise(var1[0].length, var2[0].length);
miCalc.setObservations(var1, var2);
miCalc.setDebug(true);
@ -157,10 +130,10 @@ public class MutualInfoMultiVariateTester
System.out.printf("k=%d: Average MI %.8f (expected %.8f)\n",
k, mi, expectedResults[kIndex]);
assertEquals(expectedResults[kIndex], mi, tolerance);
errors[kIndex] = 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);
}
return errors;
}
/**
@ -679,195 +652,5 @@ 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);
}
}

72
java/unittests/infodynamics/measures/continuous/kraskov/TransferEntropyTester.java
View File

@ -737,76 +737,4 @@ 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++;
}
}
}
}

66
java/unittests/infodynamics/utils/KdTreeTest.java
View File

@ -422,72 +422,6 @@ 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;

62
java/unittests/infodynamics/utils/UnivariateNearestNeighbourTest.java
View File

@ -625,66 +625,4 @@ 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));
}
}
}