ContrastMgr.cc

/*  ContrastMgr.cc

    Mark Woolrich, FMRIB Image Analysis Group

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

/*  CCOPYRIGHT  */

#include <strstream>
#include "ContrastMgr.h"
#include "miscmaths.h"
#include "Log.h"
#include "sigproc.h"
#include "gaussComparer.h"
#include "t2z.h"

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

  ContrastMgr::ContrastMgr(const Matrix& pcontrasts) : 
    contrasts(pcontrasts),
    c(),
    c_counter(0),
    numParams(pcontrasts.Ncols()),
    corrections(),
    b(),
    dof(0.0),
    sigmaSquareds(),
    varcb(),
    cb(),
    neff(),
    numTS(0)
    {}

  void ContrastMgr::Load()
    {
      Tracer ts("ContrastMgr::Load");
      // Need to read in b, sigmaSquareds, corrections and dof 
      Log& logger = Log::getInstance();
     
      // sigmaSquareds:
      sigmaSquareds.read(logger.getDir() + "/sigmasquareds");
      sigmaSquareds.threshold(0.0);
      numTS = sigmaSquareds.getVolumeSize();
     
      // b:
      Volume peVol;
      b.ReSize(numTS, numParams);
      for(int i = 1; i <= numParams; i++)
	{ 
	  // Add param number to "pe" to create filename:
	  char strc[4];
	  ostrstream osc(strc,4);
	  osc << i << '\0';
	  
	  peVol.read(logger.getDir() + "/pe" + strc);

	  peVol.setDims(sigmaSquareds.getDims());
	  peVol.setPreThresholdPositions(sigmaSquareds.getPreThresholdPositions());
	  peVol.threshold();
	  
	  b.Column(i) = peVol; 
	}

      // dof:
      ColumnVector dofVec(1);
      SIGPROC::In(logger.getDir() + "/dof", dofVec);
      dof = dofVec(1);

      // corrections:
      corrections.ReSize(numTS, numParams*numParams);
      SIGPROC::In(logger.getDir() + "/corrections", corrections);
    }

  void ContrastMgr::Save(const string& suffix)
    {
      Tracer ts("ContrastMgr::Save");
      Log& logger = Log::getInstance();

      // prepare contrast number:
      char strc[50];
      ostrstream osc(strc,50);
      osc << suffix << c_counter << '\0';

      // Write out neffs:
      neff.setDims(sigmaSquareds.getDims());
      neff.setPreThresholdPositions(sigmaSquareds.getPreThresholdPositions());
      neff.unthreshold();
      neff.writeAsFloat(logger.getDir() + "/neff" + strc);

      // Write out cope:
      cb.setDims(sigmaSquareds.getDims());
      cb.setPreThresholdPositions(sigmaSquareds.getPreThresholdPositions());
      cb.unthreshold();
      cb.writeAsFloat(logger.getDir() + "/cope" + strc);

      // Write out varcope:
      varcb.setDims(sigmaSquareds.getDims());
      varcb.setPreThresholdPositions(sigmaSquareds.getPreThresholdPositions());
      varcb.unthreshold();
      varcb.writeAsFloat(logger.getDir() + "/varcope" + strc);

      // Write out tstat:
      tstat.setDims(sigmaSquareds.getDims());
      tstat.setPreThresholdPositions(sigmaSquareds.getPreThresholdPositions());
      tstat.unthreshold();
      tstat.writeAsFloat(logger.getDir() + "/tstat" + strc);

      // Write out zstat:
      zstat.setDims(sigmaSquareds.getDims());
      zstat.setPreThresholdPositions(sigmaSquareds.getPreThresholdPositions());
      zstat.unthreshold();
      zstat.writeAsFloat(logger.getDir() + "/zstat" + strc);
    }

  void ContrastMgr::GetCorrection(Matrix& corr, const int ind)
    {
      Tracer ts("ContrastMgr::GetCorrection");

      // puts ColumnVector of length p*p from correction
      // into Matrix corr which is p*p:
      corr.ReSize(numParams, numParams);

      for (int i = 1; i <= numParams; i++)
	{
	  for (int j = 1; j <= numParams; j++)
	    {
	      corr(i,j) = corrections(ind, (i-1)*numParams + j); 
	    }
	}
    }

  void ContrastMgr::ComputeZStat()
    {
      Tracer ts("ContrastMgr::ComputeZStat");

      Log& logger = Log::getInstance();

      // calulate Zstat:
      tstat.ReSize(numTS);
      for(int i = 1; i <= numTS; i++)
	{
	  if(varcb(i) > 0 && neff(i) > 0)
	    {
	      tstat(i) = cb(i)/sqrt(varcb(i));
	    }
	  else
	    tstat(i) = 0.0;
	}
      
      // Calculate tstat:
      zstat.ReSize(numTS);
      T2z::ComputeZStats(varcb, cb, (int)dof, zstat);

      // Compare with theory:
      GaussComparer gaussComp(zstat);
      gaussComp.setup();

      ColumnVector ratios(5);
      ColumnVector probs(5);
      int co = 1;
      for(float p = 0.05; p >= 0.0005; p=p/sqrt(10))
	{
	  float temp = gaussComp.computeRatio(p, logger.str());
	  logger.str() << "p = " << p << ": ratio = " << temp << endl;
	  ratios(co) = temp;
	  probs(co) = p;
	  co++;
	}
      
      logger.out("ratios", ratios);
      logger.out("probs", probs);
      
    }

  void ContrastMgr::ComputeCope()
    { 
      Tracer ts("ContrastMgr::ComputeCope");
      cb.ReSize(numTS);

      for(int i = 1; i <= numTS; i++)
	{	
	  cb(i) = (c.t()*b.Row(i).t()).AsScalar();
	}
    }

  void ContrastMgr::ComputeNeff(bool removeEnds)
    {
      Tracer ts("ContrastMgr::ComputeNeff");
   
      Log& logger = Log::getInstance();
      Matrix corr;
      
      neff.ReSize(numTS);
      neff = 0;

      int from = 1;
      int to = numTS;
      if(removeEnds)
	{
	  from = 1000;
	  to = numTS - 1000;
	}

      int numNegs = 0;
      int maxNeff = 0;

      for(int i = from; i <= to; i++)
	{	
	  GetCorrection(corr, i);

	  neff(i) = 1/(c.t()*corr*c).AsScalar();

	  if(maxNeff < neff(i))
	    maxNeff = (int)neff(i);
	  
	  if(neff(i) < 0.0 && i > 1)
	    {
	      neff(i) = neff(i-1);
	      numNegs++;
	    }  
	}

      // Logging:
      logger.str() << "maxNeff = " << maxNeff << endl;
      logger.str() << "Number of negative neffs calculated = " << numNegs << endl;    
      logger.str() << "mean neff = " << mean(neff) << endl;
      logger.str() << "std neff = " << sqrt(var(neff)) << endl;
    }

  void ContrastMgr::ComputeVarCope()
    { 
      Tracer ts("ContrastMgr::ComputeVarCope");
      varcb.ReSize(numTS);

      for(int i = 1; i <= numTS; i++)
	{
	  if(neff(i) > 0)
	    varcb(i) = sigmaSquareds(i)/neff(i);
	  else
	    varcb(i) = 0;
	}
    }

#ifndef NO_NAMESPACE
}
#endif