This commit was manufactured by cvs2svn to create branch 'newnew'.

Cherrypick from master 2015-09-14 14:53:31 UTC Matthew Webster <mwebster@fmrib.ox.ac.uk> 'fix for tica/migp problem':
    dual_regression
    melgmix.cc
    melgmix.h
    melhlprfns.h
    melica.cc
    melica.h
    melpca.cc
    melpca.h

dual_regression

+#!/bin/sh
+
+#   dual_regression - take group-ICA maps (etc) and get subject-specific versions of them (and associated timecourses)
+#
+#   Stephen Smith and Christian Beckmann, FMRIB Image Analysis Group
+#
+#   Copyright (C) 2011-2012 University of Oxford
+#
+#   SHCOPYRIGHT
+
+
+Usage() {
+    cat <<EOF
+
+dual_regression v0.5 (beta)
+
+***NOTE*** ORDER OF COMMAND-LINE ARGUMENTS IS DIFFERENT FROM PREVIOUS VERSION
+
+Usage: dual_regression <group_IC_maps> <des_norm> <design.mat> <design.con> <n_perm> <output_directory> <input1> <input2> <input3> .........
+e.g.   dual_regression groupICA.gica/groupmelodic.ica/melodic_IC 1 design.mat design.con 500 grot \`cat groupICA.gica/.filelist\`
+
+<group_IC_maps_4D>            4D image containing spatial IC maps (melodic_IC) from the whole-group ICA analysis
+<des_norm>                    0 or 1 (1 is recommended). Whether to variance-normalise the timecourses used as the stage-2 regressors
+<design.mat>                  Design matrix for final cross-subject modelling with randomise
+<design.con>                  Design contrasts for final cross-subject modelling with randomise
+<n_perm>                      Number of permutations for randomise; set to 1 for just raw tstat output, set to 0 to not run randomise at all.
+<output_directory>            This directory will be created to hold all output and logfiles
+<input1> <input2> ...         List all subjects' preprocessed, standard-space 4D datasets
+
+<design.mat> <design.con>     can be replaced with just
+-1                            for group-mean (one-group t-test) modelling.
+If you need to add other randomise options then edit the line after "EDIT HERE" in the dual_regression script
+
+EOF
+    exit 1
+}
+
+############################################################################
+
+[ "$6" = "" ] && Usage
+
+ORIG_COMMAND=$*
+
+ICA_MAPS=`${FSLDIR}/bin/remove_ext $1` ; shift
+
+DES_NORM=--des_norm
+if [ $1 = 0 ] ; then
+  DES_NORM=""
+fi ; shift
+
+if [ $1 = "-1" ] ; then
+  DESIGN="-1"
+  shift
+else
+  dm=$1
+  dc=$2
+  DESIGN="-d $1 -t $2"
+  shift 2
+fi
+
+NPERM=$1 ; shift
+
+OUTPUT=`${FSLDIR}/bin/remove_ext $1` ; shift
+
+while [ _$1 != _ ] ; do
+  INPUTS="$INPUTS `${FSLDIR}/bin/remove_ext $1`"
+  shift
+done
+
+############################################################################
+
+mkdir $OUTPUT
+LOGDIR=${OUTPUT}/scripts+logs
+mkdir $LOGDIR
+echo $ORIG_COMMAND > $LOGDIR/command
+if [ "$DESIGN" != -1 ] ; then
+  /bin/cp $dm $OUTPUT/design.mat
+  /bin/cp $dc $OUTPUT/design.con
+fi
+
+echo "creating common mask"
+j=0
+for i in $INPUTS ; do
+  echo "$FSLDIR/bin/fslmaths $i -Tstd -bin ${OUTPUT}/mask_`${FSLDIR}/bin/zeropad $j 5` -odt char" >> ${LOGDIR}/drA
+  j=`echo "$j 1 + p" | dc -`
+done
+ID_drA=`$FSLDIR/bin/fsl_sub -T 10 -N drA -l $LOGDIR -t ${LOGDIR}/drA`
+cat <<EOF > ${LOGDIR}/drB
+#!/bin/sh
+\$FSLDIR/bin/fslmerge -t ${OUTPUT}/maskALL \`\$FSLDIR/bin/imglob ${OUTPUT}/mask_*\`
+\$FSLDIR/bin/fslmaths $OUTPUT/maskALL -Tmin $OUTPUT/mask
+\$FSLDIR/bin/imrm $OUTPUT/mask_*
+EOF
+chmod a+x ${LOGDIR}/drB
+ID_drB=`$FSLDIR/bin/fsl_sub -j $ID_drA -T 5 -N drB -l $LOGDIR ${LOGDIR}/drB`
+
+echo "doing the dual regressions"
+j=0
+for i in $INPUTS ; do
+  s=subject`${FSLDIR}/bin/zeropad $j 5`
+  echo "$FSLDIR/bin/fsl_glm -i $i -d $ICA_MAPS -o $OUTPUT/dr_stage1_${s}.txt --demean -m $OUTPUT/mask ; \
+        $FSLDIR/bin/fsl_glm -i $i -d $OUTPUT/dr_stage1_${s}.txt -o $OUTPUT/dr_stage2_$s --out_z=$OUTPUT/dr_stage2_${s}_Z --demean -m $OUTPUT/mask $DES_NORM ; \
+        $FSLDIR/bin/fslsplit $OUTPUT/dr_stage2_$s $OUTPUT/dr_stage2_${s}_ic" >> ${LOGDIR}/drC
+  j=`echo "$j 1 + p" | dc -`
+done
+ID_drC=`$FSLDIR/bin/fsl_sub -j $ID_drB -T 30 -N drC -l $LOGDIR -t ${LOGDIR}/drC`
+
+echo "sorting maps and running randomise"
+j=0
+Nics=`$FSLDIR/bin/fslnvols $ICA_MAPS`
+while [ $j -lt $Nics ] ; do
+  jj=`$FSLDIR/bin/zeropad $j 4`
+
+  RAND=""
+  if [ $NPERM -eq 1 ] ; then
+    RAND="$FSLDIR/bin/randomise -i $OUTPUT/dr_stage2_ic$jj -o $OUTPUT/dr_stage3_ic$jj -m $OUTPUT/mask $DESIGN -n 1 -V -R"
+  fi
+  if [ $NPERM -gt 1 ] ; then
+    # EDIT HERE
+    RAND="$FSLDIR/bin/randomise -i $OUTPUT/dr_stage2_ic$jj -o $OUTPUT/dr_stage3_ic$jj -m $OUTPUT/mask $DESIGN -n $NPERM -T -V"
+  fi
+
+  echo "$FSLDIR/bin/fslmerge -t $OUTPUT/dr_stage2_ic$jj \`\$FSLDIR/bin/imglob $OUTPUT/dr_stage2_subject*_ic${jj}.*\` ; \
+        $FSLDIR/bin/imrm \`\$FSLDIR/bin/imglob $OUTPUT/dr_stage2_subject*_ic${jj}.*\` ; $RAND" >> ${LOGDIR}/drD
+  j=`echo "$j 1 + p" | dc -`
+done
+ID_drD=`$FSLDIR/bin/fsl_sub -j $ID_drC -T 60 -N drD -l $LOGDIR -t ${LOGDIR}/drD`
+

melgmix.cc

+/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
+              independent components
+    
+    melgmix.cc - Gaussian Mixture Model
+
+    Christian F. Beckmann, FMRIB Analysis Group
+    
+    Copyright (C) 1999-2013 University of Oxford */
+
+/*  CCOPYRIGHT */
+
+#include "newimage/newimageall.h"
+#include "melgmix.h"
+#include "utils/log.h"
+#include "miscmaths/miscprob.h"
+#include <time.h>
+#include "libvis/miscplot.h"
+#include "libvis/miscpic.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+using namespace MISCPLOT;
+using namespace MISCPIC;
+
+string float2str(float f,int width, int prec, bool scientif){
+ 	ostringstream os;
+  int redw = int(std::abs(std::log10(std::abs(f))))+1;
+  if(width>0)
+    os.width(width);
+  if(scientif)
+    os.setf(ios::scientific);
+  os.precision(redw+std::abs(prec));
+  os.setf(ios::internal, ios::adjustfield);
+  os << f;
+  return os.str();
+} 
+  
+namespace Melodic{
+  
+  void MelGMix::setup(const RowVector& dat, 
+		const string dirname,
+		int cnum, volume<float> themask, 
+		volume<float> themean, 
+		int num_mix, float eps, bool fixit){
+    	cnumber = cnum;
+    	Mask = themask;
+    	Mean = themean;
+    	prefix = string("IC_")+num2str(cnum);
+    	
+    	fitted = false;
+    	nummix = num_mix;
+    	numdata = dat.Ncols();
+    
+    	//normalise data 
+    	datamean = mean(dat,2).AsScalar();
+    	datastdev= stdev(dat,2).AsScalar();
+    	data=(dat - datamean)/datastdev;
+
+			dbgmsg(" mapdat; mean: " << datamean << " std: " <<datastdev << endl);
+
+    	props=zeros(1,nummix);
+    	vars=zeros(1,nummix);
+    	means=zeros(1,nummix);
+    	Params=zeros(1,nummix);
+    	logprobY = 1.0;
+
+    	props = std::pow(float(nummix),float(-1.0));
+
+    	Matrix tmp1;
+    	tmp1 = data * data.t() / numdata;
+    	vars = tmp1.AsScalar();
+
+    	float Dmin, Dmax, IntSize;
+    	Dmin =  min(data).AsScalar(); Dmax = max(data).AsScalar();
+    	IntSize = Dmax / nummix;
+
+    	means(1)=mean(data,2).AsScalar(); 
+    	for (int ctr=2; ctr <= means.Ncols(); ctr++){
+      	means(ctr) =  Dmax - (ctr - 1.5) * IntSize; 
+    	}
+    	means(2)=means(1)+2*sqrt(vars(1));
+    	if(nummix>2)
+        means(3)=means(1)-2*sqrt(vars(1));
+
+    	epsilon = eps;
+    	if(epsilon >=0 && epsilon < 0.0000001)
+      	{epsilon = std::log(float(dat.Ncols()))/1000 ;}
+    	fixdim = fixit;
+			dbgmsg(" parameters; " << means << " : " << vars << " : " << props << endl);
+		}
+
+  Matrix MelGMix::threshold(const RowVector& dat, string levels){ 
+    Matrix Res;
+    Res = 1.0;
+    string tmpstr;
+    tmpstr = levels;
+    //cerr << " Levels : " << levels << endl << endl;
+    Matrix levls(1,4);
+    levls = 0;
+    Matrix fpr;
+    Matrix fdr;
+    Matrix nht;
+
+    char *p;
+    char t[1024];
+    const char *discard = ", [];{(})abceghijklmoqstuvwxyzABCEGHIJKLMNOQSTUVWXYZ~!@#$%^&*_-=+|\':></?";
+
+    strcpy(t, tmpstr.c_str());
+    p=strtok(t,discard);
+    while(p){
+      Matrix New(1,1);
+      New(1,1) = atof(p);
+      if(strchr(p,'d')){
+				levls(1,3)++;
+				if(fdr.Storage()>0){
+	  			fdr = fdr | New;
+				}else{
+	  			fdr = New;
+				}
+      }else{
+				if(strchr(p,'p')){
+	  			levls(1,2)++;
+	  			if(fpr.Storage()>0){
+	    			fpr = fpr | New;
+	  			}else{
+	    			fpr = New;
+	  			}
+				}else{
+	  			if(strchr(p,'n')){
+	    			levls(1,4)++;
+	    			if(nht.Storage()>0){
+	      			nht = nht | New;
+	    			}else{
+	      			nht = New;
+	    			}
+	  			}else{
+	    			levls(1,1)++;
+	    			levls = levls | New;
+	  			}
+				}
+      }
+      p=strtok(NULL,discard);
+    }
+
+    if(fpr.Storage()>0){levls = levls | fpr;}
+    if(fdr.Storage()>0){levls = levls | fdr;}
+    if(nht.Storage()>0){levls = levls | nht;}
+    
+  //  cerr << " levles : " << levls << endl << endl;
+    Res = threshold(data, levls);
+    set_threshmaps(Res);
+
+    return Res;
+  }
+
+  Matrix MelGMix::threshold(const RowVector& dat, Matrix& levels){  
+  	Matrix tests;
+    tests=levels;
+    Matrix Nprobs;
+
+    //if only single Gaussian: resort to nht
+    if(nummix==1||props(1)>0.999||probmap.Sum()<0.05){
+    	if(levels(1,4)==0){
+				Matrix New(1,6);
+				New(1,5)=0.05;
+				New(1,6)=0.01;
+				New(1,4)=2;New(1,1)=0;New(1,2)=0;New(1,3)=0;tests=New;
+      }else{
+				Matrix New;
+				New = levels.Columns(int(1+levels(1,1)+levels(1,2)
+				 	+levels(1,3)),levels.Ncols());
+				New(1,4) = levels(1,4);
+				New(1,1)=0;New(1,1)=0;New(1,3)=0;
+				tests=New;
+      }
+    }
+
+    int numtests = int(tests(1,1)+tests(1,2)+tests(1,3)+tests(1,4));    
+    Matrix Res(numtests,numdata);
+    Res = 0.0;
+    int next = 1;
+
+    for(int ctr1=1;ctr1<=tests(1,1);ctr1++){
+      if(4+next <= tests.Ncols()){
+				message("   alternative hypothesis test at p > " << tests(1,4+next) << endl);
+				add_infstr(" alternative hypothesis test at p > "+float2str(tests(1,4+next),0,2,false));
+				Matrix tmpNull;
+				tmpNull = dat;
+/*	
+				float cutoffpos, cutoffneg;
+				cutoffpos = means(1)+100*std::sqrt(vars(1)+0.0000001);
+				cutoffneg = means(1)-100*std::sqrt(vars(1)+0.0000001);
+	
+				for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
+	  			if( probmap(ctr) > tests(1,4+next) ){
+	    			if( dat(ctr) > means(1) )
+	      			cutoffpos = std::min(cutoffpos, float(dat(ctr)));
+	    			else
+	      			cutoffneg = std::max(cutoffneg, float(dat(ctr)));
+	 				}
+	
+				for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
+	  			if( (dat(ctr) > cutoffneg) && (dat(ctr) < cutoffpos) )
+	    			tmpNull(1,ctr)=0.0;
+*/
+				for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
+					if( probmap(ctr) < tests(1,4+next) ){	
+						tmpNull(1,ctr)=0.0;
+						}
+							
+				Res.Row(next) << tmpNull;
+      }
+      next++;
+    }
+    
+    for(int ctr1=1;ctr1<=tests(1,2);ctr1++){
+      if(4+next <=tests.Ncols()){
+				cerr << " false positives control " << tests(1,4+next)<<endl;
+				Matrix tmp1;
+				tmp1 = normcdf(dat,means(1),vars(1));
+				Matrix tmpNull;
+				tmpNull = dat; 
+				for(int ctr=1; ctr<=tmp1.Ncols(); ctr++)
+	  			if(tmp1(1,ctr) < tests(1,4+next))
+	    		tmpNull(1,ctr)=0.0;
+				Res.Row(next) << tmpNull;
+      }
+      next++;
+    }
+
+    for(int ctr1=1;ctr1<=tests(1,3);ctr1++){
+      if(4+next <=tests.Ncols()){
+				message("   Local False Discovery Rate control at p < " << tests(1,4+next) << endl);
+				add_infstr(" Local False Discovery Rate control at p < "+float2str(tests(1,4+next),0,2,false));
+				RowVector tmp=dat;
+				SortAscending(tmp);
+				RowVector newcdf(tmp);
+	  		newcdf << normcdf(tmp,means(1),vars(1));
+
+				float thrp = tmp(tmp.Ncols())+0.01;
+				float thrn = tmp(1)-0.01;
+				int ctr=tmp.Ncols();
+				do{
+					thrp = tmp(ctr);
+					ctr-=1;
+				}while(ctr>0 && ( (1.0-newcdf(ctr))*tmp.Ncols() < (tests(1,4+next)*(tmp.Ncols()-ctr+1))   ));
+
+				ctr=1;
+				do{
+					thrn = tmp(ctr);
+					ctr+=1;
+				}while(ctr<=tmp.Ncols() && ( (newcdf(ctr))*tmp.Ncols() < (tests(1,4+next)*ctr)));
+
+				tmp = dat;
+				for(ctr=1; ctr<=tmp.Ncols();ctr++)
+					if((tmp(ctr) < thrp)&&(tmp(ctr) > thrn))
+					tmp(ctr) = 0.0;
+				Res.Row(next) << tmp;
+      }
+			next++;
+    }
+
+    for(int ctr1=1;ctr1<=tests(1,4);ctr1++){
+      if(4+next <=tests.Ncols()){
+				message("   2-sided null hypothesis test at " << tests(1,4+next)<<endl);
+				add_infstr(" 2-sided null hypothesis test at "+float2str(tests(1,4+next),0,2,false));
+				double mu, sig;
+				mu  = dat.Sum()/numdata;
+				sig = var(dat,2).AsScalar();
+				Matrix tmp1;
+				tmp1 = normcdf(dat,mu,std::abs(sig));
+				Matrix tmpNull;
+				tmpNull = dat; 
+				for(int ctr=1; ctr<=tmp1.Ncols(); ctr++)
+	  			if((tmp1(1,ctr) < 1-0.5*(tests(1,4+next))&&
+	      		(tmp1(1,ctr) > 0.5*(tests(1,4+next)))))
+	    				tmpNull(1,ctr)=0.0;
+				Res.Row(next) << tmpNull;
+      }
+      next++;
+    }
+   
+    return Res;
+	}
+
+  /* GMM fitting  */
+
+  void MelGMix::gmmupdate(){
+    int it_ctr = 1;
+    bool exitloop = false;
+    float oldll;
+
+    Matrix tmp0;Matrix tmp1;Matrix prob_K__y_theta;
+    Matrix kdata;
+    RowVector prob_Y__theta;RowVector Nbar;
+    RowVector mubar;RowVector sigmabar;RowVector pibar;
+    
+    do{
+      oldll = logprobY;
+
+      //make sure all variances are acceptable
+      for(int ctr1=1; ctr1 <=vars.Ncols(); ctr1++)
+      	if(vars(ctr1)<0.0001){
+      	  vars(ctr1) = 0.0001;
+      	}
+
+      tmp0 = normpdf(data,means,vars);
+      tmp1 = SP(props.t()*ones(1,numdata),tmp0);      
+      prob_Y__theta << sum(tmp1,1);
+      logprobY = log(prob_Y__theta).Sum();
+      prob_K__y_theta = SP(tmp1,pow(ones(nummix,1)*prob_Y__theta,-1));
+      Nbar << sum(prob_K__y_theta,2).t();
+      pibar = Nbar / numdata;
+      kdata = ones(nummix,1)*data;
+      mubar <<SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1));    
+      kdata -= mubar.t()*ones(1,numdata);
+      kdata = pow(kdata,2);
+      sigmabar << SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1));
+      
+      means = mubar;
+      vars  = sigmabar;
+      props = pibar;
+
+      if(epsilon<0){exitloop = it_ctr >= -epsilon;}
+      else{exitloop = (((logprobY-oldll < epsilon)&&(it_ctr>20))
+		  	||(it_ctr>400));}      
+      it_ctr++;
+    }while(!exitloop);
+  }
+
+  void MelGMix::gmmfit(){
+    int i,j;
+
+    if(fixdim){
+      if(nummix>1){
+				gmmupdate();
+				add_params(means,vars,props,logprobY,MDL,Evi,true);
+      }else{
+				means.ReSize(1);
+				means = data.Sum()/numdata;
+				vars.ReSize(1);
+				vars = var(data,2);
+				props.ReSize(1);
+				props = 1.0;
+				gmmevidence();
+      }
+    }else{
+      RowVector Score(Params.Ncols());
+      do{
+				gmmupdate();
+				Score(nummix) = gmmevidence();    
+				int idx1,idx2;
+				RowVector pitmp = props;
+     
+				pitmp.MaximumAbsoluteValue1(idx1);
+				pitmp(idx1)=0.0;
+				pitmp.MaximumAbsoluteValue1(idx2);
+	
+				if(props.MaximumAbsoluteValue1(i)<0.9){
+	  			if((vars(idx2)>0.15)&&
+	     			(std::abs(means(idx2)-means(idx1))<0.5*vars(idx1))){
+	    				Score(nummix) = Score(nummix)/(2*(means(idx1)));
+	  				}	     
+				}
+	
+				add_params(means,vars,props,logprobY,MDL,Evi,true);
+     
+				gmmreducemm();
+				means = means.Columns(1,nummix);
+				vars  = vars.Columns(1,nummix);
+				props = props.Columns(1,nummix);
+
+      }while(nummix>1);
+      means.ReSize(1);
+      means = data.Sum()/numdata;
+      vars.ReSize(1);
+      vars = var(data,2);
+      props.ReSize(1);
+      props = 1.0;
+      Score(nummix) = gmmevidence();
+      add_params(means,vars,props,logprobY,MDL,Evi,true);
+      //identify best MM
+      Score.MinimumAbsoluteValue2(i,j);
+      means.ReSize(j);
+      vars.ReSize(j);
+      props.ReSize(j);
+      nummix = j;
+      int index; index = 3 + (j-1)*5; 
+      means = Params.SubMatrix(index,index,1,j);
+      vars  = Params.SubMatrix(index+1,index+1,1,j);
+      props = Params.SubMatrix(index+2,index+2,1,j);
+    }
+
+    props.MaximumAbsoluteValue2(i,j);
+    if(j>1){
+      float tmp;
+      tmp = means(1);means(1) = means(j);means(j)=tmp;
+      tmp = vars(1);vars(1) = vars(j);vars(j)=tmp;
+      tmp = props(1);props(1) = props(j);props(j)=tmp;
+    }
+    
+    add_params(means,vars,props,logprobY,MDL,Evi,true);
+
+    if(nummix==1)
+      probmap << normcdf(data,means(1),vars(1));
+    else{
+      Matrix Nprobs;
+      Matrix tmp0;
+      tmp0 = normpdf(data,means,vars);
+      Nprobs = SP(props.t()*ones(1,numdata),tmp0);
+      tmp0 = ones(Nprobs.Nrows(),1)*pow(sum(Nprobs,1),-1);
+      Nprobs = SP(tmp0,Nprobs);
+      probmap << SP(sum(Nprobs.Rows(2,Nprobs.Nrows()),1),
+		    pow(sum(Nprobs,1),-1));
+    } 
+  }
+
+  float MelGMix::gmmevidence(){
+    Matrix tmp0;
+    if(means.Ncols()>1){
+      tmp0 = normpdf(data,means,vars); 
+    }else{
+      tmp0 = normpdf(data,means.AsScalar(),vars.AsScalar());
+    }
+    Matrix tmp1;
+    tmp1 = SP(props.t()*ones(1,numdata),tmp0);
+    tmp0 = SP(tmp0,pow(ones(nummix,1)*sum(tmp1,1),-1));
+    tmp0 = pow(tmp0-ones(nummix,1)*tmp0.Row(nummix),2);
+    float logH = 0;
+    if(means.Ncols()>1){
+      logH = sum(log(sum(tmp0.Rows(1,nummix-1),2)),1).AsScalar();
+    }
+    logH = logH + 2*sum(log(std::sqrt(2.0)*numdata*props),2).AsScalar();
+    logH = logH - 2*sum(props,2).AsScalar();
+    
+    RowVector prob_Y__theta;
+    prob_Y__theta << sum(tmp1,1);    
+    logprobY = log(prob_Y__theta).Sum();     
+    MDL = -logprobY + (1.5*nummix-0.5)* std::log(float(numdata));   
+    Evi = -logprobY +nummix*std::log(2.0)-std::log(MISCMATHS::gamma(nummix))
+      -3*nummix/2*std::log(M_PI)+0.5*logH;
+  
+    return Evi;
+  }
+
+  void MelGMix::gmmreducemm(){
+    Matrix dlm(nummix,nummix);
+    Matrix mus(nummix,nummix);
+    Matrix rs(nummix,nummix);
+
+    for(int ctri=1;ctri<=nummix; ctri++){
+      for(int ctrj=1;ctrj<=nummix; ctrj++){
+				mus(ctri,ctrj) = (props(ctri)*means(ctri)+props(ctrj)*means(ctrj))
+	      	/(props(ctri)+props(ctrj));
+        rs(ctri,ctrj)  = (props(ctri)*(vars(ctri)+
+			  std::pow(means(ctri)-mus(ctri,ctrj),2) ) 
+        	+ props(ctrj)*(vars(ctrj) 
+         	+ std::pow(means(ctrj)-mus(ctri,ctrj),2))) 
+	        / (props(ctri)+props(ctrj));
+				dlm(ctri,ctrj) = 0.5*numdata*(
+			 		props(ctri)*std::log(
+          std::abs(rs(ctri,ctrj))/std::abs(vars(ctri))) 
+			 		+ props(ctrj)*std::log(std::abs(rs(ctri,ctrj))
+          / std::abs(vars(ctrj))));
+      }
+    }
+
+    dlm += IdentityMatrix(nummix)*dlm.Maximum();
+
+    int i,j;
+    float val;
+    val=dlm.MinimumAbsoluteValue2(i,j);
+
+    nummix--;
+    
+    RowVector newmean(nummix);
+    RowVector newvars(nummix);
+    RowVector newprop(nummix);
+
+    int ctr0 = 1;
+    for(int ctr=1; ctr<=nummix+1; ctr++){
+      if(ctr!=i&&ctr!=j){
+				newmean(ctr0) = means(ctr);
+				newvars(ctr0) = vars(ctr);
+				newprop(ctr0) = props(ctr);
+				ctr0++;
+      }
+    }
+    //cerr << "ctr0   " << ctr0 << endl;
+    if(ctr0<=nummix){
+      newmean(ctr0) = mus(i,j);
+      newvars(ctr0) = rs(i,j);
+      newprop(ctr0) = props(i)+props(j);
+      
+      means = newmean;    
+      vars=newvars;
+      props=newprop;
+    }
+  }
+
+  void MelGMix::ggmfit(){
+  	// fit a mixture of a Gaussian and multiple Gamma functions to the histogram
+  
+    float log_p_y_theta = 1.0;
+    float old_ll = 2.0;
+    float g_eps = 0.000001;
+    int it_ctr = 0;
+    double Dmax, Dmin;
+   
+    Dmax = 2 * data.Maximum();
+    Dmin = -2 * data.Minimum();
+
+    //resize means, vars and props
+    if(nummix > 3)
+      nummix = 3;
+    means = means.Columns(1,nummix);
+    vars  = vars.Columns(1,nummix);
+    props = props.Columns(1,nummix);
+
+    means(1) = -2*mean(data,2).AsScalar();
+
+    Matrix tmp1;Matrix prob_K__y_theta;
+    Matrix kdata;
+    RowVector prob_Y__theta;RowVector Nbar;
+    RowVector mubar;RowVector sigmabar;RowVector pibar;
+    Matrix p_ygx(nummix,numdata);
+    offset = 0.0;
+    float const2;
+    Matrix negdata(data);
+    negdata = -data;
+
+    while((it_ctr<30) ||
+	  ((std::abs(old_ll - log_p_y_theta)>g_eps) && (it_ctr<500))){ // fit GGM	
+ 			it_ctr++;
+			//offset = (std::min(0.2,1-props(1)))*std::sqrt(vars(1));
+
+			//make sure all variances are acceptable
+ 			for(int ctr1=1; ctr1 <=nummix; ctr1++)
+ 	  		if(vars(ctr1)<0.0001){
+ 	    		vars(ctr1) = 0.0001;
+ 	  		}
+
+ 				p_ygx = 0.0;
+ 				p_ygx.Row(1) << normpdf(data,means(1),vars(1));
+       
+ 				const2 = (2.6-props(1))*sqrt(vars(1))+means(1); //min. nht level
+ 
+				means(2) = (std::max(means(2), std::max(0.001,
+	   			0.5 * ( const2 + std::sqrt( const2*const2 + 4*vars(2))))));
+				vars(2)  = std::max(std::min(vars(2), 0.5*std::pow(means(2),2)),0.0001);
+				p_ygx.Row(2) << gammapdf(data,means(2),vars(2));
+   
+ 				if(nummix>2){
+	  			const2 = (2.6-props(1))*sqrt(vars(1))-means(1);
+	
+	  			means(3) = -(std::max(-means(3), std::max(0.001,
+	      		0.5 * ( const2 + std::sqrt( const2*const2 + 4*vars(3))))));
+	  			vars(3)  = std::max(std::min(vars(3), 0.5*std::pow(means(3),2)),0.0001);
+ 	  			p_ygx.Row(3) << gammapdf(negdata,-means(3),vars(3));
+				}
+
+ 				tmp1 = SP(props.t()*ones(1,numdata),p_ygx);
+ 				prob_Y__theta << sum(tmp1,1);
+	
+				//deal with non-modelled voxels
+				for(int ctr=1; ctr<=tmp1.Ncols(); ctr++)
+	  			if(prob_Y__theta(ctr) < 0.0001)
+	    			prob_Y__theta(ctr) = 0.0001;
+
+ 				old_ll = log_p_y_theta;
+ 				log_p_y_theta = log(prob_Y__theta).Sum();
+ 				if((it_ctr<30) ||
+	   			((std::abs(old_ll - log_p_y_theta)>g_eps) && (it_ctr<300))){//update
+	  
+ 	  			prob_K__y_theta = SP(tmp1,pow(ones(nummix,1)*prob_Y__theta,-1));
+ 	  			Nbar << sum(prob_K__y_theta,2).t();
+	  			for(int ctr=1; ctr<=Nbar.Ncols(); ctr++)
+	    			if(Nbar(ctr) < 0.0001 * numdata)
+	      			Nbar = Nbar + 0.0001;
+ 	  			pibar= Nbar / numdata;
+	  			// 	  cerr << "pibar :" << pibar << endl;
+ 	  			kdata = ones(nummix,1)*data;
+ 	  			mubar <<SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1)); 
+	  			// 	  cerr << "mubar :" << mubar << endl;
+
+ 	  			kdata -= mubar.t()*ones(1,numdata);
+ 	  			kdata = pow(kdata,2);
+ 	  			sigmabar << SP(sum(SP(kdata,prob_K__y_theta),2).t(),pow(Nbar,-1));
+      
+ 	  			means = mubar;
+ 	  			vars  = sigmabar;
+ 	  			props = pibar;
+ 					}//update
+    } //while loop
+
+    props = props / sum(props,2).AsScalar();
+    add_params(means,vars,props,logprobY,MDL,Evi,true);
+	
+    probmap << SP(sum(tmp1.Rows(2,tmp1.Nrows()),1),
+		  pow(sum(tmp1,1),-1));
+
+	dbgmsg("   mu: " << means << "  sig: " << vars << " prop: " << props << endl);
+
+   	if(props(1)<0.4 ){
+      //set up GMM
+      message("    try Gaussian Mixture Model " << endl);
+      props=zeros(1,nummix);
+      vars=zeros(1,nummix);
+      means=zeros(1,nummix);
+      Params=zeros(1,nummix);
+      logprobY = 1.0;  
+      props = std::pow(float(nummix),float(-1.0));
+
+      tmp1 = data * data.t() / numdata;
+      vars = tmp1.AsScalar();
+      float Dmin, Dmax, IntSize;
+      Dmin = min(data).AsScalar(); Dmax = max(data).AsScalar();
+      IntSize = Dmax / nummix;
+      means(1)=mean(data,2).AsScalar(); 
+      for (int ctr=2; ctr <= means.Ncols(); ctr++){
+		means(ctr) =  Dmax - (ctr - 1.5) * IntSize; 
+      }
+      means(2)=means(1)+sqrt(vars(1));
+      if(nummix>2)
+		means(3)=means(1)-sqrt(vars(1));
+      
+      fit(string("GMM"));
+    }
+
+  }
+
+  /*  INPUT / OUTPUT  */
+
+  void MelGMix::add_params(Matrix& mu, Matrix& sig, Matrix& pi, 
+		float logLH, float MDL, float Evi, bool advance){ 
+    int size = Params.Ncols();
+    if(size<2){size=2;}
+    Matrix New(5,size);
+    New = 0;
+    
+    New.SubMatrix(3,3,1,mu.Ncols())=mu;
+    New.SubMatrix(4,4,1,mu.Ncols())=sig;
+    New.SubMatrix(5,5,1,mu.Ncols())=pi;
+    New(1,1)=nummix;
+  
+    New(1,2)=-logLH;
+    New(2,1)=Evi;
+    New(2,2)=MDL;
+    if(Params.Storage()>nummix){ 
+      Params = New & Params;
+    }else{ 
+      Params =  New;
+    }
+    }
+
+  void MelGMix::get_params(int index, Matrix& mu, Matrix& sig, Matrix& pi, 
+		float logLH, float MDL, float Evi){ 
+   
+    }
+
+  void MelGMix::save(){
+    
+  }
+
+  void MelGMix::status(const string &txt){
+    cerr << txt << "epsilon " << epsilon << endl;
+    cerr << txt << "nummix  " << nummix << endl;
+    cerr << txt << "numdata " << numdata << endl;
+    cerr << txt << "means   " << means << endl;
+    cerr << txt << "vars    " << vars << endl;
+    cerr << txt << "props   " << props << endl;
+    //write_ascii_matrix(logger.appendDir(string(txt + "means")),means);
+    //write_ascii_matrix(logger.appendDir(string(txt + "vars")),vars);
+    //write_ascii_matrix(logger.appendDir(string(txt + "props")),props);
+  }
+
+  void MelGMix::create_rep(){
+ 
+  }
+  
+
+}
+
+
+

melgmix.h

+/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
+              independent components
+    
+    melgmix.h - class for Gaussian/Gamma Mixture Model
+
+    Christian F. Beckmann, FMRIB Analysis Group
+    
+    Copyright (C) 1999-2013 University of Oxford */
+
+/*  CCOPYRIGHT */
+
+#ifndef __MELGMIX_h
+#define __MELGMIX_h
+
+#include "utils/log.h"
+#include "melodic.h"
+#include "utils/options.h"
+#include "meloptions.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+  
+  class MelGMix{
+    public:
+ 
+      MelGMix(MelodicOptions &popts, Log &plogger):
+				opts(popts),
+				logger(plogger){}
+
+      ~MelGMix() { 
+				//mainhtml << endl << "<hr></CENTER></BODY></HTML>" << endl;
+	  	}
+
+      void save();
+
+      void setup(const RowVector& dat, const string dirname,
+		 		int here, volume<float> themask, 
+		 		volume<float> themean, int num_mix = 3, 
+		 		float eps = 0.0, bool fixdim = false);
+      
+      void gmmfit();
+      void ggmfit();
+
+      inline void fit(string mtype = string("GGM")){
+	  		mmtype = mtype;
+	  		if(mmtype==string("GGM")) 
+	    		this->ggmfit(); 
+	  		else 
+	    		this->gmmfit();
+
+	  		//re-insert mean and stdev
+	  		data = data*datastdev + datamean;
+	  		//threshmaps = threshmaps*datastdev + datamean;
+	  		means = means*datastdev + datamean;
+	  		vars = vars*datastdev*datastdev;
+			}
+
+      inline Matrix threshold(string levels){
+				return this->threshold(data, levels);
+			}
+      inline Matrix threshold(RowVector& levels){ 
+				return this->threshold(data, levels);
+			}
+      Matrix threshold(const RowVector& dat, Matrix& levels);
+      Matrix threshold(const RowVector& dat, string levels);
+
+      void status(const string &txt);
+ 
+      inline RowVector& get_means() {return means;}
+      inline void set_means(RowVector& Arg) {means = Arg;}
+    
+      inline RowVector& get_vars() {return vars;}
+      inline void set_vars(RowVector& Arg) {vars = Arg;}
+      
+      inline RowVector& get_pi() {return props;}
+      inline void set_pi(RowVector& Arg) {props = Arg;}
+      
+      inline RowVector& get_data() {return data;}
+      inline void set_data(RowVector& Arg) {data = Arg;}
+
+      inline RowVector& get_prob() {return probmap;}
+
+      inline float get_eps() {return epsilon;}
+      inline void set_eps(float Arg) {epsilon = Arg;}
+
+      inline Matrix& get_threshmaps() {return threshmaps;}
+      inline void set_threshmaps(Matrix& Arg) {threshmaps = Arg;}
+
+      inline bool isfitted(){return fitted;}
+
+      inline int mixtures(){return nummix;}
+     
+      inline string get_type() { return mmtype;}
+      inline void set_type(string Arg) { mmtype = Arg;}
+      
+      inline string get_prefix() { return prefix;}
+      inline void  set_prefix(string Arg) { prefix = Arg;}
+
+      inline RowVector get_probmap() {return probmap;}
+
+      inline float get_offset() {return offset;}
+      inline void set_offset(float Arg) {offset = Arg;}
+
+      inline void flipres(int num){
+				means = -means;
+				data = -data;
+				threshmaps = -threshmaps;
+				if(mmtype=="GGM"){
+	  			float tmp;
+	  			tmp= means(2);means(2)=means(3);means(3)=tmp;
+	  			tmp=vars(2);vars(2)=vars(3);vars(3)=tmp;
+	  			tmp=props(2);props(2)=props(3);props(3)=tmp;
+				}
+      }      
+
+      void create_rep();
+
+      inline void add_infstr(string what){
+				threshinfo.push_back(what);
+      }
+
+      inline string get_infstr(int num){
+				if((threshinfo.size()<(unsigned int)(num-1))||(num<1))
+	  			return string("");
+				else
+	  			return threshinfo[num-1];
+      }
+
+      inline int size_infstr(){
+      	return threshinfo.size();
+      }
+
+      inline void clear_infstr(){
+				threshinfo.clear();
+      }
+
+      inline void smooth_probs(float howmuch){
+				volume4D<float> tempVol;
+				tempVol.setmatrix(probmap,Mask);
+        tempVol[0]= smooth(tempVol[0],howmuch);
+        probmap = tempVol.matrix(Mask);
+      }
+
+	  inline Matrix get_params(){
+		Matrix tmp = zeros(3,means.Ncols());
+		tmp.Row(1) = means;
+		tmp.Row(2) = vars;
+		tmp.Row(3) = props;
+		return tmp;
+      }
+
+      double datamean;
+      double datastdev;
+
+    private:
+      __attribute__((unused)) MelodicOptions &opts;     
+      __attribute__((unused)) Log &logger; //global log file
+
+      //Log mainhtml;
+
+      void gmmupdate();
+      float gmmevidence();
+      void gmmreducemm();
+      void add_params(Matrix& mu, Matrix& sig, Matrix& pi, 
+		    float logLH, float MDL, float Evi, bool advance = false);
+      void get_params(int index, Matrix& mu, Matrix& sig, Matrix& pi, 
+		    float logLH, float MDL, float Evi);
+
+      Matrix Params;
+      Matrix threshmaps;
+
+      RowVector means;
+      RowVector vars;
+      RowVector props;
+      RowVector data;
+      RowVector probmap;
+
+      volume<float> Mean;
+      volume<float> Mask;
+
+      float epsilon;
+      float logprobY;
+      float MDL;
+      float Evi;
+      float offset;
+
+      int nummix;
+      int numdata;
+      int cnumber;
+
+      bool fitted;
+      bool fixdim;
+
+      string prefix;
+      string mmtype;
+      string dirname;
+
+      vector<string> threshinfo;
+
+  };
+}
+
+#endif

melhlprfns.h

+/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
+              independent components
+    
+    melhlprfns.cc - misc functions
+
+    Christian F. Beckmann, FMRIB Analysis Group
+    
+    Copyright (C) 1999-2013 University of Oxford */
+
+/*  CCOPYRIGHT  */
+
+#ifndef __MELODICHLPR_h
+#define __MELODICHLPR_h
+
+#include "newimage/newimageall.h"
+
+	#ifdef __APPLE__
+	#include <mach/mach.h>
+	#define mmsg(msg) { \
+	  struct task_basic_info t_info; \
+	  mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
+	  if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
+		{ \
+			cout << " MEM: " << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
+			} \
+	}
+	#else
+	#define mmsg(msg) { \
+	   cout << msg; \
+	}
+	#endif
+
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+  void update_mask(volume<float>& mask, Matrix& Data);
+  void del_vols(volume4D<float>& in, int howmany);
+
+  Matrix smoothColumns(const Matrix& inp);
+  Matrix calc_FFT(const Matrix& Mat, const bool logpwr = 0);
+
+  Matrix convert_to_pbsc(Matrix& Mat);
+
+  RowVector varnorm(Matrix& in, int dim = 30, float level = 1.6, int econ = 20000);
+       void varnorm(Matrix& in, const RowVector& vars);
+  RowVector varnorm(Matrix& in, SymmetricMatrix& Corr, int dim = 30, float level = 1.6, int econ = 20000);
+
+  Matrix SP2(const Matrix& in, const Matrix& weights, int econ = 20000);
+  void SP3(Matrix& in, const Matrix& weights);
+
+  RowVector Feta(int n1,int n2);
+  RowVector cumsum(const RowVector& Inp);
+
+  Matrix corrcoef(const Matrix& in1, const Matrix& in2);
+  Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part);
+  float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
+  float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite);
+  void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
+  void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& white, Matrix& dewhite);
+  void calc_white(const SymmetricMatrix& Corr, int dim, Matrix& white, Matrix& dewhite);
+  
+  void std_pca(const Matrix& Mat, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
+  void std_pca(const Matrix& Mat, const Matrix& weights, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ = 20000);
+  void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
+  void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
+
+  float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim = 1);
+  RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows);
+  RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows);
+  Matrix krprod(const Matrix& cols, const Matrix& rows);
+  Matrix krapprox(const Matrix& Mat, int size_col, int dim = 1);
+
+  void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2, RowVector& out3, int& out4, int num_vox, float resels);
+  void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2);
+
+  int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, SymmetricMatrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which);
+  int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, float resels, string which);
+  int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which);
+  ColumnVector ppca_select(Matrix& PPCAest, int& dim, int maxEV, string which);
+  Matrix ppca_est(const RowVector& eigenvalues, const int N1, const float N2);
+  Matrix ppca_est(const RowVector& eigenvalues, const int N);
+
+  ColumnVector acf(const ColumnVector& in, int order);
+  ColumnVector pacf(const ColumnVector& in, int maxorder = 1);
+  Matrix est_ar(const Matrix& Mat, int maxorder);
+  ColumnVector gen_ar(const ColumnVector& in, int maxorder = 1);
+  Matrix gen_ar(const Matrix& in, int maxorder);
+  Matrix gen_arCorr(const Matrix& in, int maxorder);
+ 
+	class basicGLM{
+		public:
+		
+			//constructor
+			basicGLM(){}
+		
+			//destructor
+			~basicGLM(){}
+		
+			void olsfit(const Matrix& data, const Matrix& design, 
+				const Matrix& contrasts, int DOFadjust = -1);
+
+			inline Matrix& get_t(){return t;}
+			inline Matrix& get_z(){return z;}
+			inline Matrix& get_p(){return p;}
+			inline Matrix& get_f_fmf(){return f_fmf;}
+			inline Matrix& get_pf_fmf(){return pf_fmf;}
+			inline Matrix& get_cbeta(){return cbeta;}
+			inline Matrix& get_beta(){return beta;}
+			inline Matrix& get_varcb(){return varcb;}
+			inline Matrix& get_sigsq(){return sigsq;}
+			inline Matrix& get_residu(){return residu;}
+			inline int get_dof(){return dof;}
+			
+		private:
+			Matrix beta;
+			Matrix residu;
+			Matrix sigsq;
+			Matrix varcb;
+			Matrix cbeta;
+			Matrix f_fmf, pf_fmf;
+			int dof;
+			Matrix t;
+			Matrix z;
+			Matrix p;
+  };
+//	Matrix glm_ols(const Matrix& dat, const Matrix& design);
+}
+
+#endif

melica.cc

+/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
+              independent components
+    
+    melica.cc - ICA estimation 
+
+    Christian F. Beckmann, FMRIB Analysis Group
+    
+    Copyright (C) 1999-2013 University of Oxford */
+
+/*  CCOPYRIGHT  */
+
+#include <stdlib.h>
+#include "utils/log.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+#include "newmatap.h"
+#include "newmatio.h"
+#include "melica.h"
+#include "melpca.h"
+#include "melhlprfns.h"
+#include "miscmaths/miscprob.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic {
+    
+  void MelodicICA::ica_fastica_symm(const Matrix &Data){
+    // based on Aapo Hyvrinen's fastica method
+    // see www.cis.hut.fi/projects/ica/fastica/
+    
+    //initialize matrices
+    Matrix redUMM_old, rank1_old;
+    Matrix tmpU;    
+    //srand((unsigned int)timer(NULL));
+    redUMM = melodat.get_white()*
+       unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
+    
+	if(opts.debug.value())
+		cerr << "redUMM init submatrix : " << endl << redUMM.SubMatrix(1,2,1,2) << endl;
+		
+    if(opts.guessfname.value().size()>1){
+      message("  Use columns in " << opts.guessfname.value() 
+	      << " as initial values for the mixing matrix " <<endl);
+      Matrix guess ;
+      guess  = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
+      redUMM.Columns(1,guess.Ncols()) = guess;
+    }
+    
+    symm_orth(redUMM);
+
+    int itt_ctr,itt_ctr2=1,cum_itt=0,newmaxitts = opts.maxNumItt.value(); 
+    double minAbsSin = 1.0, minAbsSin2 = 1.0;
+    if(opts.approach.value() == string("tica")) {
+      opts.maxNumItt.set_T(opts.rank1interval.value());
+      rank1_old = melodat.get_dewhite()*redUMM;
+      rank1_old = melodat.expand_dimred(rank1_old);
+      rank1_old = krapprox(rank1_old,int(rank1_old.Nrows()/melodat.get_numfiles()));
+    }
+    do{// TICA loop
+      itt_ctr = 1;
+      do{ // da loop!!!
+				redUMM_old = redUMM;      
+				//calculate IC estimates
+				tmpU = Data.t() * redUMM;
+					
+				//update redUMM depending on nonlinearity
+				if(opts.nonlinearity.value()=="pow4"){
+	  			redUMM = (Data * pow(tmpU,3.0)) / samples - 3 * redUMM;
+				}
+				if(opts.nonlinearity.value()=="pow3"){
+	  			tmpU /= opts.nlconst1.value();
+	  			redUMM = 3 * (Data * pow(tmpU,2.0)) / samples  - 
+	    			(SP(ones(dim,1)*sum(tmpU,1),redUMM))/ samples;
+				}
+				if(opts.nonlinearity.value()=="tanh"){
+	  			Matrix hyptanh;
+	  			hyptanh = tanh(opts.nlconst1.value()*tmpU);
+	  			redUMM = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
+						sum(1-pow(hyptanh,2),1),redUMM))/samples;						
+				}
+				if(opts.nonlinearity.value()=="gauss"){
+	  			Matrix tmpUsq;
+	  			Matrix tmpU2;
+	  			Matrix gauss;
+	  			Matrix dgauss;
+	  			tmpUsq = pow(tmpU,2);
+	  			tmpU2 = exp(-(opts.nlconst2.value()/2) * tmpUsq);
+	  			gauss = SP(tmpU,tmpU2);
+	  			dgauss = SP(1-opts.nlconst2.value()*tmpUsq,tmpU2);
+	  			redUMM = (Data * gauss - SP(ones(dim,1)*
+				    sum(dgauss,1),redUMM))/samples;
+				}
+           
+				// orthogonalize the unmixing-matrix 
+				symm_orth(redUMM);
+				if(opts.approach.value() == string("tica")){
+					message("");
+				}
+
+				//termination condition : angle between old and new < epsilon
+				minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_old)).Minimum();
+
+				message("  Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl);
+		//		if((abs(minAbsSin) < opts.epsilon.value())&&
+		//		  (opts.approach.value()!=string("tica"))){ break;}
+				if((abs(minAbsSin) < opts.epsilon.value())){ break;}
+	
+				itt_ctr++;
+      } while(itt_ctr < opts.maxNumItt.value());
+      cum_itt += itt_ctr;
+      itt_ctr2++;
+      if(opts.approach.value() == string("tica")){
+				message(" Rank-1 approximation of the time courses; ");
+	  		Matrix temp(melodat.get_dewhite() * redUMM);
+	  		temp = melodat.expand_dimred(temp);
+			  temp = krapprox(temp,int(temp.Nrows()/melodat.get_numfiles())); 
+				minAbsSin2 = 1 - diag(abs(corrcoef(temp,rank1_old))).Minimum();
+				rank1_old = temp;
+				temp = melodat.reduce_dimred(temp);
+				redUMM = melodat.get_white() * temp;
+
+				message(" change : " << minAbsSin2 << endl);	
+				if(abs(minAbsSin2) < opts.epsilonS.value() && abs(minAbsSin) < opts.epsilon.value()){ break;}
+			}
+    } while(
+      (itt_ctr2 < newmaxitts/opts.maxNumItt.value()) && 
+			(opts.approach.value() ==  string("tica")) && 
+			cum_itt < newmaxitts);
+
+    if((itt_ctr>=opts.maxNumItt.value() && (opts.approach.value()!=string("tica")))
+			|| (cum_itt >= newmaxitts && opts.approach.value()==string("tica"))){
+      cerr << "  No convergence after " << cum_itt  <<" steps "<<endl;
+    } else {
+      message("  Convergence after " << cum_itt  <<" steps " << endl << endl);
+      no_convergence = false;
+
+      {Matrix temp(melodat.get_dewhite() * redUMM);
+       melodat.set_mix(temp);}
+      {Matrix temp(redUMM.t()*melodat.get_white());
+      melodat.set_unmix(temp);}
+    } 
+  }
+  
+  void MelodicICA::ica_fastica_defl(const Matrix &Data){
+    if(!opts.explicitnums || opts.numICs.value()>dim){
+      opts.numICs.set_T(dim); 
+      message("  Using numICs:" << opts.numICs.value() << endl);
+    }
+     
+    //redUMM = zeros(dim); // got to start somewhere
+    redUMM = melodat.get_white()*
+      unifrnd(melodat.get_white().Ncols(),opts.numICs.value());
+    redUMM = zeros(melodat.get_white().Nrows(),opts.numICs.value());
+    Matrix guess;
+    int guesses=0;
+    if(opts.guessfname.value().size()>1){
+       message("  Use columns in " << opts.guessfname.value() << " as initial values for the mixing matrix " <<endl);
+       guess  = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
+       guesses = guess.Ncols();
+    }
+
+    int ctrIC = 1;
+    int numRestart = 0;
+    while(ctrIC<=opts.numICs.value()){   
+     	message("  Extracting IC " << ctrIC << "  ... ");
+      ColumnVector w;
+      ColumnVector w_old;   
+      ColumnVector tmpU;
+      if(ctrIC <= guesses){
+      	w = w - redUMM * redUMM.t() * w;
+      	w = w / norm2(w);  
+      	w_old = zeros(w.Nrows(),1);
+      	int itt_ctr = 1; 
+      	do{
+	 				w_old = w;
+	 				tmpU = Data.t() * w; 
+					if(opts.nonlinearity.value()=="pow4"){
+	  				w =  (Data * pow(tmpU,3.0)) / samples - 3 * w;
+					}
+					if(opts.nonlinearity.value()=="tanh"){
+	  				ColumnVector hyptanh;
+          	hyptanh = tanh(opts.nlconst1.value()*tmpU);
+ 	  				w = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
+          		sum(1-pow(hyptanh,2),1),w))/samples;
+					} 
+					if(opts.nonlinearity.value()=="pow3"){
+	  				tmpU /= opts.nlconst1.value();
+ 	  				w = 3*(Data * pow(tmpU,2.0)) / samples - 2*(SP(ones(dim,1)*
+           		sum(tmpU,1),w))/samples;
+					} 
+					if(opts.nonlinearity.value()=="gauss"){
+	  				ColumnVector tmpUsq;
+	  				ColumnVector tmpU2;
+	  				ColumnVector gauss;
+	  				ColumnVector dgauss;
+	  				tmpUsq = pow(tmpU,2);
+	  				tmpU2 = exp(-(opts.nlconst2.value()/2) * tmpUsq);
+	  				gauss = SP(tmpU,tmpU2);
+	  				dgauss = SP(1-opts.nlconst2.value()*tmpUsq,tmpU2);
+          	w = (Data * gauss - SP(ones(dim,1)*
+				 			sum(dgauss,1),w))/samples;
+					}
+
+					// orthogonalize w
+					w = w - redUMM * redUMM.t() * w;
+					w = w / norm2(w);  
+
+					//termination condition : angle between old and new < epsilon
+					if(((norm2(w-w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10)) || 
+					   ((norm2(w+w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10))){
+	 					break;
+				  	}
+        	//cout << norm2(w-w_old) << "   " << norm2(w+w_old) << endl;
+					itt_ctr++;
+     	  } while(itt_ctr <= opts.maxNumItt.value());
+
+      if(itt_ctr<opts.maxNumItt.value()){
+				redUMM.Column(ctrIC) = w;
+        message(" estimated using " << itt_ctr << " iterations " << endl);
+        ctrIC++; 
+        numRestart = 0;
+      } else{
+        if(numRestart > opts.maxRestart.value()){
+	  			message(endl << "  Estimation failed after " 
+		  			<< numRestart << " attempts " 
+		  			<< " giving up " << endl);
+	  			break;
+				}else{
+          numRestart++;
+	  			message(endl <<"  Estimation failed - restart " 
+		  			<< numRestart << endl);
+				}
+      }
+     	}
+     	if(numRestart < opts.maxRestart.value()){
+       	no_convergence = false;
+       	{Matrix temp(melodat.get_dewhite() * redUMM);
+       		melodat.set_mix(temp);}
+       	{Matrix temp(redUMM.t()*melodat.get_white());
+       		melodat.set_unmix(temp);}
+     	}
+		}
+  }
+
+  void MelodicICA::ica_maxent(const Matrix &Data){
+    // based on Aapo Hyvrinen's fastica method
+    // see www.cis.hut.fi/projects/ica/fastica/ 
+    message(" MAXENT " << endl);
+    //initialize matrices
+    Matrix redUMM_old;
+    Matrix tmpU;    
+    Matrix gtmpU;
+    double lambda = 0.015/std::log((double)melodat.get_white().Ncols());
+    
+    //srand((unsigned int)timer(NULL));
+    redUMM = melodat.get_white()*
+     	unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
+    
+    if(opts.guessfname.value().size()>1){
+      message("  Use columns in " << opts.guessfname.value() 
+	      << " as initial values for the mixing matrix " <<endl);
+      Matrix guess ;
+      guess  = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
+      redUMM.Columns(1,guess.Ncols()) = guess;
+    }
+    
+    //    symm_orth(redUMM);
+    int itt_ctr=1; 
+    double minAbsSin = 1.0;
+    Matrix Id;
+    Id = IdentityMatrix(redUMM.Ncols());
+    //cerr << " nonlinearity : " <<    opts.nonlinearity.value() << endl;
+
+    do{ // da loop!!!
+      redUMM_old = redUMM;      
+      //calculate IC estimates
+      tmpU = Data.t() * redUMM;
+      if(opts.nonlinearity.value()=="tanh"){
+				//Matrix hyptanh;
+				//hyptanh = tanh(opts.nlconst1.value()*tmpU);
+				//redUMM = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
+				//sum(1-pow(hyptanh,2),1),redUMM))/samples;
+				gtmpU = tanh(opts.nlconst1.value()*tmpU);
+				redUMM = redUMM + lambda*(Id+(1-2*gtmpU.t()*tmpU))*redUMM;
+      }
+      if(opts.nonlinearity.value()=="gauss"){
+				gtmpU = pow(1+exp(-(opts.nlconst2.value()/2) * tmpU),-1);
+				redUMM = redUMM + lambda*(Id - (gtmpU.t()-tmpU.t())*tmpU)*redUMM;
+      }
+ 
+      //termination condition : angle between old and new < epsilon
+      minAbsSin = abs(1 - diag(abs(redUMM.t()*redUMM_old)).Minimum());
+      message("  Step no. " << itt_ctr << " change : " << minAbsSin << endl);
+      if(abs(minAbsSin) < opts.epsilon.value()){ break;}
+      
+      itt_ctr++;
+    } while(itt_ctr < opts.maxNumItt.value());
+    
+    if(itt_ctr>=opts.maxNumItt.value()){
+      cerr << "  No convergence after " << itt_ctr <<" steps "<<endl;
+    } else {
+      message("  Convergence after " << itt_ctr <<" steps " << endl << endl);
+      no_convergence = false;
+      {Matrix temp(melodat.get_dewhite() * redUMM);
+       melodat.set_mix(temp);}
+      {Matrix temp(redUMM.t()*melodat.get_white());
+      melodat.set_unmix(temp);}
+    } 
+  }
+  
+  void MelodicICA::ica_jade(const Matrix &Data){ 
+    int dim_sym = (int) (dim*(dim+1))/2;  
+    int num_CM = dim_sym;
+    Matrix CM;
+    Matrix R; R = IdentityMatrix(dim);
+    Matrix Qij; Qij = zeros(dim);
+    Matrix Xim;
+    Matrix Xjm;
+    Matrix scale; scale = ones(dim,1)/samples;
+
+    for (int im =1; im <= dim; im++){
+      Xim = Data.Row(im);
+			write_ascii_matrix("Xim",Data.Row(1));
+      //Qij = SP((scale * pow(Xim,2)),Data) * Data.t();//- R - 2*R.Column(im)*R.Column(im).t();
+      Qij = (pow(Xim,2)) * Data.t();//- R - 2*R.Column(im)*R.Column(im).t();
+      if(im==1){CM = Qij; write_ascii_matrix("CM",CM);exit(2);}else{CM |= Qij;}
+      for (int jm = 1; jm < im; jm++){
+				Xjm = Data.Row(jm);
+				Qij = (SP((scale * SP(Xim,Xjm)),Data) * Data.t() - 
+					R.Column(im)*R.Column(jm).t() - R.Column(jm)*R.Column(im).t());
+				Qij *= sqrt(2);
+				CM  |= Qij;
+      }
+    }
+
+    write_ascii_matrix("CM",CM);
+    Matrix redUMM; redUMM = IdentityMatrix(dim);
+  
+    bool exitloop = false;
+    int ctr_itt = 0;
+    int ctr_updates = 0;
+    Matrix Givens; Givens = zeros(2,num_CM);
+    Matrix Givens_ip; Givens_ip = zeros(2);
+    Matrix Givens_ro; Givens_ro = zeros(2);
+    double det_on, det_off;
+    double cos_theta, sin_theta, theta;
+
+    while(!exitloop && ctr_itt <= opts.maxNumItt.value()){
+      ctr_itt++;
+      cout << "loop" <<endl;
+      for(int ctr_p = 1; ctr_p < dim; ctr_p++){
+				for(int ctr_q = ctr_p+1; ctr_q <= dim; ctr_q++){
+
+	  			for(int ctr_i = 0; ctr_i < num_CM; ctr_i++){
+	    			int Ip = ctr_p + ctr_i * dim;
+	    			int Iq = ctr_q + ctr_i * dim;
+	    			Givens(1,ctr_i + 1) = CM(ctr_p,Ip) - CM(ctr_q,Iq);
+	    			Givens(2,ctr_i + 1) = CM(ctr_p,Iq) - CM(ctr_q,Ip);
+	  			}
+	  
+	  			Givens_ip = Givens * Givens.t();
+	  			det_on = Givens_ip(1,1) - Givens_ip(2,2);
+	  			det_off = Givens_ip(1,2) + Givens_ip(2,1);
+	  			theta = 0.5 * atan2(det_off, det_on + sqrt(det_on*det_on + det_off*det_off));
+
+	  			cout << theta << endl;
+
+	  			if(abs(theta) > opts.epsilon.value()){
+	    			ctr_updates++;
+	    			message("  Step No. "<< ctr_itt << " change : " << theta << endl);
+			
+	    			//create Givens rotation matrix
+	    			cos_theta = cos(theta);
+	    			sin_theta = sin(theta);
+	    			Givens_ro(1,1) = cos_theta;
+	    			Givens_ro(1,2) = -sin_theta;
+	    			Givens_ro(2,1) = sin_theta;
+	    			Givens_ro(2,2) = cos_theta;
+
+	    			//update 2 entries of redUMM
+	    			{Matrix tmp_redUMM;
+	    				tmp_redUMM = redUMM.Column(ctr_p);
+	    				tmp_redUMM |= redUMM.Column(ctr_q);
+	    				tmp_redUMM *= Givens_ro;
+	    				redUMM.Column(ctr_p) = tmp_redUMM.Column(1);
+	    				redUMM.Column(ctr_q) = tmp_redUMM.Column(2);}
+
+	    			//update Cumulant matrix
+	    			{Matrix tmp_CM;
+	    				tmp_CM = CM.Row(ctr_p);
+	    				tmp_CM &= CM.Row(ctr_q);
+	    				tmp_CM = Givens_ro.t() * tmp_CM;
+	    				CM.Row(ctr_p) = tmp_CM.Row(1);
+	    				CM.Row(ctr_q) = tmp_CM.Row(2);}
+
+	    			//update Cumulant matrices
+	    			for(int ctr_i = 0; ctr_i < num_CM; ctr_i++){
+	      			int Ip = ctr_p + ctr_i * dim;
+	      			int Iq = ctr_q + ctr_i * dim;
+	      			CM.Column(Ip) = cos_theta*CM.Column(Ip)+sin_theta*CM.Column(Iq);
+	      			CM.Column(Iq) = cos_theta*CM.Column(Iq)-sin_theta*CM.Column(Ip);
+	    			}
+	  			}else{
+	    			exitloop = true;
+	  			}
+				}
+      }
+    }//while loop
+    if(ctr_itt > opts.maxNumItt.value()){
+    	cerr << "  No convergence after " << ctr_itt <<" steps "<<endl;
+    } else {
+     	message("  Convergence after " << ctr_itt <<" steps " << endl << endl);
+      no_convergence = false;
+      {Matrix temp(melodat.get_dewhite() * redUMM);
+       	melodat.set_mix(temp);}
+      {Matrix temp(redUMM.t()*melodat.get_white());
+       	melodat.set_unmix(temp);}
+    }
+  }
+
+  Matrix MelodicICA::sign(const Matrix &Inp){
+    Matrix Res = Inp;
+    Res = 1;
+    for(int ctr_i = 1; ctr_i <= Inp.Ncols(); ctr_i++){
+      for(int ctr_j = 1; ctr_j <= Inp.Nrows(); ctr_j++){
+				if(Inp(ctr_j,ctr_i)<0){Res(ctr_j,ctr_i)=-1;}
+      }
+    } 
+    return Res;
+  }
+
+  void MelodicICA::perf_ica(const Matrix &Data){ 
+    message("Starting ICA estimation using " << opts.approach.value() 
+	    << endl << endl);
+    dim = Data.Nrows();
+    samples = Data.Ncols();
+    no_convergence = true;
+    //switch to the chosen method   
+    if(opts.approach.value()==string("symm") ||
+       opts.approach.value()==string("tica") ||
+       opts.approach.value()==string("parafac") ||
+       opts.approach.value()==string("concat"))
+      ica_fastica_symm(Data);
+    if(opts.approach.value()==string("defl"))
+      ica_fastica_defl(Data);
+    if(opts.approach.value()==string("jade"))
+      ica_jade(Data);
+    if(opts.approach.value()==string("maxent"))
+      ica_maxent(Data);
+    
+    if(!no_convergence){//calculate the IC
+      Matrix temp(melodat.get_unmix()*melodat.get_Data());
+      //  Add the mean time course again  
+      //      temp += melodat.get_unmix()*melodat.get_meanC()*ones(1,temp.Ncols());
+
+      //re-normalise the decomposition to std(mix)=1
+      Matrix scales;
+      scales = stdev(melodat.get_mix());   
+
+      //cerr << " SCALES 1 " << scales << endl;
+      Matrix tmp, tmp2;
+      tmp = SP(melodat.get_mix(),ones(melodat.get_mix().Nrows(),1)*pow(scales,-1));
+      temp = SP(temp,scales.t()*ones(1,temp.Ncols()));
+
+      scales = scales.t();
+  
+      melodat.set_mix(tmp);
+
+      melodat.set_IC(temp);
+
+      melodat.set_ICstats(scales);
+      melodat.sort();
+
+	  //message("Calculating T- and S-modes " << endl);
+      melodat.set_TSmode();
+		
+    }
+  }
+}
+
+

melica.h

+/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
+              independent components
+    
+    melica.h - ICA estimation
+
+    Christian F. Beckmann, FMRIB Analysis Group
+    
+    Copyright (C) 1999-2013 University of Oxford */
+
+/*  CCOPYRIGHT */
+
+#ifndef __MELODICICA_h
+#define __MELODICICA_h
+#include "utils/log.h"
+#include "melpca.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+#include "newmatap.h"
+#include "newmatio.h"
+#include "melreport.h"
+#include "ggmix.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+  
+  class MelodicICA{
+    public:
+      MelodicICA(MelodicData &pmelodat, MelodicOptions &popts, Log &plogger, MelodicReport &preport):  
+				melodat(pmelodat),
+				opts(popts),
+				logger(plogger),
+				report(preport){} 
+      
+      void perf_ica(const Matrix &Data);
+      
+      bool no_convergence;
+
+    private:
+      MelodicData &melodat;
+      MelodicOptions &opts;
+      Log &logger;
+      MelodicReport &report;
+
+      int dim;
+      int samples; 
+      Matrix redUMM;
+
+      void ica_fastica_symm(const Matrix &Data);
+      void ica_fastica_defl(const Matrix &Data);
+      void ica_maxent(const Matrix &Data);
+      void ica_jade(const Matrix &Data);
+      //void tica();
+      //void tica_concat();
+      //void parafac();
+      Matrix randm(const int dim1, const int dim2);
+      void sort();
+      Matrix sign(const Matrix &Inp);
+
+  };   
+}
+
+#endif

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 "newmatap.h"
+#include "newmatio.h"
+#include "melpca.h"
+#include "melhlprfns.h"
+#include "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;
+  }
+
+}
+
+

melpca.h

+ /*  MELODIC - Multivariate exploratory linear optimized decomposition into 
+              independent components
+    
+    melpca.h - PCA and whitening 
+
+    Christian F. Beckmann, FMRIB Analysis Group
+    
+    Copyright (C) 1999-2013 University of Oxford */
+
+/*  CCOPYRIGHT */
+
+#ifndef __MELODICPCA_h
+#define __MELODICPCA_h
+#include "utils/log.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+#include "newmatap.h"
+#include "newmatio.h"
+
+using namespace Utilities;
+
+namespace Melodic{
+  
+  class MelodicReport;
+
+  class MelodicPCA{
+    public:
+      MelodicPCA(MelodicData &pmelodat, MelodicOptions &popts, Log &plogger, 
+				MelodicReport &preport):  
+					melodat(pmelodat),
+					opts(popts),
+					logger(plogger),
+					report(preport){} 
+      
+    	void perf_pca(Matrix& in, Matrix& weights);
+      inline void perf_pca(){
+				perf_pca(melodat.get_Data(),melodat.get_RXweight());
+			}
+      void perf_white();
+
+    private:
+      MelodicData &melodat;
+      MelodicOptions &opts;
+      Log &logger;
+      __attribute__((unused)) MelodicReport &report;
+
+      int pcadim();
+  };   
+}
+
+#endif