MIGP v 0.2

meldata.cc

@@ -150,76 +150,83 @@ namespace Melodic{
     return tmp;
   }
 
-  void MelodicData::set_TSmode()
+  void MelodicData::set_TSmode_depr()
   {
-   	dbgmsg(string("START: set_TSmode"));	
-	
-	Matrix tmp, tmpT, tmpS, tmpT2, tmpS2, tmpT3;
-	
-    tmp = expand_dimred(mixMatrix);
-    tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
-    tmpS = ones(numfiles, tmp.Ncols());
-
-	outMsize("tmp",tmp);
-	outMsize("tmpT",tmpT);
-	outMsize("tmpS",tmpS);
-
-  	dbgmsg(string("   approach ") << opts.approach.value() << endl);	
-
-	if(opts.approach.value()==string("tica")){
-      message("Calculating T- and S-modes " << endl);
-      explained_var = krfact(tmp,tmpT,tmpS);
-      Tmodes.clear(); Smodes.clear();
-      for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
-	    tmpT3 << reshape(tmp.Column(ctr),tmpT.Nrows(),numfiles);
-		outMsize("tmpT3", tmpT3);
-      	tmpT2 << tmpT.Column(ctr);
-      	tmpS2 << tmpS.Column(ctr);
-		tmpT3 << SP(tmpT3,pow(ones(tmpT3.Nrows(),1)*tmpS2.t(),-1));
-		if(numfiles>1)
-	      tmpT2 |= tmpT3;
-		if(mean(tmpS2,1).AsScalar()<0){
-		  tmpT2*=-1.0;
-		  tmpS2*=-1.0;
+		Matrix tmp, tmpT, tmpS, tmpT2, tmpS2, tmpT3;
+
+	    tmp = expand_dimred(mixMatrix);
+	    tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
+	    tmpS = ones(numfiles, tmp.Ncols());
+
+		outMsize("tmp",tmp);
+		outMsize("tmpT",tmpT);
+		outMsize("tmpS",tmpS);
+
+	  	dbgmsg(string("   approach ") << opts.approach.value() << endl);	
+
+		if(opts.approach.value()==string("tica")){
+	      message("Calculating T- and S-modes " << endl);
+	      explained_var = krfact(tmp,tmpT,tmpS);
+	      Tmodes.clear(); Smodes.clear();
+	      for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
+		    tmpT3 << reshape(tmp.Column(ctr),tmpT.Nrows(),numfiles);
+			outMsize("tmpT3", tmpT3);
+	      	tmpT2 << tmpT.Column(ctr);
+	      	tmpS2 << tmpS.Column(ctr);
+			tmpT3 << SP(tmpT3,pow(ones(tmpT3.Nrows(),1)*tmpS2.t(),-1));
+			if(numfiles>1)
+		      tmpT2 |= tmpT3;
+			if(mean(tmpS2,1).AsScalar()<0){
+			  tmpT2*=-1.0;
+			  tmpS2*=-1.0;
+			}
+	      	add_Tmodes(tmpT2);
+	      	add_Smodes(tmpS2);
+		  }
 		}
-      	add_Tmodes(tmpT2);
-      	add_Smodes(tmpS2);
-	  }
-	}
-	else{
-		Tmodes.clear();
-		Smodes.clear();
-		for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
-			tmpT3 << tmp.Column(ctr);
-			add_Tmodes(tmpT3);
+		else{
+			Tmodes.clear();
+			Smodes.clear();
+			for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
+				tmpT3 << tmp.Column(ctr);
+				add_Tmodes(tmpT3);
+			}
 		}
-	}
-	
-    if(opts.approach.value()!=string("concat")){
-	  //add GLM OLS fit
-	  dbgmsg(string(" GLM fitting ") << endl);
-
-	  if(Tdes.Storage()){
-	    Matrix alltcs = Tmodes.at(0).Column(1);
-	    for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
-		  alltcs|=Tmodes.at(ctr).Column(1);
-	    if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols()))
-		  glmT.olsfit(alltcs,Tdes,Tcon);
-	  }
-	  if(Sdes.Storage()){
-	    Matrix alltcs = Smodes.at(0);
-	    for(int ctr=1; ctr < (int)Smodes.size();ctr++)
-	  	  alltcs|=Smodes.at(ctr);
-	    if((alltcs.Nrows()==Sdes.Nrows())&&(Sdes.Nrows()>Sdes.Ncols()&&alltcs.Nrows()>2))
-		  glmS.olsfit(alltcs,Sdes,Scon);
-	  }
-		
-    }
-  //    else{
-//		dbgmsg(string(" Bypassing krfac ") << endl);
-//        add_Tmodes(tmp);
-//		add_Smodes(tmpS);
-//      }
+
+	    if(opts.approach.value()!=string("concat")){
+		  //add GLM OLS fit
+		  dbgmsg(string(" GLM fitting ") << endl);
+
+		  if(Tdes.Storage()){
+		    Matrix alltcs = Tmodes.at(0).Column(1);
+		    for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
+			  alltcs|=Tmodes.at(ctr).Column(1);
+		    if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols()))
+			  glmT.olsfit(alltcs,Tdes,Tcon);
+		  }
+		  if(Sdes.Storage()){
+		    Matrix alltcs = Smodes.at(0);
+		    for(int ctr=1; ctr < (int)Smodes.size();ctr++)
+		  	  alltcs|=Smodes.at(ctr);
+		    if((alltcs.Nrows()==Sdes.Nrows())&&(Sdes.Nrows()>Sdes.Ncols()&&alltcs.Nrows()>2))
+			  glmS.olsfit(alltcs,Sdes,Scon);
+		  }
+
+	    }
+	  //    else{
+	//		dbgmsg(string(" Bypassing krfac ") << endl);
+	//        add_Tmodes(tmp);
+	//		add_Smodes(tmpS);
+	//      }
+  }
+
+  void MelodicData::set_TSmode()
+  {
+   	dbgmsg(string("START: set_TSmode"));	
+	if(opts.dr.value())
+		dual_regression();
+	else
+		set_TSmode_depr();
 	
 	dbgmsg(string("END: set_TSmode"));	
   }
@@ -393,8 +400,15 @@ namespace Melodic{
   void MelodicData::setup_migp()
   {
 	dbgmsg("starting MIGP");
-    //permute input vector if desired	std::random_shuffle ( opts.inputfname.value().begin(), opts.inputfname.value().end() );
+	
+	std::vector<int> myctr;
+	for (int i=0; i< numfiles ; ++i) myctr.push_back(i); 
 
+	if(opts.migp_shuffle.value()){
+		message("Randomising input file order" << endl);
+		std::random_shuffle ( myctr.begin(), myctr.end() );
+	}
+	
 	Matrix tmpData;
 	bool tmpvarnorm = opts.varnorm.value();
 
@@ -403,7 +417,11 @@ namespace Melodic{
 	}
 		
 	for(int ctr = 0; ctr < numfiles; ctr++){
-		tmpData = process_file(opts.inputfname.value().at(ctr), numfiles) / numfiles;
+		tmpData = process_file(opts.inputfname.value().at(myctr.at(ctr)), numfiles) / numfiles;
+		
+		if (opts.migpN.value()==0){
+			opts.migpN.set_T(2*tmpData.Nrows()-1);
+		}
 		if(opts.debug.value())
 			save4D(tmpData,string("preproc_dat") + num2str(ctr+1));
 			
@@ -412,6 +430,7 @@ namespace Melodic{
 		else
   			Data &= tmpData;
 
+		outMsize("Data", Data);
 		//reduce dim down to manageable level
 		if(Data.Nrows() > opts.migpN.value()){
 			message("  Reducing data matrix to a  " << opt.migpN.value() << " dimensional subspace " << endl);
@@ -421,6 +440,8 @@ namespace Melodic{
 		    pcaE = pcaE.Columns(pcaE.Ncols()-opts.migpN.value()+1,pcaE.Ncols());
 		    Data = pcaE.t() * Data;	
 		}
+		outMsize("Data", Data);
+		
     }
 
   	//update mask
@@ -1059,7 +1080,7 @@ namespace Melodic{
     unmixMatrix = pinv(mixMatrix);
   }
 
-  void MelodicData::reregress()
+  void MelodicData::dual_regression()
   {
 	if((numfiles > 1)){
 	   

meldata.h

@@ -103,6 +103,7 @@ namespace Melodic{
 			}
       }
 
+      void set_TSmode_depr();	
       void set_TSmode();
 
       inline Matrix& get_param() {return param;} 
@@ -206,7 +207,7 @@ namespace Melodic{
       }
       
       void sort();
-	  void reregress();
+	  void dual_regression();
 
       vector<Matrix> DWM, WM;
 			basicGLM glmT, glmS;

meloptions.h

@@ -56,6 +56,8 @@ class MelodicOptions {
 	Option<bool>   dr_pca;
 	Option<bool>   migp;
 	Option<int>    migpN;
+	Option<bool>   migp_shuffle;
+	Option<bool>   dr;
 	Option<float>  vn_level;
   	Option<int>    numICs;
   	Option<string> approach;
@@ -199,9 +201,15 @@ class MelodicOptions {
    migp(string("--migp"), false,
 	   string("switch on MIGP data reduction"),
 	   false, no_argument, false),	
-   migpN(string("--migpN"), 1000,
-		string("Number of internal Eigenmaps"),
-		false, requires_argument, false),
+   migpN(string("--migpN"), 0,
+	   string("Number of internal Eigenmaps"),
+	   false, requires_argument, false),
+   migp_shuffle(string("--migp_order"), true,
+	   string("Randomise MIGP file order (default: TRUE)"),
+	   false, no_argument, false),
+   dr(string("--dr"), false,
+	   string("Dual Regression (default: TRUE)"),
+	   false, no_argument, false),	
    vn_level(string("--vn_level"), float(2.3),
 	   string("variance nomalisation threshold level (Z> value is ignored)"), 
 	   false, requires_argument, false),
@@ -408,6 +416,8 @@ class MelodicOptions {
 		options.add(dr_pca);
 		options.add(migp);
 		options.add(migpN);
+		options.add(migp_shuffle);
+		options.add(dr);
 	    options.add(vn_level);
 	    options.add(numICs);
 	    options.add(approach);