From 7ae0edbd6c036b32fa8dc909e7a67f57346d0a65 Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Mon, 8 Dec 2008 22:15:08 +0000
Subject: [PATCH] SBCA modified
---
Makefile | 11 ++++--
fsl_mvlm.cc | 99 +++++++++++++++++++++++++++++++++++++++++------------
melica.cc | 3 ++
test.cc | 61 +++++++++------------------------
4 files changed, 105 insertions(+), 69 deletions(-)
diff --git a/Makefile b/Makefile
index 552f8ad..b8ea47c 100644
--- a/Makefile
+++ b/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}
diff --git a/fsl_mvlm.cc b/fsl_mvlm.cc
index bbe57c7..cafe3e8 100644
--- a/fsl_mvlm.cc
+++ b/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);
diff --git a/melica.cc b/melica.cc
index c80b17d..cf00974 100644
--- a/melica.cc
+++ b/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);
diff --git a/test.cc b/test.cc
index b51a5b4..549caf2 100644
--- a/test.cc
+++ b/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;
--
GitLab