/*  glmrand.h

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#if !defined(__glmRand_h)
#define __glmRand_h

#include <iostream>
#include <fstream>
#define WANT_STREAM
#define WANT_MATH

#include "newmatap.h"
#include "newmatio.h"
#include "Volume.h"
#include "VolumeSeries.h"

using namespace NEWMAT;

namespace FILM{

#define BATCHSIZE 50

  class GlmRand
    {
    public:
     
      GlmRand(const int pnumrand, const int pnumts, Matrix& pcontrasts, const int contrastnum) : 
	yorig(),
	x(),
	contrasts(&pcontrasts),
	numTS(pnumts),
	numrand(pnumrand),
	randts(pnumrand*pnumts),
	datats(pnumts),
	datazs(pnumts),
	datatscount(1),
	batch_size(BATCHSIZE) {SetContrast(contrastnum);} 

      void addData(ColumnVector& p_y, Matrix& p_x);

      void SetContrast(const int p_num) {c = contrasts->Row(p_num).t();}
      const Volume& ComputeZStats();

    private:
      GlmRand(const GlmRand&);
      GlmRand& operator=(const GlmRand& p_glmRand);

      void randomise();
      void ComputeResids();
      void ComputeVar();
      void Computecb();
      void ComputeTStats();
      
      // y = bx + r
       
      // Original time series:
      ColumnVector yorig;

      VolumeSeries y;
      Matrix* x;
   
      // Contrast
      Matrix* contrasts;
      ColumnVector c;

      int numTS;
      int sizeTS;
      int numrand;

      VolumeSeries r;

      Matrix pinv_x;

      Volume cb;
      Volume var;

      Volume randts;
      Volume datats;
      Volume datazs;

      Matrix RMat;

      int datatscount;

      // batch stuff
      int batch_size;
    };

}

#endif