/*  glim.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#include <strstream>

#include "glim.h"
#include "miscmaths/miscmaths.h"
#include "utils/log.h"
#include "miscmaths/volume.h"

#ifndef NO_NAMESPACE
using namespace MISCMATHS;
using namespace Utilities;

namespace FILM {
#endif

  Glim::Glim(VolumeSeries& p_y, const Matrix& p_x):
    y(p_y),
    x(p_x),
    numTS(p_y.Ncols()),
    sizeTS(p_y.Nrows()),
    numParams(p_x.Ncols()),
    r(sizeTS, numTS, p_y.getInfo(), p_y.getPreThresholdPositions()),
    pinv_x(p_x.Ncols(), sizeTS),
    V(sizeTS,sizeTS),
    RV(sizeTS,sizeTS),
    RMat(sizeTS,sizeTS),
    batch_size(BATCHSIZE),
    corrections(numTS, numParams*numParams),
    b(numTS, numParams),
    sigmaSquareds(numTS),
    dof(sizeTS - p_x.Ncols())
    { 
    }

  void Glim::Save()
    {
      // 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(y.getInfo());
	  peVol.setPreThresholdPositions(y.getPreThresholdPositions());
	  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(y.getInfo());
      sigmaSquareds.setPreThresholdPositions(y.getPreThresholdPositions());
      sigmaSquareds.unthreshold();	
      sigmaSquareds.writeAsFloat(logger.getDir() + "/sigmasquareds");

      // dof:
      ColumnVector dofVec(1);
      dofVec = dof;
      write_ascii_matrix(logger.appendDir("dof"), dofVec);

      // corrections:
      write_ascii_matrix(logger.appendDir("corrections"), corrections);
    }

  // Called on entire data set:
  VolumeSeries& Glim::ComputeResids()
    {
      Tracer ts("ComputeResids");

      int batch_pos = 1;

      // pinv(x) = inv(x'x)x'
      //pinv_x = (x.t()*x).i()*x.t();
      pinv_x = pinv(x);

      // R = I - x*pinv(x)
      Matrix I(sizeTS, sizeTS);
      Identity(I);

      RMat = I - x*pinv_x;
      
      while(batch_pos <= numTS)
	{
	  if(batch_pos+batch_size - 1 > numTS)
	    r.Columns(batch_pos, numTS) = RMat*y.Columns(batch_pos, numTS);
	  else
	    r.Columns(batch_pos, batch_pos+batch_size-1) = RMat*y.Columns(batch_pos, batch_pos+batch_size-1);
	
	  batch_pos += batch_size;
	}
      
      return r;
    }

  // Called on entire data set:
  void Glim::ComputePes()
    { 
      Tracer ts("ComputePe");
      
      b = pinv_x*y;
    }

  void Glim::ConstructV(const ColumnVector& p_vrow)
    {
      Tracer ts("ConstructV");
      V = 0;

      for (int i = 1; i <= sizeTS; i++)
	{
	  V.SubMatrix(i,i,i,sizeTS) = p_vrow.Rows(1,sizeTS-i+1).t();
	  V.SubMatrix(i,i,1,i) = p_vrow.Rows(1,i).Reverse().t();
	}
    }

  void Glim::SetVrow(const ColumnVector& p_vrow, const int ind)
    {
      Tracer ts("SetVrow");
      
      ConstructV(p_vrow);

      // var/e = inv(x'x)x'*V*x*inv(x'x)
      Matrix corr = pinv_x*V*pinv_x.t();
      SetCorrection(corr, ind);
    }

  void Glim::SetGlobalVrow(const ColumnVector& p_vrow)
    {
      Tracer ts("Glim::SetGlobalVrow");

      ConstructV(p_vrow);
      RV = RMat*V;

      dof = Trace(RV)*Trace(RV)/Trace(RV*RV);
    }


  void Glim::UseGlobalVrow()
    {
      Tracer ts("Glim::UseGlobalVrow");

      // var/e = inv(x'x)x'*V*x*inv(x'x)
      Matrix corr = pinv_x*V*pinv_x.t();

      for(int i = 1; i <= numTS; i++)
	{
	  SetCorrection(corr, i);
	}
    }

  void Glim::ComputeSigmaSquared(const int ind)
    {
      Tracer ts("Glim::ComputeSigmaSquared");

      sigmaSquareds(ind) = (r.Column(ind).t()*r.Column(ind)/Trace(RV)).AsScalar();
    }

  void Glim::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(ind, (i-1)*p + j) = corr(i,j); 
	    }
	}
    }

#ifndef NO_NAMESPACE
}
#endif