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

Separated core MI function to different files (and hacked Makefile).

This commit is contained in:
Pedro Martinez Mediano 2017-05-19 16:53:26 +10:00
parent 4b665d525f
commit 6e9ade6542
4 changed files with 145 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
java/source/infodynamics/measures/continuous/kraskov/cuda/Makefile
View File

@ -32,7 +32,7 @@ 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
${NVCC} ${NVCCFLAGS} ${INCLUDES} -Xcompiler -fPIC -shared -o libKraskov.so gpuMILibrary.c kraskovCuda.c ${LDFLAGS} -lgpuKnnLibrary
test: libKraskov.so
${GCC} -std=c++11 -o unittest unittests.cpp -L. -lKraskov

129
java/source/infodynamics/measures/continuous/kraskov/cuda/gpuMILibrary.c Normal file
View File

@ -0,0 +1,129 @@
#include <stdlib.h>
#include <stdio.h>
#include "gpuMILibrary.h"
#include "gpuKnnLibrary.h"
#include "digamma.c"
jidt_error_t MIKraskov_C(int N, float *source, int dimx, float *dest, int dimy,
int k, int thelier, int nchunks, int returnLocals, int useMaxNorm,
int isAlgorithm1, float *result) {
int dims = dimx + dimy;
int err;
int i, j;
const int verbose = 0;
float *pointset = NULL, *distances = NULL, *radii = NULL;
int *indexes = NULL, *nx = NULL, *ny = NULL;
pointset = (float *) malloc(N * dims * sizeof(float));
for (i = 0; i < N; i++) {
for (j = 0; j < dimx; j++) {
register int idx = N*j + i;
pointset[idx] = source[idx];
}
for (j = 0; j < dimy; j++) {
register int idx = N*j + i;
pointset[N*dimx + idx] = dest[idx];
}
}
indexes = (int *) malloc(N * k * sizeof(int));
distances = (float *) malloc(N * k * sizeof(float));
if (!cudaFindKnn(indexes, distances, pointset, pointset, k,
thelier, nchunks, dims, N, useMaxNorm)) {
device_reset();
return JIDT_ERROR;
}
// 3. Find distance R from each point to its k-th neighbour in joint space
// ======================
radii = (float *) malloc(N * sizeof(float));
for (i = 0; i < N; i++) {
radii[i] = distances[N*(k-1) + i];
}
// 4. Count points strictly within R in the X-space
// ======================
// 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));
if (!cudaFindRSAll(nx, source, source, radii, thelier, nchunks, dimx, N, useMaxNorm)) {
device_reset();
return JIDT_ERROR;
}
if (!isAlgorithm1) {
findRadiiAlgorithm2(radii, dest, indexes, k, dimy, N);
}
ny = (int *) malloc(N * sizeof(int));
if (!cudaFindRSAll(ny, dest, dest, radii, thelier, nchunks, dimy, N, useMaxNorm)) {
device_reset();
return JIDT_ERROR;
}
// 6. Set locals, surrogates or digammas for return
// ======================
float sumNx, sumNy, sumDiGammas, digammaK, digammaN;
int result_size;
if (nchunks > 1) {
float digammaK = cpuDigamma(k);
float digammaN = cpuDigamma(N);
float trialLength = N/((float) nchunks);
for (int ii = 0; ii < nchunks; ii++) {
float sumDiGammas = 0;
for (i = 0; i < trialLength; i++) {
sumDiGammas += cpuDigamma(nx[i] + 1) + cpuDigamma(ny[i] + 1);
}
result[ii] = digammaK - sumDiGammas/trialLength + digammaN;
}
} else if (returnLocals) {
float digammaK = cpuDigamma(k);
float digammaN = cpuDigamma(N);
for (i = 0; i < N; i++) {
result[i] = digammaK + digammaN - cpuDigamma(nx[i] + 1) - cpuDigamma(ny[i] + 1);
}
} else {
float sumNx = 0;
float sumNy = 0;
float 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] = sumDiGammas;
result[1] = sumNx;
result[2] = sumNy;
}
err = JIDT_SUCCESS;
FREE(indexes);
FREE(distances);
FREE(nx);
FREE(ny);
FREE(distances);
FREE(radii);
FREE(pointset);
return err;
}

13
java/source/infodynamics/measures/continuous/kraskov/cuda/gpuMILibrary.h Normal file
View File

@ -0,0 +1,13 @@
#ifndef GPUMILIBRARY_H
#define GPUMILIBRARY_H
#define FREE(x) { if (x) free(x); x = NULL; }
typedef enum { JIDT_SUCCESS, JIDT_ERROR } jidt_error_t;
jidt_error_t MIKraskov_C(int N, float *source, int dimx, float *dest, int dimy,
int k, int thelier, int nchunks, int returnLocals, int useMaxNorm,
int isAlgorithm1, float *result);
#endif

111
java/source/infodynamics/measures/continuous/kraskov/cuda/kraskovCuda.c
View File

@ -2,116 +2,9 @@
#include <stdlib.h>
#include <stdio.h>
#include "digamma.c"
#include "gpuKnnLibrary.h"
#define check(ans) { _check((ans), __FILE__, __LINE__); }
typedef enum { JIDT_SUCCESS, JIDT_ERROR } jidt_error_t;
jidt_error_t MIKraskov_C(int N, float *source, int dimx, float *dest, int dimy,
int k, int thelier, int nchunks, int returnLocals, int useMaxNorm,
int isAlgorithm1, int nb_surrogates, float *result) {
int dims = dimx + dimy;
int err;
int i, j;
const int verbose = 0;
nb_surrogates = 0;
float *pointset = (float *) malloc(N * dims * sizeof(float));
for (i = 0; i < N; i++) {
for (j = 0; j < dimx; j++) {
register int idx = N*j + i;
pointset[idx] = source[idx];
}
for (j = 0; j < dimy; j++) {
register int idx = N*j + i;
pointset[N*dimx + idx] = dest[idx];
}
}
int *indexes = (int *) malloc(N * k * sizeof(int));
float *distances = (float *) malloc(N * k * sizeof(float));
if (!cudaFindKnn(indexes, distances, pointset, pointset, k,
thelier, nchunks, dims, N, useMaxNorm)) {
device_reset();
return JIDT_ERROR;
}
// 3. Find distance R from each point to its k-th neighbour in joint space
// ======================
float *radii = (float *) malloc(N * sizeof(float));
for (i = 0; i < N; i++) {
radii[i] = distances[N*(k-1) + i];
}
// 4. Count points strictly within R in the X-space
// ======================
// 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);
}
int *nx = (int *) malloc(N * sizeof(int));
if (!cudaFindRSAll(nx, source, source, radii, thelier, nchunks, dimx, N, useMaxNorm)) {
device_reset();
return JIDT_ERROR;
}
if (!isAlgorithm1) {
findRadiiAlgorithm2(radii, dest, indexes, k, dimy, N);
}
int *ny = (int *) malloc(N * sizeof(int));
if (!cudaFindRSAll(ny, dest, dest, radii, thelier, nchunks, dimy, N, useMaxNorm)) {
device_reset();
return JIDT_ERROR;
}
// 6. Set locals or digammas for return
// ======================
float sumNx, sumNy, sumDiGammas, digammaK, digammaN;
int result_size;
if (nb_surrogates > 0) {
printf("Surrogates not supported yet\n");
return JIDT_ERROR;
} else if (returnLocals) {
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 {
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] = (float) sumDiGammas;
result[1] = (float) sumNx;
result[2] = (float) sumNy;
}
return JIDT_SUCCESS;
}
#include "gpuMILibrary.h"
#ifdef __cplusplus
extern "C" {
@ -232,7 +125,7 @@ JNIEXPORT jdoubleArray JNICALL
float *result = (float *) malloc(resultSize * sizeof(float));
jidt_error_t ret = MIKraskov_C(N, source, dimx, dest, dimy,
k, thelier, nchunks, returnLocals, useMaxNorm,
isAlgorithm1, nb_surrogates, result);
isAlgorithm1, result);
if (JIDT_ERROR == ret) {
jclass Exception = (*env)->FindClass(env, "java/lang/Exception");