fsl_mvlm.cc

/*  fsl_glm -

    Christian F. Beckmann, FMRIB Analysis Group

    Copyright (C) 2006-2013 University of Oxford  */

/*  CCOPYRIGHT  */


//Header & includes
#include <vector>
#include "libvis/miscplot.h"
#include "miscmaths/miscmaths.h"
#include "miscmaths/miscprob.h"
#include "armawrap/newmat.h"
#include "newimage/newimageall.h"
#include "utils/options.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 NEWMAT;
using namespace NEWIMAGE;
using namespace MISCPLOT;
using namespace MISCMATHS;
using namespace Utilities;
using namespace std;

namespace FSL_MVLM {

  //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("\nAuthor: Christian F. Beckmann \nCopyright(C) 2006-2013 University of Oxford\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, 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)"),
                          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> 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);
  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;

  ////////////////////////////////////////////////////////////////////////////

  // 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);
    }
  }

  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(FslImageExists(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(FslImageExists(fnin.value())){//read data
      //input is 3D/4D vol
      volume4D<float> tmpdata;
      read_volume4D(tmpdata,fnin.value());

      // 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(FslImageExists(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(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){
      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") && 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());

      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);

          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);

          // 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>");

    return 0;
  }
}

using namespace FSL_MVLM;

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(nmfdim);
    options.add(nmfitt);
    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;
  }
}