glmrand.cc

/*  glmrand.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  Part of FSL - FMRIB's Software Library
    WWW:      http://www.fmrib.ox.ac.uk/fsl
    Email:    fsl@fmrib.ox.ac.uk
    
    Developed at FMRIB (Oxford Centre for Functional Magnetic Resonance
    Imaging of the Brain), Department of Clinical Neurology, Oxford
    University, Oxford, UK
    
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
(at
    your option) any later version.
    
    This program is distributed in the hope that it will be useful, but
in
    order that the University as a charitable foundation protects its
    assets for the benefit of its educational and research purposes, the
    University makes clear that no condition is made or to be implied,
nor
    is any warranty given or to be implied, as to the accuracy of FSL,
or
    that it will be suitable for any particular purpose or for use under
    any specific conditions.  Furthermore, the University disclaims all
    responsibility for the use which is made of FSL.  See the GNU
General
    Public License for more details.
    
    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
    USA */

#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);
	      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);
	      var.Rows(batch_pos, batch_pos+batch_size-1) = vartempfull;
	    }
	  batch_pos += batch_size;
	}
    }

#ifndef NO_NAMESPACE
}
#endif