00001
00003 #include <iostream>
00004 #include <string>
00005 #include <vector>
00006 #include <fstream>
00007
00008 #include <cstdio>
00009 #include <cstdlib>
00010 #include <cmath>
00011 #include <cassert>
00012
00013
00014 #include "MI.hh"
00015 using namespace std;
00016
00017 MI::MI(){}
00018
00019 MI::~MI(){}
00020
00021 void MI::Allocate(int nb_dimensions){
00022
00023 nb_dim=nb_dimensions;
00024
00025 float* AuxMeanT1=new float[nb_dim];
00026 float** AuxCovT1=new float*[nb_dim];
00027 for (int i=0;i<nb_dim;i++){
00028 AuxCovT1[i]=new float[nb_dim];
00029 }
00030 float* AuxdeterminantT1=new float[1];
00031 MeanT1.push_back(AuxMeanT1);
00032 CovT1.push_back(AuxCovT1);
00033 determinantT1.push_back(AuxdeterminantT1[0]);
00034
00035
00036 float* AuxMeanT2=new float[nb_dim];
00037 float** AuxCovT2=new float*[nb_dim];
00038 for (int i=0;i<nb_dim;i++){
00039 AuxCovT2[i]=new float[nb_dim];
00040 }
00041 float* AuxdeterminantT2=new float[1];
00042 MeanT2.push_back(AuxMeanT2);
00043 CovT2.push_back(AuxCovT2);
00044 determinantT2.push_back(AuxdeterminantT2[0]);
00045
00046 float* AuxMeanT12=new float[nb_dim*2];
00047 float** AuxCovT12=new float*[nb_dim*2];
00048 for (int i=0;i<nb_dim*2;i++){
00049 AuxCovT12[i]=new float[nb_dim*2];
00050 }
00051 float* AuxdeterminantT12=new float[1];
00052 MeanT12.push_back(AuxMeanT12);
00053 CovT12.push_back(AuxCovT12);
00054 determinantT12.push_back(AuxdeterminantT12[0]);
00055
00056
00057 }
00058
00059 void MI::Initialise1(vector<float**> *C1,vector<float*> *M1)
00060 {
00061 CovT1=*C1;
00062 MeanT1=*M1;
00063 for (int j1=0;j1<nb_dim;j1++){CovT1[0][0][j1]=0;CovT1[0][j1][0]=0;}
00064 DetCovMatrixT1();
00065
00066 }
00067
00068 void MI::Initialise2(vector<float**> *C2,vector<float*> *M2)
00069 {
00070 CovT2=*C2;
00071 MeanT2=*M2;
00072
00073 DetCovMatrixT2();
00074
00075 }
00076
00077 void MI::Initialise12(vector<float**> *C12,vector<float*> *M12)
00078 {
00079 CovT12=*C12;
00080 MeanT12=*M12;
00081
00082 DetCovMatrixT12();
00083 }
00084
00085 void MI::DetCovMatrixT1(){
00086
00087 int j1,j2,j;
00088
00089
00090
00091
00092
00093 int nb_Gdim=nb_dim;
00094 float aux;
00095 for (j1=0;j1<nb_dim;j1++){
00096 if (CovT1[0][j1][j1]==0){
00097 nb_Gdim--;
00098 }
00099 }
00100 bool again=true;
00101 while(again){
00102 for (j1=0;j1<nb_dim;j1++)
00103 {
00104 if (CovT1[0][j1][j1]==0){
00105 if (j1<nb_dim-1){
00106 for (j2=0;j2<nb_dim;j2++){
00107 aux=CovT1[0][j1][j2];
00108 CovT1[0][j1][j2]=CovT1[0][j1+1][j2];
00109 CovT1[0][j1+1][j2]=aux;
00110
00111
00112 }
00113 for (j2=0;j2<nb_dim;j2++){
00114 aux=CovT1[0][j2][j1];
00115 CovT1[0][j2][j1]=CovT1[0][j2][j1+1];
00116 CovT1[0][j2][j1+1]=aux;
00117
00118 }
00119 }
00120 }
00121 }
00122 for (j1=0;j1<nb_Gdim;j1++){
00123 if (CovT1[0][j1][j1]==0){again=true;}
00124 else {again=false;}
00125 }
00126 if (nb_Gdim==0){again=false;}
00127 }
00128
00129
00130
00131
00132
00133
00134 float **auxCov=new float*[nb_Gdim+1];
00135 for (j=0;j<nb_Gdim+1;j++)
00136 auxCov[j]=new float[nb_Gdim+1];
00137 float d;
00138 int *indx=new int[nb_Gdim+1];
00139 int n=nb_Gdim;
00140 void ludcmp(float **a,int n,int *indx,float *d);
00141
00142 for (j1=0;j1<nb_Gdim;j1++)
00143 {
00144 for (j2=0;j2<nb_Gdim;j2++){auxCov[j1+1][j2+1]=CovT1[0][j1][j2];}
00145 }
00146
00147 ludcmp(auxCov,n,indx,&d);
00148 for (j=1;j<nb_Gdim+1;j++){d=auxCov[j][j]*d;}
00149 determinantT1[0]=d;
00150
00151
00152 for (j=0;j<nb_Gdim+1;j++){delete [] auxCov[j];}
00153 delete [] auxCov;delete [] indx;
00154
00155 }
00156
00157 void MI::DetCovMatrixT2(){
00158
00159 int j1,j2,j;
00160
00161
00162
00163
00164
00165
00166 int nb_Gdim=nb_dim;
00167 float aux;
00168 for (j1=0;j1<nb_dim;j1++){
00169 if (CovT2[0][j1][j1]==0){
00170 nb_Gdim--;
00171 }
00172 }
00173 bool again=true;
00174 while(again){
00175 for (j1=0;j1<nb_dim;j1++)
00176 {
00177 if (CovT2[0][j1][j1]==0){
00178 if (j1<nb_dim-1){
00179 for (j2=0;j2<nb_dim;j2++){
00180 aux=CovT2[0][j1][j2];
00181 CovT2[0][j1][j2]=CovT2[0][j1+1][j2];
00182 CovT2[0][j1+1][j2]=aux;
00183
00184
00185 }
00186 for (j2=0;j2<nb_dim;j2++){
00187 aux=CovT2[0][j2][j1];
00188 CovT2[0][j2][j1]=CovT2[0][j2][j1+1];
00189 CovT2[0][j2][j1+1]=aux;
00190
00191 }
00192 }
00193 }
00194 }
00195 for (j1=0;j1<nb_Gdim;j1++){
00196 if (CovT2[0][j1][j1]==0){again=true;}
00197 else {again=false;}
00198 }
00199 if (nb_Gdim==0){again=false;}
00200 }
00201
00202
00203
00204
00205
00206 float **auxCov=new float*[nb_Gdim+1];
00207 for (j=0;j<nb_Gdim+1;j++)
00208 {
00209 auxCov[j]=new float[nb_Gdim+1];
00210 }
00211 float d;
00212 int *indx=new int[nb_Gdim+1];
00213 int n=nb_Gdim;
00214 void ludcmp(float **a,int n,int *indx,float *d);
00215
00216 for (j1=0;j1<nb_Gdim;j1++)
00217 {
00218 for (j2=0;j2<nb_Gdim;j2++){auxCov[j1+1][j2+1]=CovT2[0][j1][j2];}
00219
00220 }
00221 ludcmp(auxCov,n,indx,&d);
00222 for (j=1;j<nb_Gdim+1;j++){d=auxCov[j][j]*d;}
00223 determinantT2[0]=d;
00224
00225
00226 for (j=0;j<nb_Gdim+1;j++){delete [] auxCov[j];}
00227 delete [] auxCov;delete [] indx;
00228
00229
00230 }
00231
00232 void MI::DetCovMatrixT12(){
00233
00234 int nb_dimAux=nb_dim*2;
00235 int j1,j2,j;
00236
00237
00238
00239
00240
00241 int nb_Gdim=nb_dimAux;
00242 float aux;
00243 for (j1=0;j1<nb_dimAux;j1++){
00244 if (CovT12[0][j1][j1]==0){
00245 nb_Gdim--;
00246 }
00247 }
00248 bool again=true;
00249 while(again){
00250 for (j1=0;j1<nb_dimAux;j1++)
00251 {
00252 if (CovT12[0][j1][j1]==0){
00253 if (j1<nb_dimAux-1){
00254 for (j2=0;j2<nb_dimAux;j2++){
00255 aux=CovT12[0][j1][j2];
00256 CovT12[0][j1][j2]=CovT12[0][j1+1][j2];
00257 CovT12[0][j1+1][j2]=aux;
00258
00259
00260 }
00261 for (j2=0;j2<nb_dimAux;j2++){
00262 aux=CovT12[0][j2][j1];
00263 CovT12[0][j2][j1]=CovT12[0][j2][j1+1];
00264 CovT12[0][j2][j1+1]=aux;
00265
00266 }
00267 }
00268 }
00269 }
00270 for (j1=0;j1<nb_Gdim;j1++){
00271 if (CovT12[0][j1][j1]==0){again=true;}
00272 else {again=false;}
00273 }
00274 if (nb_Gdim==0){again=false;}
00275 }
00276
00277
00278
00279
00280
00281
00282 float **auxCov=new float*[nb_Gdim+1];
00283 for (j=0;j<nb_Gdim+1;j++)
00284 {
00285 auxCov[j]=new float[nb_Gdim+1];
00286 }
00287 float d;
00288 int *indx=new int[nb_Gdim+1];
00289 int n=nb_Gdim;
00290 void ludcmp(float **a,int n,int *indx,float *d);
00291
00292 for (j1=0;j1<nb_Gdim;j1++)
00293 {
00294 for (j2=0;j2<nb_Gdim;j2++){auxCov[j1+1][j2+1]=CovT12[0][j1][j2];}
00295
00296 }
00297 ludcmp(auxCov,n,indx,&d);
00298 for (j=1;j<nb_Gdim+1;j++){d=auxCov[j][j]*d;}
00299 determinantT12[0]=d;
00300
00301 for (j=0;j<nb_Gdim+1;j++){delete [] auxCov[j];}
00302 delete [] auxCov;delete [] indx;
00303
00304
00305
00306 }
00307
00308
00309
00310 float MI::MICalculus()
00311 {
00312 float res=0.5*log((determinantT1[0]*determinantT2[0])/determinantT12[0])/log(2);
00313 return res;
00314 }
00315