melpca.cc

/*  MELODIC - Multivariate exploratory linear optimized decomposition into
              independent components

    melpca.cc - PCA and whitening

    Christian F. Beckmann, FMRIB Analysis Group

    Copyright (C) 1999-2013 University of Oxford */

/*  CCOPYRIGHT  */

#include "utils/log.h"
#include "meloptions.h"
#include "meldata.h"
#include "melodic.h"
#include "armawrap/newmat.h"
#include "melpca.h"
#include "melhlprfns.h"
#include "cprob/libprob.h"

using namespace Utilities;

namespace Melodic{

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

    SymmetricMatrix Corr;
    Matrix tmpE;
    RowVector tmpD, AdjEV, PercEV;

	if(opts.paradigmfname.value().length()>0)
	{
		basicGLM tmpglm;
		tmpglm.olsfit(in,melodat.get_param(),IdentityMatrix(melodat.get_param().Ncols()));
		std_pca(tmpglm.get_residu(),weights,Corr,tmpE,tmpD);
	}
	else{
 		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;

	if (opts.paradigmfname.value().length()>0)
	    varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
    		melodat.get_EVP(),opts.pca_dim.value(),melodat.get_param(),melodat.get_paramS(),tmpWhite,tmpDeWhite);
	else
		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;
  }

}