glimGls.cc

/*  glimGls.cc

    Mark Woolrich, FMRIB Image Analysis Group

    Copyright (C) 1999-2000 University of Oxford  */

/*  CCOPYRIGHT  */

#include <sstream>

#include "glimGls.h"
#include "miscmaths/miscmaths.h"
#include "utils/log.h"
#include "miscmaths/volume.h"
#include "miscmaths/volumeseries.h"

using namespace MISCMATHS;
using namespace Utilities;

namespace FILM {

  GlimGls::GlimGls(const int pnumTS, const int psizeTS, const int pnumParams) : 
    numTS(pnumTS),
    sizeTS(psizeTS),
    numParams(pnumParams),
    corrections(numParams*numParams,numTS),
    b(numParams, numTS),
    sigmaSquareds(numTS),
    dof(sizeTS - numParams)
  {
    I.ReSize(sizeTS);
    for(int i=1; i<=sizeTS; i++)
      {
	I(i,i)=1;
      }
  }
  
  void GlimGls::CleanUp()
  {
    corrections.CleanUp();
    sigmaSquareds.CleanUp();
    b.CleanUp();
    r.CleanUp();
  }

  void GlimGls::setData(const ColumnVector& y, const Matrix& x, const int ind)
    {
      // compute b
      //Matrix inv_xx = (x.t()*x).i();
      Matrix inv_xx = pinv(x.t()*x);
      b.Column(ind) = inv_xx*x.t()*y;

      // compute r
      Matrix R = I-x*inv_xx*x.t();
      r = R*y;

      // compute sigma squareds 
      //      sigmaSquareds(ind) = (r.t()*r/sizeTS).AsScalar();
      sigmaSquareds(ind) = (r.t()*r).AsScalar()/R.Trace();

      // set corrections
      SetCorrection(inv_xx, ind);
    }

  void GlimGls::Save(const VolumeInfo& pvolinfo, const ColumnVector& prethreshpos)
    {
      // Need to save b, sigmaSquareds, corrections and dof 
      Log& logger = LogSingleton::getInstance();

      VolumeInfo volinfo = pvolinfo;

      // b:
      Volume peVol;
      for(int i = 1; i <= numParams; i++)
	{
	  peVol = b.Row(i).AsColumn();
	  volinfo.intent_code = NIFTI_INTENT_ESTIMATE;
	  peVol.setInfo(volinfo);
	  peVol.setPreThresholdPositions(prethreshpos);
	  peVol.unthreshold();
	  
	  // Add param number to "pe" to create filename:
	  ostringstream osc;
	  osc << i;
	  
	  peVol.writeAsFloat(logger.getDir() + "/pe" + osc.str().c_str());
	}

      // sigmaSquareds:
      volinfo.intent_code = NIFTI_INTENT_ESTIMATE;
      sigmaSquareds.setInfo(volinfo);
      sigmaSquareds.setPreThresholdPositions(prethreshpos);
      sigmaSquareds.unthreshold();	
      sigmaSquareds.writeAsFloat(logger.getDir() + "/sigmasquareds");

      // dof:
      ColumnVector dofVec(1);
      dofVec = dof;
      write_ascii_matrix(logger.appendDir("dof"), dofVec);
      
      // corrections (are the pes correlation matrix (x.t()*x).i() reshapen to a vector):
      VolumeInfo newvolinfo = volinfo;
      newvolinfo.v = numParams*numParams;
      newvolinfo.intent_code = NIFTI_INTENT_NONE;
      //corrections.setInfo(newvolinfo);
      //corrections.setPreThresholdPositions(prethreshpos);
      //corrections.unthresholdSeries();
      //corrections.writeAsFloat(logger.getDir() + "/corrections");
      corrections.writeThresholdedSeriesAsFloat(newvolinfo,prethreshpos,logger.getDir() + "/corrections");
    }

  void GlimGls::SetCorrection(const Matrix& corr, const int ind)
    {
      Tracer ts("SetCorrection");

      // puts Matrix corr which is p*p into Matrix correction
      // as a p*p length row:
      int p = corr.Nrows();
      
      for (int i = 1; i <= p; i++)
	{
	  for (int j = 1; j <= p; j++)
	    {
	      corrections((i-1)*p + j, ind) = corr(i,j); 
	    }
	}
    }

}