From da03bc3faa6e41adc77801365d47e0a98ccdcb01 Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Mon, 10 Nov 2008 17:58:52 +0000
Subject: [PATCH] added fsl_mvlm
---
Makefile | 9 +-
fsl_glm.cc | 2 +-
fsl_mvlm.cc | 381 ++++++++++++++++++++++++++++++++++++++++++++++++++++
3 files changed, 389 insertions(+), 3 deletions(-)
create mode 100644 fsl_mvlm.cc
diff --git a/Makefile b/Makefile
index ecf8f16..552f8ad 100644
--- a/Makefile
+++ b/Makefile
@@ -18,19 +18,21 @@ GGMIX_OBJS = ggmix.o
FSL_GLM_OBJS = melhlprfns.o fsl_glm.o
+FSL_MVLM_OBJS = melhlprfns.o fsl_mvlm.o
+
FSL_REGFILT_OBJS = melhlprfns.o fsl_regfilt.o
MELODIC_OBJS = meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o melreport.o melodic.o
TESTXFILES = test
-XFILES = fsl_glm fsl_regfilt melodic
+XFILES = fsl_glm fsl_mvlm fsl_regfilt melodic
RUNTCLS = Melodic
SCRIPTS = melodicreport
-all: ggmix fsl_regfilt fsl_glm melodic
+all: ggmix fsl_regfilt fsl_glm fsl_mvlm melodic
ggmix: ${GGMIX_OBJS}
${AR} -r libggmix.a ${GGMIX_OBJS}
@@ -44,6 +46,9 @@ test: ${TEST_OBJS}
fsl_glm: ${FSL_GLM_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_GLM_OBJS} ${LIBS}
+fsl_mvlm: ${FSL_MVLM_OBJS}
+ $(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_MVLM_OBJS} ${LIBS}
+
fsl_regfilt: ${FSL_REGFILT_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_REGFILT_OBJS} ${LIBS}
diff --git a/fsl_glm.cc b/fsl_glm.cc
index 0d78bab..bc5ff31 100644
--- a/fsl_glm.cc
+++ b/fsl_glm.cc
@@ -58,7 +58,7 @@ using namespace std;
string(" perform MELODIC variance-normalisation on data"),
false, no_argument);
Option<bool> perf_demean(string("--demean"),FALSE,
- string(" switch on de-meaning of design and data"),
+ string("switch on de-meaning of design and data"),
false, no_argument);
Option<int> help(string("-h,--help"), 0,
string("display this help text"),
diff --git a/fsl_mvlm.cc b/fsl_mvlm.cc
new file mode 100644
index 0000000..619c712
--- /dev/null
+++ b/fsl_mvlm.cc
@@ -0,0 +1,381 @@
+/* fsl_glm -
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 2006-2008 University of Oxford */
+
+/* CCOPYRIGHT */
+
+//Header & includes
+#include "libvis/miscplot.h"
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/miscprob.h"
+#include "utils/options.h"
+#include <vector>
+#include "newimage/newimageall.h"
+#include "melhlprfns.h"
+
+#define message(msg) { \
+ if(verbose.value()) \
+ { \
+ cout << msg; \
+ } \
+}
+
+#define dbgmsg(msg) { \
+ if(debug.value()) {\
+ cerr << msg << endl; } \
+}
+
+#define outMsize(msg,Mat) { \
+ if(debug.value()) \
+ cerr << " " << msg << " " <<Mat.Nrows() << " x " << Mat.Ncols() << endl; \
+}
+
+#define outM(msg,Mat) { \
+ if(verbose.value()) \
+ cout << " " << msg << " " <<Mat.Nrows() << " x " << Mat.Ncols() << endl; \
+}
+
+
+using namespace MISCPLOT;
+using namespace MISCMATHS;
+using namespace Utilities;
+using namespace std;
+
+
+//Command-line output
+// The two strings below specify the title and example usage that is
+// printed out as the help or usage message
+
+ string title=string("fsl_mvlm (Version 1.0)")+
+ string("\nCopyright(c) 2008, University of Oxford (Christian F. Beckmann)\n")+
+ string(" \n Multivariate Linear Model regression on\n")+
+ string(" time courses and/or 3D/4D imges using SVD (PCA), PLS, normalised PLS, \n")+
+ string(" CVA, SVD-CVA or MLM\n\n");
+ string examples="fsl_mvlm -i <input> -o <output> [options]";
+
+
+//Command line Options
+ Option<string> fnin(string("-i,--in"), string(""),
+ string(" input file name (text matrix or 3D/4D image file)"),
+ true, requires_argument);
+ Option<string> fnout(string("-o,--out"), string(""),
+ 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, orthPLS, CVA, SVD-CVA, MLM"),
+ false, requires_argument);
+ Option<string> fndesign(string("-d,--design"), string(""),
+ string("file name of the GLM design matrix (time courses or spatial maps)"),
+ false, requires_argument);
+ Option<string> fnmask(string("-m,--mask"), string(""),
+ string("mask image file name if input is image"),
+ false, requires_argument);
+ Option<bool> normdes(string("--des_norm"),FALSE,
+ string("switch on normalisation of the design matrix columns to unit std. deviation"),
+ false, no_argument);
+ Option<bool> perfvn(string("--vn"),FALSE,
+ string(" perform MELODIC variance-normalisation on data"),
+ false, no_argument);
+ Option<bool> perf_demean(string("--demean"),FALSE,
+ string("switch on de-meaning of design and data"),
+ false, no_argument);
+ Option<int> help(string("-h,--help"), 0,
+ string("display this help text"),
+ false,no_argument);
+ Option<bool> verbose(string("-v,--verbose"),FALSE,
+ string("switch on verbose output"),
+ false, no_argument);
+ Option<bool> debug(string("--debug"),FALSE,
+ string("switch on debug output"),
+ false, no_argument, false);
+ // Output options
+ Option<string> outres(string("--out_res"),string(""),
+ string("output file name for residuals"),
+ false, requires_argument, false);
+ Option<string> outdata(string("--out_data"),string(""),
+ string("output file name for pre-processed data"),
+ false, requires_argument);
+ Option<string> outvnscales(string("--out_vnscales"),string(""),
+ string("output file name for scaling factors for variance normalisation"),
+ false, requires_argument);
+
+//Globals
+ Melodic::basicGLM glm;
+ int voxels = 0;
+ Matrix data, tmpdata;
+ Matrix design;
+ Matrix meanR;
+ Matrix svd_X_U, svd_X_V, svd_Y_U, svd_Y_V;
+ DiagonalMatrix svd_X_D, svd_Y_D;
+ 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);
+ }
+}
+
+bool isimage(Matrix what){
+ if((voxels > 0)&&(what.Ncols()==voxels || what.Nrows()==voxels))
+ return TRUE;
+ else
+ return FALSE;
+}
+
+void saveit(Matrix what, string fname){
+ if(isimage(what))
+ save4D(what,fname);
+ else if(fsl_imageexists(fndesign.value()))
+ write_ascii_matrix(what.t(),fname);
+ else
+ write_ascii_matrix(what,fname);
+}
+
+int setup(){
+
+ dbgmsg(" In <setup>");
+ message(" Reading data " << fnin.value() << " ... ");
+ if(fsl_imageexists(fnin.value())){//read data
+ //input is 3D/4D vol
+ volume4D<float> tmpdata;
+ read_volume4D(tmpdata,fnin.value(),volinf);
+
+ // 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);
+ voxels = data.Ncols();
+ if(perf_demean.value())
+ data = remmean(data,1);
+ if(perfvn.value())
+ vnscales = Melodic::varnorm(data);
+ }
+ else
+ data = read_ascii_matrix(fnin.value());
+
+ message("done" << endl;);
+ if(fndesign.value().length()>0){
+ message(" Reading design " << fndesign.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();
+ }else{
+ design = read_ascii_matrix(fndesign.value());
+ }
+ message("done" << endl;);
+ }else
+ design = ones(data.Nrows(),1);
+
+ meanR=mean(data,1);
+ if(perf_demean.value()){
+ data = remmean(data,1);
+ design = remmean(design,1);
+ }
+
+ if(normdes.value())
+ 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(fnout.value().length() == 0){
+ string basename = fnin.value();
+ basename = make_basename(basename);
+ fnout.set_T(basename+string("_mvlm_"));
+ }
+
+ outM("Data matrix : ", data);
+ outM("Design matrix : ", design);
+
+ dbgmsg(" initial SVD : ");
+ outMsize("svd_Y_U",svd_Y_U);
+ outMsize("svd_Y_V",svd_Y_V);
+ outMsize("svd_Y_D",svd_Y_D);
+
+ outMsize("svd_X_U",svd_X_U);
+ outMsize("svd_X_V",svd_X_V);
+ outMsize("svd_X_D",svd_X_D);
+
+ dbgmsg(" Leaving <setup>");
+ return 0;
+}
+
+void write_res(){
+ dbgmsg(" In <write_res>");
+
+ message(" Writing results ... ")
+ if(isimage(svd_Y_V)){
+ saveit(svd_Y_V,fnout.value()+string("maps"));
+ saveit(svd_Y_U,fnout.value()+string("tcs"));
+ }
+ else{
+ saveit(svd_Y_V.t(),fnout.value()+string("tcs"));
+ saveit(svd_Y_U,fnout.value()+string("maps"));
+ }
+ saveit(svd_Y_D.AsColumn(),fnout.value()+string("scales"));
+ if(outres.value()>"")
+
+ if(outdata.value()>"")
+ saveit(data,outdata.value());
+ if(outvnscales.value()>"")
+ saveit(vnscales,outvnscales.value());
+
+ message("done" << endl;);
+ dbgmsg(" Leaving <write_res>");
+
+}
+
+int do_work(int argc, char* argv[]) {
+ dbgmsg(" In <do_work>");
+
+ if(setup())
+ exit(1);
+
+ //modify data
+
+ //X = svd_X_U * svd_X_D * svd_X_V.t();
+ //Y = svd_Y_U * svd_Y_D * svd_Y_V.t();
+ //X'X = svd_X_V *pow(svd_X_D,2) * svd_X_V.t();
+ //(X'X)^(-1) = svd_X_V *pow(svd_X_D,-2) * svd_X_V.t()
+ //(X'X)^(-1/2) = svd_X_V *pow(svd_X_D,-1) * svd_X_V.t()
+
+ if(approach.value()==string("PLS")) {
+ message(" Using method : " << approach.value() << endl;);
+ data = design.t() * data;
+ }
+ if(approach.value()==string("orthoPLS")) {
+ message(" Using method : " << approach.value() << endl;);
+ data = (svd_X_V * svd_X_D.i() * svd_X_V.t()) * design.t() * data;
+ }
+ if(approach.value()==string("CVA")) {
+ message(" Using method : " << approach.value() << endl;);
+ data = design.t() * svd_Y_U * svd_Y_V.t();
+ data = (svd_X_V * svd_X_D.i() * svd_X_V.t() ) * data;
+ }
+ if(approach.value()==string("SVD-CVA")) {
+ message(" Using method : " << approach.value() << endl;);
+ tmpdata = data;
+ data = design.t() * svd_Y_U;
+ data = (svd_X_V * svd_X_D.i() * svd_X_V.t() ) * data;
+ }
+ if(approach.value()==string("MLM")) {
+ message(" Using method : " << approach.value() << endl;);
+
+ Matrix RE;
+ DiagonalMatrix RD;
+ SymmetricMatrix tmp;
+ tmp << cov(data.t());
+ EigenValues(tmp,RD,RE);
+// S = RE * RD * RE.t()
+
+ tmp << sqrtm(svd_X_V * svd_X_D * svd_X_U.t() * RE * RD * RE.t() *svd_X_U * svd_X_D * svd_X_V.t());
+ data = tmp.i()*design.t() * data;
+ }
+
+ if( approach.value()!=string("MLM") && approach.value()!=string("CVA") && approach.value()!=string("PLS") &&
+ approach.value()!=string("SVD-CVA") && approach.value()!=string("orthoPLS"))
+ message(" Using method : PCA" << endl;);
+
+ //perform an SVD on data
+ outMsize(" New Data ", data);
+
+ 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;
+ write_res();
+
+ dbgmsg(" Leaving <do_work>");
+
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+int main(int argc,char *argv[]){
+ Tracer tr("main");
+ OptionParser options(title, examples);
+ try{
+ // must include all wanted options here (the order determines how
+ // the help message is printed)
+ options.add(fnin);
+ options.add(fnout);
+ options.add(approach);
+ options.add(fndesign);
+ options.add(fnmask);
+ options.add(normdes);
+ options.add(perfvn);
+ options.add(perf_demean);
+ options.add(help);
+ options.add(verbose);
+ options.add(debug);
+ options.add(outres);
+ options.add(outdata);
+ options.add(outvnscales);
+
+ options.parse_command_line(argc, argv);
+
+ // line below stops the program if the help was requested or
+ // a compulsory option was not set
+ if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
+ options.usage();
+ exit(EXIT_FAILURE);
+ }else{
+ // Call the local functions
+ return do_work(argc,argv);
+ }
+ }
+ catch(X_OptionError& e){
+ options.usage();
+ cerr << endl << e.what() << endl;
+ exit(EXIT_FAILURE);
+ }
+ catch(std::exception &e){
+ cerr << e.what() << endl;
+ }
+}
+
--
GitLab