cah_d.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */

#include <iostream>
#include <math.h>
#include <aims/math/pca_d.h>
#include <aims/math/ppca.h>
#include <aims/io/writer.h>
#include <aims/connectivity/connectivity.h> 
#include <aims/classification/individuals_d.h>
#include <bioprocessing/classification/cah.h> 

template <typename T>
BaryElement<T>::BaryElement( int card, const aims::Individuals<float>& bary, 
                             const std::set<int>& constitutedFrom, int number ) :
  Element<T>( number, card, constitutedFrom ), _bary(bary)
{}

template <typename T>
BaryElement<T>::BaryElement( const std::vector<std::list< Point3d > >& classes, 
                             const AimsData<T>& data, const std::set<int>& constitutedFrom, 
                             int number ) :
  Element<T>( number, 0, constitutedFrom )
{

  this->_card = 0 ;
  std::set<int>::const_iterator aggIt( this->_constitutedFrom.begin() ), aggLast( this->_constitutedFrom.end() ) ;
  
  while( aggIt != aggLast ) {
    this->_card += classes[*aggIt].size() ;
    ++aggIt ;
  }
  
  computeBary( classes, data ) ;    
}

template <typename T>
void
BaryElement<T>::computeBary( const std::vector<std::list< Point3d > >& classes,
                             const AimsData<T>& data )
{
  
  int indivSize = data.dimT() ;
/*   _bary.position() = Point3df(0., 0., 0.) ; */
/*   _bary.value() = std::vector<float>( indivSize ) ; */

  Point3df position(0., 0., 0.) ;
  std::vector<float> value( indivSize ) ;
  
  std::set<int>::iterator constFromIter(this->_constitutedFrom.begin()), constLast(this->_constitutedFrom.end()) ;
  
  int nbInd = 0 ;
  while( constFromIter != constLast ){
    std::list<Point3d>::const_iterator iter(classes[*constFromIter].begin()), last(classes[*constFromIter].end()) ;
    while( iter != last ){
      position += Point3df( (*iter)[0], (*iter)[1], (*iter)[2] ) ;
      for( int j = 0 ; j < indivSize ; ++j )
        value[j] += data((*iter)[0], (*iter)[1], (*iter)[2], j) ;
      ++iter ;
      ++nbInd ;
    }
    ++constFromIter ;
  }
  if( nbInd){
    float invNbInd = 1. / nbInd ;
    for( int j = 0 ; j < indivSize ; ++j )
      value[j] *= invNbInd ;
    position = position * invNbInd ;
  }
  _bary = aims::Individuals<float>(position, value) ;
}

template <typename T>
float 
BaryElement<T>::distanceTo( Element<T> * other, 
                            const std::vector< std::list< Point3d > >&,
                            const AimsData<T>&, const AimsData<float>& )
{
  BaryElement<T> * otherB = dynamic_cast<BaryElement<T>*>(other) ;
  unsigned int size = baryValue().size() ;
  if( otherB->baryValue().size() != size )
    throw std::runtime_error("Inconsistant data !") ;
  float d = 0. ;
  for( unsigned int i = 0 ; i < size ; ++i )
    d += (_bary.value()[i] - otherB->baryValue()[i])*(_bary.value()[i] - otherB->baryValue()[i]);
  d *= this->_card * otherB->card() /( this->_card + otherB->card() ) ;
  
  return d ;
}

template <typename T>
Element<T> * 
BaryElement<T>::agregateTo( Element<T> * other, 
                            const std::vector< std::list< Point3d > >&,
                            const AimsData<T>&, const AimsData<float>&, int newNb )
{
  BaryElement<T> * otherB = dynamic_cast<BaryElement<T>*>(other) ;
  unsigned int size = baryValue().size() ;
  if( otherB->baryValue().size() != size )
    throw std::runtime_error("Inconsistant data !") ;
  
  // Agglomeration
  std::set<int> constFrom( this->constitutedFrom() ) ;
  std::set<int>::const_iterator cIt(otherB->constitutedFrom().begin() ), 
    cLast(otherB->constitutedFrom().end() ) ;
  while( cIt != cLast ){
    constFrom.insert( *cIt ) ;
    ++cIt ;
  }
  
  // Calcul du nouveau barycentre des 2 classes fusionnees
  Point3df pos ;
  std::vector<float> val( size ) ;
  int newCard = this->_card + otherB->card() ;
  float invNewCard = 1./newCard ;
  for(unsigned int i=0;i<_bary.value().size();i++)
    val[i] = (this->_card*_bary.value()[i] + 
              otherB->card()*otherB->baryValue()[i]) * invNewCard ;
  pos = (float(this->_card)*_bary.position() + 
         float(otherB->card())*otherB->baryPosition()) * invNewCard ;
  
  aims::Individuals<float> newBary( pos, val );
  
  // Instanciation du nouvel element
  Element<T> *nouveau = new BaryElement<T>( newCard, newBary, constFrom, newNb );
  return nouveau ;
}

template <typename T>
PcaElement<T>::PcaElement( const std::vector<std::list< Point3d > >& classes, 
                           const AimsData<T>& data, int numberOfSignificantVps,
                           const std::set<int>& constitutedFrom, int number ) :
  Element<T>( number, 0, constitutedFrom ), _meanReconstructionError(-1.), _varReconstructionError(-1.),
_unreconstructedInertia(-1.), _numberOfSignificantVps(numberOfSignificantVps)
{
  this->_card = 0 ;
  std::set<int>::const_iterator aggIt( this->_constitutedFrom.begin() ), aggLast( this->_constitutedFrom.end() ) ;
  
  while( aggIt != aggLast ) {
    this->_card += classes[*aggIt].size() ;
    ++aggIt ;
  }

  computeReconstructionError( classes, data ) ;
}

template <typename T>
float 
PcaElement<T>::meanReconstructionError()
{
  return _meanReconstructionError ;
}

template <typename T>
float 
PcaElement<T>::varReconstructionError()
{
  return _varReconstructionError ;
}

template <typename T>
float 
PcaElement<T>::unreconstructedInertia()
{
  return _unreconstructedInertia ;
}

template <typename T>
void
PcaElement<T>::computeReconstructionError( const std::vector< std::list< Point3d > >& classes,
                                           const AimsData<T>& data )
{
  _meanReconstructionError = 0. ;
  _varReconstructionError = 0. ;
  std::set<int>::const_iterator aggIter( this->_constitutedFrom.begin() ), aggLast( this->_constitutedFrom.end() ) ;
  
  std::list<Point3d> agregatePoints ;
  while( aggIter != aggLast ){
    std::list<Point3d>::const_iterator lIter(classes[*aggIter].begin()), lLast(classes[*aggIter].end() ) ;
    while( lIter != lLast ){
      agregatePoints.push_back( *lIter ) ;
      ++lIter ;
    }
    ++aggIter ;
  }

  if( int(agregatePoints.size()) > data.dimT() ){
/*     AimsPCA pca( _numberOfSignificantVps ) ; */
/*     pca.doIt( agregatePoints, data ) ; */
/*     _unreconstructedInertia = pca.noiseInertia() ; */
    AimsPCA pca( _numberOfSignificantVps ) ;
    pca.doIt( agregatePoints, data ) ;
    AimsData<float> errorMatrix = pca.reconstructionErrorMatrix( ) ; 

    _mean = pca.mean() ;
    std::list<Point3d>::iterator iter( agregatePoints.begin() ), last( agregatePoints.end() ) ;
    AimsData<float> indiv( 1, data.dimT() ), indivTr(data.dimT()) ;
    while( iter != last ){
      float norm2 = 0. ;
      for ( int t = 0 ; t < data.dimT() ; ++t ){
        indiv(0,t) = data( (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        indivTr(t,0) = data(  (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        
        norm2 += data( (*iter)[0], (*iter)[1], (*iter)[2], t )*
          data( (*iter)[0], (*iter)[1], (*iter)[2], t ) ;
      }
      float val = indiv.cross( errorMatrix.cross(indivTr) )(0,0) / norm2 ;
      _meanReconstructionError += val ;
      _varReconstructionError += val * val ;
      ++iter ;
    }

    _meanReconstructionError /= double(agregatePoints.size()) ;
    _varReconstructionError = _varReconstructionError / double(agregatePoints.size()) - 
      _meanReconstructionError * _meanReconstructionError ;
    if( this->_constitutedFrom.size() == 1 )
      _errorMatrix = errorMatrix ;
      
    //_unreconstructedInertia = pca.noiseInertia() ;
  } else {
    _unreconstructedInertia = 1. ;
    _meanReconstructionError = 1. ;
    _varReconstructionError = 1. ;
  }
}

template <typename T>
float 
PcaElement<T>::distanceTo( Element<T> * other,
                           const std::vector< std::list< Point3d > >& classes,
                           const AimsData<T>& data, const AimsData<float>& initDistances )
{
  Element<T> * fusionEl = agregateTo( other, classes, data, initDistances ) ;
  PcaElement<T> * fusionPcaEl = dynamic_cast< PcaElement<T> * >( fusionEl ) ;
  PcaElement<T> * otherPcaEl = dynamic_cast< PcaElement<T> * >( other ) ;
  
  float distance = fusionPcaEl->card() * fusionPcaEl->meanReconstructionError() /
    ( this->card() * meanReconstructionError() + otherPcaEl->card() * otherPcaEl->meanReconstructionError() ) ;
  
  delete fusionEl ;
  return distance ;
}

template <typename T>
Element<T> * 
PcaElement<T>::agregateTo( Element<T> * other, 
                           const std::vector< std::list< Point3d > >& classes,
                           const AimsData<T>& data, const AimsData<float>&, int newNb  )
{
/*   BaryElement<T> * otherB = dynamic_cast<BaryElement<T>*>(other) ; */
  // Agglomeration
  std::set<int> constFrom( this->constitutedFrom() ) ;
  std::set<int>::const_iterator cIt(other->constitutedFrom().begin() ), 
    cLast(other->constitutedFrom().end() ) ;
  while( cIt != cLast ){
    constFrom.insert( *cIt ) ;
    ++cIt ;
  }
  
  Element<T> * newEl = new PcaElement<T>( classes, data, _numberOfSignificantVps, constFrom, newNb ) ;
  return newEl ;
}

template <typename T>
PPcaElement<T>::PPcaElement( const std::vector<std::list< Point3d > >& classes, 
                           const AimsData<T>& data, int numberOfSignificantVps,
                           const std::set<int>& constitutedFrom, int number ) :
  Element<T>( number, 0, constitutedFrom ), _numberOfSignificantVps(numberOfSignificantVps)
{
  this->_card = 0 ;
  std::set<int>::const_iterator aggIt( this->_constitutedFrom.begin() ), aggLast( this->_constitutedFrom.end() ) ;
  
  //std::cout << "Card = " ;
  while( aggIt != aggLast ) {
    this->_card += classes[*aggIt].size() ;
    //std::cout << "\t" << this->_card ;
    ++aggIt ;
  }
  
  _ppcaElement = new aims::ProbabilisticPcaElement( _numberOfSignificantVps ) ;
  
  computeReconstructionError( classes, data ) ;
}

template <typename T>
void
PPcaElement<T>::computeReconstructionError( const std::vector< std::list< Point3d > >& classes,
                                           const AimsData<T>& data )
{
  std::set<int>::const_iterator aggIter( this->_constitutedFrom.begin() ), aggLast( this->_constitutedFrom.end() ) ;
  
  std::list<Point3d> agregatePoints ;
  while( aggIter != aggLast ){
    std::list<Point3d>::const_iterator lIter(classes[*aggIter].begin()), lLast(classes[*aggIter].end() ) ;
    while( lIter != lLast ){
      agregatePoints.push_back( *lIter ) ;
      ++lIter ;
    }
    ++aggIter ;
  }

  if( int(agregatePoints.size()) > data.dimT() )
    _ppcaElement->doIt(agregatePoints, data) ;
}

template <typename T>
float 
PPcaElement<T>::distanceTo( Element<T> * other,
                            const std::vector< std::list< Point3d > >& classes,
                            const AimsData<T>& data, const AimsData<float>& initDistances )
{
  PPcaElement<T> * otherPPca = dynamic_cast< PPcaElement<T> *>(other) ;
  if( this->constitutedFrom().size() == 1 && other->constitutedFrom().size() == 1 ){
    AimsData<double> indiv(data.dimT()) ;
    std::list<Point3d>::const_iterator iter( classes[*(this->constitutedFrom().begin())].begin() ), 
      last( classes[*(this->constitutedFrom().begin())].end() ) ;
    float oneToTwoMeanError = 0. ;
    float twoToOneMeanError = 0. ;
    while( iter != last ){
      float norm2 = 0. ;
      for ( int t = 0 ; t < data.dimT() ; ++t ){
        indiv(t) = data(  (*iter)[0], (*iter)[1], (*iter)[2], t )  ;
        
        norm2 += data( (*iter)[0], (*iter)[1], (*iter)[2], t )*
          data( (*iter)[0], (*iter)[1], (*iter)[2], t ) ;
      }
      oneToTwoMeanError -= otherPPca->ppca()->posteriorProbability( indiv, 1. ) ;
      ++iter ;
    }
    
    oneToTwoMeanError /= float(classes[*(this->constitutedFrom().begin())].size()) ;
/*     std::cout << "Mean Error = " << oneToTwoMeanError << " against " << otherPcaReproc->meanReconstructionError() << std::endl  */;
    
    iter = classes[ *(otherPPca->constitutedFrom().begin())].begin() ;
    last = classes[ *(otherPPca->constitutedFrom().begin())].end() ;
    
    while( iter != last ){
      float norm2 = 0. ;
      for ( int t = 0 ; t < data.dimT() ; ++t ){
        indiv(t) = data(  (*iter)[0], (*iter)[1], (*iter)[2], t ) ;
        
        norm2 += data( (*iter)[0], (*iter)[1], (*iter)[2], t )*
          data( (*iter)[0], (*iter)[1], (*iter)[2], t ) ;
      }
      twoToOneMeanError -= _ppcaElement->posteriorProbability( indiv, 1. ) ;
      ++iter ;
    }
    twoToOneMeanError /= float(classes[ *(otherPPca->constitutedFrom().begin())].size()) ;
/*     std::cout << "Mean Error = " << twoToOneMeanError << " against " << _meanReconstructionError << std::endl ; */
    return std::max( oneToTwoMeanError, twoToOneMeanError ) ;
  }
  
  std::set<int> constFrom = otherPPca->constitutedFrom() ;
  std::set<int>::const_iterator cIter( this->constitutedFrom().begin() ), cLast( this->constitutedFrom().end() ) ;
  while( cIter != cLast ){
    constFrom.insert( *cIter ) ;
    ++cIter ;
  }

/*   std::cout << "Contains" ; */
/*   std::set<int>::const_iterator c3Iter( constFrom.begin() ), c3Last( constFrom.end() ) ; */
/*   while( c3Iter != c3Last ) */
/*     { */
/*       std::cout << "\t" << *c3Iter ;  */
/*       ++c3Iter ; */
/*     }  */

  std::set<int>::const_iterator c1Iter( constFrom.begin() ), c1Last( constFrom.end() ) ;
  
  float maxDistance = 0. ;
  while( c1Iter != c1Last )
    {
      std::set<int>::const_iterator c2Iter( c1Iter ) ;
      ++c2Iter ;
      while( c2Iter != c1Last )
        {
/*        std::cout << *c1Iter << " " << *c2Iter << std::endl ; */
/*        std::cout << "\t( " << *c1Iter << ", " << *c2Iter << " ) = " << initDistances( *c1Iter, *c2Iter )  ; */
          float initDist = initDistances( *c1Iter, *c2Iter ) ;
          if( maxDistance < initDist )
            maxDistance = initDist ;
          ++c2Iter ;
        }
      ++c1Iter ;
    }
  
/*   std::cout << " = " << maxDistance << std::endl ; */
  
  return maxDistance ;
}

template <typename T>
Element<T> * 
PPcaElement<T>::agregateTo( Element<T> * other, 
                           const std::vector< std::list< Point3d > >& classes,
                           const AimsData<T>& data, const AimsData<float>&, int newNb  )
{
/*   BaryElement<T> * otherB = dynamic_cast<BaryElement<T>*>(other) ; */
  // Agglomeration
  std::set<int> constFrom( this->constitutedFrom() ) ;
  std::set<int>::const_iterator cIt(other->constitutedFrom().begin() ), 
    cLast(other->constitutedFrom().end() ) ;
  while( cIt != cLast ){
    constFrom.insert( *cIt ) ;
    ++cIt ;
  }
  
  Element<T> * newEl = new PPcaElement<T>( classes, data, _numberOfSignificantVps, constFrom, newNb ) ;
  return newEl ;
}

template <typename T>
float
SCKBDistanceElement<T>::distanceTo( Element<T> * other, 
                                    const std::vector< std::list< Point3d > >& ,
                                    const AimsData<T>& , const AimsData<float>& )
{
  //std::cout << "Distance To :: begin" << std::endl ;
  double surf1To2 = 0., surf2To1 = 0., dist1To2 = 0., dist2To1 = 0., maxMomentRatio = 1. ;
  //double ratio ;
  for( std::set<int>::iterator itC1 = this->constitutedFrom().begin() ; itC1 != this->constitutedFrom().end() ; ++itC1 )
    for( std::set<int>::iterator itC2 = other->constitutedFrom().begin() ; itC2 != other->constitutedFrom().end() ; ++itC2 ){
      surf1To2 += (*_interactionSurf)(*itC1, *itC2) ;
      surf2To1 += (*_interactionSurf)(*itC2, *itC1) ;
      dist1To2 = std::max(dist1To2, (*_kinModel2by2Distances)(*itC1, *itC2) ) ;
      dist2To1 = std::max(dist2To1, (*_kinModel2by2Distances)(*itC2, *itC1) ) ;
//       if( (*_classMaximumMoment)[*itC1] == (*_classMaximumMoment)[*itC2] ) 
//      ratio = 1. ;
//       else if( (*_classMaximumMoment)[*itC1] == 0.  || (*_classMaximumMoment)[*itC2] == 0. )
//      ratio = 0. ;
//       else{
//      ratio = std::min( (*_classMaximumMoment)[*itC1] / (*_classMaximumMoment)[*itC2], 
//                        (*_classMaximumMoment)[*itC2] / (*_classMaximumMoment)[*itC1] ) ;
//       }
      
//       maxMomentRatio = std::min( maxMomentRatio, ratio ) ;
      
      maxMomentRatio = std::max( maxMomentRatio, (*_maxMomentValueRatios)(*itC1, *itC2) ) ;
    }
  
  double vol1 = 0, vol2 = 0 ;
  for( std::set<int>::iterator itC1 = this->constitutedFrom().begin() ; itC1 != this->constitutedFrom().end() ; ++itC1 )
    vol1 += (*_labelsVolumes)[*itC1] ;
  for( std::set<int>::iterator itC2 = other->constitutedFrom().begin() ; itC2 != other->constitutedFrom().end() ; ++itC2 )
    vol2 += (*_labelsVolumes)[*itC2] ;
  
  double maxSurfVolRatio = std::max( surf1To2 / vol1, surf2To1 / vol2 ) ;
  
  return ( 0.5 * ( dist1To2 + dist2To1 ) - _meanDist2By2 ) / sqrt(_varDist2By2) 
    - _spatialregularizationWeight * ( maxSurfVolRatio - _meanSurfVolRatio ) / sqrt(_varSurfVolRatio)
    + _timeOfMaximumWeight * ( maxMomentRatio - _meanMaxMomentRatio ) / sqrt(_varMaxMomentRatio) ;
//   return 0.5 * ( dist1To2 + dist2To1 ) - _spatialregularizationWeight * maxSurfVolRatio 
//     - _timeOfMaximumWeight * maxMomentRatio ;
}

template <typename T>
Element<T> *
SCKBDistanceElement<T>::agregateTo( Element<T> * other, 
                                    const std::vector< std::list< Point3d > >& ,
                                    const AimsData<T>& , const AimsData<float>& , int newNb )
{
  std::set<int> constFrom( this->constitutedFrom() ) ;
  std::set<int>::const_iterator cIt(other->constitutedFrom().begin() ),
    cLast(other->constitutedFrom().end() ) ;
  while( cIt != cLast ){
    constFrom.insert( *cIt ) ;
    ++cIt ;
  }
  
  return new SCKBDistanceElement<T>( _maxMomentValueRatios, _kinModel2by2Distances, _interactionSurf, _labelsVolumes,
                                     _spatialregularizationWeight, _timeOfMaximumWeight, 
                                     _meanDist2By2, _varDist2By2, _meanSurfVolRatio, _varSurfVolRatio, 
                                     _meanMaxMomentRatio, _varMaxMomentRatio,
                                     constFrom, newNb ) ;
}

template <typename T>
Max2By2DistanceElement<T>::Max2By2DistanceElement( const std::set<int>& constitutedFrom, int number ) :
  Element<T>( number, 0, constitutedFrom )
{
  //std::cout << "Max2By2DistanceElement" << std::endl ;
  this->_card = 1 ;
}

template <typename T>
float 
Max2By2DistanceElement<T>::distanceTo( Element<T> * other,
                                       const std::vector< std::list< Point3d > >&,
                                       const AimsData<T>&, const AimsData<float>& initDistances )
{

  Max2By2DistanceElement<T> * otherPPca = dynamic_cast< Max2By2DistanceElement<T> *>(other) ;
  if( this->constitutedFrom().size() == 1 && other->constitutedFrom().size() == 1 ){
    return initDistances( *( this->constitutedFrom().begin() ), *( other->constitutedFrom().begin() ) ) ;
  }
  
  std::set<int> constFrom = otherPPca->constitutedFrom() ;
  std::set<int>::const_iterator cIter( this->constitutedFrom().begin() ), cLast( this->constitutedFrom().end() ) ;
  while( cIter != cLast ){
    constFrom.insert( *cIter ) ;
    ++cIter ;
  }

/*   std::cout << "Contains" ; */
/*   std::set<int>::const_iterator c3Iter( constFrom.begin() ), c3Last( constFrom.end() ) ; */
/*   while( c3Iter != c3Last ) */
/*     { */
/*       std::cout << "\t" << *c3Iter ;  */
/*       ++c3Iter ; */
/*     }  */  

  std::set<int>::const_iterator c1Iter( constFrom.begin() ), c1Last( constFrom.end() ) ;
  
  float maxDistance = initDistances( *(constFrom.begin()), *(constFrom.rbegin()) ) ;
  while( c1Iter != c1Last ){
    std::set<int>::const_iterator c2Iter( c1Iter ) ;
    ++c2Iter ;
    while( c2Iter != c1Last )
      {
        /*        std::cout << *c1Iter << " " << *c2Iter << std::endl ; */
        /*        std::cout << "\t( " << *c1Iter << ", " << *c2Iter << " ) = " << initDistances( *c1Iter, *c2Iter )  ; */
        float initDist = initDistances( *c1Iter, *c2Iter ) ;
/*      std::cout << *c1Iter << ", " << *c2Iter << " =  " << initDist << " -> " ; */
        if( maxDistance < initDist )
          maxDistance = initDist ;
/*      std::cout << maxDistance << std::endl ; */
        ++c2Iter ;
      }
    ++c1Iter ;
  }
  
/*   std::cout << " = " << maxDistance << std::endl ; */
  
  return maxDistance ;
}

template <typename T>
Element<T> * 
Max2By2DistanceElement<T>::agregateTo( Element<T> * other, 
                           const std::vector< std::list< Point3d > >&,
                           const AimsData<T>&, const AimsData<float>&, int newNb  )
{
/*   BaryElement<T> * otherB = dynamic_cast<BaryElement<T>*>(other) ; */
  // Agglomeration
  std::set<int> constFrom( this->constitutedFrom() ) ;
  std::set<int>::const_iterator cIt(other->constitutedFrom().begin() ), 
    cLast(other->constitutedFrom().end() ) ;
  while( cIt != cLast ){
    constFrom.insert( *cIt ) ;
    ++cIt ;
  }

  //std::cout << "Max2By2DistanceElement<T>::agregateTo" << std::endl ;
  
  Element<T> * newEl = new Max2By2DistanceElement<T>( constFrom, newNb ) ;
  return newEl ;
}

template <typename T>
PcaReprocReconsErrorElement<T>::PcaReprocReconsErrorElement(const std::vector<std::list< Point3d > >& classes, 
                                                            const AimsData<T>& data, int numberOfSignificantVps,   
                                                            const std::set<int>& constitutedFrom, int number) :
  Element<T>( number, 0, constitutedFrom ), _meanReconstructionError(-1.), _varReconstructionError(-1.),
  _unreconstructedInertia(-1.), _numberOfSignificantVps(numberOfSignificantVps)
{
  this->_card = 0 ;
  std::set<int>::const_iterator aggIt( this->_constitutedFrom.begin() ), aggLast( this->_constitutedFrom.end() ) ;
  
  while( aggIt != aggLast ) {
    this->_card += classes[*aggIt].size() ;
    ++aggIt ;
  }
  if( this->_constitutedFrom.size() == 1 )
    computeReconstructionError( classes, data ) ;
}

template <typename T>
void
PcaReprocReconsErrorElement<T>::computeReconstructionError( const std::vector< std::list< Point3d > >& classes,
                                                            const AimsData<T>& data )
{
  _meanReconstructionError = 0. ;
  _varReconstructionError = 0. ;
  std::set<int>::const_iterator aggIter( this->_constitutedFrom.begin() ), aggLast( this->_constitutedFrom.end() ) ;
  
  std::list<Point3d> agregatePoints ;
  while( aggIter != aggLast ){
    std::list<Point3d>::const_iterator lIter(classes[*aggIter].begin()), lLast(classes[*aggIter].end() ) ;
    while( lIter != lLast ){
      agregatePoints.push_back( *lIter ) ;
      ++lIter ;
    }
    ++aggIter ;
  }

  if( int(agregatePoints.size()) > data.dimT() ){
/*     AimsPCA pca( _numberOfSignificantVps ) ; */
/*     pca.doIt( agregatePoints, data ) ; */
/*     _unreconstructedInertia = pca.noiseInertia() ; */
    AimsPCA pca( _numberOfSignificantVps ) ;
    pca.doIt( agregatePoints, data ) ;
    AimsData<float> errorMatrix = pca.reconstructionErrorMatrix( ) ; 

    _mean = pca.mean() ;
    std::list<Point3d>::iterator iter( agregatePoints.begin() ), last( agregatePoints.end() ) ;
    AimsData<float> indiv( 1, data.dimT() ), indivTr(data.dimT()) ;
    while( iter != last ){
      float norm2 = 0. ;
      for ( int t = 0 ; t < data.dimT() ; ++t ){
        indiv(0,t) = data( (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        indivTr(t,0) = data(  (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        
        norm2 += data( (*iter)[0], (*iter)[1], (*iter)[2], t )*
          data( (*iter)[0], (*iter)[1], (*iter)[2], t ) ;
      }
      float val = indiv.cross( errorMatrix.cross(indivTr) )(0,0) / norm2 ;
      _meanReconstructionError += val ;
      _varReconstructionError += val * val ;
      ++iter ;
    }

    _meanReconstructionError /= double(agregatePoints.size()) ;
    _varReconstructionError = _varReconstructionError / double(agregatePoints.size()) - 
      _meanReconstructionError * _meanReconstructionError ;
    _errorMatrix = errorMatrix ;
      
    //_unreconstructedInertia = pca.noiseInertia() ;
  } else {
    _unreconstructedInertia = 1. ;
    _meanReconstructionError = 1. ;
    _varReconstructionError = 1. ;
  }
}

template <typename T>
float 
PcaReprocReconsErrorElement<T>::meanReconstructionError()
{
  return _meanReconstructionError ;
}

template <typename T>
float 
PcaReprocReconsErrorElement<T>::varReconstructionError()
{
  return _varReconstructionError ;
}

template <typename T>
float 
PcaReprocReconsErrorElement<T>::unreconstructedInertia()
{
  return _unreconstructedInertia ;
}

template <typename T>
float
PcaReprocReconsErrorElement<T>::distanceTo( Element<T> * other,
                                            const std::vector< std::list< Point3d > >& classes,
                                            const AimsData<T>& data, const AimsData<float>& initDistances )
{
  PcaReprocReconsErrorElement<T> * otherPcaReproc = dynamic_cast< PcaReprocReconsErrorElement<T> *>(other) ;
  if( this->constitutedFrom().size() == 1 && other->constitutedFrom().size() == 1 ){
    AimsData<float> indiv( 1, data.dimT() ), indivTr(data.dimT()) ;
    std::list<Point3d>::const_iterator iter( classes[*(this->constitutedFrom().begin())].begin() ), 
      last( classes[*(this->constitutedFrom().begin())].end() ) ;
    AimsData<float> secondErrorMatrix = otherPcaReproc->errorMatrix() ;
    float oneToTwoMeanError = 0. ;
    float twoToOneMeanError = 0. ;
    while( iter != last ){
      float norm2 = 0. ;
      for ( int t = 0 ; t < data.dimT() ; ++t ){
        indiv(0,t) = data( (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        indivTr(t,0) = data(  (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        
        norm2 += data( (*iter)[0], (*iter)[1], (*iter)[2], t )*
          data( (*iter)[0], (*iter)[1], (*iter)[2], t ) ;
      }
      oneToTwoMeanError +=  indiv.cross( secondErrorMatrix.cross(indivTr) )(0,0) / norm2 ;
      ++iter ;
    }
    
    oneToTwoMeanError /= float(classes[*(this->constitutedFrom().begin())].size()) ;
    oneToTwoMeanError = /* 1. - exp( - */( oneToTwoMeanError - otherPcaReproc->meanReconstructionError() ) *
                                 ( oneToTwoMeanError - otherPcaReproc->meanReconstructionError() ) / 
                                 ( 2 * otherPcaReproc->varReconstructionError() /* ) */ ) ;
    
    iter = classes[ *(otherPcaReproc->constitutedFrom().begin())].begin() ;
    last = classes[ *(otherPcaReproc->constitutedFrom().begin())].end() ;
    
    while( iter != last ){
      float norm2 = 0. ;
      for ( int t = 0 ; t < data.dimT() ; ++t ){
        indiv(0,t) = data( (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        indivTr(t,0) = data(  (*iter)[0], (*iter)[1], (*iter)[2], t ) - _mean[t] ;
        
        norm2 += data( (*iter)[0], (*iter)[1], (*iter)[2], t )*
          data( (*iter)[0], (*iter)[1], (*iter)[2], t ) ;
      }
      twoToOneMeanError += indiv.cross( _errorMatrix.cross(indivTr) )(0,0) / norm2 ;
      ++iter ;
    }
    twoToOneMeanError /= float(classes[ *(otherPcaReproc->constitutedFrom().begin())].size()) ;
    twoToOneMeanError = /* 1. - exp (  -*/( twoToOneMeanError - _meanReconstructionError ) *
                                   ( twoToOneMeanError - _meanReconstructionError ) /
                                   ( 2. * _varReconstructionError ) /* ) */ ;
    
    return std::max( oneToTwoMeanError, twoToOneMeanError ) ;
  }
  
  std::set<int> constFrom = otherPcaReproc->constitutedFrom() ;
  std::set<int>::const_iterator cIter( this->constitutedFrom().begin() ), cLast( this->constitutedFrom().end() ) ;
  while( cIter != cLast ){
    constFrom.insert( *cIter ) ;
    ++cIter ;
  }

  std::set<int>::const_iterator c1Iter( constFrom.begin() ), c1Last( constFrom.end() ) ;
  
  float maxDistance = 0. ;
  while( c1Iter != c1Last )
    {
      std::set<int>::const_iterator c2Iter( c1Iter ) ;
      ++c2Iter ;
      while( c2Iter != c1Last )
        {
          float initDist = initDistances( *c1Iter, *c2Iter ) ;

          if( maxDistance < initDist )
            maxDistance = initDist ;
          ++c2Iter ;
        }
      ++c1Iter ;
    }
  
  return maxDistance ;
}

template <typename T>
Element<T> * 
PcaReprocReconsErrorElement<T>::agregateTo( Element<T> * other, 
                           const std::vector< std::list< Point3d > >& classes,
                           const AimsData<T>& data, const AimsData<float>&, int newNb  )
{
  // Agglomeration
  std::set<int> constFrom( this->constitutedFrom() ) ;
  std::set<int>::const_iterator cIt(other->constitutedFrom().begin() ), 
    cLast(other->constitutedFrom().end() ) ;
  while( cIt != cLast ){
    constFrom.insert( *cIt ) ;
    ++cIt ;
  }
  
  Element<T> * newEl = new PcaReprocReconsErrorElement<T>( classes, data, _numberOfSignificantVps, constFrom, newNb ) ;
  return newEl ;
}

template <typename T>
CahElementFactory<T>::CahElementFactory( const std::vector< std::list< Point3d > >& classes,
                                         const AimsData<T>& data ) : _classes( classes ), _data( data ) {}

template <typename T>
CahElementFactory<T>::~CahElementFactory( ) {}


template <typename T>
std::list< Element<T>* > 
CahElementFactory<T>::creator() const 
{
  std::list<Element<T>*> listeEl ;
  for( unsigned int c = 0 ; c < _classes.size() ; ++c )
    listeEl.push_back( creator(c) ) ;
  return listeEl ;
}

template <typename T>
BaryCahElementFactory<T>::BaryCahElementFactory( const std::vector< std::list< Point3d > >& classes,
                       const AimsData<T>& data ) : CahElementFactory<T>(classes, data) {}

template <typename T>
BaryCahElementFactory<T>::~BaryCahElementFactory( ) {}

template <typename T>
Element<T> * 
BaryCahElementFactory<T>::creator( int n )  const
{
  return new BaryElement<T>( this->_classes, this->_data, std::set<int>(), n ) ;
}


template <typename T>
PcaCahElementFactory<T>::PcaCahElementFactory( int nbOfSignificantVps, 
                                               const std::vector< std::list< Point3d > >& classes,
                                               const AimsData<T>& data ) : 
  CahElementFactory<T>(classes, data), _nbOfSignificantVps(nbOfSignificantVps) {}

template <typename T>
PcaCahElementFactory<T>::~PcaCahElementFactory( ) {}

template <typename T>
Element<T> * 
PcaCahElementFactory<T>::creator( int n )  const
{
  return new PcaElement<T>( this->_classes, this->_data, _nbOfSignificantVps, std::set<int>(), n ) ;
}



template <typename T>
Max2By2DistanceCahElementFactory<T>::Max2By2DistanceCahElementFactory( const AimsData<float>& initDistances ) : 
  CahElementFactory<T>( ), _initDistances( initDistances.clone() ), _nbClasses(initDistances.dimX())
{
  std::cout << "Init distances = " << _initDistances.dimX() << ", " << _initDistances.dimY() << std::endl ;
}
 
template <typename T>
Max2By2DistanceCahElementFactory<T>::~Max2By2DistanceCahElementFactory( ) {}

template <typename T>
Element<T> * 
Max2By2DistanceCahElementFactory<T>::creator( int n )  const
{
  //std::cout << "Creator " << n << std::endl ;
  return new Max2By2DistanceElement<T>( std::set<int>(), n ) ;
}

template <typename T>
std::list< Element<T>* > 
Max2By2DistanceCahElementFactory<T>::creator() const 
{
  std::list<Element<T>*> listeEl ;
  std::cout << "Nb classes = " << _nbClasses << std::endl ;
  for( int c = 0 ; c < _nbClasses ; ++c )
    listeEl.push_back( creator(c) ) ;
  return listeEl ;
}

AimsData<double> * 
kineticSampleDistance( const std::vector<aims::DiscriminantAnalysisElement*>& classes )
{
  std::cout << "CAH::kineticSampleDistance" << std::endl ;
  for( unsigned int c = 0 ; c < classes.size() ; ++c )
    if( ! classes[c]->computed() )
      std::cerr << "Not computed for class " << c << std::endl ;
  int c1, c2 ;
  double nbFrames = double( classes[0]->mean().dimX() ) ;
  AimsData<double> distance2by2( classes.size(), classes.size() ) ;
  //std::cout << "Centroids distance : " << std::endl ;
  ForEach2d( distance2by2, c1, c2 ){
    distance2by2( c1, c2 ) = classes[c1]->distance( classes[c2]->mean() ) ;
    //std::cout  << "c1 = " << c1 << "  c2 = " << c2 << " == " <<distance2by2( c1, c2 )  << std::endl ;
  }
  AimsData<double> * symDistance2by2 = new AimsData<double>( classes.size(), classes.size() ) ;
  ForEach2d( distance2by2, c1, c2 ){
    (*symDistance2by2)( c1, c2 ) = 0.5 * ( distance2by2( c1, c2 ) + distance2by2( c2, c1 ) ) / nbFrames ;
  }
//   aims::Writer< AimsData<double> > wri("symDistance2by2.ima") ;
//   wri.write( *symDistance2by2 ) ;

  return symDistance2by2 ;
}

AimsData<double>
interactionSurfaces( const AimsData<int>& locMap, const std::vector<short>& unfusionedLabels,
                     AimsData<double>& interactionSurf, std::vector<double>& labelsVolumes )
{
  std::cout << "CAH::interactionSurfaces" << std::endl ;
  int x, y, z ;
  int maxLabel = unfusionedLabels[0] ;
  for(unsigned int ind = 0 ; ind < unfusionedLabels.size() ; ++ind ){
    if( unfusionedLabels[ind] > maxLabel ){
      maxLabel = unfusionedLabels[ind] ;
    }
    //std::cout << "  " << unfusionedLabels[ind] ;
  }
  
  std::cout << std::endl << "max label computed : " << maxLabel << " on " << unfusionedLabels.size() << "indivs" <<  std::endl ;

  labelsVolumes = std::vector<double>( maxLabel + 1 ) ;
  interactionSurf = AimsData<double>( maxLabel + 1, maxLabel + 1 ) ;

  aims::Connectivity  * conn =  0 ; 
  if( locMap.dimZ() == 1 )
    conn = new aims::Connectivity( 0, 0, aims::Connectivity::CONNECTIVITY_4_XY ) ;
  else
    conn = new aims::Connectivity( 0, 0, aims::Connectivity::CONNECTIVITY_6_XYZ ) ;
  
  Point3df voxelSize( locMap.sizeX(), locMap.sizeY(), locMap.sizeZ() ) ;
  float voxelVolume = locMap.sizeX() * locMap.sizeY() * locMap.sizeZ() ;
  std::vector<float> surfaces( conn->nbNeighbors(), voxelVolume ) ;
  
  for( int n = 0 ; n < conn->nbNeighbors() ; ++n ){
    short index = 0 ;
    if( conn->xyzOffset(n)[1] != 0. )
      index = 1 ;
    else if( conn->xyzOffset(n)[2] != 0. ) 
      index = 2 ;
    surfaces[n] /= voxelSize[index] ;
    std::cout << "xyzOffset : " << conn->xyzOffset(n) << " : " << surfaces[n] << std::endl ;
  }
  
  std::cout << "surface vector computed" << std::endl ;
  
  int ind, indLabel, neigh, neighLabel ;
  Point3d p ;
  Point3d pNeigh ;
  Point3d dims( locMap.dimX(), locMap.dimY(), locMap.dimZ() ) ;
  ForEach3d( locMap, x, y, z ){
    ind = locMap( x, y, z ) ;
    indLabel = unfusionedLabels[ind] ;
    labelsVolumes[ indLabel ] += voxelVolume ;
    p = Point3d(x, y, z) ;
    for( int n = 0 ; n < conn->nbNeighbors() ; ++n ){
      pNeigh = p + conn->xyzOffset(n) ;
      if( pNeigh[0] >= 0 && pNeigh[1] >= 0 && pNeigh[2] >= 0 && 
          pNeigh[0] < locMap.dimX() && pNeigh[1] < locMap.dimY() && pNeigh[2] < locMap.dimZ() ) {
        neigh = locMap( pNeigh ) ;
        neighLabel = unfusionedLabels[neigh] ;
        if( indLabel != neighLabel ){
          interactionSurf( indLabel, neighLabel ) += surfaces[n] ;
//        std::cout << "p : " << p << " l = " << indLabel << "\tNeigh : " << pNeigh << " l == " << neighLabel 
//                  << "  IS(p, n) = " << interactionSurf( indLabel, neighLabel )<< std::endl ;
        }
      }
    }
  }
  std::cout << "interaction surfaces computed" << std::endl ;

  int c1, c2 ;
  AimsData<double> surfaceVolumeRatio( maxLabel + 1, maxLabel + 1 ) ;
  ForEach2d( interactionSurf, c1, c2 ){
    if( c1 != c2 )
      surfaceVolumeRatio( c1, c2 ) = interactionSurf( c1, c2 ) / labelsVolumes[c1] ;
  }

  AimsData<double> symetrizedSurfaceVolumeRatio( maxLabel + 1, maxLabel + 1 ) ;
  double maximum, minimum = std::max( surfaceVolumeRatio( 0, 1 ), surfaceVolumeRatio( 1, 0 ) ) ;
  maximum = minimum ;

  ForEach2d( interactionSurf, c1, c2 ){
    if( c1 != c2 ){
      symetrizedSurfaceVolumeRatio( c1, c2 ) = std::max( surfaceVolumeRatio( c1, c2 ), surfaceVolumeRatio( c2, c1 ) ) ;
      if( symetrizedSurfaceVolumeRatio( c1, c2 ) > maximum )
        maximum = symetrizedSurfaceVolumeRatio( c1, c2 ) ;
      if( symetrizedSurfaceVolumeRatio( c1, c2 ) < minimum )
        minimum = symetrizedSurfaceVolumeRatio( c1, c2 ) ;
    }
  }
  
//   aims::Writer< AimsData<double> > wri1("interactionSurfaces.ima") ;
//   wri1.write(interactionSurf) ;
//   aims::Writer< AimsData<double> > wri2("surfaceVolumeRatio.ima") ;
//   wri2.write(surfaceVolumeRatio) ;
//   aims::Writer< AimsData<double> > wri3("symetrizedSurfaceVolumeRatio.ima") ;
//   wri3.write(symetrizedSurfaceVolumeRatio) ;
  
  AimsData<double> distanceMatrix(maxLabel + 1, maxLabel + 1) ;
  ForEach2d( interactionSurf, c1, c2 ){
    if( c1 == c2 )
      distanceMatrix(c1, c2) = 0. ;
    else
      distanceMatrix(c1, c2) = (symetrizedSurfaceVolumeRatio(c1, c2) - minimum) / ( maximum - minimum ) ;
  }
  
  return distanceMatrix ;
}

template <class T>
bool
correlated( const AimsData<T>& x, const AimsData<T>& y  )
{
  if( x.dimX() == 1 && y.dimX() == 1 )
    return true ;
  
  if( x.dimX() != y.dimX() ){
    std::cerr <<"Wrong size for correlation !" << std::endl ;
    throw std::runtime_error("Wrong size for correlation !") ;
  }
  
  T xMean = 0. ;
  T yMean = 0. ;
  
  for( int t = 0 ; t < x.dimX() ; ++t ){
    xMean += x(t) ;
    yMean += y(t) ;
  }
  xMean /= T(x.dimX()) ;
  yMean /= T(y.dimX()) ;
  
  T corr = 0. ;
  for( int t = 0 ; t <  x.dimX() ; ++t )
    corr += ( x(t)  - xMean ) * ( y(t) - yMean ) ;
  return (corr > 0.) ;
}

// AimsData<double>
// momentOfMaximumDifference( const std::vector<float>& momentOfMax )
// {
//   AimsData<double> res( classesKinDescription.size(), classesKinDescription.size() ) ;
//   int c1, c2 ;
//   ForEach2d( res, c1, c2){
//     if( momentOfMax[c1] == momentOfMax[c2] )
//       res( c1, c2) == 0. ; 
//     else{
//       double minMoment = std::min(momentOfMax[c1], momentOfMax[c2]) ;
//       res( c1, c2) = 2. * std::max( momentOfMax[c1] - momentOfMax[c2], momentOfMax[c2] - momentOfMax[c1]) * 
//      exp( - (minMoment * minMoment) 
//           / ( 2. * 0.333 * classesKinDescription[0]->mean().dimX() * 0.333 * classesKinDescription[0]->mean().dimX() ) ) ;
//     }
//   }
  
//   aims::Writer< AimsData<double> > wri( "momentOfMaximumRelDiff.ima" ) ;
//   wri.write(res) ;
  
//   return res ;
// }

AimsData<double>
momentOfMaximumDifference( const std::vector<aims::DiscriminantAnalysisElement*>& classesKinDescription )
{
  std::vector<short> momentOfMax( classesKinDescription.size() ) ;
  std::vector<float> valueOfMax( classesKinDescription.size() ) ;
  int maxInd = 0 ;
  for( unsigned int c = 0 ; c < classesKinDescription.size() ; ++c ){
    for( int t = 0 ; t < classesKinDescription[c]->mean().dimX() ; ++t )
      if( classesKinDescription[c]->mean()(t) > classesKinDescription[c]->mean()(maxInd) )
        maxInd = t ;
    momentOfMax[c] = maxInd ;
  }


  AimsData<double> res( classesKinDescription.size(), classesKinDescription.size() ) ;
  int c1, c2 ;
  ForEach2d( res, c1, c2){
    if( momentOfMax[c1] == momentOfMax[c2] )
      res( c1, c2) == 0. ; 
    else{
      double minMoment = std::min(momentOfMax[c1], momentOfMax[c2]) ;
      res( c1, c2) = 2. * std::max( momentOfMax[c1] - momentOfMax[c2], momentOfMax[c2] - momentOfMax[c1]) * 
        exp( - (minMoment * minMoment) 
             / ( 2. * 0.5 * classesKinDescription[0]->mean().dimX() * 0.5 * classesKinDescription[0]->mean().dimX() ) ) ;
    }
  }
  
//   aims::Writer< AimsData<double> > wri( "momentOfMaximumRelDiff.ima" ) ;
//   wri.write(res) ;
  
  return res ;
}

AimsData<double>
momentOfMaxValueRatios( const std::vector<aims::DiscriminantAnalysisElement*>& classesKinDescription )
{
  // Getting moment of max
  std::vector<short> momentOfMax( classesKinDescription.size() ) ;
  std::vector<float> valueOfMax( classesKinDescription.size() ) ;
  int maxInd = 0 ;
  for( unsigned int c = 0 ; c < classesKinDescription.size() ; ++c ){
    for( int t = 0 ; t < classesKinDescription[c]->mean().dimX() ; ++t )
      if( classesKinDescription[c]->mean()(t) > classesKinDescription[c]->mean()(maxInd) )
        maxInd = t ;
    momentOfMax[c] = maxInd ;
  }
  AimsData<double> valueAtMaxMoment( classesKinDescription.size(), classesKinDescription.size() ) ;
  int c1, c2 ;
  ForEach2d(valueAtMaxMoment, c1, c2){
    //valueAtMaxMoment(c1, c2) = ( classesKinDescription[c1]->mean()( momentOfMax[c2] ) / classesKinDescription[c1]->mean()( momentOfMax[c1] ) ) ;
    valueAtMaxMoment(c1, c2) = 1.  - ( classesKinDescription[c2]->mean()( momentOfMax[c1] ) / classesKinDescription[c2]->mean()( momentOfMax[c2] ) ) ;
  }
  AimsData<double> res( valueAtMaxMoment.clone() );
  ForEach2d(res, c1, c2){
    res(c1, c2) = 0.5 * ( valueAtMaxMoment(c1, c2) + valueAtMaxMoment(c2, c1) ) ;
  }
  
//   aims::Writer< AimsData<double> > wri( "momentOfMaxValueRatios.ima" ) ;
//   wri.write(res) ;
  
  return res ;
}

template<typename T>
SCKBDistanceCahElementFactory<T>::SCKBDistanceCahElementFactory( const std::vector<std::list<Point3d> >& classes, 
                                                                 const std::vector<aims::DiscriminantAnalysisElement*>& classesKinDescription, 
                                                                 const AimsData<int>& locMap,
                                                                 const std::vector<short>& unfusionedLabels,
                                                                 const std::vector<float>& classMaximaMoment,
                                                                 float spatialregularizationWeight, 
                                                                 float timeOfMaximumWeight ) :
  CahElementFactory<T>( classes ), /*_classMaximumMoment( new std::vector<float>( classMaximaMoment ) ), */
  _spatialregularizationWeight(spatialregularizationWeight), 
  _timeOfMaximumWeight(timeOfMaximumWeight)
{
  _interactionSurf = new AimsData<double> ;
  _labelsVolumes = new std::vector<double> ;
  interactionSurfaces( locMap, unfusionedLabels, *_interactionSurf, *_labelsVolumes ) ;
  _maxMomentValueRatios = new AimsData<double>( momentOfMaximumDifference( classesKinDescription /*classMaximaMoment*/ ) ) ;
  _antiCorrelatedMeans = new AimsData<char>( classesKinDescription.size(), classesKinDescription.size() ) ;
  
  int c1, c2 ;
  double nbFrames = double( classesKinDescription[0]->mean().dimX() ) ;
  _kinModel2by2Distances = new AimsData<double>( classesKinDescription.size(), classesKinDescription.size() ) ;
//   std::cout << "Centroids distance : " << std::endl ;
  AimsData<double> distance2by2( classesKinDescription.size(), classesKinDescription.size() ) ;
  double max2By2Distance = 0 ;

  ForEach2d( distance2by2, c1, c2 ){
    distance2by2( c1, c2 ) = classesKinDescription[c1]->distance( classesKinDescription[c2]->mean() ) ;
    if( distance2by2( c1, c2 ) > max2By2Distance )
      max2By2Distance = distance2by2( c1, c2 ) ;
    (*_antiCorrelatedMeans)(c1, c2) = !correlated( classesKinDescription[c1]->mean(), classesKinDescription[c2]->mean()) ;
//     std::cout  << "c1 = " << c1 << "  c2 = " << c2 << " == " <<distance2by2( c1, c2 )  << std::endl ;
  }
  _meanDist2By2 = _varDist2By2 = _meanSurfVolRatio = _varSurfVolRatio =  _meanMaxMomentRatio = _varMaxMomentRatio = 0. ;
  int nbOfCases = 0 ;
  ForEach2d( distance2by2, c1, c2 ){
    (*_kinModel2by2Distances)( c1, c2 ) = 0.5 * ( distance2by2( c1, c2 ) + distance2by2( c2, c1 ) ) / nbFrames ;
    if( c1 > c2 )
      if( !(*_antiCorrelatedMeans)(c1, c2) ){
        ++nbOfCases ;
        float val = (*_kinModel2by2Distances)( c1, c2 ) ;
        _meanDist2By2 += val ;
        _varDist2By2 += val * val ;

        if( (*_labelsVolumes)[c1] == 0. || (*_labelsVolumes)[c2] == 0. )
          val = 0. ;
        else
          val = std::max( (*_interactionSurf)(c1, c2) / (*_labelsVolumes)[c1], 
                          (*_interactionSurf)(c2, c1) / (*_labelsVolumes)[c2] ) ;
        
        _meanSurfVolRatio += val ;
        _varSurfVolRatio += val * val ;
        
        val = 0.5 * (*_maxMomentValueRatios)(c1, c2) ;
        _meanMaxMomentRatio += val ;
        _varMaxMomentRatio += val * val ;
        
      }
  }
  _meanDist2By2 /= float( nbOfCases ) ;
  _meanSurfVolRatio /= float( nbOfCases ) ;
  _meanMaxMomentRatio /= float( nbOfCases ) ;
  
  _varDist2By2 = ( _varDist2By2 / float( nbOfCases ) - _meanDist2By2 * _meanDist2By2 ) 
    / float( nbOfCases - 1 ) * float( nbOfCases ) ;
  _varSurfVolRatio = ( _varSurfVolRatio / float( nbOfCases ) - _meanSurfVolRatio * _meanSurfVolRatio ) 
    / float( nbOfCases - 1 ) * float( nbOfCases ) ;
  _varMaxMomentRatio = ( _varMaxMomentRatio / float( nbOfCases ) - _meanMaxMomentRatio * _meanMaxMomentRatio ) 
    / float( nbOfCases - 1 ) * float( nbOfCases ) ;
  if( _varMaxMomentRatio == 0. )
    _varMaxMomentRatio = 1. ;
    
  //-------------------------------------
  double meanCentered1 = 0., meanCentered2 = 0. ; 

  ForEach2d( distance2by2, c1, c2 ){
    (*_kinModel2by2Distances)( c1, c2 ) = 0.5 * ( distance2by2( c1, c2 ) + distance2by2( c2, c1 ) ) / nbFrames ;
    if( c1 > c2 )
      if( !(*_antiCorrelatedMeans)(c1, c2) ){
        //std::cout << c1 << " && " << c2 << " : " << (*_kinModel2by2Distances)( c1, c2 ) - _meanDist2By2 << std::endl ;
        meanCentered1 += (*_kinModel2by2Distances)( c1, c2 ) - _meanDist2By2 ;
        float val ;
        if( (*_labelsVolumes)[c1] == 0. || (*_labelsVolumes)[c2] == 0. )
          val = 0. ;
        else
          val = std::max( (*_interactionSurf)(c1, c2) / (*_labelsVolumes)[c1], 
                          (*_interactionSurf)(c2, c1) / (*_labelsVolumes)[c2] ) ;
        
        meanCentered2 += val - _meanSurfVolRatio ;
      }
  }
  //-------------------------------------

  AimsData<double> correctedDists( classesKinDescription.size(), classesKinDescription.size() ) ;
  AimsData<double> correctedSurfVolRatio( correctedDists.clone() ), correctedMomentOfMax( correctedDists.clone() ) ;
  ForEach2d( correctedDists, c1, c2 ){
    correctedDists(c1, c2) = ( (*_kinModel2by2Distances)( c1, c2 ) - _meanDist2By2 ) / sqrt( _varDist2By2 ) ;
    correctedSurfVolRatio(c1, c2) = -( std::max( (*_interactionSurf)(c1, c2) / (*_labelsVolumes)[c1], 
                                                (*_interactionSurf)(c2, c1) / (*_labelsVolumes)[c1] ) - _meanSurfVolRatio ) / sqrt( _varSurfVolRatio ) ;
//     correctedMomentOfMax(c1, c2) = ( 1. - std::min( (*_classMaximumMoment)[c1] / (*_classMaximumMoment)[c2], 
//                                                  (*_classMaximumMoment)[c2] / (*_classMaximumMoment)[c1] ) - _meanMaxMomentRatio )
//       / sqrt( _varMaxMomentRatio ) ;
    correctedMomentOfMax(c1, c2) = ( (*_maxMomentValueRatios)(c1, c2) - _meanMaxMomentRatio )
      / sqrt( _varMaxMomentRatio ) ;
  }
  
//   aims::Writer< AimsData<double> > wri1("correctedDists.ima") ;
//   wri1.write(correctedDists) ;
  
//   aims::Writer< AimsData<double> > wri2("correctedSurfVolRatio.ima") ;
//   wri2.write(correctedSurfVolRatio) ;
  
  
//   aims::Writer< AimsData<double> > wri3("correctedMomentOfMax.ima") ;
//   wri3.write(correctedMomentOfMax) ;
  
//   aims::Writer< AimsData<char> > wri4("antiCorr.ima") ;
//   wri4.write(*_antiCorrelatedMeans) ;
  
//   aims::Writer< AimsData<double> > wri("symDistance2by2.ima") ;
//   wri.write( *_kinModel2by2Distances ) ;
  
//   std::cout << "SCKBDistanceCahElementFactory<T>::SCKBDistanceCahElementFactory end" << std::endl ;
}


template <typename T>
Element<T>*
SCKBDistanceCahElementFactory<T>::creator(int n ) const
{
  return new SCKBDistanceElement<T>( _maxMomentValueRatios, _kinModel2by2Distances, _interactionSurf, _labelsVolumes,
                                     _spatialregularizationWeight, _timeOfMaximumWeight, 
                                     _meanDist2By2, _varDist2By2, _meanSurfVolRatio, _varSurfVolRatio, _meanMaxMomentRatio, _varMaxMomentRatio,
                                     std::set<int>(), n ) ;
}

template <typename T>
PPcaCahElementFactory<T>::PPcaCahElementFactory( int nbOfSignificantVps, 
                                               const std::vector< std::list< Point3d > >& classes,
                                               const AimsData<T>& data ) : 
  CahElementFactory<T>(classes, data), _nbOfSignificantVps(nbOfSignificantVps) {}

template <typename T>
PPcaCahElementFactory<T>::~PPcaCahElementFactory( ) {}

template <typename T>
Element<T> * 
PPcaCahElementFactory<T>::creator( int n )  const
{
  return new PPcaElement<T>( this->_classes, this->_data, _nbOfSignificantVps, std::set<int>(), n ) ;
}

template <typename T>
PcaReprocReconsErrorCahElementFactory<T>::PcaReprocReconsErrorCahElementFactory( int nbOfSignificantVps, 
                                               const std::vector< std::list< Point3d > >& classes,
                                               const AimsData<T>& data ) : 
  CahElementFactory<T>(classes, data), _nbOfSignificantVps(nbOfSignificantVps) {}

template <typename T>
PcaReprocReconsErrorCahElementFactory<T>::~PcaReprocReconsErrorCahElementFactory( ) {}

template <typename T>
Element<T> * 
PcaReprocReconsErrorCahElementFactory<T>::creator( int n )  const
{
  return new PcaReprocReconsErrorElement<T>( this->_classes, this->_data, _nbOfSignificantVps, std::set<int>(), n ) ;
}

template <typename T>
Cah<T>::Cah( const CahElementFactory<T> * method ) :
  _classes( method->classes() ), _data(method->data()), _m(0), _canAddElement(true), 
  _initDistancesAlreadyComputed(false)
{
  if( method->classes().size() == 0 && method->methodType() == "Pca" )
    throw std::runtime_error("You must specify a method for ascendant hierarcical classification !") ;

  if( method->methodType() == "Max2By2Distance" )
    {
      const Max2By2DistanceCahElementFactory<T>* fact = 
        dynamic_cast< const Max2By2DistanceCahElementFactory<T>* >( method ) ;
      
      if( fact->initDistances()->dimX() != fact->initDistances()->dimY() )
        throw std::runtime_error("Incompatible pre-computed distances between element matrix size !") ;
      
      _distanceBetweenInitialElements = fact->initDistances()->clone() ;
      
      _m = fact->initDistances()->dimX() ;
      _listeN = fact->creator() ;
      _initDistancesAlreadyComputed = true ;
      _canAddElement = false ;
    } 
  else if( method->classes().size() )
    {
      _m = method->classes().size() ;    
      _listeN = method->creator() ;
      _canAddElement = false ;
    }
}

template <typename T>
bool Cah<T>::addElement(Element<T> * el )
{
  if( !_canAddElement )
    return false ;
  if( el->number() == -1 )
    el->setNumber( _m + 1 ) ;
  _listeN.push_back(el);
  _m++;
  
  return true ;
}


template <typename T>
void Cah<T>::initialisationDist()
{
  if( _initDistancesAlreadyComputed ){
    for( int i = 0 ; i < _distanceBetweenInitialElements.dimX() ; ++i )
      for( int j = i + 1 ; j < _distanceBetweenInitialElements.dimY() ; ++j )
        _mMap.insert( std::pair<float, Point2d>( _distanceBetweenInitialElements(i, j), 
                                            Point2d(i, j) ) ) ;
//     aims::Writer< AimsData<float> > wt("Distances.ima") ;
//     wt.write( _distanceBetweenInitialElements ) ;
  
    return ;
  }
  
  _distanceBetweenInitialElements = AimsData<float>(_m, _m) ;
    
  float d = -1.;
  typename std::list<Element<T>*>::iterator 
    iiter(_listeN.begin()), ilast(_listeN.end()) ;
  
  int i = 0 ;

  while( iiter != ilast )
    {
      typename std::list<Element<T>*>::iterator 
        jiter(iiter), jlast(_listeN.end()) ;
      ++jiter ; 
      int j = i + 1 ;
      while( jiter != jlast )
        {
          if( (*iiter) != (*jiter) ){
            d = (*iiter)->distanceTo( *jiter, _classes, _data, 
                                      _distanceBetweenInitialElements );
            _mMap.insert(std::pair<float, Point2d>
                         (d, Point2d((*iiter)->number(),(*jiter)->number())));
            _distanceBetweenInitialElements( (*iiter)->number(), 
                                             (*jiter)->number() ) = d ;
            _distanceBetweenInitialElements( (*jiter)->number(), 
                                             (*iiter)->number() ) = d ;
          }
          ++jiter ; ++j ;
        }
      ++iiter ; ++i ;
    }
//   aims::Writer< AimsData<float> > wt("Distances.ima") ;
//   wt.write( _distanceBetweenInitialElements ) ;
}

template <typename T>
void Cah<T>::oneStep(int etape)
{
  std::multimap<float, Point2d>::const_iterator it = _mMap.begin();
  int Na = it->second[0], Nb = it->second[1];
  Element<T> * elA = 0, *elB = 0 ;
/*   int card_a, card_b; */
  
  // On rajoute la fusion dans _result
  _result.push_back(std::pair<Point2d, float> ( Point2d(Na, Nb),it->first));
  std::cout <<Na<<"  "<<Nb<<"  "<<(it->first) << " : " << std::endl;
    
  // Recherche de Na et Nb dans _listeN, sauvegarde et suppression...
  typename std::list<Element<T>*>::iterator vectiter, tmpiter;
  vectiter=_listeN.begin();
  while(vectiter!=_listeN.end())
    {
      if( ((*vectiter)->number()) == Na )
        {
          elA = (*vectiter) ;
          tmpiter = vectiter;
          ++vectiter ;
          _listeN.erase(tmpiter);
        }
      else 
        {
          if( ((*vectiter)->number()) == Nb )
            {
              elB = (*vectiter) ;
              tmpiter = vectiter;
              ++vectiter ;
              _listeN.erase(tmpiter);
            }
          else 
            {
              ++vectiter;
            }
        }
    }
  
  Element<T> *nouveau = elA->agregateTo( elB, _classes, _data, 
                                         _distanceBetweenInitialElements, _m + etape ) ;
  
  std::set<int>::const_iterator cit(nouveau->constitutedFrom().begin()), clt(nouveau->constitutedFrom().end()) ;
  while( cit != clt ){
    ++cit ;
  } 
  // Suppression dans _mMap de toutes les distances faisant reference a Na ou Nb        
  std::multimap<float, Point2d>::iterator mapiter, tempiter;
  mapiter= _mMap.begin();
  while(mapiter!=_mMap.end())
    {
      if( (mapiter->second[0]==Na) || (mapiter->second[1]==Na) ||
          (mapiter->second[0]==Nb) || (mapiter->second[1]==Nb) )
        {
          tempiter = mapiter;
          ++mapiter;
          _mMap.erase(tempiter);
        }
      else ++mapiter;
    } 
  
  // Calcul des distances concernant le nouvel element dans _mMap
  float d = -1.;
  typename std::list<Element<T>*>::iterator 
    iiter(_listeN.begin()), ilast(_listeN.end()) ;
  while( iiter != ilast )
    {
      d = nouveau->distanceTo( *iiter, _classes, _data, _distanceBetweenInitialElements ) ;
      _mMap.insert(std::pair<float, Point2d>
                   (d, Point2d((*iiter)->number(), nouveau->number())) );
      ++iiter ;
    }
  
  
/*   delete elA, elB ; */
  // On rajoute a present le nouvel element dans _listeN
  _listeN.push_back(nouveau);
}


template <typename T>
const std::vector< std::pair<Point2d, float> >& Cah<T>::doit()
{
  initialisationDist();
  for (unsigned int i=1; i<_m;i++)
    {
      std::cout<<"Etape n°"<<i<<std::endl;
      oneStep(i);
    }
  return _result;
}




template <typename T>
std::vector<int> Cah<T>::leaves(int nc)
{
//      std::cerr<<"On entre dans leaves classes..."<<std::endl;
  std::vector<int> liste;
  for(unsigned int i=_m-nc;i<_m-1;i++)
    {
      if((unsigned int)_result[i].first[0]<(2*_m-nc+1))
        {
          liste.push_back(_result[i].first[0]);
        }
      
      if((unsigned int)_result[i].first[1]<(2*_m-nc+1))
        {
          liste.push_back(_result[i].first[1]);
        }
    }
  
  return liste;
}


template <typename T>
std::vector<int> Cah<T>::leaves(float d)
{
//      std::cerr<<"On entre dans leaves distance..."<<std::endl;
  std::vector<int> liste;
  for(int i=_m-2;(_result[i].second<d)||(i<0);i--)
        {
          if((unsigned int)_result[i].first[0]<=_m)
            {
              liste.push_back(_result[i].first[0]);
            }   
          else
            {
              if(_result[_result[i].first[0]-_m-1].second<d)
                {
                  liste.push_back(_result[i].first[0]);
                }
            }
          
          if((unsigned int)_result[i].first[1]<=_m)
            {
              liste.push_back(_result[i].first[1]);
            }   
          else
            {
              if(_result[_result[i].first[1]-_m-1].second<d)
                {
                  liste.push_back(_result[i].first[1]);
                }
            }
        }
  
  return liste;
}


//  Procédure récursive de parcours d'arbre
template <typename T>
void Cah<T>::Proc(int N, std::vector<int>& aClass )
{
  if(N <= int(_m) )
    {
//              std::cerr<<"On rentre N = "<<N<<" car il est plus petit que _m = "<<_m<<std::endl;
      aClass.push_back(N);
    }
  else
    {
      Proc(_result[N-_m-1].first[0], aClass );
      Proc(_result[N-_m-1].first[1], aClass );
    }
}


template <typename T>
std::vector<Point2d> Cah<T>::cut( const std::vector< std::pair<Point2d, float> >& res_cah, int nc )
{
   
//      std::cerr<<"On entre dans cut classes..."<<std::endl;
  if( res_cah.size() ) 
    _result = res_cah ;

  _m = _result.size() + 1;
  //    std::cerr<<"_m vaut : "<<_m<<std::endl;
  
  std::vector<int> liste = leaves(nc);
  
  //    std::cerr<<"leaves est constituée de : ";
//      for(unsigned i=0;i<liste.size();i++)    std::cerr<<liste[i]<<" ";
//      std::cerr<<std::endl;

//      std::vector<std::vector<int> > res;

  std::vector<Point2d> res;
  std::map<int,int> t;
  for(unsigned int c=0;c<liste.size();c++)
    {
      std::vector<int> aClasse ;
//              std::cerr<<"On lance Proc avec "<<liste[c]<<std::endl;
      Proc(liste[c], aClasse);
//              res.push_back(_uneClasse);
      for(unsigned int i=0;i<aClasse.size();i++)
        {
          t[aClasse[i]] = (c+1);
        }
    }

  std::map<int,int>::iterator mapit;
  Point2d point;
  for(mapit=t.begin();mapit!=t.end();mapit++)
    {
      point[0] = mapit->first;
      point[1] = mapit->second;
      res.push_back(point);
    }
  
  return res;
}


template <typename T>
std::vector<Point2d> Cah<T>::distanceCut( const std::vector< std::pair<Point2d, float> >& res_cah, float d )
{
/*   if( res_cah.size() )  */
/*     _result = res_cah ; */
  
  int nc= 0 ;
  while( (unsigned int)(nc + 1) < res_cah.size() && res_cah[nc + 1].second < d  )
    ++nc ;

  return cut( res_cah, _m + 1 - nc ) ;

}

template <typename T>
std::vector<Point2d> Cah<T>::maxCurvatureCut( const std::vector< std::pair<Point2d, float> >& res_cah )
{
/*   if( res_cah.size() )  */
/*     _result = res_cah ; */
/*   _result = res_cah ; */
  
  AimsData<float> derivate( res_cah.size() ) ;
  AimsData<float> cost( res_cah.size() ) ;
  AimsData<float> curvature( res_cah.size() ) ;
  
  for( unsigned int i = 1 ; i < res_cah.size() - 1 ; ++i ){
    cost[i] = res_cah[i].second ;
    derivate[i] = ( res_cah[i+1].second - res_cah[i-1].second) * 0.5 ;
  }
  
  int maxCurvInd = 0 ;
  float maxCurv = 0. ;
  for( unsigned int i = 1 ; i < res_cah.size() - 1 ; ++i ){
    curvature[i] = ( res_cah[i+1].second - 2 * res_cah[i].second + res_cah[i-1].second) * 0.25 ;
    
    if( curvature[i] > maxCurv ){
      maxCurv = curvature[i] ;
      maxCurvInd = i ;
    }
  }
  
//   aims::Writer< AimsData <float> > wri1("cost.ima") ;
//   wri1.write(cost) ;
//   aims::Writer< AimsData <float> > wri2("derivate.ima") ;
//   wri2.write(derivate) ;
//   aims::Writer< AimsData <float> > wri3("curv.ima") ;
//   wri3.write(curvature) ;
  
  return cut( res_cah, _m + 1 - maxCurvInd ) ;
}

template <typename T>
void Cah<T>::save( const std::string& fileout )
{
  std::string car = fileout + "_cah" ;

  std::ofstream os( car.c_str() ) ; 
  
  for(unsigned i=0;i<_result.size();i++)
    {
      os << _result[i].first[0] << " " << _result[i].first[1]
         << " " << _result[i].second << std::endl ;
    }
  
  os.close();

  AimsData<float> derivate( _result.size() ) ;
  AimsData<float> cost( _result.size() ) ;
  AimsData<float> curvature( _result.size() ) ;
  
  for( unsigned int i = 1 ; i < _result.size() - 3 ; ++i ){
    cost[i] = _result[i].second ;
    derivate[i] = ( _result[i+1].second - _result[i-1].second) * 0.5 ;
  }
  
  for( unsigned int i = 1 ; i < _result.size() - 1 ; ++i )
    curvature[i] = ( _result[i+1].second - 2 * _result[i].second + _result[i-1].second) * 0.25 ;
  
//   aims::Writer< AimsData <float> > wri1(fileout + "_cost.ima") ;
//   wri1.write(cost) ;
//   aims::Writer< AimsData <float> > wri2(fileout + "_derivative.ima") ;
//   wri2.write(derivate) ;
//   aims::Writer< AimsData <float> > wri3(fileout + "_curvature.ima") ;
//   wri3.write(curvature) ;
}

template <typename T>
const std::vector< std::pair<Point2d, float> >& Cah<T>::load( const std::string& filein )
{
  // Lecture fichier cah
  _result.clear();
  unsigned int Na = 0, Nb = 0;
  float d = -1.;

  std::ifstream is( filein.c_str() ) ;
  
  is >> Na ;
  is >> Nb ;
  is >> d ;
  _result.push_back(std::pair<Point2d, float> (Point2d(Na, Nb), d) );

  while( !is.eof() ) 
    {
      is >> Na ;
      is >> Nb ;
      is >> d ;
      
      _result.push_back(std::pair<Point2d, float> (Point2d(Na, Nb), d) );
    }

  is.close();
  return _result ;
}