From 809d386866a2d4b53d64ae9d866b13bcaf55e0cc Mon Sep 17 00:00:00 2001
From: Mark Woolrich <woolrich@fmrib.ox.ac.uk>
Date: Wed, 25 Aug 2004 11:06:56 +0000
Subject: [PATCH] fixed bug introduced by switching film to miscmaths
volumeseries
---
Makefile | 4 +-
film_gls.cc | 6 +-
film_gls.ccbak | 216 ---------------------------------
film_gls_res.cc | 10 +-
glim.baccc | 210 --------------------------------
glim.bach | 89 --------------
glmrand.cc | 213 ---------------------------------
glmrand.h | 99 ----------------
ols.ccbak | 238 -------------------------------------
ols.hbak | 96 ---------------
sigprocbak.ccbak | 302 -----------------------------------------------
sigprocbak.hbak | 38 ------
12 files changed, 10 insertions(+), 1511 deletions(-)
delete mode 100755 film_gls.ccbak
delete mode 100755 glim.baccc
delete mode 100755 glim.bach
delete mode 100644 glmrand.cc
delete mode 100644 glmrand.h
delete mode 100755 ols.ccbak
delete mode 100755 ols.hbak
delete mode 100755 sigprocbak.ccbak
delete mode 100755 sigprocbak.hbak
diff --git a/Makefile b/Makefile
index 51902fe..e2533bc 100644
--- a/Makefile
+++ b/Makefile
@@ -15,8 +15,8 @@ OBJS = ContrastMgrOptions.o ContrastMgr.o BandPassOptions.o gaussComparer.o Aut
all: ${XFILES}
-film_ols:${OBJS} film_ols.o
- ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} film_ols.o ${LIBS}
+#film_ols:${OBJS} film_ols.o
+# ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} film_ols.o ${LIBS}
contrast_mgr:${OBJS} contrast_mgr.o
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} contrast_mgr.o ${LIBS}
diff --git a/film_gls.cc b/film_gls.cc
index 9a9a2ff..1430ef6 100644
--- a/film_gls.cc
+++ b/film_gls.cc
@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
for(int i = 1; i <= numTS; i++)
{
glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
- residuals.getSeries(i) = glimGls.getResiduals();
+ residuals.setSeries(glimGls.getResiduals(),i);
}
cout << "Completed" << endl;
@@ -175,7 +175,7 @@ int main(int argc, char *argv[])
co++;
}
- x.getSeries(i) = xw;
+ x.setSeries(xw,i);
glimGls.setData(x.getSeries(i), designmattw, i);
}
@@ -189,7 +189,7 @@ int main(int argc, char *argv[])
}
glimGls.setData(x.getSeries(i), parad.getDesignMatrix(), i);
- residuals.getSeries(i) = glimGls.getResiduals();
+ residuals.setSeries(glimGls.getResiduals(),i);
}
}
diff --git a/film_gls.ccbak b/film_gls.ccbak
deleted file mode 100755
index 7cb5e6d..0000000
--- a/film_gls.ccbak
+++ /dev/null
@@ -1,216 +0,0 @@
-/* film_gls.cc
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include <iostream>
-#include <fstream>
-#include <strstream>
-#define WANT_STREAM
-#define WANT_MATH
-
-#include "newmatap.h"
-#include "newmatio.h"
-#include "VolumeSeries.h"
-#include "Volume.h"
-#include "glim.h"
-#include "sigproc.h"
-#include "miscmaths.h"
-#include "gaussComparer.h"
-#include "Log.h"
-#include "AutoCorrEstimator.h"
-#include "paradigm.h"
-#include "FilmGlsOptions.h"
-#include "glimGls.h"
-
-#include <string>
-
-#ifndef NO_NAMESPACE
-using namespace NEWMAT;
-using namespace SIGPROC;
-using namespace TACO;
-using namespace UTILS;
-#endif
-
-int main(int argc, char *argv[])
-{
- try{
- rand();
- // parse command line to find out directory name for logging:
- ofstream out2;
- FilmGlsOptions& globalopts = FilmGlsOptions::getInstance();
- globalopts.parse_command_line(argc, argv, out2);
-
- // Setup logging:
- Log& logger = Log::getInstance();
- logger.setLogFile("glslogfile");
- logger.establishDir(globalopts.datadir);
-
- // parse command line again to output arguments to logfile
- globalopts.parse_command_line(argc, argv, logger.str());
-
- // load non-temporally filtered data
- VolumeSeries x;
- x.read(globalopts.inputfname);
-
- // if needed output the 12th volume for use later
- Volume epivol;
- if(globalopts.smoothACEst)
- {
- epivol = x.getVolume(12).AsColumn();
- epivol.setDims(x.getDims());
-
- epivol.writeAsInt(logger.getDir() + "/" + globalopts.epifname);
- }
-
- // This also removes the mean from each of the time series:
- x.thresholdSeries(globalopts.thresh, true);
-
- // if needed later also threshold the epi volume
- if(globalopts.smoothACEst)
- {
- epivol.setPreThresholdPositions(x.getPreThresholdPositions());
- epivol.threshold();
- }
-
- int sizeTS = x.getNumVolumes();
- int numTS = x.getNumSeries();
-
- // Load paradigm:
- Paradigm parad;
- parad.load(globalopts.paradigmfname, "", false, sizeTS);
-
- // Sort out detrending:
- if(globalopts.detrend)
- {
- SIGPROC::Detrend(x, false);
- }
-
- if(globalopts.verbose)
- {
- logger.out("Gc", parad.getDesignMatrix());
- }
-
- // Set up OLS GLM for non-whitened data
- Glim glim(x, parad.getDesignMatrix());
-
- cerr << "Computing residuals for non-whitened data... ";
- VolumeSeries& rnotw = glim.ComputeResids();
- cerr << "Completed" << endl;
-
- // Estimate Autocorrelations:
- AutoCorrEstimator acEst(rnotw);
-
- if(globalopts.fitAutoRegressiveModel)
- {
- acEst.fitAutoRegressiveModel();
- if(globalopts.verbose)
- {
- AutoCorrEstimator acEstForLogging(rnotw);
- acEstForLogging.calcRaw();
- logger.out("rawac", acEstForLogging.getEstimates());
- logger.out("autoregac", acEst.getEstimates());
- }
- }
- else
- {
- acEst.calcRaw();
-
- if(globalopts.verbose)
- {
- //logger.out("rawac", acEst.getEstimates());
- }
-
- // Smooth raw estimates:
- if(globalopts.smoothACEst)
- {
- acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, globalopts.epith);
-
- }
-
- // Apply constraints to estimate autocorr:
- acEst.pava();
-
- }
-
- ColumnVector I(sizeTS);
- I = 0;
- I(1) = 1;
-
- // Setup OLS GLM for temporally filtered data:
- int numParams = parad.getDesignMatrix().Ncols();
- GlimGls glimGls(numTS, sizeTS, numParams);
-
- ColumnVector xw(sizeTS);
- ColumnVector xprew(sizeTS);
-
- acEst.setDesignMatrix(parad.getDesignMatrix());
- int co = 1;
-
- // Loop through voxels:
- cerr << "Calculating params..." << endl;
-
- for(int i = 1; i <= numTS; i++)
- {
- // Use autocorr estimate to prewhiten data:
- xprew = x.getSeries(i);
-
- Matrix designmattw;
-
- acEst.preWhiten(xprew, xw, i, designmattw);
-
- if(co > 1000 || i == 5618 || i == 5582)
- {
- co = 1;
- cerr << (float)i/(float)numTS << ",";
- }
- else
- co++;
-
- glimGls.setData(xw, designmattw, i);
- }
- cerr << "Completed" << endl;
-
-
- if(globalopts.verbose)
- {
- // Save E
- VolumeSeries& E = acEst.getE();
- ColumnVector& countLargeE = acEst.getCountLargeE();
-
- logger.out("countLargeE", countLargeE);
- VolumeSeries::Dims dims = x.getDims();
- dims.v = acEst.getZeroPad();
-
- E.setDims(dims);
- E.setPreThresholdPositions(x.getPreThresholdPositions());
- E.unthresholdSeries();
- E.writeAsFloat(logger.getDir() + "/E");
-
- VolumeSeries& threshac = acEst.getEstimates();
- threshac.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
- threshac.writeAsFloat(logger.getDir() + "/threshac");
-
- rnotw.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
- rnotw.writeAsFloat(logger.getDir() + "/rnotw");
- }
-
- // Write out necessary data:
- cerr << "Saving results... ";
- glimGls.Save(x.getDims(), x.getPreThresholdPositions());
- cerr << "Completed" << endl;
- }
- catch(Exception p_excp)
- {
- cerr << p_excp.what() << endl;
- }
- catch(...)
- {
- cerr << "Image error" << endl;
- }
- return 0;
-}
-
diff --git a/film_gls_res.cc b/film_gls_res.cc
index 02dbcf8..1d89f90 100644
--- a/film_gls_res.cc
+++ b/film_gls_res.cc
@@ -1,4 +1,4 @@
-/* film_gls.cc
+/* film_gls_res.cc
Mark Woolrich, FMRIB Image Analysis Group
@@ -118,7 +118,7 @@ int main(int argc, char *argv[])
for(int i = 1; i <= numTS; i++)
{
glimGls.setData(x.getSeries(i), dms[i-1], i);
- residuals.getSeries(i) = glimGls.getResiduals();
+ residuals.setSeries(glimGls.getResiduals());
}
cerr << "Completed" << endl;
@@ -199,7 +199,7 @@ int main(int argc, char *argv[])
else
co++;
- x.getSeries(i) = xw;
+ x.setSeries(xw,i);
dms[i-1] = designmattw;
}
@@ -249,7 +249,7 @@ int main(int argc, char *argv[])
if(globalopts.verbose)
{
- residuals.getSeries(i) = glimGls.getResiduals();
+ residuals.setSeries(glimGls.getResiduals(),i);
}
}
cerr << "Completed" << endl;
@@ -274,7 +274,7 @@ int main(int argc, char *argv[])
for(int i = 1; i <= numTS; i++)
{
- dmsmat.getSeries(i) = dms[i-1].Column(j);
+ dmsmat.setSeries(dms[i-1].Column(j),i);
}
dmsmat.setInfo(vinfo);
dmsmat.setPreThresholdPositions(x.getPreThresholdPositions());
diff --git a/glim.baccc b/glim.baccc
deleted file mode 100755
index e8e87d2..0000000
--- a/glim.baccc
+++ /dev/null
@@ -1,210 +0,0 @@
-/* glim.cc
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include <strstream>
-
-#include "glim.h"
-#include "miscmaths.h"
-#include "Log.h"
-#include "Volume.h"
-
-#ifndef NO_NAMESPACE
-using namespace MISCMATHS;
-using namespace UTILS;
-using namespace TACO;
-namespace SIGPROC {
-#endif
-
- Glim::Glim(const VolumeSeries& p_y, const Matrix& p_x):
- y(p_y),
- x(p_x),
- numTS(p_y.Ncols()),
- sizeTS(p_y.Nrows()),
- numParams(p_x.Ncols()),
- r(sizeTS, numTS, p_y.getDims(), p_y.getPreThresholdPositions()),
- pinv_x(p_x.Ncols(), sizeTS),
- V(sizeTS,sizeTS),
- RV(sizeTS,sizeTS),
- RMat(sizeTS,sizeTS),
- batch_size(BATCHSIZE),
- corrections(numTS, numParams*numParams),
- b(numTS, numParams),
- sigmaSquareds(numTS),
- dofs(numTS),
- traceRVs(numTS)
- {
- }
-
- void Glim::Save()
- {
- // Need to save b, sigmaSquareds, corrections and dof
- Log& logger = Log::getInstance();
-
- // b:
- Volume peVol;
- for(int i = 1; i <= numParams; i++)
- {
- peVol = b.Row(i).AsColumn();
- peVol.setDims(y.getDims());
- peVol.setPreThresholdPositions(y.getPreThresholdPositions());
- peVol.unthreshold();
-
- // Add param number to "pe" to create filename:
- char strc[4];
- ostrstream osc(strc,4);
- osc << i << '\0';
-
- peVol.writeAsFloat(logger.getDir() + "/pe" + strc);
- }
-
- // sigmaSquareds:
- sigmaSquareds.setDims(y.getDims());
- sigmaSquareds.setPreThresholdPositions(y.getPreThresholdPositions());
- sigmaSquareds.unthreshold();
- sigmaSquareds.writeAsFloat(logger.getDir() + "/sigmasquareds");
-
- // dof:
- dofs.setDims(y.getDims());
- dofs.setPreThresholdPositions(y.getPreThresholdPositions());
- dofs.unthreshold();
- dofs.writeAsFloat(logger.getDir() + "/dofs");
-
- // traceRVs:
- traceRVs.setDims(y.getDims());
- traceRVs.setPreThresholdPositions(y.getPreThresholdPositions());
- traceRVs.unthreshold();
- traceRVs.writeAsFloat(logger.getDir() + "/traceRVs");
-
- // corrections:
- logger.out("corrections", corrections);
- }
-
- // Called on entire data set:
- const VolumeSeries& Glim::ComputeResids()
- {
- Tracer ts("ComputeResids");
-
- int batch_pos = 1;
-
- // pinv(x) = inv(x'x)x'
- pinv_x = (x.t()*x).i()*x.t();
-
- // R = I - x*pinv(x)
- Matrix I(sizeTS, sizeTS);
- Identity(I);
-
- RMat = I - x*pinv_x;
-
- while(batch_pos <= numTS)
- {
- if(batch_pos+batch_size - 1 > numTS)
- r.Columns(batch_pos, numTS) = RMat*y.Columns(batch_pos, numTS);
- else
- r.Columns(batch_pos, batch_pos+batch_size-1) = RMat*y.Columns(batch_pos, batch_pos+batch_size-1);
-
- batch_pos += batch_size;
- }
-
- return r;
- }
-
- // Called on entire data set:
- void Glim::ComputePes()
- {
- Tracer ts("ComputePe");
-
- b = pinv_x*y;
- }
-
- void Glim::ConstructV(const ColumnVector& p_vrow)
- {
- Tracer ts("ConstructV");
- V = 0;
-
- for (int i = 1; i <= sizeTS; i++)
- {
- V.SubMatrix(i,i,i,sizeTS) = p_vrow.Rows(1,sizeTS-i+1).t();
- V.SubMatrix(i,i,1,i) = p_vrow.Rows(1,i).Reverse().t();
- }
- }
-
- void Glim::SetVrow(const ColumnVector& p_vrow, const int ind)
- {
- Tracer ts("SetVrow");
-
- ConstructV(p_vrow);
-
- // var/e = inv(x'x)x'*V*x*inv(x'x)
- Matrix corr = pinv_x*V*pinv_x.t();
- SetCorrection(corr, ind);
- RV = RMat*V;
- traceRVs(ind) = Trace(RV);
-
- dofs(ind) = traceRVs(ind)*traceRVs(ind)/Trace(RV*RV);
- }
-
- void Glim::SetGlobalVrow(const ColumnVector& p_vrow)
- {
- Tracer ts("Glim::SetGlobalVrow");
-
- ConstructV(p_vrow);
- RV = RMat*V;
-
- traceRVs = Trace(RV);
-
- dofs = Trace(RV)*Trace(RV)/Trace(RV*RV);
- }
-
-
- void Glim::UseGlobalVrow()
- {
- Tracer ts("Glim::UseGlobalVrow");
-
- // var/e = inv(x'x)x'*V*x*inv(x'x)
- Matrix corr = pinv_x*V*x*(x.t()*x).i();
-
- for(int i = 1; i <= numTS; i++)
- {
- SetCorrection(corr, i);
- }
- }
-
- void Glim::ComputeSigmaSquared(const int ind)
- {
- Tracer ts("Glim::ComputeSigmaSquared");
-
- sigmaSquareds(ind) = (r.Column(ind).t()*r.Column(ind)/traceRVs(ind)).AsScalar();
- }
-
- void Glim::SetCorrection(const Matrix& corr, const int ind)
- {
- Tracer ts("SetCorrection");
-
- // puts Matrix corr which is p*p into Matrix correction
- // as a p*p length row:
- int p = corr.Nrows();
-
- for (int i = 1; i <= p; i++)
- {
- for (int j = 1; j <= p; j++)
- {
- corrections(ind, (i-1)*p + j) = corr(i,j);
- }
- }
- }
-
-#ifndef NO_NAMESPACE
-}
-#endif
-
-
-
-
-
-
-
diff --git a/glim.bach b/glim.bach
deleted file mode 100755
index 41cfe36..0000000
--- a/glim.bach
+++ /dev/null
@@ -1,89 +0,0 @@
-/* glim.h
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include <iostream>
-#include <fstream>
-#define WANT_STREAM
-#define WANT_MATH
-
-#include "newmatap.h"
-#include "newmatio.h"
-#include "VolumeSeries.h"
-#include "Volume.h"
-
-#ifndef NO_NAMESPACE
-using namespace NEWMAT;
-using namespace TACO;
-namespace SIGPROC{
-#endif
-
-#if !defined(__glim_h)
-#define __glim_h
-
-#define FALSE 0
-#define TRUE 1
-
-#define BATCHSIZE 50
-
- class Glim
- {
- public:
- Glim(const VolumeSeries& p_y, const Matrix& p_x);
-
- void Save();
-
- const VolumeSeries& ComputeResids();
- void ComputePes();
- void SetVrow(const ColumnVector& p_vrow, const int ind);
- void SetGlobalVrow(const ColumnVector& p_vrow);
- void ComputeSigmaSquared(const int ind);
- void UseGlobalVrow();
-
- private:
- Glim();
- Glim(const Glim&);
- Glim& operator=(const Glim& p_glim);
-
- void SetCorrection(const Matrix& corr, const int ind);
- void ConstructV(const ColumnVector& p_vrow);
-
- // y = bx + r
- const VolumeSeries& y;
- const Matrix& x;
-
- int numTS;
- int sizeTS;
- int numParams;
-
- VolumeSeries r;
- Matrix pinv_x;
- Matrix V;
- Matrix RV;
- Matrix RMat;
-
- int batch_size;
-
- // Data to be saved:
- Matrix corrections;
- Matrix b;
- Volume sigmaSquareds;
- Volume dofs;
- Volume traceRVs;
-
- };
-
-
-#ifndef NO_NAMESPACE
- }
-#endif
-
-#endif
-
-
-
-
diff --git a/glmrand.cc b/glmrand.cc
deleted file mode 100644
index 96e535b..0000000
--- a/glmrand.cc
+++ /dev/null
@@ -1,213 +0,0 @@
-/* glmrand.cc
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "glmrand.h"
-#include "miscmaths/miscmaths.h"
-#include "utils/log.h"
-#include "histogram.h"
-#include "miscmaths/t2z.h"
-#define __STL_NO_DRAND48
-
-#include <vector.h>
-#include <algo.h>
-
-#ifndef NO_NAMESPACE
-using namespace MISCMATHS;
-using namespace FILM;
-using namespace Utilities;
-namespace FILM {
-#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
-
-
-
-
-
-
-
-
-
diff --git a/glmrand.h b/glmrand.h
deleted file mode 100644
index 32677e1..0000000
--- a/glmrand.h
+++ /dev/null
@@ -1,99 +0,0 @@
-/* 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 "miscmaths/volume.h"
-#include "miscmaths/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
-
-
-
diff --git a/ols.ccbak b/ols.ccbak
deleted file mode 100755
index 0cfe973..0000000
--- a/ols.ccbak
+++ /dev/null
@@ -1,238 +0,0 @@
-/* ols.cc
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "ols.h"
-#include "miscmaths.h"
-#include "Log.h"
-
-#ifndef NO_NAMESPACE
-using namespace MISCMATHS;
-namespace FILM {
-#endif
-
- Ols::Ols(const Matrix& p_y, const Matrix& p_x, const Matrix& p_contrasts):
- y(p_y),
- x(p_x),
- contrasts(p_contrasts),
- numTS(p_y.Ncols()),
- sizeTS(p_y.Nrows()),
- r(sizeTS,numTS),
- pinv_x(p_x.Ncols(), sizeTS),
- var_on_e(0.0),
- cb(numTS),
- b(p_x.Ncols()),
- var(numTS),
- dof(sizeTS - p_x.Ncols()),
- V(sizeTS,sizeTS),
- RV(sizeTS,sizeTS),
- RMat(sizeTS,sizeTS),
- batch_size(BATCHSIZE)
- {
- SetContrast(1);
- }
-
-// void Ols::ConstructV(const ColumnVector& p_vrow)
-// {
-// Tracer ts("ConstructV");
-// V = 0;
-
-// for (int i = 1; i <= sizeTS; i++)
-// {
-// V.SubMatrix(i,i,i,sizeTS) = p_vrow.Rows(1,sizeTS-i+1).t();
-// V.SubMatrix(i,i,1,i) = p_vrow.Rows(1,i).Reverse().t();
-// }
-// }
-
-// float Ols::SetupWithV(const ColumnVector& p_vrow, bool p_justvarone)
-// {
-// Tracer ts("SetupWithV");
-
-// ConstructV(p_vrow);
-
-// // var/e = c'inv(x'x)x'*V*x*inv(x'x)*c
-// var_on_e = (c.t()*pinv_x*V*x*(x.t()*x).i()*c).AsScalar();
-
-// if(!p_justvarone)
-// {
-// // dof = 2*trace(RV)^2/trace(R*V*R*V);
-// RV = RMat*V;
-// dof = Trace(RV)*Trace(RV)/Trace(RV*RV);
-// }
-
-// return dof;
-// }
-
-// float Ols::SetupWithV(const ColumnVector& p_vrow, const ColumnVector& p_kfft, ColumnVector& vrow, bool p_justvarone, const int zeropad)
-// {
-// Tracer ts("SetupWithV");
-
-// // make sure p_vrow is cyclic (even function)
-// //ColumnVector vrow(zeropad);
-// vrow.ReSize(zeropad);
-
-// vrow = 0;
-// vrow.Rows(1,sizeTS/2) = p_vrow.Rows(1,sizeTS/2);
-// vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = p_vrow.Rows(2, sizeTS/2).Reverse();
-
-// // fft vrow
-// ColumnVector fft_real;
-// ColumnVector fft_im;
-// ColumnVector dummy(zeropad);
-// dummy = 0;
-
-// ColumnVector realifft(zeropad);
-
-// FFT(vrow, dummy, fft_real, fft_im);
-
-// FFTI(SP(fft_real, p_kfft), dummy, realifft, dummy);
-
-// vrow = realifft.Rows(1,sizeTS);
-
-// // Normalise vrow:
-// vrow = vrow/vrow(1);
-
-// ConstructV(vrow);
-
-// // var/e = c'inv(x'x)x'*V*x*inv(x'x)*c
-// var_on_e = (c.t()*pinv_x*V*x*(x.t()*x).i()*c).AsScalar();
-
-// if(!p_justvarone)
-// {
-// // dof = 2*trace(RV)^2/trace(R*V*R*V);
-// RV = RMat*V;
-// dof = Trace(RV)*Trace(RV)/Trace(RV*RV);
-// }
-
-// return dof;
-// }
-
- const Matrix& Ols::ComputeResids()
- {
- Tracer ts("ComputeResids");
-
- int batch_pos = 1;
-
- // pinv(x) = inv(x'x)x'
- pinv_x = (x.t()*x).i()*x.t();
-
- // R = I - x*pinv(x)
- Matrix I(sizeTS, sizeTS);
- Identity(I);
-
- RMat = I - x*pinv_x;
-
- while(batch_pos <= numTS)
- {
- if(batch_pos+batch_size - 1 > numTS)
- r.Columns(batch_pos, numTS) = RMat*y.Columns(batch_pos, numTS);
- else
- r.Columns(batch_pos, batch_pos+batch_size-1) = RMat*y.Columns(batch_pos, batch_pos+batch_size-1);
-
- batch_pos += batch_size;
- }
-
- return r;
- }
-
-// const ColumnVector& Ols::Computecb()
-// {
-// Tracer ts("Computecb");
-
-// // cerr << "Computing cbs";
-// int batch_pos = 1;
-
-// while(batch_pos <= numTS)
-// {
-// if(batch_pos+batch_size - 1 > numTS)
-// cb.Rows(batch_pos, numTS) = (c.t()*pinv_x*y.Columns(batch_pos, numTS)).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;
-// // cerr << ".";
-// }
-// //cerr << endl;
-// return cb;
-// }
-
-// const float Ols::Computecb(const int ind)
-// {
-// Tracer ts("Computecb");
-
-// cb(ind) = ((c.t()*pinv_x*y.Column(ind)).t()).AsScalar();
-// return cb(ind);
-// }
-
-// const ColumnVector& Ols::Computeb(const int ind)
-// {
-// Tracer ts("Computeb");
-
-// b = pinv_x*y.Column(ind);
-// return b;
-// }
-
-// const ColumnVector& Ols::ComputeVar()
-// {
-// Tracer ts("ComputeVar");
-
-// // cerr << "Computing Vars";
-// int batch_pos = 1;
-
-// Matrix varmatfull(batch_size, batch_size);
-// ColumnVector vartempfull(batch_size);
-
-// while(batch_pos <= numTS)
-// {
-// if(batch_pos+batch_size - 1 > numTS)
-// {
-// // var = e*var_on_e
-// // e is the estimate of the variance of the timeseries, sigma^2
-// Matrix varmat = (r.Columns(batch_pos, numTS).t()*r.Columns(batch_pos, numTS)/Trace(RV))*var_on_e;
-// ColumnVector vartemp;
-// //getdiag(vartemp, varmat); // obsolete fn
-// vartemp = diag(varmat); // MJ NOTE: new fn
-// var.Rows(batch_pos, numTS) = vartemp;
-// }
-// else
-// {
-// varmatfull = (r.Columns(batch_pos, batch_pos+batch_size-1).t()*r.Columns(batch_pos, batch_pos+batch_size-1)/Trace(RV))*var_on_e;
-// //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;
-// // cerr << ".";
-// }
-
-// // cerr << endl;
-// return var;
-// }
-
-// const float Ols::ComputeVar(const int ind)
-// {
-// Tracer ts("ComputeVar");
-
-// // var = e*var_on_e
-// // e is the estimate of the variance of the timeseries, sigma^2
-// var(ind) = ((r.Column(ind).t()*r.Column(ind)/Trace(RV))*var_on_e).AsScalar();
-// // var(ind) = (r.Column(ind).t()*r.Column(ind)*var_on_e).AsScalar();
-
-// return var(ind);
-// }
-
-#ifndef NO_NAMESPACE
-}
-#endif
-
-
-
-
-
-
-
-
-
diff --git a/ols.hbak b/ols.hbak
deleted file mode 100755
index 6755f72..0000000
--- a/ols.hbak
+++ /dev/null
@@ -1,96 +0,0 @@
-/* ols.h
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#if !defined(__ols_h)
-#define __ols_h
-
-#include <iostream>
-#include <fstream>
-#define WANT_STREAM
-#define WANT_MATH
-
-#include "newmatap.h"
-#include "newmatio.h"
-
-using namespace NEWMAT;
-namespace FILM {
-
-#define BATCHSIZE 50
-
- class Ols
- {
- public:
-
- Ols(const Matrix& p_y, const Matrix& p_x, const Matrix& p_contrasts);
-
- const Matrix& ComputeResids();
- const ColumnVector& ComputeVar();
- const ColumnVector& Computecb();
- const float ComputeVar(const int ind);
- const float Computecb(const int ind);
- const ColumnVector& Computeb(const int ind);
-
- float SetupWithV(const ColumnVector& p_vrow, bool p_justvarone = false);
- float SetupWithV(const ColumnVector& p_vrow, const ColumnVector& p_krowspec, ColumnVector& p_res, bool p_justvarone = false, const int zeropad = 512);
-
- const Matrix& Getr() const {return r;}
-
- const ColumnVector& Getvar() const {return var;}
- const ColumnVector& Getcb() const {return cb;}
- float Getdof() const {return dof;}
- float Getvarone() const {return var_on_e;}
-
- void SetContrast(const int p_num) {c = contrasts.Row(p_num).t();}
- void Setvarone(const float p_varone) {var_on_e = p_varone;}
-
- private:
- Ols();
- Ols(const Ols&);
- Ols& operator=(const Ols& p_ols);
-
- void ConstructV(const ColumnVector& p_vrow);
-
- // y = bx + r
- const Matrix& y;
- const Matrix& x;
-
- // Contrast
- const Matrix& contrasts;
- ColumnVector c;
-
- int numTS;
- int sizeTS;
-
- Matrix r;
-
- Matrix pinv_x;
- float var_on_e;
-
- ColumnVector cb;
-
- ColumnVector b;
-
- ColumnVector var;
- float dof;
-
- // Covariance matrix is sigma^2*V:
- Matrix V;
- Matrix RV;
-
- Matrix RMat;
-
- // batch stuff
- int batch_size;
- };
-
-}
-
-#endif
-
-
-
diff --git a/sigprocbak.ccbak b/sigprocbak.ccbak
deleted file mode 100755
index fed675c..0000000
--- a/sigprocbak.ccbak
+++ /dev/null
@@ -1,302 +0,0 @@
-/* sigproc.cc
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include "ols.h"
-#include "sigproc.h"
-#include "miscmaths.h"
-#include "glm.h"
-
-#ifndef NO_NAMESPACE
-namespace SIGPROC {
- using namespace MISCMATHS;
-#endif
-
- void Out(string p_fname, const Matrix& p_mat)
- {
- ofstream out;
- out.open(p_fname.c_str(), ios::out);
- out << "/Matrix" << endl;
- out << p_mat;
- out.close();
- }
-
- void Out(string p_fname, const ColumnVector& p_mat)
- {
- ofstream out;
- out.open(p_fname.c_str(), ios::out);
- out << p_mat;
- out.close();
- }
-
- void Out(string p_fname, const RowVector& p_mat)
- {
- ofstream out;
- out.open(p_fname.c_str(), ios::out);
- out << p_mat;
- out.close();
- }
-
- void In(string p_fname, ColumnVector& p_mat)
- {
- int size = p_mat.Nrows();
-
- ifstream in;
- in.open(p_fname.c_str(), ios::in);
-
- if(!in)
- {
- string err("Unable to open ");
- err = err + p_fname;
- throw Exception(err.c_str());
- }
-
- for(int i = 1; i <= size; i++)
- {
- in >> p_mat(i);
- }
-
- in.close();
- }
-
- void In(string p_fname, Matrix& p_mat)
- {
- ifstream in;
- in.open(p_fname.c_str(), ios::in);
-
- if(!in)
- {
- cerr << "Unable to open " << p_fname << endl;
- throw;
- }
-
- int numWaves = 0;
- int numPoints = 0;
-
- string str;
-
- while(true)
- {
- in >> str;
- if(str == "/Matrix")
- break;
- else if(str == "/NumWaves")
- {
- in >> numWaves;
- }
- else if(str == "/NumPoints" || str == "/NumContrasts")
- {
- in >> numPoints;
- }
- }
-
- if(numWaves != 0)
- {
- p_mat.ReSize(numPoints, numWaves);
- }
- else
- {
- numWaves = p_mat.Ncols();
- numPoints = p_mat.Nrows();
- }
-
- for(int i = 1; i <= numPoints; i++)
- {
- for(int j = 1; j <= numWaves; j++)
- {
- in >> p_mat(i,j);
- }
- }
-
- in.close();
- }
-
- int EstablishZeroPadding(int sizeTS) {
-
- Tracer tr("EstablishZeroPadding");
-
- return (int)pow(2,ceil(log(sizeTS)/log(2)));
- }
-
- void BandPass(Matrix& p_ts, int lowcut, int highcut)
- {
- Tracer tr("BandPass");
-
- int sizeTS = p_ts.Nrows();
- int numTS = p_ts.Ncols();
-
- int zeropad = EstablishZeroPadding(sizeTS);
- ColumnVector x(zeropad);
- x = 0;
- ColumnVector fft_real;
- ColumnVector fft_im;
- ColumnVector dummy(zeropad);
- ColumnVector dummy2;
- dummy = 0;
- ColumnVector realifft(zeropad);
-
- // convert cutofs to zeropad cutoffs:
- lowcut = (int)(lowcut*zeropad/(float)sizeTS);
- highcut = (int)(highcut*zeropad/(float)sizeTS);
-
- for(int i = 1; i <= numTS; i++)
- {
- x.Rows(1,sizeTS) = p_ts.Column(i);
-
- FFT(x, dummy, fft_real, fft_im);
-
- // Perform cutoff here:
- // everything greater than highcut and lower than lowcut is removed
- if(lowcut > 0 && lowcut < sizeTS){
- fft_real.Rows(1, lowcut) = 0;
- fft_im.Rows(1, lowcut) = 0;
- fft_real.Rows(zeropad - lowcut, zeropad) = 0;
- fft_im.Rows(zeropad - lowcut, zeropad) = 0;
- }
- if(highcut < sizeTS && highcut > 0){
- fft_real.Rows(highcut, zeropad - highcut) = 0;
- fft_im.Rows(highcut, zeropad - highcut) = 0;
- }
- //////////////////////
-
- FFTI(fft_real, fft_im, realifft, dummy2);
- p_ts.Column(i) = realifft.Rows(1,sizeTS);
- }
- }
-
- void AutoCorr(const Matrix& p_ts, Matrix& p_ret, Matrix& fts, int p_zeropad)
- {
- Tracer tr("AutoCorr");
-
- int sizeTS = p_ts.Nrows();
- int numTS = p_ts.Ncols();
-
- ColumnVector x(p_zeropad);
- x = 0;
- ColumnVector fft_real;
- ColumnVector fft_im;
- ColumnVector dummy(p_zeropad);
- ColumnVector dummy2;
- dummy = 0;
- ColumnVector realifft(p_zeropad);
- p_ret.ReSize(sizeTS,numTS);
- fts.ReSize(p_zeropad,numTS);
-
- for(int i = 1; i <= numTS; i++)
- {
- x.Rows(1,sizeTS) = p_ts.Column(i);
- FFT(x, dummy, fft_real, fft_im);
-
- for(int j = 1; j <= p_zeropad; j++)
- {
- // (x+iy)(x-iy) = x^2 + y^2
- fft_real(j) = fft_real(j)*fft_real(j) + fft_im(j)*fft_im(j);
- fft_im(j) = 0;
- fts(j,i) = fft_real(j)*fft_real(j) + fft_im(j)*fft_im(j);
- }
-
- // normalise fts:
- fts.Column(i) = fts.Column(i)/fts.Column(i).Sum();
-
- FFTI(fft_real, fft_im, realifft, dummy2);
-
- float varx = var(ColumnVector(x.Rows(1,sizeTS)));
- p_ret.Column(i) = realifft.Rows(1,sizeTS);
-
- for(int j = 1; j <= sizeTS-1; j++)
- {
- // Correction to make autocorr unbiased and normalised
- p_ret(j,i) = p_ret(j,i)/((sizeTS-j)*varx);
- }
- }
-
- }
-
- int Pacf(const ColumnVector& x, int minorder, int maxorder, ColumnVector& betas)
- {
- Tracer ts("pacf");
- int order = -1;
- int sizeTS = x.Nrows();
-
- // Set c
- Matrix c(1,1);
- c(1,1) = 1;
-
- Glm glm;
-
- for(int i = minorder+1; i <= maxorder+1; i++)
- {
- ColumnVector y = x.Rows(i+1,sizeTS);
-
- // Setup design matrix
- Matrix X(sizeTS-i, i);
- X = 0;
-
- for(int j = 1; j <= i; j++)
- {
- X.Column(j) = x.Rows(i+1-j,sizeTS-j).AsColumn();
- }
-
- glm.setData(y, X, c);
-
- glm.ComputeResids();
- betas = glm.Getb();
-
- if((abs(betas(i)) < (1/sizeTS) + (2/pow(sizeTS,0.5)) && order == -1)
- || i == maxorder+1)
- {
- order = i-1;
- break;
- }
- }
- return order;
- }
-
- void Detrend(Matrix& p_ts, bool p_justremovemean)
- {
- Tracer trace("SIGPROC::Detrend");
-
- cerr << "Detrending...";
- int sizeTS = p_ts.Nrows();
-
- // Set c
- Matrix c(1,1);
- c(1,1) = 1;
-
- // p_ts = b*a + p_ret OLS regression
- Matrix a(sizeTS, 3);
-
- if(p_justremovemean) a.ReSize(sizeTS,1);
-
- // Create a
- for(int i = 1; i <= sizeTS; i++)
- {
- a(i,1) = 1;
- if(!p_justremovemean)
- {
- a(i,2) = (float)i/sizeTS;
- a(i,3) = pow((float)i/sizeTS,2);
- }
- }
-
- // Do OLS for all TS:
- Ols ols(p_ts, a, c);
-
- p_ts = ols.ComputeResids();
-
- cerr << " Completed" << endl;
- }
-
-#ifndef NO_NAMESPACE
-}
-#endif
-
-
-
-
-
-
diff --git a/sigprocbak.hbak b/sigprocbak.hbak
deleted file mode 100755
index 68a1354..0000000
--- a/sigprocbak.hbak
+++ /dev/null
@@ -1,38 +0,0 @@
-/* sigproc.h
-
- Mark Woolrich, FMRIB Image Analysis Group
-
- Copyright (C) 1999-2000 University of Oxford */
-
-/* CCOPYRIGHT */
-
-#include <iostream>
-#include <fstream>
-#define WANT_STREAM
-#define WANT_MATH
-#include <string>
-
-#include "newmatap.h"
-#include "newmatio.h"
-
-#if !defined(__sigproc_h)
-#define __sigproc_h
-
-using namespace NEWMAT;
-namespace SIGPROC{
-
- int EstablishZeroPadding(int num);
- void BandPass(Matrix& p_ts, int lowcut, int highcut);
- void AutoCorr(const Matrix& p_ts, Matrix& p_ret, Matrix& fts, int p_zeropad);
- void Detrend(Matrix& p_ts, bool p_justremovemean = false);
- int Pacf(const ColumnVector& x, int minorder, int maxorder, ColumnVector& betas);
- void Out(string p_fname, const Matrix& out);
- void Out(string p_fname, const ColumnVector& out);
- void Out(string p_fname, const RowVector& out);
- void In(string p_fname, ColumnVector& p_mat);
- void In(string p_fname, Matrix& p_mat);
-
-}
-#endif
-
-
--
GitLab