diff --git a/Makefile b/Makefile
index 51902fe811bc8029711b9bf6df6c36640e203da1..e2533bc352fcf50ea147e41a5739a3a90ccbbf2d 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 9a9a2ff04be36795f45715a21aee5a6ea355f335..1430ef6117a1bf3410656dbcd672628f7a879e7a 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 7cb5e6d56b45b1b860887e7a846bbb9ed0dfa3b0..0000000000000000000000000000000000000000
--- 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 02dbcf87ca5a342bb5524b6c2cf0ac8fbd764a76..1d89f90999e02a55c47c826f1f943d0f5e4f4b81 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 e8e87d201cfbceb91faadce83169956f053177c5..0000000000000000000000000000000000000000
--- 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 41cfe361a7002a73964743d2380c7347204518fc..0000000000000000000000000000000000000000
--- 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 96e535b058ec8d30ae7c02b188ccf27cd50abe08..0000000000000000000000000000000000000000
--- 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 32677e1bf882d923d7496ae8948390c056c702ce..0000000000000000000000000000000000000000
--- 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 0cfe9730faba4812d91eac5dd217c3dd6a3cf463..0000000000000000000000000000000000000000
--- 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 6755f729ec3dc2e706323e143c6756696c1cf881..0000000000000000000000000000000000000000
--- 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 fed675c4ee315dd83dedb1bf874e2a17f25a0fa0..0000000000000000000000000000000000000000
--- 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 68a13540c7c9bb5e4e5386c931ed85c4740392f7..0000000000000000000000000000000000000000
--- 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
-
-