melpca.cc

/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
              independent components
    
    melpca.cc - PCA and whitening 

    Christian F. Beckmann, FMRIB Image Analysis Group
    
    Copyright (C) 1999-2004 University of Oxford */

/*  CCOPYRIGHT  */

#include "newimage/newimageall.h"
#include "utils/log.h"
#include "meloptions.h"
#include "meldata.h"
#include "melodic.h"
#include "newmatap.h"
#include "newmatio.h"
#include "melpca.h"
#include "melhlprfns.h"
#include "libprob.h"

using namespace Utilities;
using namespace NEWIMAGE;

namespace Melodic{

  void MelodicPCA::perf_pca(Matrix& in, Matrix& weights)
   {    
     message("Starting PCA  ... ");

     Matrix Corr;
     Matrix tmpE;
     RowVector tmpD, AdjEV, PercEV;
   
     std_pca(in,weights,Corr,tmpE,tmpD); 
     if(opts.tsmooth.value()){
       message(endl << "  temporal smoothing of Eigenvectors " << endl);
       tmpE=smoothColumns(tmpE);
     }
   
     adj_eigspec(tmpD,AdjEV,PercEV);
     melodat.set_pcaE(tmpE);
     melodat.set_pcaD(tmpD);
     melodat.set_EVP(PercEV); 
     melodat.set_EV((AdjEV));
     write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV);
   
     message("done" << endl);
   }
  
  void MelodicPCA::perf_white()
  {
    int N = melodat.get_pcaE().Ncols();    
    if(opts.pca_dim.value() > N){
      message("dimensionality too large - using -dim " << N << 
              " instead " << endl);
      opts.pca_dim.set_T(N);
    }
    if(opts.pca_dim.value() < 0){
      if(opts.remove_meanvol.value()){
	opts.pca_dim.set_T(N-2);
      }else{
   	opts.pca_dim.set_T(N-1);
      }
    }
    if(opts.pca_dim.value() ==0){
      opts.pca_dim.set_T(pcadim());
      if(melodat.get_Data().Nrows()<20){
	opts.pca_dim.set_T(N-2);
	message("too few data points for full estimation, using "
		<< opts.pca_dim.value() << " instead"<< endl);
      }
    }
    if(opts.approach.value()==string("jade") && opts.pca_dim.value() > 30){
      message("dimensionality too large for jade estimation - using --dim 30 instead" << endl);
      opts.pca_dim.set_T(30);
    }

    message("Start whitening using  "<< opts.pca_dim.value()<<" dimensions ... " << endl);

    Matrix tmpWhite;
    Matrix tmpDeWhite;

    float varp = 1.0;
    varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(), 
    	      melodat.get_EVP() ,opts.pca_dim.value(),
    	      tmpWhite,tmpDeWhite);

    melodat.set_white(tmpWhite);
    melodat.set_dewhite(tmpDeWhite);
    message("  retaining "<< varp*100 <<" percent of the variability " << endl);
    message(" ... done"<< endl << endl);
   }

  int MelodicPCA::pcadim()
  { 
    message("Estimating the number of sources from the data (PPCA) ..." << endl);

    ColumnVector PPCA; RowVector AdjEV, PercEV;   
    int res = ppca_dim(melodat.get_Data(),melodat.get_RXweight(), PPCA, AdjEV, PercEV, melodat.get_resels(), opts.pca_est.value());
     
    write_ascii_matrix(logger.appendDir("PPCA"),PPCA);			      
    write_ascii_matrix(logger.appendDir("eigenvalues_adjusted"),AdjEV.t());
    write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV.t());

    melodat.set_EVP(PercEV);
    melodat.set_EV(AdjEV);
    melodat.set_PPCA(PPCA);
    
    return res;
  }

}