Computes the differential multi-information of a given multivariate set of @@ -46,8 +48,11 @@ import java.util.Vector; * *
* - *- * TODO Add fast nearest neighbour searches to the child classes + *
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: 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 AIS
- * calculations. 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. The method returns: The method returns: Computes the differential multi-information of two given multivariate
@@ -45,272 +47,79 @@ import infodynamics.utils.MatrixUtils;
*
* @author Joseph Lizier (email,
* www)
+ * @author Ipek Özdemir
*/
public class MultiInfoCalculatorKraskov1
extends MultiInfoCalculatorKraskov {
- @Override
- public double computeAverageLocalOfObservations() throws Exception {
- if (miComputed) {
- return mi;
- }
- return computeAverageLocalOfObservations(null);
+ public MultiInfoCalculatorKraskov1() {
+ super();
+ isAlgorithm1 = true;
}
@Override
- public double computeAverageLocalOfObservations(int[][] reordering) throws Exception {
- if (V == 1) {
- miComputed = true;
- return 0.0;
- }
- if (!tryKeepAllPairsNorms || (data.length * V > MAX_DATA_SIZE_FOR_KEEP_ALL_PAIRS_NORM)) {
- double[][] originalData = data;
- if (reordering != null) {
- // Generate a new re-ordered data
- data = MatrixUtils.reorderDataForVariables(originalData, reordering);
- }
- // Compute the MI
- double newMI = computeAverageLocalOfObservationsWhileComputingDistances();
- // restore data
- data = originalData;
- return newMI;
+ 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];
}
- if (norms == null) {
- computeNorms();
- }
-
- // Count the average number of points within eps_x and eps_y
- double averageDiGammas = 0;
- double[] avNx = new double[V];
-
- int cutoffForKthMinLinear = (int) (Math.log(N) / Math.log(2.0));
-
- for (int t = 0; t < N; t++) {
- // Compute eps for this time step:
- // First grab marginal norms to all neighbours
- // (note that norm of point t to itself will be set to infinity).
-
- // Create storage for the reordered time steps for the variables
- int[] tForEachMarginal = reorderedTimeStepsForEachMarginal(reordering, t);
-
- double[] jointNorm = new double[N];
- for (int t2 = 0; t2 < N; t2++) {
- int[] t2ForEachMarginal = reorderedTimeStepsForEachMarginal(reordering, t2);
- // Find the max marginal norm between the vector at t and the reordered vector
- // at t2:
- jointNorm[t2] = 0;
- for (int v = 0; v < V; v++) {
- double normForThisVar = norms[v][tForEachMarginal[v]][t2ForEachMarginal[v]];
- if (normForThisVar > jointNorm[t2]) {
- jointNorm[t2] = normForThisVar;
- }
- }
- }
- // Then find the kth closest neighbour:
- double epsilon = 0.0;
- if (k <= cutoffForKthMinLinear) {
- // just do a linear search for the minimum
- epsilon = MatrixUtils.kthMin(jointNorm, k);
- } else {
- // Sort the array of joint norms first
- java.util.Arrays.sort(jointNorm);
- // And find the distance to it's kth closest neighbour
- // (we subtract one since the array is indexed from zero)
- epsilon = jointNorm[k-1];
- }
-
- // Count the number of points (in each marginal variable)
- // whose marginal distance is less than eps
- int[] n_x = new int[V];
- for (int t2 = 0; t2 < N; t2++) {
- int[] t2ForEachMarginal = reorderedTimeStepsForEachMarginal(reordering, t2);
- for (int v = 0; v < V; v++) {
- if (norms[v][tForEachMarginal[v]][t2ForEachMarginal[v]] < epsilon) {
- n_x[v]++;
- }
- }
- }
- // Track the averages, and take the digamma before adding into the
- // average:
- for (int v = 0; v < V; v++) {
- avNx[v] += n_x[v];
- averageDiGammas += MathsUtils.digamma(n_x[v]+1);
- }
- }
- averageDiGammas /= (double) N;
- if (debug) {
- for (int v = 0; v < V; v++) {
- avNx[v] /= (double)N;
- System.out.print(String.format("Average n_x[%d]=%.3f, ", v, avNx[v]));
- }
- System.out.println();
- }
-
- mi = MathsUtils.digamma(k) - averageDiGammas + (double) (V - 1) * MathsUtils.digamma(N);
- miComputed = true;
- return mi;
- }
-
- /**
- * This method correctly computes the average multi-info, but recomputes the
- * marginal distances between all tuples in time.
- * Kept here for cases where we have too many observations
- * to keep the norm between all pairs, and for testing purposes.
- *
- * @see #computeAverageLocalOfObservations()
- * @return average multi-info
- * @throws Exception
- */
- public double computeAverageLocalOfObservationsWhileComputingDistances() throws Exception {
-
- if (V == 1) {
- miComputed = true;
- return 0.0;
- }
- // Count the average number of points within eps for each marginal variable
- double averageDiGammas = 0;
- double[] avNx = new double[V];
- int cutoffForKthMinLinear = (int) (Math.log(N) / Math.log(2.0));
-
- for (int t = 0; t < N; t++) {
- // Compute eps for this time step:
- // First get marginal norms to all neighbours
- // (note that norm of point t to itself will be set to infinity).
-
- double[][] normsForT = EuclideanUtils.computeNorms(data, t);
- double[] jointNorm = new double[N];
- for (int t2 = 0; t2 < N; t2++) {
- jointNorm[t2] = MatrixUtils.max(normsForT[t2]);
- }
- // Then find the kth closest neighbour:
- double epsilon = 0.0;
- if (k <= cutoffForKthMinLinear) {
- // just do a linear search for the minimum
- epsilon = MatrixUtils.kthMin(jointNorm, k);
- } else {
- // Sort the array of joint norms first
- java.util.Arrays.sort(jointNorm);
- // And find the distance to it's kth closest neighbour
- // (we subtract one since the array is indexed from zero)
- epsilon = jointNorm[k-1];
- }
-
- // Count the number of points (in each marginal variable)
- // whose marginal distance is less than eps
- int[] n_x = new int[V];
- for (int t2 = 0; t2 < N; t2++) {
- for (int v = 0; v < V; v++) {
- if (normsForT[t2][v] < epsilon) {
- n_x[v]++;
- }
- }
- }
- for (int v = 0; v < V; v++) {
- avNx[v] += n_x[v];
- }
- // And take the digamma before adding into the
- // average:
- for (int v = 0; v < V; v++) {
- averageDiGammas += MathsUtils.digamma(n_x[v]+1);
- }
- }
- averageDiGammas /= (double) N;
- if (debug) {
- for (int v = 0; v < V; v++) {
- avNx[v] /= (double)N;
- System.out.print(String.format("Average n_x[%d]=%.3f, ", v, avNx[v]));
- }
- System.out.println();
- }
-
- mi = MathsUtils.digamma(k) - averageDiGammas + (double) (V - 1) * MathsUtils.digamma(N);
- miComputed = true;
- return mi;
- }
-
- @Override
- public double[] computeLocalOfPreviousObservations() throws Exception {
- double[] localMi = new double[N];
- int cutoffForKthMinLinear = (int) (Math.log(N) / Math.log(2.0));
-
- if (V == 1) {
- miComputed = true;
- return localMi;
- }
-
// Constants:
- double digammaK = MathsUtils.digamma(k);
- double Vminus1TimesdigammaN = (double) (V - 1) * MathsUtils.digamma(N);
-
- // Count the average number of points within eps_x[v] for each marginal v
- double averageDiGammas = 0;
- double[] avNx = new double[V];
-
- for (int t = 0; t < N; t++) {
- // Compute eps for this time step:
- // First get marginal norms to all neighbours
- // (note that norm of point t to itself will be set to infinity).
-
- double[][] norms = EuclideanUtils.computeNorms(data, t);
- double[] jointNorm = new double[N];
- for (int t2 = 0; t2 < N; t2++) {
- jointNorm[t2] = MatrixUtils.max(norms[t2]);
- }
-
- // Then find the kth closest neighbour:
- double epsilon = 0.0;
- if (k <= cutoffForKthMinLinear) {
- // just do a linear search for the minimum
- epsilon = MatrixUtils.kthMin(jointNorm, k);
- } else {
- // Sort the array of joint norms first
- java.util.Arrays.sort(jointNorm);
- // And find the distance to it's kth closest neighbour
- // (we subtract one since the array is indexed from zero)
- epsilon = jointNorm[k-1];
- }
-
- // Count the number of points (in each marginal variable)
- // whose marginal distance is less than eps
- int[] n_x = new int[V];
- for (int t2 = 0; t2 < N; t2++) {
- for (int v = 0; v < V; v++) {
- if (norms[t2][v] < epsilon) {
- n_x[v]++;
- }
- }
- }
+ double dimensionsMinus1TimesDiGammaN = (double) (dimensions - 1) * digammaN;
- // And take the digammas, and add into the local
- localMi[t] = digammaK + Vminus1TimesdigammaN;
- for (int v = 0; v < V; v++) {
- double digammaNxPlusOne = MathsUtils.digamma(n_x[v]+1);
- localMi[t] -= digammaNxPlusOne;
- // And keep track of the averages
- averageDiGammas += digammaNxPlusOne;
- avNx[v] += n_x[v];
+ // Count the average number of points within eps_x for each marginal x of each point
+ double sumDiGammas = 0;
+ double[] sumNMarginals = new double[dimensions];
+
+ for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
+ // Compute eps for this time step by
+ // finding the kth closest neighbour for point t:
+ PriorityQueue Computes the differential multi-information of two given multivariate
@@ -53,315 +54,84 @@ import infodynamics.utils.MatrixUtils;
public class MultiInfoCalculatorKraskov2
extends MultiInfoCalculatorKraskov {
- protected static final int JOINT_NORM_VAL_COLUMN = 0;
- protected static final int JOINT_NORM_TIMESTEP_COLUMN = 1;
-
- @Override
- public double computeAverageLocalOfObservations() throws Exception {
- if (miComputed) {
- return mi;
- }
- return computeAverageLocalOfObservations(null);
+ public MultiInfoCalculatorKraskov2() {
+ super();
+ isAlgorithm1 = false;
}
-
- @Override
- public double computeAverageLocalOfObservations(int[][] reordering) throws Exception {
- if (V == 1) {
- miComputed = true;
- return 0.0;
- }
- if (!tryKeepAllPairsNorms || (data.length * V > MAX_DATA_SIZE_FOR_KEEP_ALL_PAIRS_NORM)) {
- double[][] originalData = data;
- // Generate a new re-ordered data
- if (reordering != null) {
- // Generate a new re-ordered data
- data = MatrixUtils.reorderDataForVariables(originalData, reordering);
- }
- // Compute the MI
- double newMI = computeAverageLocalOfObservationsWhileComputingDistances();
- // restore data
- data = originalData;
- return newMI;
- }
-
- // Otherwise we will use the norms we've already computed, and use a "virtual"
- // reordered data2.
-
- if (norms == null) {
- computeNorms();
- }
-
- // Count the average number of points within eps_x[v]
- double averageDiGammas = 0;
- double[] avNx = new double[V];
-
- for (int t = 0; t < N; t++) {
- // Compute eps for each marginal for this time step:
- // First grab marginal norms to all neighbours
- // (note that norm of point t to itself will be set to infinity).
-
- // Get the reordered time steps for the variables
- int[] tForEachMarginal = reorderedTimeStepsForEachMarginal(reordering, t);
-
- double[][] jointNorm = new double[N][2];
- for (int t2 = 0; t2 < N; t2++) {
- int[] t2ForEachMarginal = reorderedTimeStepsForEachMarginal(reordering, t2);
- // Find the max marginal norm between the vector at t and the reordered vector
- // at t2:
- jointNorm[t2][JOINT_NORM_VAL_COLUMN] = 0;
- for (int v = 0; v < V; v++) {
- double normForThisVar = norms[v][tForEachMarginal[v]][t2ForEachMarginal[v]];
- if (normForThisVar > jointNorm[t2][JOINT_NORM_VAL_COLUMN]) {
- jointNorm[t2][JOINT_NORM_VAL_COLUMN] = normForThisVar;
- }
- }
- // And store the time step for back reference after the
- // array is sorted.
- jointNorm[t2][JOINT_NORM_TIMESTEP_COLUMN] = t2;
- }
- // Then find the k closest neighbours:
- double[] eps_x = new double[V];
- if (k == 1) {
- // just do a linear search for the minimum epsilon value
- int timeStepOfMin = MatrixUtils.minIndex(jointNorm, JOINT_NORM_VAL_COLUMN);
- int[] timeStepOfMinForEachMarginal =
- reorderedTimeStepsForEachMarginal(reordering, timeStepOfMin);
- for (int v = 0; v < V; v++) {
- eps_x[v] = norms[v][tForEachMarginal[v]][timeStepOfMinForEachMarginal[v]];
- }
- } else {
- // Sort the array of joint norms
- java.util.Arrays.sort(jointNorm, FirstIndexComparatorDouble.getInstance());
- // and now we have the closest k points.
- // Find eps_{x,y} as the maximum x and y norms amongst this set:
- for (int j = 0; j < k; j++) {
- int timeStepOfJthPoint = (int)jointNorm[j][JOINT_NORM_TIMESTEP_COLUMN];
- int[] timeStepOfJthPointForEachMarginal =
- reorderedTimeStepsForEachMarginal(reordering, timeStepOfJthPoint);
- for (int v = 0; v < V; v++) {
- if (norms[v][tForEachMarginal[v]][timeStepOfJthPointForEachMarginal[v]] > eps_x[v]) {
- eps_x[v] = norms[v][tForEachMarginal[v]][timeStepOfJthPointForEachMarginal[v]];
- }
- }
- }
- }
-
- // Count the number of points (in each marginal variable)
- // whose marginal distance is less than eps in that marginal dimension
- int[] n_x = new int[V];
- for (int t2 = 0; t2 < N; t2++) {
- int[] t2ForEachMarginal = reorderedTimeStepsForEachMarginal(reordering, t2);
- for (int v = 0; v < V; v++) {
- if (norms[v][tForEachMarginal[v]][t2ForEachMarginal[v]] <= eps_x[v]) {
- n_x[v]++;
- }
- }
- }
- // Track the averages, and take the digamma before adding into the
- // average:
- for (int v = 0; v < V; v++) {
- avNx[v] += n_x[v];
- averageDiGammas += MathsUtils.digamma(n_x[v]);
- }
- }
- averageDiGammas /= (double) N;
- if (debug) {
- for (int v = 0; v < V; v++) {
- avNx[v] /= (double)N;
- System.out.print(String.format("Average n_x[%d]=%.3f, ", v, avNx[v]));
- }
- System.out.println();
- }
-
- mi = MathsUtils.digamma(k) - (double) (V - 1) /(double)k - averageDiGammas +
- (double) (V - 1) * MathsUtils.digamma(N);
- miComputed = true;
- return mi;
- }
-
- /**
- * This method correctly computes the average multi-info, but recomputes the
- * marginal distances between all tuples in time.
- * Kept here for cases where we have too many observations
- * to keep the norm between all pairs, and for testing purposes.
- *
- * @see #computeAverageLocalOfObservations()
- * @return average multi-info
- * @throws Exception
- */
- public double computeAverageLocalOfObservationsWhileComputingDistances() throws Exception {
-
- if (V == 1) {
- miComputed = true;
- return 0.0;
- }
- // Count the average number of points within eps for each marginal variable
- double averageDiGammas = 0;
- double[] avNx = new double[V];
- for (int t = 0; t < N; t++) {
- // Compute eps_x (for each marginal) for this time step:
- // First get the marginal norms to all neighbours
- // (note that norm of point t to itself will be set to infinity).
-
- double[][] normsForT = EuclideanUtils.computeNorms(data, t);
- double[][] jointNorm = new double[N][2];
- for (int t2 = 0; t2 < N; t2++) {
- jointNorm[t2][JOINT_NORM_VAL_COLUMN] = MatrixUtils.max(normsForT[t2]);
- // And store the time step for back reference after the
- // array is sorted.
- jointNorm[t2][JOINT_NORM_TIMESTEP_COLUMN] = t2;
- }
-
- // Then find the k closest neighbours:
- double[] eps_x = new double[V];
- if (k == 1) {
- // just do a linear search for the minimum epsilon value
- int timeStepOfMin = MatrixUtils.minIndex(jointNorm, JOINT_NORM_VAL_COLUMN);
- for (int v = 0; v < V; v++) {
- eps_x[v] = normsForT[timeStepOfMin][v];
- }
- } else {
- // Sort the array of joint norms
- java.util.Arrays.sort(jointNorm, FirstIndexComparatorDouble.getInstance());
- // and now we have the closest k points.
- // Find eps_{x,y} as the maximum x and y norms amongst this set:
- for (int j = 0; j < k; j++) {
- int timeStepOfJthPoint = (int)jointNorm[j][JOINT_NORM_TIMESTEP_COLUMN];
- for (int v = 0; v < V; v++) {
- if (normsForT[timeStepOfJthPoint][v] > eps_x[v]) {
- eps_x[v] = normsForT[timeStepOfJthPoint][v];
- }
- }
- }
- }
-
- // Count the number of points (in each marginal variable)
- // whose marginal distance is less than eps in that marginal dimension
- int[] n_x = new int[V];
- for (int t2 = 0; t2 < N; t2++) {
- for (int v = 0; v < V; v++) {
- if (normsForT[t2][v] <= eps_x[v]) {
- n_x[v]++;
- }
- }
- }
- // Track the averages, and take the digamma before adding into the
- // average:
- for (int v = 0; v < V; v++) {
- avNx[v] += n_x[v];
- averageDiGammas += MathsUtils.digamma(n_x[v]);
- }
- }
- averageDiGammas /= (double) N;
- if (debug) {
- for (int v = 0; v < V; v++) {
- avNx[v] /= (double)N;
- System.out.print(String.format("Average n_x[%d]=%.3f, ", v, avNx[v]));
- }
- System.out.println();
+ 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];
}
- mi = MathsUtils.digamma(k) - (double) (V - 1) /(double)k - averageDiGammas +
- (double) (V - 1) * MathsUtils.digamma(N);
- miComputed = true;
- return mi;
- }
-
- @Override
- public double[] computeLocalOfPreviousObservations() throws Exception {
- double[] localMi = new double[N];
- if (V == 1) {
- miComputed = true;
- return localMi;
- }
-
// Constants:
- double digammaK = MathsUtils.digamma(k);
- double Vminus1TimesDigammaN = (double) (V - 1) * MathsUtils.digamma(N);
- double Vminus1TimesInvK = (double) (V - 1) / (double)k;
+ double dimensionsMinus1DivK = (double) (dimensions - 1) / (double)k;
+ double dimensionsMinus1TimesDiGammaN = (double) (dimensions - 1) * digammaN;
- // Count the average number of points within eps_x[v] for each marginal v
- double averageDiGammas = 0;
- double[] avNx = new double[V];
-
- for (int t = 0; t < N; t++) {
- // Compute eps_x (for each marginal) for this time step:
- // First get the marginal norms to all neighbours
- // (note that norm of point t to itself will be set to infinity).
-
- double[][] normsForT = EuclideanUtils.computeNorms(data, t);
- double[][] jointNorm = new double[N][2];
- for (int t2 = 0; t2 < N; t2++) {
- jointNorm[t2][JOINT_NORM_VAL_COLUMN] = MatrixUtils.max(normsForT[t2]);
- // And store the time step for back reference after the
- // array is sorted.
- jointNorm[t2][JOINT_NORM_TIMESTEP_COLUMN] = t2;
- }
+ // Count the average number of points within eps_x for each marginal x of each point
+ double sumDiGammas = 0;
+ double[] sumNMarginals = new double[dimensions];
+
+ for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
+ // Compute eps_x for each marginal x for this time step by
+ // finding the kth closest neighbours for point t:
+ PriorityQueue
@@ -58,57 +63,22 @@ import java.util.Vector;
*
* @author Joseph Lizier (email,
* www)
+ * @author Ipek Özdemir
*/
-public abstract class MultiInfoCalculatorKraskov implements
- MultiInfoCalculator {
+public abstract class MultiInfoCalculatorKraskov
+ extends MultiInfoCalculatorCommon
+ implements Cloneable { // See comments on clonability above
/**
* we compute distances to the kth nearest neighbour
*/
protected int k = 4;
+
/**
- * Cached observations
+ * The norm type in use (see {@link #PROP_NORM_TYPE})
*/
- protected double[][] data;
- /**
- * Whether we are in debug mode
- */
- protected boolean debug;
- /**
- * Last average multi-info computed
- */
- protected double mi;
- protected boolean miComputed;
-
- /**
- * Set of individually supplied observations
- */
- private Vector
+ *
+ *
+ * @param returnLocals whether to return an array or local values, or else
+ * sums of these values
+ * @return either the average multi-info, or array of local multi-info value,
+ * in nats not bits
+ * @throws Exception
+ */
+ protected double[] computeFromObservations(boolean returnLocals) throws Exception {
+
+ double[] returnValues = null;
+
+ // 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.
+ double[][][] separateMarginals = null;
+ int[] dimensionsArray = null;
+ if (kdTreeJoint == null) {
+ // We need to pull out 2D time series (of only one variable)
+ // for each marginal variable here
+ separateMarginals = new double[dimensions][][];
+ dimensionsArray = new int[dimensions];
+ for (int d = 0; d < dimensions; d++) {
+ separateMarginals[d] =
+ MatrixUtils.selectColumns(observations, new int[] {d});
+ dimensionsArray[d] = 1;
+ }
+ kdTreeJoint = new KdTree(dimensionsArray, separateMarginals);
+ kdTreeJoint.setNormType(normType);
+ }
+ if (rangeSearchersInMarginals == null) {
+ rangeSearchersInMarginals = new UnivariateNearestNeighbourSearcher[dimensions];
+ for (int d = 0; d < dimensions; d++) {
+ rangeSearchersInMarginals[d] =
+ new UnivariateNearestNeighbourSearcher(
+ MatrixUtils.selectColumn(observations, d));
+ rangeSearchersInMarginals[d].setNormType(normType);
+ }
+ }
+
+ if (numThreads == 1) {
+ // Single-threaded implementation:
+ returnValues = partialComputeFromObservations(0, totalObservations, returnLocals);
+
+ } else {
+ // We're going multithreaded:
+ if (returnLocals) {
+ // We're computing local MI
+ returnValues = new double[totalObservations];
+ } else {
+ // We're computing average MI
+ returnValues = new double[1 + dimensions];
+ }
+
+ // 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];
+ MultiInfoKraskovThreadRunner[] runners = new MultiInfoKraskovThreadRunner[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 MultiInfoKraskovThreadRunner(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 multi-info; 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 MI, keep the running sums of digammas and counts
+ MatrixUtils.addInPlace(returnValues, runners[t].getReturnValues());
+ }
+ }
+ }
+
+ // Finalise the results:
+ if (returnLocals) {
+ return returnValues;
+ } else {
+ // Compute the average number of points within eps_x and eps_y
+ double averageDiGammas = returnValues[MultiInfoKraskovThreadRunner.INDEX_SUM_DIGAMMAS] / (double) totalObservations;
+ double[] avNMarginals = new double[dimensions];
+ for (int d = 0; d < dimensions; d++) {
+ avNMarginals[d] = returnValues[1 + d] / (double) totalObservations;
+ if (debug) {
+ System.out.printf("Average n_%d=%.3f, ", d, avNMarginals[d]);
+ }
+ }
+ if (debug) {
+ System.out.println();
+ }
+
+ // Finalise the average result, depending on which algorithm we are implementing:
+ if (isAlgorithm1) {
+ return new double[] { digammaK - averageDiGammas + (double) (dimensions - 1) * digammaN };
+ } else {
+ return new double[] { digammaK - ((double) (dimensions - 1) / (double)k) - averageDiGammas +
+ (double) (dimensions - 1) * digammaN };
+ }
+ }
}
/**
- * Utility function used for debugging, printing digamma constants
+ * 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.
*
- * @param N
- * @return
+ *
+ *
+ *
+ * @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
*/
- public abstract String printConstants(int N) throws Exception;
+ protected abstract double[] partialComputeFromObservations(
+ int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception;
/**
- * Utility to take a reordering matrix and return the array of reordered time indices from
- * which to find the reordered data to be inserted at timeStep.
+ * Private class to handle multi-threading of the Kraskov algorithms.
+ * Each instance calls partialComputeFromObservations()
+ * to compute nearest neighbours for a part of the data.
*
- * @param reordering the specific new orderings to use. First index is the variable number
- * (can be for all variables, or one less than all if the first is not to be reordered),
- * second index is the time step, the value is the reordered time step to use
- * for that variable at the given time step.
- * If null, no reordering is performed.
- * @param timeStep
- * @return array of reordered time indices from
- * which to find the reordered data to be inserted at timeStep
+ *
+ * @author Joseph Lizier (email,
+ * www)
+ * @author Ipek Özdemir
*/
- protected int[] reorderedTimeStepsForEachMarginal(int[][] reordering, int timeStep) {
- // Create storage for the reordered time steps for the variables
- int[] tForEachMarginal = new int[V];
- if (reordering == null) {
- // We're not reordering
- for (int v = 0; v < V; v++) {
- tForEachMarginal[v] = timeStep;
+ private class MultiInfoKraskovThreadRunner implements Runnable {
+ protected MultiInfoCalculatorKraskov miCalc;
+ 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 MultiInfoKraskovThreadRunner(
+ MultiInfoCalculatorKraskov miCalc,
+ int myStartTimePoint, int numberOfTimePoints,
+ boolean computeLocals) {
+ this.miCalc = miCalc;
+ 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;
}
- } else {
- boolean reorderingFirstColumn = (reordering.length == V);
- int reorderIndex = 0;
- // Handle the first column
- if (reorderingFirstColumn) {
- tForEachMarginal[0] = reordering[reorderIndex++][timeStep];
- } else {
- tForEachMarginal[0] = timeStep;
- }
- // Handle subsequent columns
- for (int v = 1; v < V; v++) {
- tForEachMarginal[v] = reordering[reorderIndex++][timeStep];
+ return returnValues;
+ }
+
+ /**
+ * Start the thread for the given parameters
+ */
+ public void run() {
+ try {
+ returnValues = miCalc.partialComputeFromObservations(
+ myStartTimePoint, numberOfTimePoints, computeLocals);
+ } catch (Exception e) {
+ // Store the exception for later retrieval
+ problem = e;
+ return;
}
}
- return tForEachMarginal;
}
+ // end class MultiInfoKraskovThreadRunner
}
diff --git a/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov1.java b/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov1.java
index 81aa300..3483ca7 100755
--- a/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov1.java
+++ b/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov1.java
@@ -18,10 +18,12 @@
package infodynamics.measures.continuous.kraskov;
+import java.util.Calendar;
+import java.util.PriorityQueue;
+
import infodynamics.measures.continuous.MultiInfoCalculator;
-import infodynamics.utils.EuclideanUtils;
import infodynamics.utils.MathsUtils;
-import infodynamics.utils.MatrixUtils;
+import infodynamics.utils.NeighbourNodeData;
/**
*