/*  glm.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#include "glm.h"
#include "miscmaths/miscmaths.h"
#include "utils/log.h"

using namespace MISCMATHS;

namespace FILM {

  void Glm::setData(Matrix& p_y, Matrix& p_x, Matrix& p_contrasts)
    {
      y = &p_y;
      x = &p_x;
      numTS = p_y.Ncols();
      sizeTS = p_y.Nrows();
      contrasts = &p_contrasts;
      SetContrast(1);
    }
  
  const Matrix& Glm::ComputeResids()
    {
      Tracer ts("ComputeResids");

      Matrix& d = *x;
      Computeb();

      // r = (I - x(x'x)-1x')y
      Matrix I(sizeTS, sizeTS);
      Identity(I);      
      r = (I-d*inv_xx*d.t())*(y->Column(1));
	
      return r;
    }

  const float Glm::Computecb()
    { 
      Tracer ts("Computecb");

      cb = (c.t()*b).AsScalar();
      return cb;
    }

  const ColumnVector& Glm::Computeb()
    {
      Tracer ts("Computeb");

      Matrix& d = *x;
      
      // inv_xx = inv(x'x)
      // inv_xx = (d.t()*d).i();
      inv_xx = pinv(d.t()*d);

      b = inv_xx*d.t()*(y->Column(1));

      return b;
    }

  const float Glm::ComputeVar()
    {
      Tracer ts("ComputeVar");
   
      // var = e*var_on_e
      // e is the estimate of the variance of the timeseries, sigma^2
      var = (r.Column(1).t()*r.Column(1)*c.t()*inv_xx*c).AsScalar()/sizeTS;

      return var;
    }

#ifndef NO_NAMESPACE
}
#endif