GLM reporting / unequal length fix

meldata.cc

@@ -144,7 +144,7 @@ namespace Melodic{
 	
     tmp = expand_dimred(mixMatrix);
     tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
-    tmpS = zeros(numfiles, tmp.Ncols());
+    tmpS = ones(numfiles, tmp.Ncols());
 
 	outMsize("tmp",tmp);
 	outMsize("tmpT",tmpT);
@@ -153,49 +153,60 @@ namespace Melodic{
   	dbgmsg(string("   approach ") << opts.approach.value() << endl);	
 
 	if(opts.approach.value()!=string("concat")){
-		message("Calculating T- and S-modes " << endl);
-    	explained_var = krfact(tmp,tmpT,tmpS);
-    	if(opts.approach.value()==string("tica")){		
-    		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);
+      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 GLM OLS fit
-	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);
+      	add_Tmodes(tmpT2);
+      	add_Smodes(tmpS2);
+	  }
 	}
-	
-    }
 	else{
-		add_Tmodes(tmp);
+		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);
+//      }
 	
-	dbgmsg(string("END: set_TSmode"));
-	}
+	dbgmsg(string("END: set_TSmode"));	
   }
 
   void MelodicData::setup()
@@ -233,6 +244,11 @@ namespace Melodic{
 			if(tmpData.Ncols() == alldat.Ncols() && tmpData.Nrows() == alldat.Nrows())
       			alldat += tmpData;	
 			else{
+				    if(opts.approach.value()==string("tica")){
+						cerr << "ERROR:: data dimensions do not match, TICA not possible \n\n";
+						exit(2); 
+					}
+					
 					if(tmpData.Ncols() == alldat.Ncols()){
 						int mindim = min(alldat.Nrows(),tmpData.Nrows());
 						alldat = alldat.Rows(1,mindim);
@@ -511,6 +527,7 @@ namespace Melodic{
       }
      
     //Output T- & S-modes
+ 
     save_Tmodes();
     save_Smodes();
 

meldata.h

@@ -90,8 +90,11 @@ namespace Melodic{
       inline void save_Tmodes(){
 			if(Tmodes.size()>0){
 				Matrix tmp = Tmodes.at(0); 
-				for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++)
-	  			tmp |= Tmodes.at(ctr);
+				outMsize("tmp",tmp);
+				for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++){
+					outMsize("Tmodes ",Tmodes.at(ctr));
+     	  			tmp |= Tmodes.at(ctr);
+				}
 				saveascii(tmp,opts.outputfname.value() + "_Tmodes");
 			}
       }