MM report option added

meldata.cc

@@ -413,16 +413,21 @@ namespace Melodic{
 	  			// ICadjust = IC;
 				}	
 				else{
+					
 					Matrix resids = stdev(Data - mixMatrix * IC);
 					for(int ctr=1;ctr<=resids.Ncols();ctr++)
 						if(resids(1,ctr) < 0.05)
 							resids(1,ctr)=1;
-	  			stdNoisei = pow(stdev(Data - mixMatrix * IC)*
+	  //			stdNoisei = pow(stdev(Data - mixMatrix * IC)*
+		//				std::sqrt((float)(Data.Nrows()-1))/
+		//				std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
+	  			stdNoisei = pow(resids*
 						std::sqrt((float)(Data.Nrows()-1))/
 						std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
-	  
+	  	
 	  			ColumnVector diagvals;
 	  			diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5);
+	
 	  			save4D(SP(IC,diagvals*stdNoisei),opts.outputfname.value() + "_IC");
 				}
 

melica.cc

@@ -101,9 +101,11 @@ namespace Melodic{
 				}
 				//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())&&
+		//		  (opts.approach.value()!=string("tica"))){ break;}
+				if((abs(minAbsSin) < opts.epsilon.value())){ break;}
 	
 				itt_ctr++;
       } while(itt_ctr < opts.maxNumItt.value());
@@ -119,7 +121,7 @@ namespace Melodic{
 				temp = melodat.reduce_dimred(temp);
 				redUMM = melodat.get_white() * temp;
 
-				message(" change : " << abs(minAbsSin2) << endl);	
+				message(" change : " << minAbsSin2 << endl);	
 				if(abs(minAbsSin2) < opts.epsilonS.value() && abs(minAbsSin) < opts.epsilon.value()){ break;}
 			}
     } while(
@@ -467,6 +469,7 @@ namespace Melodic{
       melodat.set_mix(tmp);
 
       melodat.set_IC(temp);
+
       melodat.set_ICstats(scales);
       melodat.sort();
 

melodic.cc

@@ -114,7 +114,7 @@ int main(int argc, char *argv[]){
 	      				opts.pca_dim.set_T(opts.pca_dim.value()-retry*opts.retrystep);
 	    				}
 	    			else{
-	      			if(opts.pca_dim.value()-retry*opts.retrystep <  melodat.data_dim()){
+	      			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 
@@ -234,7 +234,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
     message(" done" << endl);
   }
 
-  message("Reading mixing matrix " << opts.filtermix.value());
+  message("Reading mixing matrix " << opts.filtermix.value() << " ... ");
   mixMatrix = read_ascii_matrix(opts.filtermix.value());
   if (mixMatrix.Storage()<=0) {
     cerr <<" Please specify the mixing matrix correctly" << endl;
@@ -257,13 +257,12 @@ void mmonly(Log& logger, MelodicOptions& opts,
     }
   }
 
-
   melodat.tempInfo = ICvolInfo;
   melodat.set_mask(Mask);
   melodat.set_mean(Mean);
   melodat.set_IC(ICs);
   melodat.set_mix(mixMatrix);
-  fmixMatrix = calc_FFT(melodat.expand_mix(), opts.logPower.value());
+  fmixMatrix = calc_FFT(mixMatrix, opts.logPower.value());
   melodat.set_fmix(fmixMatrix);
   fmixMatrix = pinv(mixMatrix);
   melodat.set_unmix(fmixMatrix);
@@ -272,7 +271,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
   
   Matrix mmres;
   Matrix pmaps;//(ICs);
-  if(opts.perf_mm.value())
+ if(opts.perf_mm.value())
     mmres = mmall(logger,opts,melodat,report,pmaps);
 	}