*** empty log message ***

test.cc

@@ -30,242 +30,57 @@ using namespace std;
   Option<string> fnin(string("-i,--in"), string(""),
 		string("input file name (matrix 3D or 4D image)"),
 		true, requires_argument);
-  Option<string> fnout(string("-o,--out"), string(""),
-		string(""),
-		true, requires_argument);
-  Option<string> fndesign(string("-d,--design"), string(""),
-		string("file name of the GLM design matrix (time courses or spatial maps)"),
-		true, requires_argument);
-  Option<string> fnmask(string("-m,--mask"), string(""),
-		string("mask image"),
-		false, requires_argument);
-  Option<string> fncontrasts(string("-c,--contrasts"), string(""),
-		string("matrix of t-statistics contrasts"),
-		false, requires_argument);
-  Option<string> fnftest(string("-f,--ftests"), string(""),
-		string("matrix of F-tests on contrasts"),
-		false, requires_argument);
-	Option<int> dofset(string("--dof"),0,
-		string("set degrees-of-freedom explicitly"),
-		false, requires_argument);
-	Option<string> filter(string("--filter"),string(""),
-		string("filter out part of the regression model"),
-		false, requires_argument);
-	Option<bool> perfvn(string("--vn"),FALSE,
-		string("perfrom variance-normalisation on data"),
-		false, requires_argument);
 	Option<int> help(string("-h,--help"), 0,
 		string("display this help text"),
 		false,no_argument);
-	// Output options	
-	Option<string> outcope(string("--out_cope"),string(""),
-		string("output COPEs"),
-		false, requires_argument);
-	Option<string> outz(string("--out_z"),string(""),
-		string("output Z-stats"),
-		false, requires_argument);
-	Option<string> outt(string("--out_t"),string(""),
-		string("output t-stats"),
-		false, requires_argument);
-	Option<string> outp(string("--out_p"),string(""),
-		string("output p-values of Z-stats"),
-		false, requires_argument);
-	Option<string> outf(string("--out_f"),string(""),
-		string("output F-value of full model fit"),
-		false, requires_argument);
-	Option<string> outpf(string("--out_pf"),string(""),
-		string("output p-value for full model fit"),
-		false, requires_argument);
-	Option<string> outfilt(string("--out_filt"),string(""),
-		string("output filtered data"),
-		false, requires_argument);
-	Option<string> outres(string("--out_res"),string(""),
-		string("output residuals"),
-		false, requires_argument);
-	Option<string> outvarcb(string("--out_varcb"),string(""),
-		string("output variance of COPEs"),
-		false, requires_argument);
-	Option<string> outsigsq(string("--out_sigsq"),string(""),
-		string("output residual noise variance sigma-square"),
-		false, requires_argument);
-	Option<string> outdata(string("--out_data"),string(""),
-		string("output data"),
-		false, requires_argument);
-	Option<string> outvnscales(string("--out_vnscales"),string(""),
-		string("output scaling factors for variance normalisation"),
-		false, requires_argument);
 		/*
 }
 */
 //Globals {
-	Melodic::basicGLM glm;
-	bool isimage = FALSE;
-	Matrix data;
-	Matrix design;
-	Matrix contrasts;
-	Matrix fcontrasts;
-	Matrix meanR;
-	RowVector vnscales;
-	volume<float> mask;
+	
 	volumeinfo volinf;  /*
 }
 */
 ////////////////////////////////////////////////////////////////////////////
 
 // Local functions
-void save4D(Matrix what, string fname){
-		if(what.Ncols()==data.Ncols()||what.Nrows()==data.Nrows()){
-			volume4D<float> tempVol;
-			if(what.Nrows()>what.Ncols())
-				tempVol.setmatrix(what.t(),mask);
-			else
-				tempVol.setmatrix(what,mask);
-			save_volume4D(tempVol,fname,volinf);
-		}
-}
-
-void saveit(Matrix what, string fname){
-	if(isimage)
-		save4D(what,fname);
-	else
-		write_ascii_matrix(what,fname);
-}
-
-int dofilter(){
-	if(!isimage){
-		cerr << "ERROR: need to specify 4D input to use filtering" << endl;
-		return 1;
-	}
-  Matrix unmixMatrix = pinv(design);
-  Matrix maps = unmixMatrix * data;
 
-  Matrix noisedes;
-  Matrix noisemaps;
-
-  int ctr=0;    
-  char *p;
-  char t[1024];
-  const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
-  
-  strcpy(t, filter.value().c_str());
-  p=strtok(t,discard);
-  ctr = atoi(p);
-  if(ctr>0 && ctr<=design.Ncols()){
-    noisedes = design.Column(ctr);
-    noisemaps  = maps.Row(ctr).t();    
-  }
-
-  do{
-    p=strtok(NULL,discard);
-    if(p){
-			ctr = atoi(p);
-      if(ctr>0 && ctr<=design.Ncols()){
-  			noisedes |= design.Column(ctr);
-  			noisemaps  |= maps.Row(ctr).t();
-			}
-    }
-  }while(p);
-  Matrix newData;
-  newData = data - noisedes * noisemaps.t();
-  newData = newData + ones(newData.Nrows(),1)*meanR;
-  
-	save4D(newData,outfilt.value());
-  return 0;	
-}
-
-int setup(){
-	if(fsl_imageexists(fnin.value())){//read data
-		//input is 3D/4D vol
-		isimage = TRUE;
-		volume4D<float> tmpdata;
-		read_volume4D(tmpdata,fnin.value(),volinf);
+int do_work(int argc, char* argv[]) {
+	Matrix design;
+	design=read_ascii_matrix(fnin.value());
+		Matrix joined;
+		Matrix weights;
 		
-		// create mask
-		if(fnmask.value()>""){
-			read_volume(mask,fnmask.value());
-			if(!samesize(tmpdata[0],mask)){
-				cerr << "ERROR: Mask image does not match input image" << endl;
-				return 1;
-			};
-		}else{
-			mask = tmpdata[0]*0.0+1.0;	
-		}
 		
-		data = tmpdata.matrix(mask);
-	}
-	else
-		data = read_ascii_matrix(fnin.value());	
-
-	if(fsl_imageexists(fndesign.value())){//read design
-		volume4D<float> tmpdata;
-		read_volume4D(tmpdata,fndesign.value());
-		if(!samesize(tmpdata[0],mask)){
-			cerr << "ERROR: GLM design does not match input image in size" << endl;
-			return 1;
-		}
-		design = tmpdata.matrix(mask).t();
-		data = data.t();
-		isimage = FALSE;
-	}else{
-		design = read_ascii_matrix(fndesign.value());
-	}
-
-	meanR=mean(data,1);
-	data = remmean(data,1);
-	design = remmean(design,1);
-	if(perfvn.value())
-		vnscales = Melodic::varnorm(data);
-	if(fncontrasts.value()>""){//read contrast		
-		contrasts = read_ascii_matrix(fncontrasts.value());
-		if(!(contrasts.Ncols()==design.Ncols())){
-			cerr << "ERROR: contrast matrix GLM design does not match GLM design" << endl;
-			return 1;
+	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;
 		}
-	}else{
-		contrasts = Identity(design.Ncols());
-		contrasts &= -1.0 * contrasts;
+			
+		
 	}
-	return 0;	
-}
+	
+	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;
 
-void write_res(){	
-	if(fnout.value()>"")
-		saveit(glm.get_beta(),fnout.value());
-	if(outcope.value()>"")
-		saveit(glm.get_cbeta(),outcope.value());
-	if(outz.value()>"")
-		saveit(glm.get_z(),outz.value());
-	if(outt.value()>"")
-		saveit(glm.get_t(),outt.value());
-	if(outp.value()>"")
-		saveit(glm.get_p(),outp.value());
-	if(outf.value()>"")
-		saveit(glm.get_f_fmf(),outf.value());
-	if(outpf.value()>"")
-		saveit(glm.get_pf_fmf(),outpf.value());
-	if(outres.value()>"")
-		saveit(glm.get_residu(),outres.value());
-	if(outvarcb.value()>"")
-		saveit(glm.get_varcb(),outvarcb.value());
-	if(outsigsq.value()>"")
-		saveit(glm.get_sigsq(),outsigsq.value());
-	if(outdata.value()>"")
-		saveit(data,outdata.value());
-	if(outvnscales.value()>"")
-		saveit(vnscales,outvnscales.value());
-}
+	cerr << weights << endl <<endl << Rows << endl;
 
-int do_work(int argc, char* argv[]) {
-  if(setup())
-		exit(1);
-	
-	if(filter.value()>"")
-		if(dofilter())
-		exit(1);	
-	else{
-		glm.olsfit(data,design,contrasts,dofset.value());
-		write_res();
-  }
+	cerr << SP(weights,pow(Rows,-1)) << endl;
 	return 0;
 }
 
@@ -278,27 +93,9 @@ int do_work(int argc, char* argv[]) {
 	    // must include all wanted options here (the order determines how
 	    //  the help message is printed)
 			options.add(fnin);
-			options.add(fnout);
-			options.add(fndesign);
-			options.add(fnmask);
-			options.add(fncontrasts);
-			options.add(fnftest);
-			options.add(dofset);
-			options.add(filter);
-			options.add(perfvn);
+		
 			options.add(help);
-			options.add(outcope);
-			options.add(outz);
-			options.add(outt);
-			options.add(outp);
-			options.add(outf);
-			options.add(outpf);
-			options.add(outfilt);
-			options.add(outres);
-			options.add(outvarcb);
-			options.add(outsigsq);
-			options.add(outdata);
-			options.add(outvnscales);
+		
 	    options.parse_command_line(argc, argv);
 
 	    // line below stops the program if the help was requested or