Skip to content
Snippets Groups Projects
Commit 809d3868 authored by Mark Woolrich's avatar Mark Woolrich
Browse files

fixed bug introduced by switching film to miscmaths volumeseries

parent c9938709
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Makefile
+ 2
− 2
View file @ 809d3868
......@@ -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}
......
film_gls.cc
+ 3
− 3
View file @ 809d3868
......@@ -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);
}
}
......
film_gls.ccbak deleted 100755 → 0
+ 0
− 216
View file @ c9938709
/* 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;
}
film_gls_res.cc
+ 5
− 5
View file @ 809d3868
/* 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());
......
glim.baccc deleted 100755 → 0
+ 0
− 210
View file @ c9938709
/* 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
glim.bach deleted 100755 → 0
+ 0
− 89
View file @ c9938709
/* 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
glmrand.cc deleted 100644 → 0
+ 0
− 213
View file @ c9938709
/* 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
glmrand.h deleted 100644 → 0
+ 0
− 99
View file @ c9938709
/* 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
ols.ccbak deleted 100755 → 0
+ 0
− 238
View file @ c9938709
/* 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
ols.hbak deleted 100755 → 0
+ 0
− 96
View file @ c9938709
/* 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
sigprocbak.ccbak deleted 100755 → 0
+ 0
− 302
View file @ c9938709
/* 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
sigprocbak.hbak deleted 100755 → 0
+ 0
− 38
View file @ c9938709
/* 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment