/*  glmrand.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#include "glmrand.h"
#include "miscmaths.h"
#include "ols.h"
#include "Log.h"
#include "histogram.h"
#include "t2z.h"
#define __STL_NO_DRAND48

#include <vector.h>
#include <algo.h>

#ifndef NO_NAMESPACE
using namespace MISCMATHS;
using namespace TACO;
using namespace UTILS;
namespace SIGPROC {
#endif

  void GlmRand::addData(ColumnVector& p_y, Matrix& p_x)
    {
      Tracer ts("GlmRand::addData");

      yorig = p_y;
      x = &p_x;
      sizeTS = p_y.Nrows();
      randomise();
    }

  void GlmRand::randomise()
  {
    Tracer ts("GlmRand::randomise");

    y.ReSize(sizeTS, numrand+1);
    vector<float> yorigvec;
    
    // put in origy:
    y.getSeries(1) = yorig.AsColumn();
 
    //// void columnVector2Vector(const ColumnVector& cvec, vector<float>& vec)
    ColumnVector cvec = yorig;
    vector<float>& vec = yorigvec;
    for(int j=1; j<=sizeTS; j++)
      {
	vec.push_back(cvec(j));	
      }
    ////////
   
    // put in num randomised versions of yorig:
    for(int i=1; i<=numrand; i++)
      {
	random_shuffle(yorigvec.begin(), yorigvec.end());
	
	//// void vector2ColumnVector(const vector<float>& vec, ColumnVector& cvec)
	vec = yorigvec;
	for(int j=1; j<=sizeTS; j++)
	{
	  cvec(j) = vec[j-1];
	}		  
	////////
	y.getSeries(i+1) = cvec.AsColumn();
      }
   
    ComputeResids();
    Computecb();
    ComputeVar();
    ComputeTStats();
  }

  void GlmRand::ComputeTStats()
  {
    Tracer ts("GlmRand::ComputeTStats");

    datats(datatscount) = cb(1)/sqrt(var(1));

    for(int i=1; i<=numrand; i++)
      {
	randts(((datatscount-1)*numrand)+i) = cb(i+1)/sqrt(var(i+1));
      } 

    datatscount++;
  }

  const Volume& GlmRand::ComputeZStats()
  {
    Tracer ts("GlmRand::ComputeZStats");

    Log::getInstance().out("randts",randts);
    Log::getInstance().out("datats",datats);

    Histogram hist(randts, randts.getVolumeSize());
    hist.generate();

    Volume logprob(numTS);
    float logtotal = log((float)hist.integrateAll());
    cerr << logtotal << endl;

    T2z& t2z = T2z::getInstance();
    for(int i=1; i<=numTS; i++)
      {
	float numtoinf = (float)hist.integrateToInf(datats(i));

	if(!(numtoinf>0))
	  numtoinf = 1;
	
	logprob(i) = log(numtoinf)-logtotal; 
	datazs(i) = t2z.convertlogp2z(logprob(i));
      }
    
    Log::getInstance().out("logprob",logprob);
    Log::getInstance().out("datazs",datazs);

    return datazs;
  }

  void GlmRand::ComputeResids()
    {
      Tracer ts("GlmRand::ComputeResids");

      int batch_pos = 1;
      Matrix& d = *x;

      pinv_x = (d.t()*d).i()*d.t();

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

      RMat = I-d*pinv_x;

      r.ReSize(sizeTS, numrand+1);
      while(batch_pos <= numrand+1)
	{
	  if(batch_pos+batch_size - 1 >  numrand+1)
	    r.Columns(batch_pos, numrand+1) = RMat*y.Columns(batch_pos,  numrand+1);
	  else
	    r.Columns(batch_pos, batch_pos+batch_size-1) = RMat*y.Columns(batch_pos, batch_pos+batch_size-1);
	
	  batch_pos += batch_size;
	}
    }

  void GlmRand::Computecb()
    { 
      Tracer ts("Computecb");
      
      int batch_pos = 1;
      cb.ReSize(numrand+1);
      while(batch_pos <= numrand+1)
	{
	  if(batch_pos+batch_size - 1 > numrand+1)
	    cb.Rows(batch_pos, numrand+1) = (c.t()*pinv_x*y.Columns(batch_pos, numrand+1)).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;
	}
    }

  void GlmRand::ComputeVar()
    { 
      Tracer ts("ComputeVar");

      int batch_pos = 1;
      var.ReSize(numrand+1);
      Matrix varmatfull(batch_size, batch_size);
      ColumnVector vartempfull(batch_size);

      Matrix& d = *x;
      
      // inv_xx = inv(x'x)
      float var_on_e = (c.t()*((d.t()*d).i())*c).AsScalar();

      while(batch_pos <= numrand+1)
	{
	  if(batch_pos+batch_size - 1 > numrand+1)
	    {
	      // var = e*var_on_e
	      // e is the estimate of the variance of the timeseries, sigma^2
	      Matrix varmat = (r.Columns(batch_pos, numrand+1).t()*r.Columns(batch_pos, numrand+1))*var_on_e/sizeTS;
	      ColumnVector vartemp;
	      //getdiag(vartemp, varmat);  // obsolete fn
	      vartemp = diag(varmat);  // MJ NOTE: new fn
	      var.Rows(batch_pos, numrand+1) = vartemp;
	    }      
	  else
	    {
	      varmatfull = (r.Columns(batch_pos, batch_pos+batch_size-1).t()*r.Columns(batch_pos, batch_pos+batch_size-1))*var_on_e/sizeTS;
	      //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;
	}
    }

#ifndef NO_NAMESPACE
}
#endif