diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..45fe1c0 Binary files /dev/null and b/.gitignore differ diff --git a/build.xml b/build.xml index 49e3283..18b7bbe 100755 --- a/build.xml +++ b/build.xml @@ -1,5 +1,5 @@ - + Build file for the Java Information Dynamics Toolkit @@ -19,18 +19,25 @@ + + + - + + + + + - + - + + + + + + + @@ -55,17 +68,19 @@ - + + - + + @@ -138,6 +153,21 @@ *********************************** --> + + + + + + + + + + + + + + + diff --git a/demos/java/example10GPUBenchmark.sh b/demos/java/example10GPUBenchmark.sh new file mode 100644 index 0000000..0a2b143 --- /dev/null +++ b/demos/java/example10GPUBenchmark.sh @@ -0,0 +1,18 @@ +#!/bin/bash + +# Make sure the latest example source file is compiled: +# Note: this does not nulify our claim that you can compile the code +# once and then only change the properties file - we only compile here +# every time this script is run so that we can capture any changes the user +# made to the demo source code. The demo as written should be compiled once +# then one can make dynamic changes to the props file and simply +# run the class file without recompiling. +javac -classpath "../../infodynamics.jar" "infodynamics/demos/Example10GPUBenchmark.java" + +# Run the example, feeding in the file names as the command line argument +file1=WhiteNoise.txt +file2=Correlated2D.txt +java -classpath ".:../../infodynamics.jar" infodynamics.demos.Example10GPUBenchmark $file1 $file2 + +# Plot the results +python plotExample10BenchmarkResults.py $file1 $file2 diff --git a/demos/java/infodynamics/demos/Example10GPUBenchmark.java b/demos/java/infodynamics/demos/Example10GPUBenchmark.java new file mode 100755 index 0000000..d3d5269 --- /dev/null +++ b/demos/java/infodynamics/demos/Example10GPUBenchmark.java @@ -0,0 +1,216 @@ +/* + * 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.demos; + +import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate; +import infodynamics.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov1; +import infodynamics.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov2; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; +import infodynamics.utils.ParsedProperties; + +import java.io.FileWriter; +import java.io.PrintWriter; +import java.io.IOException; + +import java.lang.Math; + +/** + * = Example 10 - GPU benchmark script = + * + * This class is used to demonstrate how the GPU module is activated and + * includes a benchmark to test how fast it is, in comparison with its CPU + * counterpart. + * + * If run from Linux, the script plotExample10BenchmarkResults.py + * will use the Python library matplotlib to plot the results of the benchmark. + * + * @author Pedro AM Mediano + * + */ +public class Example10GPUBenchmark { + + + /** + * Calculate MI for the given multivariate time series using the Kraskov + * algorithm, either with CPU or GPU implementation. + * + * @param src Source variable for MI calculation + * @param tgt Target variable for MI calculation + * @param useGPU Whether or not to use the GPU for this calculation + * @return An array with two doubles. The first component is the time in ms + * it took to run the calculation and the second component is the actual + * value of the calculated MI. + */ + public static double[] runEstimator(double[][] src, double[][] tgt, boolean useGPU) throws Exception { + MutualInfoCalculatorMultiVariate miCalc = new MutualInfoCalculatorMultiVariateKraskov1(); + if (useGPU) + miCalc.setProperty("USE_GPU", "true"); + miCalc.setProperty("NOISE_LEVEL_TO_ADD", "0"); + miCalc.setProperty("NORMALISE", "false"); + miCalc.setProperty("k", "4"); + miCalc.initialise(src[0].length, tgt[0].length); + miCalc.setObservations(src, tgt); + double[] timeAndValue = new double[2]; + long startTime = System.nanoTime(); + timeAndValue[1] = miCalc.computeAverageLocalOfObservations(); + timeAndValue[0] = (System.nanoTime() - startTime)/1000000.0; + return timeAndValue; + } + + /** + * Benchmark on a small predefined dataset, for debugging purposes. + */ + public static void Benchmark0() throws Exception { + double[][] src = new double[][] {{0}, {1}, {2}, {3}, {4}, {5}, {6}}; + double[][] tgt = new double[][] {{0}, {0}, {0}, {0}, {0}, {0}, {0}}; + System.out.printf("CPU value: %f\n", runEstimator(src, tgt, false)[1]); + System.out.printf("GPU value: %f\n", runEstimator(src, tgt, true)[1]); + + } + + /** + * Benchmark on white Gaussian data of arbitrary dimension. + */ + public static void Benchmark1(String filename) throws Exception { + + boolean append_to_file = false; + FileWriter write = new FileWriter(filename, append_to_file); + PrintWriter print_line = new PrintWriter(write); + + int[] N_vec = new int[] {500, 1000, 2000, 4000}; + int[] D_vec = new int[] {1, 3, 5}; + + RandomGenerator rg = new RandomGenerator(); + int nb_repetitions = 5; + + for (int i = 0; i < N_vec.length; i++) { + for (int j = 0; j < D_vec.length; j++) { + + int N = N_vec[i]; + int D = D_vec[j]; + + double[][] src = rg.generateNormalData(N, D, 0, 1); + double[][] tgt = rg.generateNormalData(N, D, 0, 1); + + double cpuTime = 0, gpuTime = 0; + for (int r = 0; r < nb_repetitions; r++) { + // CPU calculation + double[] cpuVals = runEstimator(src, tgt, false); + cpuTime += cpuVals[0]; + + // GPU calculation + double[] gpuVals = runEstimator(src, tgt, true); + gpuTime += gpuVals[0]; + + if (Math.abs(cpuVals[1] - gpuVals[1]) > 1e-4) { + System.out.printf("Values differ. CPU: %f, GPU: %f\n", cpuVals[1], gpuVals[1]); + } + + } + + print_line.printf("%d\t%d\t%f\t%f\n", N, D, + cpuTime/nb_repetitions, gpuTime/nb_repetitions); + + } + } + + print_line.close(); + + return; + } + + + + /** + * Benchmark on 2D correlated Gaussian. + * + * Reference: + * + * http://math.stackexchange.com/questions/446093/generate-correlated-normal-random-variables + * + */ + public static void Benchmark2(String filename) throws Exception { + + boolean append_to_file = false; + FileWriter write = new FileWriter(filename, append_to_file); + PrintWriter print_line = new PrintWriter(write); + + int[] N_vec = new int[] {500, 1000, 2000, 4000}; + double[] corr_vec = new double[] {0, 0.2, 0.4, 0.6}; + + RandomGenerator rg = new RandomGenerator(); + int nb_repetitions = 5; + + for (int i = 0; i < N_vec.length; i++) { + for (int j = 0; j < corr_vec.length; j++) { + + int N = N_vec[i]; + double corr = corr_vec[j]; + + double[] r1 = rg.generateNormalData(N, 0, 1); + double[] r2 = rg.generateNormalData(N, 0, 1); + double[][] src = new double[N][1]; + double[][] tgt = new double[N][1]; + for (int t = 0; t < N; t++) { + src[t][0] = r1[t]; + tgt[t][0] = corr*src[t][0] + Math.sqrt(1 - corr*corr)*r2[t]; + } + + double cpuTime = 0, gpuTime = 0; + for (int r = 0; r < nb_repetitions; r++) { + // CPU calculation + double[] cpuVals = runEstimator(src, tgt, false); + cpuTime += cpuVals[0]; + + // GPU calculation + double[] gpuVals = runEstimator(src, tgt, true); + gpuTime += gpuVals[0]; + + if (Math.abs(cpuVals[1] - gpuVals[1]) > 1e-4) { + System.out.printf("Values differ. CPU: %f, GPU: %f\n", cpuVals[1], gpuVals[1]); + } + } + + // Save time in milliseconds + print_line.printf("%d\t%f\t%f\t%f\n", N, corr, + cpuTime/nb_repetitions, gpuTime/nb_repetitions); + + } + } + + print_line.close(); + + return; + } + + /** + * Run several of the benchmarks above and write the results to a + * .txt file. + * + * @param args List of file names to save the results of the benchmark. + */ + public static void main(String[] args) throws Exception { + // Benchmark0(); + Benchmark1(args[0]); + Benchmark2(args[1]); + return; + } + +} diff --git a/demos/java/plotExample10BenchmarkResults.py b/demos/java/plotExample10BenchmarkResults.py new file mode 100644 index 0000000..f127ec3 --- /dev/null +++ b/demos/java/plotExample10BenchmarkResults.py @@ -0,0 +1,47 @@ +#!/usr/bin/python + +""" +plotExample10BenchmarkResults.py: Show 3D plots with the result of benchmarking +the CPU and GPU implementations of the KSG algorithm for mutual information. +""" + +__author__ = "Pedro AM Mediano" +__email__ = "pmediano@imperial.ac.uk" +__license__ = "GPL" + +import sys +import numpy as np +import matplotlib.pyplot as pl +from mpl_toolkits.mplot3d import Axes3D + +## Load files +filename1 = sys.argv[1] +filename2 = sys.argv[2] +res1 = np.loadtxt(filename1) +res2 = np.loadtxt(filename2) + +## There is a bug with the legends of 3D scatterplots, so we have to save +# the plot style for the legend separately +img1 = pl.Line2D([0],[0], linestyle="none", c='r', marker='o') +img2 = pl.Line2D([0],[0], linestyle="none", c='b', marker='^') + +## Plot times for first benchmark +fig1 = pl.figure() +ax1 = fig1.add_subplot(111, projection='3d') +ax1.scatter(np.log10(res1[:,0]), res1[:,1], res1[:,2], c='r', marker='o') +ax1.scatter(np.log10(res1[:,0]), res1[:,1], res1[:,3], c='b', marker='^') +ax1.legend([img1, img2], ['CPU', 'GPU'], numpoints = 1) +ax1.set_xlabel(r'$\log_{10} N$'); ax1.set_ylabel('D'); ax1.set_zlabel('Time [ms]'); +ax1.set_title(filename1.partition(".")[0]) + +## Plot times for second benchmark +fig2 = pl.figure() +ax2 = fig2.add_subplot(111, projection='3d') +ax2.scatter(np.log10(res2[:,0]), res2[:,1], res2[:,2], c='r', marker='o') +ax2.scatter(np.log10(res2[:,0]), res2[:,1], res2[:,3], c='b', marker='^') +ax2.legend([img1, img2], ['CPU', 'GPU'], numpoints = 1) +ax2.set_xlabel(r'$\log_{10} N$'); ax2.set_ylabel(r'$\rho$'); ax2.set_zlabel('Time [ms]'); +ax2.set_title(filename2.partition(".")[0]) + +pl.show() + diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java index e2a269f..85f42cb 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java @@ -21,6 +21,8 @@ package infodynamics.measures.continuous.kraskov; import java.util.Calendar; import java.util.PriorityQueue; import java.util.Random; +import java.nio.file.Path; +import java.nio.file.Paths; import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate; import infodynamics.measures.continuous.MutualInfoMultiVariateCommon; @@ -46,7 +48,7 @@ import infodynamics.utils.NeighbourNodeData; *

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

    - *
  • For constructors see the child classes.
    • + *
  • For constructors see the child classes.
    • *
  • Further properties are defined in {@link #setProperty(String, String)}.
    • *
  • Computed values are in nats, not bits!
    • *
@@ -54,7 +56,7 @@ import infodynamics.utils.NeighbourNodeData; * *

References:
*

    - *
  • Kraskov, A., Stoegbauer, H., Grassberger, P., + *
  • Kraskov, A., Stoegbauer, H., Grassberger, P., * "Estimating mutual information", * Physical Review E 69, (2004) 066138.
    • *
@@ -64,822 +66,907 @@ import infodynamics.utils.NeighbourNodeData; * @author Ipek Özdemir */ public abstract class MutualInfoCalculatorMultiVariateKraskov - extends MutualInfoMultiVariateCommon - implements MutualInfoCalculatorMultiVariate { + extends MutualInfoMultiVariateCommon + implements MutualInfoCalculatorMultiVariate { - /** - * 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 in the full joint space (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 whether to normalise the incoming data to - * mean 0, standard deviation 1 (default true) - */ - public static final String PROP_NORMALISE = "NORMALISE"; - /** - * Property name for an amount of random Gaussian noise to be - * added to the data (default is 1e-8, matching the MILCA toolkit). - */ - public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD"; - /** - * Property name for 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 normalise the incoming data - */ - protected boolean normalise = true; - /** - * 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(); - /** - * Private variable to record which KSG algorithm number - * this instance is implementing - */ - protected boolean isAlgorithm1 = false; - - /** - * protected k-d tree data structure (for fast nearest neighbour searches) - * representing the joint source-dest space - */ - protected KdTree kdTreeJoint; - /** - * protected data structure (for fast nearest neighbour searches) - * representing the source space - */ - protected NearestNeighbourSearcher nnSearcherSource; - /** - * protected data structure (for fast nearest neighbour searches) - * representing the dest space - */ - protected NearestNeighbourSearcher nnSearcherDest; - - /** - * 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; - - /** - * Construct an instance of the KSG MI calculator - */ - public MutualInfoCalculatorMultiVariateKraskov() { - super(); - } + /** + * 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 in the full joint space (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 whether to normalise the incoming data to + * mean 0, standard deviation 1 (default true) + */ + public static final String PROP_NORMALISE = "NORMALISE"; + /** + * Property name for an amount of random Gaussian noise to be + * added to the data (default is 1e-8, matching the MILCA toolkit). + */ + public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD"; + /** + * Property name for 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"; + /** + * Property name for the flag to enable or disable the GPU module. + */ + public static final String PROP_USE_GPU = "USE_GPU"; + + /** + * Whether to normalise the incoming data + */ + protected boolean normalise = true; + /** + * 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(); + /** + * Private variable to record which KSG algorithm number + * this instance is implementing + */ + protected boolean isAlgorithm1 = false; + /** + * Whether to enable the GPU module + */ + protected boolean useGPU = false; + /** + * Check whether C native code has been loaded + */ + protected boolean cudaLibraryLoaded = false; + + /** + * protected k-d tree data structure (for fast nearest neighbour searches) + * representing the joint source-dest space + */ + protected KdTree kdTreeJoint; + /** + * protected data structure (for fast nearest neighbour searches) + * representing the source space + */ + protected NearestNeighbourSearcher nnSearcherSource; + /** + * protected data structure (for fast nearest neighbour searches) + * representing the dest space + */ + protected NearestNeighbourSearcher nnSearcherDest; + + /** + * 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; + + /** + * Construct an instance of the KSG MI calculator + */ + public MutualInfoCalculatorMultiVariateKraskov() { + super(); + } - @Override - public void initialise(int sourceDimensions, int destDimensions) { - kdTreeJoint = null; - nnSearcherSource = null; - nnSearcherDest = null; - super.initialise(sourceDimensions, destDimensions); - } + @Override + public void initialise(int sourceDimensions, int destDimensions) { + kdTreeJoint = null; + nnSearcherSource = null; + nnSearcherDest = null; + super.initialise(sourceDimensions, destDimensions); + } - /** - * Sets properties for the KSG MI 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_NORMALISE} -- whether to normalise the incoming individual - * variables to mean 0 and standard deviation 1 (true by default)
    • - *
  • {@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.
    • - *
  • {@link #PROP_NUM_THREADS} -- the integer number of parallel threads - * to use in the computation. Can be passed as a string "USE_ALL" - * to use all available processors on the machine. - * Default is "USE_ALL". - *
  • any valid properties for {@link MutualInfoMultiVariateCommon#setProperty(String, String)}.
    • - *
- * - *

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_NORMALISE)) { - normalise = Boolean.parseBoolean(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); - } - } + /** + * Sets properties for the KSG MI 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_NORMALISE} -- whether to normalise the incoming individual + * variables to mean 0 and standard deviation 1 (true by default)
    • + *
  • {@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.
    • + *
  • {@link #PROP_NUM_THREADS} -- the integer number of parallel threads + * to use in the computation. Can be passed as a string "USE_ALL" + * to use all available processors on the machine. + * Default is "USE_ALL". + *
  • any valid properties for {@link MutualInfoMultiVariateCommon#setProperty(String, String)}.
    • + *
+ * + *

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_NORMALISE)) { + normalise = Boolean.parseBoolean(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 if (propertyName.equalsIgnoreCase(PROP_USE_GPU)) { + useGPU = Boolean.parseBoolean(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); + } + } - /** - * Get property values for the calculator. - * - *

Valid property names, and what their - * values should represent, are the same as those for - * {@link #setProperty(String, String)}

- * - *

Unknown property values are responded to with a null return value.

- * - * @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); - } - } + /** + * Get property values for the calculator. + * + *

Valid property names, and what their + * values should represent, are the same as those for + * {@link #setProperty(String, String)}

+ * + *

Unknown property values are responded to with a null return value.

+ * + * @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() - */ - @Override - public void finaliseAddObservations() throws Exception { - // Allow the parent to generate the data for us first - super.finaliseAddObservations(); - - if (dynCorrExcl && addedMoreThanOneObservationSet) { - // We have not properly implemented dynamic correlation exclusion for - // multiple observation sets, so throw an error - throw new RuntimeException("Addition of multiple observation sets is not currently " + - "supported with property " + PROP_DYN_CORR_EXCL_TIME + " set"); - } - - if (totalObservations <= k + 2*dynCorrExclTime) { - throw new Exception("There are less observations provided (" + - totalObservations + - ") 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(sourceObservations); - MatrixUtils.normalise(destObservations); - } - - if (addNoise) { - Random random = new Random(); - // Add Gaussian noise of std dev noiseLevel to the data - for (int r = 0; r < sourceObservations.length; r++) { - for (int c = 0; c < dimensionsSource; c++) { - sourceObservations[r][c] += - random.nextGaussian()*noiseLevel; - } - for (int c = 0; c < dimensionsDest; c++) { - destObservations[r][c] += - random.nextGaussian()*noiseLevel; - } - } - } + /* (non-Javadoc) + * @see infodynamics.measures.continuous.MutualInfoMultiVariateCommon#finaliseAddObservations() + */ + @Override + public void finaliseAddObservations() throws Exception { + // Allow the parent to generate the data for us first + super.finaliseAddObservations(); + + if (dynCorrExcl && addedMoreThanOneObservationSet) { + // We have not properly implemented dynamic correlation exclusion for + // multiple observation sets, so throw an error + throw new RuntimeException("Addition of multiple observation sets is not currently " + + "supported with property " + PROP_DYN_CORR_EXCL_TIME + " set"); + } + + if (totalObservations <= k + 2*dynCorrExclTime) { + throw new Exception("There are less observations provided (" + + totalObservations + + ") 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(sourceObservations); + MatrixUtils.normalise(destObservations); + } + + if (addNoise) { + Random random = new Random(); + // Add Gaussian noise of std dev noiseLevel to the data + for (int r = 0; r < sourceObservations.length; r++) { + for (int c = 0; c < dimensionsSource; c++) { + sourceObservations[r][c] += + random.nextGaussian()*noiseLevel; + } + for (int c = 0; c < dimensionsDest; c++) { + destObservations[r][c] += + random.nextGaussian()*noiseLevel; + } + } + } - // Set the constants: - digammaK = MathsUtils.digamma(k); - digammaN = MathsUtils.digamma(totalObservations); - } + // Set the constants: + digammaK = MathsUtils.digamma(k); + digammaN = MathsUtils.digamma(totalObservations); + } - /** - * {@inheritDoc} - * - * @return the average MI in nats (not bits!) - */ - public double computeAverageLocalOfObservations() throws Exception { - // Compute the MI - double startTime = Calendar.getInstance().getTimeInMillis(); - lastAverage = computeFromObservations(false)[0]; - miComputed = 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 average MI in nats (not bits!) + */ + public double computeAverageLocalOfObservations() throws Exception { + // Compute the MI + double startTime = Calendar.getInstance().getTimeInMillis(); + lastAverage = computeFromObservations(false)[0]; + miComputed = 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 MI under the new ordering, in nats (not bits!). - */ - public double computeAverageLocalOfObservations(int[] reordering) throws Exception { - if (reordering == null) { - return computeAverageLocalOfObservations(); - } - - // Save internal variables pertaining to the original order for - // later reinstatement: - // (don't need to save and reinstate kdTreeSource, since we're not - // altering the source data order). - KdTree originalKdTreeJoint = kdTreeJoint; - kdTreeJoint = null; // So that it is rebuilt for the new ordering - NearestNeighbourSearcher originalKdTreeDest = nnSearcherDest; - nnSearcherDest = null; // So that it is rebuilt for the new ordering - double[][] originalData2 = destObservations; - - // Generate a new re-ordered data2 - destObservations = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData2, reordering); - // Compute the MI - double newMI = computeFromObservations(false)[0]; - - // restore original variables: - destObservations = originalData2; - kdTreeJoint = originalKdTreeJoint; - nnSearcherDest = originalKdTreeDest; - - return newMI; - } + /** + * {@inheritDoc} + * + * @return the MI under the new ordering, in nats (not bits!). + */ + public double computeAverageLocalOfObservations(int[] reordering) throws Exception { + if (reordering == null) { + return computeAverageLocalOfObservations(); + } + + // Save internal variables pertaining to the original order for + // later reinstatement: + // (don't need to save and reinstate kdTreeSource, since we're not + // altering the source data order). + KdTree originalKdTreeJoint = kdTreeJoint; + kdTreeJoint = null; // So that it is rebuilt for the new ordering + NearestNeighbourSearcher originalKdTreeDest = nnSearcherDest; + nnSearcherDest = null; // So that it is rebuilt for the new ordering + double[][] originalData2 = destObservations; + + // Generate a new re-ordered data2 + destObservations = MatrixUtils.extractSelectedTimePointsReusingArrays(originalData2, reordering); + // Compute the MI + double newMI = computeFromObservations(false)[0]; + + // restore original variables: + destObservations = originalData2; + kdTreeJoint = originalKdTreeJoint; + nnSearcherDest = originalKdTreeDest; + + return newMI; + } - /** - *

Computes the local values of the MI, - * for each valid observation in the previously supplied observations - * (with PDFs computed using all of the previously supplied observation sets).

- * - *

If the samples were supplied via a single call such as - * {@link #setObservations(double[][], double[][])}, - * then the return value is a single time-series of local - * channel measure values corresponding to these samples.

- * - *

Otherwise where disjoint time-series observations were supplied using several - * calls such as {@link #addObservations(double[][], double[][])} - * then the local values for each disjoint observation set will be appended here - * to create a single "time-series" return array.

- * - * @return the "time-series" of local MIs in nats - * @throws Exception - */ - public double[] computeLocalOfPreviousObservations() throws Exception { - double[] localValues = computeFromObservations(true); - lastAverage = MatrixUtils.mean(localValues); - miComputed = true; - return localValues; - } + /** + *

Computes the local values of the MI, + * for each valid observation in the previously supplied observations + * (with PDFs computed using all of the previously supplied observation sets).

+ * + *

If the samples were supplied via a single call such as + * {@link #setObservations(double[][], double[][])}, + * then the return value is a single time-series of local + * channel measure values corresponding to these samples.

+ * + *

Otherwise where disjoint time-series observations were supplied using several + * calls such as {@link #addObservations(double[][], double[][])} + * then the local values for each disjoint observation set will be appended here + * to create a single "time-series" return array.

+ * + * @return the "time-series" of local MIs in nats + * @throws Exception + */ + public double[] computeLocalOfPreviousObservations() throws Exception { + double[] localValues = computeFromObservations(true); + lastAverage = MatrixUtils.mean(localValues); + miComputed = true; + return localValues; + } - /** - * This method, specified in {@link MutualInfoCalculatorMultiVariate} - * is not implemented yet here. - */ - public double[] computeLocalUsingPreviousObservations(double[][] states1, double[][] states2) throws Exception { - // TODO If implemented, will need to incorporate any time difference here. - // Will also need to handle normalisation of the incoming data - // appropriately - throw new Exception("Local method not implemented yet"); - } - - /** - * This protected method handles the multiple threads which - * computes either the average or local MI (over parts of the total - * observations), computing the x and y - * distances between all tuples in time. - * - *

The method returns:

    - *
  1. for (returnLocals == false), an array of size 1, - * containing the average MI
    2. - *
  2. for local MIs (returnLocals == true), the array of local MI values
    3. - *
- * - * @param returnLocals whether to return an array or local values, or else - * sums of these values - * @return either the average MI, or array of local MI value, in nats not bits - * @throws Exception - */ - protected double[] computeFromObservations(boolean returnLocals) throws Exception { - int N = sourceObservations.length; // number of observations - - double[] returnValues = null; - - ensureKdTreesConstructed(); - - if (numThreads == 1) { - // Single-threaded implementation: - returnValues = partialComputeFromObservations(0, N, returnLocals); - - } else { - // We're going multithreaded: - if (returnLocals) { - // We're computing local MI - returnValues = new double[N]; - } else { - // We're computing average MI - returnValues = new double[MiKraskovThreadRunner.RETURN_ARRAY_LENGTH]; - } - - // 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 Kraskov MI with %d threads (%d timesteps each, plus %d residual)\n", - numThreads, lTimesteps, res); - } - Thread[] tCalculators = new Thread[numThreads]; - MiKraskovThreadRunner[] runners = new MiKraskovThreadRunner[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 MiKraskovThreadRunner(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 MI; 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[MiKraskovThreadRunner.INDEX_SUM_DIGAMMAS] / (double) N; - double avNx = returnValues[MiKraskovThreadRunner.INDEX_SUM_NX] / (double) N; - double avNy = returnValues[MiKraskovThreadRunner.INDEX_SUM_NY] / (double) N; - if (debug) { - System.out.println(String.format("Average n_x=%.3f, Average n_y=%.3f", avNx, avNy)); - } - - // Finalise the average result, depending on which algorithm we are implementing: - if (isAlgorithm1) { - return new double[] { digammaK - averageDiGammas + digammaN }; - } else { - return new double[] { digammaK - (1.0 / (double)k) - averageDiGammas + digammaN }; - } - } - } - - /** - * 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. - * - *

The method returns:

    - *
  1. for average MIs (returnLocals == false), the relevant sums of digamma(n_x+1) and digamma(n_y+1) - * for a partial set of the observations
    2. - *
  2. for local MIs (returnLocals == true), the array of local MI 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) and digamma(n_y+1), then - * sum of n_x and finally sum of n_y (these latter two are for debugging purposes). - * @throws Exception - */ - protected abstract double[] partialComputeFromObservations( - int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception; + /** + * This method, specified in {@link MutualInfoCalculatorMultiVariate} + * is not implemented yet here. + */ + public double[] computeLocalUsingPreviousObservations(double[][] states1, double[][] states2) throws Exception { + // TODO If implemented, will need to incorporate any time difference here. + // Will also need to handle normalisation of the incoming data + // appropriately + throw new Exception("Local method not implemented yet"); + } + + /** + * This protected method handles the multiple threads which + * computes either the average or local MI (over parts of the total + * observations), computing the x and y + * distances between all tuples in time. + * + *

The method returns:

    + *
  1. for (returnLocals == false), an array of size 1, + * containing the average MI
    2. + *
  2. for local MIs (returnLocals == true), the array of local MI values
    3. + *
+ * + * @param returnLocals whether to return an array or local values, or else + * sums of these values + * @return either the average MI, or array of local MI value, in nats not bits + * @throws Exception + */ + protected double[] computeFromObservations(boolean returnLocals) throws Exception { + int N = sourceObservations.length; // number of observations + + double[] returnValues = null; + + ensureKdTreesConstructed(); - /** - * 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() - * to compute nearest neighbours for a part of the data. - * - * - * @author Joseph Lizier (email, - * www) - * @author Ipek Özdemir - */ - private class MiKraskovThreadRunner implements Runnable { - protected MutualInfoCalculatorMultiVariateKraskov 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 static final int INDEX_SUM_NX = 1; - public static final int INDEX_SUM_NY = 2; - public static final int RETURN_ARRAY_LENGTH = 3; + if (useGPU) { + returnValues = gpuComputeFromObservations(0, N, returnLocals); + } else if (numThreads == 1) { + // Single-threaded implementation: + returnValues = partialComputeFromObservations(0, N, returnLocals); + + } else { + // We're going multithreaded: + if (returnLocals) { + // We're computing local MI + returnValues = new double[N]; + } else { + // We're computing average MI + returnValues = new double[MiKraskovThreadRunner.RETURN_ARRAY_LENGTH]; + } + + // 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 Kraskov MI with %d threads (%d timesteps each, plus %d residual)\n", + numThreads, lTimesteps, res); + } + Thread[] tCalculators = new Thread[numThreads]; + MiKraskovThreadRunner[] runners = new MiKraskovThreadRunner[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 MiKraskovThreadRunner(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 MI; 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[MiKraskovThreadRunner.INDEX_SUM_DIGAMMAS] / (double) N; + double avNx = returnValues[MiKraskovThreadRunner.INDEX_SUM_NX] / (double) N; + double avNy = returnValues[MiKraskovThreadRunner.INDEX_SUM_NY] / (double) N; + if (debug) { + System.out.println(String.format("Average n_x=%.3f, Average n_y=%.3f", avNx, avNy)); + } + + // Finalise the average result, depending on which algorithm we are implementing: + if (isAlgorithm1) { + return new double[] { digammaK - averageDiGammas + digammaN }; + } else { + return new double[] { digammaK - (1.0 / (double)k) - averageDiGammas + digammaN }; + } + } + } + + /** + * 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. + * + *

The method returns:

    + *
  1. for average MIs (returnLocals == false), the relevant sums of digamma(n_x+1) and digamma(n_y+1) + * for a partial set of the observations
    2. + *
  2. for local MIs (returnLocals == true), the array of local MI 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) and digamma(n_y+1), then + * sum of n_x and finally sum of n_y (these latter two are for debugging purposes). + * @throws Exception + */ + protected abstract double[] partialComputeFromObservations( + int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception; - public MiKraskovThreadRunner( - MutualInfoCalculatorMultiVariateKraskov 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 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) { - throw problem; - } - 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; - } - } - } - // 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 (email, - * www) - * @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 GPU implementations. + * This method serves the same purpose as partialComputeFromObservations, + * but for GPU computation. Each algorithm must override this method + * and implement a GPU routine to calculate all values in a single call + * to the GPU code. + */ + protected double[] gpuComputeFromObservations( + int startTimePoint, int numTimePoints, boolean returnLocals) throws Exception { - /** - * 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(); + ensureCudaLibraryLoaded(); - 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 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 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; - } + boolean useMaxNorm; + + if ( normType == EuclideanUtils.NORM_MAX_NORM) { + useMaxNorm = true; + } else if ( normType == EuclideanUtils.NORM_EUCLIDEAN || normType == EuclideanUtils.NORM_EUCLIDEAN_SQUARED) { + useMaxNorm = false; + } else { + throw new Exception("Only max and square norms are implemented. Abort."); + } + + double[] res; + + try { + res = MIKraskov(totalObservations, sourceObservations, dimensionsSource, + destObservations, dimensionsDest, k, returnLocals, useMaxNorm, isAlgorithm1); + } catch (Throwable e) { + System.out.println("WARNING. Error in GPU code. Reverting back to CPU."); + e.printStackTrace(); + res = partialComputeFromObservations(0, totalObservations, returnLocals); + } + + return res; + } + + + /** + * Native method to calculate MI in GPU. + */ + private native double[] MIKraskov( + int N, double[][] source, int dimx, double[][] dest, int dimy, + int k, boolean returnLocals, boolean useMaxNorm, boolean isAlgorithm1); + + /** + * 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); + } + } + + /** + * Internal method to ensure that the Cuda native library has been correctly + * loaded. + */ + protected void ensureCudaLibraryLoaded() throws Exception { + + if (!cudaLibraryLoaded) { + + Path jarPath = Paths.get(MutualInfoCalculatorMultiVariateKraskov.class.getProtectionDomain().getCodeSource().getLocation().toURI()); + String fileSep = System.getProperty("file.separator", ""); + if (fileSep.length() == 0) { + throw new Exception("Unable to load GPU library."); + } + String jarFolder = jarPath.toString().substring(0, jarPath.toString().lastIndexOf(fileSep)); + String relPath = fileSep + "java" + fileSep + "source" + fileSep + + "infodynamics" + fileSep + "measures" + fileSep + + "continuous" + fileSep + "kraskov" + fileSep + "cuda" + + fileSep + "libKraskov.so"; + String fullPath = jarFolder + relPath; + System.load(fullPath); + cudaLibraryLoaded = true; + + } + } + + /** + * 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() + * to compute nearest neighbours for a part of the data. + * + * + * @author Joseph Lizier (email, + * www) + * @author Ipek Özdemir + */ + private class MiKraskovThreadRunner implements Runnable { + protected MutualInfoCalculatorMultiVariateKraskov 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 static final int INDEX_SUM_NX = 1; + public static final int INDEX_SUM_NY = 2; + public static final int RETURN_ARRAY_LENGTH = 3; + + public MiKraskovThreadRunner( + MutualInfoCalculatorMultiVariateKraskov 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 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) { + throw problem; + } + 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; + } + } + } + // 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 (email, + * www) + * @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 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 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; + } } diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java index 363d16c..878514a 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java @@ -23,7 +23,9 @@ import java.util.PriorityQueue; import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate; import infodynamics.utils.MathsUtils; +import infodynamics.utils.MatrixUtils; import infodynamics.utils.NeighbourNodeData; +import infodynamics.utils.EuclideanUtils; /** *

Computes the differential mutual information of two given multivariate sets of @@ -80,7 +82,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1 // First element in the PQ is the kth NN, // and epsilon = kthNnData.distance NeighbourNodeData kthNnData = nnPQ.poll(); - + // Count the number of points whose x distance is less // than eps, and whose y distance is less than // epsilon = kthNnData.distance @@ -88,7 +90,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1 t, kthNnData.distance, dynCorrExclTime); int n_y = nnSearcherDest.countPointsStrictlyWithinR( t, kthNnData.distance, dynCorrExclTime); - + sumNx += n_x; sumNy += n_y; // And take the digammas: @@ -103,7 +105,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1 localMi[t-startTimePoint] = digammaK - digammaNxPlusOne - digammaNyPlusOne + digammaN; } } - + if (debug) { Calendar rightNow2 = Calendar.getInstance(); long endTime = rightNow2.getTimeInMillis(); @@ -119,4 +121,5 @@ public class MutualInfoCalculatorMultiVariateKraskov1 return new double[] {sumDiGammas, sumNx, sumNy}; } } + } diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java index f42f726..1d5a12d 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java @@ -129,4 +129,5 @@ public class MutualInfoCalculatorMultiVariateKraskov2 return new double[] {sumDiGammas, sumNx, sumNy}; } } + } diff --git a/java/source/infodynamics/measures/continuous/kraskov/cuda/Makefile b/java/source/infodynamics/measures/continuous/kraskov/cuda/Makefile new file mode 100644 index 0000000..7bf2501 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/cuda/Makefile @@ -0,0 +1,39 @@ + +## Set these if the CUDA Toolkit is not in your PATH +# Location of the CUDA Toolkit binaries and libraries +# CUDA_PATH ?= /usr/local/cuda +# CUDA_INC_PATH ?= $(CUDA_PATH)/include +# CUDA_LIB_PATH ?= $(CUDA_PATH)/lib +# CUDA_BIN_PATH ?= $(CUDA_PATH)/bin + +# Common binaries and flags. TODO: the architecture flag may not be necessary +NVCC ?= nvcc +GCC ?= g++ +CCFLAGS ?= -O3 +# NVCCFLAGS ?= -O3 -arch=sm_30 +NVCCFLAGS ?= -O3 + +ifdef DEBUG + CCFLAGS += -g + NVCCFLAGS += -g -G +endif + +# Common includes and paths for CUDA. This assumes the CUDA toolkit is in PATH +INCLUDES := -I. -I/usr/lib/jvm/java-7-openjdk-amd64/include +LDFLAGS += -L. -lcuda -lcudart + +# Target rules +all: libKraskov.so + +gpuKnnLibrary.o: gpuKnnLibrary.cu helperfunctions.cu gpuKnnBF_kernel.cu + ${NVCC} ${NVCCFLAGS} ${INCLUDES} -Xcompiler -fPIC -c gpuKnnLibrary.cu + +libgpuKnnLibrary.a: gpuKnnLibrary.o + ${AR} -r libgpuKnnLibrary.a gpuKnnLibrary.o + +libKraskov.so: libgpuKnnLibrary.a kraskovCuda.c + ${NVCC} ${NVCCFLAGS} ${INCLUDES} -Xcompiler -fPIC -shared -o libKraskov.so kraskovCuda.c ${LDFLAGS} -lgpuKnnLibrary + +clean: + rm -f *.o *.so *.a + diff --git a/java/source/infodynamics/measures/continuous/kraskov/cuda/digamma.c b/java/source/infodynamics/measures/continuous/kraskov/cuda/digamma.c new file mode 100644 index 0000000..a74e3e1 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/cuda/digamma.c @@ -0,0 +1,120 @@ +/************************************* + An ANSI-C implementation of the digamma-function for real arguments based + on the Chebyshev expansion proposed in appendix E of + http://arXiv.org/abs/math.CA/0403344 . This is identical to the implementation + by Jet Wimp, Math. Comp. vol 15 no 74 (1961) pp 174 (see Table 1). + For other implementations see + the GSL implementation for Psi(Digamma) in + http://www.gnu.org/software/gsl/manual/html_node/Psi-_0028Digamma_0029-Function.html + + Richard J. Mathar, 2005-11-24 + ***************************************/ +#include + +#ifndef M_PIl +/** The constant Pi in high precision */ +#define M_PIl 3.1415926535897932384626433832795029L +#endif +#ifndef M_GAMMAl +/** Euler's constant in high precision */ +#define M_GAMMAl 0.5772156649015328606065120900824024L +#endif +#ifndef M_LN2l +/** the natural logarithm of 2 in high precision */ +#define M_LN2l 0.6931471805599453094172321214581766L +#endif + +/** The digamma function in long double precision. + * @param x the real value of the argument + * @return the value of the digamma (psi) function at that point + * @author Richard J. Mathar + * @since 2005-11-24 + */ +#ifdef __cplusplus +extern "C" { +#endif +long double cpuDigamma(long double x) +{ + int n; + /* force into the interval 1..3 */ + if( x < 0.0L ) + return cpuDigamma(1.0L-x)+M_PIl/tanl(M_PIl*(1.0L-x)) ; /* reflection formula */ + else if( x < 1.0L ) + return cpuDigamma(1.0L+x)-1.0L/x ; + else if ( x == 1.0L) + return -M_GAMMAl ; + else if ( x == 2.0L) + return 1.0L-M_GAMMAl ; + else if ( x == 3.0L) + return 1.5L-M_GAMMAl ; + else if ( x > 3.0L) + /* duplication formula */ + return 0.5L*(cpuDigamma(x/2.0L)+cpuDigamma((x+1.0L)/2.0L))+M_LN2l ; + else { + /* Just for your information, the following lines contain + * the Maple source code to re-generate the table that is + * eventually becoming the Kncoe[] array below + * interface(prettyprint=0) : + * Digits := 63 : + * r := 0 : + * + * for l from 1 to 60 do + * d := binomial(-1/2,l) : + * r := r+d*(-1)^l*(Zeta(2*l+1) -1) ; + * evalf(r) ; + * print(%,evalf(1+Psi(1)-r)) ; + *o d : + * + * for N from 1 to 28 do + * r := 0 : + * n := N-1 : + * + * for l from iquo(n+3,2) to 70 do + * d := 0 : + * for s from 0 to n+1 do + * d := d+(-1)^s*binomial(n+1,s)*binomial((s-1)/2,l) : + * od : + * if 2*l-n > 1 then + * r := r+d*(-1)^l*(Zeta(2*l-n) -1) : + * fi : + * od : + * print(evalf((-1)^n*2*r)) ; + *od : + *quit : + */ + static long double Kncoe[] = { .30459198558715155634315638246624251L, + .72037977439182833573548891941219706L, -.12454959243861367729528855995001087L, + .27769457331927827002810119567456810e-1L, -.67762371439822456447373550186163070e-2L, + .17238755142247705209823876688592170e-2L, -.44817699064252933515310345718960928e-3L, + .11793660000155572716272710617753373e-3L, -.31253894280980134452125172274246963e-4L, + .83173997012173283398932708991137488e-5L, -.22191427643780045431149221890172210e-5L, + .59302266729329346291029599913617915e-6L, -.15863051191470655433559920279603632e-6L, + .42459203983193603241777510648681429e-7L, -.11369129616951114238848106591780146e-7L, + .304502217295931698401459168423403510e-8L, -.81568455080753152802915013641723686e-9L, + .21852324749975455125936715817306383e-9L, -.58546491441689515680751900276454407e-10L, + .15686348450871204869813586459513648e-10L, -.42029496273143231373796179302482033e-11L, + .11261435719264907097227520956710754e-11L, -.30174353636860279765375177200637590e-12L, + .80850955256389526647406571868193768e-13L, -.21663779809421233144009565199997351e-13L, + .58047634271339391495076374966835526e-14L, -.15553767189204733561108869588173845e-14L, + .41676108598040807753707828039353330e-15L, -.11167065064221317094734023242188463e-15L } ; + + register long double Tn_1 = 1.0L ; /* T_{n-1}(x), started at n=1 */ + register long double Tn = x-2.0L ; /* T_{n}(x) , started at n=1 */ + register long double resul = Kncoe[0] + Kncoe[1]*Tn ; + + x -= 2.0L ; + + for(n = 2 ; n < sizeof(Kncoe)/sizeof(long double) ;n++) + { + const long double Tn1 = 2.0L * x * Tn - Tn_1 ; /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */ + resul += Kncoe[n]*Tn1 ; + Tn_1 = Tn ; + Tn = Tn1 ; + } + return resul ; + } +} +#ifdef __cplusplus +} +#endif + diff --git a/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnBF_kernel.cu b/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnBF_kernel.cu new file mode 100644 index 0000000..bb0768f --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnBF_kernel.cu @@ -0,0 +1,455 @@ +#ifndef _TEMPLATE_KERNEL_H_ +#define _TEMPLATE_KERNEL_H_ + +//#include +#ifndef INFINITY +#define INFINITY 0x7F800000 +#endif + +/** + * Typedef function pointer for arbitrary norm functions. + */ +typedef float (*normFunction_t)(const float* g_uquery, const float* g_vpoint, + int pointdim, int signallength); + +/** + * Calculate max norm (L-inf) between two points. + */ +__device__ __host__ float +maxMetricPoints(const float* g_uquery, const float* g_vpoint, int pointdim, int signallength){ + float r_u1; + float r_v1; + float r_d1,r_dim=0; + + r_dim=0; + for(int d=0; d*(kdistances+k)) && (kk;k2--){ + *(kdistances+k2)=*(kdistances+k2-1); + *(kindexes+k2)=*(kindexes+k2-1); + } + //Replace + *(kdistances+k)=distance; + *(kindexes+k)=indexv; + + //printf("\n -> Modificacion pila: %.f %.f. New max distance: %.f", *kdistances, *(kdistances+1), *(kdistances+kth-1)); + return *(kdistances+kth-1); +} + + +/* + * Main KNN kernel. Find the nearest k neighbours to each point according + * to the supplied norm function. + */ +__global__ void +kernelKNNshared(const float* g_uquery, const float* g_vpointset, + int *g_indexes, float* g_distances, const int pointdim, + const int triallength, const int signallength, const int kth, + const int exclude, normFunction_t *normFunction) { + + // Shared memory + extern __shared__ char array[]; + float *kdistances; + int *kindexes; + kdistances = (float*)array; + kindexes = (int*)array+kth*blockDim.x; + + const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x; + const unsigned int itrial = tid / triallength; // indextrial + +if(tidcondition2)){ + float temp_dist = normFunction[0](g_uquery+indexu, g_vpointset+indexv,pointdim, signallength); + if(temp_dist <= r_kdist){ + r_kdist = insertPointKlist(kth,temp_dist,t,kdistances+threadIdx.x*kth,kindexes+threadIdx.x*kth); + } + } + //printf("tid:%d indexes: %d, %d distances: %.f %.f\n",tid, *kindexes, *(kindexes+1), *kdistances, *(kdistances+1)); + } + + __syncthreads(); + // COPY TO GLOBAL MEMORY + for (int k = 0; k < kth; k++) { + g_indexes[tid+k*signallength] = kindexes[threadIdx.x*kth+k]; + g_distances[tid+k*signallength] = kdistances[threadIdx.x*kth+k];//*(kdistances+k); + } +} + +} + + +/* + * Range search for one data point using bruteforce. + */ +__global__ void +kernelBFRSshared(const float* g_uquery, const float* g_vpointset, int *g_npoints, int pointdim, int triallength, int signallength, int exclude, float radius) +{ + + // Shared memory + extern __shared__ char array[]; + int *s_npointsrange; + s_npointsrange = (int*)array; + + const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x; + const unsigned int itrial = tid / triallength; // indextrial + +if(tidcondition2)){ + float temp_dist = maxMetricPoints(g_uquery+indexu, g_vpointset+indexv,pointdim, signallength); + if(temp_dist <= radius){ + s_npointsrange[threadIdx.x]++; + } + } + + } + + __syncthreads(); + + // COPY TO GLOBAL MEMORY + g_npoints[tid] = s_npointsrange[threadIdx.x]; + +} +} + +/* + * Range search using bruteforce in multiple GPUs + */ +__global__ void +kernelBFRSMultishared(const float* g_uquery, const float* g_vpointset, int *g_npoints, int pointdim, int triallength, int signallength, int exclude, const float* vecradius) +{ + + // shared memory + extern __shared__ char array[]; + int *s_npointsrange; + s_npointsrange = (int*)array; + float radius=0; + const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x; + const unsigned int itrial = tid / triallength; // indextrial + +if(tidcondition2)){ + float temp_dist = maxMetricPoints(g_uquery+indexu, g_vpointset+indexv,pointdim, signallength); + if(temp_dist <= radius){ + s_npointsrange[threadIdx.x]++; + } + } + + } + + __syncthreads(); + //printf("\ntid:%d npoints: %d\n",tid, s_npointsrange[threadIdx.x]); + //COPY TO GLOBAL MEMORY + g_npoints[tid] = s_npointsrange[threadIdx.x]; + +} +} + + +/* + * Range search for all data points using bruteforce. + */ +__global__ void +kernelBFRSAllshared(const float* g_uquery, const float* g_vpointset, + int *g_npoints, int pointdim, int triallength, int signallength, + int exclude, const float* vecradius, normFunction_t *normFunction) { + + // Shared memory + extern __shared__ char array[]; + int *s_npointsrange; + s_npointsrange = (int *) array; + float radius = 0; + const unsigned int tid = threadIdx.x + blockDim.x*blockIdx.x; + const unsigned int itrial = tid / triallength; // indextrial + + if(tidcondition2)){ + float temp_dist = normFunction[0](g_uquery+indexu, g_vpointset+indexv,pointdim, signallength); + // PEDRO: For KSG algorithm 1 this should be strictly less than R, and in TRENTOOL code it's less or equal. It's a float comparison, so I don't think it matters anyway. + if(temp_dist < radius){ + s_npointsrange[threadIdx.x]++; + } + } + + } + + __syncthreads(); + + //COPY TO GLOBAL MEMORY + g_npoints[tid] = s_npointsrange[threadIdx.x]; + + } +} + + + +__device__ void digammaXp1(double *y); + +/** + * Optimized reduction kernel, obtained from Mark Harris' lecture on + * GPU optimization: + * https://docs.nvidia.com/cuda/samples/6_Advanced/reduction/doc/reduction.pdf + */ +template +// __device__ void warpReduce(volatile int *sdata, unsigned int tid) { +__device__ void warpReduce(volatile float *sdata, unsigned int tid) { + if (blockSize >= 64) sdata[tid] += sdata[tid + 32]; + if (blockSize >= 32) sdata[tid] += sdata[tid + 16]; + if (blockSize >= 16) sdata[tid] += sdata[tid + 8]; + if (blockSize >= 8) sdata[tid] += sdata[tid + 4]; + if (blockSize >= 4) sdata[tid] += sdata[tid + 2]; + if (blockSize >= 2) sdata[tid] += sdata[tid + 1]; +} + +// template +// __global__ void reduce6(int *g_idata, int *g_odata, unsigned int n) { +__global__ void reduce6(int *g_nx, int *g_ny, float *g_odata, unsigned int n) { +// extern __shared__ int sdata[]; +extern __shared__ float sdata[]; +unsigned int tid = threadIdx.x; +const unsigned int blockSize = 512; +unsigned int i = blockIdx.x*(blockSize*2) + tid; +unsigned int gridSize = blockSize*2*gridDim.x; +sdata[tid] = 0; + +while (i= 512) { if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads(); } +if (blockSize >= 256) { if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads(); } +if (blockSize >= 128) { if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads(); } +if (tid < 32) warpReduce(sdata, tid); +if (tid == 0) g_odata[blockIdx.x] = sdata[0]; +} + + +__global__ void reduce1(int *g_nx, int *g_ny, float *g_odata, unsigned int N) { + extern __shared__ float sdata2[]; + // each thread loads one element from global to shared mem + unsigned int tid = threadIdx.x; + unsigned int i = blockIdx.x*blockDim.x + threadIdx.x; + + if (i < N) { + double dgX = g_nx[i]; + double dgY = g_ny[i]; + digammaXp1(&dgX); + digammaXp1(&dgY); + sdata2[tid] = dgX + dgY; + } else { + sdata2[tid] = 0; + } + __syncthreads(); + + // do reduction in shared mem + for (unsigned int s=1; s < blockDim.x; s *= 2) { + if (tid % (2*s) == 0) { + sdata2[tid] += sdata2[tid + s]; + } + __syncthreads(); + } + // write result for this block to global mem + if (tid == 0) g_odata[blockIdx.x] = sdata2[0]; +} + + + + + + +/** + * Taken from NVIDIA dev forum: + * https://devtalk.nvidia.com/default/topic/516516/kernel-launch-failure-in-matlab/ + */ +__device__ void digammaXp1(double *y) { + + // double x = *y; + double x = *y + 1; + double neginf = -INFINITY; + const double c = 12, + digamma1 = -0.57721566490153286, + trigamma1 = 1.6449340668482264365, /* pi^2/6 */ + s = 1e-6, + s3 = 1./12, + s4 = 1./120, + s5 = 1./252, + s6 = 1./240, + s7 = 1./132; + // s8 = 691./32760, + // s9 = 1./12, + // s10 = 3617./8160; + double result; + + /* Illegal arguments */ + if((x == neginf) || isnan(x)) { + *y = NAN; + return; + } + + /* Singularities */ + if((x <= 0) && (floor(x) == x)) { + *y = neginf; + return; + } + + /* Negative values */ + + /* Use the reflection formula (Jeffrey 11.1.6): + * digamma(-x) = digamma(x+1) + pi*cot(pi*x) + * + * This is related to the identity + * digamma(-x) = digamma(x+1) - digamma(z) + digamma(1-z) + * where z is the fractional part of x + * For example: + + * digamma(-3.1) = 1/3.1 + 1/2.1 + 1/1.1 + 1/0.1 + digamma(1-0.1) + * = digamma(4.1) - digamma(0.1) + digamma(1-0.1) + * Then we use + * digamma(1-z) - digamma(z) = pi*cot(pi*z) + * + + if(x < 0) { + *p = digamma(p,1-x) + M_PI/tan(-M_PI*x); + return; + } + */ + + /* Use Taylor series if argument <= S */ + if(x <= s) { + *y = digamma1 - 1/x + trigamma1*x; + return; + } + + /* Reduce to digamma(X + N) where (X + N) >= C */ + result = 0; + while(x < c) { + result -= 1/x; + x++; + } + + /* Use de Moivre's expansion if argument >= C */ + /* This expansion can be computed in Maple via asympt(Psi(x),x) */ + if(x >= c) { + double r = 1/x, t; + result += log(x) - 0.5*r; + r *= r; +#if 0 + result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7)))); +#else + /* this version for lame compilers */ + t = (s5 - r * (s6 - r * s7)); + result -= r * (s3 - r * (s4 - r * t)); +#endif + } + + /* assign the result to the pointer*/ + *y = result; + return; + +} + + + + + + +#endif // #ifndef _TEMPLATE_KERNEL_H_ + + diff --git a/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnLibrary.cu b/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnLibrary.cu new file mode 100644 index 0000000..ec3acb4 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnLibrary.cu @@ -0,0 +1,278 @@ +#include "gpuKnnBF_kernel.cu" +#include +#include "helperfunctions.cu" + +#ifdef __cplusplus +extern "C" { +#endif +int cudaFindKnn(int* h_bf_indexes, float* h_bf_distances, float* h_pointset, + float* h_query, int kth, int thelier, int nchunks, int pointdim, + int signallength, unsigned int useMaxNorm) { + float *d_bf_pointset, *d_bf_query; + int *d_bf_indexes; + float *d_bf_distances; + + unsigned int meminputsignalquerypointset= pointdim * signallength * sizeof(float); + unsigned int mem_bfcl_outputsignaldistances= kth * signallength * sizeof(float); + unsigned int mem_bfcl_outputsignalindexes = kth * signallength * sizeof(int); + + checkCudaErrors( cudaMalloc( (void**) &(d_bf_query), meminputsignalquerypointset)); + checkCudaErrors( cudaMalloc( (void**) &(d_bf_pointset), meminputsignalquerypointset)); + //GPU output + checkCudaErrors( cudaMalloc( (void**) &(d_bf_distances), mem_bfcl_outputsignaldistances )); + checkCudaErrors( cudaMalloc( (void**) &(d_bf_indexes), mem_bfcl_outputsignalindexes )); + + cudaError_t error = cudaGetLastError(); + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + + //Upload input data + checkCudaErrors( cudaMemcpy(d_bf_query, h_query, meminputsignalquerypointset, cudaMemcpyHostToDevice )); + checkCudaErrors( cudaMemcpy(d_bf_pointset, h_pointset, meminputsignalquerypointset, cudaMemcpyHostToDevice )); + error = cudaGetLastError(); + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + // Kernel parameters + dim3 threads(1,1,1); + dim3 grid(1,1,1); + threads.x = 512; + grid.x = (signallength-1)/threads.x + 1; + int memkernel = kth*sizeof(float)*threads.x+\ + kth*sizeof(int)*threads.x; + int triallength = signallength / nchunks; + + // Pointer to the function used to calculate norms + normFunction_t *normFunction; + cudaMalloc( (void **) &normFunction, sizeof(normFunction_t) ); + if (useMaxNorm) { + cudaMemcpyFromSymbol(normFunction, pMaxNorm, sizeof(normFunction_t)); + } else { + cudaMemcpyFromSymbol(normFunction, pSquareNorm, sizeof(normFunction_t)); + } + + // Launch kernel + kernelKNNshared<<>>(d_bf_query, d_bf_pointset, + d_bf_indexes, d_bf_distances, pointdim, triallength, signallength, kth, + thelier, normFunction); + + checkCudaErrors( cudaDeviceSynchronize() ); + + //Download result + checkCudaErrors( cudaMemcpy( h_bf_distances, d_bf_distances, mem_bfcl_outputsignaldistances, cudaMemcpyDeviceToHost) ); + checkCudaErrors( cudaMemcpy( h_bf_indexes, d_bf_indexes, mem_bfcl_outputsignalindexes, cudaMemcpyDeviceToHost) ); + error = cudaGetLastError(); + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + + //Free resources + checkCudaErrors(cudaFree(d_bf_query)); + checkCudaErrors(cudaFree(d_bf_pointset)); + checkCudaErrors(cudaFree(d_bf_distances)); + checkCudaErrors(cudaFree(d_bf_indexes)); + cudaFree(normFunction); + // cudaDeviceReset(); + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + + return 1; +} +#ifdef __cplusplus +} +#endif + +#ifdef __cplusplus +extern "C" { +#endif +int cudaFindKnnSetGPU(int* h_bf_indexes, float* h_bf_distances, + float* h_pointset, float* h_query, int kth, int thelier, int nchunks, + int pointdim, int signallength, unsigned int useMaxNorm, int deviceid) { + cudaSetDevice(deviceid); + return cudaFindKnn(h_bf_indexes, h_bf_distances, h_pointset, h_query, + kth, thelier, nchunks, pointdim, signallength, useMaxNorm); +} +#ifdef __cplusplus +} +#endif + +/* + * Range search being radius a vector of length number points in queryset/pointset + */ +#ifdef __cplusplus +extern "C" { +#endif +int cudaFindRSAll(int* h_bf_npointsrange, float* h_pointset, float* h_query, + float* h_vecradius, int thelier, int nchunks, int pointdim, + int signallength, unsigned int useMaxNorm) { + + float *d_bf_pointset, *d_bf_query, *d_bf_vecradius; + int *d_bf_npointsrange; + + unsigned int meminputsignalquerypointset= pointdim * signallength * sizeof(float); + unsigned int mem_bfcl_outputsignalnpointsrange= signallength * sizeof(int); + unsigned int mem_bfcl_inputvecradius = signallength * sizeof(float); + + checkCudaErrors( cudaMalloc( (void**) &(d_bf_query), meminputsignalquerypointset)); + checkCudaErrors( cudaMalloc( (void**) &(d_bf_pointset), meminputsignalquerypointset)); + checkCudaErrors( cudaMalloc( (void**) &(d_bf_npointsrange), mem_bfcl_outputsignalnpointsrange )); + checkCudaErrors( cudaMalloc( (void**) &(d_bf_vecradius), mem_bfcl_inputvecradius )); + + cudaError_t error = cudaGetLastError(); + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + //Upload input data + checkCudaErrors( cudaMemcpy(d_bf_query, h_query, meminputsignalquerypointset, cudaMemcpyHostToDevice )); + checkCudaErrors( cudaMemcpy(d_bf_pointset, h_pointset, meminputsignalquerypointset, cudaMemcpyHostToDevice )); + checkCudaErrors( cudaMemcpy(d_bf_vecradius, h_vecradius, mem_bfcl_inputvecradius, cudaMemcpyHostToDevice )); + + error = cudaGetLastError(); + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + + + // Kernel parameters + dim3 threads(1,1,1); + dim3 grid(1,1,1); + threads.x = 512; + grid.x = (signallength-1)/threads.x + 1; + int memkernel = sizeof(int)*threads.x; + int triallength = signallength / nchunks; + + // Pointer to the function used to calculate norms + normFunction_t *normFunction; + cudaMalloc( (void **) &normFunction, sizeof(normFunction_t) ); + if (useMaxNorm) { + cudaMemcpyFromSymbol(normFunction, pMaxNorm, sizeof(normFunction_t)); + } else { + cudaMemcpyFromSymbol(normFunction, pSquareNorm, sizeof(normFunction_t)); + } + + // Launch kernel + kernelBFRSAllshared<<< grid.x, threads.x, memkernel>>>( + d_bf_query, d_bf_pointset, d_bf_npointsrange, pointdim, + triallength, signallength, thelier, d_bf_vecradius, normFunction); + + checkCudaErrors(cudaDeviceSynchronize()); + + checkCudaErrors( cudaMemcpy( h_bf_npointsrange, d_bf_npointsrange,mem_bfcl_outputsignalnpointsrange, cudaMemcpyDeviceToHost) ); + + + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + + // Free resources + checkCudaErrors(cudaFree(d_bf_query)); + checkCudaErrors(cudaFree(d_bf_pointset)); + checkCudaErrors(cudaFree(d_bf_npointsrange)); + checkCudaErrors(cudaFree(d_bf_vecradius)); + checkCudaErrors(cudaFree(normFunction)); + cudaDeviceReset(); + if(error!=cudaSuccess){ + fprintf(stderr,"%s",cudaGetErrorString(error)); + return 0; + } + + return 1; +} +#ifdef __cplusplus +} +#endif + +#ifdef __cplusplus +extern "C" { +#endif +int cudaFindRSAllSetGPU(int* h_bf_npointsrange, float* h_pointset, + float* h_query, float* h_vecradius, int thelier, int nchunks, + int pointdim, int signallength, unsigned int useMaxNorm, int deviceid) { + cudaSetDevice(deviceid); + return cudaFindRSAll(h_bf_npointsrange, h_pointset, h_query, h_vecradius, + thelier, nchunks, pointdim, signallength, useMaxNorm); +} +#ifdef __cplusplus +} +#endif + + +#ifdef __cplusplus +extern "C" { +#endif +int findRadiiAlgorithm2(float *radii, const float *data, const int *indexes, + unsigned int k, unsigned int dim, unsigned int N) { + + unsigned int i, j; + + for (j = 0; j < N; j++) { + radii[j] = 0.0f; + for (i = 0; i < k; i++) { + float d = maxMetricPoints(data + j, data + indexes[j + i*N], dim, N); + if (d > radii[j]) { + radii[j] = d; + } + } + } + + return 1; + +} +#ifdef __cplusplus +} +#endif + + + +#ifdef __cplusplus +extern "C" { +#endif +int computeSumDigammas(float *sumDiGammas, int *nx, int *ny, unsigned int N) { + + int *d_nx, *d_ny; + float *d_sumDiGammas; + + checkCudaErrors( cudaMalloc((void **) &d_nx, N * sizeof(int)) ); + checkCudaErrors( cudaMalloc((void **) &d_ny, N * sizeof(int)) ); + checkCudaErrors( cudaMalloc((void **) &d_sumDiGammas, sizeof(float)) ); + + checkCudaErrors( cudaMemcpy(d_nx, nx, N*sizeof(int), cudaMemcpyHostToDevice) ); + checkCudaErrors( cudaMemcpy(d_ny, ny, N*sizeof(int), cudaMemcpyHostToDevice) ); + + unsigned int threads_per_block = 512; + dim3 n_blocks, n_threads; + n_blocks.x = ((N % threads_per_block) != 0) ? (N / threads_per_block + 1) : (N / threads_per_block); + n_threads.x = threads_per_block; + + printf("blocks = %i, threads = %i\n", n_blocks.x, n_threads.x); + + // reduce6<<>>(d_nx, d_ny, d_sumDiGammas, N); + reduce1<<>>(d_nx, d_ny, d_sumDiGammas, N); + + checkCudaErrors( cudaDeviceSynchronize() ); + + checkCudaErrors( cudaMemcpy(sumDiGammas, d_sumDiGammas, sizeof(float), cudaMemcpyDeviceToHost) ); + + cudaFree(d_nx); + cudaFree(d_ny); + cudaFree(d_sumDiGammas); + + *sumDiGammas = 0.0; + + return 1; + +} +#ifdef __cplusplus +} +#endif + + diff --git a/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnLibrary.h b/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnLibrary.h new file mode 100644 index 0000000..13b6ae7 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/cuda/gpuKnnLibrary.h @@ -0,0 +1,26 @@ +#ifndef GPUKNNLIBRARY_H +#define GPUKNNLIBRARY_H + +int cudaFindKnn(int* h_bf_indexes, float* h_bf_distances, float* h_pointset, + float* h_query, int kth, int thelier, int nchunks, int pointdim, + int signallength, unsigned int useMaxNorm); + +int cudaFindKnnSetGPU(int* h_bf_indexes, float* h_bf_distances, + float* h_pointset, float* h_query, int kth, int thelier, int nchunks, + int pointdim, int signallength, unsigned int useMaxNorm, int deviceid); + +int cudaFindRSAll(int* h_bf_npointsrange, float* h_pointset, + float* h_query, float* h_vecradius, int thelier, int nchunks, + int pointdim, int signallength, unsigned int useMaxNorm); + +int cudaFindRSAllSetGPU(int* h_bf_npointsrange, float* h_pointset, + float* h_query, float* h_vecradius, int thelier, int nchunks, + int pointdim, int signallength, unsigned int useMaxNorm, int deviceid); + +int findRadiiAlgorithm2(float *radii, const float *data, const int *indexes, + unsigned int k, unsigned int dim, unsigned int N); + +int computeSumDigammas(float *sumDiGammas, int *nx, int *ny, unsigned int N); + +#endif + diff --git a/java/source/infodynamics/measures/continuous/kraskov/cuda/helperfunctions.cu b/java/source/infodynamics/measures/continuous/kraskov/cuda/helperfunctions.cu new file mode 100644 index 0000000..bfc286c --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/cuda/helperfunctions.cu @@ -0,0 +1,230 @@ +//////////////////////////////////////////////////////////////////////////////// +// These are CUDA Helper functions + +#include +#include + +// This will output the proper CUDA error strings in the event that a CUDA host call returns an error +#define checkCudaErrors(err) __checkCudaErrors (err, __FILE__, __LINE__) + +inline void __checkCudaErrors(cudaError err, const char *file, const int line) +{ + if (cudaSuccess != err) + { + fprintf(stderr, "%s(%i) : CUDA Runtime API error %d: %s.\n",file, line, (int)err, cudaGetErrorString(err)); + exit(-1); + } +} + +// This will output the proper error string when calling cudaGetLastError +#define getLastCudaError(msg) __getLastCudaError (msg, __FILE__, __LINE__) + +inline void __getLastCudaError(const char *errorMessage, const char *file, const int line) +{ + cudaError_t err = cudaGetLastError(); + + if (cudaSuccess != err) + { + fprintf(stderr, "%s(%i) : getLastCudaError() CUDA error : %s : (%d) %s.\n", + file, line, errorMessage, (int)err, cudaGetErrorString(err)); + exit(-1); + } +} + +// General GPU Device CUDA Initialization +int gpuDeviceInit(int devID) +{ + int deviceCount; + checkCudaErrors(cudaGetDeviceCount(&deviceCount)); + + if (deviceCount == 0) + { + fprintf(stderr, "gpuDeviceInit() CUDA error: no devices supporting CUDA.\n"); + exit(-1); + } + + if (devID < 0) + { + devID = 0; + } + + if (devID > deviceCount-1) + { + fprintf(stderr, "\n"); + fprintf(stderr, ">> %d CUDA capable GPU device(s) detected. <<\n", deviceCount); + fprintf(stderr, ">> gpuDeviceInit (-device=%d) is not a valid GPU device. <<\n", devID); + fprintf(stderr, "\n"); + return -devID; + } + + cudaDeviceProp deviceProp; + checkCudaErrors(cudaGetDeviceProperties(&deviceProp, devID)); + + if(!deviceProp.deviceOverlap){ + printf("\n Device will not handle overlaps, so no speed up from streams \n"); + } + if (deviceProp.major < 1) + { + fprintf(stderr, "gpuDeviceInit(): GPU device does not support CUDA.\n"); + exit(-1); + } + + checkCudaErrors(cudaSetDevice(devID)); + printf("gpuDeviceInit() CUDA Device [%d]: \"%s\n", devID, deviceProp.name); + + + return devID; +} + +#ifndef MAX +#define MAX(a,b) (a > b ? a : b) +#endif + +// Beginning of GPU Architecture definitions +inline int _ConvertSMVer2Cores(int major, int minor) +{ + // Defines for GPU Architecture types (using the SM version to determine the # of cores per SM + typedef struct + { + int SM; // 0xMm (hexidecimal notation), M = SM Major version, and m = SM minor version + int Cores; + } sSMtoCores; + + sSMtoCores nGpuArchCoresPerSM[] = + { + { 0x10, 8 }, // Tesla Generation (SM 1.0) G80 class + { 0x11, 8 }, // Tesla Generation (SM 1.1) G8x class + { 0x12, 8 }, // Tesla Generation (SM 1.2) G9x class + { 0x13, 8 }, // Tesla Generation (SM 1.3) GT200 class + { 0x20, 32 }, // Fermi Generation (SM 2.0) GF100 class + { 0x21, 48 }, // Fermi Generation (SM 2.1) GF10x class + { 0x30, 192}, // Fermi Generation (SM 3.0) GK10x class + { -1, -1 } + }; + + int index = 0; + + while (nGpuArchCoresPerSM[index].SM != -1) + { + if (nGpuArchCoresPerSM[index].SM == ((major << 4) + minor)) + { + return nGpuArchCoresPerSM[index].Cores; + } + + index++; + } + + printf("MapSMtoCores undefined SM %d.%d is undefined (please update to the latest SDK)!\n", major, minor); + return -1; +} +// end of GPU Architecture definitions + + +// This function returns the best GPU (with maximum GFLOPS) +int gpuGetMaxGflopsDeviceId() +{ + int current_device = 0, sm_per_multiproc = 0; + int max_compute_perf = 0, max_perf_device = 0; + int device_count = 0, best_SM_arch = 0; + cudaDeviceProp deviceProp; + cudaGetDeviceCount(&device_count); + + // Find the best major SM Architecture GPU device + while (current_device < device_count) + { + cudaGetDeviceProperties(&deviceProp, current_device); + + if (deviceProp.major > 0 && deviceProp.major < 9999) + { + best_SM_arch = MAX(best_SM_arch, deviceProp.major); + } + + current_device++; + } + + // Find the best CUDA capable GPU device + current_device = 0; + + while (current_device < device_count) + { + cudaGetDeviceProperties(&deviceProp, current_device); + + if (deviceProp.major == 9999 && deviceProp.minor == 9999) + { + sm_per_multiproc = 1; + } + else + { + sm_per_multiproc = _ConvertSMVer2Cores(deviceProp.major, deviceProp.minor); + } + + int compute_perf = deviceProp.multiProcessorCount * sm_per_multiproc * deviceProp.clockRate; + + if (compute_perf > max_compute_perf) + { + // If we find GPU with SM major > 2, search only these + if (best_SM_arch > 2) + { + // If our device==dest_SM_arch, choose this, or else pass + if (deviceProp.major == best_SM_arch) + { + max_compute_perf = compute_perf; + max_perf_device = current_device; + } + } + else + { + max_compute_perf = compute_perf; + max_perf_device = current_device; + } + } + + ++current_device; + } + + return max_perf_device; +} + + +// Initialization code to find the best CUDA Device +int findCudaDevice(int argc, const char **argv) +{ + cudaDeviceProp deviceProp; + int devID = 0; + + // If the command-line has a device number specified, use it + /*if (checkCmdLineFlag(argc, argv, "device")) + { + devID = getCmdLineArgumentInt(argc, argv, "device="); + + if (devID < 0) + { + printf("Invalid command line parameter\n "); + exit(-1); + } + else + { + devID = gpuDeviceInit(devID); + + if (devID < 0) + { + printf("exiting...\n"); + //shrQAFinishExit(argc, (const char **)argv, QA_FAILED); + exit(-1); + } + } + } + else + {*/ + // Otherwise pick the device with highest Gflops/s + devID = gpuGetMaxGflopsDeviceId(); + checkCudaErrors(cudaSetDevice(devID)); + checkCudaErrors(cudaGetDeviceProperties(&deviceProp, devID)); + printf("GPU Device %d: \"%s\" with compute capability %d.%d\n\n", devID, deviceProp.name, deviceProp.major, deviceProp.minor); + //} + + return devID; +} +// end of CUDA Helper Functions + + diff --git a/java/source/infodynamics/measures/continuous/kraskov/cuda/kraskovCuda.c b/java/source/infodynamics/measures/continuous/kraskov/cuda/kraskovCuda.c new file mode 100644 index 0000000..439499a --- /dev/null +++ b/java/source/infodynamics/measures/continuous/kraskov/cuda/kraskovCuda.c @@ -0,0 +1,316 @@ +#include +#include +#include + +#include "digamma.c" +#include "gpuKnnLibrary.h" + +#define check(ans) { _check((ans), __FILE__, __LINE__); } + +#ifdef __cplusplus +extern "C" { +#endif +/* + * Class: infodynamics_measures_continuous_kraskov_MutualInfoCalculatorMultiVariateKraskov + * Method: MIKraskov + * Signature: (I[DI[DIIZZZ)[D + */ +JNIEXPORT jdoubleArray JNICALL + Java_infodynamics_measures_continuous_kraskov_MutualInfoCalculatorMultiVariateKraskov_MIKraskov( + JNIEnv *env, jobject thisObj, jint j_N, + jobjectArray j_sourceArray, jint j_dimx, + jobjectArray j_destArray, jint j_dimy, + jint j_k, jboolean j_returnLocals, jboolean j_useMaxNorm, jboolean j_isAlgorithm1) { + + // Declare variables + // Variables starting with j_ are Java types + float *source, *dest, *pointset, *radii, *distances; + int *nx, *ny, *indexes; + jdouble *jd_source, *jd_dest, *result; + unsigned int N, dimx, dimy, dims, k, returnLocals, useMaxNorm, isAlgorithm1, result_size; + float sumNx, sumNy, sumDiGammas; + unsigned int i, j; + long double digammaK, digammaN; + int thelier = 0, nchunks = 1, err; + const int verbose = 0; + + if (verbose) printf("Verbose is on.\n"); + + // Check that incoming data has correct size + jsize sourceLength = (*env)->GetArrayLength(env, j_sourceArray); + jsize destLength = (*env)->GetArrayLength(env, j_destArray); + + if (sourceLength != j_N || destLength != j_N) { + cudaDeviceReset(); + jclass Exception = (*env)->FindClass(env, "java/lang/Exception"); + (*env)->ThrowNew(env, Exception, "Data has wrong length."); + } + + // 1. Load data into GPU + // ===================== + // 1a. Copy atomic variables from Java + N = j_N; + k = j_k; + dimx = j_dimx; + dimy = j_dimy; + dims = dimx + dimy; + returnLocals = j_returnLocals ? 1 : 0; + useMaxNorm = j_useMaxNorm ? 1 : 0; + isAlgorithm1 = j_isAlgorithm1 ? 1 : 0; + + if (verbose > 0) { + printf("Loading data from Java to C.\n"); + printf("N = %i, dimx = %i, dimy = %i\n", N, dimx, dimy); + } + + // 1b. Get arrays from Java + source = (float *) malloc(N * dimx * sizeof(float)); + dest = (float *) malloc(N * dimy * sizeof(float)); + pointset = (float *) malloc(N * dims * sizeof(float)); + + if (NULL == source || NULL == dest || NULL == pointset) { + printf("Error allocating data.\n"); + cudaDeviceReset(); + jclass Exception = (*env)->FindClass(env, "java/lang/Exception"); + (*env)->ThrowNew(env, Exception, "Error allocating data."); + } + + for (i = 0; i < N; i++) { + jdoubleArray j_sourceRow = (*env)->GetObjectArrayElement(env, j_sourceArray, i); + jdoubleArray j_destRow = (*env)->GetObjectArrayElement(env, j_destArray, i); + + jdouble *sourceRow = (*env)->GetDoubleArrayElements(env, j_sourceRow, NULL); + jdouble *destRow = (*env)->GetDoubleArrayElements(env, j_destRow, NULL); + + // Data in java are doubles, but GPUs need floats. + // We have to cast them manually, so we can't memcopy + + // The following two for-loops get two matrices in T-by-D indexing (i.e. + // first dimension is time, second is variable) and return the data in + // column-major form + for (j = 0; j < dimx; j++) { + register float x = (float) sourceRow[j]; + source[N*j + i] = x; + pointset[N*j + i] = x; + } + + for (j = 0; j < dimy; j++) { + register float x = (float) destRow[j]; + dest[N*j + i] = x; + pointset[N*dimx + N*j + i] = x; + } + + // // The following two for-loops get the same matrices but return arrays + // // in row-major form. Code is left here for future development purposes. + // for (j = 0; j < dimx; j++) { + // register float x = (float) sourceRow[j]; + // source[dimx*i + j] = x; + // pointset[dims*i + j] = x; + // } + + // for (j = 0; j < dimy; j++) { + // register float x = (float) destRow[j]; + // dest[dimy*i + j] = x + // pointset[i*dims + dimx + j] = x; + // } + + + (*env)->ReleaseDoubleArrayElements(env, j_sourceRow, sourceRow, 0); + (*env)->ReleaseDoubleArrayElements(env, j_destRow, destRow, 0); + + (*env)->DeleteLocalRef(env, j_sourceRow); + (*env)->DeleteLocalRef(env, j_destRow); + + // FIXME: I'm not entirely sure I'm freeing all the memory here. I should + // check for memory leaks more carefully. + + } + + if (verbose > 1) { + printf("Flat source array: "); + for (i = 0; i < N*dimx; i++) printf("%f, ", source[i]); + printf("\n"); + printf("Flat dest array: "); + for (i = 0; i < N*dimy; i++) printf("%f, ", dest[i]); + printf("\n"); + printf("Flat pointset: "); + for (i = 0; i < N*dims; i++) printf("%f, ", pointset[i]); + printf("\n"); + } + + + // 2. Get indices and distances to the k NN of each point + // ====================== + + indexes = (int *) malloc(N * k * sizeof(int)); + distances = (float *) malloc(N * k * sizeof(float)); + + err = cudaFindKnn(indexes, distances, pointset, pointset, k, + thelier, nchunks, dims, N, useMaxNorm); + if (err != 1) { + cudaDeviceReset(); + jclass Exception = (*env)->FindClass(env, "java/lang/Exception"); + (*env)->ThrowNew(env, Exception, "Cuda error finding nearest neighbours."); + } + + if (verbose > 1) { + printf("Distances: \n"); + for (j = 0; j < N; j++) { + for (i = 0; i < k; i++) { + printf("%f\t", distances[j + i*N]); + } + printf("\n"); + } + printf("Indexes: \n"); + for (j = 0; j < N; j++) { + for (i = 0; i < k; i++) { + printf("%i\t", indexes[j + i*N]); + } + printf("\n"); + } + } + + + // 3. Find distance R from each point to its k-th neighbour in joint space + // ====================== + + if (verbose > 0) { + printf("Getting radii from calculated dists.\n"); + } + + radii = (float *) malloc(N * sizeof(float)); + for (i = 0; i < N; i++) { + radii[i] = distances[N*(k-1) + i]; + } + + if (verbose > 1) { + for (i = 0; i < N; i++) { + printf("%f\t", radii[i]); + } + printf("\n"); + } + + + // 4. Count points strictly within R in the X-space + // ====================== + + if (verbose > 0) { + printf("Counting points within radius.\n"); + } + + // If we're using algorithm 2 then we need the radius in the marginal space, + // not the joint (as calculated above) + if (!isAlgorithm1) { + findRadiiAlgorithm2(radii, source, indexes, k, dimx, N); + } + + nx = (int *) malloc(N * sizeof(int)); + err = cudaFindRSAll(nx, source, source, radii, thelier, nchunks, dimx, N, useMaxNorm); + if (err != 1) { + cudaDeviceReset(); + jclass Exception = (*env)->FindClass(env, "java/lang/Exception"); + (*env)->ThrowNew(env, Exception, "Cuda error during source range search."); + } + + if (verbose > 0) { + printf("X count: "); + for (i = 0; i < N; i++) printf("%i, ", nx[i]); + printf("\n\n"); + } + + // 4. Count points strictly within R in the X-space + // ====================== + + if (!isAlgorithm1) { + findRadiiAlgorithm2(radii, dest, indexes, k, dimy, N); + } + + ny = (int *) malloc(N * sizeof(int)); + err = cudaFindRSAll(ny, dest, dest, radii, thelier, nchunks, dimy, N, useMaxNorm); + if (err != 1) { + cudaDeviceReset(); + jclass Exception = (*env)->FindClass(env, "java/lang/Exception"); + (*env)->ThrowNew(env, Exception, "Cuda error during dest range search."); + } + + if (verbose > 0) { + printf("Y count: "); + for (i = 0; i < N; i++) printf("%i, ", ny[i]); + printf("\n\n"); + } + + + // 6. Set locals or digammas for return + // ====================== + if (verbose > 0) { + printf("Start calculating digammas\n"); + } + + if (returnLocals) { + result = (jdouble *) malloc(N * sizeof(jdouble)); + result_size = N; + digammaK = cpuDigamma(k); + digammaN = cpuDigamma(N); + for (i = 0; i < N; i++) { + result[i] = digammaK + digammaN - cpuDigamma(nx[i] + 1) - cpuDigamma(ny[i] + 1); + } + + } else { + + result = (jdouble *) malloc(3 * sizeof(jdouble)); + result_size = 3; + + sumNx = 0; + sumNy = 0; + sumDiGammas = 0; + for (i = 0; i < N; i++) { + sumNx += nx[i]; + sumNy += ny[i]; + sumDiGammas += cpuDigamma(nx[i] + 1) + cpuDigamma(ny[i] + 1); + + } + + result[0] = (jdouble) sumDiGammas; + result[1] = (jdouble) sumNx; + result[2] = (jdouble) sumNy; + + } + + + // 7. Free memory and return + // ========================= + + if (verbose > 0) { + printf("Free memory\n"); + } + + if (source) free(source); + if (dest) free(dest); + if (pointset) free(pointset); + if (radii) free(radii); + if (distances) free(distances); + if (indexes) free(indexes); + if (nx) free(nx); + if (ny) free(ny); + + // Set Java array for return + jdoubleArray outJNIArray = (*env)->NewDoubleArray(env, result_size); // allocate + if (NULL == outJNIArray) { + cudaDeviceReset(); + jclass Exception = (*env)->FindClass(env, "java/lang/Exception"); + (*env)->ThrowNew(env, Exception, "Error creating return array."); + } + (*env)->SetDoubleArrayRegion(env, outJNIArray, 0 , result_size, result); // copy + + if (result) free(result); + + return outJNIArray; + + +} // End of function MIKraskov1 + + +#ifdef __cplusplus +} +#endif +