glimGls.cc

/*  glimGls.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#include <strstream>

#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 {

  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();
      b.Column(ind) = inv_xx*x.t()*y;

      // compute r
      Matrix I(sizeTS, sizeTS);
      Identity(I);
      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& volinfo, const ColumnVector& prethreshpos)
    {
      // Need to save b, sigmaSquareds, corrections and dof 
      Log& logger = LogSingleton::getInstance();

      // b:
      Volume peVol;
      for(int i = 1; i <= numParams; i++)
	{
	  peVol = b.Row(i).AsColumn();
	  peVol.setInfo(volinfo);
	  peVol.setPreThresholdPositions(prethreshpos);
	  peVol.unthreshold();
	  
	  // Add param number to "pe" to create filename:
	  char strc[4];
	  ostrstream osc(strc,4);
	  osc << i << '\0';
	  
	  peVol.writeAsFloat(logger.getDir() + "/pe" + strc);
	}

      // sigmaSquareds:
      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;
      corrections.setInfo(newvolinfo);
      corrections.setPreThresholdPositions(prethreshpos);
      corrections.unthresholdSeries();
      corrections.writeAsFloat(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); 
	    }
	}
    }

}