ols.ccbak

/*  ols.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#include "ols.h"
#include "miscmaths.h"
#include "Log.h"

#ifndef NO_NAMESPACE
using namespace MISCMATHS;
namespace FILM {
#endif

  Ols::Ols(const Matrix& p_y, const Matrix& p_x, const Matrix& p_contrasts):
    y(p_y),
    x(p_x),
    contrasts(p_contrasts),
    numTS(p_y.Ncols()),
    sizeTS(p_y.Nrows()),
    r(sizeTS,numTS),
    pinv_x(p_x.Ncols(), sizeTS),
    var_on_e(0.0),
    cb(numTS),
    b(p_x.Ncols()),
    var(numTS),
    dof(sizeTS - p_x.Ncols()),
    V(sizeTS,sizeTS),
    RV(sizeTS,sizeTS),
    RMat(sizeTS,sizeTS),
    batch_size(BATCHSIZE)
    {
      SetContrast(1);
    }
  
//    void Ols::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();
//         }
//      }

//    float Ols::SetupWithV(const ColumnVector& p_vrow, bool p_justvarone)
//      {
//        Tracer ts("SetupWithV");
      
//        ConstructV(p_vrow);

//        // var/e = c'inv(x'x)x'*V*x*inv(x'x)*c
//        var_on_e = (c.t()*pinv_x*V*x*(x.t()*x).i()*c).AsScalar();

//        if(!p_justvarone)
//  	{
//  	  // dof = 2*trace(RV)^2/trace(R*V*R*V);
//  	  RV = RMat*V;
//  	  dof = Trace(RV)*Trace(RV)/Trace(RV*RV);
//  	}

//        return dof;
//      }

//    float Ols::SetupWithV(const ColumnVector& p_vrow, const ColumnVector& p_kfft, ColumnVector& vrow, bool p_justvarone, const int zeropad)
//      {
//        Tracer ts("SetupWithV");

//        // make sure p_vrow is cyclic (even function)
//        //ColumnVector vrow(zeropad);
//        vrow.ReSize(zeropad);

//        vrow = 0;
//        vrow.Rows(1,sizeTS/2) = p_vrow.Rows(1,sizeTS/2);
//        vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = p_vrow.Rows(2, sizeTS/2).Reverse();

//        // fft vrow
//        ColumnVector fft_real;
//        ColumnVector fft_im;
//        ColumnVector dummy(zeropad);
//        dummy = 0;
      
//        ColumnVector realifft(zeropad);

//        FFT(vrow, dummy, fft_real, fft_im);

//        FFTI(SP(fft_real, p_kfft), dummy, realifft, dummy);
      
//        vrow = realifft.Rows(1,sizeTS);

//        // Normalise vrow:
//        vrow = vrow/vrow(1);

//        ConstructV(vrow);     

//        // var/e = c'inv(x'x)x'*V*x*inv(x'x)*c
//        var_on_e = (c.t()*pinv_x*V*x*(x.t()*x).i()*c).AsScalar();
      
//        if(!p_justvarone)
//  	{
//  	  // dof = 2*trace(RV)^2/trace(R*V*R*V);
//  	  RV = RMat*V;
//  	  dof = Trace(RV)*Trace(RV)/Trace(RV*RV);
//  	}

//        return dof;
//      }

  const Matrix& Ols::ComputeResids()
    {
      Tracer ts("ComputeResids");

      int batch_pos = 1;

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

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

//    const ColumnVector& Ols::Computecb()
//      { 
//        Tracer ts("Computecb");
      
//        //     cerr << "Computing cbs";
//        int batch_pos = 1;
      
//        while(batch_pos <= numTS)
//  	{
//  	  if(batch_pos+batch_size - 1 > numTS)
//  	    cb.Rows(batch_pos, numTS) = (c.t()*pinv_x*y.Columns(batch_pos, numTS)).t();
//  	  else
//  	    cb.Rows(batch_pos, batch_pos+batch_size-1) = (c.t()*pinv_x*y.Columns(batch_pos, batch_pos+batch_size-1)).t();
//  	  batch_pos += batch_size;
//  	  //	  cerr << ".";
//  	}
//        //cerr << endl;
//        return cb;
//      }

//    const float Ols::Computecb(const int ind)
//      { 
//        Tracer ts("Computecb");

//        cb(ind) = ((c.t()*pinv_x*y.Column(ind)).t()).AsScalar();
//        return cb(ind);
//      }

//    const ColumnVector& Ols::Computeb(const int ind)
//      { 
//        Tracer ts("Computeb");
      
//        b = pinv_x*y.Column(ind);
//        return b;
//      }

//    const ColumnVector& Ols::ComputeVar()
//      { 
//        Tracer ts("ComputeVar");

//        //      cerr << "Computing Vars";
//        int batch_pos = 1;

//        Matrix varmatfull(batch_size, batch_size);
//        ColumnVector vartempfull(batch_size);

//        while(batch_pos <= numTS)
//  	{
//  	  if(batch_pos+batch_size - 1 > numTS)
//  	    {
//  	      // var = e*var_on_e
//  	      // e is the estimate of the variance of the timeseries, sigma^2
//  	      Matrix varmat = (r.Columns(batch_pos, numTS).t()*r.Columns(batch_pos, numTS)/Trace(RV))*var_on_e;
//  	      ColumnVector vartemp;
//  	      //getdiag(vartemp, varmat);  // obsolete fn
//  	      vartemp = diag(varmat);  // MJ NOTE: new fn
//  	      var.Rows(batch_pos, numTS) = vartemp;
//  	    }      
//  	  else
//  	    {
//  	      varmatfull = (r.Columns(batch_pos, batch_pos+batch_size-1).t()*r.Columns(batch_pos, batch_pos+batch_size-1)/Trace(RV))*var_on_e;
//  	      //getdiag(vartempfull, varmatfull);  // obsolete fn
//  	      vartempfull = diag(varmatfull);  // MJ NOTE: new fn
//  	      var.Rows(batch_pos, batch_pos+batch_size-1) = vartempfull;
//  	    }
//  	  batch_pos += batch_size;
//  	  //  cerr << ".";
//  	}
      
//        // cerr << endl;
//        return var;
//      }

//    const float Ols::ComputeVar(const int ind)
//      { 
//        Tracer ts("ComputeVar");
   
//        // var = e*var_on_e
//        // e is the estimate of the variance of the timeseries, sigma^2
//        var(ind) = ((r.Column(ind).t()*r.Column(ind)/Trace(RV))*var_on_e).AsScalar();
//        // var(ind) = (r.Column(ind).t()*r.Column(ind)*var_on_e).AsScalar();

//        return var(ind);
//      }

#ifndef NO_NAMESPACE
}
#endif