aimsalgo  5.1.2
Neuroimaging image processing
mixtureOfPpca_d.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 #ifndef MIXTURE_OF_PPCA_D_H
36 #define MIXTURE_OF_PPCA_D_H
37 
38 #include <cstdlib>
39 #include <aims/classification/mixtureOfPpca.h>
40 #include <aims/io/finder.h>
41 #include <aims/io/writer.h>
42 #include <aims/vector/vector.h>
43 #include <aims/math/gausslu.h>
44 #include <aims/bucket/bucketMap.h>
45 #include <cartobase/type/string_conversion.h>
46 #include <aims/math/svd.h>
47 #include <vector>
48 #include <iostream>
49 #include <iomanip>
50 #include <math.h>
51 #include <aims/math/random.h>
52 
53 /* This code seems to be used... nowhere...
54  Anyway it is not currently working: namespace problems will prevent
55  compiling.
56 */
57 
58 namespace aims
59 {
60 
61 double * PpcaAnalyserElement::_exp = 0 ;
62 
63 PpcaAnalyserElement::PpcaAnalyserElement( int significantNumberOfVp, bool useOnlyCorrIndiv ) :
64  _useOnlyCorrIndiv( useOnlyCorrIndiv ), _significantNumberOfVp( significantNumberOfVp ), _computed( false ), _energy(0.), _valid(true)
65 {
66  if( !_exp ){
67  _exp = new double[75010] ;
68  for( int i = 0 ; i < 7510 ; ++i )
69  _exp[i] = exp( -i/100. ) ;
70  }
71 }
72 
73 
74 template <class T>
75 void
76 PpcaAnalyserElement::init(
77  const std::list<int>& selectedIndividuals, double initialPi,
78  const carto::rc_ptr<carto::Volume<T> >& individuals, double noiseRef )
79 {
80  carto::VolumeRef<T> indivMatrix( max(int(selectedIndividuals.size()), 1), individuals->getSizeY() ) ;
81  list<int>::const_iterator iter( selectedIndividuals.begin() ),
82  last( selectedIndividuals.end() ) ;
83  int i = 0 ;
84 
85  while( iter != last ){
86  for( int t = 0 ; t < individuals->getSizeY() ; ++t )
87  indivMatrix( i, t ) = individuals( *iter, t ) ;
88  ++i ;
89  ++iter ;
90  }
91 
92  _Pi = initialPi ;
93  doIt( indivMatrix, individuals->getSizeX(), noiseRef ) ;
94 }
95 
96 
97 template <class T>
98 void
99 PpcaAnalyserElement::doIt( const carto::rc_ptr<carto::Volume<T> >& indivMatrix,
100  int totalNbOfIndividuals, double noiseRef )
101 {
102  _computed = true ;
103  int nbInd = indivMatrix->getSizeX() ;
104  int nbFrame = indivMatrix->getSizeY() ;
105 
106  carto::VolumeRef<double> centeredIndivMatrix( nbInd, nbFrame ) ;
107  _mean = carto::VolumeRef<double>( nbFrame ) ;
108 
109  for( int ind = 0 ; ind < nbInd ; ++ind )
110  for( int t = 0 ; t < nbFrame ; ++t )
111  _mean[t] += indivMatrix( ind, t ) ;
112  for( int t = 0 ; t < nbFrame ; ++t ){
113  _mean[t] /= centeredIndivMatrix->getSizeX() ;
114  for( int ind = 0 ; ind < centeredIndivMatrix->getSizeX() ; ++ind )
115  centeredIndivMatrix( ind, t ) = indivMatrix( ind, t ) - _mean[t] ;
116  }
117 // cout << "mean = " << _mean << endl ;
118 
119  // Matrice des correlations
120  carto::VolumeRef<double> matVarCov( nbFrame, nbFrame ) ;
121  int x, y ;
122  ForEach2d( matVarCov, x, y ){
123  for( int k = 0 ; k < centeredIndivMatrix->getSizeX() ; ++k )
124  matVarCov(x, y) += centeredIndivMatrix(k, x) * centeredIndivMatrix(k, y) ;
125  matVarCov(x, y) /= centeredIndivMatrix->getSizeX() - 1 ;
126  }
127 
128  // Decomposition SVD
129  AimsSVD<double> svd;
130  svd.setReturnType( AimsSVD<double>::VectorOfSingularValues ) ;
131  carto::VolumeRef<double> eigenValues = svd.doit( matVarCov ) ;
132  svd.sort( matVarCov, eigenValues ) ;
133 // cout << "valeurs propres = " << eigenValues << endl ;
134 
135  // Calcul de "noise variance"(somme des valeurs propres non significatives/par leur nb)=sigma^2
136  _sigma2 = 0. ;
137  for( int i = _significantNumberOfVp ; i < nbFrame ; ++i )
138  _sigma2 += eigenValues[i] ;
139  _sigma2 /= (double)( nbFrame - _significantNumberOfVp ) ;
140  cout << endl << "sigma2 = " << setprecision(3) << _sigma2 << endl ;
141 
142  carto::VolumeRef<double> selectedEigenVectors( nbFrame, _significantNumberOfVp ) ; //dim dxq
143  for( int i = 0 ; i < _significantNumberOfVp ; ++i ){
144  for( int t = 0 ; t < nbFrame ; ++t ){
145  selectedEigenVectors(t, i) = matVarCov( t, i ) ;
146  }
147  }
148 // cout << "selected EV = " << endl << setprecision(3) << selectedEigenVectors << endl ;
149 
150  // Calcul de la "weight" matrice Wi (dim. dxq)
151  carto::VolumeRef<double> tempWi( _significantNumberOfVp, _significantNumberOfVp ) ;
152 
153  for( int i = 0 ; i < _significantNumberOfVp ; ++i )
154  tempWi(i, i) = sqrt( eigenValues[i] - _sigma2 ) ;
155  _Wi = selectedEigenVectors.cross( tempWi ) ;
156 
157  // calcul de la matrice de covariance Ci = sigma^2I+WWT
158  carto::VolumeRef<double> Ci = _Wi.cross( _Wi.clone().transpose() ) ;
159  for( int i = 0 ; i < nbFrame ; ++i )
160  Ci(i, i) += _sigma2 ;
161 
162  // calcul de detCi et invCi
163  double s = 0.000001 ;
164  carto::VolumeRef< double > u = Ci.clone(); // because it is replaced on SVD output
165  carto::VolumeRef< double > v( Ci->getSizeY(), Ci->getSizeY() );
166  AimsSVD< double > svd2;
167  svd2.setReturnType( AimsSVD<double>::MatrixOfSingularValues ) ;
168  carto::VolumeRef< double > w = svd2.doit( u, &v );
169 
170  // Set the singular values lower than threshold to zero
171  double ts = s * w.maximum();
172  int i, n = w->getSizeX();
173  double detCi = 1. ;
174  for ( i = 0 ; i < n; i++ ){
175  if ( w( i, i ) < 0. ){
176  cout << endl << "pour i = " << i << ", valeur w = " << w( i, i ) << " ! VALEUR NEGATIVE CHANGEE !" << endl ;
177  w( i, i ) = - w( i, i ) ;
178  }
179  detCi *= w( i, i ) / _sigma2 ;
180  if ( w( i, i ) > ts ) w( i, i ) = 1.0f / w( i, i );
181  else w( i, i ) = 0.0f;
182  }
183  // Get the SVD inverse of h ( h-1 = V.W-2.V'.H' = V.W-1.U' )
184  // for solution of H'.H.x = H'.b <=> min(||H.x-b||^2)
185  _invCi = v.cross( w.cross( u.transpose() ) ) ;
186 // cout << "invCi = " << endl << _invCi << endl ;
187 
188 // cout << "det Ci = " << detCi << endl ;
189  _normFactor = 1 / ( pow( detCi, 0.5 ) * pow( _sigma2 / noiseRef, nbFrame / 2. ) ) ;
190 // cout << "normFactor avant = " << 1 / pow( detCi, 0.5 ) << endl ;
191  cout << "norm factor avec noiseRef = " << _normFactor << endl ;
192 
193  // calcul de la matrice de covariance posterieure Mi = sigma^2I+WTW et de invMi
194  carto::VolumeRef<double> Mi = ( _Wi.clone().transpose() ).cross( _Wi ) ;
195  for( int i = 0 ; i < _significantNumberOfVp ; ++i )
196  Mi(i, i) += _sigma2 ;
197 
198  u = Mi.clone(); // because it is replaced on SVD output
199  v( Mi->getSizeY(), Mi->getSizeY() );
200  AimsSVD< double > svd3;
201  w = svd3.doit( u, &v );
202 
203  ts = s * w.maximum();
204  n = w->getSizeX();
205  for ( i = 0 ; i < n; i++ ){
206  if ( w( i, i ) > ts ) w( i, i ) = 1.0f / w( i, i );
207  else w( i, i ) = 0.0f;
208  }
209  _invMi = v.cross( w.cross( u.transpose() ) ) ;
210 // cout << "invMi = " << _invMi << endl ;
211 
212  _Rn = carto::VolumeRef<double>( totalNbOfIndividuals ) ;
213 }
214 
215 
216 template <class T>
217 bool
218 PpcaAnalyserElement::newStep1(
219  const carto::rc_ptr<carto::Volume<T> >& indivMatrix, bool useOnlyCorrIndiv )
220 {
221  int nbInd = indivMatrix->getSizeX() ;
222  int nbFrame = indivMatrix->getSizeY() ;
223 
225 // calcul pour tous les individus de p(tn|i)
226 // + CALCUL DE LA CORRELATION ENTRE CHAQUE INDIV ET LE VECTEUR MOYEN
227 
228  carto::VolumeRef<double> indivn( nbFrame ) ;
229  carto::VolumeRef<double> centeredIndivn( nbFrame ) ;
230  _dist = carto::VolumeRef<double>( nbInd ) ;
231  _pTni = carto::VolumeRef<double>( nbInd ) ;
232 
233  double sumInd, sumMean, indMean, meanMean, corr = 0. ;
234  int count = 0 ;
235  carto::VolumeRef<double> antiCorr( nbInd, nbFrame ) ;
236  carto::VolumeRef<double> Corr( nbInd, nbFrame ) ;
237  int nbIndTaken = 0 ;
238 
239  //float invDimNormFactor = 1. / ( nbFrame - _significantNumberOfVp ) ;
240 
241  for( int n = 0 ; n < nbInd ; ++n ){
242  sumInd = 0. ;
243  sumMean = 0. ;
244  indMean = 0. ;
245  meanMean = 0. ;
246  corr = 0. ;
247 
248  for( int t = 0 ; t < nbFrame ; ++t ){
249  indivn[t] = indivMatrix(n, t) ;
250  centeredIndivn[t] = indivMatrix(n, t) - _mean[t] ;
251  sumInd += indivn[t] ;
252  sumMean += _mean[t] ;
253  }
254  // AJOUT 0411 calcul de la corr�lation entre chaque individu et le vecteur moyen
255  if( useOnlyCorrIndiv == 1 ){
256  indMean = sumInd / (double)nbFrame ;
257  meanMean = sumMean / (double)nbFrame ;
258  for( int t = 0 ; t < nbFrame ; ++t )
259  corr += ( indivn[t] - indMean ) * ( _mean[t] - meanMean ) ;
260 
261  if( corr < 0. ){
262  // probabilit� nulle si anti-corr�l�
263 /* for( int t = 0 ; t < nbFrame ; ++t ) */
264 /* antiCorr( nbAntiCorr, t ) = indivn[t] ; */
265 /* ++nbAntiCorr ; */
266  _pTni[n] = 0. ;
267  ++count ;
268  } else {
269  // calcul de p(tn|i)
270 /* for( int t = 0 ; t < nbFrame ; ++t ) */
271 /* Corr( nbCorr, t ) = indivn[t] ; */
272 /* ++nbCorr ; */
273  _dist[n] = ( centeredIndivn.clone().transpose() ).cross( _invCi.cross( centeredIndivn ) )(0, 0) ;
274  if( _dist[n] > 1500 )
275  _pTni[n] = 0 ;
276  else
277  _pTni[n] = /*pow( 2*M_PI, - nbFrame/2. ) * */ _normFactor * _exp[ int(0.5 * _dist[n] * 100. + 0.000001 ) ] ;
278  }
279 
280  } else {
281  _dist[n] = ( centeredIndivn.clone().transpose() ).cross( _invCi.cross( centeredIndivn ) )(0, 0) ;
282  if( _dist[n] > 1500 )
283  _pTni[n] = 0 ;
284  else
285  _pTni[n] = /*pow( 2*M_PI, - nbFrame/2. ) * */ _normFactor * _exp[ int(0.5 * _dist[n] * 100. + 0.000001) ] ;
286  }
287  if( _pTni[n] > 0. )
288  ++nbIndTaken ;
289  }
290 
291 /* Writer < carto::VolumeRef<double> > wr( "correlated.ima" ) ; */
292 /* wr.write( Corr ) ; */
293 /* Writer < carto::VolumeRef<double> > wr2( "anti-correlated.ima" ) ; */
294 /* wr2.write( antiCorr ) ; */
295 /* Writer < carto::VolumeRef<double> > wr3( "mean.ima" ) ; */
296 /* wr3.write( _mean ) ; */
297 
298 /* double sumPtni = 0., prodPtni = 0., meanPtni = 0. , varPtni = 0. ; */
299 /* // moyenne et variance des pTni */
300 /* for( int n = 0 ; n < nbInd ; ++n ){ */
301 /* sumPtni += _pTni[n] ; */
302 /* prodPtni += _pTni[n] * _pTni[n] ; */
303 /* } */
304 /* meanPtni = sumPtni / nbInd ; */
305 /* varPtni = prodPtni / nbInd ; */
306 /* cout << ( _useOnlyCorrIndiv ? " anticorreles = " : "" ) << ( _useOnlyCorrIndiv ? toString(count) : "" ) */
307 /* << " - moyenne et variance des p(tn|i) = " */
308 /* << meanPtni << " - " << varPtni << endl ; */
309  if( nbIndTaken < indivMatrix->getSizeY() * 2 )
310  _valid = false ;
311  else
312  _valid = true ;
313  return _valid ;
314 }
315 
316 
317 template <class T>
318 double
319 PpcaAnalyserElement::newStep2(
320  carto::rc_ptr<carto::Volume<double> > pTn,
321  const carto::rc_ptr<carto::Volume<T> >& indivMatrix,
322  double noiseRef )
323 {
324  bool explicitComputing = true ;
325 
326  int nbInd = indivMatrix->getSizeX() ;
327  int nbFrame = indivMatrix->getSizeY() ;
328  _energy = 0. ;
329 
331  // _Rn[i] = probabilityOfTnKnownI * Pi[i] / pTn ;
332 
333 
334 // _Rn = carto::VolumeRef<double>( nbInd ) ;
335  double previousRn ;
336  int nbIndTemp = 0 ; // nb d'individus dont p(tn)!=0
337 // double previousRn = 0. ;
338  _sumDiff2Rni = 0. ;
339 
340  //cout << "Rn computing begin ------ " ;
341  for( int n = 0 ; n < nbInd ; ++n ){
342  previousRn = _Rn[n] ;
343  if( pTn[n] != 0. ){ // on ignore les points problematiques
344  _Rn[n] = _pTni[n] * _Pi / pTn[n] ;
345  //_energy += 1 - 4. * (_Rn[n] - 0.5) * (_Rn[n] - 0.5) ;
346  ++nbIndTemp ;
347  }
348  else
349  _Rn[n] = 0. ;
350 
351  // if( n < 10 && _Rn[n] != 0.) cout << "n = " << n << " - nouveau/ancien Rn = " << _Rn[n] << " " << previousRn << endl ;
352  //_sumDiff2Rni += pow( _Rn[n] - previousRn, 2. ) ;
353  }
354  //_energy /= nbInd ;
355  //_sumDiff2Rni /= nbInd ;
356 
357  // cout << "Rn computing end" << endl ;
358 
359  cout << /*endl <<*/ "critere pour la classe = " << setprecision(9) << _sumDiff2Rni << " - " ;
360 
361 
363  // _Pi = 1/nbInd * sum(Rni)
364 
365  _Pi = 0. ;
366  for( int n = 0 ; n < nbInd ; ++n )
367  _Pi += _Rn[n] ;
368  _Pi /= nbIndTemp ;
369  cout << "Pi = " << setprecision(3) << _Pi << endl ;
370 
371  double nbIndClass = _Pi * nbIndTemp ;
372  cout << "'nb d'individus appartenant' a la classe = " << nbIndClass << endl ;
373 
374 
376  // _mean = sum(Rni*tn)/sum(Rni)
377 
378  double sumRn = 0. ;
379  carto::VolumeRef<double> sumRnTn( nbFrame ) ;
380  int nullRniNb = 0, nearNullRni = 0 ;
381  for( int n = 0 ; n < nbInd ; ++n ){
382  sumRn += _Rn[n] ;
383  if( _Rn[n] == 0. )
384  ++nullRniNb ;
385  if( _Rn[n] < 10e-30 )
386  ++nearNullRni ;
387  for( int t = 0 ; t < nbFrame ; ++t )
388  sumRnTn[t] += _Rn[n] * indivMatrix(n, t) ;
389  }
390 
391  cout << "Nombre de RNi nulls : " << nullRniNb << "et de Rni presque nulls" << nearNullRni << endl ;
392  for( int t = 0 ; t < nbFrame ; ++t )
393  _mean[t] = sumRnTn[t] / sumRn ;
394 // cout << "mean = " << _mean << endl ;
395 
396 
398  // Si = 1/Pi*nbInd sum( Rni (tn-mui) (tn-mui)T )
399 
400  carto::VolumeRef<double> centeredIndiv( nbFrame ) ;
401  carto::VolumeRef<double> Si( nbFrame, nbFrame ) ;
402 
403 /* for( int n = 0 ; n < nbInd ; ++n ) */
404 /* for( int t = 0 ; t < nbFrame ; ++t ) */
405 /* centeredIndivMatrix(n, t) = indivMatrix(n, t) - _mean[t] ; */
406 
407  int x, y ;
408  for( int n = 0 ; n < nbInd ; ++n )
409  if( _Rn[n] != 0. ){
410  for( int t = 0 ; t < nbFrame ; ++t )
411  centeredIndiv(t) = indivMatrix(n, t) - _mean[t] ;
412  ForEach2d( Si, x, y )
413  Si(x, y) += _Rn[n] * centeredIndiv(x) * centeredIndiv(y) ;
414  }
415  ForEach2d( Si, x, y )
416  Si(x, y) /= _Pi * nbIndTemp ;
417 
418 
419  //------------------- CALCUL EXPLICITE ---------------------------------
420 
421  double detCi = 1. ;
422 
423  if( explicitComputing ){
424  // Decomposition SVD
425 
426  //cout << "Compute Si svd begin" << endl ;
427  AimsSVD<double> svd;
428  svd.setReturnType( AimsSVD<double>::VectorOfSingularValues ) ;
429  carto::VolumeRef<double> eigenValues = svd.doit( Si ) ;
430  svd.sort( Si, eigenValues ) ;
431 
432  //cout << "Compute Si svd end" << endl ;
433 
434  // Calcul de "noise variance"(somme des valeurs propres non significatives/par leur nb)=sigma^2
435  _sigma2 = 0. ;
436  for( int i = _significantNumberOfVp ; i < nbFrame ; ++i )
437  _sigma2 += eigenValues[i] ;
438  _sigma2 /= (double)( nbFrame - _significantNumberOfVp ) ;
439  cout << "sigma2 = " << _sigma2 << " (explicite)" << " - " ;
440 
441  carto::VolumeRef<double> selectedEigenVectors( nbFrame, _significantNumberOfVp ) ; //dim dxq
442  for( int i = 0 ; i < _significantNumberOfVp ; ++i )
443  for( int t = 0 ; t < nbFrame ; ++t )
444  selectedEigenVectors(t, i) = Si( t, i ) ;
445 
446  // Calcul de la "weight" matrice Wi (dim. dxq)
447  carto::VolumeRef<double> tempWi( _significantNumberOfVp, _significantNumberOfVp ) ;
448 
449  for( int i = 0 ; i < _significantNumberOfVp ; ++i )
450  tempWi(i, i) = sqrt( eigenValues[i] - _sigma2 ) ;
451  _Wi = selectedEigenVectors.cross( tempWi ) ;
452 
453  // calcul de la matrice de covariance Ci = sigma^2I+WWT
454  carto::VolumeRef<double> Ci = _Wi.cross( _Wi.clone().transpose() ) ;
455  for( int i = 0 ; i < nbFrame ; ++i )
456  Ci(i, i) += _sigma2 ;
457 
458  // calcul de detCi et invCi
459  double s = 0.000001 ;
460  carto::VolumeRef< double > u = Ci.clone(); // because it is replaced on SVD output
461  carto::VolumeRef< double > v( Ci->getSizeY(), Ci->getSizeY() );
462  AimsSVD< double > svd2;
463  svd2.setReturnType( AimsSVD<double>::MatrixOfSingularValues ) ;
464  carto::VolumeRef< double > w = svd2.doit( u, &v );
465 
466  // Set the singular values lower than threshold to zero
467  double ts = s * w.maximum();
468  int i, n = w->getSizeX();
469  for ( i = 0 ; i < n; i++ ){
470  if ( w( i, i ) < 0. ){
471  cout << endl << "pour i = " << i << ",valeur w = " << w( i, i ) << "! VALEUR NEGATIVE CHANGEE!" << endl ;
472  w( i, i ) = - w( i, i ) ;
473  }
474  detCi *= w( i, i ) / _sigma2 ;
475  if ( w( i, i ) > ts ) w( i, i ) = 1.0f / w( i, i );
476  else w( i, i ) = 0.0f;
477  }
478  // Get the SVD inverse of h ( h-1 = V.W-2.V'.H' = V.W-1.U' )
479  // for solution of H'.H.x = H'.b <=> min(||H.x-b||^2)
480  _invCi = v.cross( w.cross( u.transpose() ) ) ;
481 
482  } else {
483 
484  //---------------------Iterative computing----------------------------------
485 
487  // _Wi = SiWi (sigma2I + Mi-1 WiT Si Wi)-1 avec Mi = sigma2I + WiT*Wi
488 
489  carto::VolumeRef<double> newWi = _invMi.cross( ( _Wi.clone().transpose() ).cross( Si.cross( _Wi ) ) ) ;
490  for( int i = 0 ; i < newWi->getSizeX() ; ++i )
491  newWi(i, i) += _sigma2 ;
492 
493  double s = 0.000001 ;
494  carto::VolumeRef<double> u = newWi.clone(); // because it is replaced on SVD output
495  carto::VolumeRef<double> v( newWi->getSizeY(), newWi->getSizeY() );
496  AimsSVD<double> svd1;
497  svd1.setReturnType( AimsSVD<double>::MatrixOfSingularValues ) ;
498  carto::VolumeRef<double> w = svd1.doit( u, &v );
499  double ts = s * w.maximum();
500  int i, n = w->getSizeX();
501 
502  for ( i = 0 ; i < n ; i++ ){
503  if ( w( i, i ) < 0. ){
504  cout << endl << "pour i = " << i << ", w = " << w( i, i ) << " ! VALEUR NEGATIVE CHANGEE !" << endl ;
505  w( i, i ) = - w( i, i ) ;
506  }
507  if ( w( i, i ) > ts ) w( i, i ) = 1.0f / w( i, i );
508  else w( i, i ) = 0.0f;
509  }
510  newWi = v.cross( w.cross( u.transpose() ) ) ;
511 
512  newWi = Si.cross( _Wi.cross( newWi ) ) ;
513 // _Wi = newWi ;
514 
515 
517  // Mi = _sigmai2I + WiT Wi & invMi
518 
519 /* carto::VolumeRef<double> Mi = ( newWi.clone().transpose() ).cross( newWi ) ; */
520 /* for( int i = 0 ; i < Mi->getSizeX() ; ++i ) */
521 /* Mi(i, i) += _sigma2 ; */
522 
523 /* carto::VolumeRef<double> u2 = Mi.clone(); */
524 /* carto::VolumeRef<double> v2( Mi->getSizeY(), Mi->getSizeY() ) ; */
525 /* AimsSVD< double > svd2 ; */
526 /* svd2.setReturnType( AimsSVD<double>::MatrixOfSingularValues ) ; */
527 /* carto::VolumeRef<double> w2 = svd2.doit( u2, &v2 ) ; */
528 /* ts = s * w2.maximum(); */
529 /* n = w2->getSizeX(); */
530 
531 /* for( i = 0 ; i < n ; i++ ){ */
532 /* if ( w2( i, i ) < 0. ){ */
533 /* cout << endl << "pour i = " << i << ", w2 = " << w2( i, i ) << " ! VALEUR NEGATIVE CHANGEE !" << endl ; */
534 /* w2( i, i ) = - w2( i, i ) ; */
535 /* } */
536 /* if ( w2( i, i ) > ts ) w2( i, i ) = 1.0f / w2( i, i ) ; */
537 /* else w2( i, i ) = 0.0f; */
538 /* } */
539 /* _invMi = v2.cross( w2.cross( u2.transpose() ) ) ; */
540 
541 
543  // _sigmai2 = 1/d tr(Si - Si Wi Mi-1 WiT)
544 
545  _sigma2 = 0. ;
546  carto::VolumeRef<double> temp = Si.cross( _Wi.cross( _invMi.cross( newWi.clone().transpose() ) ) ) ;
547 
548  for( int i = 0 ; i < temp->getSizeX() ; ++i )
549  _sigma2 += ( Si(i, i) - temp(i, i) ) ;
550  _sigma2 /= nbFrame ;
551  cout << "sigma2 = " << _sigma2 << " - " ;
552  if( _sigma2 < 0. ){
553  cout << "VALEUR NEGATIVE CHANGEE ! " << endl ;
554  _sigma2 = - _sigma2 ;
555  }
556 
557  // essai
558  _Wi = newWi ;
559 
560 
562  // Mi = _sigmai2I + WiT Wi & invMi
563 
564  carto::VolumeRef<double> Mi = ( _Wi.clone().transpose() ).cross( _Wi ) ;
565  for( int i = 0 ; i < Mi->getSizeX() ; ++i )
566  Mi(i, i) += _sigma2 ;
567 
568  carto::VolumeRef<double> u2 = Mi.clone();
569  carto::VolumeRef<double> v2( Mi->getSizeY(), Mi->getSizeY() ) ;
570  AimsSVD< double > svd2 ;
571  svd2.setReturnType( AimsSVD<double>::MatrixOfSingularValues ) ;
572  carto::VolumeRef<double> w2 = svd2.doit( u2, &v2 ) ;
573  ts = s * w2.maximum();
574  n = w2->getSizeX();
575  for( i = 0 ; i < n ; i++ ){
576  if ( w2( i, i ) < 0. ){
577  cout << endl << "pour i = " << i << ", w2 = " << w2( i, i ) << " ! VALEUR NEGATIVE CHANGEE !" << endl ;
578  w2( i, i ) = - w2( i, i ) ;
579  }
580  if ( w2( i, i ) > ts ) w2( i, i ) = 1.0f / w2( i, i ) ;
581  else w2( i, i ) = 0.0f;
582  }
583  _invMi = v2.cross( w2.cross( u2.transpose() ) ) ;
584 
585 
587  // Ci = sigmai2 + Wi WiT --> Ci-1 & det(Ci)
588 
589  carto::VolumeRef<double> Ci = _Wi.cross( _Wi.clone().transpose() ) ;
590  for( int i = 0 ; i < nbFrame ; ++i )
591  Ci(i, i) += _sigma2 ;
592 
593  // verification Ci pas de valeurs negatives & symetrique
594  bool sym = true ;
595  for( i = 0 ; i < Ci->getSizeX() ; i++ ){
596  if( Ci(i, i) < 0 ){
597  cout << "pour i = " << i << ", Ci = " << Ci( i, i ) << " ! VALEUR NEGATIVE CHANGEE !" << endl ;
598  Ci( i, i ) = - Ci( i, i ) ;
599  }
600  for( int j = 0 ; j < Ci->getSizeY() ; j++ ){
601  if( Ci(i,j) != Ci(j,i) ){
602  sym = false ;
603  cout << "pour i,j = " << i << "," << j << " Ci(i,j) = " << Ci(i,j) << ",Ci(j,i) = " << Ci(j,i) << endl ;
604  }
605  }
606  }
607  if( sym == false ) cout << "--> Ci pas symetrique !" << endl ;
608 
609 
610  carto::VolumeRef<double> u3 = Ci.clone();
611  carto::VolumeRef<double> v3( Ci->getSizeY(), Ci->getSizeY() );
612  AimsSVD< double > svd3 ;
613  svd3.setReturnType( AimsSVD<double>::MatrixOfSingularValues ) ;
614  carto::VolumeRef<double> w3 = svd3.doit( u3, &v3 );
615  ts = s * w3.maximum();
616  n = w3->getSizeX();
617  bool saveCi = false ;
618 
619 // cout << endl << "valeurs diagonale w3 = " ;
620  for( i = 0 ; i < n ; i++ ){
621 // cout << w3( i, i ) << " " ;
622  if ( w3( i, i ) < 0. ){
623  saveCi = true ;
624  cout << endl << "pour i = " << i << ", valeur w3 = " << w3( i, i ) << " ! VALEUR NEGATIVE CHANGEE !" << endl ;
625  w3( i, i ) = - w3( i, i ) ;
626  }
627  if( _sigma2 == 0. )
628  detCi *= w3( i, i ) ;
629  else
630  detCi *= w3( i, i ) / _sigma2 ;
631  if ( w3( i, i ) > ts ) w3( i, i ) = 1.0f / w3( i, i );
632  else w3( i, i ) = 0.0f;
633  }
634  _invCi = v3.cross( w3.cross( u3.transpose() ) ) ;
635 // cout << "_invCi = " << endl << _invCi << endl ;
636 
637  // sauvegarde de Ci dans un fichier
638 /* if( saveCi == true ){ */
639 /* ofstream ci_file, diff_file ; */
640 /* ci_file.open( "Ci_file.txt" ) ; */
641 /* w3_file.open( "w3_file.txt" ) ; */
642 /* diff_file.open( "diff_file.txt" ) ; */
643 
644 /* bool newLine ; */
645 /* cout << "dimensions de Ci = " << Ci->getSizeX() << " " << Ci->getSizeY() << endl ; */
646 /* for( int j = 0 ; j < Ci->getSizeY() ; j++ ){ */
647 /* newLine = true ; */
648 /* for( int i = 0 ; i < Ci->getSizeX() ; i++ ){ */
649 /* if( newLine == true && j != 0 ) ci_file << endl ; */
650 /* ci_file << Ci( i, j ) << " " ; */
651 /* newLine = false ; */
652 /* } */
653 /* } */
654 /* ci_file << endl ; */
655 /* ci_file.close() ; */
656 
657 /* // verification */
658 /* carto::VolumeRef<double> prod = u3.cross( w3.cross( v3.clone().transpose() ) ) ; */
659 /* for( int j = 0 ; j < Ci->getSizeY() ; j++ ){ */
660 /* newLine = true ; */
661 /* for( int i = 0 ; i < Ci->getSizeX() ; i++ ){ */
662 /* if( newLine == true && j != 0 ) diff_file << endl ; */
663 /* diff_file << Ci(i,j) - prod(i,j) << " " ; */
664 /* newLine = false ; */
665 /* } */
666 /* } */
667 /* diff_file << endl ; */
668 /* diff_file.close() ; */
669 /* } */
670 
671  }
672 
673  // AJOUT 12/11: calcul de l'�nergie pour la mixture i Ei = nbOfInd * sigma2
674  //_energy = nbIndTemp * _sigma2 ;
675  cout << "energy = " << _energy ;
676 
677  if( _sigma2 == 0. )
678  _normFactor = 1 / pow( detCi, 0.5 ) ;
679  else
680  _normFactor = 1 / ( pow( detCi, 0.5 ) * pow( _sigma2 / noiseRef, nbFrame / 2. ) ) ;
681  cout << " - detCi = " << detCi /*<< " - normFactor sans noiseRef = " << 1 / pow( detCi, 0.5 )*/
682  << " - normFactor = " << _normFactor << endl ;
683 
684  return _sigma2 ;
685 }
686 
687 template <class T>
688 double
689 MixtureOfPPCA<T>::pTnComputation()
690 {
691  int nbInd = _individuals->getSizeX() ;
692  int nbIndTemp = 0 ;
693  double sumPtn = 0., prodPtn = 0., meanPtn = 0. , varPtn = 0. ;
694  carto::VolumeRef<double> pTni( _individuals->getSizeX() ) ;
695  double Pi = 0. ;
696 
697  _nullPtnIndiv.clear() ;
698 
699  _pTn = 0 ;
700  cout << endl ;
701  for( int n = 0 ; n < nbInd ; ++n ){
702  for( int c = 0 ; c < _nbOfClasses ; ++c ){
703  Pi = _elements[c].getPi() ;
704  pTni = _elements[c].getPtni() ;
705  _pTn[n] += pTni[n] * Pi ;
706  }
707  if( _pTn[n] != 0 )
708  ++nbIndTemp ;
709  else
710  _nullPtnIndiv.push_back( n ) ;
711  }
712 
713  _nbOfRejected = nbInd - nbIndTemp ;
714 
715  cout << "nb d'individus pris en compte = " << nbIndTemp << " - nb d'individus rejetes = "
716  << nbInd - nbIndTemp << endl ;
717 
718  // moyenne et variance des pTn
719  double oldSumLnPTn = _logLikelihood ;
720  _logLikelihood = 0. ;
721  for( int n = 0 ; n < nbInd ; ++n ){
722  if( _pTn[n] != 0. )
723  _logLikelihood += log( _pTn[n] ) ;
724  else
725  _logLikelihood += -100000000. ;
726  sumPtn += _pTn[n] ;
727  prodPtn += _pTn[n] * _pTn[n] ;
728  }
729  meanPtn = sumPtn / nbIndTemp ;
730  varPtn = prodPtn / nbIndTemp ;
731  cout << endl << "CRITERIUM TO MAXIMIZE : LOG LIKELIHOOD = " << _logLikelihood << " and difference = "
732  << (_logLikelihood - oldSumLnPTn)/oldSumLnPTn << " for "
733  << nbInd - nbIndTemp << " rejections" << endl ;
734  cout << "moyenne et variance des p(tn) = " << meanPtn << " - " << varPtn << endl ;
735 
736  return (_logLikelihood - oldSumLnPTn)/oldSumLnPTn ;
737 }
738 
739 template <class T>
740 double
741 MixtureOfPPCA<T>::distComputation()
742 {
743  carto::VolumeRef<double> mean( _individuals->getSizeY() ) ;
744  carto::VolumeRef<double> invCi( _individuals->getSizeY(), _individuals->getSizeY() ) ;
745  carto::VolumeRef<double> centeredIndivn( _individuals->getSizeY() ) ;
746  carto::VolumeRef<double> dist(1, 1) ;
747  carto::VolumeRef<double> Rn( _individuals->getSizeX() ) ;
748  carto::VolumeRef<double> Rni( _individuals->getSizeX(), _nbOfClasses ) ;
749  carto::VolumeRef<int> results( _individuals->getSizeX() ) ;
750  double max = 0., global_dist = 0. ;
751 
752  for( int c = 0 ; c < _nbOfClasses ; ++c ){
753  Rn = _elements[c].getRn() ;
754  for( int ind = 0 ; ind < _individuals->getSizeX() ; ++ind )
755  Rni(ind, c) = Rn[ind] ;
756  }
757 
758  for( int ind = 0 ; ind < _individuals->getSizeX() ; ++ind ){
759  max = Rni( ind, 0 ) ;
760  results[ind] = 1 ;
761  for( int c = 1 ; c < _nbOfClasses ; ++c ){
762  if( Rni(ind, c) > max ){
763  max = Rni( ind, c ) ;
764  results[ind] = c + 1 ;
765  }
766  }
767  mean = _elements[results[ind]-1].getMean() ;
768  invCi = _elements[results[ind]-1].getInvCi() ;
769  for( int t = 0 ; t < _individuals->getSizeY() ; ++t )
770  centeredIndivn[t] = _individuals(ind, t) - mean[t] ;
771  dist = ( centeredIndivn.clone().transpose() ).cross( invCi.cross( centeredIndivn ) ) ;
772  global_dist += dist(0,0) ;
773  }
774 
775  return global_dist ;
776 }
777 
778 
779 template <class T>
780 bool
781 MixtureOfPPCA<T>::classesVisualisation( int nbOfIterations, const string & fileOut, bool theEnd )
782 {
783  bool res = true ;
784 
785 // POUR TOUS LES INDIVIDUS
786  double max = 0. ;
787  carto::VolumeRef<double> score( _individuals->getSizeX(), _nbOfClasses ) ;
788  carto::VolumeRef<int> results( _individuals->getSizeX() ) ;
789  vector< list< int > > tempList( _nbOfClasses ) ;
790 
791  carto::VolumeRef<double> Rni = getRni() ;
792 
793 /* for( int ind = 0 ; ind < _individuals->getSizeX() ; ++ind ){ */
794 /* sumRn = 0. ; */
795 /* for( int c = 0 ; c < _nbOfClasses ; ++c ) */
796 /* sumRn += Rni( ind, c ) ; */
797 /* for( int c = 0 ; c < _nbOfClasses ; ++c ) */
798 /* score( ind, c) = Rni(ind, c) / sumRn ; */
799 /* } */
800 
801 // CALCUL DES CLASSES FINALES
802  for( int ind = 0 ; ind < _individuals->getSizeX() ; ++ind ){
803  max = 0. ;
804  results[ind] = 0 ;
805  for( int c = 0 ; c < _nbOfClasses ; ++c ){
806  if( Rni(ind, c) > max ){
807  max = Rni( ind, c ) ;
808  results[ind] = c + 1 ;
809  }
810  }
811  if( results[ind] > 0 )
812  _finalClasses[results[ind] - 1].push_back( ind ) ;
813  }
814 
815 
816 // SAUVEGARDE DES CLASSES toutes les 10 iterations
817 // carto::VolumeRef< short > resultImage( data->getSizeX(), data->getSizeY(), data->getSizeZ() ) ; // data inconnu ici
818 // resultImage.setVoxelSize( data->getVoxelSize() ) ;
819  if( nbOfIterations% 30 == 0 || theEnd ){
820  carto::VolumeRef< short > resultImage( 128, 128, 111 ) ;
821  carto::VolumeRef< short > rejected( 128, 128, 111 ) ;
822  carto::VolumeRef< float > rniImage( 128, 128, 111, _nbOfClasses ) ;
823  rejected.setVoxelSize( 0.858086, 0.858086, 2.425 ) ;
824  rejected = 0 ;
825  rniImage.setVoxelSize( 0.858086, 0.858086, 2.425 ) ;
826  rniImage = 0. ;
827  resultImage.setVoxelSize( 0.858086, 0.858086, 2.425 ) ;
828  resultImage = 0 ;
829  Point3d theVoxel ;
830  // recuperer le x,y,z du voxel dont la valeur est *iter
831  for( short c = 0 ; c < _nbOfClasses ; ++c ){
832  list<int>::const_iterator iter( _finalClasses[c].begin() ),
833  last( _finalClasses[c].end() ) ;
834  while( iter != last ){
835  theVoxel = _indPosVector[*iter] ;
836  resultImage( theVoxel ) = c + 1 ;
837  rniImage(theVoxel[0], theVoxel[1], theVoxel[2], c) = Rni( *iter, c ) ;
838  ++iter ;
839  }
840  }
841 
842  list<int>::const_iterator it( _nullPtnIndiv.begin() ),
843  lt( _nullPtnIndiv.end() ) ;
844  while( it != lt ){
845  theVoxel = _indPosVector[*it] ;
846  rejected( theVoxel ) = 1 ;
847  ++it ;
848  }
849 
850  string iter_nb = carto::toString( nbOfIterations ) ;
851  Writer< carto::VolumeRef< short > > writer( "iter" + iter_nb + "_classes_" + fileOut ) ;
852  writer.write( resultImage ) ;
853 
854  Writer< carto::VolumeRef< short > > writer3( "iter" + iter_nb + "_rejected__" + fileOut ) ;
855  writer3.write( rejected ) ;
856 
857  Writer< carto::VolumeRef< float > > writer2( "iter" + iter_nb + "_rni.ima" ) ;
858  writer2.write( rniImage ) ;
859  }
860 
861  cout << endl << "Nb d'individus a l'iteration " << nbOfIterations << " ds chaque classe = " ;
862  for( int i = 0 ; i < _nbOfClasses ; ++i )
863  cout << _finalClasses[i].size() << " " ;
864  cout << endl ;
865 
866  // A LA DERNIERE ITERATION
867  if( theEnd ){
868 
869  // sauvegarder la matrice des individus _individuals (carto::VolumeRef) et les classes finales (1 bucketMap)
870  Writer< carto::VolumeRef<T> > writerMatrix( "indivMatrix" ) ;
871  Finder f ;
872  string fformat = f.format() ;
873  if( !writerMatrix.write( _individuals, false, &fformat ) )
874  return( false ) ;
875 
876 // AimsBucket< int > bClasses ;
877  AimsBucket< short > bClasses ;
878  int count ;
879 
880  for( int c = 0 ; c < _nbOfClasses ; ++c ){
881  count = 0 ;
882  list<int>::const_iterator iter( _finalClasses[c].begin() ),
883  last( _finalClasses[c].end() ) ;
884  while( iter != last ){
885 // AimsBucketItem< int > item ;
886  AimsBucketItem< short > item ;
887  item.location() = Point3d( count, count, count ) ;
888  item.value() = *iter ;
889  bClasses[c].push_back( item ) ;
890  ++iter ;
891  ++count ;
892  }
893  }
894 // Writer< AimsBucket< int > > writerClasses( fileOut + "_finalClasses" ) ;
895  Writer< AimsBucket< short > > writerClasses( fileOut + "_finalClasses" ) ;
896  writerClasses.write( bClasses, true, &fformat ) ;
897  bClasses.clear() ;
898 
899  carto::VolumeRef<double> classCorrelationMatrix( _nbOfClasses, _nbOfClasses ) ;
900  vector<double> classesNorm( _nbOfClasses, 0. ) ;
901  multimap<double, Point2d> fusion ;
902  multimap<double, int> split ;
903 
904  for( int i = 0 ; i < _nbOfClasses ; ++i )
905  for( int n = 0 ; n < _individuals->getSizeX() ; ++n ) {
906  for( int j = i ; j < _nbOfClasses ; ++j )
907  classCorrelationMatrix(i, j) += Rni( n, i ) * Rni(n, j ) ;
908 
909  classesNorm[i] += Rni( n, i ) * Rni( n, i ) ;
910  }
911 
912  for( int i = 0 ; i < _nbOfClasses ; ++i ){
913  classesNorm[i] = sqrt( classesNorm[i] ) ;
914  double meanNorm = 0. ;
915  for( int t = 0 ; t < _elements[i].getMean()->getSizeX() ; ++t )
916  meanNorm += _elements[i].getMean()[t] ;
917 
918  split.insert( pair<double, int>( _elements[i].getSigma2() / meanNorm, i) ) ;
919  }
920  for( int i = 0 ; i < _nbOfClasses ; ++i )
921  for( int j = i ; j < _nbOfClasses ; ++j ){
922  if ( classesNorm[i] * classesNorm[j] > 0 )
923  classCorrelationMatrix(i, j) /= classesNorm[i] * classesNorm[j] ;
924  else
925  classCorrelationMatrix(i, j) = 0. ;
926 
927  // To get a cost, we take the opposite of correlation matrix
928  classCorrelationMatrix(j, i) = classCorrelationMatrix(i, j) ;
929  if( i != j )
930  fusion.insert( pair<double, Point2d> ( classCorrelationMatrix(i, j), Point2d(i, j) ) ) ;
931  }
932 
933  cout << "Classes to fusion : " << endl ;
934 
935  count = 0 ;
936  for( multimap<double, Point2d>::reverse_iterator rit = fusion.rbegin() ; rit != fusion.rend() && count < 5 ; ++rit, ++count )
937  cout << (rit->second)[0] << " avec " << (rit->second)[1] << " : " << rit->first << endl ;
938 
939  cout << "Classes to separe : " << endl ;
940  count = 0 ;
941  for( multimap<double, int>::reverse_iterator rit = split.rbegin() ; rit != split.rend() && count < 5 ; ++rit, ++count )
942  cout << rit->second << " : " << rit->first << endl ;
943  }
944  return res ;
945 }
946 
947 
948 template <class T>
949 bool
950 MixtureOfPPCA<T>::distMatrixComputation()
951 {
952  bool res = true ;
953  int count ;
954  double sumDist ;
955  carto::VolumeRef<double> meanDistMatrix( _nbOfClasses, _nbOfClasses ) ;
956  multimap< float, Point2d > fScore ;
957  multimap< float, Point2d > meanFScore ;
958 
959  for( int fromC = 0 ; fromC < _nbOfClasses ; ++fromC ){
960  for( int toC = 0 ; toC < _nbOfClasses ; ++toC ){
961  count = 0 ;
962  sumDist = 0. ;
963  list<int>::const_iterator iter( _finalClasses[fromC].begin() ),
964  last( _finalClasses[fromC].end() ) ;
965  while( iter != last ){
966  sumDist += _distToClasses( *iter, toC ) ;
967  ++count ;
968  ++iter ;
969  }
970  meanDistMatrix( fromC, toC ) = sumDist / count ;
971  }
972  }
973 /* cout << endl << "matrice des distances moyennes inter classes = " */
974 /* << endl << meanDistMatrix << endl ; */
975 
976  // remplissage de la map de fusion fScore
977  for( int i = 0 ; i < _nbOfClasses ; ++i )
978  for( int j = 0 ; j < _nbOfClasses ; ++j )
979  fScore.insert( pair< float, Point2d >( meanDistMatrix(i, j), Point2d(i, j) ) ) ;
980 
981  // remplissage de la map meanFScore
982  for( int i = 0 ; i < _nbOfClasses ; ++i )
983  for( int j = i+1 ; j < _nbOfClasses ; ++j )
984  meanFScore.insert( pair< float, Point2d >( ( meanDistMatrix(i, j) + meanDistMatrix(j, i) )/2.,
985  Point2d(i, j) ) ) ;
986 
987  cout << endl << "elements de la map des scores = " << endl;
988  multimap< float, Point2d >::iterator iter( fScore.begin() ), last( fScore.end() ) ;
989  while( iter != last ){
990  if( (*iter).first > 1.5 && (*iter).first < 100. ){
991  cout << " [" << (*iter).first << ", " << (*iter).second << "]" ;
992  cout << " sigmas = " << _sigma2Matrix[ (*iter).second[0] ] << " et "
993  << _sigma2Matrix[ (*iter).second[1] ] << endl;
994  }
995  ++iter ;
996  }
997 
998  cout << endl << "elements de la map des scores moyens 2 a 2 = " << endl;
999  multimap< float, Point2d >::iterator iterM( meanFScore.begin() ), lastM( meanFScore.end() ) ;
1000  while( iterM != lastM ){
1001  if( exp( -0.5 * (*iterM).first ) > 0.01 ){
1002  cout << " [" << (*iterM).first << ", " << (*iterM).second << "]" ;
1003  cout << " exp(-score/2) = " << exp( -0.5 * (*iterM).first ) << endl ;
1004  }
1005  ++iterM ;
1006  }
1007 
1008 
1009  return res;
1010 }
1011 
1012 
1013 template <class T>
1014 MixtureOfPPCA<T>::MixtureOfPPCA(
1015  int nbOfClasses, int significantNumberOfVp,
1016  int maxNbOfIterations,
1017  const carto::rc_ptr<carto::Volume<T> >& individuals, const std::vector< Point3d > indPosVector,
1018  const std::vector< std::list <int> >& initialClasses,
1019  const std::string & fileOut, int runNb, int iterationToUseOnlyCorrelatedIndiv ) :
1020  _nbOfClasses( nbOfClasses ), _valid( true ), _significantNumberOfEigenValues( significantNumberOfVp ), _maxNbOfIterations( maxNbOfIterations ), _individuals( individuals ), _indPosVector( indPosVector ), _fileOut( fileOut ), _runNb( runNb ), _itToUseOnlyCorrelatedIndiv(iterationToUseOnlyCorrelatedIndiv)
1021 {
1022  cout << "Entering MixtureOfPPCA Constructor" << endl ;
1023 
1024  double Pi_init = 1. / (double)_nbOfClasses ;
1025  _elements.reserve( _nbOfClasses ) ;
1026  _pTn = carto::VolumeRef<double>( _individuals->getSizeX() ) ;
1027  _pTn = 1 ;
1028  _noiseRef = 1. ;
1029  _sigma2init = vector<double>( _nbOfClasses ) ;
1030 
1031  cout << endl ;
1032  std::vector< std::list<int> > initClasses( initialClasses ) ;
1033 
1034  if( initClasses.size() == 0 ){
1035  initClasses = std::vector< std::list<int> >( _nbOfClasses ) ;
1036  for( int i = 0 ; i < _individuals->getSizeX() ; ++i ){
1037  double aux = UniformRandom( 0., 0.99999 ) ;
1038  int c = int( aux * _nbOfClasses ) ;
1039  initClasses[c].push_back( i ) ;
1040  }
1041  std::cout << std::endl ;
1042  }
1043 
1044  for( int c = 0 ; c < _nbOfClasses ; ++c ){
1045  std::cout << "Class " << c << " : ";
1046  for( std::list<int>::iterator iter = initClasses[c].begin() ; iter < initClasses[c].end() ; ++iter )
1047  std::cout << "\t" << (*iter) ;
1048  std::cout << std::endl ;
1049  }
1050 
1051  for( int c = 0 ; c < _nbOfClasses ; ++c ){
1052  if( int( initClasses[c].size() ) > _individuals->getSizeY() ) {
1053  PpcaAnalyserElement el( _significantNumberOfEigenValues ) ;
1054  cout << endl << "Initialization class " << c << " ..." ;
1055 
1056 /* cout << endl << "numeros des indivs de la classe " << c << " : " ; */
1057 /* count = 0 ; */
1058 /* list<int>::const_iterator iter( initClasses[c].begin() ), */
1059 /* last( initClasses[c].end() ) ; */
1060 /* while( iter != last && count < 10 ){ */
1061 /* cout << *iter << " " ; */
1062 /* ++count ; */
1063 /* ++iter ; */
1064 /* } */
1065 
1066  // AJOUT: calcul d'un facteur noiseRef = bruit de la classe 0 pour normaliser normFactor
1067  if( c == 0 ){
1068  el.init( initClasses[c], Pi_init, individuals, _noiseRef ) ;
1069  _noiseRef = el.getSigma2() ;
1070  cout << "Noise Ref = " << setprecision(3) << _noiseRef << endl ;
1071  }
1072  el.init( initClasses[c], Pi_init, individuals, _noiseRef ) ;
1073  cout << "Class " << c << " initialized !" << endl ;
1074  _elements.push_back( el ) ;
1075  }
1076  else
1077  cerr << "Not enough data for element " << c << endl ;
1078  }
1079 
1080 }
1081 
1082 
1083 template <class T>
1084 bool
1085 //const vector< vector< list< int > > >&
1086 MixtureOfPPCA<T>::doIt()
1087 {
1088  int nbOfIterations = 0 ;
1089  bool res ;
1090 // carto::VolumeRef<double> distToClass( _individuals->getSizeX() ) ;
1091 // _distToClasses = carto::VolumeRef<double>( _individuals->getSizeX(), _nbOfClasses ) ;
1092  _finalClasses = vector< list< int > >( _nbOfClasses ) ;
1093  _sigma2Matrix = carto::VolumeRef<double>( _nbOfClasses ) ;
1094  double sumOfEnergies, criteria ;
1095 
1096  carto::VolumeRef<double> imageOfMean( _nbOfClasses, _maxNbOfIterations, _individuals->getSizeY() ) ;
1097  carto::VolumeRef<double> imageOfSigma2( _nbOfClasses, _maxNbOfIterations ) ;
1098  carto::VolumeRef<double> imageOfPi( _nbOfClasses, _maxNbOfIterations ) ;
1099  carto::VolumeRef<double> imageOfRniDiff( _nbOfClasses, _maxNbOfIterations ) ;
1100  carto::VolumeRef<double> imageOfEnergy( _maxNbOfIterations ) ;
1101  carto::VolumeRef<double> imageOfLogLikelihood( _maxNbOfIterations ) ;
1102  carto::VolumeRef<double> imageOfLogLikelihoodProgression( _maxNbOfIterations ) ;
1103 
1104  _nbOfRejected = 0 ;
1105  int oldNbOfRejected = 0 ;
1106  bool notYetTheEnd = true ;
1107 
1108  do{
1109  for( int i = 0 ; i < _nbOfClasses ; ++i )
1110  _finalClasses[i].clear() ;
1111  cout << endl << "Iteration n�: " << nbOfIterations /*<< endl*/ ;
1112 
1113  // 1ERE PARTIE
1114  cout << endl << "Calcul des p(tn|i) pour chaque classe... " << endl ;
1115  for( int i = 0 ; i < _nbOfClasses ; ++i ){
1116  cerr << "Class " << i << ": " ;
1117  if( _elements[i].isValid() )
1118  _elements[i].newStep1( _individuals, 0 /*(_itToUseOnlyCorrelatedIndiv < 0 ? -1 : nbOfIterations > _itToUseOnlyCorrelatedIndiv )*/ ) ;
1119  }
1120 
1121  oldNbOfRejected = _nbOfRejected ;
1122 
1123  // INTERMEDIAIRE - calcul de p(tn) = sum( p(tn/i)*p(i) )
1124  cout << endl << "Calcul des p(tn)..." ;
1125  imageOfLogLikelihoodProgression(nbOfIterations) = pTnComputation() ;
1126  imageOfLogLikelihood(nbOfIterations) = _logLikelihood ;
1127 
1128  // 2EME PARTIE
1129  cout << endl << "Calcul des autres parametres pour chaque classe..." ;
1130 
1131  //_noiseRef = _elements[0].newStep2( _pTn, _individuals, 1. ) ;
1132  criteria = 0. ;
1133  double oldSumOfEnergies = sumOfEnergies ;
1134  sumOfEnergies = 0. ;
1135  for( int i = 0 ; i < _nbOfClasses ; ++i ){
1136  cout << "Class " << i << ": " ;
1137  if( _elements[i].isValid() ){
1138  _elements[i].newStep2( _pTn, _individuals, /*_noiseRef*/1000000. ) ;
1139  //_sigma2Matrix[i] = _elements[i].getSigma2() ;
1140  //cout << "ok !" << endl ;
1141  sumOfEnergies += _elements[i].getEnergy() ;
1142  criteria += _elements[i].getSumDiff2Rni() ;
1143 
1144  //------ AJOUT 12/11 pour visualiser les param�tres sous forme d'carto::VolumeRef -----
1145  for( int t = 0 ; t < _individuals->getSizeY() ; ++t )
1146  imageOfMean( i, nbOfIterations, t ) = ( _elements[i].getMean() )[t] ;
1147  imageOfSigma2( i, nbOfIterations ) = _elements[i].getSigma2() ;
1148  imageOfPi( i, nbOfIterations ) = _elements[i].getPi() ;
1149  imageOfRniDiff( i, nbOfIterations ) = _elements[i].getSumDiff2Rni() ;
1150  }
1151  }
1152 
1153  cout << endl << "Critere d'arret = " << criteria << " or " << (oldSumOfEnergies - sumOfEnergies)/_elements.size() << endl ;
1154  imageOfEnergy( nbOfIterations ) = sumOfEnergies ;
1155 
1156 /* Writer < carto::VolumeRef<double> > wr( "Mean3D.ima" ) ; */
1157 /* wr.write( imageOfMean ) ; */
1158 /* Writer < carto::VolumeRef<double> > wr2( "Sigma2D.ima" ) ; */
1159 /* wr2.write( imageOfSigma2 ) ; */
1160 /* Writer < carto::VolumeRef<double> > wr3( "Pi2D.ima" ) ; */
1161 /* wr3.write( imageOfPi ) ; */
1162 /* Writer < carto::VolumeRef<double> > wr4( "RniDiff2D.ima" ) ; */
1163 /* wr4.write( imageOfRniDiff ) ; */
1164 /* Writer < carto::VolumeRef<double> > wr5( "Energy1D.ima" ) ; */
1165 /* wr5.write( imageOfEnergy ) ; */
1166 /* Writer < carto::VolumeRef<double> > wr6( "logLikelihood.ima" ) ; */
1167 /* wr6.write( imageOfLogLikelihood ) ; */
1168 /* Writer < carto::VolumeRef<double> > wr7( "logLikelihoodProgression.ima" ) ; */
1169 /* wr7.write( imageOfLogLikelihoodProgression ) ; */
1170 
1171  notYetTheEnd = (nbOfIterations < _maxNbOfIterations) &&
1172  ( (nbOfIterations < 10 ) ||
1173  !( (oldNbOfRejected == _nbOfRejected) && (fabs( imageOfLogLikelihoodProgression(nbOfIterations) ) < 0.0000001) ) ) ;
1174 
1175  // fonction pour calculer et visualiser les classes
1176  res = classesVisualisation( nbOfIterations, _fileOut, !notYetTheEnd ) ;
1177 
1178  ++nbOfIterations ;
1179  } while( notYetTheEnd ) ;
1180 
1181  // matrice des distances de chaque individu � chaque classe
1182  /* for( int i = 0 ; i < _nbOfClasses ; ++i ){ */
1183  /* distToClass = _elements[i].getDist() ; */
1184  /* for( int n = 0 ; n < distToClass->getSizeX() ; ++n ){ */
1185  /* // distToClass d�j� au carr� - division par la dimension */
1186  /* _distToClasses(n, i) = distToClass[n] / _individuals->getSizeY() ; */
1187  /* } */
1188  /* } */
1189 
1190  // _valid = distMatrixComputation() ;
1191 
1192  // return _finalClasses ;
1193  // return _nullPtnIndiv ;
1194  return _valid ;
1195 }
1196 
1197 
1198 
1199 template <class T>
1200 carto::VolumeRef<double>
1201 MixtureOfPPCA<T>::getRni()
1202 {
1203  carto::VolumeRef<double> Rn( _individuals->getSizeX() ) ;
1204  carto::VolumeRef<double> RnAll( _individuals->getSizeX(), _nbOfClasses ) ;
1205  for( int c = 0 ; c < _nbOfClasses ; ++c ){
1206  Rn = _elements[c].getRn() ;
1207  for( int ind = 0 ; ind < _individuals->getSizeX() ; ++ind )
1208  RnAll(ind, c) = Rn[ind] ;
1209  }
1210 
1211  return RnAll ;
1212 }
1213 
1214 }
1215 
1216 #endif
ForEach2d
#define ForEach2d(thing, x, y)
AimsBucketItem::location
const AimsVector< short, 3 > & location() const
AimsBucketItem::value
const T & value() const
AimsBucket::push_back
void push_back(const AimsBucketItem< T > &item)
AimsSVD
Definition: svd.h:59
AimsSVD::sort
void sort(carto::VolumeRef< T > &, carto::VolumeRef< T > &, carto::VolumeRef< T > *v=NULL)
sort the U and V matrices and the W vector in decreasing order
AimsSVD::doit
carto::VolumeRef< T > doit(carto::VolumeRef< T > &, carto::VolumeRef< T > *v=NULL)
Singular Value Decomposition.
AimsSVD::setReturnType
void setReturnType(SVDReturnType rt)
Definition: svd.h:76
aims::Finder::format
std::string format() const
aims::MixtureOfPPCA::getRni
carto::VolumeRef< double > getRni()
Definition: mixtureOfPpca_d.h:1201
aims::MixtureOfPPCA::classesVisualisation
bool classesVisualisation(int nbOfIterations, const std::string &fileOut, bool theEnd)
Definition: mixtureOfPpca_d.h:781
aims::MixtureOfPPCA::MixtureOfPPCA
MixtureOfPPCA(int nbOfClasses, int significantNumberOfVp, int maxNbOfIterations, const carto::rc_ptr< carto::Volume< T > > &individuals, const std::vector< Point3d > indPosVector, const std::vector< std::list< int > > &initialClasses, const std::string &fileOut, int runNb, int iterationToUseOnlyCorrelatedIndiv=false)
Definition: mixtureOfPpca_d.h:1014
aims::PpcaAnalyserElement::getSigma2
double getSigma2() const
Definition: mixtureOfPpca.h:72
aims::PpcaAnalyserElement::init
void init(const std::list< int > &selectedIndividuals, double initialPi, const carto::rc_ptr< carto::Volume< T > > &individuals, double noiseRef=1.)
Definition: mixtureOfPpca_d.h:76
aims::PpcaAnalyserElement::newStep1
bool newStep1(const carto::rc_ptr< carto::Volume< T > > &indivMatrix, bool useOnlyCorrIndiv)
Definition: mixtureOfPpca_d.h:218
aims::PpcaAnalyserElement::newStep2
double newStep2(carto::rc_ptr< carto::Volume< double > > pTn, const carto::rc_ptr< carto::Volume< T > > &indivMatrix, double noiseRef=1.)
Definition: mixtureOfPpca_d.h:319
aims::PpcaAnalyserElement::PpcaAnalyserElement
PpcaAnalyserElement(int significantNumberOfVp, bool useOnlyCorrIndiv=false)
Definition: mixtureOfPpca_d.h:63
aims::Writer::write
virtual bool write(const T &obj, bool ascii=false, const std::string *format=0)
carto::VolumeRef::setVoxelSize
void setVoxelSize(float vx, float vy=1., float vz=1., float vt=1.)
carto::VolumeRef::getSizeY
int getSizeY() const
carto::VolumeRef::getSizeX
int getSizeX() const
carto::rc_ptr
aims::meshdistance::max
float max(float x, float y)
Definition: thickness.h:98
aims::meshdistance::cross
Point3df cross(const Point3df &A, const Point3df &B)
Definition: projection.h:241
split
std::vector< std::string > split(const std::string &text, const std::string &sep)
carto::toString
std::string toString(const T &object)
UniformRandom
double UniformRandom()
Uniform distribution between 0.0 and 1.0.
string_conversion.h