convergence

meldata.cc

@@ -63,7 +63,7 @@ namespace Melodic{
     }
 
     meanC = mean(tmpData,2);
-		tmpData = remmean(tmpData,2);
+	//	tmpData = remmean(tmpData,2);
 
     //convert to power spectra
     if(opts.pspec.value()){
@@ -86,7 +86,7 @@ namespace Melodic{
     //variance - normalisation
     if(opts.varnorm.value()){
       message("  Normalising by voxel-wise variance ..."); 
-      stdDev = varnorm(tmpData,tmpData.Nrows(),2.3);
+      stdDev = varnorm(tmpData,std::min(30,tmpData.Nrows()-1),2.3);
       stdDevi = pow(stdDev,-1); 
       message(" done" << endl);
     }
@@ -131,9 +131,10 @@ namespace Melodic{
 
   void MelodicData::set_TSmode()
   {
-    message("Calculating T- and S-modes " << endl);
     Matrix tmp, tmpT, tmpS, tmpT2, tmpS2;
     tmp = expand_dimred(mixMatrix);
+		saveascii(tmp,"exp");
+		saveascii(mixMatrix,"mix");
     tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
     tmpS = zeros(numfiles, tmp.Ncols());
     krfact(tmp,tmpT,tmpS);
@@ -175,11 +176,15 @@ namespace Melodic{
 //		bool tmpvarnorm = opts.varnorm.value();
 //		opts.varnorm.set_T(false);
     alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
-    for(int ctr = 1; ctr < numfiles; ctr++){
-      tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
+ 		if(opts.debug.value())
+			save4D(alldat,string("preproc_dat") + num2str(1));   
+		for(int ctr = 1; ctr < numfiles; ctr++){
+    	tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
+			if(opts.debug.value())
+				save4D(tmpData /numfiles,string("preproc_dat") + num2str(ctr+1));
       alldat += tmpData / numfiles;
     }
-
+	
 	//	opts.varnorm.set_T(tmpvarnorm);
 		//variance - normalisation
 	 // if(opts.varnorm.value()){
@@ -199,25 +204,37 @@ namespace Melodic{
 		if(numfiles>1)
     	message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);
 
+		if(opts.debug.value())
+			save4D(alldat,"alldat");
     //estimate model order
     ColumnVector tmpPPCA;
     RowVector AdjEV, PercEV;
     Matrix Corr, tmpE;
     int order;
 
-    order = ppca_dim(alldat, RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());	  
+    order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());	  
 
-    if(opts.pca_dim.value() == 0){
+	  if(opts.pca_dim.value() == 0){
       opts.pca_dim.set_T(order);
 			PPCA=tmpPPCA;
   	}
     order = opts.pca_dim.value();
     calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
-   
+
+    if(opts.debug.value()){
+			outMsize("pcaE",pcaE); saveascii(pcaE,"pcaE");
+			outMsize("pcaD",pcaD); saveascii(pcaD,"pcaD");
+			outMsize("AdjEV",AdjEV); saveascii(AdjEV,"AdjEV");
+			outMsize("PercEV",PercEV); saveascii(PercEV,"PercEV");
+			outMsize("tmpPPCA",tmpPPCA); saveascii(tmpPPCA,"tmpPPCA");
+			outMsize("whiteMatrix",whiteMatrix); saveascii(whiteMatrix,"whiteMatrix");
+			outMsize("dewhiteMatrix",dewhiteMatrix); saveascii(dewhiteMatrix,"dewhiteMatrix");
+			cerr << "Order: " << order << endl;
+		}
+
 		EV = AdjEV;
 		EVP = PercEV;
 
-
     if(numfiles == 1){
       Data = alldat;
       Matrix tmp = Identity(Data.Nrows());
@@ -230,15 +247,21 @@ namespace Melodic{
 //				if(opts.varnorm.value())
 //					varnorm(tmpData,stdDev);
 					
-				//  whiten (separate / joint) 
+				//  whiten (separate / joint)
+				Matrix newWM,newDWM; 
 				if(!opts.joined_whiten.value()){	  
 					message("  Individual whitening in a " << order << " dimensional subspace " << endl);
       	  std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
-	  			calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
+	  			calc_white(pcaE, pcaD, order, newWM, newDWM);
+				}else{
+					std_pca(whiteMatrix*tmpData, RXweight, Corr, pcaE, pcaD);
+					calc_white(pcaE, pcaD, order, newWM, newDWM);		
+					newDWM=(dewhiteMatrix*newDWM);
+					newWM=(newWM*whiteMatrix);
 				}
-				tmpData = whiteMatrix * tmpData;
-				DWM.push_back(dewhiteMatrix);
-				WM.push_back(whiteMatrix);
+				DWM.push_back(newDWM);
+				WM.push_back(newWM);
+				tmpData = newWM * tmpData;
 				
 				//concatenate Data
 				if(Data.Storage() == 0)
@@ -251,7 +274,8 @@ namespace Melodic{
     message("  Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl);
     
     meanC=mean(Data,2);
-        
+		if(opts.debug.value())
+			save4D(Data,"concat_data");    
     //save the mean & mask
     save_volume(Mask,logger.appendDir("mask"));
     save_volume(Mean,logger.appendDir("mean"));
@@ -295,17 +319,17 @@ namespace Melodic{
 
 		//read in post-proc design matrices etc
 		if(opts.fn_Tdesign.value().length()>0)
-			Tdes = read_vest(opts.fn_Tdesign.value());
+			Tdes = read_ascii_matrix(opts.fn_Tdesign.value());
 		if(opts.fn_Sdesign.value().length()>0)
-			Sdes = read_vest(opts.fn_Sdesign.value());
+			Sdes = read_ascii_matrix(opts.fn_Sdesign.value());
 		if(opts.fn_Tcon.value().length()>0)
-			Tcon = read_vest(opts.fn_Tcon.value());
+			Tcon = read_ascii_matrix(opts.fn_Tcon.value());
 		if(opts.fn_Scon.value().length()>0)
-			Scon = read_vest(opts.fn_Scon.value());
+			Scon = read_ascii_matrix(opts.fn_Scon.value());
 		if(opts.fn_TconF.value().length()>0)
-			TconF = read_vest(opts.fn_TconF.value());
+			TconF = read_ascii_matrix(opts.fn_TconF.value());
 		if(opts.fn_SconF.value().length()>0)
-			SconF = read_vest(opts.fn_SconF.value());
+			SconF = read_ascii_matrix(opts.fn_SconF.value());
 		
 		Tdes = remmean(Tdes,1);
 		Sdes = remmean(Sdes,1);	

melhlprfns.cc

@@ -132,7 +132,7 @@ namespace Melodic{
     Matrix tmpE, white, dewhite;
     RowVector tmpD, tmpD2;
 
-    std_pca(in, Corr, tmpE, tmpD);
+    std_pca(remmean(in,2), Corr, tmpE, tmpD);
     calc_white(tmpE,tmpD, dim, white, dewhite);
     
     Matrix ws = white * in;

melica.cc

@@ -32,7 +32,7 @@ namespace Melodic{
     // see www.cis.hut.fi/projects/ica/fastica/
     
     //initialize matrices
-    Matrix redUMM_old;
+    Matrix redUMM_old, redUMM_veryold;
     Matrix tmpU;    
     
     //srand((unsigned int)timer(NULL));
@@ -50,12 +50,13 @@ namespace Melodic{
     
     symm_orth(redUMM);
 
-    int itt_ctr,itt_ctr2=1,cum_itt=0; 
+    int itt_ctr,itt_ctr2=1,cum_itt=0,newmaxitts = opts.maxNumItt.value(); 
     double minAbsSin = 1.0;
     if(opts.approach.value() == string("tica"))
       opts.maxNumItt.set_T(opts.rank1interval.value());
     do{// TICA loop
       itt_ctr = 1;
+			redUMM_veryold = redUMM;
       do{ // da loop!!!
 				redUMM_old = redUMM;      
 				//calculate IC estimates
@@ -92,6 +93,9 @@ namespace Melodic{
 				// orthogonalize the unmixing-matrix 
 				symm_orth(redUMM);
 
+				if(opts.approach.value() == string("tica")){
+					message(" TICA:");
+				}
 				//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);
@@ -113,9 +117,12 @@ namespace Melodic{
 	  		outMsize(" 2nd unmmixing matrix ", temp);
 	  		redUMM = melodat.get_white() * temp;
 	  		outMsize(" redUMM ", redUMM);
+				minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_veryold)).Minimum();
+				message("  Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl);	
+				if(abs(minAbsSin) < opts.epsilon.value()){ break;}
 			}
     } while((itt_ctr2 < 3 || itt_ctr >= opts.maxNumItt.value()) && 
-			(opts.approach.value() ==  string("tica")) && cum_itt < 100*opts.maxNumItt.value());
+			(opts.approach.value() ==  string("tica")) && cum_itt < newmaxitts);
 
     if(itt_ctr>=opts.maxNumItt.value()){
       cerr << "  No convergence after " << cum_itt  <<" steps "<<endl;
@@ -459,6 +466,7 @@ namespace Melodic{
       melodat.set_IC(temp);
       melodat.set_ICstats(scales);
       melodat.sort();
+	    message("Calculating T- and S-modes " << endl);
       melodat.set_TSmode();
 		
     }

meloptions.h

@@ -251,7 +251,7 @@ class MelodicOptions {
 		 false, requires_argument),
    fn_TconF(string("--Tconf"), string(""),
 		 string("        F-contrast matrix across time-domain"),
-		 false, requires_argument),
+		 false, requires_argument, false),
    fn_Sdesign(string("--Sdes"), string(""),
 		 string("        design matrix across subject-domain"),
 		 false, requires_argument),

melreport.cc

@@ -692,7 +692,7 @@ namespace Melodic{
 				newplot.add_label("dim. estimate");
       }
 
-      //CB's latest BUG!!!			newplot.grid_swapdefault();
+      newplot.grid_swapdefault();
       newplot.timeseries(what,
 			 	report.appendDir("EVplot.png"),
 			 	string("Eigenspectrum Analysis"),