SBCA modified

Makefile

@@ -12,12 +12,14 @@ USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB} -L${LIB_GD} -L${LIB_GDC} -L${LIB_PNG}
 
 LIBS = -lutils -lnewimage -lmiscplot -lmiscpic -lmiscmaths -lfslio -lniftiio -lznz -lnewmat -lprob -lm  -lgdc -lgd -lpng -lz
 
-TEST_OBJS = test.o test2.o
+TEST_OBJS = test.o 
 
 GGMIX_OBJS = ggmix.o
 
 FSL_GLM_OBJS = melhlprfns.o fsl_glm.o
 
+FSL_SBCA_OBJS = melhlprfns.o fsl_sbca.o
+
 FSL_MVLM_OBJS = melhlprfns.o fsl_mvlm.o
 
 FSL_REGFILT_OBJS = melhlprfns.o fsl_regfilt.o
@@ -26,13 +28,13 @@ MELODIC_OBJS =  meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o
 
 TESTXFILES = test
 
-XFILES = fsl_glm fsl_mvlm fsl_regfilt melodic
+XFILES = fsl_glm fsl_sbca fsl_mvlm fsl_regfilt melodic
 
 RUNTCLS = Melodic
 
 SCRIPTS = melodicreport
 
-all: ggmix fsl_regfilt fsl_glm fsl_mvlm melodic
+all: ggmix fsl_regfilt fsl_glm fsl_mvlm fsl_sbca melodic
 
 ggmix: ${GGMIX_OBJS}
 	${AR} -r libggmix.a ${GGMIX_OBJS} 
@@ -46,6 +48,9 @@ test: ${TEST_OBJS}
 fsl_glm: ${FSL_GLM_OBJS}
 	$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_GLM_OBJS} ${LIBS}
 
+fsl_sbca: ${FSL_SBCA_OBJS}
+		$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_SBCA_OBJS} ${LIBS}
+
 fsl_mvlm: ${FSL_MVLM_OBJS}
 		$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_MVLM_OBJS} ${LIBS}
 

fsl_mvlm.cc

@@ -64,7 +64,7 @@ using namespace std;
 		string("basename for output files "),
 		true, requires_argument);
 	Option<string> approach(string("-a,--alg"), string("PCA"),
-		string("algorithm for decomposition: PCA (or SVD; default), PLS, orthoPLS, CVA, SVD-CVA, MLM"),
+		string("algorithm for decomposition: PCA (or SVD; default), PLS, orthoPLS, CVA, SVD-CVA, MLM, NMF"),
 		false, requires_argument);	
     Option<string> fndesign(string("-d,--design"), string(""),
 		string("file name of the GLM design matrix (time courses or spatial maps)"),
@@ -81,6 +81,12 @@ using namespace std;
 	Option<bool> perf_demean(string("--demean"),FALSE,
 		string("switch on de-meaning of design and data"),
 		false, no_argument);
+	Option<int> nmfdim(string("--nmf_dim"), 0,
+		string(" Number of underlying factors for NMF"),
+		false,requires_argument);
+	Option<int> nmfitt(string("--nmfitt"), 100,
+		string("number of NMF itterations (default 100)"),
+		false,requires_argument);
 	Option<int> help(string("-h,--help"), 0,
 		string("display this help text"),
 		false,no_argument);
@@ -202,10 +208,12 @@ int setup(){
 		design =  SP(design,ones(design.Nrows(),1)*pow(stdev(design,1),-1));
 	
 	SVD( design, svd_X_D, svd_X_U, svd_X_V );
-	if(data.Nrows()>=data.Ncols())
-		SVD ( data, svd_Y_D, svd_Y_U, svd_Y_V );
-	else{
-		SVD ( data.t(), svd_Y_D, svd_Y_V, svd_Y_U );
+	if(approach.value()!=string("NMF")){
+		if(data.Nrows()>=data.Ncols())
+			SVD ( data, svd_Y_D, svd_Y_U, svd_Y_V );
+		else{
+			SVD ( data.t(), svd_Y_D, svd_Y_V, svd_Y_U );
+		}
 	}		
 
 	if(fnout.value().length() == 0){
@@ -303,30 +311,75 @@ int do_work(int argc, char* argv[]) {
 	}
 
     if( approach.value()!=string("MLM") && approach.value()!=string("CVA") && approach.value()!=string("PLS") &&
-		approach.value()!=string("SVD-CVA") && approach.value()!=string("orthoPLS"))
+		approach.value()!=string("SVD-CVA") && approach.value()!=string("orthoPLS") && approach.value()!=string("NMF"))
 		message(" Using method : PCA" << endl;);
 	
 	//perform an SVD on data	
 	outMsize(" New Data ", data);
+	
+	if(approach.value()!=string("NMF")){	
+		if(data.Nrows()>=data.Ncols())
+			SVD ( data, svd_Y_D, svd_Y_U, svd_Y_V );
+		else{
+			SVD ( data.t(), svd_Y_D, svd_Y_V, svd_Y_U );
+			svd_Y_U = svd_Y_U.t();
+			svd_Y_V = svd_Y_V.t();		
+		}	
+
+		dbgmsg(" Finished SVD : ");
+		outMsize("svd_Y_U",svd_Y_U);
+		outMsize("svd_Y_V",svd_Y_V);
+		outMsize("svd_Y_D",svd_Y_D);
+	
+		svd_Y_V = sqrtm(svd_Y_D) * svd_Y_V;	
+		svd_Y_U = svd_Y_U * sqrtm(svd_Y_D);	
+
+		if(approach.value()==string("SVD-CVA"))
+			svd_Y_V = svd_Y_V *tmpdata;
+	}
+	else{ //NMF
+		float err, err_old;
+		Matrix Ratio, Diff;
+
+		if(nmfdim.value()==0)
+			nmfdim.set_T(data.Nrows());
 		
-	if(data.Nrows()>=data.Ncols())
-		SVD ( data, svd_Y_D, svd_Y_U, svd_Y_V );
-	else{
-		SVD ( data.t(), svd_Y_D, svd_Y_V, svd_Y_U );
-		svd_Y_U = svd_Y_U.t();
-		svd_Y_V = svd_Y_V.t();		
-	}	
+		message("Using "<< nmfdim.value() << " dimensions" << endl;); 	
+		svd_Y_U =  unifrnd(data.Nrows(), nmfdim.value());
+		svd_Y_V =  unifrnd(nmfdim.value(), data.Ncols());
+		// re-scale columns of svd_Y_U to unit amplitude
+		Ratio = pow(stdev(svd_Y_U),-1);
+		svd_Y_U = SP(svd_Y_U,ones(svd_Y_U.Nrows(),1)*Ratio);
+				
+		Diff = data - svd_Y_U * svd_Y_V;
+		err = Diff.SumAbsoluteValue()/(data.Ncols()*data.Nrows());
+			
+		for(int k=1; k< nmfitt.value(); k++)
+		{
+		//	Ratio = SP(data,pow(svd_Y_U * svd_Y_V,-1));
+		//	svd_Y_U  =  SP(svd_Y_U, Ratio * svd_Y_V.t());
+		//	svd_Y_U  =  SP(svd_Y_U, pow( ones(svd_Y_U.Nrows(),1) * sum(svd_Y_U,1),-1));
+		//	svd_Y_V  =  SP(svd_Y_V, svd_Y_U.t()* Ratio);
+		//	
+		// Lee & Seung multiplicatice updates
+		
+			Ratio = SP(svd_Y_U.t() * data, pow(svd_Y_U.t() * svd_Y_U * svd_Y_V ,-1));
+			svd_Y_V  = SP(svd_Y_V,Ratio);
 
-	dbgmsg(" Finished SVD : ");
-	outMsize("svd_Y_U",svd_Y_U);
-	outMsize("svd_Y_V",svd_Y_V);
-	outMsize("svd_Y_D",svd_Y_D);
-	
-	svd_Y_V = sqrtm(svd_Y_D) * svd_Y_V;	
-	svd_Y_U = svd_Y_U * sqrtm(svd_Y_D);	
+			Ratio = SP(data * svd_Y_V.t(),pow(svd_Y_U * (svd_Y_V * svd_Y_V.t()),-1));
+			svd_Y_U  = SP(svd_Y_U,Ratio);
 
-	if(approach.value()==string("SVD-CVA"))
-		svd_Y_V = svd_Y_V *tmpdata;
+			// re-scale columns of svd_Y_U to unit amplitude
+			Ratio = pow(stdev(svd_Y_U),-1);
+			svd_Y_U = SP(svd_Y_U,ones(svd_Y_U.Nrows(),1)*Ratio);
+			
+			Diff = data - svd_Y_U * svd_Y_V;
+			err_old = err;
+			err = Diff.SumSquare()/(data.Ncols()*data.Nrows());			
+			message(" Error " << err << endl;);
+		}	
+	}
+	
 	write_res();
 	
 	dbgmsg(" Leaving <do_work>");
@@ -350,6 +403,8 @@ int main(int argc,char *argv[]){
 			options.add(normdes);
 			options.add(perfvn);
 			options.add(perf_demean);
+			options.add(nmfdim);
+			options.add(nmfitt);
 			options.add(help);
 			options.add(verbose);
 			options.add(debug);

melica.cc

@@ -39,6 +39,9 @@ namespace Melodic{
     redUMM = melodat.get_white()*
        unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
     
+	if(opts.debug.value())
+		cerr << "redUMM init submatrix : " << endl << redUMM.SubMatrix(1,2,1,2) << endl;
+		
     if(opts.guessfname.value().size()>1){
       message("  Use columns in " << opts.guessfname.value() 
 	      << " as initial values for the mixing matrix " <<endl);

test.cc

@@ -9,6 +9,7 @@
 #include "miscmaths/miscprob.h"
 #include "utils/options.h"
 #include <vector>
+#include <time.h>
 #include "newimage/newimageall.h"
 #include "melhlprfns.h"
 
@@ -20,11 +21,11 @@ using namespace std;
 // The two strings below specify the title and example usage that is
 // printed out as the help or usage message
 
-  string title=string("fsl_glm (Version 1.0)")+
-		string("\nCopyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
-		string(" \n Simple GLM usign ordinary least-squares regression on\n")+
-		string(" time courses and/or 3D/4D volumes\n\n");
-  string examples="fsl_glm <input> -d <design> [options]";
+  string title=string("fsl_BLAH")+
+		string("\nCopyright(c) 2008, University of Oxford (Christian F. Beckmann)\n")+
+		string(" \n \n")+
+		string(" \n");
+  string examples="fsl_BLAH [options]";
 
 //Command line Options {
   Option<string> fnin(string("-i,--in"), string(""),
@@ -46,43 +47,7 @@ using namespace std;
 // Local functions
 
 int do_work(int argc, char* argv[]) {
-	Matrix design;
-	design=read_ascii_matrix(fnin.value());
-		Matrix joined;
-		Matrix weights;
-		
-		
-	cerr << " Design : "  << design.Nrows() << " x " << design.Ncols() << endl <<endl;
-	for(int i=1; i <10; i++){
-		weights = normrnd(10,design.Ncols());
-		
-	
-		joined=SP(design,ones(design.Nrows(),1)*weights.Row(1));
-		
-		for(int j=2;j<=weights.Nrows();j++){
-			Matrix tmp;
-			tmp=SP(design,ones(design.Nrows(),1)*weights.Row(j));
-			joined &= tmp;
-		}
-			
-		
-	}
-	
-	cerr << " joined : " << joined.Nrows() << " x " << joined.Ncols() << endl << endl;
-	
-	Matrix Cols,Rows;
-	
-	Melodic::krfact(joined,design.Nrows(),Cols,Rows);
-	
-	cerr << " Cols : " << Cols.Nrows() << " x " << Cols.Ncols() << endl << endl;
-	
-	cerr << " Rows : " << Rows.Nrows() << " x " << Rows.Ncols() << endl << endl;
-
-	cerr << weights << endl <<endl << Rows << endl;
 
-	cerr << remmean(SP(weights,pow(Rows,-1))) << endl;
-	
-	
 	return 0;
 }
 
@@ -94,8 +59,16 @@ int do_work(int argc, char* argv[]) {
 	  try{
 	    // must include all wanted options here (the order determines how
 	    //  the help message is printed)
-			options.add(fnin);
-		
+	
+		  double tmptime = time(NULL);
+		  srand((unsigned int) tmptime);
+
+		cerr << (unsigned int) tmptime << endl << endl;
+
+		cerr << unifrnd(2,2) << endl;
+		exit(1);
+	/*
+			options.add(fnin);		
 			options.add(help);
 		
 	    options.parse_command_line(argc, argv);
@@ -108,7 +81,7 @@ int do_work(int argc, char* argv[]) {
 	    }else{
 	  		// Call the local functions
 	  		return do_work(argc,argv);
-			}
+			}*/
 		}catch(X_OptionError& e) {
 			options.usage();
 	  	cerr << endl << e.what() << endl;