MELODIC 3.0 beta

meldata.cc

@@ -181,10 +181,10 @@ namespace Melodic{
     }
     else{
       order = opts.pca_dim.value();
-      std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
+      std_pca(alldat, RXweight, Corr, pcaE, pcaD);
       calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
     }  
-      
+   
     if(numfiles < 2){
       Data = alldat;
       Matrix tmp = Identity(Data.Nrows());
@@ -193,7 +193,7 @@ namespace Melodic{
     } else {
       for(int ctr = 0; ctr < numfiles; ctr++){
 	tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
-
+	
 	//  whiten (separate / joint) 
 	if(!opts.joined_whiten.value()){	  
       	  std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
@@ -248,6 +248,7 @@ namespace Melodic{
     char callRMstr[1000];
     ostrstream osc(callRMstr,1000);
     osc  << "rm " << string(Mean_fname) <<"*  " << '\0';
+
     system(callRMstr);
  
     if(!samesize(Mean,Mask)){
@@ -266,9 +267,24 @@ namespace Melodic{
     //seed the random number generator
     double tmptime = time(NULL);
     srand((unsigned int) tmptime);
-  }
-
 
+		//read in post-proc design matrices etc
+		if(opts.fn_Tdesign.value().length()>0)
+			Tdes = read_vest(opts.fn_Tdesign.value());
+		if(opts.fn_Sdesign.value().length()>0)
+			Sdes = read_vest(opts.fn_Sdesign.value());
+		if(opts.fn_Tcon.value().length()>0)
+			Tcon = read_vest(opts.fn_Tcon.value());
+		if(opts.fn_Scon.value().length()>0)
+			Scon = read_vest(opts.fn_Scon.value());
+		if(opts.fn_TconF.value().length()>0)
+			TconF = read_vest(opts.fn_TconF.value());
+		if(opts.fn_SconF.value().length()>0)
+			SconF = read_vest(opts.fn_SconF.value());
+		
+		Tdes = remmean(Tdes,1);
+		Sdes = remmean(Sdes,1);	
+  }
 
 // {{{ Save
   void MelodicData::save()
@@ -334,7 +350,7 @@ namespace Melodic{
     //Output mixMatrix
     if(mixMatrix.Storage()>0){
       saveascii(mixMatrix, opts.outputfname.value() + "_mix");
-      mixFFT=calc_FFT(mixMatrix, opts.logPower.value());
+      mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
       saveascii(mixFFT,opts.outputfname.value() + "_FTmix");      
     }
 
@@ -746,7 +762,7 @@ namespace Melodic{
       ICstats |= copeN;
     }
     
-    mixFFT=calc_FFT(mixMatrix, opts.logPower.value());
+    mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
     unmixMatrix = pinv(mixMatrix);
 
     //if(ICstats.Storage()>0){cout << "ICstats: " << ICstats.Nrows() <<"x" << ICstats.Ncols() << endl;}else{cout << "ICstats empty " <<endl;}

meldata.h

@@ -29,32 +29,31 @@ namespace Melodic{
 
       //constructor
       MelodicData(MelodicOptions &popts, Log &plogger):  
-	opts(popts),
-	logger(plogger)
-	{	
-	  after_mm = false;
-	  Resels = 0;
-	}  
+				opts(popts),logger(plogger)
+			{	
+	  		after_mm = false;
+	  		Resels = 0;
+			}  
  
       void save();
 
       Matrix process_file(string fname, int numfiles);
 
       inline void save4D(Matrix what, string fname){
-	 volume4D<float> tempVol;
-	 tempVol.setmatrix(what,Mask);
-	 save_volume4D(tempVol,logger.appendDir(fname),tempInfo);
-	 message("  " << logger.appendDir(fname) << endl);
+	 			volume4D<float> tempVol;
+	 			tempVol.setmatrix(what,Mask);
+	 			save_volume4D(tempVol,logger.appendDir(fname),tempInfo);
+	 			message("  " << logger.appendDir(fname) << endl);
       }
       
       inline void saveascii(Matrix what, string fname){
-	 write_ascii_matrix(logger.appendDir(fname),what);   
-	 message("  " << logger.appendDir(fname) << endl);   
+	 			write_ascii_matrix(logger.appendDir(fname),what);   
+	 			message("  " << logger.appendDir(fname) << endl);   
       }
  
       inline void savebinary(Matrix what, string fname){
-         write_binary_matrix(what,logger.appendDir(fname));  
-	 message("  " << logger.appendDir(fname) << endl);    
+      	write_binary_matrix(what,logger.appendDir(fname));  
+	 			message("  " << logger.appendDir(fname) << endl);    
       }
 
       int  remove_components();
@@ -78,20 +77,20 @@ namespace Melodic{
       inline Matrix& get_Smodes(int what) {return Smodes.at(what);}
       inline void add_Smodes(Matrix& Arg) {Smodes.push_back(Arg);}      
       inline void save_Smodes(){
-	Matrix tmp = Smodes.at(0); 
-	for(unsigned int ctr = 1; ctr < Smodes.size(); ctr++)
-	  tmp |= Smodes.at(ctr);
-	saveascii(tmp,opts.outputfname.value() + "_Smodes");
+				Matrix tmp = Smodes.at(0); 
+				for(unsigned int ctr = 1; ctr < Smodes.size(); ctr++)
+	  			tmp |= Smodes.at(ctr);
+				  saveascii(tmp,opts.outputfname.value() + "_Smodes");
       }
 
       inline vector<Matrix>& get_Tmodes() {return Tmodes;}
       inline Matrix& get_Tmodes(int what) {return Tmodes.at(what);}
       inline void add_Tmodes(Matrix& Arg) {Tmodes.push_back(Arg);}
       inline void save_Tmodes(){
-	Matrix tmp = Tmodes.at(0); 
-	for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++)
-	  tmp |= Tmodes.at(ctr);
-	saveascii(tmp,opts.outputfname.value() + "_Tmodes");
+				Matrix tmp = Tmodes.at(0); 
+				for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++)
+	  			tmp |= Tmodes.at(ctr);
+				saveascii(tmp,opts.outputfname.value() + "_Tmodes");
       }
 
       void set_TSmode();
@@ -111,12 +110,12 @@ namespace Melodic{
       inline Matrix& get_mix() {return mixMatrix;}
 
       inline void set_mix(Matrix& Arg) {
-	mixMatrix = Arg;
-	if (Tmodes.size() < 1)
-	  for (int ctr = 1; ctr <= Arg.Ncols(); ctr++){
-	    Matrix tmp = Arg.Column(ctr);
-	    add_Tmodes(tmp);
-	  }
+				mixMatrix = Arg;
+				if (Tmodes.size() < 1)
+	  			for (int ctr = 1; ctr <= Arg.Ncols(); ctr++){
+	    			Matrix tmp = Arg.Column(ctr);
+	    			add_Tmodes(tmp);
+	  			}
       }
 
       Matrix expand_mix(); 
@@ -167,24 +166,24 @@ namespace Melodic{
       inline void set_after_mm(bool val) {after_mm = val;}
 
       inline void flipres(int num){
-	IC.Row(num) = -IC.Row(num);
-	mixMatrix.Column(num) = -mixMatrix.Column(num);
-	mixFFT=calc_FFT(mixMatrix);
-	unmixMatrix = pinv(mixMatrix);
+				IC.Row(num) = -IC.Row(num);
+				mixMatrix.Column(num) = -mixMatrix.Column(num);
+				mixFFT=calc_FFT(mixMatrix);
+				unmixMatrix = pinv(mixMatrix);
 
         if(ICstats.Storage()>0&&ICstats.Ncols()>3){
-	  double tmp;
-	  tmp = ICstats(num,3);
-	  ICstats(num,3) = -1.0*ICstats(num,4);
-	  ICstats(num,4) = -1.0*tmp;
-	}
+	  			double tmp;
+	  			tmp = ICstats(num,3);
+	  			ICstats(num,3) = -1.0*ICstats(num,4);
+	  			ICstats(num,4) = -1.0*tmp;
+				}
       }
       
       void sort();
 
       volumeinfo tempInfo;
-
       vector<Matrix> DWM, WM;
+			basicGLM glmT, glmS;
 
     private:
       MelodicOptions &opts;     
@@ -204,6 +203,7 @@ namespace Melodic{
       Matrix mixFFT;
       Matrix IC;
       Matrix ICstats;
+			Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
 
       vector<Matrix> Tmodes;
       vector<Matrix> Smodes;

melhlprfns.cc

@@ -12,6 +12,8 @@
 #include "melhlprfns.h"
 #include "libprob.h"
 #include "miscmaths/miscprob.h"
+#include "miscmaths/t2z.h"
+#include "miscmaths/f2z.h"
 
 namespace Melodic{
 
@@ -781,4 +783,52 @@ namespace Melodic{
     return res;
   }  //Matrix gen_arCorr
 
+	Matrix basicGLM::olsfit(const Matrix& data, const Matrix& design, 
+		const Matrix& contrasts, int DOFadjust)
+	{
+		beta = zeros(design.Ncols(),1); 
+		residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1); 
+		t = zeros(1); z = zeros(1); p=-1.0*ones(1);
+		dof = (int)-1; cbeta = -1.0*ones(1); 
+
+		if(data.Nrows()==design.Nrows()){
+			Matrix dat = remmean(data,1);
+			Matrix pinvdes = pinv(design);
+			
+			beta = pinvdes * dat;
+			residu = dat - design*beta;
+			Matrix R = Identity(design.Nrows()) - design * pinvdes;
+			Matrix R2 = R*R;
+			float tR = R.Trace();
+			sigsq = sum(SP(residu,residu))/tR;
+			
+			dof = (int)(tR*tR/R2.Trace() - DOFadjust);
+			
+			float fact = float(design.Nrows() - 1 - design.Ncols()) / design.Ncols();
+			f_fmf = SP(var(design*beta),pow(var(residu),-1))* fact;
+			pf_fmf = f_fmf.Row(1); pf_fmf &= pf_fmf;
+			for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
+				pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),
+				int(design.Nrows() -1 -design.Ncols()),f_fmf.Column(ctr1).AsScalar());
+			pf_fmf.Row(2) = pf_fmf.Row(1) * pf_fmf.Ncols();
+				
+			if(contrasts.Storage()>0 && contrasts.Nrows()==beta.Ncols()){
+				cbeta = contrasts.t()*beta;
+				Matrix tmp = contrasts.t()*pinvdes*pinvdes.t()*contrasts;
+				varcb = diag(tmp)*sigsq;
+				t = SP(cbeta,pow(varcb,-0.5));
+				z = t; p=t; 
+				for(int ctr1=1;ctr1<=t.Ncols();ctr1++){
+					ColumnVector tmp = t.Column(ctr1);
+					T2z::ComputeZStats(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
+					z.Column(ctr1) << tmp;
+					T2z::ComputePs(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
+					p.Column(ctr1) << exp(tmp);
+				}
+			}
+		}	
+		return beta;
+	}
+
+
 }

melhlprfns.h

@@ -70,7 +70,44 @@ namespace Melodic{
   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(){}
+		
+			Matrix olsfit(const Matrix& data, const Matrix& design, 
+				const Matrix& contrasts, int DOFadjust = 0);
+
+			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_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

@@ -499,7 +499,6 @@ write_ascii_matrix("Xim",Data.Row(1));
       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());
@@ -511,18 +510,20 @@ write_ascii_matrix("Xim",Data.Row(1));
       scales = stdev(melodat.get_mix());   
 
       //      cerr << " SCALES 1 " << scales << endl;
-      Matrix tmp;
+      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();
 
       melodat.set_TSmode();
+		
     }
   }
 }

melodic.cc

@@ -56,7 +56,6 @@ void mmonly(Log& logger, MelodicOptions& opts,
 
 int main(int argc, char *argv[])
 {
-
   try{
     // Setup logging:
     Log& logger  =   LogSingleton::getInstance();
@@ -67,16 +66,16 @@ int main(int argc, char *argv[])
 
     //set up data object
     MelodicData melodat(opts,logger);
-      
+    //set up report object     
     MelodicReport report(melodat,opts,logger);
    
     if (opts.filtermode || opts.filtermix.value().length()>0){
       if(opts.filtermode){ // just filter out some noise from a previous run
-	melodat.setup();
-	melodat.remove_components();
+    		melodat.setup();
+				melodat.remove_components();
       }
       else
-	mmonly(logger,opts,melodat,report);
+				mmonly(logger,opts,melodat,report);
     } 
     else
     {  // standard PICA now
@@ -84,105 +83,100 @@ int main(int argc, char *argv[])
       bool no_conv;
       bool leaveloop = false;
 
-      //      melodat.setup2();
       melodat.setup();
 
       do{
-	//do PCA pre-processing
-	MelodicPCA pcaobj(melodat,opts,logger,report);
-	pcaobj.perf_pca();
-	pcaobj.perf_white();
-
-	//do ICA
-	MelodicICA icaobj(melodat,opts,logger,report);
-	icaobj.perf_ica(melodat.get_white()*melodat.get_Data());
+				//do PCA pre-processing
+				MelodicPCA pcaobj(melodat,opts,logger,report);
+				pcaobj.perf_pca();
+				pcaobj.perf_white();
+
+				//do ICA
+				MelodicICA icaobj(melodat,opts,logger,report);
+				icaobj.perf_ica(melodat.get_white()*melodat.get_Data());
     
-	no_conv = icaobj.no_convergence;
-
-	opts.maxNumItt.set_T(500);
-	if((opts.approach.value()=="symm")&&
-	   (retry > std::min(opts.retrystep,3))){
-	  if(no_conv){
-	    retry++;
-	    opts.approach.set_T("defl");
-	    message(endl << "Restarting MELODIC using deflation approach" 
-		    << endl << endl);
-	  }
-	  else{
-	    leaveloop = true;
-	  }
-	}
-	else{
-	  if(no_conv){
-	    retry++;
-	    if(opts.pca_dim.value()-retry*opts.retrystep > 
-	       0.1*melodat.data_dim()){
-	      opts.pca_dim.set_T(opts.pca_dim.value()-retry*opts.retrystep);
-	    }
-	    else{
-	      if(opts.pca_dim.value()-retry*opts.retrystep <  melodat.data_dim()){
-		opts.pca_dim.set_T(opts.pca_dim.value()+retry*opts.retrystep);
-	      }else{
-		leaveloop = true; //stupid, but break does not compile 
-		                  //on all platforms
-	      }
-	    }
-	    if(!leaveloop){
+				no_conv = icaobj.no_convergence;
+
+				opts.maxNumItt.set_T(500);
+				if((opts.approach.value()=="symm")&&
+	   		(retry > std::min(opts.retrystep,3))){
+	  			if(no_conv){
+	    			retry++;
+	    			opts.approach.set_T("defl");
+	    			message(endl << "Restarting MELODIC using deflation approach" 
+		    			<< endl << endl);
+	  			}
+	  			else{
+	    			leaveloop = true;
+	  			}
+				}
+				else{
+	  			if(no_conv){
+	    			retry++;
+	    			if(opts.pca_dim.value()-retry*opts.retrystep > 
+	       			0.1*melodat.data_dim()){
+	      				opts.pca_dim.set_T(opts.pca_dim.value()-retry*opts.retrystep);
+	    				}
+	    			else{
+	      			if(opts.pca_dim.value()-retry*opts.retrystep <  melodat.data_dim()){
+								opts.pca_dim.set_T(opts.pca_dim.value()+retry*opts.retrystep);
+	      			}else{
+								leaveloop = true; //stupid, but break does not compile 
+		                  						//on all platforms
+	      			}
+	    			}
+	    			if(!leaveloop){
 	      message(endl << "Restarting MELODIC using -d " 
 		      << opts.pca_dim.value() 
 		      << endl << endl);
 	    }
-	  }
-	}
+	  			}
+				}
       } while (no_conv && retry<opts.maxRestart.value() && !leaveloop);	
      
       if(!no_conv){
-	//save raw IC results
-	melodat.save();
-
-	Matrix pmaps;//(melodat.get_IC());
-	Matrix mmres;
-
-	if(opts.perf_mm.value())
-	  mmres = mmall(logger,opts,melodat,report,pmaps);
-	else{
-	  if( bool(opts.genreport.value()) ){
-	    message(endl 
-		    << "Creating web report in " << report.getDir() 
-		    << " " << endl);
-	    for(int ctr=1; ctr<= melodat.get_IC().Nrows(); ctr++){
-	      string prefix = "IC_"+num2str(ctr);
-	      message("  " << ctr);
-	      report.IC_simplerep(prefix,ctr,melodat.get_IC().Nrows());
-	    }
-
-	    
-	    
-	    message(endl << endl <<
-		    " To view the output report point your web browser at " <<
-		    report.getDir() + "/00index.html" << endl<< endl); 
-	  }
-	}		 
-
-	if( bool(opts.genreport.value()) ){
-	  report.analysistxt();
-	  report.PPCA_rep();
-	}
-
-	message("finished!" << endl << endl);
-      } else { 
-	message(endl <<"No convergence -- giving up " << endl << endl);
+				//save raw IC results
+				melodat.save();
+				Matrix pmaps;//(melodat.get_IC());
+				Matrix mmres;
+
+				if(opts.perf_mm.value())
+	  			mmres = mmall(logger,opts,melodat,report,pmaps);
+				else{
+	  			if( bool(opts.genreport.value()) ){
+	    			message(endl 
+		    			<< "Creating web report in " << report.getDir() 
+		    			<< " " << endl);
+	    			for(int ctr=1; ctr<= melodat.get_IC().Nrows(); ctr++){
+	      			string prefix = "IC_"+num2str(ctr);
+	      			message("  " << ctr);
+	      			report.IC_simplerep(prefix,ctr,melodat.get_IC().Nrows());
+	    			}
+
+	    			message(endl << endl <<
+		    			" To view the output report point your web browser at " <<
+		    			report.getDir() + "/00index.html" << endl<< endl); 
+	  			}
+				}		 
+
+				if( bool(opts.genreport.value()) ){
+	  			report.analysistxt();
+	  			report.PPCA_rep();
+				}
+
+				message("finished!" << endl << endl);
+      } 
+			else { 
+				message(endl <<"No convergence -- giving up " << endl << endl);
       }	     
     }
   }
-  catch(Exception e) 
-    {
+  catch(Exception e) {
       cerr << endl << e.what() << endl;
-    }
-  catch(X_OptionError& e) 
-    {
+  }
+  catch(X_OptionError& e) {
       cerr << endl << e.what() << endl;
-    }
+  }
 
   return 0;
 }
@@ -271,7 +265,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
   melodat.set_mean(Mean);
   melodat.set_IC(ICs);
   melodat.set_mix(mixMatrix);
-  fmixMatrix = calc_FFT(mixMatrix, opts.logPower.value());
+  fmixMatrix = calc_FFT(melodat.expand_mix(), opts.logPower.value());
   melodat.set_fmix(fmixMatrix);
   fmixMatrix = pinv(mixMatrix);
   melodat.set_unmix(fmixMatrix);

melodic.h

@@ -34,7 +34,7 @@
 
 namespace Melodic{
 
-const string version = "3.0 alpha 1";  
+const string version = "3.0 beta";  
 
 }
 

meloptions.cc

@@ -64,11 +64,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
     } 
   if(! options.check_compulsory_arguments())
     {
-      print_version();
-      cout << endl;
-      cout << "Usage: melodic <options> -i <filename>" << endl <<
-	"Please specify input file name correctly or try melodic --help" 
-	   << endl << endl;
+			print_usage(argc, argv);
       exit(2);
     }     
 
@@ -173,8 +169,8 @@ MelodicOptions* MelodicOptions::gopt = NULL;
     while (!fs.eof()) {
       getline(fs,cline);
       if(cline.length()>0)
-	tmpfnames.push_back(cline);
-    }
+				tmpfnames.push_back(cline);
+    }	
     fs.close();
     inputfname.set_T(tmpfnames);
   }
@@ -211,7 +207,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
 
 void MelodicOptions::print_usage(int argc, char *argv[])
 {
-  print_version();
+  print_copyright();
   options.usage();
   /* cout << "Usage: " << argv[0] << " ";
 
@@ -221,14 +217,14 @@ void MelodicOptions::print_usage(int argc, char *argv[])
 
 void MelodicOptions::print_version()
 {
-  cout << "MELODIC Version " << version;
+  cout << endl <<"MELODIC Version " << version;
   cout.flush();
 }
 
 void MelodicOptions::print_copyright()
 {
   print_version();
-  cout << endl << "Copyright (C) 1999-2002 University of Oxford " << endl;
+  cout << endl << "Copyright (C) 1999-2007 University of Oxford " << endl;
   cout.flush();
 }
 

meloptions.h

@@ -80,6 +80,13 @@ class MelodicOptions {
   Option<float>  tr;
   Option<bool>   logPower;
 
+	Option<string> fn_Tdesign;
+	Option<string> fn_Tcon;
+	Option<string> fn_TconF;
+	Option<string> fn_Sdesign;
+	Option<string> fn_Scon;	
+	Option<string> fn_SconF;	
+	
   Option<bool>   output_all;
   Option<bool>   output_unmix;
   Option<bool>   output_MMstats;
@@ -103,7 +110,6 @@ class MelodicOptions {
   Option<float> nlconst2;
   Option<float> smooth_probmap;
 
-
   Option<bool> remove_meanvol;
   Option<bool> remove_endslices;
   Option<bool> rescale_nht;
@@ -186,7 +192,7 @@ class MelodicOptions {
    varnorm(string("--vn,--varnorm"), true,
 	   string("switch off variance normalisation"), 
 	   false, no_argument),
-   pbsc(string("--pbsc"), true,
+   pbsc(string("--pbsc"), false,
 	   string("switch off conversion to percent BOLD signal change"), 
 	   false, no_argument),
    pspec(string("--pspec"), false,
@@ -207,7 +213,7 @@ class MelodicOptions {
    maxRestart(string("--maxrestart"),  6,
 	   string("maximum number of restarts\n"), 
 	   false, requires_argument),
-   rank1interval(string("--rank1interval"),  30,
+   rank1interval(string("--rank1interval"),  5,
 	   string("number of iterations between rank-1 approximation (TICA)\n"), 
 		 false, requires_argument,false),
     mmthresh(string("--mmthresh"), string("0.5"),
@@ -237,6 +243,24 @@ class MelodicOptions {
    logPower(string("--logPower"),  false,
 	   string("calculate log of power for frequency spectrum\n"), 
 	    false, no_argument),
+	 fn_Tdesign(string("--Tdes"), string(""),
+	   string("design matrix across time-domain"),
+		 false, requires_argument),
+	 fn_Tcon(string("--Tcon"), string(""),
+		 string("t-contrast matrix across time-domain"),
+		 false, requires_argument),
+   fn_TconF(string("--Tconf"), string(""),
+		 string("F-contrast matrix across time-domain"),
+		 false, requires_argument),
+   fn_Sdesign(string("--Sdes"), string(""),
+		 string("design matrix across subject-domain"),
+		 false, requires_argument),
+	 fn_Scon(string("--Scon"), string(""),
+		 string("t-contrast matrix across subject-domain"),
+		 false, requires_argument),	
+	 fn_SconF(string("--Sconf"), string(""),
+		 string("F-contrast matrix across subject-domain"),
+		 false, requires_argument),	
    output_all(string("--Oall"),  false,
 	   string("output everything"), 
 	   false, no_argument),
@@ -356,6 +380,12 @@ class MelodicOptions {
 	    options.add(genreport);
 	    options.add(tr);
 	    options.add(logPower);
+			options.add(fn_Tdesign);
+			options.add(fn_Tcon);
+			options.add(fn_TconF);
+			options.add(fn_Sdesign);
+			options.add(fn_Scon);
+			options.add(fn_SconF);
 	    options.add(output_all);
 	    options.add(output_unmix);
 	    options.add(output_MMstats);

test.cc

@@ -37,18 +37,23 @@ int do_work(int argc, char* argv[])
 {
  
 	Matrix mat;
-	mat = normrnd(300,(int)num.value())+2;
 	
 	cout << "BLAH " << num.value() << endl;
-  	for (int i=1; i <= (int)num.value();i++){
+	mat=normrnd(300,1);
+	miscplot::Timeseries(mat.t(),string("test0.png"),string("TEST"));
+	
+  for (int i=1; i <= (int)num.value();i++){
 		cout << "Processing " << i << endl;
-	  	miscplot newplot;
-		ColumnVector col;
-       	col = mat.Column(i);
-      	newplot.add_bpdata(col);
-     	newplot.boxplot(string("test")+num2str(i)+string(".png"),string("TEST"));
+	  miscplot newplot;
+		newplot.GDCglobals_set();
+		mat = normrnd(300,3)+2;
+		Matrix col;	
+    col = mat;
+    newplot.add_bpdata(col);
+    //	newplot.add_bpdata(col);
+    newplot.boxplot(string("test")+num2str(i)+string(".png"),string("TEST"));
 	}	
-  	return 0;
+  return 0;
 }
 
 ////////////////////////////////////////////////////////////////////////////