mirror of https://github.com/jlizier/jidt
Added auto-embedding with Ragwitz criteria to AIS Kraskov calculator, partially addressing Issue 38. Also adds getProperty() method to AIS calculators and MI calculators.
This commit is contained in:
joseph.lizier
parent
0f5bf692ac
commit
687978103c
|
|
@ -150,6 +150,21 @@ public interface ActiveInfoStorageCalculator {
|
|||
*/
|
||||
public void setProperty(String propertyName, String propertyValue) throws Exception;
|
||||
|
||||
/**
|
||||
* Get current property values for the calculator.
|
||||
*
|
||||
* <p>Valid property names, and what their
|
||||
* values should represent, are the same as those for
|
||||
* {@link #setProperty(String, String)}</p>
|
||||
*
|
||||
* <p>Unknown property values are responded to with a null return value.</p>
|
||||
*
|
||||
* @param propertyName name of the property
|
||||
* @return current value of the property
|
||||
* @throws Exception for invalid property values
|
||||
*/
|
||||
public String getProperty(String propertyName) throws Exception;
|
||||
|
||||
/**
|
||||
* Sets a single time-series from which to compute the PDF for the AIS.
|
||||
* Cannot be called in conjunction with other methods for setting/adding
|
||||
|
|
|
|||
|
|
@ -78,6 +78,12 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
*/
|
||||
protected boolean debug = false;
|
||||
|
||||
/**
|
||||
* Storage for source observations supplied via {@link #addObservations(double[])}
|
||||
* type calls
|
||||
*/
|
||||
protected Vector<double[]> vectorOfObservationTimeSeries;
|
||||
|
||||
/**
|
||||
* Construct using an instantiation of the named MI calculator
|
||||
*
|
||||
|
|
@ -158,7 +164,6 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
public void initialise(int k, int tau) throws Exception {
|
||||
this.k = k;
|
||||
this.tau = tau;
|
||||
miCalc.initialise(k, 1);
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -210,20 +215,31 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
}
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getProperty(String propertyName)
|
||||
throws Exception {
|
||||
|
||||
if (propertyName.equalsIgnoreCase(K_PROP_NAME)) {
|
||||
return Integer.toString(k);
|
||||
} else if (propertyName.equalsIgnoreCase(TAU_PROP_NAME)) {
|
||||
return Integer.toString(tau);
|
||||
} else {
|
||||
// No property was set on this class, assume it is for the underlying
|
||||
// MI calculator, even if it is for
|
||||
// MutualInfoCalculatorMultiVariate.PROP_TIME_DIFF which
|
||||
// is not a valid property for the AIS calculator:
|
||||
return miCalc.getProperty(propertyName);
|
||||
}
|
||||
}
|
||||
|
||||
/* (non-Javadoc)
|
||||
* @see infodynamics.measures.continuous.ActiveInfoStorageCalculator#setObservations(double[])
|
||||
*/
|
||||
@Override
|
||||
public void setObservations(double[] observations) throws Exception {
|
||||
if (observations.length - (k-1)*tau - 1 <= 0) {
|
||||
// There are no observations to add here
|
||||
throw new Exception("Not enough observations to set here given k and tau");
|
||||
}
|
||||
double[][] currentDestPastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(observations, k, tau, (k-1)*tau, observations.length - (k-1)*tau - 1);
|
||||
double[][] currentDestNextVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(observations, 1, (k-1)*tau + 1, observations.length - (k-1)*tau - 1);
|
||||
miCalc.setObservations(currentDestPastVectors, currentDestNextVectors);
|
||||
startAddObservations();
|
||||
addObservations(observations);
|
||||
finaliseAddObservations();
|
||||
}
|
||||
|
||||
/* (non-Javadoc)
|
||||
|
|
@ -231,7 +247,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
*/
|
||||
@Override
|
||||
public void startAddObservations() {
|
||||
miCalc.startAddObservations();
|
||||
vectorOfObservationTimeSeries = new Vector<double[]>();
|
||||
}
|
||||
|
||||
/* (non-Javadoc)
|
||||
|
|
@ -239,6 +255,20 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
*/
|
||||
@Override
|
||||
public void addObservations(double[] observations) throws Exception {
|
||||
// Store these observations in our vector for now
|
||||
vectorOfObservationTimeSeries.add(observations);
|
||||
}
|
||||
|
||||
/**
|
||||
* Protected method to internally parse and submit observations through
|
||||
* to the underlying MI calculator once any internal parameter settings
|
||||
* have been finalised (in the case of automatically determining the embedding
|
||||
* parameters)
|
||||
*
|
||||
* @param observations time series of observations
|
||||
* @throws Exception
|
||||
*/
|
||||
protected void addObservationsAfterParamsDetermined(double[] observations) throws Exception {
|
||||
if (observations.length - (k-1)*tau - 1 <= 0) {
|
||||
// There are no observations to add here
|
||||
// Don't throw an exception, do nothing since more observations
|
||||
|
|
@ -251,7 +281,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
MatrixUtils.makeDelayEmbeddingVector(observations, 1, (k-1)*tau + 1, observations.length - (k-1)*tau - 1);
|
||||
miCalc.addObservations(currentDestPastVectors, currentDestNextVectors);
|
||||
}
|
||||
|
||||
|
||||
/* (non-Javadoc)
|
||||
* @see infodynamics.measures.continuous.ActiveInfoStorageCalculator#addObservations(double[], int, int)
|
||||
*/
|
||||
|
|
@ -261,11 +291,37 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
addObservations(MatrixUtils.select(observations, startTime, numTimeSteps));
|
||||
}
|
||||
|
||||
/**
|
||||
* Hook in case child implementations need to perform any processing on the
|
||||
* observation time series prior to their being processed and supplied
|
||||
* to the underlying MI calculator.
|
||||
* Primarily this is to allow the child implementation to automatically determine
|
||||
* embedding parameters if desired.
|
||||
* Child implementations do not need to override this default empty implementation
|
||||
* if no new functionality is required.
|
||||
*/
|
||||
public void preFinaliseAddObservations() throws Exception {
|
||||
// Empty implementation supplied by default.
|
||||
}
|
||||
|
||||
/* (non-Javadoc)
|
||||
* @see infodynamics.measures.continuous.ActiveInfoStorageCalculator#finaliseAddObservations()
|
||||
*/
|
||||
@Override
|
||||
public void finaliseAddObservations() throws Exception {
|
||||
// Auto embed if required
|
||||
preFinaliseAddObservations();
|
||||
|
||||
// Initialise the MI calculator, including any auto-embedding length
|
||||
miCalc.initialise(k, 1);
|
||||
miCalc.startAddObservations();
|
||||
// Send all of the observations through:
|
||||
for (double[] observations : vectorOfObservationTimeSeries) {
|
||||
addObservationsAfterParamsDetermined(observations);
|
||||
}
|
||||
vectorOfObservationTimeSeries = null; // No longer required
|
||||
|
||||
// TODO do we need to throw an exception if there are no observations to add?
|
||||
miCalc.finaliseAddObservations();
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -102,4 +102,19 @@ public interface MutualInfoCalculatorMultiVariate extends ChannelCalculatorMulti
|
|||
* @throws Exception
|
||||
*/
|
||||
public double computeAverageLocalOfObservations(int[] newOrdering) throws Exception;
|
||||
|
||||
/**
|
||||
* Get current the value for a given property for the calculator.
|
||||
*
|
||||
* <p>Valid property names, and what their
|
||||
* values should represent, are the same as those for
|
||||
* {@link #setProperty(String, String)}</p>
|
||||
*
|
||||
* <p>Unknown property values are responded to with a null return value.</p>
|
||||
*
|
||||
* @param propertyName name of the property
|
||||
* @return current value of the property
|
||||
* @throws Exception for invalid property values
|
||||
*/
|
||||
public String getProperty(String propertyName) throws Exception;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -175,6 +175,18 @@ public abstract class MutualInfoMultiVariateCommon implements
|
|||
}
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getProperty(String propertyName)
|
||||
throws Exception {
|
||||
|
||||
if (propertyName.equalsIgnoreCase(PROP_TIME_DIFF)) {
|
||||
return Integer.toString(timeDiff);
|
||||
} else {
|
||||
// No property was recognised here
|
||||
return null;
|
||||
}
|
||||
}
|
||||
|
||||
public void setObservations(double[][] source, double[][] destination) throws Exception {
|
||||
startAddObservations();
|
||||
addObservations(source, destination);
|
||||
|
|
|
|||
|
|
@ -688,6 +688,25 @@ public class MutualInfoCalculatorMultiVariateKernel
|
|||
}
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getProperty(String propertyName)
|
||||
throws Exception {
|
||||
|
||||
if (propertyName.equalsIgnoreCase(KERNEL_WIDTH_PROP_NAME) ||
|
||||
propertyName.equalsIgnoreCase(EPSILON_PROP_NAME)) {
|
||||
return Double.toString(kernelWidth);
|
||||
} else if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
|
||||
return Boolean.toString(normalise);
|
||||
} else if (propertyName.equalsIgnoreCase(DYN_CORR_EXCL_TIME_NAME)) {
|
||||
return Integer.toString(dynCorrExclTime);
|
||||
} else if (propertyName.equalsIgnoreCase(FORCE_KERNEL_COMPARE_TO_ALL)) {
|
||||
return Boolean.toString(forceCompareToAll);
|
||||
} else {
|
||||
// try the superclass:
|
||||
return super.getProperty(propertyName);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
*
|
||||
* @return the kernel width in use in the calculator
|
||||
|
|
|
|||
|
|
@ -21,6 +21,7 @@ package infodynamics.measures.continuous.kraskov;
|
|||
import infodynamics.measures.continuous.ActiveInfoStorageCalculator;
|
||||
import infodynamics.measures.continuous.ActiveInfoStorageCalculatorViaMutualInfo;
|
||||
import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
|
||||
import infodynamics.utils.MatrixUtils;
|
||||
|
||||
/**
|
||||
* An Active Information Storage (AIS) calculator (implementing {@link ActiveInfoStorageCalculator})
|
||||
|
|
@ -43,7 +44,10 @@ import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
|
|||
* <li>{@link #setProperty(String, String)} allowing properties for
|
||||
* {@link MutualInfoCalculatorMultiVariateKraskov#setProperty(String, String)}
|
||||
* (except {@link MutualInfoCalculatorMultiVariate#PROP_TIME_DIFF} as outlined
|
||||
* in {@link ActiveInfoStorageCalculatorViaMutualInfo#setProperty(String, String)}).</li>
|
||||
* in {@link ActiveInfoStorageCalculatorViaMutualInfo#setProperty(String, String)}).
|
||||
* Embedding parameters may be automatically determined as per the Ragwitz criteria
|
||||
* by setting the property {@link #PROP_AUTO_EMBED_METHOD} to {@link #AUTO_EMBED_METHOD_RAGWITZ}
|
||||
* (plus additional parameter settings for this).</li>
|
||||
* <li>Computed values are in <b>nats</b>, not bits!</li>
|
||||
* </ul>
|
||||
* </p>
|
||||
|
|
@ -54,6 +58,8 @@ import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
|
|||
* <a href="http://dx.doi.org/10.1016/j.ins.2012.04.016">
|
||||
* "Local measures of information storage in complex distributed computation"</a>,
|
||||
* Information Sciences, vol. 208, pp. 39-54, 2012.</li>
|
||||
* <li>Ragwitz and Kantz, "Markov models from data by simple nonlinear time series
|
||||
* predictors in delay embedding spaces", Physical Review E, vol 65, 056201 (2002).</li>
|
||||
* </ul>
|
||||
*
|
||||
* @author Joseph Lizier (<a href="joseph.lizier at gmail.com">email</a>,
|
||||
|
|
@ -74,7 +80,63 @@ public class ActiveInfoStorageCalculatorKraskov
|
|||
* Class name for KSG MI estimator via KSG algorithm 2
|
||||
*/
|
||||
public static final String MI_CALCULATOR_KRASKOV2 = MutualInfoCalculatorMultiVariateKraskov2.class.getName();
|
||||
|
||||
|
||||
/**
|
||||
* Property name for the auto-embedding method. Defaults to {@link #AUTO_EMBED_METHOD_NONE}
|
||||
*/
|
||||
public static final String PROP_AUTO_EMBED_METHOD = "AUTO_EMBED_METHOD";
|
||||
/**
|
||||
* Valid value for the property {@link #PROP_AUTO_EMBED_METHOD} indicating that
|
||||
* no auto embedding should be done (i.e. to use manually supplied parameters)
|
||||
*/
|
||||
public static final String AUTO_EMBED_METHOD_NONE = "NONE";
|
||||
/**
|
||||
* Valid value for the property {@link #PROP_AUTO_EMBED_METHOD} indicating that
|
||||
* the Ragwitz optimisation technique should be used for automatic embedding
|
||||
*/
|
||||
public static final String AUTO_EMBED_METHOD_RAGWITZ = "RAGWITZ";
|
||||
/**
|
||||
* Internal variable tracking what type of auto embedding (if any)
|
||||
* we are using
|
||||
*/
|
||||
protected String autoEmbeddingMethod = AUTO_EMBED_METHOD_NONE;
|
||||
|
||||
/**
|
||||
* Property name for maximum k (embedding length) for the auto-embedding search. Default to 1
|
||||
*/
|
||||
public static final String PROP_K_SEARCH_MAX = "AUTO_EMBED_K_SEARCH_MAX";
|
||||
/**
|
||||
* Internal variable for storing the maximum embedding length to search up to for
|
||||
* automating the parameters.
|
||||
*/
|
||||
protected int k_search_max = 1;
|
||||
|
||||
/**
|
||||
* Property name for maximum tau (embedding delay) for the auto-embedding search. Default to 1
|
||||
*/
|
||||
public static final String PROP_TAU_SEARCH_MAX = "AUTO_EMBED_TAU_SEARCH_MAX";
|
||||
/**
|
||||
* Internal variable for storing the maximum embedding delay to search up to for
|
||||
* automating the parameters.
|
||||
*/
|
||||
protected int tau_search_max = 1;
|
||||
|
||||
/**
|
||||
* Property name for the number of nearest neighbours to use for the auto-embedding search (Ragwitz criteria).
|
||||
* Defaults to match the value in use for {@link MutualInfoCalculatorMultiVariateKraskov#PROP_K}
|
||||
*/
|
||||
public static final String PROP_RAGWITZ_NUM_NNS = "AUTO_EMBED_RAGWITZ_NUM_NNS";
|
||||
/**
|
||||
* Internal variable for storing the number of nearest neighbours to use for the
|
||||
* auto embedding search (Ragwitz criteria)
|
||||
*/
|
||||
protected int ragwitz_num_nns = 1;
|
||||
/**
|
||||
* Internal variable to track whether the property {@link #PROP_RAGWITZ_NUM_NNS} has been
|
||||
* set yet
|
||||
*/
|
||||
protected boolean ragwitz_num_nns_set = false;
|
||||
|
||||
/**
|
||||
* Creates a new instance of the Kraskov-Stoegbauer-Grassberger style AIS calculator.
|
||||
*
|
||||
|
|
@ -126,4 +188,181 @@ public class ActiveInfoStorageCalculatorKraskov
|
|||
throw new ClassNotFoundException("Algorithm must be 1 or 2");
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Sets properties for the AIS calculator.
|
||||
* New property values are not guaranteed to take effect until the next call
|
||||
* to an initialise method.
|
||||
*
|
||||
* <p>Valid property names, and what their
|
||||
* values should represent, include:</p>
|
||||
* <ul>
|
||||
* <li>{@link #PROP_AUTO_EMBED_METHOD} -- method by which the calculator
|
||||
* automatically determines the embedding history length ({@link #K_PROP_NAME})
|
||||
* and embedding delay ({@link #TAU_PROP_NAME}). Default is {@link #AUTO_EMBED_METHOD_NONE} meaning
|
||||
* values are set manually; other accepted values include: {@link #AUTO_EMBED_METHOD_RAGWITZ} for use
|
||||
* of the Ragwitz criteria (searching up to {@link #PROP_K_SEARCH_MAX} and
|
||||
* {@link #PROP_TAU_SEARCH_MAX})</li>
|
||||
* <li>{@link #PROP_K_SEARCH_MAX} -- maximum embedded history length to search
|
||||
* up to if automatically determining the embedding parameters (as set by
|
||||
* {@link #PROP_AUTO_EMBED_METHOD}); default is 1</li>
|
||||
* <li>{@link #PROP_TAU_SEARCH_MAX} -- maximum embedded history length to search
|
||||
* up to if automatically determining the embedding parameters (as set by
|
||||
* {@link #PROP_AUTO_EMBED_METHOD}); default is 1</li>
|
||||
* <li>{@link #PROP_RAGWITZ_NUM_NNS} -- number of nearest neighbours to use
|
||||
* in the auto-embedding if the property {@link #PROP_AUTO_EMBED_METHOD}
|
||||
* has been set to {@link #AUTO_EMBED_METHOD_RAGWITZ}. Defaults to the property value
|
||||
* set for {@link MutualInfoCalculatorMultiVariateKraskov.PROP_K}</li>
|
||||
* <li>Any properties accepted by {@link super#setProperty(String, String)}</li>
|
||||
* <li>Or properties accepted by the underlying
|
||||
* {@link MutualInfoCalculatorMultiVariateKraskov#setProperty(String, String)} implementation.</li>
|
||||
* </ul>
|
||||
* <p>One should set {@link MutualInfoCalculatorMultiVariateKraskov#PROP_K} here, the number
|
||||
* of neighbouring points one should count up to in determining the joint kernel size.</p>
|
||||
*
|
||||
* @param propertyName name of the property
|
||||
* @param propertyValue value of the property.
|
||||
* @throws Exception if there is a problem with the supplied value).
|
||||
*/
|
||||
public void setProperty(String propertyName, String propertyValue)
|
||||
throws Exception {
|
||||
boolean propertySet = true;
|
||||
if (propertyName.equalsIgnoreCase(PROP_AUTO_EMBED_METHOD)) {
|
||||
// New method set for determining the embedding parameters
|
||||
autoEmbeddingMethod = propertyValue;
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_K_SEARCH_MAX)) {
|
||||
// Set max embedding history length for auto determination of embedding
|
||||
k_search_max = Integer.parseInt(propertyValue);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_TAU_SEARCH_MAX)) {
|
||||
// Set maximum embedding delay for auto determination of embedding
|
||||
tau_search_max = Integer.parseInt(propertyValue);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_RAGWITZ_NUM_NNS)) {
|
||||
// Set the number of nearest neighbours to use in case of Ragwitz auto embedding:
|
||||
ragwitz_num_nns = Integer.parseInt(propertyValue);
|
||||
ragwitz_num_nns_set = true;
|
||||
} else {
|
||||
propertySet = false;
|
||||
// Assume it was a property for the parent class or underlying MI calculator
|
||||
super.setProperty(propertyName, propertyValue);
|
||||
}
|
||||
if (debug && propertySet) {
|
||||
System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
|
||||
" to " + propertyValue);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Get property values for the calculator.
|
||||
*
|
||||
* <p>Valid property names, and what their
|
||||
* values should represent, are the same as those for
|
||||
* {@link #setProperty(String, String)}</p>
|
||||
*
|
||||
* <p>Unknown property values are responded to with a null return value.</p>
|
||||
*
|
||||
* @param propertyName name of the property
|
||||
* @return current value of the property
|
||||
* @throws Exception for invalid property values
|
||||
*/
|
||||
public String getProperty(String propertyName)
|
||||
throws Exception {
|
||||
|
||||
if (propertyName.equalsIgnoreCase(PROP_AUTO_EMBED_METHOD)) {
|
||||
return autoEmbeddingMethod;
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_K_SEARCH_MAX)) {
|
||||
return Integer.toString(k_search_max);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_TAU_SEARCH_MAX)) {
|
||||
return Integer.toString(tau_search_max);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_RAGWITZ_NUM_NNS)) {
|
||||
if (ragwitz_num_nns_set) {
|
||||
return Integer.toString(ragwitz_num_nns);
|
||||
} else {
|
||||
return miCalc.getProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_K);
|
||||
}
|
||||
} else {
|
||||
// Assume it was a property for the parent class or underlying MI calculator
|
||||
return super.getProperty(propertyName);
|
||||
}
|
||||
}
|
||||
|
||||
@Override
|
||||
public void preFinaliseAddObservations() throws Exception {
|
||||
// Automatically determine the embedding parameters for the given time series
|
||||
|
||||
if (autoEmbeddingMethod.equalsIgnoreCase(AUTO_EMBED_METHOD_NONE)) {
|
||||
return;
|
||||
}
|
||||
// Else we need to auto embed
|
||||
|
||||
double bestPredictionError = Double.POSITIVE_INFINITY;
|
||||
int k_candidate_best = 1;
|
||||
int tau_candidate_best = 1;
|
||||
|
||||
if (autoEmbeddingMethod.equalsIgnoreCase(AUTO_EMBED_METHOD_RAGWITZ)) {
|
||||
if (debug) {
|
||||
System.out.printf("Beginning Ragwitz auto-embedding with k_max=%d, tau_max=%d\n",
|
||||
k_search_max, tau_search_max);
|
||||
}
|
||||
|
||||
for (int k_candidate = 1; k_candidate <= k_search_max; k_candidate++) {
|
||||
for (int tau_candidate = 1; tau_candidate <= tau_search_max; tau_candidate++) {
|
||||
// Use our internal MI calculator in case it has any particular
|
||||
// properties we need to have been set already
|
||||
miCalc.initialise(k_candidate, 1);
|
||||
miCalc.startAddObservations();
|
||||
// Send all of the observations through:
|
||||
for (double[] observations : vectorOfObservationTimeSeries) {
|
||||
double[][] currentDestPastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(observations, k_candidate,
|
||||
tau_candidate, (k_candidate-1)*tau_candidate,
|
||||
observations.length - (k_candidate-1)*tau_candidate - 1);
|
||||
double[][] currentDestNextVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(observations, 1,
|
||||
(k_candidate-1)*tau_candidate + 1,
|
||||
observations.length - (k_candidate-1)*tau_candidate - 1);
|
||||
miCalc.addObservations(currentDestPastVectors, currentDestNextVectors);
|
||||
}
|
||||
miCalc.finaliseAddObservations();
|
||||
// Now grab the prediction errors of the next value from the required number of
|
||||
// nearest neighbours of the previous state: (array is of only one term)
|
||||
double[] predictionError;
|
||||
if (ragwitz_num_nns_set) {
|
||||
predictionError =
|
||||
((MutualInfoCalculatorMultiVariateKraskov) miCalc).
|
||||
computePredictionErrorsFromObservations(false, ragwitz_num_nns);
|
||||
} else {
|
||||
predictionError =
|
||||
((MutualInfoCalculatorMultiVariateKraskov) miCalc).
|
||||
computePredictionErrorsFromObservations(false);
|
||||
}
|
||||
if (debug) {
|
||||
System.out.printf("Embedding prediction error (dim=%d) for k=%d,tau=%d is %.3f\n",
|
||||
predictionError.length, k_candidate, tau_candidate,
|
||||
predictionError[0] / (double) miCalc.getNumObservations());
|
||||
}
|
||||
if ((predictionError[0] / (double) miCalc.getNumObservations())
|
||||
< bestPredictionError) {
|
||||
// This parameter setting is the best so far:
|
||||
// (Note division by number of observations to normalise
|
||||
// for less observations for larger k and tau)
|
||||
bestPredictionError = predictionError[0] / (double) miCalc.getNumObservations();
|
||||
k_candidate_best = k_candidate;
|
||||
tau_candidate_best = tau_candidate;
|
||||
}
|
||||
if (k_candidate == 1) {
|
||||
// tau is irrelevant, so no point testing other values
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Make sure the embedding length and delay are set here
|
||||
k = k_candidate_best;
|
||||
tau = tau_candidate_best;
|
||||
if (debug) {
|
||||
System.out.printf("Embedding parameters set to k=%d,tau=%d (for prediction error %.3f)\n",
|
||||
k, tau, bestPredictionError);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@
|
|||
package infodynamics.measures.continuous.kraskov;
|
||||
|
||||
import java.util.Calendar;
|
||||
import java.util.PriorityQueue;
|
||||
import java.util.Random;
|
||||
|
||||
import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
|
||||
|
|
@ -28,6 +29,7 @@ import infodynamics.utils.KdTree;
|
|||
import infodynamics.utils.MathsUtils;
|
||||
import infodynamics.utils.MatrixUtils;
|
||||
import infodynamics.utils.NearestNeighbourSearcher;
|
||||
import infodynamics.utils.NeighbourNodeData;
|
||||
|
||||
/**
|
||||
* <p>Computes the differential mutual information of two given multivariate sets of
|
||||
|
|
@ -256,6 +258,40 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
|
|||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Get property values for the calculator.
|
||||
*
|
||||
* <p>Valid property names, and what their
|
||||
* values should represent, are the same as those for
|
||||
* {@link #setProperty(String, String)}</p>
|
||||
*
|
||||
* <p>Unknown property values are responded to with a null return value.</p>
|
||||
*
|
||||
* @param propertyName name of the property
|
||||
* @return current value of the property
|
||||
* @throws Exception for invalid property values
|
||||
*/
|
||||
public String getProperty(String propertyName)
|
||||
throws Exception {
|
||||
|
||||
if (propertyName.equalsIgnoreCase(PROP_K)) {
|
||||
return Integer.toString(k);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_NORM_TYPE)) {
|
||||
return KdTree.convertNormTypeToString(normType);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) {
|
||||
return Boolean.toString(normalise);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_DYN_CORR_EXCL_TIME)) {
|
||||
return Integer.toString(dynCorrExclTime);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
|
||||
return Double.toString(noiseLevel);
|
||||
} else if (propertyName.equalsIgnoreCase(PROP_NUM_THREADS)) {
|
||||
return Integer.toString(numThreads);
|
||||
} else {
|
||||
// try the superclass:
|
||||
return super.getProperty(propertyName);
|
||||
}
|
||||
}
|
||||
|
||||
/* (non-Javadoc)
|
||||
* @see infodynamics.measures.continuous.MutualInfoMultiVariateCommon#finaliseAddObservations()
|
||||
*/
|
||||
|
|
@ -415,25 +451,7 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
|
|||
|
||||
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.
|
||||
if (kdTreeJoint == null) {
|
||||
kdTreeJoint = new KdTree(new int[] {dimensionsSource, dimensionsDest},
|
||||
new double[][][] {sourceObservations, destObservations});
|
||||
kdTreeJoint.setNormType(normType);
|
||||
}
|
||||
if (nnSearcherSource == null) {
|
||||
nnSearcherSource = NearestNeighbourSearcher.create(sourceObservations);
|
||||
nnSearcherSource.setNormType(normType);
|
||||
}
|
||||
if (nnSearcherDest == null) {
|
||||
nnSearcherDest = NearestNeighbourSearcher.create(destObservations);
|
||||
nnSearcherDest.setNormType(normType);
|
||||
}
|
||||
ensureKdTreesConstructed();
|
||||
|
||||
if (numThreads == 1) {
|
||||
// Single-threaded implementation:
|
||||
|
|
@ -534,6 +552,119 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
|
|||
protected abstract double[] partialComputeFromObservations(
|
||||
int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception;
|
||||
|
||||
/**
|
||||
* 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(new int[] {dimensionsSource, dimensionsDest},
|
||||
new double[][][] {sourceObservations, destObservations});
|
||||
kdTreeJoint.setNormType(normType);
|
||||
}
|
||||
if (nnSearcherSource == null) {
|
||||
nnSearcherSource = NearestNeighbourSearcher.create(sourceObservations);
|
||||
nnSearcherSource.setNormType(normType);
|
||||
}
|
||||
if (nnSearcherDest == null) {
|
||||
nnSearcherDest = NearestNeighbourSearcher.create(destObservations);
|
||||
nnSearcherDest.setNormType(normType);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the prediction error in one variable from the k nearest neighbours (kNNs) of the observation
|
||||
* of the other variable.
|
||||
* The kNNs of the variable to predict from are formed from the supplied norm type within that variable.
|
||||
* The number of kNNs to use here is the current property value set for {@link #PROP_K}.
|
||||
* The prediction error is a sum of absolute errors for each dimension within the variable to predict
|
||||
*
|
||||
* @param predictFirstVariable true for predicting the first variable (Source) or
|
||||
* false for predicting the second variable (destination)
|
||||
* @return array of prediction errors for each dimension of the predicted variable
|
||||
* @throws Exception
|
||||
*/
|
||||
public double[] computePredictionErrorsFromObservations(boolean predictFirstVariable) throws Exception {
|
||||
return computePredictionErrorsFromObservations(predictFirstVariable, k);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the prediction error in one variable from the k nearest neighbours (kNNs) of the observation
|
||||
* of the other variable.
|
||||
* The kNNs of the variable to predict from are formed from the supplied norm type within that variable.
|
||||
* The prediction error is a sum of absolute errors for each dimension within the variable to predict
|
||||
*
|
||||
* @param predictFirstVariable true for predicting the first variable (Source) or
|
||||
* false for predicting the second variable (destination)
|
||||
* @param kNNs number of nearest neighbours to use
|
||||
* @return array of prediction errors for each dimension of the predicted variable
|
||||
* @throws Exception
|
||||
*/
|
||||
public double[] computePredictionErrorsFromObservations(boolean predictFirstVariable, int kNNs) throws Exception {
|
||||
int N = sourceObservations.length; // number of observations
|
||||
|
||||
double[] totalErrors = null;
|
||||
|
||||
ensureKdTreesConstructed();
|
||||
|
||||
if (numThreads == 1) {
|
||||
// Single-threaded implementation:
|
||||
totalErrors = partialComputePredictionErrorFromObservations(0, N, kNNs, predictFirstVariable);
|
||||
} else {
|
||||
// We're going multithreaded:
|
||||
totalErrors = new double[predictFirstVariable ? dimensionsSource : dimensionsDest];
|
||||
|
||||
// Distribute the observations to the threads for the parallel processing
|
||||
int lTimesteps = N / numThreads; // each thread gets the same amount of data
|
||||
int res = N % numThreads; // the first thread gets the residual data
|
||||
if (debug) {
|
||||
System.out.printf("Computing prediction errors for variable %d from variable %d with %d threads (%d timesteps each, plus %d residual)\n",
|
||||
predictFirstVariable ? 1 : 2, predictFirstVariable ? 2 : 1,
|
||||
numThreads, lTimesteps, res);
|
||||
}
|
||||
Thread[] tCalculators = new Thread[numThreads];
|
||||
MiKraskovPredictionThreadRunner[] runners = new MiKraskovPredictionThreadRunner[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 MiKraskovPredictionThreadRunner(this, startTime,
|
||||
numTimesteps, kNNs, predictFirstVariable);
|
||||
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) {
|
||||
tCalculators[t].join();
|
||||
}
|
||||
// Now we add in the data from this completed thread:
|
||||
MatrixUtils.addInPlace(totalErrors, runners[t].getReturnValues());
|
||||
}
|
||||
}
|
||||
|
||||
// Finalise the results:
|
||||
if (debug) {
|
||||
System.out.printf("Total prediction error from variable %d to variable %d=",
|
||||
predictFirstVariable ? 2 : 1, predictFirstVariable ? 1 : 2);
|
||||
MatrixUtils.printArray(System.out, 3, totalErrors);
|
||||
}
|
||||
return totalErrors;
|
||||
}
|
||||
|
||||
/**
|
||||
* Private class to handle multi-threading of the Kraskov algorithms.
|
||||
* Each instance calls partialComputeFromObservations()
|
||||
|
|
@ -573,8 +704,8 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
|
|||
* or throw any exception that was encountered by the
|
||||
* thread.
|
||||
*
|
||||
* @return an exception previously encountered by this thread.
|
||||
* @throws Exception
|
||||
* @return the relevant return values from this part of the data
|
||||
* @throws Exception an exception previously encountered by this thread.
|
||||
*/
|
||||
public double[] getReturnValues() throws Exception {
|
||||
if (problem != null) {
|
||||
|
|
@ -598,4 +729,149 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
|
|||
}
|
||||
}
|
||||
// end class MiKraskovThreadRunner
|
||||
|
||||
/**
|
||||
* Private class to handle multi-threading of the prediction from
|
||||
* k nearest neighbours.
|
||||
* Each instance calls partialComputePredictionErrorFromObservations()
|
||||
* to compute nearest neighbours for a part of the data.
|
||||
*
|
||||
*
|
||||
* @author Joseph Lizier (<a href="joseph.lizier at gmail.com">email</a>,
|
||||
* <a href="http://lizier.me/joseph/">www</a>)
|
||||
* @author Ipek Özdemir
|
||||
*/
|
||||
private class MiKraskovPredictionThreadRunner implements Runnable {
|
||||
protected MutualInfoCalculatorMultiVariateKraskov miCalc;
|
||||
protected int myStartTimePoint;
|
||||
protected int numberOfTimePoints;
|
||||
protected int kNNs;
|
||||
protected boolean predictFirstVariable;
|
||||
|
||||
protected double[] returnValues = null;
|
||||
protected Exception problem = null;
|
||||
|
||||
public MiKraskovPredictionThreadRunner(
|
||||
MutualInfoCalculatorMultiVariateKraskov miCalc,
|
||||
int myStartTimePoint, int numberOfTimePoints,
|
||||
int kNNs, boolean predictFirstVariable) {
|
||||
this.miCalc = miCalc;
|
||||
this.myStartTimePoint = myStartTimePoint;
|
||||
this.numberOfTimePoints = numberOfTimePoints;
|
||||
this.kNNs = kNNs;
|
||||
this.predictFirstVariable = predictFirstVariable;
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the sum of prediction errors from this part of the data,
|
||||
* or throw any exception that was encountered by the
|
||||
* thread.
|
||||
*
|
||||
* @return sum of prediction errors from this part of the data, for
|
||||
* each of the dimensions of the relevant variable
|
||||
* @throws Exception an exception previously encountered by this thread.
|
||||
*/
|
||||
public double[] getReturnValues() throws Exception {
|
||||
if (problem != null) {
|
||||
throw problem;
|
||||
}
|
||||
return returnValues;
|
||||
}
|
||||
|
||||
/**
|
||||
* Start the thread for the given parameters
|
||||
*/
|
||||
public void run() {
|
||||
try {
|
||||
returnValues = miCalc.partialComputePredictionErrorFromObservations(
|
||||
myStartTimePoint, numberOfTimePoints, kNNs, predictFirstVariable);
|
||||
} catch (Exception e) {
|
||||
// Store the exception for later retrieval
|
||||
problem = e;
|
||||
return;
|
||||
}
|
||||
}
|
||||
}
|
||||
// end class MiKraskovPredictionThreadRunner
|
||||
|
||||
/**
|
||||
* Protected method to be used internally for threaded implementations.
|
||||
* This method implements the guts of examining prediction errors in one variable
|
||||
* from the k nearest neighbours of the other.
|
||||
* It is intended to be called by one thread to work on that specific
|
||||
* sub-set of the data.
|
||||
*
|
||||
* @param startTimePoint start time for the partial set we examine
|
||||
* @param numTimePoints number of time points (including startTimePoint to examine)
|
||||
* @param kNNs number of nearest neighbours to use
|
||||
* @param predictFirstVariable whether to use the second variable to predict the first (true)
|
||||
* or first variable to predict the second (false)
|
||||
* @return an array of the sum of square prediction errors for each dimension within the predicted
|
||||
* variable.
|
||||
* @throws Exception
|
||||
*/
|
||||
protected double[] partialComputePredictionErrorFromObservations(
|
||||
int startTimePoint, int numTimePoints, int kNNs, boolean predictFirstVariable) throws Exception {
|
||||
|
||||
double startTime = Calendar.getInstance().getTimeInMillis();
|
||||
|
||||
double[] totalErrors = new double[predictFirstVariable ? dimensionsSource : dimensionsDest];
|
||||
|
||||
for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) {
|
||||
// Find the k nearest neighbours for the relevant predictor variable
|
||||
if (predictFirstVariable) {
|
||||
// First variable value to predict:
|
||||
double[] sourceValueToPredict = sourceObservations[t];
|
||||
// Find kNNs of second variable
|
||||
PriorityQueue<NeighbourNodeData> nnPQ =
|
||||
nnSearcherDest.findKNearestNeighbours(kNNs, t, dynCorrExclTime);
|
||||
double[] predictedValue = new double[dimensionsSource];
|
||||
for (NeighbourNodeData kthNnData : nnPQ) {
|
||||
// Retrieve the source value this corresponds to
|
||||
double[] neighbourSourceValue = sourceObservations[kthNnData.sampleIndex];
|
||||
// And include it's contribution in the prediction
|
||||
for (int d = 0; d < dimensionsSource; d++) {
|
||||
predictedValue[d] += neighbourSourceValue[d];
|
||||
}
|
||||
}
|
||||
// Now add in the square prediction errors from the prediction:
|
||||
for (int d = 0; d < dimensionsSource; d++) {
|
||||
predictedValue[d] /= (double) kNNs;
|
||||
totalErrors[d] += (sourceValueToPredict[d] - predictedValue[d]) *
|
||||
(sourceValueToPredict[d] - predictedValue[d]);
|
||||
}
|
||||
} else { // predict second variable
|
||||
// Second variable value to predict:
|
||||
double[] destValueToPredict = destObservations[t];
|
||||
// Find kNNs of first variable
|
||||
PriorityQueue<NeighbourNodeData> nnPQ =
|
||||
nnSearcherSource.findKNearestNeighbours(kNNs, t, dynCorrExclTime);
|
||||
double[] predictedValue = new double[dimensionsDest];
|
||||
for (NeighbourNodeData kthNnData : nnPQ) {
|
||||
// Retrieve the dest value this corresponds to
|
||||
double[] neighbourDestValue = destObservations[kthNnData.sampleIndex];
|
||||
// And include it's contribution in the prediction
|
||||
for (int d = 0; d < dimensionsDest; d++) {
|
||||
predictedValue[d] += neighbourDestValue[d];
|
||||
}
|
||||
}
|
||||
// Now add in the square prediction errors from the prediction:
|
||||
for (int d = 0; d < dimensionsDest; d++) {
|
||||
predictedValue[d] /= (double) kNNs;
|
||||
totalErrors[d] += (destValueToPredict[d] - predictedValue[d]) *
|
||||
(destValueToPredict[d] - predictedValue[d]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (debug) {
|
||||
Calendar rightNow2 = Calendar.getInstance();
|
||||
long endTime = rightNow2.getTimeInMillis();
|
||||
System.out.println("Subset " + startTimePoint + ":" +
|
||||
(startTimePoint + numTimePoints) + " Calculation time: " +
|
||||
((endTime - startTime)/1000.0) + " sec" );
|
||||
}
|
||||
|
||||
return totalErrors;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -3709,6 +3709,13 @@ public class MatrixUtils {
|
|||
out.println();
|
||||
}
|
||||
|
||||
public static void printArray(PrintStream out, int decimalPlaces, double[] array) {
|
||||
for (int r = 0; r < array.length; r++) {
|
||||
out.printf(String.format("%%.%df ", decimalPlaces), array[r]);
|
||||
}
|
||||
out.println();
|
||||
}
|
||||
|
||||
public static void printArray(PrintStream out, int[] array) {
|
||||
for (int r = 0; r < array.length; r++) {
|
||||
out.print(array[r] + " ");
|
||||
|
|
|
|||
|
|
@ -120,6 +120,22 @@ public abstract class NearestNeighbourSearcher {
|
|||
normTypeToUse = validateNormType(normType);
|
||||
}
|
||||
|
||||
/**
|
||||
*
|
||||
* @return the norm type in use
|
||||
*/
|
||||
public int getNormType() {
|
||||
return normTypeToUse;
|
||||
}
|
||||
|
||||
/**
|
||||
*
|
||||
* @return the norm type in use as a String
|
||||
*/
|
||||
public String getNormTypeAsString() {
|
||||
return convertNormTypeToString(normTypeToUse);
|
||||
}
|
||||
|
||||
/**
|
||||
* Validate whether a specified norm type is supported,
|
||||
* and return the int corresponding to that type,
|
||||
|
|
@ -150,6 +166,24 @@ public abstract class NearestNeighbourSearcher {
|
|||
" is not supported in NearestNeighbourSearcher");
|
||||
}
|
||||
|
||||
/**
|
||||
* @param an identifier for a norm type
|
||||
* @return representation of that norm type as a String
|
||||
*/
|
||||
public static String convertNormTypeToString(int normType) {
|
||||
if (normType == EuclideanUtils.NORM_EUCLIDEAN_SQUARED) {
|
||||
return EuclideanUtils.NORM_EUCLIDEAN_SQUARED_STRING;
|
||||
} else if (normType == EuclideanUtils.NORM_EUCLIDEAN) {
|
||||
return EuclideanUtils.NORM_EUCLIDEAN_STRING;
|
||||
} else if (normType == EuclideanUtils.NORM_MAX_NORM) {
|
||||
return EuclideanUtils.NORM_MAX_NORM_STRING;
|
||||
}
|
||||
// Execution should never reach this point as we control
|
||||
// what normTypeToUse gets set to; it might be possible though
|
||||
// if a child class mis-handles the value
|
||||
throw new Error("normTypeToUse set to an invalid value: " + normType);
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the node which is the nearest neighbour for a given
|
||||
* sample index in the data set. The node itself is
|
||||
|
|
|
|||
Loading…
Reference in New Issue