25 #include <queso/InfoTheory.h>
34 void distANN_XY(
const ANNpointArray dataX,
const ANNpointArray dataY,
36 unsigned int dimX,
unsigned int dimY,
37 unsigned int xN,
unsigned int yN,
38 unsigned int k,
double eps )
46 nnIdx =
new ANNidx[ k+1 ];
47 nnDist =
new ANNdist[ k+1 ];
48 kdTree =
new ANNkd_tree( dataY, yN, dimY );
51 for(
unsigned int i = 0; i < xN ; i++ )
53 kdTree->annkSearch( dataX[ i ], k+1, nnIdx, nnDist, eps );
55 double my_dist = nnDist[ k ];
61 ANNdistArray nnDist_tmp =
new ANNdist[ yN ];
62 ANNidxArray nnIdx_tmp =
new ANNidx[ yN ];
63 kdTree->annkSearch( dataX[ i ], yN, nnIdx_tmp, nnDist_tmp, eps );
65 for(
unsigned int my_k = k + 1; my_k < yN; ++my_k )
66 if( nnDist_tmp[ my_k ] > 0.0 )
68 my_dist = nnDist_tmp[ my_k ];
75 distsXY[ i ] = my_dist;
92 void normalizeANN_XY( ANNpointArray dataXY,
unsigned int dimXY,
93 ANNpointArray dataX,
unsigned int dimX,
94 ANNpointArray dataY,
unsigned int dimY,
98 ANNpoint meanXY, stdXY;
100 meanXY = annAllocPt(dimXY);
101 stdXY = annAllocPt(dimXY);
104 for(
unsigned int i = 0; i < N; i++ ) {
105 for(
unsigned int j = 0; j < dimX; j++ ) {
106 meanXY[ j ] += dataXY[ i ][ j ];
108 for(
unsigned int j = 0; j < dimY; j++ ) {
109 meanXY[ dimX + j ] += dataXY[ i ][ dimX + j ];
112 for(
unsigned int j = 0; j < dimXY; j++ ) {
113 meanXY[ j ] = meanXY[ j ] / (double)N;
117 for(
unsigned int i = 0; i < N; i++ ) {
118 for(
unsigned int j = 0; j < dimXY; j++ ) {
119 stdXY[ j ] += pow( dataXY[ i ][ j ] - meanXY[ j ], 2.0 );
122 for(
unsigned int j = 0; j < dimXY; j++ ) {
123 stdXY[ j ] = sqrt( stdXY[ j ] / ((
double)N-1.0) );
135 for(
unsigned int i = 0; i < N; i++ ) {
137 for(
unsigned int j = 0; j < dimXY; j++ ) {
138 dataXY[ i ][ j ] = ( dataXY[ i ][ j ] - meanXY[ j ] ) / stdXY[ j ];
142 for(
unsigned int j = 0; j < dimX; j++ ) {
143 dataX[ i ][ j ] = dataXY[ i ][ j ];
145 for(
unsigned int j = 0; j < dimY; j++ ) {
146 dataY[ i ][ j ] = dataXY[ i ][ dimX + j ];
155 double computeMI_ANN( ANNpointArray dataXY,
156 unsigned int dimX,
unsigned int dimY,
157 unsigned int k,
unsigned int N,
double eps )
160 ANNpointArray dataX, dataY;
166 unsigned int dimXY = dimX + dimY;
169 dataX = annAllocPts(N,dimX);
170 dataY = annAllocPts(N,dimY);
171 distsXY =
new double[N];
174 normalizeANN_XY( dataXY, dimXY, dataX, dimX, dataY, dimY, N);
177 kdTreeX =
new ANNkd_tree( dataX, N, dimX );
178 kdTreeY =
new ANNkd_tree( dataY, N, dimY );
179 distANN_XY( dataXY, dataXY, distsXY, dimXY, dimXY, N, N, k, eps );
182 double marginal_contrib = 0.0;
183 for(
unsigned int i = 0; i < N; i++ ) {
185 int no_pts_X = kdTreeX->annkFRSearch( dataX[ i ], distsXY[ i ], 0, NULL, NULL, eps);
186 int no_pts_Y = kdTreeY->annkFRSearch( dataY[ i ], distsXY[ i ], 0, NULL, NULL, eps);
188 marginal_contrib += gsl_sf_psi_int( no_pts_X+1 ) + gsl_sf_psi_int( no_pts_Y+1 );
190 MI_est = gsl_sf_psi_int( k ) + gsl_sf_psi_int( N ) - marginal_contrib / (double)N;
196 annDeallocPts( dataX );
197 annDeallocPts( dataY );
207 template<
template <
class P_V,
class P_M>
class RV,
class P_V,
class P_M>
208 double estimateMI_ANN(
const RV<P_V,P_M>& jointRV,
209 const unsigned int xDimSel[],
unsigned int dimX,
210 const unsigned int yDimSel[],
unsigned int dimY,
211 unsigned int k,
unsigned int N,
double eps )
213 ANNpointArray dataXY;
216 unsigned int dimXY = dimX + dimY;
219 dataXY = annAllocPts(N,dimXY);
222 P_V smpRV( jointRV.imageSet().vectorSpace().zeroVector() );
223 for(
unsigned int i = 0; i < N; i++ ) {
225 jointRV.realizer().realization( smpRV );
228 for(
unsigned int j = 0; j < dimX; j++ ) {
229 dataXY[ i ][ j ] = smpRV[ xDimSel[j] ];
231 for(
unsigned int j = 0; j < dimY; j++ ) {
232 dataXY[ i ][ dimX + j ] = smpRV[ yDimSel[j] ];
237 MI_est = computeMI_ANN( dataXY,
242 annDeallocPts( dataXY );
251 template<
class P_V,
class P_M,
252 template <
class P_V,
class P_M>
class RV_1,
253 template <
class P_V,
class P_M>
class RV_2>
254 double estimateMI_ANN(
const RV_1<P_V,P_M>& xRV,
255 const RV_2<P_V,P_M>& yRV,
256 const unsigned int xDimSel[],
unsigned int dimX,
257 const unsigned int yDimSel[],
unsigned int dimY,
258 unsigned int k,
unsigned int N,
double eps )
260 ANNpointArray dataXY;
263 unsigned int dimXY = dimX + dimY;
266 dataXY = annAllocPts(N,dimXY);
269 P_V smpRV_x( xRV.imageSet().vectorSpace().zeroVector() );
270 P_V smpRV_y( yRV.imageSet().vectorSpace().zeroVector() );
272 for(
unsigned int i = 0; i < N; i++ ) {
274 xRV.realizer().realization( smpRV_x );
275 yRV.realizer().realization( smpRV_y );
278 for(
unsigned int j = 0; j < dimX; j++ ) {
279 dataXY[ i ][ j ] = smpRV_x[ xDimSel[j] ];
281 for(
unsigned int j = 0; j < dimY; j++ ) {
282 dataXY[ i ][ dimX + j ] = smpRV_y[ yDimSel[j] ];
287 MI_est = computeMI_ANN( dataXY,
292 annDeallocPts( dataXY );
301 template <
class P_V,
class P_M,
302 template <
class P_V,
class P_M>
class RV_1,
303 template <
class P_V,
class P_M>
class RV_2>
304 double estimateKL_ANN( RV_1<P_V,P_M>& xRV,
306 unsigned int xDimSel[],
unsigned int dimX,
307 unsigned int yDimSel[],
unsigned int dimY,
308 unsigned int xN,
unsigned int yN,
309 unsigned int k,
double eps )
323 dataX = annAllocPts( xN, dimX );
324 dataY = annAllocPts( yN, dimY );
325 distsX =
new double[xN];
326 distsXY =
new double[xN];
329 P_V xSmpRV( xRV.imageSet().vectorSpace().zeroVector() );
330 for(
unsigned int i = 0; i < xN; i++ ) {
332 xRV.realizer().realization( xSmpRV );
334 for(
unsigned int j = 0; j < dimX; j++ ) {
335 dataX[ i ][ j ] = xSmpRV[ xDimSel[j] ];
340 P_V ySmpRV( yRV.imageSet().vectorSpace().zeroVector() );
341 for(
unsigned int i = 0; i < yN; i++ ) {
343 yRV.realizer().realization( ySmpRV );
345 for(
unsigned int j = 0; j < dimY; j++ ) {
346 dataY[ i ][ j ] = ySmpRV[ yDimSel[j] ];
351 distANN_XY( dataX, dataX, distsX, dimX, dimX, xN, xN, k+1, eps );
352 distANN_XY( dataX, dataY, distsXY, dimX, dimY, xN, yN, k, eps );
355 double sum_log_ratio = 0.0;
356 for(
unsigned int i = 0; i < xN; i++ )
358 sum_log_ratio += log( distsXY[i] / distsX[i] );
360 KL_est = (double)dimX/(
double)xN * sum_log_ratio + log( (
double)yN / ((
double)xN-1.0 ) );
363 annDeallocPts( dataX );
364 annDeallocPts( dataY );
376 template <
class P_V,
class P_M,
377 template <
class P_V,
class P_M>
class RV_1,
378 template <
class P_V,
class P_M>
class RV_2>
379 double estimateCE_ANN( RV_1<P_V,P_M>& xRV,
381 unsigned int xDimSel[],
unsigned int dimX,
382 unsigned int yDimSel[],
unsigned int dimY,
383 unsigned int xN,
unsigned int yN,
384 unsigned int k,
double eps )
398 dataX = annAllocPts( xN, dimX );
399 dataY = annAllocPts( yN, dimY );
400 distsXY =
new double[xN];
401 kdTree =
new ANNkd_tree( dataY, yN, dimY );
404 P_V xSmpRV( xRV.imageSet().vectorSpace().zeroVector() );
405 for(
unsigned int i = 0; i < xN; i++ ) {
407 xRV.realizer().realization( xSmpRV );
409 for(
unsigned int j = 0; j < dimX; j++ ) {
410 dataX[ i ][ j ] = xSmpRV[ xDimSel[j] ];
415 P_V ySmpRV( yRV.imageSet().vectorSpace().zeroVector() );
416 for(
unsigned int i = 0; i < yN; i++ ) {
418 yRV.realizer().realization( ySmpRV );
420 for(
unsigned int j = 0; j < dimY; j++ ) {
421 dataY[ i ][ j ] = ySmpRV[ yDimSel[j] ];
426 distANN_XY( dataX, dataY, distsXY, dimX, dimY, xN, yN, k, eps );
427 kdTree =
new ANNkd_tree( dataY, yN, dimY );
430 double sum_log = 0.0;
431 for(
unsigned int i = 0; i < xN; i++ )
433 sum_log += log( distsXY[i] );
435 CE_est = (double)dimX/(
double)xN * sum_log + log( (
double)yN ) - gsl_sf_psi_int( k );
438 annDeallocPts( dataX );
439 annDeallocPts( dataY );
447 #endif // QUESO_HAS_ANN
#define queso_error_msg(msg)