Commits

Zoltan Szabo committed 8b71da6

Edgeworth expansion based Shannon entropy estimator: accelerated (C++); see 'Edgeworth_t1_t2_t3.cpp'.

Comments (0)

Files changed (19)

+v0.17 (Nov 6, 2012):
+-Edgeworth expansion based Shannon entropy estimator: accelerated (C++); see 'Edgeworth_t1_t2_t3.cpp'.
 -'Tsallis entropy <- Renyi entropy' meta estimator: added; see 'HTsallis_HRenyi_initialization.m', 'HTsallis_HRenyi_estimation.m'. 
 
 v0.16 (Nov 2, 2012):

code/H_I_D/base_estimators/DRenyi_kNN_k_initialization.m

             %co.k = 3;%k-nearest neighbors
             %co.epsi = 0; %=0: exact kNN; >0: approximate kNN, the true (not squared) distances can not exceed the real distance more than a factor of (1+epsi).
             
-    co.alpha = 0.99; %The Renyi divergence equals to the Kullback-Leibler divergence in limit, i.e., D_{R,alpha} -> KL, provided that alpha ->1.
+    co.alpha = 0.99; %The Renyi divergence (D_{R,alpha}) equals to the Kullback-Leibler divergence (D) in limit: D_{R,alpha} -> D, as alpha -> 1.
 
 %initialize the ann wrapper in Octave, if needed:
     initialize_Octave_ann_wrapper_if_needed(co.kNNmethod);

code/H_I_D/base_estimators/DTsallis_kNN_k_initialization.m

             %co.k = 3;%k-nearest neighbors
             %co.epsi = 0; %=0: exact kNN; >0: approximate kNN, the true (not squared) distances can not exceed the real distance more than a factor of (1+epsi).
 				
-    co.alpha = 0.99; %The Tsallis divergence equals to the Kullback-Leibler divergence in limit, i.e., D_{T,alpha} -> KL, provided that alpha ->1.
+    co.alpha = 0.99; %The Tsallis divergence (D_{T,alpha}) equals to the Kullback-Leibler divergence (D) in limit: D_{T,alpha} -> D, as alpha -> 1.
 
 %initialize the ann wrapper in Octave, if needed:
     initialize_Octave_ann_wrapper_if_needed(co.kNNmethod);

code/H_I_D/base_estimators/HRenyi_GSF_initialization.m

     %Method to compute the geodesic spanning forest:
         co.GSFmethod = 'MatlabBGL_Kruskal';
 
-    co.alpha = 0.99; %The Renyi entropy equals to the Shannon differential entropy, in limit, i.e., Renyi=H_{R,alpha} -> Shannon=H, provided that alpha ->1.
+    co.alpha = 0.99; %The Renyi entropy (H_{R,alpha}) equals to the Shannon differential entropy (H) in limit: H_{R,alpha} -> H, as alpha -> 1.
     co.additive_constant_is_relevant = 0; %1:additive constant is relevant (you can precompute it via 'estimate_HRenyi_constant.m'), 0:not relevant
     
 %initialize the ann wrapper in Octave, if needed:

code/H_I_D/base_estimators/HRenyi_MST_initialization.m

     co.mult = mult;
     
 %other fields:
-    co.alpha = 0.99; %The Renyi entropy equals to the Shannon differential entropy, in limit, i.e., Renyi=H_{R,alpha} -> Shannon=H, provided that alpha ->1.
+    co.alpha = 0.99; %The Renyi entropy (H_{R,alpha}) equals to the Shannon differential entropy (H) in limit: H_{R,alpha} -> H, as alpha -> 1.
     %Possibilites for the MST (minimum spanning tree) method:
         co.MSTmethod = 'MatlabBGL_Prim';
         %co.MSTmethod = 'MatlabBGL_Kruskal';

code/H_I_D/base_estimators/HRenyi_kNN_1tok_estimation.m

     L = compute_length_HRenyi_kNN_1tok(Y,co);
 
 %estimation:
-    if co.additive_constant_is_relevant %the additive constant is relevant in the Rényi entropy estimation
+    if co.additive_constant_is_relevant %the additive constant is relevant in the Renyi entropy estimation
         FN = filename_of_HRenyi_constant(d,co);
         if exist(FN)
             load(FN,'constant');

code/H_I_D/base_estimators/HRenyi_kNN_1tok_initialization.m

             %co.k = 3;%k-nearest neighbors
             %co.epsi = 0; %=0: exact kNN; >0: approximate kNN, the true (not squared) distances can not exceed the real distance more than a factor of (1+epsi).
 
-    co.alpha = 0.99; %The Renyi entropy equals to the Shannon differential entropy, in limit, i.e., Renyi=H_{R,alpha} -> Shannon=H, provided that alpha ->1.
+    co.alpha = 0.99; %The Renyi entropy (H_{R,alpha}) equals to the Shannon differential entropy (H) in limit: H_{R,alpha} -> H, as alpha -> 1.
     co.additive_constant_is_relevant = 0; %1:additive constant is relevant (you can precompute it via 'estimate_HRenyi_constant.m'), 0:not relevant
     
 %initialize the ann wrapper in Octave, if needed:

code/H_I_D/base_estimators/HRenyi_kNN_S_estimation.m

     L = compute_length_HRenyi_kNN_S(Y,co);
 
 %estimation:
-    if co.additive_constant_is_relevant %the additive constant is relevant in the Rényi entropy estimation
+    if co.additive_constant_is_relevant %the additive constant is relevant in the Renyi entropy estimation
         FN = filename_of_HRenyi_constant(d,co);
         if exist(FN)
             load(FN,'constant');

code/H_I_D/base_estimators/HRenyi_kNN_S_initialization.m

             %co.k = [1,2,4];%=S: nearest neighbor set
             %co.epsi = 0; %=0: exact kNN; >0: approximate kNN, the true (not squared) distances can not exceed the real distance more than a factor of (1+epsi).
 				
-    co.alpha = 0.99; %The Renyi entropy equals to the Shannon differential entropy, in limit, i.e., Renyi=H_{R,alpha} -> Shannon=H, provided that alpha ->1.
+    co.alpha = 0.99; %The Renyi entropy (H_{R,alpha) equals to the Shannon differential entropy (H) in limit: H_{R,alpha} -> H, as alpha -> 1.
     co.additive_constant_is_relevant = 0; %1:additive constant is relevant (you can precompute it via 'estimate_HRenyi_constant.m'), 0:not relevant    
     
 %initialize the ann wrapper in Octave, if needed:

code/H_I_D/base_estimators/HRenyi_kNN_k_initialization.m

             %co.k = 3;%k-nearest neighbors
             %co.epsi = 0; %=0: exact kNN; >0: approximate kNN, the true (not squared) distances can not exceed the real distance more than a factor of (1+epsi).
             
-    co.alpha = 0.95; %The Renyi entropy equals to the Shannon differential entropy, in limit, i.e., Renyi=H_{R,alpha} -> Shannon=H, provided that alpha ->1.
+    co.alpha = 0.95; %The Renyi entropy (H_{R,alpha}) equals to the Shannon differential entropy (H) in limit: H_{R,alpha} -> H, as alpha -> 1.
     
 %initialize the ann wrapper in Octave, if needed:
     initialize_Octave_ann_wrapper_if_needed(co.kNNmethod);

code/H_I_D/base_estimators/HRenyi_weightedkNN_initialization.m

     co.mult = mult;
     
 %other fields:
-    co.alpha = 0.95; %The Renyi entropy equals to the Shannon differential entropy, in limit, i.e., Renyi=H_{R,alpha} -> Shannon=H, provided that alpha ->1.
+    co.alpha = 0.95; %The Renyi entropy (H_{R,alpha}) equals to the Shannon differential entropy (H) in limit: H_{R,alpha} -> Shannon=H, as alpha -> 1.
     %Possibilities for 'co.kNNmethod' (see 'kNN_squared_distances.m'): 
         %I: 'knnFP1': fast pairwise distance computation and C++ partial sort; parameter: co.k.
         %II: 'knnFP2': fast pairwise distance computation; parameter: co.k. 												 		

code/H_I_D/base_estimators/HTsallis_kNN_k_initialization.m

             %co.k = 3;%k-nearest neighbors
             %co.epsi = 0; %=0: exact kNN; >0: approximate kNN, the true (not squared) distances can not exceed the real distance more than a factor of (1+epsi).
             
-    co.alpha = 0.95; 
+    co.alpha = 0.95; %The Tsallis entropy (H_{T,alpha}) equals to the Shannon differential (H) entropy in limit: H_{T,alpha} -> H, as alpha -> 1.
     
 %initialize the ann wrapper in Octave, if needed:
     initialize_Octave_ann_wrapper_if_needed(co.kNNmethod);

code/H_I_D/meta_estimators/DJdistance_initialization.m

     co.mult = mult;
     
 %other fields:
-    co.member_name = 'Renyi_kNN_k'; %you can change it to any Kullback-Leibler entropy estimator, the Renyi divergence converges to the Shannon's as \alpha->1.
+    co.member_name = 'Renyi_kNN_k'; %you can change it to any Kullback-Leibler entropy estimator; the Renyi divergence (D_{R,alpha}) converges to the Shannon's (D): D_{R,alpha} -> D, as alpha-> 1.
     co.member_co = D_initialization(co.member_name,mult);

code/H_I_D/meta_estimators/HRPensemble_initialization.m

 %other fields:
     co.group_size = 500;%assumption: co.group_size <= number of samples
 	co.dim_RP = 2; %assumption: co.dim_RP <= dimension of the samples
-    co.member_name = 'Shannon_kNN_k'; %you can change it to any entropy estimator.
+    co.member_name = 'Shannon_kNN_k'; %you can change it to any entropy estimator
     co.member_co = H_initialization(co.member_name,mult);
-    
+    

code/H_I_D/meta_estimators/HTsallis_HRenyi_estimation.m

 %co.mult:OK.
 
 H =  H_estimation(Y,co.member_co); %Renyi entropy
-H = (exp(H * (1-co.alpha)) - 1) / (1-co.alpha); %Renyi->Tsallis entropy
+H = (exp(H * (1-co.alpha)) - 1) / (1-co.alpha); %transform Renyi entropy to Tsallis entropy
 

code/H_I_D/meta_estimators/IShannon_HShannon_initialization.m

     
 %other fields:
     co.member_name = 'Shannon_kNN_k'; %you can change it to any Shannon entropy estimator capable of 'multiplicative 
-	%factor correct' H estimation (mult=1) below. Note: Renyi entropy (H_{R,alpha}) also gives in limit (alpha->1) the Shannon entropy (H).
+	%factor correct' H estimation (mult=1) below. Note: Renyi entropy (H_{R,alpha}) also gives in limit the Shannon entropy (H): H_{R,alpha} -> H, as alpha -> 1.
     co.member_co = H_initialization(co.member_name,1);%'1': since we use the relation '(*)' (multiplicative factors in entropy estimations are NOT allowed).
    
     

code/H_I_D/utilities/Edgeworth_t1_t2_t3.cpp

+#include "mex.h"
+#include "math.h"
+
+/* .cpp version of 'Edgeworth_t1_t2_t3.m' */
+
+/* Copyright (C) 2012 Zoltan Szabo ("http://nipg.inf.elte.hu/szzoli", "szzoli (at) cs (dot) elte (dot) hu")
+%
+%This file is part of the ITE (Information Theoretical Estimators) Matlab/Octave toolbox.
+%
+%ITE 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 software 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 ITE. If not, see <http://www.gnu.org/licenses/>. */
+
+void mexFunction(
+    int nlhs, mxArray *plhs[],
+    int nrhs, const mxArray *prhs[])
+{
+
+    double *Y;    /* pointer to the input matrix */    
+    mwSize d,T;   /* d:dimension, T:number of samples */
+    
+    Y = mxGetPr(prhs[0]); /* create a pointer to the input matrix */
+    d = mxGetM(prhs[0]);  /* number of rows */
+    T = mxGetN(prhs[0]);  /* number of columns */
+
+    double *t1,*t2,*t3;    
+    /* Allocating array for t1,t2,t3 */
+    plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
+    t1 = mxGetPr(plhs[0]);    
+    plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
+    t2 = mxGetPr(plhs[1]);    
+    plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL);
+    t3 = mxGetPr(plhs[2]);    
+    
+    mwSize i,j,k,u;
+    double kappa_iii,kappa_iij,kappa_ijk;
+    
+    *t1 = 0;
+    for (i=1; i<d+1; i++) {      /* i: in Matlab sense */
+        kappa_iii = 0;
+        for (u=1; u<T+1; u++) {  /* u: in Matlab sense */
+            kappa_iii = kappa_iii + pow(Y[(i-1)+(u-1)*d],3);
+        }
+        *t1 = *t1 + pow(kappa_iii/T,2);
+    }
+    
+    *t2 = 0;
+    /* i<j */
+    for (i=1; i<(d-1)+1; i++) {   /* i: in Matlab sense */
+    for (j=i+1; j<d+1; j++) {     /* j: in Matlab sense */
+         kappa_iij = 0;
+         for (u=1; u<T+1; u++) {  /* u: in Matlab sense */
+             kappa_iij = kappa_iij + pow(Y[(i-1)+(u-1)*d],2) * Y[(j-1)+(u-1)*d];
+         }
+         *t2 = *t2 + pow(kappa_iij/T,2);           
+    }        
+    }
+    
+    /* j<i */
+    for (j=1; j<(d-1)+1; j++) {  /* i: in Matlab sense */
+    for (i=j+1; i<d+1; i++) {    /* j: in Matlab sense */
+        kappa_iij = 0;
+        for (u=1; u<T+1; u++) {  /* u: in Matlab sense */
+            kappa_iij = kappa_iij + pow(Y[(i-1)+(u-1)*d],2) * Y[(j-1)+(u-1)*d];
+        }
+        *t2 = *t2 + pow(kappa_iij/T,2);           
+    }        
+    }
+    *t2 = 3 * *t2;
+    
+    *t3 = 0;
+    /* i<j<k */
+    for (i=1; i<(d-2)+1; i++) {   /* i: in Matlab sense */
+    for (j=i+1; j<(d-1)+1; j++) { /* j: in Matlab sense */
+    for (k=j+1; k<d+1; k++) {     /* k: in Matlab sense */
+        kappa_ijk = 0;
+        for (u=1; u<T+1; u++) {   /* u: in Matlab sense */
+            kappa_ijk = kappa_ijk + Y[(i-1)+(u-1)*d] * Y[(j-1)+(u-1)*d] *Y[(k-1)+(u-1)*d];
+        }
+        *t3 = *t3 + pow(kappa_ijk/T,2);
+    }
+    }
+    }
+    *t3 = *t3 / 6;
+}

code/ITE_install.m

 %You should have received a copy of the GNU General Public License along with ITE. If not, see <http://www.gnu.org/licenses/>.
 
 %parameters:
-    compile_NCut = 1;%1=compile, 0=do not compile
-    compile_ANN = 1;%1=compile, 0=do not compile
-    compile_knn = 1;%1=compile, 0=do not compile
-	compile_TCA = 1;%1=compile, 0=do not compile (chol_gauss); not necessary, but can further speed-up the computations; the package also contains the purely Matlab/Octave 'chol_gauss.m'
-    compile_SWICA = 1; %1=compile, 0=do not compile; not necessary, but can accelerate computations; the package also contains the purely Matlab/Octave 'SW_kappa.m' and 'SW_sigma.m'
-    compile_Hoeffding_term1 = 1; %1=compile, 0=do not compile; not necessary, but can be more ecomical in terms of memory used + accelerate computations; the package also contains the purely Matlab/Octave 'Hoeffding_term1.m'
+	%compile:
+		compile_NCut = 1;%1=compile, 0=do not compile
+		compile_ANN = 1;%1=compile, 0=do not compile
+		compile_knn = 1;%1=compile, 0=do not compile
+		compile_TCA = 1;%1=compile, 0=do not compile (chol_gauss); not necessary, but can further speed-up the computations; the package also contains the purely Matlab/Octave 'chol_gauss.m'
+		compile_SWICA = 1; %1=compile, 0=do not compile; not necessary, but can accelerate computations; the package also contains the purely Matlab/Octave 'SW_kappa.m' and 'SW_sigma.m'
+		compile_Hoeffding_term1 = 1; %1=compile, 0=do not compile; not necessary, but can be more ecomical in terms of memory used + accelerate computations; the package also contains the purely Matlab/Octave 'Hoeffding_term1.m'
+		compile_Edgeworth_t1_t2_t3 = 1; %1=compile, 0=do not compile; not necessary, but can accelerate computations; the package also contains the purely Matlab/Octave 'Edgeworth_t1_t2_t3.m'
+    %download:
+		download_ARfit = 1;%1=download+extract,0=do not download
 	delete_Ncut_in_Octave = 0; %1=yes, 0=no
 
 disp('Installation: started.');
         rmdir(strcat(ITE_code_dir,'/shared/embedded/ann_wrapperM/'),'s');%delete 'ann_wrapperM' with subfolders.    
         disp('We are working in Octave environment => ann_wrapper for Matlab: deleted.');%the 'ann' file may otherwise cause conflicts
     end
-    
-%download and extract the ARfit package to '/shared/embedded/ARfit':
+	
+if download_ARfit  %download and extract the ARfit package to '/shared/embedded/ARfit':
     disp('ARfit package: downloading, extraction: started.');
     [FN,status] = urlwrite('http://www.gps.caltech.edu/~tapio/arfit/arfit.zip','arfit.zip');
     if status %downloading: successful
     else
         disp('Error: ARfit package could not been downloaded.');
     end    
-  
+end
+
 %add 'ITE_code_dir' with its subfolders to the Matlab PATH:
     addpath(genpath(ITE_code_dir));
     if environment_Matlab
         mex Hoeffding_term1.cpp;
         disp('Hoeffding_term1.cpp compilation: ready.');
     end
+
+%compile 'Edgeworth_t1_t2_t3', if needed:
+    if compile_Edgeworth_t1_t2_t3
+        cd(strcat(ITE_code_dir,'/H_I_D/utilities'));
+        disp('Edgeworth_t1_t2_t3.cpp compilation: started.');
+        mex Edgeworth_t1_t2_t3.cpp;
+        disp('Edgeworth_t1_t2_t3.cpp compilation: ready.');
+	end
     
 %compile 'knn', if needed:
     if compile_knn%if needed
         disp('knn (top.cpp) compilation: started.');
         build; %compile the knn package (top.cpp)
         disp('knn (top.cpp) compilation: ready.');
-    end	
-    
+    end
+	
 %change back to the stored working directory:
     cd(pd);
     
 			end
     end        
     %ARfit:
-        d = 2; A = 0.9 * random_orthogonal(d); T = 1000;   Ls = 1;
-        v = arsim(zeros(d,1),A,eye(d),T);
-        [w_temp, Fx_hat, C, SBCs] = arfit(v, Ls, Ls);
-        disp('ARfit quick test: successful.');
+		if download_ARfit %verify only in case of downloading
+			d = 2; A = 0.9 * random_orthogonal(d); T = 1000;   Ls = 1;
+			v = arsim(zeros(d,1),A,eye(d),T);
+			[w_temp, Fx_hat, C, SBCs] = arfit(v, Ls, Ls);
+			disp('ARfit quick test: successful.');
+		end
 	%knn:
 		if compile_knn%verify only in case of compilation
 			Y = rand(3,100); Q = rand(3,200);
 			knn(Q,Y,1);
 			disp('knn quick test: successful.');
-		end
+		end
Add a comment to this file

doc/ITE_documentation.pdf

Binary file modified.

Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.