aimsalgo 6.0.0
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
58namespace aims
59{
60
61double * PpcaAnalyserElement::_exp = 0 ;
62
63PpcaAnalyserElement::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
74template <class T>
75void
76PpcaAnalyserElement::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
97template <class T>
98void
99PpcaAnalyserElement::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
216template <class T>
217bool
218PpcaAnalyserElement::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
317template <class T>
318double
319PpcaAnalyserElement::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
687template <class T>
688double
689MixtureOfPPCA<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
739template <class T>
740double
741MixtureOfPPCA<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
779template <class T>
780bool
781MixtureOfPPCA<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
948template <class T>
949bool
950MixtureOfPPCA<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
1013template <class T>
1014MixtureOfPPCA<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
1083template <class T>
1084bool
1085//const vector< vector< list< int > > >&
1086MixtureOfPPCA<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
1199template <class T>
1200carto::VolumeRef<double>
1201MixtureOfPPCA<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::value
const T & value() const
AimsBucketItem::location
const AimsVector< short, 3 > & location() 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::MatrixOfSingularValues
@ MatrixOfSingularValues
Definition svd.h:64
AimsSVD::VectorOfSingularValues
@ VectorOfSingularValues
Definition svd.h:65
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::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::transpose
AimsData< T > transpose(const AimsData< T > &thing)
carto::toString
std::string toString(const T &object)
max
T max(const Volume< T > &vol)
UniformRandom
double UniformRandom()
Uniform distribution between 0.0 and 1.0.
Point3d
AimsVector< int16_t, 3 > Point3d
Point2d
AimsVector< int16_t, 2 > Point2d