Skip to content
Snippets Groups Projects
Commit fc7593dd authored by Christian Beckmann's avatar Christian Beckmann
Browse files

added -p option

parent 89457170
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
meldata.cc
+ 82
− 51
View file @ fc7593dd
......@@ -168,21 +168,21 @@ namespace Melodic{
add_Smodes(tmpS2);
}
//add GLM OLS fit
if(Tdes.Storage()){
Matrix alltcs = Tmodes.at(0).Column(1);
for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
alltcs|=Tmodes.at(ctr).Column(1);
if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols()))
glmT.olsfit(alltcs,Tdes,Tcon);
}
if(Sdes.Storage()){
Matrix alltcs = Smodes.at(0);
for(int ctr=1; ctr < (int)Smodes.size();ctr++)
alltcs|=Smodes.at(ctr);
if((alltcs.Nrows()==Sdes.Nrows())&&(Sdes.Nrows()>Sdes.Ncols()&&alltcs.Nrows()>2))
glmS.olsfit(alltcs,Sdes,Scon);
}
//add GLM OLS fit
if(Tdes.Storage()){
Matrix alltcs = Tmodes.at(0).Column(1);
for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
alltcs|=Tmodes.at(ctr).Column(1);
if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols()))
glmT.olsfit(alltcs,Tdes,Tcon);
}
if(Sdes.Storage()){
Matrix alltcs = Smodes.at(0);
for(int ctr=1; ctr < (int)Smodes.size();ctr++)
alltcs|=Smodes.at(ctr);
if((alltcs.Nrows()==Sdes.Nrows())&&(Sdes.Nrows()>Sdes.Ncols()&&alltcs.Nrows()>2))
glmS.olsfit(alltcs,Sdes,Scon);
}
}
void MelodicData::setup()
......@@ -256,9 +256,17 @@ namespace Melodic{
PPCA=tmpPPCA;
}
order = opts.pca_dim.value();
calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
if(opts.debug.value()){
if (opts.paradigmfname.value().length()>0){
Matrix tmpPscales;
tmpPscales = param.t() * alldat;
paramS = stdev(tmpPscales.t());
calc_white(pcaE, pcaD, order, param, paramS, whiteMatrix, dewhiteMatrix);
}else
calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
if(opts.debug.value()){
outMsize("pcaE",pcaE); saveascii(pcaE,"pcaE");
outMsize("pcaD",pcaD); saveascii(pcaD,"pcaD");
outMsize("AdjEV",AdjEV); saveascii(AdjEV,"AdjEV");
......@@ -332,16 +340,16 @@ namespace Melodic{
//create mask
create_mask(Mask);
//setup background image
if(opts.bgimage.value()>""){
read_volume(background,opts.bgimage.value());
if(!samesize(Mean,background)){
cerr << "ERROR:: background image and data have different dimensions \n\n";
exit(2);
}
}else{
background = Mean;
}
//setup background image
if(opts.bgimage.value()>""){
read_volume(background,opts.bgimage.value());
if(!samesize(Mean,background)){
cerr << "ERROR:: background image and data have different dimensions \n\n";
exit(2);
}
}else{
background = Mean;
}
if(!samesize(Mean,Mask)){
cerr << "ERROR:: mask and data have different dimensions \n\n";
......@@ -360,28 +368,46 @@ namespace Melodic{
double tmptime = time(NULL);
srand((unsigned int) tmptime);
//read in post-proc design matrices etc
if(opts.fn_Tdesign.value().length()>0)
Tdes = read_ascii_matrix(opts.fn_Tdesign.value());
if(opts.fn_Sdesign.value().length()>0)
Sdes = read_ascii_matrix(opts.fn_Sdesign.value());
if(opts.fn_Tcon.value().length()>0)
Tcon = read_ascii_matrix(opts.fn_Tcon.value());
if(opts.fn_Scon.value().length()>0)
Scon = read_ascii_matrix(opts.fn_Scon.value());
if(opts.fn_TconF.value().length()>0)
TconF = read_ascii_matrix(opts.fn_TconF.value());
if(opts.fn_SconF.value().length()>0)
SconF = read_ascii_matrix(opts.fn_SconF.value());
if(opts.paradigmfname.value().length()>0){
message(" Use columns in " << opts.paradigmfname.value()
<< " for PCA initialisation" <<endl);
param = read_ascii_matrix(opts.paradigmfname.value());
if(numfiles>1 && Sdes.Storage() == 0){
Sdes = ones(numfiles,1);
if(Scon.Storage() == 0){
Scon = ones(1,1);
Scon &= -1*Scon;
}
Matrix tmpPU, tmpPV;
DiagonalMatrix tmpPD;
SVD(param, tmpPD, tmpPU, tmpPV);
param = tmpPU;
opts.pca_dim.set_T(std::max(opts.pca_dim.value(), param.Ncols()+3));
if(opts.debug.value()){
outMsize("Paradigm",param); saveascii(param,"param");
outMsize("ParadigmS",paramS); saveascii(param,"paramS");
}
Tdes = remmean(Tdes,1);
opts.guessfname.set_T(opts.paradigmfname.value());
}
//read in post-proc design matrices etc
if(opts.fn_Tdesign.value().length()>0)
Tdes = read_ascii_matrix(opts.fn_Tdesign.value());
if(opts.fn_Sdesign.value().length()>0)
Sdes = read_ascii_matrix(opts.fn_Sdesign.value());
if(opts.fn_Tcon.value().length()>0)
Tcon = read_ascii_matrix(opts.fn_Tcon.value());
if(opts.fn_Scon.value().length()>0)
Scon = read_ascii_matrix(opts.fn_Scon.value());
if(opts.fn_TconF.value().length()>0)
TconF = read_ascii_matrix(opts.fn_TconF.value());
if(opts.fn_SconF.value().length()>0)
SconF = read_ascii_matrix(opts.fn_SconF.value());
if(numfiles>1 && Sdes.Storage() == 0){
Sdes = ones(numfiles,1);
if(Scon.Storage() == 0){
Scon = ones(1,1);
Scon &= -1*Scon;
}
}
Tdes = remmean(Tdes,1);
}
void MelodicData::save()
......@@ -585,7 +611,8 @@ namespace Melodic{
RXweight = tmpRX.matrix(Mask);
}
void MelodicData::est_smoothness(){
void MelodicData::est_smoothness()
{
if(Resels == 0){
string SM_path = opts.binpath + "smoothest";
string Mask_fname = logger.appendDir("mask");
......@@ -620,7 +647,8 @@ namespace Melodic{
}
}
unsigned long MelodicData::standardise(volume<float>& mask, volume4D<float>& R){
unsigned long MelodicData::standardise(volume<float>& mask, volume4D<float>& R)
{
unsigned long count = 0;
int M=R.tsize();
......@@ -809,7 +837,8 @@ namespace Melodic{
}
} //void create_mask()
void MelodicData::sort(){
void MelodicData::sort()
{
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
numTime = mixMatrix.Nrows(), i,j;
......@@ -891,11 +920,13 @@ namespace Melodic{
unmixMatrix = pinv(mixMatrix);
}
void MelodicData::reregress(){
void MelodicData::reregress()
{
if((numfiles > 1)){
}
}
void MelodicData::status(const string &txt)
{
cout << "MelodicData Object " << txt << endl;
......
meldata.h
+ 7
− 1
View file @ fc7593dd
......@@ -94,6 +94,12 @@ namespace Melodic{
void set_TSmode();
inline Matrix& get_param() {return param;}
inline void set_param(Matrix& Arg) {param = Arg;}
inline Matrix& get_paramS() {return paramS;}
inline void set_paramS(Matrix& Arg) {paramS = Arg;}
inline Matrix& get_white() {return whiteMatrix;}
inline void set_white(Matrix& Arg) {whiteMatrix = Arg;}
......@@ -194,7 +200,7 @@ namespace Melodic{
volumeinfo tempInfo;
vector<Matrix> DWM, WM;
basicGLM glmT, glmS;
Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF, param, paramS;
RowVector explained_var;
private:
......
melhlprfns.cc
+ 53
− 1
View file @ fc7593dd
......@@ -169,7 +169,7 @@ namespace Melodic{
return Res;
} //Matrix SP
Matrix corrcoef(const Matrix& in1, const Matrix& in2){
Matrix corrcoef(const Matrix& in1, const Matrix& in2){
Matrix tmp = in1;
tmp |= in2;
Matrix out;
......@@ -219,6 +219,50 @@ namespace Melodic{
return Res;
} //Matrix calc_corr
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite)
{
// tmpD2= tmpD | tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols());
// cerr << tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols()) << endl;
//
// Matrix tmpPE;
// tmpPE = SP(param,ones(param.Nrows(),1)*pow(stdev(param,1)*std::sqrt((float)param.Nrows()),-1));
// RE |= tmpPE;
// RowVector tmpD2;
// tmpD2 = tmpD | stdev(param,1).AsRow()*std::sqrt((float)param.Nrows());
// RD << abs(diag(tmpD2.t()));
// RD << RD.SymSubMatrix(N-dim+1,N);
Matrix RE;
DiagonalMatrix RD;
int N = tmpE.Ncols();
dim = std::min(dim,N);
// cerr << stdev(param) << endl;
RE = tmpE.Columns(std::min(N-dim+1+param.Ncols(),N-2),N);
RE |= param;
// cerr << paramS.Nrows() << " x " << paramS.Ncols() << endl;
// cerr << paramS << endl;
RowVector tmpD2;
tmpD2 = tmpD | pow(paramS,2).AsRow();
RD << abs(diag(tmpD2.t()));
// cerr << " " <<tmpD2.Ncols() << " " << N << " " << dim << endl;
RD << RD.SymSubMatrix(N-dim+1+param.Ncols(),N+param.Ncols());
float res = 1.0;
white = sqrt(abs(RD.i()))*RE.t();
dewhite = RE*sqrt(RD);
if(dim < N)
res = PercEV(dim);
return res;
} //Matrix calc_white
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite)
{
Matrix RE;
......@@ -245,6 +289,13 @@ namespace Melodic{
tmp2 = calc_white(tmpE,tmpD, tmp, dim, white, dewhite);
} //Matrix calc_white
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite)
{
RowVector tmp(tmpE.Ncols());
float tmp2;
tmp2 = calc_white(tmpE,tmpD, tmp, dim, param, paramS, white, dewhite);
} //Matrix calc_white
void calc_white(const Matrix& Corr, int dim, Matrix& white, Matrix& dewhite)
{
Matrix RE;
......@@ -257,6 +308,7 @@ namespace Melodic{
calc_white(RE,tmp2, dim, white, dewhite);
} //Matrix calc_white
void std_pca(const Matrix& Mat, const Matrix& weights, Matrix& Corr, Matrix& evecs, RowVector& evals)
{
if(weights.Storage()>0)
......
melhlprfns.h
+ 4
− 2
View file @ fc7593dd
......@@ -29,7 +29,7 @@ namespace Melodic{
Matrix convert_to_pbsc(Matrix& Mat);
RowVector varnorm(Matrix& in, int dim = 30, float level = 1.6);
void varnorm(Matrix& in, const RowVector& vars);
void varnorm(Matrix& in, const RowVector& vars);
RowVector varnorm(Matrix& in, Matrix& Corr, int dim = 30, float level = 1.6);
Matrix SP2(const Matrix& in, const Matrix& weights, bool econ = 0);
......@@ -37,11 +37,13 @@ namespace Melodic{
RowVector Feta(int n1,int n2);
RowVector cumsum(const RowVector& Inp);
Matrix corrcoef(const Matrix& in1, const Matrix& in2);
Matrix corrcoef(const Matrix& in1, const Matrix& in2);
Matrix calc_corr(const Matrix& in, bool econ = 0);
Matrix calc_corr(const Matrix& in, const Matrix& weights, bool econ = 0);
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite);
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& white, Matrix& dewhite);
void calc_white(const Matrix& Corr, int dim, Matrix& white, Matrix& dewhite);
......
melodic.cc
+ 5
− 6
View file @ fc7593dd
......@@ -96,8 +96,8 @@ int main(int argc, char *argv[]){
no_conv = icaobj.no_convergence;
opts.maxNumItt.set_T(500);
if((opts.approach.value()=="symm")&&
(retry > std::min(opts.retrystep,3))){
if((opts.approach.value()=="symm")&&(retry > std::min(opts.retrystep,3)))
{
if(no_conv){
retry++;
opts.approach.set_T("defl");
......@@ -166,7 +166,7 @@ int main(int argc, char *argv[]){
message("finished!" << endl << endl);
}
else {
else {
message(endl <<"No convergence -- giving up " << endl << endl);
}
}
......@@ -277,8 +277,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
mmres = mmall(logger,opts,melodat,report,pmaps);
}
Matrix mmall(Log& logger, MelodicOptions& opts,
MelodicData& melodat, MelodicReport& report, Matrix& pmaps){
Matrix mmall(Log& logger, MelodicOptions& opts,MelodicData& melodat, MelodicReport& report, Matrix& pmaps){
Matrix mmpars(5*melodat.get_IC().Nrows(),5);
mmpars = 0;
......@@ -307,7 +306,7 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
ICmap = SP(melodat.get_IC().Row(ctr),diagvals(ctr)*melodat.get_stdNoisei());
}else{
ICmap = melodat.get_IC().Row(ctr);
}
}
string wherelog;
if(opts.genreport.value())
wherelog = report.getDir();
......
meloptions.h
+ 66
− 66
View file @ fc7593dd
......@@ -155,11 +155,11 @@ class MelodicOptions {
inline MelodicOptions::MelodicOptions() :
logdir(string("-o,--outdir"), string("log.txt"),
string("output directory name\n"),
false, requires_argument),
string("output directory name\n"),
false, requires_argument),
inputfname(string("-i,--in"), std::vector<string>(),
string("input file names (either single file name or comma-separated list or text file)"),
true, requires_argument),
string("input file names (either single file name or comma-separated list or text file)"),
true, requires_argument),
outputfname(string("-O,--out"), string("melodic"),
string("output file name"),
false, requires_argument,false),
......@@ -187,12 +187,12 @@ class MelodicOptions {
joined_whiten(string("--sep_whiten"), true,
string("switch on separate whitening"),
false, no_argument),
joined_vn(string("--sep_vn"), true,
string("switch off joined variance nomalisation"),
false, no_argument),
vn_level(string("--vn_level"), float(2.3),
string("variance nomalisation threshold level (Z> value is ignored)"),
false, requires_argument, false),
joined_vn(string("--sep_vn"), true,
string("switch off joined variance nomalisation"),
false, no_argument),
vn_level(string("--vn_level"), float(2.3),
string("variance nomalisation threshold level (Z> value is ignored)"),
false, requires_argument, false),
numICs(string("-n,--numICs"), -1,
string("numer of IC's to extract (for deflation approach)"),
false, requires_argument),
......@@ -220,9 +220,9 @@ class MelodicOptions {
epsilon(string("--eps,--epsilon"), 0.0005,
string("minimum error change"),
false, requires_argument),
epsilonS(string("--epsS,--epsilonS"), 0.03,
string("minimum error change for rank-1 approximation in TICA"),
false, requires_argument),
epsilonS(string("--epsS,--epsilonS"), 0.03,
string("minimum error change for rank-1 approximation in TICA"),
false, requires_argument),
maxNumItt(string("--maxit"), 500,
string("\tmaximum number of iterations before restart"),
false, requires_argument),
......@@ -247,54 +247,54 @@ class MelodicOptions {
smodename(string("--smode"), string(""),
string("\tmatrix of session modes for report generation"),
false, requires_argument),
filter(string("-f,--filter"), string(""),
string("component numbers to remove\n "),
false, requires_argument),
filter(string("-f,--filter"), string(""),
string("component numbers to remove\n "),
false, requires_argument),
genreport(string("--report"), false,
string("generate Melodic web report"),
false, no_argument),
guireport(string("--guireport"), string(""),
string("modify report for GUI use"),
false, requires_argument, false),
bgimage(string("--bgimage"), string(""),
string("specify background image for report (default: mean image)\n "),
false, requires_argument),
guireport(string("--guireport"), string(""),
string("modify report for GUI use"),
false, requires_argument, false),
bgimage(string("--bgimage"), string(""),
string("specify background image for report (default: mean image)\n "),
false, requires_argument),
tr(string("--tr"), 0.0,
string("\tTR in seconds"),
false, requires_argument),
logPower(string("--logPower"), false,
string("calculate log of power for frequency spectrum\n"),
false, no_argument),
addsigchng(string("--sigchng"), false,
string("add signal change estimates to report pages\n"),
false, no_argument, false),
allPPCA(string("--allPPCA"), false,
string("add all PPCA plots\n"),
false, no_argument, false),
varplots(string("--varplots"), false,
string("add std error envelopes to time course plots\n"),
false, no_argument, false),
varvals(string("--varvals"), false,
string("add rank1 values after plots\n"),
false, no_argument, false),
fn_Tdesign(string("--Tdes"), string(""),
logPower(string("--logPower"), false,
string("calculate log of power for frequency spectrum\n"),
false, no_argument),
addsigchng(string("--sigchng"), false,
string("add signal change estimates to report pages\n"),
false, no_argument, false),
allPPCA(string("--allPPCA"), false,
string("add all PPCA plots\n"),
false, no_argument, false),
varplots(string("--varplots"), false,
string("add std error envelopes to time course plots\n"),
false, no_argument, false),
varvals(string("--varvals"), false,
string("add rank1 values after plots\n"),
false, no_argument, false),
fn_Tdesign(string("--Tdes"), string(""),
string(" design matrix across time-domain"),
false, requires_argument),
fn_Tcon(string("--Tcon"), string(""),
string(" t-contrast matrix across time-domain"),
false, requires_argument),
false, requires_argument),
fn_Tcon(string("--Tcon"), string(""),
string(" t-contrast matrix across time-domain"),
false, requires_argument),
fn_TconF(string("--Tconf"), string(""),
string(" F-contrast matrix across time-domain"),
false, requires_argument, false),
string(" F-contrast matrix across time-domain"),
false, requires_argument, false),
fn_Sdesign(string("--Sdes"), string(""),
string(" design matrix across subject-domain"),
false, requires_argument),
fn_Scon(string("--Scon"), string(""),
string(" t-contrast matrix across subject-domain"),
false, requires_argument),
fn_SconF(string("--Sconf"), string(""),
string(" F-contrast matrix across subject-domain"),
false, requires_argument,false),
string(" design matrix across subject-domain"),
false, requires_argument),
fn_Scon(string("--Scon"), string(""),
string(" t-contrast matrix across subject-domain"),
false, requires_argument),
fn_SconF(string("--Sconf"), string(""),
string(" F-contrast matrix across subject-domain"),
false, requires_argument,false),
output_all(string("--Oall"), false,
string(" output everything"),
false, no_argument),
......@@ -335,7 +335,7 @@ class MelodicOptions {
string("file name of guess of mixing matrix"),
false, requires_argument, false),
paradigmfname(string("-p,--paradigm"), string(""),
string("file name of FEAT paradigm file"),
string("file name of FEAT paradigm file (design.mat)"),
false, requires_argument, false),
dummy(string("--dummy"), 0,
string("number of dummy volumes"),
......@@ -355,21 +355,21 @@ class MelodicOptions {
remove_meanvol(string("--keep_meanvol"), true,
string("do not subtract mean volume"),
false, no_argument, false),
remove_meantc(string("--keep_meantc"), true,
string("do not remove mean time course"),
false, no_argument, false),
remove_meantc(string("--remove_meantc"), false,
string("remove mean time course"),
false, no_argument, false),
remove_endslices(string("--remEndslices"), false,
string("delete end slices (motion correction artefacts)"),
false, no_argument,false),
rescale_nht(string("--rescale_nht"), true,
string("switch off map rescaling after mixture-modelling"),
false, no_argument,false),
guess_remderiv(string("--remderiv"), false,
guess_remderiv(string("--remove_deriv"), false,
string("removes every second entry in paradigm file (EV derivatives)"),
false, no_argument, false),
false, no_argument),
temporal(string("--temporal"), false,
string("perform temporal ICA"),
false, no_argument, false),
false, no_argument, false),
retrystep(3),
options(title, usageexmpl)
{
......@@ -409,19 +409,19 @@ class MelodicOptions {
options.add(filter);
options.add(genreport);
options.add(guireport);
options.add(bgimage);
options.add(bgimage);
options.add(tr);
options.add(logPower);
options.add(addsigchng);
options.add(allPPCA);
options.add(varplots);
options.add(varvals);
options.add(fn_Tdesign);
options.add(fn_Tcon);
options.add(fn_TconF);
options.add(fn_Sdesign);
options.add(fn_Scon);
options.add(fn_SconF);
options.add(fn_Tdesign);
options.add(fn_Tcon);
options.add(fn_TconF);
options.add(fn_Sdesign);
options.add(fn_Scon);
options.add(fn_SconF);
options.add(output_all);
options.add(output_unmix);
options.add(output_MMstats);
......
melpca.cc
+ 18
− 3
View file @ fc7593dd
......@@ -32,7 +32,17 @@ namespace Melodic{
Matrix tmpE;
RowVector tmpD, AdjEV, PercEV;
std_pca(in,weights,Corr,tmpE,tmpD);
if(opts.paradigmfname.value().length()>0)
{
basicGLM tmpglm;
tmpglm.olsfit(in,melodat.get_param(),Identity(melodat.get_param().Ncols()));
std_pca(tmpglm.get_residu(),weights,Corr,tmpE,tmpD);
// std_pca(in,weights,Corr,tmpE,tmpD,melodat.get_param());
}
else{
std_pca(in,weights,Corr,tmpE,tmpD);
}
if(opts.tsmooth.value()){
message(endl << " temporal smoothing of Eigenvectors " << endl);
tmpE=smoothColumns(tmpE);
......@@ -80,8 +90,13 @@ namespace Melodic{
Matrix tmpDeWhite;
float varp = 1.0;
varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
melodat.get_EVP() ,opts.pca_dim.value(),tmpWhite,tmpDeWhite);
if (opts.paradigmfname.value().length()>0)
varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
melodat.get_EVP(),opts.pca_dim.value(),melodat.get_param(),melodat.get_paramS(),tmpWhite,tmpDeWhite);
else
varp = calc_white(melodat.get_pcaE(),melodat.get_pcaD(),
melodat.get_EVP(),opts.pca_dim.value(),tmpWhite,tmpDeWhite);
melodat.set_white(tmpWhite);
melodat.set_dewhite(tmpDeWhite);
......
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