From b3192e225ffbcbb7aa87deb7dde8e840f9b0d5d4 Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Sun, 24 Jan 2021 21:54:09 +0000 Subject: [PATCH 1/8] Added convenience methods to multiply an array times a constant value. --- .../infodynamics/utils/MatrixUtils.java | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/java/source/infodynamics/utils/MatrixUtils.java b/java/source/infodynamics/utils/MatrixUtils.java index 98cc551..bebc867 100755 --- a/java/source/infodynamics/utils/MatrixUtils.java +++ b/java/source/infodynamics/utils/MatrixUtils.java @@ -819,6 +819,51 @@ public class MatrixUtils { return returnValues; } + /** + * Multiplies all items in an array times a constant value + * + * @param array + * @param value + * @return array * constant value + */ + public static int[] multiply(int[] array, int value) throws Exception { + int[] returnValues = new int[array.length]; + for (int i = 0; i < returnValues.length; i++) { + returnValues[i] = array[i] * value; + } + return returnValues; + } + + /** + * Multiplies all items in an array times a constant value + * + * @param array + * @param value + * @return array * constant value + */ + public static double[] multiply(int[] array, double value) throws Exception { + double[] returnValues = new double[array.length]; + for (int i = 0; i < returnValues.length; i++) { + returnValues[i] = array[i] * value; + } + return returnValues; + } + + /** + * Multiplies all items in an array times a constant value + * + * @param array + * @param value + * @return array * constant value + */ + public static double[] multiply(double[] array, double value) throws Exception { + double[] returnValues = new double[array.length]; + for (int i = 0; i < returnValues.length; i++) { + returnValues[i] = array[i] * value; + } + return returnValues; + } + /** * Return the matrix product A x B * From 06cc6d649de7a68a693b6bd4ccc948f82be267d0 Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Sun, 24 Jan 2021 21:54:42 +0000 Subject: [PATCH 2/8] Added `select*` methods for `int[][]` matrices. --- .../infodynamics/utils/MatrixUtils.java | 78 +++++++++++++++++++ 1 file changed, 78 insertions(+) diff --git a/java/source/infodynamics/utils/MatrixUtils.java b/java/source/infodynamics/utils/MatrixUtils.java index bebc867..5b14833 100755 --- a/java/source/infodynamics/utils/MatrixUtils.java +++ b/java/source/infodynamics/utils/MatrixUtils.java @@ -1580,6 +1580,84 @@ public class MatrixUtils { return data; } + /** + * Extract the required columns from the matrix + * + * @param matrix + * @param fromCol + * @param cols + * @return + */ + public static int[][] selectColumns(int matrix[][], + int fromCol, int cols) { + int[][] data = new int[matrix.length][cols]; + for (int r = 0; r < matrix.length; r++) { + for (int cIndex = 0; cIndex < cols; cIndex++) { + data[r][cIndex] = matrix[r][cIndex + fromCol]; + } + } + return data; + } + + /** + * Extract the required columns from the matrix + * + * @param matrix + * @param columns + * @return + */ + public static int[][] selectColumns(int matrix[][], int columns[]) { + int[][] data = new int[matrix.length][columns.length]; + for (int r = 0; r < matrix.length; r++) { + for (int cIndex = 0; cIndex < columns.length; cIndex++) { + data[r][cIndex] = matrix[r][columns[cIndex]]; + } + } + return data; + } + + /** + * Extract the required columns from the matrix + * + * @param matrix + * @param includeColumnFlags + * @return + */ + public static int[][] selectColumns(int matrix[][], boolean includeColumnFlags[]) { + Vector v = new Vector(); + + for (int i = 0; i < includeColumnFlags.length; i++) { + if (includeColumnFlags[i]) { + v.add(i); + } + } + int[][] data = new int[matrix.length][v.size()]; + for (int r = 0; r < matrix.length; r++) { + for (int outputColumnIndex = 0; outputColumnIndex < v.size(); outputColumnIndex++) { + int outputColumn = v.get(outputColumnIndex); + data[r][outputColumnIndex] = matrix[r][outputColumn]; + } + } + return data; + } + + /** + * Extract the required columns from the matrix + * + * @param matrix + * @param columns + * @return + */ + public static int[][] selectColumns(int matrix[][], List columns) { + int[][] data = new int[matrix.length][columns.size()]; + for (int r = 0; r < matrix.length; r++) { + for (int cIndex = 0; cIndex < columns.size(); cIndex++) { + data[r][cIndex] = matrix[r][columns.get(cIndex)]; + } + } + return data; + } + /** * Extract the required rows from the matrix * From 403745800e1da8a35037eb4ef23e76e1a09285ae Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Sun, 24 Jan 2021 21:56:26 +0000 Subject: [PATCH 3/8] Added abstract parent class for multivariate IT measures on continuous data. --- ...ltiVariateInfoMeasureCalculatorCommon.java | 575 ++++++++++++++++++ 1 file changed, 575 insertions(+) create mode 100644 java/source/infodynamics/measures/continuous/MultiVariateInfoMeasureCalculatorCommon.java diff --git a/java/source/infodynamics/measures/continuous/MultiVariateInfoMeasureCalculatorCommon.java b/java/source/infodynamics/measures/continuous/MultiVariateInfoMeasureCalculatorCommon.java new file mode 100644 index 0000000..9e9af2f --- /dev/null +++ b/java/source/infodynamics/measures/continuous/MultiVariateInfoMeasureCalculatorCommon.java @@ -0,0 +1,575 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous; + +import infodynamics.utils.EmpiricalMeasurementDistribution; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; + +import java.util.Arrays; +import java.util.Random; +import java.util.Vector; + +/** + * Implements a base class with common functionality for child class + * implementations of multivariate information measures via various estimators. + * + *

Multivariate information measures are functionals of probability + * distributions over R^n, and typical examples include multi-information + * (a.k.a. total correlation), dual total correlation, O-information, and connected + * information.

+ * + *

These measures can be computed via different kinds of estimators, such as + * linear-gaussian, KSG estimators, etc (see the child classes linked above). + *

+ * + *

+ * Usage of the child classes is intended to follow this paradigm: + *

+ *
    + *
  1. Construct the calculator;
    2. + *
  2. Set properties using {@link #setProperty(String, String)};
    3. + *
  3. Initialise the calculator using {@link #initialise(int)};
    4. + *
  4. Provide the observations/samples for the calculator + * to set up the PDFs, using: + *
      + *
    • {@link #setObservations(double[][])} + * for calculations based on single time-series, OR
      • + *
    • The following sequence:
        + *
      1. {@link #startAddObservations()}, then
        2. + *
      2. One or more calls to {@link #addObservations(double[][])} or + * {@link #addObservation(double[])}, then
        3. + *
      3. {@link #finaliseAddObservations()};
        4. + *
      • + *
    + *
  5. Compute the required quantities, being one or more of: + *
      + *
    • the average measure: + * {@link #computeAverageLocalOfObservations()};
      • + *
    • the local values for these samples: + * {@link #computeLocalOfPreviousObservations()}
      • + *
    • local values for a specific set of samples: + * {@link #computeLocalUsingPreviousObservations(double[][])}
      • + *
    + *
    6. + *
  6. + * Return to step 2 or 3 to re-use the calculator on a new data set. + *
    7. + *
+ *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public abstract class MultiVariateInfoMeasureCalculatorCommon + implements InfoMeasureCalculatorContinuous { + + /** + * Number of joint variables to consider + */ + protected int dimensions = 1; + /** + * Number of samples supplied + */ + protected int totalObservations = 0; + /** + * Whether we are in debug mode + */ + protected boolean debug = false; + /** + * Cached supplied observations + */ + protected double[][] observations; + /** + * Set of individually supplied observations + */ + protected Vector individualObservations; + /** + * Whether the user has supplied more than one (disjoint) set of samples + */ + protected boolean addedMoreThanOneObservationSet; + /** + * Whether the measure has been computed for the latest supplied data + */ + protected boolean isComputed = false; + /** + * Cached last the measure value calculated + */ + protected double lastAverage; + + /** + * Whether to normalise incoming values + */ + protected boolean normalise = true; + + /** + * Property name for whether to normalise incoming values to mean 0, + * standard deviation 1 (default true) + */ + public static final String PROP_NORMALISE = "NORMALISE"; + + + /** + * Initialise the calculator for (re-)use, with the existing + * (or default) values of parameters, with number of + * joint variables specified. + * Clears any PDFs of previously supplied observations. + * + * @param dimensions the number of joint variables to consider + */ + public void initialise() { + initialise(dimensions); + } + + /** + * Initialise the calculator for (re-)use, with the existing + * (or default) values of parameters, with number of + * joint variables specified. + * Clears an PDFs of previously supplied observations. + * + * @param dimensions the number of joint variables to consider + */ + public void initialise(int dimensions) { + this.dimensions = dimensions; + lastAverage = 0.0; + totalObservations = 0; + isComputed = false; + observations = null; + addedMoreThanOneObservationSet = false; + } + + /** + * Set properties for the calculator. + * New property values are not guaranteed to take effect until the next call + * to an initialise method. + * + *

Valid property names, and what their + * values should represent, include:

+ *
    + *
  • {@link #PROP_NORMALISE} -- whether to normalise the incoming variables + * to mean 0, standard deviation 1, or not (default false).
    • + *
+ * + *

Unknown property values are ignored.

+ * + * @param propertyName name of the property + * @param propertyValue value of the property + * @throws Exception for invalid property values + */ + public void setProperty(String propertyName, String propertyValue) + throws Exception { + boolean propertySet = true; + if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) { + normalise = Boolean.parseBoolean(propertyValue); + } else { + // No property was set here + propertySet = false; + } + if (debug && propertySet) { + System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName + + " to " + propertyValue); + } + } + + public String getProperty(String propertyName) throws Exception { + if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) { + return Boolean.toString(normalise); + } else { + return null; + } + } + + /** + * Sets a single series from which to compute the PDF for the instantiated measure.. + * Cannot be called in conjunction with other methods for setting/adding + * observations. + * + *

The supplied series may be a time-series, or may be simply + * a set of separate observations + * without a time interpretation.

+ * + *

Should only be called once, the last call contains the + * observations that are used (they are not accumulated).

+ * + * @param observations series of multivariate observations + * (first index is time or observation index, second is variable number) + * @throws Exception + */ + public void setObservations(double[][] observations) throws Exception { + startAddObservations(); + addObservations(observations); + finaliseAddObservations(); + addedMoreThanOneObservationSet = false; + } + + /** + * Signal that we will add in the samples for computing the PDF + * from several disjoint time-series or trials via calls to + * {@link #addObservation(double[])} or {@link #addObservations(double[][])} + * rather than {@link #setDebug(boolean)}. + */ + public void startAddObservations() { + individualObservations = new Vector(); + } + + /** + *

Adds a new (single) observation 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.

+ * + *

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).

+ * + * @param observation a single multivariate observation + * (index is variable number) + */ + public void addObservation(double[] observation) { + if (individualObservations.size() > 0) { + addedMoreThanOneObservationSet = true; + } + individualObservations.add(observation); + } + + /** + *

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.

+ * + *

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).

+ * + * @param observations series of multivariate observations + * (first index is time or observation index, second is variable number) + */ + public void addObservations(double[][] observations) { + // This implementation is not particularly efficient, + // however it will suffice for now. + for (int s = 0; s < observations.length; s++) { + addObservation(observations[s]); + } + } + + /** + * {@inheritDoc} + * + * This class provides a basic implementation, generating + * the internal set of samples in observations; child classes + * should then process these observations as required. + * + */ + public void finaliseAddObservations() throws Exception { + observations = new double[individualObservations.size()][]; + for (int t = 0; t < observations.length; t++) { + observations[t] = individualObservations.elementAt(t); + } + // Allow vector to be reclaimed + individualObservations = null; + + if (observations[0].length != dimensions) { + throw new Exception("Incorrect number of dimensions " + observations[0].length + + " in supplied observations (expected " + dimensions + ")"); + } + totalObservations = observations.length; + } + + /** + * Generate a resampled distribution of what the measure would look like, + * under a null hypothesis that the individual values of each + * variable in the + * samples have no relation to each other. + * That is, we destroy the p(x,y,z,..) correlations, while + * retaining the p(x), p(y),.. marginals, to check how + * significant this measure actually was. + * + *

See Section II.E "Statistical significance testing" of + * the JIDT paper below for a description of how this is done for MI, + * we are extending that here. + *

+ * + *

Note that if several disjoint time-series have been added + * as observations using {@link #addObservations(double[])} etc., + * then these separate "trials" will be mixed up in the generation + * of surrogates here.

+ * + *

This method (in contrast to {@link #computeSignificance(int[][][])}) + * creates random shufflings of the next values for the surrogate MultiInfo + * calculations.

+ * + * @param numPermutationsToCheck number of surrogate samples to permute + * to generate the distribution. + * @return the distribution of surrogate measure values under this null hypothesis. + * @see "J.T. Lizier, 'JIDT: An information-theoretic + * toolkit for studying the dynamics of complex systems', 2014." + * @throws Exception + */ + public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck) throws Exception { + // Generate the re-ordered indices: + RandomGenerator rg = new RandomGenerator(); + int[][][] newOrderings = new int[numPermutationsToCheck][][]; + // Generate numPermutationsToCheck * (dimensions-1) permutations of 0 .. data.length-1 + for (int n = 0; n < numPermutationsToCheck; n++) { + // (Not necessary to check for distinct random perturbations) + newOrderings[n] = rg.generateRandomPerturbations(totalObservations, dimensions-1); + } + return computeSignificance(newOrderings); + } + + /** + * Generate a resampled distribution of what the measure would look like, + * under a null hypothesis that the individual values of each + * variable in the + * samples have no relation to eachother. + * That is, we destroy the p(x,y,z,..) correlations, while + * retaining the p(x), p(y),.. marginals, to check how + * significant this measure actually was. + * + *

See Section II.E "Statistical significance testing" of + * the JIDT paper below for a description of how this is done for MI, + * we are extending that here. + *

+ * + *

Note that if several disjoint time-series have been added + * as observations using {@link #addObservations(double[])} etc., + * then these separate "trials" will be mixed up in the generation + * of surrogates here.

+ * + *

This method (in contrast to {@link #computeSignificance(int)}) + * allows the user to specify how to construct the surrogates, + * such that repeatable results may be obtained.

+ * + * @param newOrderings a specification of how to shuffle the values + * to create the surrogates to generate the distribution with. The first + * index is the permutation number (i.e. newOrderings.length is the number + * of surrogate samples we use to bootstrap to generate the distribution here.) + * The second index is the variable number (minus 1, since we don't reorder + * the first variable), + * Each array newOrderings[i][v] should be an array of length N (where + * would be the value returned by {@link #getNumObservations()}), + * containing a permutation of the values in 0..(N-1). + * @return the distribution of surrogate measure values under this null hypothesis. + * @see "J.T. Lizier, 'JIDT: An information-theoretic + * toolkit for studying the dynamics of complex systems', 2014." + * @throws Exception where the length of each permutation in newOrderings + * is not equal to the number N samples that were previously supplied. + */ + public EmpiricalMeasurementDistribution computeSignificance(int[][][] newOrderings) throws Exception { + + int numPermutationsToCheck = newOrderings.length; + if (!isComputed) { + computeAverageLocalOfObservations(); + } + + // Store the real observations and their measure value: + double actualMeasure = lastAverage; + + EmpiricalMeasurementDistribution measDistribution = new EmpiricalMeasurementDistribution(numPermutationsToCheck); + + int countWhereSurrogateIsMoreSignificantThanOriginal = 0; + for (int i = 0; i < numPermutationsToCheck; i++) { + // Compute the measure under this reordering + double newMeasure = computeAverageLocalOfObservations(newOrderings[i]); + measDistribution.distribution[i] = newMeasure; + if (debug){ + System.out.println("New measure value was " + newMeasure); + } + if (newMeasure >= actualMeasure) { + countWhereSurrogateIsMoreSignificantThanOriginal++; + } + } + + // Restore the actual measure and the observations + lastAverage = actualMeasure; + + // And return the significance + measDistribution.pValue = (double) countWhereSurrogateIsMoreSignificantThanOriginal / (double) numPermutationsToCheck; + measDistribution.actualValue = actualMeasure; + return measDistribution; + } + + /** + * Compute what the measure would look like were all time series (bar the first) + * reordered as per the array of time indices in newOrdering. + * + *

The reordering array contains the reordering for each marginal variable + * (first index). The user should ensure that all values 0..N-1 are + * represented exactly once in the array reordering and that no other values + * are included here.

+ * + *

Note that if several disjoint time-series have been added as + * observations using {@link #addObservations(double[])} etc., then these + * separate "trials" will be mixed up in the generation of a shuffled source + * series here.

+ * + *

This method is primarily intended for use in {@link + * #computeSignificance(int[][])} however has been made public in case users + * wish to access it.

+ * + * @param newOrdering the specific permuted new orderings to use. First index + * is the variable number (minus 1, since we don't reorder the first + * variable), second index is the time step, the value is the reordered time + * step to use for that variable at the given time step. The values must be + * an array of length N (where + * would be the value returned by {@link #getNumObservations()}), containing + * a permutation of the values in 0..(N-1). If null, no reordering is + * performed. + * @return what the average measure would look like under this reordering + * @throws Exception + */ + public double computeAverageLocalOfObservations(int[][] newOrdering) + throws Exception { + + if (newOrdering == null) { + return computeAverageLocalOfObservations(); + } + + // Take a clone of the object to compute the measure of the surrogates: + // (this is a shallow copy, it doesn't make new copies of all + // the arrays) + MultiVariateInfoMeasureCalculatorCommon surrogateCalculator = + (MultiVariateInfoMeasureCalculatorCommon) this.clone(); + + // Generate a new re-ordered source data + double[][] shuffledData = + MatrixUtils.reorderDataForVariables( + observations, newOrdering); + // Perform new initialisations + surrogateCalculator.initialise(dimensions); + // Set new observations + surrogateCalculator.setObservations(shuffledData); + // Compute the MI + return surrogateCalculator.computeAverageLocalOfObservations(); + } + + /** + * Calculates the local measure at every sample provided since the last time the + * calculator was initialised. + * + * @return the "time-series" of local measure values in nats (not bits!) + * @throws Exception + */ + public double[] computeLocalOfPreviousObservations() throws Exception { + // Cannot do if observations haven't been set + if (observations == null) { + throw new Exception("Cannot compute local values of previous observations " + + "if they have not been set!"); + } + + return computeLocalUsingPreviousObservations(observations); + } + + /** + * Compute the local measure values for each of the + * supplied samples in states. + * + *

PDFs are computed using all of the previously supplied + * observations, but not those in states + * (unless they were + * some of the previously supplied samples).

+ * + * @param states series of multivariate observations + * (first index is time or observation index, second is variable number) + * @return the series of local measure values. + * @throws Exception + */ + public abstract double[] computeLocalUsingPreviousObservations(double states[][]) + throws Exception; + + /** + * Shortcut method to initialise the calculator, set observations and compute + * the average measure in one line. + * + * @param new_observations series of multivariate observations + * (first index is time or observation index, second is variable number) + */ + public double compute(double[][] new_observations) throws Exception { + initialise(new_observations[0].length); + setObservations(new_observations); + return computeAverageLocalOfObservations(); + } + + /** + * Shortcut method to initialise the calculator, set observations and compute + * the local measure in one line. + * + * @param new_observations series of multivariate observations + * (first index is time or observation index, second is variable number) + */ + public double[] computeLocals(double[][] new_observations) throws Exception { + initialise(new_observations[0].length); + setObservations(new_observations); + return computeLocalOfPreviousObservations(); + } + + public int getNumObservations() throws Exception { + return totalObservations; + } + + public void setDebug(boolean debug) { + this.debug = debug; + } + + public double getLastAverage() { + return lastAverage; + } + + /** + * Returns an int[] array with all integers from 0 to + * N-1, except idx. + * + *

This method is primarily intended for internal use (to extract blocks + * of covariance matrices excluding one variable).

+ * + * @param idx index of integer to omit + * @param N upper limit of the integer array + */ + protected int[] allExcept(int idx, int N) { + boolean[] v = new boolean[N]; + Arrays.fill(v, true); + v[idx] = false; + + int[] v2 = new int[N - 1]; + int counter = 0; + for (int i = 0; i < N; i++) { + if (v[i]) { + v2[counter] = i; + counter++; + } + } + + return v2; + } + +} + From f67e2748192bb557f2da4a7bdafa46b31946fdff Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Sun, 24 Jan 2021 22:06:39 +0000 Subject: [PATCH 4/8] Added Gaussian implementation of various multivariate IT measures and unit tests. --- ...ualTotalCorrelationCalculatorGaussian.java | 122 ++++++++++++++ ...iVariateInfoMeasureCalculatorGaussian.java | 152 ++++++++++++++++++ .../gaussian/OInfoCalculatorGaussian.java | 134 +++++++++++++++ .../gaussian/SInfoCalculatorGaussian.java | 120 ++++++++++++++ ...alCorrelationCalculatorGaussianTester.java | 133 +++++++++++++++ .../OInfoCalculatorGaussianTester.java | 150 +++++++++++++++++ .../SInfoCalculatorGaussianTester.java | 129 +++++++++++++++ 7 files changed, 940 insertions(+) create mode 100644 java/source/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussian.java create mode 100644 java/source/infodynamics/measures/continuous/gaussian/MultiVariateInfoMeasureCalculatorGaussian.java create mode 100644 java/source/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussian.java create mode 100644 java/source/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussian.java create mode 100644 java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java create mode 100644 java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java create mode 100644 java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java diff --git a/java/source/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussian.java new file mode 100644 index 0000000..7d3c572 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussian.java @@ -0,0 +1,122 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2017, Joseph T. Lizier, Ipek Oezdemir and Pedro Mediano + * + * 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 . + */ + +package infodynamics.measures.continuous.gaussian; + +import infodynamics.utils.MatrixUtils; + +/** + *

Computes the differential dual total correlation (DTC) of a given multivariate + * double[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorGaussian}), + * assuming that the probability distribution function for these observations is + * a multivariate Gaussian distribution.

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorCommon}. + *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class DualTotalCorrelationCalculatorGaussian + extends MultiVariateInfoMeasureCalculatorGaussian { + + /** + * Constructor. + */ + public DualTotalCorrelationCalculatorGaussian() { + // Nothing to do + } + + /** + * {@inheritDoc} + * + * @return the average DTC in nats (not bits!) + * @throws Exception if not sufficient data have been provided, or if the + * supplied covariance matrix is invalid. + */ + public double computeAverageLocalOfObservations() throws Exception { + + if (covariance == null) { + throw new Exception("Cannot calculate DTC without having " + + "a covariance either supplied or computed via setObservations()"); + } + + if (!isComputed) { + double dtc = - (dimensions - 1)*Math.log(MatrixUtils.determinantSymmPosDefMatrix(covariance)); + for (int i = 0; i < dimensions; i++) { + int[] idx = allExcept(i, dimensions); + double[][] marginal_cov = MatrixUtils.selectRowsAndColumns(covariance, idx, idx); + dtc += Math.log(MatrixUtils.determinantSymmPosDefMatrix(marginal_cov)); + } + // This "0.5" comes from the entropy formula for Gaussians: h = 0.5*logdet(2*pi*e*Sigma) + lastAverage = 0.5*dtc;; + isComputed = true; + } + + return lastAverage; + } + + /** + * {@inheritDoc} + * + * @return the "time-series" of local DTC values in nats (not bits!) + * for the supplied states. + * @throws Exception if not sufficient data have been provided, or if the + * supplied covariance matrix is invalid. + */ + public double[] computeLocalUsingPreviousObservations(double[][] states) throws Exception { + + if ((means == null) || (covariance == null)) { + throw new Exception("Cannot compute local values without having means " + + "and covariance either supplied or computed via setObservations()"); + } + + EntropyCalculatorMultiVariateGaussian hCalc = new EntropyCalculatorMultiVariateGaussian(); + hCalc.initialise(dimensions); + hCalc.setCovarianceAndMeans(covariance, means); + double[] localValues = MatrixUtils.multiply(hCalc.computeLocalUsingPreviousObservations(states), -(dimensions - 1)); + + for (int i = 0; i < dimensions; i++) { + int[] idx = allExcept(i, dimensions); + double[][] marginal_cov = MatrixUtils.selectRowsAndColumns(covariance, idx, idx); + double[] marginal_means = MatrixUtils.select(means, idx); + double[][] marginal_state = MatrixUtils.selectColumns(states, idx); + + hCalc.initialise(dimensions - 1); + hCalc.setCovarianceAndMeans(marginal_cov, marginal_means); + double[] thisLocals = hCalc.computeLocalUsingPreviousObservations(marginal_state); + + localValues = MatrixUtils.add(localValues, thisLocals); + + } + + return localValues; + + } + +} + diff --git a/java/source/infodynamics/measures/continuous/gaussian/MultiVariateInfoMeasureCalculatorGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/MultiVariateInfoMeasureCalculatorGaussian.java new file mode 100644 index 0000000..f8d0933 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/gaussian/MultiVariateInfoMeasureCalculatorGaussian.java @@ -0,0 +1,152 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2017, Joseph T. Lizier, Ipek Oezdemir and Pedro Mediano + * + * 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 . + */ + +package infodynamics.measures.continuous.gaussian; + +import infodynamics.measures.continuous.MultiVariateInfoMeasureCalculatorCommon; +import infodynamics.utils.MatrixUtils; + +/** + *

Base class with common functionality for child class implementations of + * multivariate information measures on a given multivariate + * double[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorCommon}), + * assuming that the probability distribution function for these observations is + * a multivariate Gaussian distribution.

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorCommon}, + * with: + *

    + *
  • For constructors see the child classes.
    • + *
  • Further properties are defined in {@link #setProperty(String, String)}.
    • + *
  • Computed values are in nats, not bits!
    • + *
+ *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public abstract class MultiVariateInfoMeasureCalculatorGaussian + extends MultiVariateInfoMeasureCalculatorCommon { + + /** + * Covariance of the system. Can be calculated from supplied observations + * or supplied directly by the user. + */ + double[][] covariance = null; + + /** + * Means of the system. Can be calculated from supplied observations + * or supplied directly by the user. + */ + double[] means = null; + + /** + * Whether the current covariance matrix has been determined from data or + * supplied directly. This changes the approach to local measures and + * significance testing. + */ + boolean covFromObservations; + + @Override + public void finaliseAddObservations() throws Exception { + super.finaliseAddObservations(); + + this.means = MatrixUtils.means(observations); + setCovariance(MatrixUtils.covarianceMatrix(observations), true); + + return; + } + + /** + *

Set the covariance of the distribution for which we will compute the + * measure directly, without supplying observations.

+ * + *

See {@link #setCovariance(double[][], boolean)}. + * + * @param covariance covariance matrix of the system + * @throws Exception for covariance matrix not matching the expected dimensions, + * being non-square, asymmetric or non-positive definite + */ + public void setCovariance(double[][] covariance) throws Exception { + setCovariance(covariance, false); + } + + /** + *

Set the covariance of the distribution for which we will compute the + * measure.

+ * + *

This is an alternative to sequences of calls to {@link #setObservations(double[][])} or + * {@link #addObservations(double[][])} etc. + * Note that without setting any observations, you cannot later + * call {@link #computeLocalOfPreviousObservations()}.

+ * + * @param covariance covariance matrix of the system + * @param means mean of the system + */ + public void setCovarianceAndMeans(double[][] covariance, double[] means) + throws Exception { + this.means = means; + setCovariance(covariance, false); + } + + /** + *

Set the covariance of the distribution for which we will compute the + * measure.

+ * + *

This is an alternative to sequences of calls to {@link #setObservations(double[][])} or + * {@link #addObservations(double[][])} etc. + * Note that without setting any observations, you cannot later + * call {@link #computeLocalOfPreviousObservations()}, and without + * providing the means of the variables, you cannot later call + * {@link #computeLocalUsingPreviousObservations(double[][])}.

+ * + * @param covariance covariance matrix of the system + * @param covFromObservations whether the covariance matrix + * was determined internally from observations or not + * @throws Exception for covariance matrix not matching the expected dimensions, + * being non-square, asymmetric or non-positive definite + */ + public void setCovariance(double[][] cov, boolean covFromObservations) throws Exception { + + if (!covFromObservations) { + // Make sure we're not keeping any observations + observations = null; + } + if (cov.length != dimensions) { + throw new Exception("Supplied covariance matrix does not match initialised number of dimensions"); + } + if (cov.length != cov[0].length) { + throw new Exception("Covariance matrices must be square"); + } + + this.covFromObservations = covFromObservations; + this.covariance = cov; + + } + +} + diff --git a/java/source/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussian.java new file mode 100644 index 0000000..83c75dc --- /dev/null +++ b/java/source/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussian.java @@ -0,0 +1,134 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2017, Joseph T. Lizier, Ipek Oezdemir and Pedro Mediano + * + * 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 . + */ + +package infodynamics.measures.continuous.gaussian; + +import infodynamics.utils.MatrixUtils; + +/** + *

Computes the differential O-information of a given multivariate + * double[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorGaussian}), + * assuming that the probability distribution function for these observations is + * a multivariate Gaussian distribution.

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorCommon}. + *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class OInfoCalculatorGaussian + extends MultiVariateInfoMeasureCalculatorGaussian { + + /** + * Constructor. + */ + public OInfoCalculatorGaussian() { + // Nothing to do + } + + /** + * {@inheritDoc} + * + * @return the average O-info in nats (not bits!) + * @throws Exception if not sufficient data have been provided, or if the + * supplied covariance matrix is invalid. + */ + public double computeAverageLocalOfObservations() throws Exception { + + if (covariance == null) { + throw new Exception("Cannot calculate O-Info without having " + + "a covariance either supplied or computed via setObservations()"); + } + + if (!isComputed) { + double oinfo = (dimensions - 2)*Math.log(MatrixUtils.determinantSymmPosDefMatrix(covariance)); + for (int i = 0; i < dimensions; i++) { + int[] idx = allExcept(i, dimensions); + double[][] marginal_cov = MatrixUtils.selectRowsAndColumns(covariance, idx, idx); + oinfo += Math.log(covariance[i][i]) - Math.log(MatrixUtils.determinantSymmPosDefMatrix(marginal_cov)); + } + // This "0.5" comes from the entropy formula for Gaussians: h = 0.5*logdet(2*pi*e*Sigma) + lastAverage = 0.5*oinfo;; + isComputed = true; + } + + return lastAverage; + } + + /** + * {@inheritDoc} + * + * @return the "time-series" of local O-info values in nats (not bits!) + * for the supplied states. + * @throws Exception if not sufficient data have been provided, or if the + * supplied covariance matrix is invalid. + */ + public double[] computeLocalUsingPreviousObservations(double[][] states) throws Exception { + + if ((means == null) || (covariance == null)) { + throw new Exception("Cannot compute local values without having means " + + "and covariance either supplied or computed via setObservations()"); + } + + EntropyCalculatorMultiVariateGaussian hCalc = new EntropyCalculatorMultiVariateGaussian(); + hCalc.initialise(dimensions); + hCalc.setCovarianceAndMeans(covariance, means); + double[] localValues = MatrixUtils.multiply(hCalc.computeLocalUsingPreviousObservations(states), dimensions - 2); + + for (int i = 0; i < dimensions; i++) { + int[] idx = allExcept(i, dimensions); + + // Local entropy of this variable (i) only + double[][] this_cov = MatrixUtils.selectRowsAndColumns(covariance, i, 1, i, 1); + double[] this_means = MatrixUtils.select(means, i, 1); + double[][] this_state = MatrixUtils.selectColumns(states, i, 1); + + hCalc.initialise(1); + hCalc.setCovarianceAndMeans(this_cov, this_means); + double[] thisLocals = hCalc.computeLocalUsingPreviousObservations(this_state); + + // Local entropy of the rest of the variables (0, ... i-1, i+1, ... D) + double[][] rest_cov = MatrixUtils.selectRowsAndColumns(covariance, idx, idx); + double[] rest_means = MatrixUtils.select(means, idx); + double[][] rest_state = MatrixUtils.selectColumns(states, idx); + + hCalc.initialise(dimensions - 1); + hCalc.setCovarianceAndMeans(rest_cov, rest_means); + double[] restLocals = hCalc.computeLocalUsingPreviousObservations(rest_state); + + localValues = MatrixUtils.add(localValues, MatrixUtils.subtract(thisLocals, restLocals)); + + } + + return localValues; + + } + +} + + diff --git a/java/source/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussian.java new file mode 100644 index 0000000..02c5627 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussian.java @@ -0,0 +1,120 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2017, Joseph T. Lizier, Ipek Oezdemir and Pedro Mediano + * + * 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 . + */ + +package infodynamics.measures.continuous.gaussian; + +import infodynamics.utils.MatrixUtils; + +/** + *

Computes the differential S-information of a given multivariate + * double[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorGaussian}), + * assuming that the probability distribution function for these observations is + * a multivariate Gaussian distribution.

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorCommon}. + *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class SInfoCalculatorGaussian + extends MultiVariateInfoMeasureCalculatorGaussian { + + /** + * Constructor. + */ + public SInfoCalculatorGaussian() { + // Nothing to do + } + + /** + * {@inheritDoc} + * + * @return the average S-info in nats (not bits!) + * @throws Exception if not sufficient data have been provided, or if the + * supplied covariance matrix is invalid. + */ + public double computeAverageLocalOfObservations() throws Exception { + + if (covariance == null) { + throw new Exception("Cannot calculate O-Info without having " + + "a covariance either supplied or computed via setObservations()"); + } + + if (!isComputed) { + + MultiInfoCalculatorGaussian tcCalc = new MultiInfoCalculatorGaussian(); + tcCalc.initialise(dimensions); + tcCalc.setCovariance(covariance); + double tc = tcCalc.computeAverageLocalOfObservations(); + + DualTotalCorrelationCalculatorGaussian dtcCalc = new DualTotalCorrelationCalculatorGaussian(); + dtcCalc.initialise(dimensions); + dtcCalc.setCovariance(covariance); + double dtc = dtcCalc.computeAverageLocalOfObservations(); + + lastAverage = tc + dtc; + isComputed = true; + } + + return lastAverage; + } + + /** + * {@inheritDoc} + * + * @return the "time-series" of local S-info values in nats (not bits!) + * for the supplied states. + * @throws Exception if not sufficient data have been provided, or if the + * supplied covariance matrix is invalid. + */ + public double[] computeLocalUsingPreviousObservations(double[][] states) throws Exception { + + if ((means == null) || (covariance == null)) { + throw new Exception("Cannot compute local values without having means " + + "and covariance either supplied or computed via setObservations()"); + } + + MultiInfoCalculatorGaussian tcCalc = new MultiInfoCalculatorGaussian(); + tcCalc.initialise(dimensions); + tcCalc.setCovarianceAndMeans(covariance, means); + double[] localTC = tcCalc.computeLocalUsingPreviousObservations(states); + + DualTotalCorrelationCalculatorGaussian dtcCalc = new DualTotalCorrelationCalculatorGaussian(); + dtcCalc.initialise(dimensions); + dtcCalc.setCovarianceAndMeans(covariance, means); + double[] localDTC = dtcCalc.computeLocalUsingPreviousObservations(states); + + double[] localValues = MatrixUtils.add(localTC, localDTC); + + return localValues; + + } + +} + + diff --git a/java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java b/java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java new file mode 100644 index 0000000..a5cd702 --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java @@ -0,0 +1,133 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.gaussian; + +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; +import junit.framework.TestCase; + +public class DualTotalCorrelationCalculatorGaussianTester extends TestCase { + + /** + * For two variables, DTC is equal to mutual information. + */ + public void testTwoVariables() throws Exception { + double[][] cov = new double[][] {{1, 0.5}, {0.5, 1}}; + + DualTotalCorrelationCalculatorGaussian dtcCalc = new DualTotalCorrelationCalculatorGaussian(); + dtcCalc.initialise(2); + dtcCalc.setCovariance(cov); + double dtc = dtcCalc.computeAverageLocalOfObservations(); + + MutualInfoCalculatorMultiVariateGaussian miCalc = new MutualInfoCalculatorMultiVariateGaussian(); + miCalc.initialise(1,1); + miCalc.setCovariance(cov, 1); + double mi = miCalc.computeAverageLocalOfObservations(); + + assertEquals(dtc, mi, 1e-6); + } + + /** + * Compare against the direct calculation of DTC as a sum of entropies using the + * entropy calculator. + */ + public void testCompareWithEntropy() throws Exception { + double[][] cov = new double[][] {{1, 0.4, 0.3}, {0.4, 1, 0.2}, {0.3, 0.2, 1}}; + + DualTotalCorrelationCalculatorGaussian dtcCalc = new DualTotalCorrelationCalculatorGaussian(); + dtcCalc.initialise(3); + dtcCalc.setCovariance(cov); + double dtc = dtcCalc.computeAverageLocalOfObservations(); + + // Calculate using an entropy calculator and picking submatrices manually + EntropyCalculatorMultiVariateGaussian hCalc = new EntropyCalculatorMultiVariateGaussian(); + hCalc.initialise(3); + hCalc.setCovariance(cov); + double dtc_hCalc = -2 * hCalc.computeAverageLocalOfObservations(); + + hCalc.initialise(2); + hCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {0,1}, new int[] {0,1})); + dtc_hCalc += hCalc.computeAverageLocalOfObservations(); + hCalc.initialise(2); + hCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {0,2}, new int[] {0,2})); + dtc_hCalc += hCalc.computeAverageLocalOfObservations(); + hCalc.initialise(2); + hCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {1,2}, new int[] {1,2})); + dtc_hCalc += hCalc.computeAverageLocalOfObservations(); + + assertEquals(dtc, dtc_hCalc, 1e-6); + + } + + /** + * Confirm that the local values average correctly back to the average value + */ + public void testLocalsAverageCorrectly() throws Exception { + + int dimensions = 4; + int timeSteps = 1000; + DualTotalCorrelationCalculatorGaussian dtcCalc = new DualTotalCorrelationCalculatorGaussian(); + dtcCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + dtcCalc.setObservations(data); + + double dtc = dtcCalc.computeAverageLocalOfObservations(); + double[] dtcLocal = dtcCalc.computeLocalOfPreviousObservations(); + + System.out.printf("Average was %.5f\n", dtc); + + assertEquals(dtc, MatrixUtils.mean(dtcLocal), 0.00001); + } + + /** + * Confirm that for 2D the local values equal the local MI values + * + */ + public void testLocalsEqualMI() throws Exception { + + int dimensions = 2; + int timeSteps = 100; + DualTotalCorrelationCalculatorGaussian dtcCalc = new DualTotalCorrelationCalculatorGaussian(); + dtcCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + dtcCalc.setObservations(data); + double[] dtcLocal = dtcCalc.computeLocalOfPreviousObservations(); + + MutualInfoCalculatorMultiVariateGaussian miCalc = new MutualInfoCalculatorMultiVariateGaussian(); + miCalc.initialise(1, 1); + miCalc.setObservations(MatrixUtils.selectColumn(data, 0), MatrixUtils.selectColumn(data, 1)); + double[] miLocal = miCalc.computeLocalOfPreviousObservations(); + + for (int t = 0; t < timeSteps; t++) { + assertEquals(dtcLocal[t], miLocal[t], 0.00001); + } + } + +} + diff --git a/java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java b/java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java new file mode 100644 index 0000000..15e3c46 --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java @@ -0,0 +1,150 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.gaussian; + +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; +import junit.framework.TestCase; + +public class OInfoCalculatorGaussianTester extends TestCase { + + /** + * For two variables, O-info is zero + */ + public void testTwoVariables() throws Exception { + double[][] cov = new double[][] {{1, 0.5}, {0.5, 1}}; + + OInfoCalculatorGaussian oCalc = new OInfoCalculatorGaussian(); + oCalc.initialise(2); + oCalc.setCovariance(cov); + double oinfo = oCalc.computeAverageLocalOfObservations(); + + assertEquals(oinfo, 0, 1e-6); + } + + + /** + * For factorisable pairwise interactions, O-info is zero + */ + public void testPairwise() throws Exception { + double[][] cov = new double[][] {{ 1, 0.5, 0, 0}, + {0.5, 1, 0, 0}, + { 0, 0, 1, 0.5}, + { 0, 0, 0.5, 1}}; + + OInfoCalculatorGaussian oCalc = new OInfoCalculatorGaussian(); + oCalc.initialise(4); + oCalc.setCovariance(cov); + double oinfo = oCalc.computeAverageLocalOfObservations(); + + assertEquals(oinfo, 0, 1e-6); + } + + /** + * Compare against the direct calculation of O-info as a sum of entropies using the + * entropy calculator. + */ + public void testCompareWithEntropy() throws Exception { + double[][] cov = new double[][] {{1, 0.4, 0.3}, {0.4, 1, 0.2}, {0.3, 0.2, 1}}; + + OInfoCalculatorGaussian oCalc = new OInfoCalculatorGaussian(); + oCalc.initialise(3); + oCalc.setCovariance(cov); + double oinfo = oCalc.computeAverageLocalOfObservations(); + + // Calculate using an entropy calculator and picking submatrices manually + EntropyCalculatorMultiVariateGaussian hCalc = new EntropyCalculatorMultiVariateGaussian(); + hCalc.initialise(3); + hCalc.setCovariance(cov); + double oinfo_hCalc = hCalc.computeAverageLocalOfObservations(); + + hCalc.initialise(2); + hCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {0,1}, new int[] {0,1})); + oinfo_hCalc -= hCalc.computeAverageLocalOfObservations(); + hCalc.initialise(2); + hCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {0,2}, new int[] {0,2})); + oinfo_hCalc -= hCalc.computeAverageLocalOfObservations(); + hCalc.initialise(2); + hCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {1,2}, new int[] {1,2})); + oinfo_hCalc -= hCalc.computeAverageLocalOfObservations(); + + hCalc.initialise(1); + hCalc.setCovariance(new double[][] {{cov[0][0]}}); + oinfo_hCalc += hCalc.computeAverageLocalOfObservations(); + hCalc.initialise(1); + hCalc.setCovariance(new double[][] {{cov[1][1]}}); + oinfo_hCalc += hCalc.computeAverageLocalOfObservations(); + hCalc.initialise(1); + hCalc.setCovariance(new double[][] {{cov[2][2]}}); + oinfo_hCalc += hCalc.computeAverageLocalOfObservations(); + + assertEquals(oinfo, oinfo_hCalc, 1e-6); + + } + + /** + * Confirm that the local values average correctly back to the average value + */ + public void testLocalsAverageCorrectly() throws Exception { + + int dimensions = 4; + int timeSteps = 1000; + OInfoCalculatorGaussian oCalc = new OInfoCalculatorGaussian(); + oCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + oCalc.setObservations(data); + + double oinfo = oCalc.computeAverageLocalOfObservations(); + double[] oLocal = oCalc.computeLocalOfPreviousObservations(); + + System.out.printf("Average was %.5f\n", oinfo); + + assertEquals(oinfo, MatrixUtils.mean(oLocal), 0.00001); + } + + /** + * Confirm that for 2D all local values equal zero + */ + public void testLocalsEqualZero() throws Exception { + + int dimensions = 2; + int timeSteps = 100; + OInfoCalculatorGaussian oCalc = new OInfoCalculatorGaussian(); + oCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + oCalc.setObservations(data); + double[] oLocal = oCalc.computeLocalOfPreviousObservations(); + + for (int t = 0; t < timeSteps; t++) { + assertEquals(oLocal[t], 0, 0.00001); + } + } + +} + diff --git a/java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java b/java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java new file mode 100644 index 0000000..bc33ad7 --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java @@ -0,0 +1,129 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.gaussian; + +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; +import junit.framework.TestCase; + +public class SInfoCalculatorGaussianTester extends TestCase { + + /** + * For two variables, S-info is twice the MI between them + */ + public void testTwoVariables() throws Exception { + double[][] cov = new double[][] {{1, 0.5}, {0.5, 1}}; + + SInfoCalculatorGaussian sCalc = new SInfoCalculatorGaussian(); + sCalc.initialise(2); + sCalc.setCovariance(cov); + double sinfo = sCalc.computeAverageLocalOfObservations(); + + MutualInfoCalculatorMultiVariateGaussian miCalc = new MutualInfoCalculatorMultiVariateGaussian(); + miCalc.initialise(1,1); + miCalc.setCovariance(cov, false); + double mi = miCalc.computeAverageLocalOfObservations(); + + assertEquals(sinfo, 2*mi, 1e-6); + } + + + /** + * Compare against the direct calculation of S-info as a sum of mutual informations. + */ + public void testCompareWithEntropy() throws Exception { + double[][] cov = new double[][] {{1, 0.4, 0.3}, {0.4, 1, 0.2}, {0.3, 0.2, 1}}; + + SInfoCalculatorGaussian sCalc = new SInfoCalculatorGaussian(); + sCalc.initialise(3); + sCalc.setCovariance(cov); + double sinfo = sCalc.computeAverageLocalOfObservations(); + + // Calculate using a mutual info calculator and picking submatrices manually + MutualInfoCalculatorMultiVariateGaussian miCalc = new MutualInfoCalculatorMultiVariateGaussian(); + miCalc.initialise(2,1); + miCalc.setCovariance(cov, false); + double sinfo_miCalc = miCalc.computeAverageLocalOfObservations(); + + miCalc.initialise(2,1); + miCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {2,0,1}, new int[] {2,0,1}), false); + sinfo_miCalc += miCalc.computeAverageLocalOfObservations(); + miCalc.initialise(2,1); + miCalc.setCovariance(MatrixUtils.selectRowsAndColumns(cov, new int[] {1,2,0}, new int[] {1,2,0}), false); + sinfo_miCalc += miCalc.computeAverageLocalOfObservations(); + + assertEquals(sinfo, sinfo_miCalc, 1e-6); + + } + + /** + * Confirm that the local values average correctly back to the average value + */ + public void testLocalsAverageCorrectly() throws Exception { + + int dimensions = 4; + int timeSteps = 1000; + SInfoCalculatorGaussian sCalc = new SInfoCalculatorGaussian(); + sCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + sCalc.setObservations(data); + + double sinfo = sCalc.computeAverageLocalOfObservations(); + double[] sLocal = sCalc.computeLocalOfPreviousObservations(); + + System.out.printf("Average was %.5f\n", sinfo); + + assertEquals(sinfo, MatrixUtils.mean(sLocal), 0.00001); + } + + /** + * Confirm that for 2D all local values equal zero + */ + public void testLocalsEqualMI() throws Exception { + + int dimensions = 2; + int timeSteps = 100; + SInfoCalculatorGaussian sCalc = new SInfoCalculatorGaussian(); + sCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + sCalc.setObservations(data); + double[] sLocal = sCalc.computeLocalOfPreviousObservations(); + + MutualInfoCalculatorMultiVariateGaussian miCalc = new MutualInfoCalculatorMultiVariateGaussian(); + miCalc.initialise(1, 1); + miCalc.setObservations(MatrixUtils.selectColumn(data, 0), MatrixUtils.selectColumn(data, 1)); + double[] miLocal = miCalc.computeLocalOfPreviousObservations(); + + for (int t = 0; t < timeSteps; t++) { + assertEquals(sLocal[t], 2*miLocal[t], 0.00001); + } + } + +} + From b3abd382aa63246c24909eb5bc975d918ee3563e Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Sun, 24 Jan 2021 22:10:14 +0000 Subject: [PATCH 5/8] Added abstract parent class for multivariate IT measures on discrete data. --- ...iVariateInfoMeasureCalculatorDiscrete.java | 279 ++++++++++++++++++ 1 file changed, 279 insertions(+) create mode 100755 java/source/infodynamics/measures/discrete/MultiVariateInfoMeasureCalculatorDiscrete.java diff --git a/java/source/infodynamics/measures/discrete/MultiVariateInfoMeasureCalculatorDiscrete.java b/java/source/infodynamics/measures/discrete/MultiVariateInfoMeasureCalculatorDiscrete.java new file mode 100755 index 0000000..1d3f826 --- /dev/null +++ b/java/source/infodynamics/measures/discrete/MultiVariateInfoMeasureCalculatorDiscrete.java @@ -0,0 +1,279 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.discrete; + +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +/** + * Implements a base class with common functionality for child class + * implementations of multivariate information measures. + * + *

Multivariate information measures are functionals of probability + * distributions over R^n, and typical examples include multi-information + * (a.k.a. total correlation), dual total correlation, O-information, and connected + * information.

+ * + *

Usage of child classes is intended to follow this paradigm:

+ *
    + *
  1. Construct the calculator;
    2. + *
  2. Initialise the calculator using {@link #initialise()};
    3. + *
  3. Provide the observations/samples for the calculator + * to set up the PDFs, using one or more calls to + * sets of {@link #addObservations(int[][], int[])} methods, then
    4. + *
  4. Compute the required quantities, being one or more of: + *
      + *
    • the average measure: {@link #computeAverageLocalOfObservations()};
      • + *
    + *
    5. + *
  5. + * Return to step 2 to re-use the calculator on a new data set. + *
    6. + *
+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public abstract class MultiVariateInfoMeasureCalculatorDiscrete + extends InfoMeasureCalculatorDiscrete { + + /** + * Count of occurrences of each joint state in the provided observations. + */ + protected int[] jointCount = null; + + /** + * Count of occurrences of each state of each variable in the provided + * observations. + * + * For a given variable v and state i, + * smallMarginalCounts[v][i] counts how many times variable + * v was observed in state i, + */ + protected int[][] smallMarginalCounts = null; // marginalCounts[marginalIndex][state] + + /** + * Count of occurrences of each state of each (D-1)-dimensional marginal in + * the provided observations. + * + * For a given variable v and state i, + * bigMarginalCounts[v][i] counts how many times the _rest_ of + * the system, excluding variable v, was observed in state i, + */ + protected int[][] bigMarginalCounts = null; + + /** + * Number of variables in the system. + */ + protected int numVars; + + /** + * Number of possible states of the whole system. + */ + protected int jointStates; + + /** + * Whether the first local value has been checked. (Used to initialise some variables + * related to computation of local values. + */ + protected boolean checkedFirst = false; + + /** + * Abstract constructor (to be called by child classes). + * + * @param base number of symbols for each variable. + * E.g. binary variables are in base-2. + * @param numVars numbers of joint variables that the measure + * will be computed over. + */ + protected MultiVariateInfoMeasureCalculatorDiscrete(int base, int numVars) { + super(base); + this.numVars = numVars; + jointStates = MathsUtils.power(base, numVars); + try { + jointCount = new int[jointStates]; + smallMarginalCounts = new int[numVars][base]; + bigMarginalCounts = new int[numVars][jointStates]; + } catch (OutOfMemoryError e) { + // Allow any Exceptions to be thrown, but catch and wrap + // Error as a RuntimeException + throw new RuntimeException("Requested memory for the base " + + base + " with " + numVars + + " variables is too large for the JVM at this time", e); + } + } + + @Override + public void initialise(){ + super.initialise(); + MatrixUtils.fill(jointCount, 0); + MatrixUtils.fill(smallMarginalCounts, 0); + MatrixUtils.fill(bigMarginalCounts, 0); + } + + /** + * Given multiple time samples of a homogeneous array of variables (states), + * add the observations of all sets of numVars of these + * Do this for every time point + * + * @param states 2D array of values of an array of variables + * at many observations (first index is time, second is variable index) + */ + public void addObservations(int[][] states) throws Exception { + int[] jointStates = MatrixUtils.computeCombinedValues(states, base); + for (int t = 0; t < states.length; t++) { + for (int i = 0; i < numVars; i++) { + // Extract values of the 1D and the (N-1)D marginals + int thisValue = states[t][i]; + int bigMarginalState = computeBigMarginalState(jointStates[t], i, thisValue); + + // Update counts + bigMarginalCounts[i][bigMarginalState]++; + smallMarginalCounts[i][thisValue]++; + } + jointCount[jointStates[t]]++; + observations++; + } + } + + @Override + public double computeAverageLocalOfObservations() { + + int[] jointTuple = new int[numVars]; + checkedFirst = false; + try { + average = computeForGivenTupleFromVarIndex(jointTuple, 0); + } catch (Exception e) { + System.out.println("Something went wrong during the calculation."); + average = -1; + } + + return average; + } + + /** + * Private utility to compute the contribution to the measure for all tuples + * starting with tuple[0..(fromIndex-1)]. + * + * @param tuple + * @param fromIndex + * @return + */ + public double computeForGivenTupleFromVarIndex(int[] tuple, int fromIndex) throws Exception { + double miCont = 0; + if (fromIndex == numVars) { + // The whole tuple is filled in, so compute the contribution to the MI from this tuple + int jointValue = MatrixUtils.computeCombinedValues(new int[][] {tuple}, base)[0]; + + if (jointCount[jointValue] == 0) { + // This joint state does not occur, so it makes no contribution here + return 0; + } + + double jointProb = (double) jointCount[jointValue] / (double) observations; + double localValue = computeLocalValueForTuple(tuple, jointValue); + miCont = jointProb * localValue; + + } else { + // Fill out the next part of the tuple and make the recursive calls + for (int v = 0; v < base; v++) { + tuple[fromIndex] = v; + miCont += computeForGivenTupleFromVarIndex(tuple, fromIndex + 1); + } + } + return miCont; + } + + /** + * Shortcut method to initialise the calculator, add observations and compute + * the average measure in one line. + * + * @param state series of multivariate observations + * (first index is time or observation index, second is variable number) + */ + public double compute(int[][] states) throws Exception { + initialise(); + addObservations(states); + return computeAverageLocalOfObservations(); + } + + /** + * Internal method to update maximum and minimum values of local information + * measures. + * + * @param localValue instance of computed local information measure + */ + protected void checkLocals(double localValue) { + if (!checkedFirst) { + max = localValue; + min = localValue; + checkedFirst = true; + } else { + if (localValue > max) { + max = localValue; + } + if (localValue < min) { + min = localValue; + } + } + } + + /** + * Method to be implemented by all child classes to compute the local value + * of the measure for a given tuple. + * + * @param tuple state of the system at a given time (index is variable number) + * @param jointValue int representing the state of the system + */ + protected abstract double computeLocalValueForTuple(int[] tuple, int jointValue) + throws Exception; + + /** + * Method to be implemented by all child classes to compute the local value + * of the measure for a given tuple. + * + * @param tuple state of the system at a given time (index is variable number) + */ + protected double computeLocalValueForTuple(int[] tuple) throws Exception { + int jointValue = MatrixUtils.computeCombinedValues(new int[][] {tuple}, base)[0]; + return computeLocalValueForTuple(tuple, jointValue); + } + + /** + * Small utility function to compute the state of the system excluding one variable. + * + * @param jointState state of the full system + * @param varIdx index of the variable to be excluded + * @param varValue value of the variable in question in the system state + */ + protected int computeBigMarginalState(int jointState, int varIdx, int varValue) { + int bigMarginalState = jointState - varValue*MathsUtils.power(base, numVars - varIdx - 1); + return bigMarginalState; + } + +} + From a58a01fbfd29c5978147c14596a92ed8218cf51a Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Sun, 24 Jan 2021 22:13:15 +0000 Subject: [PATCH 6/8] Added discrete implementation of various multivariate IT measures and unit tests. --- ...ualTotalCorrelationCalculatorDiscrete.java | 84 ++++++++++++++ .../discrete/OInfoCalculatorDiscrete.java | 85 ++++++++++++++ .../discrete/SInfoCalculatorDiscrete.java | 95 +++++++++++++++ .../discrete/DualTotalCorrelationTester.java | 108 ++++++++++++++++++ .../measures/discrete/OInfoTester.java | 87 ++++++++++++++ .../measures/discrete/SInfoTester.java | 79 +++++++++++++ 6 files changed, 538 insertions(+) create mode 100755 java/source/infodynamics/measures/discrete/DualTotalCorrelationCalculatorDiscrete.java create mode 100755 java/source/infodynamics/measures/discrete/OInfoCalculatorDiscrete.java create mode 100755 java/source/infodynamics/measures/discrete/SInfoCalculatorDiscrete.java create mode 100755 java/unittests/infodynamics/measures/discrete/DualTotalCorrelationTester.java create mode 100755 java/unittests/infodynamics/measures/discrete/OInfoTester.java create mode 100755 java/unittests/infodynamics/measures/discrete/SInfoTester.java diff --git a/java/source/infodynamics/measures/discrete/DualTotalCorrelationCalculatorDiscrete.java b/java/source/infodynamics/measures/discrete/DualTotalCorrelationCalculatorDiscrete.java new file mode 100755 index 0000000..8bcbc93 --- /dev/null +++ b/java/source/infodynamics/measures/discrete/DualTotalCorrelationCalculatorDiscrete.java @@ -0,0 +1,84 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.discrete; + +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +/** + *

Computes the dual total correlation (DTC) of a given multivariate + * int[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorDiscrete}).

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorDiscrete}. + *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class DualTotalCorrelationCalculatorDiscrete + extends MultiVariateInfoMeasureCalculatorDiscrete { + + /** + * Construct an instance. + * + * @param base number of symbols for each variable. + * E.g. binary variables are in base-2. + * @param numVars numbers of joint variables that DTC + * will be computed over. + */ + public DualTotalCorrelationCalculatorDiscrete(int base, int numVars) { + super(base, numVars); + } + + protected double computeLocalValueForTuple(int[] tuple, int jointValue) { + + if (jointCount[jointValue] == 0) { + // This joint state does not occur, so it makes no contribution here + return 0; + } + + double jointProb = (double) jointCount[jointValue] / (double) observations; + double logValue = (numVars - 1) * Math.log(jointProb); + + for (int i = 0; i < numVars; i++) { + int marginalState = computeBigMarginalState(jointValue, i, tuple[i]); + double marginalProb = (double) bigMarginalCounts[i][marginalState] / (double) observations; + logValue -= Math.log(marginalProb); + } + + double localValue = logValue / log_2; + + if (jointProb > 0.0) { + checkLocals(localValue); + } + + return localValue; + } + +} + diff --git a/java/source/infodynamics/measures/discrete/OInfoCalculatorDiscrete.java b/java/source/infodynamics/measures/discrete/OInfoCalculatorDiscrete.java new file mode 100755 index 0000000..118df47 --- /dev/null +++ b/java/source/infodynamics/measures/discrete/OInfoCalculatorDiscrete.java @@ -0,0 +1,85 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.discrete; + +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +/** + *

Computes the O-information of a given multivariate + * int[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorDiscrete}).

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorDiscrete}. + *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class OInfoCalculatorDiscrete + extends MultiVariateInfoMeasureCalculatorDiscrete { + + /** + * Construct an instance. + * + * @param base number of symbols for each variable. + * E.g. binary variables are in base-2. + * @param numVars numbers of joint variables that DTC + * will be computed over. + */ + public OInfoCalculatorDiscrete(int base, int numVars) { + super(base, numVars); + } + + protected double computeLocalValueForTuple(int[] tuple, int jointValue) { + + if (jointCount[jointValue] == 0) { + // This joint state does not occur, so it makes no contribution here + return 0; + } + + double jointProb = (double) jointCount[jointValue] / (double) observations; + double logValue = (2 - numVars) * Math.log(jointProb); + + for (int i = 0; i < numVars; i++) { + int bigMarginalState = computeBigMarginalState(jointValue, i, tuple[i]); + double bigMarginalProb = (double) bigMarginalCounts[i][bigMarginalState] / (double) observations; + double smallMarginalProb = (double) smallMarginalCounts[i][tuple[i]] / (double) observations; + logValue += Math.log(bigMarginalProb) - Math.log(smallMarginalProb); + } + + double localValue = logValue / log_2; + + if (jointProb > 0.0) { + checkLocals(localValue); + } + + return localValue; + } + +} + diff --git a/java/source/infodynamics/measures/discrete/SInfoCalculatorDiscrete.java b/java/source/infodynamics/measures/discrete/SInfoCalculatorDiscrete.java new file mode 100755 index 0000000..150118d --- /dev/null +++ b/java/source/infodynamics/measures/discrete/SInfoCalculatorDiscrete.java @@ -0,0 +1,95 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.discrete; + +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +/** + *

Computes the S-information of a given multivariate + * int[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorDiscrete}).

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorDiscrete}. + *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class SInfoCalculatorDiscrete + extends MultiVariateInfoMeasureCalculatorDiscrete { + + /** + * Construct an instance. + * + * @param base number of symbols for each variable. + * E.g. binary variables are in base-2. + * @param numVars numbers of joint variables that DTC + * will be computed over. + */ + public SInfoCalculatorDiscrete(int base, int numVars) { + super(base, numVars); + } + + protected double computeLocalValueForTuple(int[] tuple, int jointValue) { + + if (jointCount[jointValue] == 0) { + // This joint state does not occur, so it makes no contribution here + return 0; + } + + double jointProb = (double) jointCount[jointValue] / (double) observations; + + // Local TC value + double localTC = Math.log(jointProb); + for (int i = 0; i < numVars; i++) { + int marginalState = tuple[i]; + double marginalProb = (double) smallMarginalCounts[i][marginalState] / (double) observations; + localTC -= Math.log(marginalProb); + } + + // Local DTC value + double localDTC = (numVars - 1) * Math.log(jointProb); + for (int i = 0; i < numVars; i++) { + int marginalState = computeBigMarginalState(jointValue, i, tuple[i]); + double marginalProb = (double) bigMarginalCounts[i][marginalState] / (double) observations; + localDTC -= Math.log(marginalProb); + } + + // Combine local TC and DTC into S-info + double logValue = localTC + localDTC; + double localValue = logValue / log_2; + + if (jointProb > 0.0) { + checkLocals(localValue); + } + + return localValue; + } + +} + diff --git a/java/unittests/infodynamics/measures/discrete/DualTotalCorrelationTester.java b/java/unittests/infodynamics/measures/discrete/DualTotalCorrelationTester.java new file mode 100755 index 0000000..5557e6b --- /dev/null +++ b/java/unittests/infodynamics/measures/discrete/DualTotalCorrelationTester.java @@ -0,0 +1,108 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.discrete; + +import infodynamics.utils.RandomGenerator; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.MathsUtils; + +import java.util.Arrays; +import junit.framework.TestCase; + + +public class DualTotalCorrelationTester extends TestCase { + + public void testIndependent() throws Exception { + DualTotalCorrelationCalculatorDiscrete dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, 3); + double dtc = dtcCalc.compute(new int[][] {{0,0,0},{0,0,1},{0,1,0},{0,1,1},{1,0,0},{1,0,1},{1,1,0},{1,1,1}}); + assertEquals(0.0, dtc, 0.000001); + } + + public void testXor() throws Exception { + // 3 variables + DualTotalCorrelationCalculatorDiscrete dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, 3); + double dtc = dtcCalc.compute(new int[][] {{0,0,1},{0,1,0},{1,0,0},{1,1,1}}); + assertEquals(2.0, dtc, 0.000001); + + // 4 variables + dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, 4); + dtc = dtcCalc.compute(new int[][] {{0,0,0,1},{0,0,1,0},{0,1,0,0},{1,0,0,0},{1,1,1,0},{1,1,0,1},{1,0,1,1},{0,1,1,1}}); + assertEquals(3.0, dtc, 0.000001); + } + + public void testCopy() throws Exception { + // 3 variables + DualTotalCorrelationCalculatorDiscrete dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, 3); + double dtc = dtcCalc.compute(new int[][] {{0,0,0},{1,1,1}}); + assertEquals(1.0, dtc, 0.000001); + + // 4 variables + dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, 4); + dtc = dtcCalc.compute(new int[][] {{0,0,0,0},{1,1,1,1}}); + assertEquals(1.0, dtc, 0.000001); + } + + public void testCompareEntropy() throws Exception { + // Generate random data and check that it matches the explicit computation + // using entropy calculators + RandomGenerator rg = new RandomGenerator(); + int D = 4; + int[][] data = rg.generateRandomInts(10, D, 2); + + // DTC calculator + DualTotalCorrelationCalculatorDiscrete dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, D); + double dtc_direct = dtcCalc.compute(data); + + // Entropy calculators + EntropyCalculatorDiscrete hCalc = new EntropyCalculatorDiscrete(MathsUtils.power(2, D)); + hCalc.initialise(); + hCalc.addObservations(MatrixUtils.computeCombinedValues(data, 2)); + double dtc_test = (1 - D) * hCalc.computeAverageLocalOfObservations(); + + hCalc = new EntropyCalculatorDiscrete(MathsUtils.power(2, D-1)); + for (int i = 0; i < 4; i++) { + hCalc.initialise(); + hCalc.addObservations(MatrixUtils.computeCombinedValues(MatrixUtils.selectColumns(data, allExcept(i, D)), 2)); + dtc_test += hCalc.computeAverageLocalOfObservations(); + } + + assertEquals(dtc_direct, dtc_test, 0.000001); + + } + + protected int[] allExcept(int idx, int N) { + boolean[] v = new boolean[N]; + Arrays.fill(v, true); + v[idx] = false; + + int[] v2 = new int[N - 1]; + int counter = 0; + for (int i = 0; i < N; i++) { + if (v[i]) { + v2[counter] = i; + counter++; + } + } + + return v2; + } + + +} + diff --git a/java/unittests/infodynamics/measures/discrete/OInfoTester.java b/java/unittests/infodynamics/measures/discrete/OInfoTester.java new file mode 100755 index 0000000..9d9c47b --- /dev/null +++ b/java/unittests/infodynamics/measures/discrete/OInfoTester.java @@ -0,0 +1,87 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.discrete; + +import infodynamics.utils.RandomGenerator; + +import junit.framework.TestCase; + +public class OInfoTester extends TestCase { + + public void testIndependent() throws Exception { + OInfoCalculatorDiscrete oCalc = new OInfoCalculatorDiscrete(2, 3); + double oinfo = oCalc.compute(new int[][] {{0,0,0},{0,0,1},{0,1,0},{0,1,1},{1,0,0},{1,0,1},{1,1,0},{1,1,1}}); + assertEquals(0.0, oinfo, 0.000001); + } + + public void testXor() throws Exception { + // 3 variables + OInfoCalculatorDiscrete oCalc = new OInfoCalculatorDiscrete(2, 3); + double oinfo = oCalc.compute(new int[][] {{0,0,1},{0,1,0},{1,0,0},{1,1,1}}); + assertEquals(-1.0, oinfo, 0.000001); + + // 4 variables + oCalc = new OInfoCalculatorDiscrete(2, 4); + oinfo = oCalc.compute(new int[][] {{0,0,0,1},{0,0,1,0},{0,1,0,0},{1,0,0,0},{1,1,1,0},{1,1,0,1},{1,0,1,1},{0,1,1,1}}); + assertEquals(-2.0, oinfo, 0.000001); + } + + public void testCopy() throws Exception { + // 3 variables + OInfoCalculatorDiscrete oCalc = new OInfoCalculatorDiscrete(2, 3); + double oinfo = oCalc.compute(new int[][] {{0,0,0},{1,1,1}}); + assertEquals(1.0, oinfo, 0.000001); + + // 4 variables + oCalc = new OInfoCalculatorDiscrete(2, 4); + oinfo = oCalc.compute(new int[][] {{0,0,0,0},{1,1,1,1}}); + assertEquals(2.0, oinfo, 0.000001); + } + + public void testPairwise() throws Exception { + // Variables 0 and 1 are correlated and independent from 2 and 3, that are + // also correlated + OInfoCalculatorDiscrete oCalc = new OInfoCalculatorDiscrete(2, 4); + double oinfo = oCalc.compute(new int[][] {{0,0,0,0},{0,0,1,1},{1,1,0,0},{1,1,1,1}}); + assertEquals(0.0, oinfo, 0.000001); + } + + public void testCompareTCAndDTC() throws Exception { + // Generate random data and check that it matches the explicit computation + // using TC and DTC calculators + RandomGenerator rg = new RandomGenerator(); + int[][] data = rg.generateRandomInts(10, 4, 2); + + // O-info calculator + OInfoCalculatorDiscrete oCalc = new OInfoCalculatorDiscrete(2, 4); + double oinfo_direct = oCalc.compute(data); + + // TC and DTC calculators + MultiInformationCalculatorDiscrete tcCalc = new MultiInformationCalculatorDiscrete(2, 4); + DualTotalCorrelationCalculatorDiscrete dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, 4); + tcCalc.initialise(); + tcCalc.addObservations(data); + double oinfo_test = tcCalc.computeAverageLocalOfObservations() - dtcCalc.compute(data); + + assertEquals(oinfo_direct, oinfo_test, 0.000001); + + } + +} + diff --git a/java/unittests/infodynamics/measures/discrete/SInfoTester.java b/java/unittests/infodynamics/measures/discrete/SInfoTester.java new file mode 100755 index 0000000..959c18f --- /dev/null +++ b/java/unittests/infodynamics/measures/discrete/SInfoTester.java @@ -0,0 +1,79 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.discrete; + +import infodynamics.utils.RandomGenerator; + +import junit.framework.TestCase; + +public class SInfoTester extends TestCase { + + public void testIndependent() throws Exception { + SInfoCalculatorDiscrete sCalc = new SInfoCalculatorDiscrete(2, 3); + double sinfo = sCalc.compute(new int[][] {{0,0,0},{0,0,1},{0,1,0},{0,1,1},{1,0,0},{1,0,1},{1,1,0},{1,1,1}}); + assertEquals(0.0, sinfo, 0.000001); + } + + public void testXor() throws Exception { + // 3 variables + SInfoCalculatorDiscrete sCalc = new SInfoCalculatorDiscrete(2, 3); + double sinfo = sCalc.compute(new int[][] {{0,0,1},{0,1,0},{1,0,0},{1,1,1}}); + assertEquals(3.0, sinfo, 0.000001); + + // 4 variables + sCalc = new SInfoCalculatorDiscrete(2, 4); + sinfo = sCalc.compute(new int[][] {{0,0,0,1},{0,0,1,0},{0,1,0,0},{1,0,0,0},{1,1,1,0},{1,1,0,1},{1,0,1,1},{0,1,1,1}}); + assertEquals(4.0, sinfo, 0.000001); + } + + public void testCopy() throws Exception { + // 3 variables + SInfoCalculatorDiscrete sCalc = new SInfoCalculatorDiscrete(2, 3); + double sinfo = sCalc.compute(new int[][] {{0,0,0},{1,1,1}}); + assertEquals(3.0, sinfo, 0.000001); + + // 4 variables + sCalc = new SInfoCalculatorDiscrete(2, 4); + sinfo = sCalc.compute(new int[][] {{0,0,0,0},{1,1,1,1}}); + assertEquals(4.0, sinfo, 0.000001); + } + + public void testCompareTCAndDTC() throws Exception { + // Generate random data and check that it matches the explicit computation + // using TC and DTC calculators + RandomGenerator rg = new RandomGenerator(); + int[][] data = rg.generateRandomInts(10, 4, 2); + + // O-info calculator + SInfoCalculatorDiscrete sCalc = new SInfoCalculatorDiscrete(2, 4); + double sinfo_direct = sCalc.compute(data); + + // TC and DTC calculators + MultiInformationCalculatorDiscrete tcCalc = new MultiInformationCalculatorDiscrete(2, 4); + DualTotalCorrelationCalculatorDiscrete dtcCalc = new DualTotalCorrelationCalculatorDiscrete(2, 4); + tcCalc.initialise(); + tcCalc.addObservations(data); + double sinfo_test = tcCalc.computeAverageLocalOfObservations() + dtcCalc.compute(data); + + assertEquals(sinfo_direct, sinfo_test, 0.000001); + + } + +} + From bc8c233e68939842ce96100efe6beeb6fdb51a9f Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Mon, 25 Jan 2021 12:41:01 +0000 Subject: [PATCH 7/8] Added KSG implementation of various multivariate IT measures and unit tests. --- ...DualTotalCorrelationCalculatorKraskov.java | 132 +++++ ...tiVariateInfoMeasureCalculatorKraskov.java | 538 ++++++++++++++++++ .../kraskov/OInfoCalculatorKraskov.java | 157 +++++ .../kraskov/SInfoCalculatorKraskov.java | 146 +++++ ...talCorrelationCalculatorKraskovTester.java | 143 +++++ .../kraskov/OInfoCalculatorKraskovTester.java | 62 ++ .../kraskov/SInfoCalculatorKraskovTester.java | 62 ++ 7 files changed, 1240 insertions(+) create mode 100755 java/source/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskov.java create mode 100755 java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java create mode 100755 java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java create mode 100755 java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java create mode 100644 java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java create mode 100644 java/unittests/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskovTester.java create mode 100644 java/unittests/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskovTester.java diff --git a/java/source/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskov.java new file mode 100755 index 0000000..97eb5fe --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskov.java @@ -0,0 +1,132 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.kraskov; + +import infodynamics.utils.EuclideanUtils; +import infodynamics.utils.NeighbourNodeData; +import infodynamics.utils.KdTree; +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +import java.util.PriorityQueue; +import java.util.Calendar; +import java.util.Random; + +/** + *

Computes the differential dual total correlation of a given multivariate + * set of observations, using Kraskov-Stoegbauer-Grassberger (KSG) estimation + * (see Kraskov et al., below).

+ * + *

Usage is as per the paradigm outlined for + * {@link MultiVariateInfoMeasureCalculatorCommon}.

+ * + *

Finally, note that {@link Cloneable} is implemented allowing clone() + * to produce only an automatic shallow copy, which is fine + * for the statistical significance calculation it is intended for + * (none of the array + * data will be changed there). + *

+ * + *

References:
+ *

+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class DualTotalCorrelationCalculatorKraskov + extends MultiVariateInfoMeasureCalculatorKraskov + implements Cloneable { // See comments on clonability above + + + protected double[] partialComputeFromObservations( + int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception { + + double startTime = Calendar.getInstance().getTimeInMillis(); + + double[] localMi = null; + if (returnLocals) { + localMi = new double[numTimePoints]; + } + + // Constants: + double dimensionsMinus1TimesDiGammaN = (double) (dimensions - 1) * digammaN; + + // Count the average number of points within eps_x for each marginal x of each point + double totalSumF = 0.0; + + for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) { + // Compute eps for this time step by + // finding the kth closest neighbour for point t: + PriorityQueue nnPQ = + kdTreeJoint.findKNearestNeighbours(k, t, dynCorrExclTime); + // First element in the PQ is the kth NN, + // and epsilon = kthNnData.distance + NeighbourNodeData kthNnData = nnPQ.poll(); + + // Distance to kth neighbour in joint space + double eps = kthNnData.distance; + + double sumF = 0.0; + + sumF += (digammaK - digammaN); + + for (int d = 0; d < dimensions; d++) { + int n = rangeSearchersInBigMarginals[d].countPointsStrictlyWithinR( + t, eps, dynCorrExclTime); + sumF -= (MathsUtils.digamma(n + 1) - digammaN)/(dimensions - 1); + } + + sumF *= (dimensions - 1); + + totalSumF += sumF; + + if (returnLocals) { + localMi[t-startTimePoint] = sumF; + } + } + + if (debug) { + Calendar rightNow2 = Calendar.getInstance(); + long endTime = rightNow2.getTimeInMillis(); + System.out.println("Subset " + startTimePoint + ":" + + (startTimePoint + numTimePoints) + " Calculation time: " + + ((endTime - startTime)/1000.0) + " sec" ); + } + + // Select what to return: + if (returnLocals) { + return localMi; + } else { + double[] returnArray = new double[] {totalSumF/((double) totalObservations)}; + return returnArray; + } + + } + +} + diff --git a/java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java new file mode 100755 index 0000000..df8bcf4 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java @@ -0,0 +1,538 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.kraskov; + +import infodynamics.measures.continuous.MultiVariateInfoMeasureCalculatorCommon; +import infodynamics.utils.EuclideanUtils; +import infodynamics.utils.NeighbourNodeData; +import infodynamics.utils.KdTree; +import infodynamics.utils.UnivariateNearestNeighbourSearcher; +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +import java.util.PriorityQueue; +import java.util.Calendar; +import java.util.Random; + +/** + *

Base class with common functionality for child class implementations of + * multivariate information measures on a given multivariate + * double[][] set of + * observations (extending {@link MultiVariateInfoMeasureCalculatorCommon}), + * using Kraskov-Stoegbauer-Grassberger (KSG) estimation + * (see Kraskov et al., below).

+ * + *

Usage is as per the paradigm outlined for {@link MultiVariateInfoMeasureCalculatorCommon}, + * with: + *

    + *
  • For constructors see the child classes.
    • + *
  • Further properties are defined in {@link #setProperty(String, String)}.
    • + *
  • Computed values are in nats, not bits!
    • + *
+ *

+ * + *

Finally, note that {@link Cloneable} is implemented allowing clone() + * to produce only an automatic shallow copy, which is fine + * for the statistical significance calculation it is intended for + * (none of the array data will be changed there). + *

+ * + *

References:
+ *

+ * @author Pedro A.M. Mediano (email, + * www) + */ +public abstract class MultiVariateInfoMeasureCalculatorKraskov + extends MultiVariateInfoMeasureCalculatorCommon + implements Cloneable { // See comments on clonability above + + /** + * we compute distances to the kth nearest neighbour + */ + protected int k = 4; + + /** + * The norm type in use (see {@link #PROP_NORM_TYPE}) + */ + protected int normType = EuclideanUtils.NORM_MAX_NORM; + + /** + * Property name for the number of K nearest neighbours used in + * the KSG algorithm (default 4). + */ + public final static String PROP_K = "k"; + /** + * Property name for what type of norm to use between data points + * for each marginal variable -- Options are defined by + * {@link KdTree#setNormType(String)} and the + * default is {@link EuclideanUtils#NORM_MAX_NORM}. + */ + public final static String PROP_NORM_TYPE = "NORM_TYPE"; + /** + * Property name for an amount of random Gaussian noise to be + * added to the data (default 1e-8 to match the noise order in MILCA toolkit.). + */ + public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD"; + /** + * Property name for a dynamics exclusion time window + * otherwise known as Theiler window (see Kantz and Schreiber). + * Default is 0 which means no dynamic exclusion window. + */ + public static final String PROP_DYN_CORR_EXCL_TIME = "DYN_CORR_EXCL"; + /** + * Property name for the number of parallel threads to use in the + * computation (default is to use all available) + */ + public static final String PROP_NUM_THREADS = "NUM_THREADS"; + /** + * Valid property value for {@link #PROP_NUM_THREADS} to indicate + * that all available processors should be used. + */ + public static final String USE_ALL_THREADS = "USE_ALL"; + + /** + * 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; + /** + * Whether we use dynamic correlation exclusion + */ + protected boolean dynCorrExcl = false; + /** + * Size of dynamic correlation exclusion window. + */ + protected int dynCorrExclTime = 0; + /** + * Number of parallel threads to use in the computation; + * defaults to use all available. + */ + protected int numThreads = Runtime.getRuntime().availableProcessors(); + /** + * Protected k-d tree data structure (for fast nearest neighbour searches) + * representing the joint space + */ + protected KdTree kdTreeJoint; + // /** + // * protected data structures (for fast nearest neighbour searches) + // * representing the marginal spaces + // */ + // protected KdTree[] rangeSearchersInMarginals; + /** + * Protected data structures (for fast nearest neighbour searches) + * representing the marginal spaces of each individual variable + */ + protected UnivariateNearestNeighbourSearcher[] rangeSearchersInSmallMarginals; + /** + * Protected data structures (for fast nearest neighbour searches) + * representing the marginal spaces of each set of (D-1) variables + */ + protected KdTree[] rangeSearchersInBigMarginals; + /** + * Constant for digamma(k), with k the number of nearest neighbours selected + */ + protected double digammaK; + /** + * Constant for digamma(N), with N the number of samples. + */ + protected double digammaN; + + + public void initialise(int dimensions) { + this.dimensions = dimensions; + lastAverage = 0.0; + totalObservations = 0; + isComputed = false; + observations = null; + kdTreeJoint = null; + rangeSearchersInSmallMarginals = null; + rangeSearchersInBigMarginals = null; + } + + /** + * Sets properties for the KSG multivariate measure calculator. + * New property values are not guaranteed to take effect until the next call + * to an initialise method. + * + *

Valid property names, and what their + * values should represent, include:

+ *
    + *
  • {@link #PROP_K} -- number of k nearest neighbours to use in joint kernel space + * in the KSG algorithm (default is 4).
    • + *
  • {@link #PROP_NORM_TYPE} -- normalization type to apply to + * working out the norms between the points in each marginal space. + * Options are defined by {@link KdTree#setNormType(String)} - + * default is {@link EuclideanUtils#NORM_MAX_NORM}.
    • + *
  • {@link #PROP_DYN_CORR_EXCL_TIME} -- a dynamics exclusion time window, + * also known as Theiler window (see Kantz and Schreiber); + * default is 0 which means no dynamic exclusion window.
    • + *
  • {@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. + * (Recommended by Kraskov. MILCA uses 1e-8; but adds in + * a random amount of noise in [0,noiseLevel) ). + * Default 1e-8 to match the noise order in MILCA toolkit..
    • + *
+ * + *

Unknown property values are ignored.

+ * + * @param propertyName name of the property + * @param propertyValue value of the property + * @throws Exception for invalid property values + */ + public void setProperty(String propertyName, String propertyValue) throws Exception { + boolean propertySet = true; + if (propertyName.equalsIgnoreCase(PROP_K)) { + k = Integer.parseInt(propertyValue); + } else if (propertyName.equalsIgnoreCase(PROP_NORM_TYPE)) { + normType = KdTree.validateNormType(propertyValue); + } else if (propertyName.equalsIgnoreCase(PROP_DYN_CORR_EXCL_TIME)) { + dynCorrExclTime = Integer.parseInt(propertyValue); + dynCorrExcl = (dynCorrExclTime > 0); + } 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_NUM_THREADS)) { + if (propertyValue.equalsIgnoreCase(USE_ALL_THREADS)) { + numThreads = Runtime.getRuntime().availableProcessors(); + } else { // otherwise the user has passed in an integer: + numThreads = Integer.parseInt(propertyValue); + } + } else { + // No property was set here + propertySet = false; + // try the superclass: + super.setProperty(propertyName, propertyValue); + } + if (debug && propertySet) { + System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName + + " to " + propertyValue); + } + } + + @Override + public void finaliseAddObservations() throws Exception { + super.finaliseAddObservations(); + + if ((observations == null) || (observations[0].length == 0)) { + throw new Exception("Computing measure with a null set of data"); + } + if (observations.length <= k + 2*dynCorrExclTime) { + throw new Exception("There are less observations provided (" + + observations.length + + ") than required for the number of nearest neighbours parameter (" + + k + ") and any dynamic correlation exclusion (" + dynCorrExclTime + ")"); + } + + // Normalise the data if required + if (normalise) { + // We can overwrite these since they're already + // a copy of the users' data. + MatrixUtils.normalise(observations); + } + + // Add small random noise if required + if (addNoise) { + Random random = new Random(); + // Add Gaussian noise of std dev noiseLevel to the data + for (int r = 0; r < observations.length; r++) { + for (int c = 0; c < dimensions; c++) { + observations[r][c] += random.nextGaussian()*noiseLevel; + } + } + } + + // Set the constants: + digammaK = MathsUtils.digamma(k); + digammaN = MathsUtils.digamma(totalObservations); + } + + /** + * Internal method to ensure that the Kd-tree data structures to represent the + * observational data have been constructed (should be called prior to attempting + * to use these data structures) + */ + protected void ensureKdTreesConstructed() throws Exception { + + // We need to construct the k-d trees for use by the child + // classes. We check each tree for existence separately + // since source can be used across original and surrogate data + // TODO can parallelise these -- best done within the kdTree -- + // though it's unclear if there's much point given that + // the tree construction itself afterwards can't really be well parallelised. + if (kdTreeJoint == null) { + kdTreeJoint = new KdTree(observations); + kdTreeJoint.setNormType(normType); + } + if (rangeSearchersInSmallMarginals == null) { + rangeSearchersInSmallMarginals = new UnivariateNearestNeighbourSearcher[dimensions]; + for (int d = 0; d < dimensions; d++) { + rangeSearchersInSmallMarginals[d] = new UnivariateNearestNeighbourSearcher( + MatrixUtils.selectColumn(observations, d)); + rangeSearchersInSmallMarginals[d].setNormType(normType); + } + } + if (rangeSearchersInBigMarginals == null) { + rangeSearchersInBigMarginals = new KdTree[dimensions]; + for (int d = 0; d < dimensions; d++) { + rangeSearchersInBigMarginals[d] = new KdTree( + MatrixUtils.selectColumns(observations, allExcept(d, dimensions))); + rangeSearchersInBigMarginals[d].setNormType(normType); + } + } + + } + + /** + * {@inheritDoc} + * + * @return the average measure in nats (not bits!) + */ + public double computeAverageLocalOfObservations() throws Exception { + // Compute the measure + double startTime = Calendar.getInstance().getTimeInMillis(); + lastAverage = computeFromObservations(false)[0]; + isComputed = true; + if (debug) { + Calendar rightNow2 = Calendar.getInstance(); + long endTime = rightNow2.getTimeInMillis(); + System.out.println("Calculation time: " + ((endTime - startTime)/1000.0) + " sec" ); + } + return lastAverage; + } + + /** + * {@inheritDoc} + * + * @return the "time-series" of local measure values in nats (not bits!) + * @throws Exception + */ + public double[] computeLocalOfPreviousObservations() throws Exception { + double[] localValues = computeFromObservations(true); + lastAverage = MatrixUtils.mean(localValues); + isComputed = true; + return localValues; + } + + /** + * This method, specified in {@link MultiVariateInfoMeasureCalculatorCommon} + * is not implemented yet here. + */ + public double[] computeLocalUsingPreviousObservations(double[][] states) throws Exception { + // TODO If this is implemented, will need to normalise the incoming + // observations the same way that previously supplied ones were + // normalised (if they were normalised, that is) + throw new Exception("Local method not implemented yet"); + } + + /** + * This protected method handles the multiple threads which + * computes either the average or local measure (over parts of the total + * observations), computing the + * distances between all tuples in time. + * + *

The method returns:

    + *
  1. for (returnLocals == false), an array of size 1, + * containing the average measure
    2. + *
  2. for (returnLocals == true), the array of local + * measure values
    3. + *
+ * + * @param returnLocals whether to return an array or local values, or else + * sums of these values + * @return either the average measure, or array of local measure value, + * in nats not bits + * @throws Exception + */ + protected double[] computeFromObservations(boolean returnLocals) throws Exception { + + double[] returnValues = null; + + ensureKdTreesConstructed(); + + if (numThreads == 1) { + // Single-threaded implementation: + returnValues = partialComputeFromObservations(0, totalObservations, returnLocals); + + } else { + // We're going multithreaded: + if (returnLocals) { + // We're computing locals + returnValues = new double[totalObservations]; + } else { + // We're computing average + returnValues = new double[1]; + } + + // Distribute the observations to the threads for the parallel processing + int lTimesteps = totalObservations / numThreads; // each thread gets the same amount of data + int res = totalObservations % numThreads; // the first thread gets the residual data + if (debug) { + System.out.printf("Computing Kraskov Multi-Info with %d threads (%d timesteps each, plus %d residual)\n", + numThreads, lTimesteps, res); + } + Thread[] tCalculators = new Thread[numThreads]; + KraskovThreadRunner[] runners = new KraskovThreadRunner[numThreads]; + for (int t = 0; t < numThreads; t++) { + int startTime = (t == 0) ? 0 : lTimesteps * t + res; + int numTimesteps = (t == 0) ? lTimesteps + res : lTimesteps; + if (debug) { + System.out.println(t + ".Thread: from " + startTime + + " to " + (startTime + numTimesteps)); // Trace Message + } + runners[t] = new KraskovThreadRunner(this, startTime, numTimesteps, returnLocals); + tCalculators[t] = new Thread(runners[t]); + tCalculators[t].start(); + } + + // Here, we should wait for the termination of the all threads + // and collect their results + for (int t = 0; t < numThreads; t++) { + if (tCalculators[t] != null) { // TODO Ipek: can you comment on why we're checking for null here? + tCalculators[t].join(); + } + // Now we add in the data from this completed thread: + if (returnLocals) { + // We're computing local measure; copy these local values + // into the full array of locals + System.arraycopy(runners[t].getReturnValues(), 0, + returnValues, runners[t].myStartTimePoint, runners[t].numberOfTimePoints); + } else { + // We're computing the average measure, keep the running sums of digammas and counts + MatrixUtils.addInPlace(returnValues, runners[t].getReturnValues()); + } + } + } + + return returnValues; + + } + + /** + * Protected method to be used internally for threaded implementations. This + * method implements the guts of each Kraskov algorithm, computing the number + * of nearest neighbours in each dimension for a sub-set of the data points. + * It is intended to be called by one thread to work on that specific sub-set + * of the data. + * + *

Child classes should implement the computation of any specific measure + * using this method.

+ * + *

The method returns:

    + *
  1. for average measures (returnLocals == false), the relevant sums of + * digamma(n_x+1) in each marginal + * for a partial set of the observations
    2. + *
  2. for local measures (returnLocals == true), the array of local values
    3. + *
+ * + * @param startTimePoint start time for the partial set we examine + * @param numTimePoints number of time points (including startTimePoint to examine) + * @param returnLocals whether to return an array or local values, or else + * sums of these values + * @return an array of sum of digamma(n_x+1) for each marginal x, then + * sum of n_x for each marginal x (these latter ones are for debugging purposes). + * @throws Exception + */ + protected abstract double[] partialComputeFromObservations( + int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception; + + /** + * Private class to handle multi-threading of the Kraskov algorithms. + * Each instance calls partialComputeFromObservations() + * to compute nearest neighbours for a part of the data. + * + * + * @author Joseph Lizier (email, + * www) + * @author Ipek Özdemir + */ + private class KraskovThreadRunner implements Runnable { + protected MultiVariateInfoMeasureCalculatorKraskov calc; + protected int myStartTimePoint; + protected int numberOfTimePoints; + protected boolean computeLocals; + + protected double[] returnValues = null; + protected Exception problem = null; + + public static final int INDEX_SUM_DIGAMMAS = 0; + + public KraskovThreadRunner( + MultiVariateInfoMeasureCalculatorKraskov calc, + int myStartTimePoint, int numberOfTimePoints, + boolean computeLocals) { + this.calc = calc; + this.myStartTimePoint = myStartTimePoint; + this.numberOfTimePoints = numberOfTimePoints; + this.computeLocals = computeLocals; + } + + /** + * Return the values from this part of the data, + * or throw any exception that was encountered by the + * thread. + * + * @return an exception previously encountered by this thread. + * @throws Exception + */ + public double[] getReturnValues() throws Exception { + if (problem != null) { + throw problem; + } + return returnValues; + } + + /** + * Start the thread for the given parameters + */ + public void run() { + try { + returnValues = calc.partialComputeFromObservations( + myStartTimePoint, numberOfTimePoints, computeLocals); + } catch (Exception e) { + // Store the exception for later retrieval + problem = e; + return; + } + } + } + // end class KraskovThreadRunner + +} diff --git a/java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java new file mode 100755 index 0000000..dc411f1 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java @@ -0,0 +1,157 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.kraskov; + +import infodynamics.utils.EuclideanUtils; +import infodynamics.utils.NeighbourNodeData; +import infodynamics.utils.KdTree; +import infodynamics.utils.UnivariateNearestNeighbourSearcher; +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +import java.util.PriorityQueue; +import java.util.Calendar; +import java.util.Random; +import java.util.Arrays; + +/** + *

Computes the differential O-information of a given multivariate + * set of observations, using Kraskov-Stoegbauer-Grassberger (KSG) estimation + * (see Kraskov et al., below).

+ * + *

Usage is as per the paradigm outlined for + * {@link MultiVariateInfoMeasureCalculatorCommon}.

+ * + *

Finally, note that {@link Cloneable} is implemented allowing clone() + * to produce only an automatic shallow copy, which is fine + * for the statistical significance calculation it is intended for + * (none of the array + * data will be changed there). + *

+ * + *

References:
+ *

    + *
  • Rosas, F., Mediano, P., Gastpar, M., Jensen, H., + * "Quantifying high-order effects via multivariate extensions of the + * mutual information".
    • + * + *
  • Kraskov, A., Stoegbauer, H., Grassberger, P., + * "Estimating mutual information", + * Physical Review E 69, (2004) 066138.
    • + *
+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class OInfoCalculatorKraskov + extends MultiVariateInfoMeasureCalculatorKraskov + implements Cloneable { // See comments on clonability above + + + protected double[] partialComputeFromObservations( + int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception { + + // If data is 2D, return 0 before doing any computation + if (dimensions == 2) { + if (returnLocals) { + double[] localMi = new double[numTimePoints]; + Arrays.fill(localMi, 0); + return localMi; + } else { + return new double[] {0}; + } + } + + double startTime = Calendar.getInstance().getTimeInMillis(); + + double[] localMi = null; + if (returnLocals) { + localMi = new double[numTimePoints]; + } + + // Constants: + double dimensionsMinus1TimesDiGammaN = (double) (dimensions - 1) * digammaN; + + // Count the average number of points within eps_x for each marginal x of each point + double totalSumF = 0.0; + + for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) { + // Compute eps for this time step by + // finding the kth closest neighbour for point t: + PriorityQueue nnPQ = + kdTreeJoint.findKNearestNeighbours(k, t, dynCorrExclTime); + // First element in the PQ is the kth NN, + // and epsilon = kthNnData.distance + NeighbourNodeData kthNnData = nnPQ.poll(); + + // Distance to kth neighbour in joint space + double eps = kthNnData.distance; + + double sumF = 0.0; + + sumF += (digammaK - digammaN); + + for (int d = 0; d < dimensions; d++) { + int n_small = rangeSearchersInSmallMarginals[d].countPointsStrictlyWithinR( + t, eps, dynCorrExclTime); + + int n_big = rangeSearchersInBigMarginals[d].countPointsStrictlyWithinR( + t, eps, dynCorrExclTime); + + + sumF -= (MathsUtils.digamma(n_big + 1) - digammaN)/(dimensions - 2); + sumF += (MathsUtils.digamma(n_small + 1) - digammaN)/(dimensions - 2); + + if (debug) { + // Only tracking this for debugging purposes: + System.out.printf("t=%d, d=%d, n_small=%d, n_big=%d, sumF=%.3f\n", + t, d, n_small, n_big, sumF); + } + + } + + sumF *= (2 - dimensions); + + totalSumF += sumF; + + if (returnLocals) { + localMi[t-startTimePoint] = sumF; + } + } + + if (debug) { + Calendar rightNow2 = Calendar.getInstance(); + long endTime = rightNow2.getTimeInMillis(); + System.out.println("Subset " + startTimePoint + ":" + + (startTimePoint + numTimePoints) + " Calculation time: " + + ((endTime - startTime)/1000.0) + " sec" ); + } + + // Select what to return: + if (returnLocals) { + return localMi; + } else { + double[] returnArray = new double[] {totalSumF/((double) totalObservations)}; + return returnArray; + } + + } + +} + diff --git a/java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java new file mode 100755 index 0000000..2bfb800 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java @@ -0,0 +1,146 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.kraskov; + +import infodynamics.utils.EuclideanUtils; +import infodynamics.utils.NeighbourNodeData; +import infodynamics.utils.KdTree; +import infodynamics.utils.UnivariateNearestNeighbourSearcher; +import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; + +import java.util.PriorityQueue; +import java.util.Calendar; +import java.util.Random; +import java.util.Arrays; + +/** + *

Computes the differential S-information of a given multivariate + * set of observations, using Kraskov-Stoegbauer-Grassberger (KSG) estimation + * (see Kraskov et al., below).

+ * + *

Usage is as per the paradigm outlined for + * {@link MultiVariateInfoMeasureCalculatorCommon}.

+ * + *

Finally, note that {@link Cloneable} is implemented allowing clone() + * to produce only an automatic shallow copy, which is fine + * for the statistical significance calculation it is intended for + * (none of the array + * data will be changed there). + *

+ * + *

References:
+ *

    + *
  • Rosas, F., Mediano, P., Gastpar, M., Jensen, H., + * "Quantifying high-order effects via multivariate extensions of the + * mutual information".
    • + * + *
  • Kraskov, A., Stoegbauer, H., Grassberger, P., + * "Estimating mutual information", + * Physical Review E 69, (2004) 066138.
    • + *
+ * + * @author Pedro A.M. Mediano (email, + * www) + */ +public class SInfoCalculatorKraskov + extends MultiVariateInfoMeasureCalculatorKraskov + implements Cloneable { // See comments on clonability above + + + protected double[] partialComputeFromObservations( + int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception { + + double startTime = Calendar.getInstance().getTimeInMillis(); + + double[] localMi = null; + if (returnLocals) { + localMi = new double[numTimePoints]; + } + + // Constants: + double dimensionsMinus1TimesDiGammaN = (double) (dimensions - 1) * digammaN; + + // Count the average number of points within eps_x for each marginal x of each point + double totalSumF = 0.0; + + for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) { + // Compute eps for this time step by + // finding the kth closest neighbour for point t: + PriorityQueue nnPQ = + kdTreeJoint.findKNearestNeighbours(k, t, dynCorrExclTime); + // First element in the PQ is the kth NN, + // and epsilon = kthNnData.distance + NeighbourNodeData kthNnData = nnPQ.poll(); + + // Distance to kth neighbour in joint space + double eps = kthNnData.distance; + + double sumF = 0.0; + + sumF += (digammaK - digammaN); + + for (int d = 0; d < dimensions; d++) { + int n_small = rangeSearchersInSmallMarginals[d].countPointsStrictlyWithinR( + t, eps, dynCorrExclTime); + + int n_big = rangeSearchersInBigMarginals[d].countPointsStrictlyWithinR( + t, eps, dynCorrExclTime); + + + sumF -= (MathsUtils.digamma(n_big + 1) - digammaN)/dimensions; + sumF -= (MathsUtils.digamma(n_small + 1) - digammaN)/dimensions; + + if (debug) { + // Only tracking this for debugging purposes: + System.out.printf("t=%d, d=%d, n_small=%d, n_big=%d, sumF=%.3f\n", + t, d, n_small, n_big, sumF); + } + + } + + sumF *= dimensions; + + totalSumF += sumF; + + if (returnLocals) { + localMi[t-startTimePoint] = sumF; + } + } + + if (debug) { + Calendar rightNow2 = Calendar.getInstance(); + long endTime = rightNow2.getTimeInMillis(); + System.out.println("Subset " + startTimePoint + ":" + + (startTimePoint + numTimePoints) + " Calculation time: " + + ((endTime - startTime)/1000.0) + " sec" ); + } + + // Select what to return: + if (returnLocals) { + return localMi; + } else { + double[] returnArray = new double[] {totalSumF/((double) totalObservations)}; + return returnArray; + } + + } + +} + diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java new file mode 100644 index 0000000..fbaf88b --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java @@ -0,0 +1,143 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.kraskov; + +import infodynamics.measures.continuous.gaussian.DualTotalCorrelationCalculatorGaussian; + +import infodynamics.utils.ArrayFileReader; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; +import junit.framework.TestCase; + +public class DualTotalCorrelationCalculatorKraskovTester extends TestCase { + + /** + * For two variables, DTC is equal to mutual information. + */ + public void testTwoVariables() throws Exception { + + double[][] data; + ArrayFileReader afr = new ArrayFileReader("demos/data/2randomCols-1.txt"); + data = afr.getDouble2DMatrix(); + + DualTotalCorrelationCalculatorKraskov dtcCalc = new DualTotalCorrelationCalculatorKraskov(); + dtcCalc.setProperty("NOISE_LEVEL_TO_ADD", "0"); + dtcCalc.initialise(2); + dtcCalc.setObservations(data); + double dtc = dtcCalc.computeAverageLocalOfObservations(); + + + MutualInfoCalculatorMultiVariateKraskov1 miCalc = new MutualInfoCalculatorMultiVariateKraskov1(); + miCalc.setProperty("NOISE_LEVEL_TO_ADD", "0"); + miCalc.initialise(1,1); + miCalc.setObservations(MatrixUtils.selectColumn(data, 0), + MatrixUtils.selectColumn(data, 1)); + double mi = miCalc.computeAverageLocalOfObservations(); + + assertEquals(dtc, mi, 1e-6); + } + + /** + * Compare against the values obtained by the Gaussian DTC calculator when the data + * is actually Gaussian. + */ + public void testCompareWithGaussian() throws Exception { + + int N = 100000, D = 3; + RandomGenerator rg = new RandomGenerator(); + rg.setSeed(1); + double[][] data = rg.generateNormalData(N, D, 0, 1); + for (int i = 0; i < N; i++) { + data[i][2] = data[i][2] + 0.1*data[i][0]; + data[i][1] = data[i][1] + 0.5*data[i][0]; + } + + DualTotalCorrelationCalculatorKraskov dtcCalc_ksg = new DualTotalCorrelationCalculatorKraskov(); + dtcCalc_ksg.setProperty("NOISE_LEVEL_TO_ADD", "0"); + dtcCalc_ksg.initialise(3); + dtcCalc_ksg.setObservations(data); + double dtc_ksg = dtcCalc_ksg.computeAverageLocalOfObservations(); + + DualTotalCorrelationCalculatorGaussian dtcCalc_gau = new DualTotalCorrelationCalculatorGaussian(); + dtcCalc_gau.initialise(3); + dtcCalc_gau.setObservations(data); + double dtc_gau = dtcCalc_gau.computeAverageLocalOfObservations(); + + assertEquals(dtc_ksg, dtc_gau, 0.001); + + } + + /** + * Confirm that the local values average correctly back to the average value + * + */ + public void testLocalsAverageCorrectly() throws Exception { + + int dimensions = 4; + int timeSteps = 1000; + DualTotalCorrelationCalculatorKraskov dtcCalc = new DualTotalCorrelationCalculatorKraskov(); + dtcCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + dtcCalc.setObservations(data); + + double dtc = dtcCalc.computeAverageLocalOfObservations(); + double[] dtcLocal = dtcCalc.computeLocalOfPreviousObservations(); + + System.out.printf("Average was %.5f\n", dtc); + + assertEquals(dtc, MatrixUtils.mean(dtcLocal), 0.00001); + } + + /** + * Confirm that for 2D the local values equal the local MI values + * + */ + public void testLocalsEqualMI() throws Exception { + + int dimensions = 2; + int timeSteps = 100; + DualTotalCorrelationCalculatorKraskov dtcCalc = new DualTotalCorrelationCalculatorKraskov(); + dtcCalc.initialise(dimensions); + + // generate some random data + RandomGenerator rg = new RandomGenerator(); + double[][] data = rg.generateNormalData(timeSteps, dimensions, + 0, 1); + + dtcCalc.setObservations(data); + double[] dtcLocal = dtcCalc.computeLocalOfPreviousObservations(); + + MutualInfoCalculatorMultiVariateKraskov1 miCalc = new MutualInfoCalculatorMultiVariateKraskov1(); + miCalc.initialise(1, 1); + miCalc.setObservations(MatrixUtils.selectColumn(data, 0), MatrixUtils.selectColumn(data, 1)); + double[] miLocal = miCalc.computeLocalOfPreviousObservations(); + + for (int t = 0; t < timeSteps; t++) { + assertEquals(dtcLocal[t], miLocal[t], 0.00001); + } + } + + +} + diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskovTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskovTester.java new file mode 100644 index 0000000..6b82f55 --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskovTester.java @@ -0,0 +1,62 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.kraskov; + +import infodynamics.measures.continuous.gaussian.OInfoCalculatorGaussian; + +import infodynamics.utils.ArrayFileReader; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; +import junit.framework.TestCase; + +public class OInfoCalculatorKraskovTester extends TestCase { + + /** + * Compare against the values obtained by the Gaussian O-info calculator when the data + * is actually Gaussian. + */ + public void testCompareWithGaussian() throws Exception { + + int N = 100000, D = 3; + RandomGenerator rg = new RandomGenerator(); + rg.setSeed(3); + double[][] data = rg.generateNormalData(N, D, 0, 1); + for (int i = 0; i < N; i++) { + data[i][2] = data[i][2] + 0.4*data[i][0]; + data[i][1] = data[i][1] + 0.5*data[i][0]; + } + + OInfoCalculatorKraskov oCalc_ksg = new OInfoCalculatorKraskov(); + oCalc_ksg.setProperty("NOISE_LEVEL_TO_ADD", "0"); + oCalc_ksg.initialise(3); + oCalc_ksg.setObservations(data); + double o_ksg = oCalc_ksg.computeAverageLocalOfObservations(); + + OInfoCalculatorGaussian oCalc_gau = new OInfoCalculatorGaussian(); + oCalc_gau.initialise(3); + oCalc_gau.setObservations(data); + double o_gau = oCalc_gau.computeAverageLocalOfObservations(); + + assertEquals(o_ksg, o_gau, 0.001); + + } + +} + + diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskovTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskovTester.java new file mode 100644 index 0000000..442c03e --- /dev/null +++ b/java/unittests/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskovTester.java @@ -0,0 +1,62 @@ +/* + * Java Information Dynamics Toolkit (JIDT) + * Copyright (C) 2012, 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 . + */ + +package infodynamics.measures.continuous.kraskov; + +import infodynamics.measures.continuous.gaussian.SInfoCalculatorGaussian; + +import infodynamics.utils.ArrayFileReader; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; +import junit.framework.TestCase; + +public class SInfoCalculatorKraskovTester extends TestCase { + + /** + * Compare against the values obtained by the Gaussian S-info calculator when the data + * is actually Gaussian. + */ + public void testCompareWithGaussian() throws Exception { + + int N = 100000, D = 3; + RandomGenerator rg = new RandomGenerator(); + rg.setSeed(3); + double[][] data = rg.generateNormalData(N, D, 0, 1); + for (int i = 0; i < N; i++) { + data[i][2] = data[i][2] + 0.4*data[i][0]; + data[i][1] = data[i][1] + 0.5*data[i][0]; + } + + SInfoCalculatorKraskov sCalc_ksg = new SInfoCalculatorKraskov(); + sCalc_ksg.setProperty("NOISE_LEVEL_TO_ADD", "0"); + sCalc_ksg.initialise(3); + sCalc_ksg.setObservations(data); + double s_ksg = sCalc_ksg.computeAverageLocalOfObservations(); + + SInfoCalculatorGaussian sCalc_gau = new SInfoCalculatorGaussian(); + sCalc_gau.initialise(3); + sCalc_gau.setObservations(data); + double s_gau = sCalc_gau.computeAverageLocalOfObservations(); + + assertEquals(s_ksg, s_gau, 0.001); + + } + +} + + From 07883af73430b7760d1174564323f2e517c0fd9f Mon Sep 17 00:00:00 2001 From: Pedro Mediano Date: Tue, 9 Feb 2021 10:23:30 +0000 Subject: [PATCH 8/8] Updated printed messages with cross-platform newline character. --- .../kraskov/MultiVariateInfoMeasureCalculatorKraskov.java | 2 +- .../measures/continuous/kraskov/OInfoCalculatorKraskov.java | 2 +- .../measures/continuous/kraskov/SInfoCalculatorKraskov.java | 2 +- .../gaussian/DualTotalCorrelationCalculatorGaussianTester.java | 2 +- .../continuous/gaussian/OInfoCalculatorGaussianTester.java | 2 +- .../continuous/gaussian/SInfoCalculatorGaussianTester.java | 2 +- .../kraskov/DualTotalCorrelationCalculatorKraskovTester.java | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java index df8bcf4..ad925f3 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MultiVariateInfoMeasureCalculatorKraskov.java @@ -405,7 +405,7 @@ public abstract class MultiVariateInfoMeasureCalculatorKraskov int lTimesteps = totalObservations / numThreads; // each thread gets the same amount of data int res = totalObservations % numThreads; // the first thread gets the residual data if (debug) { - System.out.printf("Computing Kraskov Multi-Info with %d threads (%d timesteps each, plus %d residual)\n", + System.out.printf("Computing Kraskov Multi-Info with %d threads (%d timesteps each, plus %d residual)%n", numThreads, lTimesteps, res); } Thread[] tCalculators = new Thread[numThreads]; diff --git a/java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java index dc411f1..3410d73 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/OInfoCalculatorKraskov.java @@ -120,7 +120,7 @@ public class OInfoCalculatorKraskov if (debug) { // Only tracking this for debugging purposes: - System.out.printf("t=%d, d=%d, n_small=%d, n_big=%d, sumF=%.3f\n", + System.out.printf("t=%d, d=%d, n_small=%d, n_big=%d, sumF=%.3f%n", t, d, n_small, n_big, sumF); } diff --git a/java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java index 2bfb800..165db44 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/SInfoCalculatorKraskov.java @@ -109,7 +109,7 @@ public class SInfoCalculatorKraskov if (debug) { // Only tracking this for debugging purposes: - System.out.printf("t=%d, d=%d, n_small=%d, n_big=%d, sumF=%.3f\n", + System.out.printf("t=%d, d=%d, n_small=%d, n_big=%d, sumF=%.3f%n", t, d, n_small, n_big, sumF); } diff --git a/java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java b/java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java index a5cd702..0bcf7c6 100644 --- a/java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java +++ b/java/unittests/infodynamics/measures/continuous/gaussian/DualTotalCorrelationCalculatorGaussianTester.java @@ -95,7 +95,7 @@ public class DualTotalCorrelationCalculatorGaussianTester extends TestCase { double dtc = dtcCalc.computeAverageLocalOfObservations(); double[] dtcLocal = dtcCalc.computeLocalOfPreviousObservations(); - System.out.printf("Average was %.5f\n", dtc); + System.out.printf("Average was %.5f%n", dtc); assertEquals(dtc, MatrixUtils.mean(dtcLocal), 0.00001); } diff --git a/java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java b/java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java index 15e3c46..def82ad 100644 --- a/java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java +++ b/java/unittests/infodynamics/measures/continuous/gaussian/OInfoCalculatorGaussianTester.java @@ -118,7 +118,7 @@ public class OInfoCalculatorGaussianTester extends TestCase { double oinfo = oCalc.computeAverageLocalOfObservations(); double[] oLocal = oCalc.computeLocalOfPreviousObservations(); - System.out.printf("Average was %.5f\n", oinfo); + System.out.printf("Average was %.5f%n", oinfo); assertEquals(oinfo, MatrixUtils.mean(oLocal), 0.00001); } diff --git a/java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java b/java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java index bc33ad7..79cae70 100644 --- a/java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java +++ b/java/unittests/infodynamics/measures/continuous/gaussian/SInfoCalculatorGaussianTester.java @@ -92,7 +92,7 @@ public class SInfoCalculatorGaussianTester extends TestCase { double sinfo = sCalc.computeAverageLocalOfObservations(); double[] sLocal = sCalc.computeLocalOfPreviousObservations(); - System.out.printf("Average was %.5f\n", sinfo); + System.out.printf("Average was %.5f%n", sinfo); assertEquals(sinfo, MatrixUtils.mean(sLocal), 0.00001); } diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java index fbaf88b..046cd7e 100644 --- a/java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java +++ b/java/unittests/infodynamics/measures/continuous/kraskov/DualTotalCorrelationCalculatorKraskovTester.java @@ -104,7 +104,7 @@ public class DualTotalCorrelationCalculatorKraskovTester extends TestCase { double dtc = dtcCalc.computeAverageLocalOfObservations(); double[] dtcLocal = dtcCalc.computeLocalOfPreviousObservations(); - System.out.printf("Average was %.5f\n", dtc); + System.out.printf("Average was %.5f%n", dtc); assertEquals(dtc, MatrixUtils.mean(dtcLocal), 0.00001); }