MIGP with fixed TS plotting and Dual Regression

meldata.cc

@@ -220,15 +220,64 @@ namespace Melodic{
 	//      }
   }
 
+void MelodicData::dual_regression()
+  {
+	dbgmsg(string("START: dual_regression") << endl);	
+
+	Tmodes.clear();
+	Smodes.clear();
+
+	bool tmpvarnorm = opts.varnorm.value();
+	opts.varnorm.set_T(false);
+
+	Log drO;
+
+	if(opts.dr_out.value())
+		drO.makeDir(logger.appendDir("dr"),"dr.log");
+
+	Matrix tmpData, alltcs, tmp, pinvdes;
+	for(int ctr = 0; ctr < numfiles; ctr++){
+		tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
+		outMsize("IC",IC);
+		//may want to remove the spatial means first
+
+		
+		tmp = IC*IC.t();   		outMsize("tmp",tmp);
+		pinvdes = tmp.i()*IC;		  
+		tmp = tmpData*pinvdes.t();      outMsize(string("tmp ")+opts.inputfname.value().at(ctr),tmp);
+		if(alltcs.Storage()==0)
+			alltcs=tmp;
+		else
+			alltcs&=tmp;
+			
+		// output DR
+		if(opts.dr_out.value()){
+			write_ascii_matrix(drO.appendDir("dr_stage1_subject"+num2str(ctr)+".txt"),tmp);
+			pinvdes = tmp.t()*tmp;
+			pinvdes = pinvdes.i();
+			tmp = pinvdes*tmp.t()*tmpData;
+			save4D(tmp,string("dr/dr_stage2_subject"+num2str(ctr)));
+		}
+    }
+
+	for(int ctr = 1; ctr <= alltcs.Ncols(); ctr++){
+		tmp << alltcs.Column(ctr);
+		add_Tmodes(tmp);
+	}
+
+    opts.varnorm.set_T(tmpvarnorm);	
+	dbgmsg(string("END: dual_regression") << endl);	
+  }
+
   void MelodicData::set_TSmode()
   {
-   	dbgmsg(string("START: set_TSmode"));	
+   	dbgmsg(string("START: set_TSmode")<< endl);	
 	if(opts.dr.value())
 		dual_regression();
 	else
 		set_TSmode_depr();
 	
-	dbgmsg(string("END: set_TSmode"));	
+	dbgmsg(string("END: set_TSmode")<< endl);	
   }
 
   void MelodicData::setup_classic()
@@ -432,7 +481,7 @@ namespace Melodic{
 
 		outMsize("Data", Data);
 		//reduce dim down to manageable level
-		if(Data.Nrows() > opts.migpN.value()){
+		if(Data.Nrows() > opts.migp_factor.value()*opts.migpN.value()){
 			message("  Reducing data matrix to a  " << opt.migpN.value() << " dimensional subspace " << endl);
 			Matrix pcaE, Corr;
 			RowVector pcaD;
@@ -1080,13 +1129,6 @@ namespace Melodic{
     unmixMatrix = pinv(mixMatrix);
   }
 
-  void MelodicData::dual_regression()
-  {
-	if((numfiles > 1)){
-	   
-	}
-  }
-  
   void MelodicData::status(const string &txt)
   {
     cout << "MelodicData Object " << txt << endl;

meloptions.h

@@ -57,7 +57,9 @@ class MelodicOptions {
 	Option<bool>   migp;
 	Option<int>    migpN;
 	Option<bool>   migp_shuffle;
+	Option<int>	   migp_factor;
 	Option<bool>   dr;
+	Option<bool>   dr_out;
 	Option<float>  vn_level;
   	Option<int>    numICs;
   	Option<string> approach;
@@ -207,8 +209,14 @@ class MelodicOptions {
    migp_shuffle(string("--migp_order"), true,
 	   string("Randomise MIGP file order (default: TRUE)"),
 	   false, no_argument, false),
+   migp_factor(string("--migp_factor"), 2,
+	   string("Internal Factor of mem-threshold relative to number of Eigenmaps (default: 2)"),
+       false, requires_argument, false),
    dr(string("--dr"), false,
 	   string("Dual Regression (default: TRUE)"),
+	   false, no_argument, false),
+   dr_out(string("--dr_out"), false,
+	   string("Dual Regression output for MIGP/concat ICA"),
 	   false, no_argument, false),	
    vn_level(string("--vn_level"), float(2.3),
 	   string("variance nomalisation threshold level (Z> value is ignored)"), 
@@ -417,7 +425,9 @@ class MelodicOptions {
 		options.add(migp);
 		options.add(migpN);
 		options.add(migp_shuffle);
+		options.add(migp_factor);
 		options.add(dr);
+		options.add(dr_out);
 	    options.add(vn_level);
 	    options.add(numICs);
 	    options.add(approach);