diff --git a/meldata.cc b/meldata.cc
index 128050dfe2434e308c973547f9cf74a90e882939..28569991919c879677d3b9bf9a8ecdef28768ea3 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -568,7 +568,7 @@ namespace Melodic{
dbgmsg(string("START: setup") << endl);
numfiles = (int)opts.inputfname.value().size();
- setup_misc();
+ setup_misc();
if(opts.debug.value())
memmsg(" after setup_misc ");
@@ -576,12 +576,14 @@ namespace Melodic{
Data = process_file(opts.inputfname.value().at(0));
}
else{
- if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
- opts.approach.set_T("concat");
- if(opts.migp.value())
- setup_migp();
- else
- setup_classic();
+ if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
+ opts.approach.set_T("concat");
+ if (opts.approach.value()==string("tica"))
+ opts.migp.set_T(false);
+ if( opts.approach.value()==string("concat") && opts.migp.value() )
+ setup_migp();
+ else
+ setup_classic();
}
message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
diff --git a/melhlprfns.cc b/melhlprfns.cc
new file mode 100644
index 0000000000000000000000000000000000000000..07c32a43a753e15d8a7024322365b30269833847
--- /dev/null
+++ b/melhlprfns.cc
@@ -0,0 +1,912 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melhlprfns.cc - misc functions
+
+ Christian F. Beckmann, FMRIB Analysis Group
+
+ Copyright (C) 1999-2013 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "melhlprfns.h"
+#include "libprob.h"
+#include "miscmaths/miscprob.h"
+#include "miscmaths/t2z.h"
+#include "miscmaths/f2z.h"
+
+namespace Melodic{
+
+ void update_mask(volume<float>& mask, Matrix& Data)
+ {
+ Matrix DStDev=stdev(Data);
+ volume4D<float> tmpMask, RawData;
+ tmpMask.setmatrix(DStDev,mask);
+
+ float tMmax;
+ volume<float> tmpMask2;
+ tmpMask2 = tmpMask[0];
+
+ tMmax = tmpMask2.max();
+ double st_mean = DStDev.Sum()/DStDev.Ncols();
+ double st_std = stdev(DStDev.t()).AsScalar();
+
+ volume4D<float> newmask;
+ newmask = binarise(tmpMask2,(float) max((float) st_mean-3*st_std,(float) 0.01*st_mean),tMmax);
+
+ Matrix newmaskM,newData;
+ newmaskM = newmask.matrix(mask);
+ int N = Data.Nrows();
+
+ if(int(newmaskM.Row(1).SumAbsoluteValue() + 0.3) < newmaskM.Ncols()){
+ RawData.setmatrix(Data.Row(1),mask);
+ newData = RawData.matrix(newmask[0]);
+ for(int r=2; r <= N; r++){
+ RawData.setmatrix(Data.Row(r),mask);
+ newData &= RawData.matrix(newmask[0]);
+ }
+ Data = newData;
+ mask = newmask[0];
+ }
+ }
+
+ void del_vols(volume4D<float>& in, int howmany)
+ {
+ for(int ctr=1; ctr<=howmany; ctr++){
+ in.deletevolume(ctr);
+ }
+ }
+
+ Matrix calc_FFT(const Matrix& Mat, const bool logpwr)
+ {
+ Matrix res;
+ for(int ctr=1; ctr <= Mat.Ncols(); ctr++){
+ ColumnVector tmpCol;
+ tmpCol=Mat.Column(ctr);
+ ColumnVector FtmpCol_real;
+ ColumnVector FtmpCol_imag;
+ ColumnVector tmpPow;
+ if(tmpCol.Nrows()%2 != 0){
+ Matrix empty(1,1); empty=0;
+ tmpCol &= empty;
+ }
+ RealFFT(tmpCol,FtmpCol_real,FtmpCol_imag);
+ tmpPow = pow(FtmpCol_real,2)+pow(FtmpCol_imag,2);
+ tmpPow = tmpPow.Rows(2,tmpPow.Nrows());
+ if(logpwr) tmpPow = log(tmpPow);
+ if(res.Storage()==0){res= tmpPow;}else{res|=tmpPow;}
+ }
+ return res;
+ } //Matrix calc_FFT()
+
+ Matrix smoothColumns(const Matrix& inp)
+ {
+ Matrix temp(inp);
+ int ctr1 = temp.Nrows();
+ Matrix temp2(temp);
+ temp2=0;
+
+ temp = temp.Row(4) & temp.Row(3) & temp.Row(2) & temp & temp.Row(ctr1-1)
+ & temp.Row(ctr1-2) &temp.Row(ctr1-3);
+
+ double kern[] ={0.0045 , 0.055, 0.25, 0.4, 0.25, 0.055, 0.0045};
+ double fac = 0.9090909;
+
+
+ for(int cc=1;cc<=temp2.Ncols();cc++){
+ for(int cr=1;cr<=temp2.Nrows();cr++){
+ temp2(cr,cc) = fac*( kern[0] * temp(cr,cc) + kern[1] * temp(cr+1,cc) +
+ kern[2] * temp(cr+2,cc) + kern[3] * temp(cr+3,cc) +
+ kern[4] * temp(cr+4,cc) + kern[5] * temp(cr+5,cc) +
+ kern[6] * temp(cr+6,cc));
+ }
+ }
+ return temp2;
+ } //Matrix smoothColumns()
+
+ Matrix convert_to_pbsc(Matrix& inp)
+ {
+ Matrix meanimg;
+ meanimg = mean(inp);
+ float eps = 0.00001;
+
+ for(int ctr=1; ctr <= inp.Ncols(); ctr++){
+ if(meanimg(1,ctr) < eps)
+ meanimg(1,ctr) = eps;
+ }
+
+ for(int ctr=1; ctr <= inp.Nrows(); ctr++){
+ Matrix tmp;
+ tmp << inp.Row(ctr);
+ inp.Row(ctr) << 100 * SP((tmp - meanimg),pow(meanimg,-1));
+ }
+
+ inp = remmean(inp);
+ return meanimg;
+ } //void convert_to_pbsc
+
+ RowVector varnorm(Matrix& in, int dim, float level, int econ)
+ {
+ SymmetricMatrix Corr(cov_r(in,false,econ));
+ RowVector out;
+ out = varnorm(in,Corr,dim,level, econ);
+ return out;
+ } //RowVector varnorm
+
+ void varnorm(Matrix& in, const RowVector& vars)
+ {
+ for(int ctr=1; ctr <=in.Nrows();ctr++)
+ in.Row(ctr) = SD(in.Row(ctr),vars);
+ }
+
+ RowVector varnorm(Matrix& in, SymmetricMatrix& Corr, int dim, float level, int econ)
+ {
+
+ Matrix tmpE, white, dewhite;
+ RowVector tmpD, tmpD2;
+
+ std_pca(remmean(in,2), Corr, tmpE, tmpD, econ);
+ calc_white(tmpE,tmpD, dim, white, dewhite);
+
+ Matrix ws = white * in;
+ for(int ctr1 = 1; ctr1<=ws.Ncols(); ctr1++)
+ for(int ctr2 = 1; ctr2<=ws.Nrows(); ctr2++)
+ if(std::abs(ws(ctr2,ctr1)) < level)
+ ws(ctr2,ctr1)=0;
+ tmpD = stdev(in - dewhite*ws);
+ for(int ctr1 = 1; ctr1<=tmpD.Ncols(); ctr1++)
+ if(tmpD(ctr1) < 0.01){
+ tmpD(ctr1) = 1.0;
+ in.Column(ctr1) = 0.0*in.Column(ctr1);
+ }
+ varnorm(in,tmpD);
+
+ return tmpD;
+ } //RowVector varnorm
+
+
+
+ Matrix SP2(const Matrix& in, const Matrix& weights, int econ)
+ {
+ Matrix Res;
+ Res = in;
+ if(econ>0){
+ ColumnVector tmp;
+ for(int ctr=1; ctr <= in.Ncols(); ctr++){
+ tmp = in.Column(ctr);
+ tmp = tmp * weights(1,ctr);
+ Res.Column(ctr) = tmp;
+ }
+ }
+ else{
+ Res = ones(in.Nrows(),1)*weights.Row(1);
+ Res = SP(in,Res);
+ }
+ return Res;
+ } //Matrix SP2
+
+ void SP3(Matrix& in, const Matrix& weights)
+ {
+ for(int ctr=1; ctr <= in.Nrows(); ctr++){
+ in.Row(ctr) << SP(in.Row(ctr),weights.AsRow());
+ }
+ }
+
+ Matrix corrcoef(const Matrix& in1, const Matrix& in2){
+ Matrix tmp = in1;
+ tmp |= in2;
+ Matrix out;
+ out=MISCMATHS::corrcoef(tmp,0);
+ return out.SubMatrix(1,in1.Ncols(),in1.Ncols()+1,out.Ncols());
+ }
+
+ Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part){
+ Matrix tmp1 = in1, tmp2 = in2, out;
+ if(part.Storage()){
+ tmp1 = tmp1 - part * pinv(part) * tmp1;
+ tmp2 = tmp2 - part * pinv(part) * tmp2;
+ }
+
+ out = corrcoef(tmp1,tmp2);
+ return out;
+ }
+
+ 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;
+ DiagonalMatrix RD;
+ int N = tmpE.Ncols();
+ dim = std::min(dim,N);
+ RE = tmpE.Columns(N-dim+1,N);
+ RD << abs(diag(tmpD.t()));
+ RD << RD.SymSubMatrix(N-dim+1,N);
+
+ 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
+
+ void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& white, Matrix& dewhite)
+ {
+ RowVector tmp(tmpE.Ncols());
+ float tmp2;
+ 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 SymmetricMatrix& Corr, int dim, Matrix& white, Matrix& dewhite)
+ {
+ Matrix RE;
+ DiagonalMatrix RD;
+ RowVector tmp2;
+ EigenValues(Corr,RD,RE);
+ tmp2 = diag(RD).t();
+ calc_white(RE,tmp2, dim, white, dewhite);
+ } //Matrix calc_white
+
+
+ void std_pca(const Matrix& Mat, const Matrix& weights, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ)
+ {
+ if(weights.Storage()>0)
+ Corr = cov_r(Mat, weights, econ);
+ else
+ Corr = cov_r(Mat,false,econ);
+
+ DiagonalMatrix tmpD;
+ EigenValues(Corr,tmpD,evecs);
+ evals = tmpD.AsRow();
+ } //void std_pca
+
+ void std_pca(const Matrix& Mat, SymmetricMatrix& Corr, Matrix& evecs, RowVector& evals, int econ)
+ {
+ Matrix weights;
+ std_pca(Mat,weights,Corr,evecs,evals, econ);
+ } //void std_pca
+
+ void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc, int iter)
+ {
+ Matrix guess;
+ guess = normrnd(Mat.Nrows(),num_pc);
+ em_pca(Mat,guess,evecs,evals,num_pc,iter);
+ } //void em_pca
+
+ void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc, int iter)
+ {
+ Matrix C;
+ if(guess.Ncols() < num_pc){
+ C=normrnd(Mat.Nrows(),num_pc);
+ C.Columns(1,guess.Ncols()) = guess;
+ }
+ else
+ C = guess;
+
+ Matrix tmp, tmp2;
+ for(int ctr=1; ctr <= iter; ctr++){
+ // E-Step
+ tmp = C.t()*C;
+ tmp = tmp.i();
+ tmp = tmp * C.t();
+ tmp = tmp * Mat;
+ // M-Step
+ tmp2 = tmp * tmp.t();
+ tmp2 = tmp2.i();
+ tmp2 = Mat*tmp.t()*tmp2;
+ C = tmp2;
+ }
+
+ symm_orth(C);
+ Matrix Evc;
+ SymmetricMatrix tmpC;
+ RowVector Evl;
+ tmp = C.t() * Mat;
+ std_pca(tmp,tmpC,Evc,Evl);
+ evals = Evl;
+ evecs = C * Evc;
+ } //void em_pca
+
+ float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
+ {
+ SymmetricMatrix Corr;
+ Matrix Evecs, tmpWM, tmpDWM, tmp;
+ RowVector Evals;
+ std_pca(Mat.t(), Corr, Evecs, Evals);
+ calc_white(Corr, dim, tmpWM, tmpDWM);
+ tmp = tmpWM * Mat.t();
+ cols = tmp.t();
+ rows << tmpDWM;
+ float res;
+ Evals=fliplr(Evals);
+ res = sum(Evals.Columns(1,dim),2).AsScalar()/sum(Evals,2).AsScalar()*100;
+ return res;
+ } // rankapprox
+
+ RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows)
+ {
+ Matrix out; RowVector res(Mat.Ncols());
+ for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++){
+ Matrix tmpVals(cols.Nrows(),rows.Nrows());
+ for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
+ tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() *
+ (ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1);
+
+ Matrix tmpcols, tmprows;
+ res(ctr1) =rankapprox(tmpVals, tmpcols, tmprows);
+ cols.Column(ctr1) = tmpcols;
+ rows.Column(ctr1) = tmprows;
+ }
+ return res;
+ } // krfact
+
+ RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows)
+ {
+ RowVector res;
+ cols = zeros(colnum,Mat.Ncols());
+ rows = zeros(int(Mat.Nrows() / colnum),Mat.Ncols());
+ res = krfact(Mat,cols,rows);
+ return res;
+ } // krfact
+
+ Matrix krprod(const Matrix& cols, const Matrix& rows)
+ {
+ Matrix out;
+ out = zeros(cols.Nrows()*rows.Nrows(),cols.Ncols());
+ for(int ctr1 = 1; ctr1 <= cols.Ncols(); ctr1++)
+ for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
+ {
+ out.SubMatrix(cols.Nrows()*(ctr2-1)+1,cols.Nrows()*ctr2,ctr1,ctr1) << cols.Column(ctr1) * rows(ctr2,ctr1);
+ }
+ return out;
+ } // krprod
+
+ Matrix krapprox(const Matrix& Mat, int size_cols, int dim)
+ {
+ Matrix out, cols, rows;
+ out = zeros(Mat.Nrows(), Mat.Ncols());
+ cols = zeros(size_cols,Mat.Ncols());
+ rows = zeros(int(Mat.Nrows() / size_cols), Mat.Ncols());
+ krfact(Mat,cols,rows);
+ out = krprod(cols, rows);
+ return out;
+ } // krapprox
+
+ void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2, RowVector& out3, int& out4, int num_vox, float resels)
+ {
+ RowVector AdjEV;
+ AdjEV << in.AsRow();
+ AdjEV = AdjEV.Columns(3,AdjEV.Ncols());
+ AdjEV = AdjEV.Reverse();
+
+ RowVector CircleLaw;
+ int NumVox = (int) floor(num_vox/(2.5*resels));
+
+ CircleLaw = Feta(int(AdjEV.Ncols()), NumVox);
+
+ for(int ctr=1;ctr<=CircleLaw.Ncols(); ctr++){
+ if(CircleLaw(ctr)<5*10e-10){CircleLaw(ctr) = 5*10e-10;}
+ }
+
+ /* float slope;
+ slope = CircleLaw.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
+ int(AdjEV.Ncols()/4)).Sum() /
+ AdjEV.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
+ int(AdjEV.Ncols()/4)).Sum();*/
+
+ RowVector PercEV(AdjEV);
+ PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
+
+ AdjEV << SP(AdjEV,pow(CircleLaw.Columns(1,AdjEV.Ncols()),-1));
+
+ SortDescending(AdjEV);
+ int maxEV = 1;
+ float threshold = 0.98;
+ for(int ctr_i = 1; ctr_i < PercEV.Ncols(); ctr_i++){
+ if((PercEV(ctr_i)<threshold)&&(PercEV(ctr_i+1)>=threshold)){maxEV=ctr_i;}
+ }
+
+ if(maxEV<3){maxEV=PercEV.Ncols()/2;}
+ RowVector NewEV;
+ Matrix temp1;
+ temp1 = abs(AdjEV);
+ NewEV << temp1.SubMatrix(1,1,1,maxEV);
+
+ AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
+
+ out1 = AdjEV;
+ out2 = PercEV;
+ out3 = NewEV;
+ out4 = maxEV;
+ } //adj_eigspec
+
+ void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2)
+ {
+ RowVector AdjEV, PercEV;
+ AdjEV = in.Reverse();
+ SortDescending(AdjEV);
+
+ PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
+ AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
+ out1 = AdjEV;
+ out2 = PercEV;
+ } //adj_eigspec
+
+ RowVector Feta(int n1, int n2)
+ {
+ float nu = (float) n1/n2;
+ float bm = pow((1-sqrt(nu)),2.0);
+ float bp = pow((1+sqrt(nu)),2.0);
+
+ float lrange = 0.9*bm;
+ float urange = 1.1*bp;
+
+ RowVector eta(30*n1);
+ float rangestepsize = (urange - lrange) / eta.Ncols();
+ for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
+ eta(ctr_i) = lrange + rangestepsize * (ctr_i);
+ }
+
+ RowVector teta(10*n1);
+ teta = 0;
+ float stepsize = (bp - bm) / teta.Ncols();
+ for(int ctr_i = 1; ctr_i <= teta.Ncols(); ctr_i++){
+ teta(ctr_i) = stepsize*(ctr_i);
+ }
+
+ RowVector feta(teta);
+ feta = SP(pow(2*M_PI*nu*(teta + bm),-1), pow(SP(teta, bp-bm-teta),0.5));
+
+ teta = teta + bm;
+
+ RowVector claw(eta);
+ claw = 0;
+ double tmpval(0);
+ int ctr_j(1);
+ for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
+ for(; ctr_j <= teta.Ncols() && teta(ctr_j) < eta(ctr_i) ; ctr_j++)
+ tmpval += feta(ctr_j);
+ claw(ctr_i) = max(n1*(1-stepsize*tmpval),0.0);
+ }
+
+ RowVector Res(n1); //invert the CDF
+ Res = 0;
+ for(int ctr_i = 1; ctr_i < eta.Ncols(); ctr_i++){
+ if(floor(claw(ctr_i))>floor(claw(ctr_i+1))){
+ Res(int(floor(claw(ctr_i)))) = eta(ctr_i);
+ }
+ }
+
+ return Res;
+ } //RowVector Feta
+
+ RowVector cumsum(const RowVector& Inp)
+ {
+ UpperTriangularMatrix UT(Inp.Ncols());
+ UT=1.0;
+ RowVector Res;
+ Res = Inp * UT;
+ return Res;
+ } //RowVector cumsum
+
+ int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, SymmetricMatrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which)
+ {
+ std_pca(in,weights,Corr,tmpE,tmpD);
+
+ int maxEV = 1;
+ RowVector NewEV;
+ adj_eigspec(tmpD.AsRow(),AdjEV,PercEV,NewEV,maxEV,in.Ncols(),resels);
+
+ int res;
+ PPCA = ppca_est(NewEV, in.Ncols(),resels);
+ ColumnVector tmp = ppca_select(PPCA, res, maxEV, which);
+
+ PPCA = tmp | PPCA;
+ return res;
+ } //int ppca_dim
+
+ int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, float resels, string which)
+ {
+ RowVector tmpD;
+ Matrix tmpE;
+ SymmetricMatrix Corr;
+
+ int res = ppca_dim(in, weights, PPCA, AdjEV, PercEV, Corr, tmpE, tmpD, resels, which);
+ return res;
+ } //int ppca_dim
+
+ int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which)
+ {
+ ColumnVector PPCA;
+ RowVector AdjEV, PercEV;
+ int res = ppca_dim(in,weights,PPCA,AdjEV,PercEV,resels,which);
+ return res;
+ } //int ppca_dim
+
+ ColumnVector ppca_select(Matrix& PPCAest, int& dim, int maxEV, string which)
+ {
+ RowVector estimators(5);
+ estimators = 1.0;
+
+ for(int ctr=1; ctr<=PPCAest.Ncols(); ctr++){
+ PPCAest.Column(ctr) = (PPCAest.Column(ctr) -
+ min(PPCAest.Column(ctr)).AsScalar()) /
+ ( max(PPCAest.Column(ctr)).AsScalar() -
+ min(PPCAest.Column(ctr)).AsScalar());
+ }
+
+ int ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,2) < PPCAest(ctr_i+1,2))&&(ctr_i<maxEV))
+ {estimators(1)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,3) < PPCAest(ctr_i+1,3))&&(ctr_i<maxEV))
+ {estimators(2)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,4) < PPCAest(ctr_i+1,4))&&(ctr_i<maxEV))
+ {estimators(3)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,5) < PPCAest(ctr_i+1,5))&&(ctr_i<maxEV))
+ {estimators(4)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,6) < PPCAest(ctr_i+1,6))&&(ctr_i<maxEV))
+ {estimators(5)=ctr_i+1;ctr_i++;}
+
+ int res = 0;
+ ColumnVector PPCA;
+ RowVector PercEV(PPCAest.Column(1).t());
+ PercEV = cumsum(PercEV / sum(PercEV,2).AsScalar());
+
+ if(which == string("aut")) {
+ if(int(estimators(2)) < int(estimators(1)) &&
+ float(PercEV(int(estimators(2))))>0.8){
+ res=int(estimators(2));
+ PPCA << PPCAest.Column(3);
+ }else{
+ res = int(estimators(1));
+ PPCA << PPCAest.Column(2);
+ }
+ }
+ if(which == string("lap")){
+ res = int(estimators(1));
+ PPCA << PPCAest.Column(2);
+ }
+ if(which == string("bic")){
+ res = int(estimators(2));
+ PPCA << PPCAest.Column(3);
+ }
+ if(which == string("mdl")){
+ res = int(estimators(3));
+ PPCA << PPCAest.Column(4);
+ }
+ if(which == string("rrn")){
+ res = int(estimators(4));
+ PPCA << PPCAest.Column(5);
+ }
+ if(which == string("aic")){
+ res = int(estimators(5));
+ PPCA << PPCAest.Column(6);
+ }
+ if(which == string("median")){
+ RowVector unsorted(estimators);
+ SortAscending(unsorted);
+ ctr_i=1;
+ res=int(unsorted(3));
+ while(res != int(estimators(ctr_i)))
+ ctr_i++;
+ PPCA << PPCAest.Column(ctr_i);
+ }
+ if(res==0 || which == string("mean")){//median estimator
+ PPCA = mean(PPCAest.Columns(2,6),2);
+ res=int(mean(estimators,2).AsScalar());
+ }
+
+ dim = res;
+ return PPCA;
+ } //RowVector ppca_select
+
+ Matrix ppca_est(const RowVector& eigenvalues, const int N1, const float N2)
+ {
+ Matrix Res;
+ Res = ppca_est(eigenvalues, (int) floor(N1/(2.5*N2)));
+ return Res;
+ } //Matrix ppca_est
+
+ Matrix ppca_est(const RowVector& eigenvalues, const int N)
+ {
+ RowVector logLambda(eigenvalues);
+ logLambda = log(eigenvalues);
+
+ int d = logLambda.Ncols();
+
+ RowVector k(d);
+ for(int ctr_i = 1; ctr_i <=d; ctr_i++){
+ k(ctr_i)=ctr_i;
+ }
+
+ RowVector m(d);
+ m=d*k-0.5*SP(k,k+1);
+
+ RowVector loggam(d);
+ loggam = 0.5*k.Reverse();
+ for(int ctr_i = 1; ctr_i <=d; ctr_i++){
+ loggam(ctr_i)=lgam(loggam(ctr_i));
+ }
+ loggam = cumsum(loggam);
+
+ RowVector l_probU(d);
+ l_probU = -log(2)*k + loggam - cumsum(0.5*log(M_PI)*k.Reverse());
+
+ RowVector tmp1;
+ tmp1 = -cumsum(eigenvalues).Reverse()+sum(eigenvalues,2).AsScalar();
+ tmp1(1) = 0.95*tmp1(2);
+ tmp1=tmp1.Reverse();
+
+ RowVector tmp2;
+ tmp2 = -cumsum(logLambda).Reverse()+sum(logLambda,2).AsScalar();
+ tmp2(1)=tmp2(2);
+ tmp2=tmp2.Reverse();
+
+ RowVector tmp3;
+ tmp3 = d-k;
+ tmp3(d) = 1.0;
+
+ RowVector tmp4;
+ tmp4 = SP(tmp1,pow(tmp3,-1));
+ for(int ctr_i = 1; ctr_i <=d; ctr_i++){
+ if(tmp4(ctr_i)<0.01){tmp4(ctr_i)=0.01;}
+ if(tmp3(ctr_i)<0.01){tmp3(ctr_i)=0.01;}
+ if(tmp1(ctr_i)<0.01){tmp1(ctr_i)=0.01;}
+ }
+
+ RowVector l_nu;
+ l_nu = SP(-N/2*(d-k),log(tmp4));
+ l_nu(d) = 0;
+
+ RowVector l_lam;
+ l_lam = -(N/2)*cumsum(logLambda);
+
+ RowVector l_lhood;
+ l_lhood = SP(pow(tmp3,-1),tmp2)-log(SP(pow(tmp3,-1),tmp1));
+
+ Matrix t1,t2, t3;
+ UpperTriangularMatrix triu(d);
+ triu = 1.0;
+ for(int ctr_i = 1; ctr_i <= triu.Ncols(); ctr_i++){
+ triu(ctr_i,ctr_i)=0.0;
+ }
+ t1 = (ones(d,1) * eigenvalues);
+ t1 = SP(triu,t1.t() - t1);
+ t2 = pow(tmp4,-1).t()*ones(1,d);
+ t3 = ones(d,1)*pow(eigenvalues,-1);
+ t2 = SP(triu, t2.t()-t3.t());
+ for(int ctr_i = 1; ctr_i <= t1.Ncols(); ctr_i++){
+ for(int ctr_j = 1; ctr_j <= t1.Nrows(); ctr_j++){
+ if(t1(ctr_j,ctr_i)<=0){t1(ctr_j,ctr_i)=1;}
+ }
+ }
+ for(int ctr_i = 1; ctr_i <= t2.Ncols(); ctr_i++){
+ for(int ctr_j = 1; ctr_j <= t2.Nrows(); ctr_j++){
+ if(t2(ctr_j,ctr_i)<=0){t2(ctr_j,ctr_i)=1;}
+ }
+ }
+ t1 = cumsum(sum(log(t1),2).AsRow());
+ t2 = cumsum(sum(log(t2),2).AsRow());
+
+ RowVector l_Az(d);
+ l_Az << (t1+t2);
+
+ RowVector l_lap;
+ l_lap = l_probU + l_nu +l_Az + l_lam + 0.5*log(2*M_PI)*(m+k)-0.5*log(N)*k;
+
+ RowVector l_BIC;
+ l_BIC = l_lam + l_nu - 0.5*log(N)*(m+k);
+
+ RowVector l_RRN;
+ l_RRN = -0.5*N*SP(k,log(SP(cumsum(eigenvalues),pow(k,-1))))+l_nu;
+
+ RowVector l_AIC;
+ l_AIC = -2*N*SP(tmp3,l_lhood)+ 2*(1+d*k+0.5*(k-1));
+ l_AIC = -l_AIC;
+
+ RowVector l_MDL;
+ l_MDL = -N*SP(tmp3,l_lhood)+ 0.5*(1+d*k+0.5*(k-1))*log(N);
+ l_MDL = -l_MDL;
+
+ Matrix Res;
+
+ Res = eigenvalues.t();
+ Res |= l_lap.t();
+ Res |= l_BIC.t();
+ Res |= l_MDL.t();
+ Res |= l_RRN.t();
+ Res |= l_AIC.t();
+
+
+ return Res;
+ } //Matrix ppca_est
+
+ ColumnVector acf(const ColumnVector& in, int order)
+ {
+ double meanval;
+ meanval = mean(in).AsScalar();
+ int tpoints = in.Nrows();
+
+ ColumnVector y, res;
+ Matrix X, tmp;
+
+ y = in.Rows(order+1,tpoints) - meanval;
+ X = zeros(order + 1, order);
+ for(int ctr1 = 1; ctr1 <= order; ctr1++)
+ X.Column(ctr1) = in.Rows(order + 1 - ctr1, tpoints - ctr1) - meanval;
+ tmp = X.t()*X;
+ tmp = tmp.i();
+ tmp = tmp * X.t();
+ res << tmp * y;
+ return res;
+ } //ColumnVector acf
+
+ ColumnVector pacf(const ColumnVector& in, int maxorder)
+ {
+ int tpoint = in.Nrows();
+ ColumnVector res;
+ res = acf(in, maxorder);
+ for(int ctr1 = 1; ctr1 <= maxorder; ctr1++)
+ if ( res.Column(ctr1).AsScalar() < (1/tpoint) + 2/(float)std::pow(tpoint,0.5))
+ res.Column(ctr1) = 0;
+ return res;
+ } //ColumnVector pacf
+
+ Matrix est_ar(const Matrix& Mat, int maxorder)
+ {
+ Matrix res;
+ res = zeros(maxorder, Mat.Ncols());
+ ColumnVector tmp;
+ for (int ctr = 1; ctr <= Mat.Ncols(); ctr++){
+ tmp = pacf(Mat.Column(ctr));
+ res.Column(ctr) = tmp;
+ }
+ return res;
+ } //Matrix est_ar
+
+ ColumnVector gen_ar(const ColumnVector& in, int maxorder)
+ {
+ float sdev;
+ sdev = stdev(in).AsScalar();
+ ColumnVector x, arcoeff, scaled;
+ scaled = in / sdev;
+ arcoeff = pacf( scaled, maxorder);
+ x = normrnd(in.Nrows(),1).AsColumn() * sdev;
+ for(int ctr1=2; ctr1 <= in.Nrows(); ctr1++)
+ for(int ctr2 = 1; ctr2 <= maxorder; ctr2++)
+ x(ctr1) = arcoeff(ctr2) * x(std::max(1,int(ctr1-ctr2))) + x(ctr1);
+ return x;
+ } //ColumnVector gen_ar
+
+ Matrix gen_ar(const Matrix& in, int maxorder)
+ {
+ Matrix res;
+ res = in;
+ ColumnVector tmp;
+ for(int ctr=1; ctr <= in.Ncols(); ctr++){
+ tmp = in.Column(ctr);
+ res.Column(ctr) = gen_ar(tmp, maxorder);
+ }
+ return res;
+ } //Matrix gen_ar
+
+ Matrix gen_arCorr(const Matrix& in, int maxorder)
+ {
+ Matrix res;
+ res = zeros(in.Nrows(), in.Nrows());
+ ColumnVector tmp;
+ for(int ctr=1; ctr<= in.Ncols(); ctr++){
+ tmp = in.Column(ctr);
+ tmp = gen_ar(tmp, maxorder);
+ res += tmp * tmp.t();
+ }
+ return res;
+ } //Matrix gen_arCorr
+
+ void basicGLM::olsfit(const Matrix& data, const Matrix& design,
+ const Matrix& contrasts, int requestedDOF)
+ {
+ beta = zeros(design.Ncols(),1);
+ residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1);
+ t = zeros(1); z = zeros(1); p=-1.0*ones(1);
+ dof = (int)-1; cbeta = -1.0*ones(1);
+
+ if(data.Nrows()==design.Nrows()){
+ Matrix dat = data;
+ Matrix tmp = design.t()*design;
+ Matrix pinvdes = tmp.i()*design.t();
+
+ beta = pinvdes * dat;
+ residu = dat - design*beta;
+ dof = ols_dof(design);
+ if ( requestedDOF>0)
+ dof = requestedDOF;
+ sigsq = sum(SP(residu,residu))/dof;
+
+ float fact = float(dof) / design.Ncols();
+ f_fmf = SP(sum(SP(design*beta,design*beta)),pow(sum(SP(residu,residu)),-1)) * fact;
+
+ pf_fmf = f_fmf.Row(1);
+ for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
+ pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),dof,f_fmf.Column(ctr1).AsScalar());
+
+ if(contrasts.Storage()>0 && contrasts.Ncols()==beta.Nrows()){
+ cbeta = contrasts*beta;
+ Matrix tmp = contrasts*pinvdes*pinvdes.t()*contrasts.t();
+ varcb = diag(tmp)*sigsq;
+ t = SP(cbeta,pow(varcb,-0.5));
+ z = t; p=t;
+ for(int ctr1=1;ctr1<=t.Ncols();ctr1++){
+ ColumnVector tmp = t.Column(ctr1);
+ T2z::ComputeZStats(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
+ z.Column(ctr1) << tmp;
+ T2z::ComputePs(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
+ p.Column(ctr1) << exp(tmp);
+ }
+ }
+ }
+
+ }
+
+
+}
diff --git a/melodic.h b/melodic.h
new file mode 100644
index 0000000000000000000000000000000000000000..6721b42328a51c7dcb8bcfbb80ed8b3190a039a3
--- /dev/null
+++ b/melodic.h
@@ -0,0 +1,88 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melodic.h - main program header
+
+ Christian F. Beckmann and Matthew Webster, FMRIB Analysis Group
+
+ Copyright (C) 1999-2013 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#ifndef __MELODIC_h
+#define __MELODIC_h
+
+#include <strstream>
+
+#ifdef __APPLE__
+#include <mach/mach.h>
+#define memmsg(msg) { \
+ MelodicOptions&opt = MelodicOptions::getInstance(); \
+ struct task_basic_info t_info; \
+ mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
+ if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
+ { \
+ if(opt.debug.value()) {\
+ cout << " MEM: " << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
+ cout.flush(); \
+ } \
+ } \
+}
+#else
+#define memmsg(msg) { \
+ MelodicOptions&opt = MelodicOptions::getInstance(); \
+ if(opt.debug.value()) {\
+ cout << msg; \
+ cout.flush(); \
+ } \
+}
+#endif
+
+
+
+// a simple message macro that takes care of cout and log
+#define message(msg) { \
+ MelodicOptions& opt = MelodicOptions::getInstance(); \
+ if(opt.verbose.value()) \
+ { \
+ cout << msg; \
+ } \
+ Log& logger = LogSingleton::getInstance(); \
+ logger.str() << msg; \
+ cout.flush(); \
+}
+
+#define dbgmsg(msg) { \
+ MelodicOptions&opt = MelodicOptions::getInstance(); \
+ if(opt.debug.value()) {\
+ cout << msg; } \
+}
+
+#define outMsize(msg,Mat) { \
+ MelodicOptions& opt = MelodicOptions::getInstance(); \
+ if(opt.debug.value()) \
+ cerr << " " << msg << " " <<Mat.Nrows() << " x " << Mat.Ncols() << endl; \
+}
+
+namespace Melodic{
+
+const string version = "3.14";
+
+// The two strings below specify the title and example usage that is
+// printed out as the help or usage message
+const string title=string("MELODIC (Version ")+version+")"+
+ string("\n Multivariate Exploratory Linear Optimised Decomposition into Independent Components\n")+
+ string("\nAuthor: Christian F. Beckmann \n Copyright(c) 2001-2013 University of Oxford");
+
+const string usageexmpl=string(" melodic -i <filename> <options>")+
+ string("\n \t \t to run melodic")+
+ // string("\n melodic -i <filename> --mix=melodic_mix")+
+ // string(" --filter=\"string of component numbers\"")+
+ // string("\n \t \t to remove estimated ICs from input")+
+ string("\n melodic -i <filename> --ICs=melodic_IC")+
+ string(" --mix=melodic_mix <options>")+
+ string("\n \t \t to run Mixture Model based inference on estimated ICs")+
+ string("\n melodic --help ");
+}
+
+#endif
diff --git a/meloptions.cc b/meloptions.cc
index b75604f5ad4d5b82b8f9a95461c0108b4f609d14..0cb87bc5990eee242a032558438df57a1385bc46 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -162,16 +162,20 @@ MelodicOptions* MelodicOptions::gopt = NULL;
//in the case of indirect inputs, create the vector of input names here
if(!FslFileExists(inputfname.value().at(0))){
std::vector< string > tmpfnames;
- ifstream fs(inputfname.value().at(0).c_str());
+ ifstream fs(inputfname.value().at(0).c_str() );
+ if (!fs.good()) {
+ cerr << "Error reading input file: " << inputfname.value().at(0).c_str() << endl;
+ exit(1);
+ }
string cline;
- while (!fs.eof()) {
+ while (fs.good()) {
getline(fs,cline);
if(cline.length()>0) {
while ( cline.find(' ') != cline.npos )
cline.erase( cline.find(' '),1);
tmpfnames.push_back(cline);
}
- }
+ }
fs.close();
inputfname.set_T(tmpfnames);
}
diff --git a/meloptions.h b/meloptions.h
new file mode 100644
index 0000000000000000000000000000000000000000..4b15d17bffb1f5362fe3eae7709fc08a321c15f0
--- /dev/null
+++ b/meloptions.h
@@ -0,0 +1,527 @@
+ /* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ meloptions.h - class for command line options
+
+ Christian F. Beckmann, FMRIB Analysis Group
+
+ Copyright (C) 1999-2013 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#ifndef __MELODICOPTIONS_h
+#define __MELODICOPTIONS_h
+
+#include <string>
+#include <strstream>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <stdlib.h>
+#include <stdio.h>
+#include "utils/options.h"
+#include "utils/log.h"
+#include "melodic.h"
+
+using namespace Utilities;
+
+namespace Melodic {
+
+class MelodicOptions {
+ public:
+ static MelodicOptions& getInstance();
+ ~MelodicOptions() { delete gopt; }
+
+ string version;
+ string binpath;
+ string logfname;
+ bool filtermode;
+ bool explicitnums;
+
+ Option<string> logdir;
+ Option< std::vector<string> > inputfname;
+
+ Option<string> outputfname;
+
+ Option<string> maskfname;
+ Option<bool> use_mask;
+ Option<bool> update_mask;
+ Option<bool> perf_bet;
+ Option<float> threshold;
+
+ Option<int> pca_dim;
+ Option<string> pca_est;
+ Option<bool> joined_whiten;
+ Option<bool> joined_vn;
+ Option<bool> dr_pca;
+ Option<bool> migp;
+ Option<int> migpN;
+ Option<bool> migp_shuffle;
+ Option<int> migp_factor;
+ Option<bool> dr;
+ Option<bool> dr_out;
+ Option<float> vn_level;
+ Option<int> numICs;
+ Option<string> approach;
+ Option<string> nonlinearity;
+
+ Option<bool> varnorm;
+ Option<bool> varnorm2;
+ Option<bool> pbsc;
+ Option<bool> pspec;
+ Option<string> segment;
+ Option<bool> tsmooth;
+ Option<float> epsilon;
+ Option<float> epsilonS;
+ Option<int> maxNumItt;
+ Option<int> maxRestart;
+ Option<int> rank1interval;
+
+ Option<string> mmthresh;
+ Option<bool> perf_mm;
+ Option<string> ICsfname;
+ Option<string> filtermix;
+ Option<string> smodename;
+ Option<string> filter;
+
+ Option<bool> genreport;
+ Option<string> guireport;
+ Option<string> bgimage;
+ Option<float> tr;
+ Option<bool> logPower;
+ Option<bool> addsigchng;
+ Option<bool> allPPCA;
+ Option<bool> varplots;
+ Option<bool> varvals;
+
+ Option<string> fn_Tdesign;
+ Option<string> fn_Tcon;
+ Option<string> fn_TconF;
+ Option<string> fn_Sdesign;
+ Option<string> fn_Scon;
+ Option<string> fn_SconF;
+
+ Option<bool> output_all;
+ Option<bool> output_unmix;
+ Option<bool> output_MMstats;
+ Option<bool> output_pca;
+ Option<bool> output_white;
+ Option<bool> output_origIC;
+ Option<bool> output_mean;
+
+ Option<bool> verbose;
+ Option<bool> vers;
+ Option<bool> copyright;
+ Option<bool> help;
+ Option<bool> debug;
+
+ Option<string> guessfname;
+ Option<string> paradigmfname;
+ Option<string> axials_str;
+
+ Option<int> dummy;
+ Option<int> repeats;
+ Option<int> seed;
+ Option<float> nlconst1;
+ Option<float> nlconst2;
+ Option<float> smooth_probmap;
+ Option<string> insta_fn;
+
+ Option<bool> remove_meanvol;
+ Option<bool> remove_meantc;
+ Option<bool> remove_endslices;
+ Option<bool> rescale_nht;
+
+ Option<bool> guess_remderiv;
+ Option<bool> temporal;
+
+ Option<float> retryfactor;
+ Option<int> econ;
+
+ void parse_command_line(int argc, char** argv, Log& logger, const string &p_version);
+
+ private:
+ MelodicOptions();
+ const MelodicOptions& operator=(MelodicOptions&);
+ MelodicOptions(MelodicOptions&);
+
+ OptionParser options;
+
+ static MelodicOptions* gopt;
+
+ void print_usage(int argc, char *argv[]);
+ void print_version();
+ void print_copyright();
+ void status();
+
+};
+
+ inline MelodicOptions& MelodicOptions::getInstance(){
+ if(gopt == NULL)
+ gopt = new MelodicOptions();
+
+ return *gopt;
+ }
+
+ inline MelodicOptions::MelodicOptions() :
+ logdir(string("-o,--outdir"), string("log.txt"),
+ 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),
+ outputfname(string("-O,--out"), string("melodic"),
+ string("output file name"),
+ false, requires_argument,false),
+ maskfname(string("-m,--mask"), string(""),
+ string("file name of mask for thresholding"),
+ false, requires_argument),
+ use_mask(string("--nomask"), true,
+ string("switch off masking"),
+ false, no_argument),
+ update_mask(string("--update_mask"), true,
+ string("switch off mask updating"),
+ false, no_argument),
+ perf_bet(string("--nobet"), true,
+ string("\tswitch off BET"),
+ false, no_argument),
+ threshold(string("--bgthreshold"), 0.01,
+ string("brain / non-brain threshold (only if --nobet selected)\n"),
+ false, requires_argument),
+ pca_dim(string("-d,--dim"), 0,
+ string("dimensionality reduction into #num dimensions (default: automatic estimation)"),
+ false, requires_argument),
+ pca_est(string("--dimest"), string("lap"),
+ string("use specific dim. estimation technique: lap, bic, mdl, aic, mean (default: lap)"),
+ false, requires_argument),
+ joined_whiten(string("--sep_whiten"), false,
+ string("switch on separate whitening"),
+ false, no_argument, false),
+ joined_vn(string("--sep_vn"), true,
+ string("switch on separate variance nomalisation (as opposed to separate VN)"),
+ false, no_argument),
+ dr_pca(string("--mod_pca"), true,
+ string("switch off modified PCA for concat ICA"),
+ false, no_argument, false),
+ migp(string("--disableMigp"), true,
+ string("switch off MIGP data reduction"),
+ false, no_argument),
+ migpN(string("--migpN"), 0,
+ string("Number of internal Eigenmaps"),
+ false, requires_argument),
+ migp_shuffle(string("--migp_shuffle"), true,
+ string("Randomise MIGP file order (default: TRUE)"),
+ false, no_argument),
+ migp_factor(string("--migp_factor"), 2,
+ string("Internal Factor of mem-threshold relative to number of Eigenmaps (default: 2)"),
+ false, requires_argument),
+ dr(string("--dr"), false,
+ string("Dual Regression (default: false)"),
+ false, no_argument, false),
+ dr_out(string("--dr_out"), false,
+ string("Dual Regression output for MIGP/concat ICA"),
+ false, no_argument, false),
+ 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),
+ approach(string("-a,--approach"), string("symm"),
+ string("approach for decomposition, 2D: defl, symm (default), 3D: tica, concat (default)"),
+ false, requires_argument),
+ nonlinearity(string("--nl"), string("pow3"),
+ string("\tnonlinearity: gauss, tanh, pow3, pow4"),
+ false, requires_argument),
+ varnorm(string("--vn,--varnorm"), true,
+ string("switch off variance normalisation"),
+ false, no_argument),
+ varnorm2(string("--vn2"), true,
+ string("switch off 2nd level variance normalisation"),
+ false, no_argument, false),
+ pbsc(string("--pbsc"), false,
+ string(" switch on conversion to percent BOLD signal change"),
+ false, no_argument, false),
+ pspec(string("--pspec"), false,
+ string(" switch on conversion to powerspectra"),
+ false, no_argument, false),
+ segment(string("--covarweight"), string(""),
+ string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"),
+ false, requires_argument),
+ tsmooth(string("--spca"), false,
+ string("smooth the eigenvectors prior to IC decomposition"),
+ false, no_argument, false),
+ epsilon(string("--eps"), 0.00005,
+ string("minimum error change"),
+ false, requires_argument),
+ epsilonS(string("--epsS"), 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),
+ maxRestart(string("--maxrestart"), -1,
+ string("maximum number of restarts\n"),
+ false, requires_argument),
+ rank1interval(string("--rank1interval"), 10,
+ string("number of iterations between rank-1 approximation (TICA)\n"),
+ false, requires_argument,false),
+ mmthresh(string("--mmthresh"), string("0.5"),
+ string("threshold for Mixture Model based inference"),
+ false, requires_argument),
+ perf_mm(string("--no_mm"), true,
+ string("\tswitch off mixture modelling on IC maps\n "),
+ false, no_argument),
+ ICsfname(string("--ICs"), string(""),
+ string("\tinput filename of the IC components file for mixture modelling"),
+ false, requires_argument),
+ filtermix(string("--mix"), string(""),
+ string("\tinput filename of mixing matrix for mixture modelling / filtering"),
+ false, requires_argument),
+ smodename(string("--smode"), string(""),
+ string("\tinput filename of matrix of session modes for report generation"),
+ false, requires_argument),
+ filter(string("-f,--filter"), string(""),
+ string("list of 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),
+ 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(""),
+ string(" design matrix across time-domain"),
+ 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),
+ 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),
+ output_all(string("--Oall"), false,
+ string(" output everything"),
+ false, no_argument),
+ output_unmix(string("--Ounmix"), false,
+ string("output unmixing matrix"),
+ false, no_argument),
+ output_MMstats(string("--Ostats"), false,
+ string("output thresholded maps and probability maps"),
+ false, no_argument),
+ output_pca(string("--Opca"), false,
+ string("\toutput PCA results"),
+ false, no_argument),
+ output_white(string("--Owhite"), false,
+ string("output whitening/dewhitening matrices"),
+ false, no_argument),
+ output_origIC(string("--Oorig"), false,
+ string("\toutput the original ICs"),
+ false, no_argument),
+ output_mean(string("--Omean"), false,
+ string("\toutput mean volume\n"),
+ false, no_argument),
+ verbose(string("-v,--verbose"), false,
+ string("switch on diagnostic messages"),
+ false, no_argument),
+ vers(string("-V,--version"), false,
+ string("prints version information"),
+ false, no_argument),
+ copyright(string("--copyright"), false,
+ string("prints copyright information"),
+ false, no_argument),
+ help(string("-h,--help"), false,
+ string("prints this help message"),
+ false, no_argument),
+ debug(string("--debug"), false,
+ string(" switch on debug messages"),
+ false, no_argument),
+ guessfname(string("--init_ica"), string(""),
+ string("file name of FEAT paradigm file (design.mat) for ICA initialisation"),
+ false, requires_argument, false),
+ paradigmfname(string("--init_pca"), string(""),
+ string("file name of FEAT paradigm file (design.mat) for PCA initialisation"),
+ false, requires_argument, false),
+ axials_str(string("--report_maps"), string(" -s 2 -A 950 "),
+ string("control string for spatial map images (see slicer)"),
+ false, requires_argument),
+ dummy(string("--dummy"), 0,
+ string("number of dummy volumes"),
+ false, requires_argument,false),
+ repeats(string("--repeats"), 1,
+ string("number of repeats (multistart)"),
+ false, requires_argument, false),
+ seed(string("--seed"), -1,
+ string("integer seed for random number generator within melodic"),
+ false, requires_argument, false),
+ nlconst1(string("--nl1,--nlconst1"), 1.0,
+ string("nonlinear constant 1"),
+ false, requires_argument, false),
+ nlconst2(string("--nl2,--nlconst2"), 1.0,
+ string("nonlinear constant 2"),
+ false, requires_argument, false),
+ smooth_probmap(string("--smooth_pm"), 0.0,
+ string("width of smoothing kernel for probability maps"),
+ false, requires_argument, false),
+ insta_fn(string("--insta_fn"), string(""),
+ string(" mask file name for instacorr calculation"),
+ false, requires_argument, false),
+ remove_meanvol(string("--keep_meanvol"), true,
+ string("do not subtract mean volume"),
+ 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("--remove_deriv"), false,
+ string("removes every second entry in paradigm file (EV derivatives)"),
+ false, no_argument, false),
+ temporal(string("--temporal"), false,
+ string("perform temporal ICA"),
+ false, no_argument, false),
+ retryfactor(string("--retryfactor"), float(0.95),
+ string("multiplicative factor for determining new dim if estimated dim fails to converge"),
+ false, requires_argument, false),
+ econ(string("--econ"), 20000,
+ string("set ctrl parameter for helperfns econ mode"),
+ false, requires_argument, false),
+ options(title, usageexmpl)
+ {
+ try {
+ options.add(logdir);
+ options.add(inputfname);
+ options.add(outputfname);
+ options.add(guessfname);
+ options.add(maskfname);
+ options.add(use_mask);
+ options.add(update_mask);
+ options.add(perf_bet);
+ options.add(threshold);
+ options.add(pca_dim);
+ options.add(pca_est);
+ options.add(joined_whiten);
+ options.add(joined_vn);
+ options.add(dr_pca);
+ options.add(migp);
+ options.add(migpN);
+ options.add(migp_shuffle);
+ options.add(migp_factor);
+ options.add(dr);
+ options.add(dr_out);
+ options.add(vn_level);
+ options.add(numICs);
+ options.add(approach);
+ options.add(nonlinearity);
+ options.add(varnorm);
+ options.add(varnorm2);
+ options.add(pbsc);
+ options.add(pspec);
+ options.add(segment);
+ options.add(tsmooth);
+ options.add(epsilon);
+ options.add(epsilonS);
+ options.add(maxNumItt);
+ options.add(maxRestart);
+ options.add(rank1interval);
+ options.add(mmthresh);
+ options.add(perf_mm);
+ options.add(ICsfname);
+ options.add(filtermix);
+ options.add(smodename);
+ options.add(filter);
+ options.add(genreport);
+ options.add(guireport);
+ 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(output_all);
+ options.add(output_unmix);
+ options.add(output_MMstats);
+ options.add(output_pca);
+ options.add(output_white);
+ options.add(output_origIC);
+ options.add(output_mean);
+ options.add(verbose);
+ options.add(vers);
+ options.add(copyright);
+ options.add(help);
+ options.add(debug);
+
+ options.add(guessfname);
+ options.add(paradigmfname);
+ options.add(axials_str);
+ options.add(dummy);
+ options.add(repeats);
+ options.add(seed);
+ options.add(nlconst1);
+ options.add(nlconst2);
+ options.add(smooth_probmap);
+ options.add(insta_fn);
+ options.add(remove_meanvol);
+ options.add(remove_meantc);
+ options.add(remove_endslices);
+ options.add(rescale_nht);
+ options.add(guess_remderiv);
+ options.add(temporal);
+ options.add(retryfactor);
+ options.add(econ);
+ }
+ catch(X_OptionError& e) {
+ options.usage();
+ cerr << endl << e.what() << endl;
+ }
+ catch(std::exception &e) {
+ cerr << e.what() << endl;
+ }
+
+ }
+}
+
+#endif
+
+
+
diff --git a/melreport.cc b/melreport.cc
new file mode 100644
index 0000000000000000000000000000000000000000..417f3ef78497e2d93826c2f13e58526a5f8e276e
--- /dev/null
+++ b/melreport.cc
@@ -0,0 +1,807 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melreport.cc - report generation
+
+ Christian F. Beckmann, FMRIB Analysis Group
+
+ Copyright (C) 1999-2013 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "newimage/newimageall.h"
+#include "utils/log.h"
+#include "melreport.h"
+#include "melhlprfns.h"
+#include "miscmaths/miscprob.h"
+
+namespace Melodic{
+
+ void MelodicReport::IC_rep(MelGMix &mmodel, int cnum, int dim, Matrix ICstats){
+
+ if( bool(opts.genreport.value()) ){
+ addlink(mmodel.get_prefix()+".html",num2str(cnum));
+ IChtml.setDir(report.getDir(),mmodel.get_prefix()+".html");
+
+ {//start IC page
+ IChtml << "<HTML><HEAD><link REL=stylesheet TYPE=text/css href=file:" +
+ (string) getenv("FSLDIR") +"/doc/fsl.css>" << endl
+ << "<style type=\"text/css\">OBJECT { width: 100% }</style>"
+ << "<TITLE>FSL</TITLE></HEAD>" << endl
+ << "<IFRAME height=" << int(melodat.get_numfiles()/30 + 1)*50
+ <<"px width=100% src=nav.html frameborder=0></IFRAME><BR>"<< endl
+ << "<Center>" << endl;
+
+ if(cnum>1)
+ IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
+ <<".html\"><</a> - ";
+ else
+ IChtml << "<a href=\"" << string("IC_")+num2str(melodat.get_mix().Ncols())
+ <<".html\"><</a> - ";
+
+ if(cnum<dim)
+ IChtml << "<a href=\"" << string("IC_")+num2str(cnum+1)
+ <<".html\">></a>";
+ else
+ IChtml << "<a href=\"" << string("IC_")+num2str(1)
+ <<".html\">></a>";
+ IChtml << "<p><H3>MELODIC Component " << num2str(cnum)
+ << "<br></b></H3><hr><p>" << endl;
+ }
+ {//output IC stats
+ if(ICstats.Storage()>0&&ICstats.Nrows()>=cnum){
+ IChtml << fixed << setprecision(2) << std::abs(ICstats(cnum,1)) << " % of explained variance";
+ if(ICstats.Ncols()>1)
+ IChtml << "; " << std::abs(ICstats(cnum,2)) << " % of total variance";
+ if(ICstats.Ncols()>2&&opts.addsigchng.value()){
+ IChtml << "<p>" <<endl;
+ IChtml << " " << ICstats(cnum,3) << " % signal change (pos peak voxel); " << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;
+ }
+ IChtml << "<hr><p>" << endl;
+ }
+ }
+
+ volume4D<float> tempVol;
+ if(mmodel.get_threshmaps().Storage()>0&&
+ (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols()))
+ {//Output thresholded IC map
+ tempVol.setmatrix(mmodel.get_threshmaps().Row(1),melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
+
+ volume<float> newvol;
+
+ miscpic newpic;
+
+ float map1min = std::max((map1 + binarise(tempVol[0],tempVol[0].min(),
+ float(0.0)) * map1.max()).min(),float(0.001));
+ float map1max = std::max(map1.max(),float(0.001));
+ float map2min = std::min(map2.min(),float(-0.001));
+ float map2max = std::min((map2 + binarise(tempVol[0],float(0.0),
+ tempVol[0].max()) * map2.min()).max(),float(-0.001));
+
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
+ melodat.get_bg().percentile(0.02),
+ melodat.get_bg().percentile(0.98),
+ map1min, map1max, map2max, map2min,
+ 0, 0);
+ char instr[10000];
+
+ sprintf(instr," ");
+ strcat(instr,axials_instr.c_str());
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+
+ "_thresh.png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - thresholded IC map ") + mmodel.get_infstr(1));
+ newpic.set_cbar(string("ysb"));
+
+ if((std::abs(map1.max()-map1.min())>0.01) || (std::abs(map2.max()-map2.min())>0.01))
+ newpic.slicer(newvol, instr);
+ else
+ newpic.slicer(newvol, instr);
+ IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
+ IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
+ +"_thresh.png\" ALT=\"MMfit\"></A><p>" << endl;
+ }
+
+ {//plot time course
+ IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
+ miscplot newplot;
+
+ Matrix tmptc = melodat.get_Tmodes(cnum-1).Column(1).t();
+
+ newplot.col_replace(0,0xFF0000);
+
+ newplot.add_label(string("IC ")+num2str(cnum)+" time course");
+ //add GLM OLS fit
+ if(melodat.Tdes.Storage()>0 &&
+ melodat.glmT.get_beta().Nrows() == melodat.Tdes.Ncols()){
+ tmptc &= melodat.glmT.get_beta().Column(cnum).t() * melodat.Tdes.t();
+ newplot.add_label("full model fit");
+ }
+
+ //add deviation around time course
+ if(melodat.get_Tmodes(cnum-1).Ncols()>1 && opts.varplots.value()){
+ Matrix tmp = stdev(melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t(),1);
+ tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()+tmp;
+ tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()-tmp;
+ newplot.add_label("std error across subjects");
+ newplot.col_replace(tmptc.Nrows()-1,0x808080);
+ newplot.col_replace(tmptc.Nrows()-2,0x808080);
+ }
+
+ if(opts.tr.value()>0.0)
+ newplot.add_xlabel(string("Time (seconds); TR = ")+
+ float2str(opts.tr.value(),0,2,0)+" s");
+ else
+ newplot.add_xlabel(string("Time (TRs)"));
+
+ newplot.add_ylabel("Normalised Response");
+ newplot.set_yrange(tmptc.Row(1).Minimum()-0.05*(tmptc.Row(1).Maximum() -
+ tmptc.Row(1).Minimum()),tmptc.Row(1).Maximum()+
+ 0.05*(tmptc.Row(1).Maximum()-tmptc.Row(1).Minimum()));
+ newplot.grid_swapdefault();
+
+ newplot.set_xysize(750,150);
+
+ newplot.timeseries(tmptc,
+ report.appendDir(string("t")+num2str(cnum)+".png"),
+ string("Timecourse No. ")+num2str(cnum),
+ opts.tr.value(),150,12,1,false);
+
+ if(melodat.get_Tmodes(cnum-1).Ncols()>1)
+ tmptc &= melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t();
+ write_ascii_matrix(report.appendDir(string("t")
+ +num2str(cnum)+".txt"),tmptc.t());
+ IChtml << "<A HREF=\"" << string("t")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
+
+ if(melodat.get_numfiles()>1 && melodat.explained_var.Storage()>0
+ && melodat.explained_var.Ncols()>=cnum && opts.varvals.value())
+ IChtml << "Rank-1 approximation of individual time courses explains "
+ << std::abs(melodat.explained_var(cnum)) << "% of variance.<p>" << endl;
+ }//time series plot
+
+ if(!opts.pspec.value())
+ {//plot frequency
+ miscplot newplot;
+ RowVector empty(1);
+ empty = 0.0;
+ int fact = int(std::pow(10.0,int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
+ if(opts.logPower.value())
+ newplot.add_ylabel(string("log-Power"));
+ else
+ newplot.add_ylabel(string("Power"));
+
+ Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1).Column(1), opts.logPower.value());
+
+ newplot.set_Ylabel_fmt("%.0f");
+ newplot.set_yrange(0.0,1.02*fmixtc.Maximum());
+ newplot.grid_swapdefault();
+ newplot.set_xysize(750,150);
+
+ if(opts.tr.value()>0.0){
+ newplot.add_xlabel(string("Frequency (in Hz / ")+num2str(fact)+ " )");
+ newplot.timeseries(empty | fmixtc.t(),
+ report.appendDir(string("f")+
+ num2str(cnum)+".png"),
+ string("Powerspectrum of timecourse"),
+ fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()), 150,0,2);
+ }else{
+ newplot.add_xlabel(string("Frequency (in cycles); ")
+ +"frequency(Hz)=cycles/("
+ +num2str(melodat.get_Tmodes(0).Nrows())
+ +"* TR); period(s)=("
+ +num2str(melodat.get_Tmodes(0).Nrows())
+ +"* TR)/cycles"
+ );
+ newplot.timeseries(fmixtc.t(),
+ report.appendDir(string("f")+num2str(cnum)+".png"),
+ string("Powerspectrum of timecourse"));
+ }
+ write_ascii_matrix(report.appendDir(string("f")
+ +num2str(cnum)+".txt"), fmixtc);
+ IChtml << "<A HREF=\"" << string("f")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
+ }//frequency plot
+ {//add T-mode GLM F-stats for full model fit & contrasts
+ if(melodat.Tdes.Storage() > 0 &&
+ melodat.glmT.get_beta().Nrows() == melodat.Tdes.Ncols()){
+ IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
+ "<CAPTION><EM> <b>GLM (OLS) on time series </b></EM></CAPTION>" << endl
+ << "<TR valign=middle><TH ><EM>GLM β's</EM> <TH> <EM> F-test on <br> full model fit </em>";
+ if(melodat.Tcon.Storage() > 0)
+ IChtml << "<TH ><EM>Contrasts</EM>"<<endl;
+ IChtml << "<TR><TD><TABLE border=0><TR><TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Tdes.Ncols();ctr++)
+ IChtml << " PE(" <<num2str(ctr)+"): <br>" << endl;
+ IChtml << "<TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Tdes.Ncols();ctr++)
+ IChtml << melodat.glmT.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
+ IChtml << "</TABLE>" <<
+ " <TD align=center> F = "<< melodat.glmT.get_f_fmf().Column(cnum) <<
+ " <BR> dof1 = " << melodat.Tdes.Ncols() << "; dof2 = "
+ << melodat.glmT.get_dof() << "<BR>" <<endl;
+ if(melodat.glmT.get_pf_fmf().Column(cnum).AsScalar() < 0.05)
+ IChtml << fixed << setprecision(5) <<"<b> p < " << melodat.glmT.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)<b></TD>" << endl;
+ else
+ IChtml << fixed << setprecision(5) << " p < " <<
+ melodat.glmT.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)</TD>" << endl;
+ if(melodat.Tcon.Storage() > 0 &&
+ melodat.Tdes.Ncols() == melodat.Tcon.Ncols()){
+ IChtml << fixed << setprecision(2) << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ IChtml << "COPE(" << num2str(ctr) << "): <br>" << endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ IChtml <<" z = <BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ IChtml << melodat.glmT.get_z().Column(cnum).Row(ctr) <<";<BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ if(melodat.glmT.get_p().Column(cnum).Row(ctr).AsScalar() < 0.05)
+ IChtml << fixed << setprecision(5) << "<b> p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
+ "</b><BR>" << endl;
+ else
+ IChtml << fixed << setprecision(5) <<" p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
+ "<BR>" << endl;
+ IChtml << "</TABLE></td></tr>" << endl;
+ }
+ }
+ IChtml << "</TABLE><p>" << endl;
+ }
+
+ if(cnum <= (int)melodat.get_Smodes().size())
+ {//plot subject mode
+
+ Matrix smode;
+ smode = melodat.get_Smodes(cnum-1);
+
+ if(smode.Nrows() > 1){
+ IChtml << "<hr><H3> Sessions/Subjects mode </H3><p>" << endl <<endl;
+ miscplot newplot;
+
+ //add GLM OLS fit
+ //newplot.setscatter(smode,2);
+
+ if(melodat.Sdes.Storage() > 0&&
+ melodat.glmS.get_beta().Nrows() == melodat.Sdes.Ncols()){
+ smode |= melodat.Sdes * melodat.glmS.get_beta().Column(cnum);
+ newplot.add_label(string("IC ")+num2str(cnum)+" subject/session-mode");
+ newplot.add_label("full model fit");
+ }
+ newplot.grid_swapdefault();
+ newplot.set_Ylabel_fmt("%.2f");
+ newplot.add_xlabel(" Subject Number");
+
+
+ newplot.set_xysize(750,150);
+ newplot.timeseries(smode.t(),
+ report.appendDir(string("s")+num2str(cnum)+".png"),
+ string("Subject/Session mode No. ") + num2str(cnum));
+ newplot.clear_xlabel();
+ newplot.clear_labels();
+ newplot.set_xysize(120,200);
+ newplot.set_minmaxscale(1.1);
+
+ newplot.boxplot((Matrix)smode.Column(1),
+ report.appendDir(string("b")+num2str(cnum)+".png"),
+ string("Subject/Session mode"));
+ write_ascii_matrix(report.appendDir(string("s")
+ +num2str(cnum)+".txt"), smode);
+ IChtml << "<A HREF=\"" << string("s")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("s")+num2str(cnum)+".png\"></A>" << endl;
+ IChtml << "<A HREF=\"" << string("s")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("b")+num2str(cnum)+".png\"></A><p>" << endl;
+ }
+ {//add S-mode GLM F-stats for full model fit & contrasts
+ if(melodat.Sdes.Storage() > 0 &&
+ melodat.glmS.get_beta().Nrows() == melodat.Sdes.Ncols()){
+ IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
+ "<CAPTION><EM> <b>GLM (OLS) on subject/session-mode </b></EM></CAPTION>" << endl
+ << "<TR valign=middle><TH ><EM>GLM β's</EM> <TH> <EM> F-test on <br> full model fit </em>";
+ if(melodat.Scon.Storage() > 0)
+ IChtml << "<TH ><EM>Contrasts</EM>"<<endl;
+ IChtml << "<TR><TD><TABLE border=0><TR><TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
+ IChtml << " PE(" <<num2str(ctr)+"): <br>" << endl;
+ IChtml << "<TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
+ IChtml << melodat.glmS.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
+ IChtml << "</TABLE>" <<
+ " <TD align=center> F = "<< melodat.glmS.get_f_fmf().Column(cnum) <<
+ " <BR> dof1 = " << melodat.Sdes.Ncols() << "; dof2 = "
+ << melodat.glmS.get_dof() << "<BR>" <<endl;
+ if(melodat.glmS.get_pf_fmf().Column(cnum).AsScalar() < 0.05)
+ IChtml << fixed << setprecision(5) <<"<b> p < " << melodat.glmS.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)<b></TD>" << endl;
+ else
+ IChtml << fixed << setprecision(5) << " p < " <<
+ melodat.glmS.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)</TD>" << endl;
+ if(melodat.Scon.Storage() > 0 && melodat.Sdes.Storage() > 0 &&
+ melodat.Sdes.Ncols() == melodat.Scon.Ncols()){
+ IChtml << fixed << setprecision(2) << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml << "COPE(" << num2str(ctr) << "): <br>" << endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml <<" z = <BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml << melodat.glmS.get_z().Column(cnum).Row(ctr) <<";<BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ if(melodat.glmS.get_p().Column(cnum).Row(ctr).AsScalar() < 0.05)
+ IChtml << fixed << setprecision(5) << "<b> p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
+ "</b><BR>" << endl;
+ else
+ IChtml << fixed << setprecision(5) <<" p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
+ "<BR>" << endl;
+ IChtml << "</TABLE></td></tr>" << endl;
+ }
+ }
+ IChtml << "</TABLE><p>" << endl;
+ }
+ }//subject mode plot
+
+ if(mmodel.get_threshmaps().Storage()>0&&
+ (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
+ (mmodel.get_threshmaps().Nrows()>1))
+ {//Output other thresholded IC map
+
+ for(int tctr=2; tctr<=mmodel.get_threshmaps().Nrows(); tctr++){
+ tempVol.setmatrix(mmodel.get_threshmaps().Row(tctr),melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
+
+ volume<float> newvol;
+ miscpic newpic;
+
+ float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
+ float(0.0)) * map1.max()).min();
+ float map1max = map1.max();
+ float map2min = map2.min();
+ float map2max = (map2 + binarise(tempVol[0],float(0.0),
+ tempVol[0].max()) * map2.min()).max();
+
+ //cerr << endl << map1min << " " << map1max << endl
+ // << map2min << " " << map2max << endl;
+
+ // if(map1.max()-map1.min()>0.01)
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
+ melodat.get_bg().percentile(0.02),
+ melodat.get_bg().percentile(0.98),
+ map1min, map1max, map2max, map2min,
+ 0, 0);
+
+ char instr[10000];
+
+ sprintf(instr," ");
+ strcat(instr,axials_instr.c_str());
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+"_thresh"+
+ num2str(tctr)+".png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - thresholded IC map ("+
+ num2str(tctr)+") ")+
+ mmodel.get_infstr(tctr));
+ newpic.set_cbar(string("ysb"));
+ //cerr << instr << endl;
+ newpic.slicer(newvol, instr);
+
+ IC_rep_det(mmodel, cnum, dim);
+ IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
+ IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
+ +"_thresh"+num2str(tctr)+".png\" ALT=\"MMfit\"></A><p>" << endl;
+ }
+ }
+
+ { //finish IC page
+ IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
+ << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "FMRIB Software Library</A>.</FONT></Center>" << endl
+ << "</BODY></HTML>" << endl;
+ } //finish IC page
+ IC_rep_det(mmodel, cnum, dim);
+ }
+ }
+
+ void MelodicReport::IC_rep_det(MelGMix &mmodel, int cnum, int dim){
+ if( bool(opts.genreport.value()) ){
+
+ {//start IC2 page
+ IChtml2.setDir(report.getDir(),mmodel.get_prefix()+"_MM.html");
+ IChtml2 << "<HTML><HEAD><link REL=stylesheet TYPE=text/css href=file:" +
+ (string) getenv("FSLDIR") +"/doc/fsl.css>" << endl
+ << "<style type=\"text/css\">OBJECT { width: 100% }</style>"
+ << "<TITLE>FSL</TITLE></HEAD>" << endl
+ << "<IFRAME height="<< int(melodat.get_numfiles()/30 + 1)*50
+ <<"px width=100% src=nav.html frameborder=0></IFRAME><p>"<< endl
+ << "<CENTER>";
+ if(cnum>1)
+ IChtml2 << "<b><a href=\"" << string("IC_")+num2str(cnum-1)
+ <<"_MM.html\"><</a> - ";
+ else
+ IChtml2 << "<b><a href=\"" << string("IC_")+num2str(melodat.get_mix().Ncols())
+ <<"_MM.html\"><</a> - ";
+ // IChtml << "<a href=\"00index.html\"> index </a>" ;
+
+ if(cnum<dim)
+ IChtml2 << "<a href=\"" << string("IC_")+num2str(cnum+1)
+ <<"_MM.html\">></a>";
+ else
+ IChtml2 << "<a href=\"" << string("IC_")+num2str(1)
+ <<"_MM.html\">></a>";
+ IChtml2 << "<p><H3>Component " << num2str(cnum)
+ << " Mixture Model fit <br></b></H3><hr><p>" << endl;
+ }
+
+ volume4D<float> tempVol;
+
+ if(melodat.get_IC().Storage()>0)
+ {//Output raw IC map
+
+ // tempVol.setmatrix(melodat.get_IC().Row(cnum),
+ // melodat.get_mask());
+ tempVol.setmatrix(mmodel.get_data(),
+ melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0),
+ tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(),
+ float(-0.0));
+
+ volume<float> newvol;
+ miscpic newpic;
+
+ // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
+ // float(0.0)) * map1.max()).robustmin();
+ float map1max = map1.percentile(0.99);
+ float map2min = map2.percentile(0.01);
+ //float map2max = (map2 + binarise(tempVol[0],float(0.0),
+ // tempVol[0].max()) * map2.min()).robustmax();
+
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
+ float(0.0),
+ float(0.0),
+ float(0.01), map1max, float(-0.01), map2min,
+ 0, 0);
+
+ char instr[10000];
+
+ sprintf(instr," ");
+ strcat(instr,axials_instr.c_str());
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+
+ ".png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - raw Z transformed IC map (1 - 99 percentile)"));
+ newpic.set_cbar(string("ysb"));
+
+ newpic.slicer(newvol, instr);
+ }
+ IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
+ IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
+ ".png\"><A><p>" << endl;
+
+ if(mmodel.get_probmap().Storage()>0&&
+ (mmodel.get_probmap().Ncols() == mmodel.get_data().Ncols())&&
+ (mmodel.get_probmap().Nrows() == mmodel.get_data().Nrows()))
+ {//Output probmap
+ tempVol.setmatrix(mmodel.get_probmap(),melodat.get_mask());
+
+ volume<float> map;
+ map = tempVol[0];
+
+ volume<float> newvol;
+ miscpic newpic;
+
+ newpic.overlay(newvol, melodat.get_bg(), map, map,
+ melodat.get_bg().percentile(0.02),
+ melodat.get_bg().percentile(0.98),
+ float(0.1), float(1.0), float(0.0), float(0.0),
+ 0, 0);
+
+ char instr[10000];
+
+ sprintf(instr," ");
+ strcat(instr,"-l render1 ");
+ strcat(instr,axials_instr.c_str());
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+
+ "_prob.png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - Mixture Model probability map"));
+
+ newpic.set_cbar(string("y"));
+ newpic.slicer(newvol, instr);
+
+ IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
+ IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
+ "_prob.png\">" << endl;
+
+ IChtml2 << "</A><p>" << endl;
+ }
+
+ RowVector dat = mmodel.get_data().Row(1);
+ if(dat.Maximum()>dat.Minimum())
+ {//Output GGM/GMM fit
+ miscplot newplot;
+
+ if(mmodel.get_type()=="GGM"){
+ newplot.add_label("IC map histogram");
+ newplot.add_label("full GGM fit");
+ newplot.add_label("background Gaussian");
+ newplot.add_label("Gamma distributions");
+ newplot.gmmfit(mmodel.get_data().Row(1),
+ mmodel.get_means(),
+ mmodel.get_vars(),
+ mmodel.get_pi(),
+ report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
+ string(mmodel.get_prefix() +
+ " Gaussian/Gamma Mixture Model("+num2str(mmodel.mixtures())+") fit"),
+ true, float(0.0), float(0.0));
+
+ }
+ else{
+ newplot.add_label("IC map histogram");
+ newplot.add_label("full GMM fit");
+ newplot.add_label("individual Gaussians");
+ newplot.gmmfit(mmodel.get_data().Row(1),
+ mmodel.get_means(),
+ mmodel.get_vars(),
+ mmodel.get_pi(),
+ report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
+ string(mmodel.get_prefix() +
+ " Gaussian Mixture Model("+num2str(mmodel.mixtures())+") fit"),
+ false, float(0.0), float(2.0));
+
+ }
+
+ // IChtml2 << "<A HREF=\"" +mmodel.get_prefix()+"_MMfit_detail.png\"> ";
+ IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
+ "_MMfit.png\"><p>" << endl;
+ } //GGM/GMM plot
+
+ {//MM parameters
+ IChtml2 << "<br> " << mmodel.get_prefix()
+ << " Mixture Model fit <br>" << endl
+ << "<br> Means : " << mmodel.get_means() << endl
+ << "<br> Vars : " << mmodel.get_vars() << endl
+ << "<br> Prop. : " << mmodel.get_pi() << endl;
+ }
+
+ { //finish IC2 page
+ IChtml2<< "<HR><FONT SIZE=1>This page produced automatically by "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
+ << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "FMRIB Software Library</A>.</FONT></CENTER>" << endl
+ << "</BODY></HTML>" << endl;
+ } //finish IC2 page
+ }
+ }
+
+ void MelodicReport::IC_simplerep(string prefix, int cnum, int dim){
+ if( bool(opts.genreport.value()) ){
+ addlink(prefix+".html",num2str(cnum));
+ IChtml.setDir(report.getDir(),prefix+".html");
+
+ {//start IC page
+
+ IChtml << "<HTML> " << endl
+ << "<TITLE>MELODIC Component " << num2str(cnum)
+ << "</TITLE>" << endl
+ << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
+ << "/doc/images/fsl-bg.jpg\">" << endl
+ << "<hr><CENTER><H1>MELODIC Component " << num2str(cnum)
+ << "</H1>"<< endl;
+
+ if(cnum>1)
+ IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
+ <<".html\">previous</a> - ";
+
+ IChtml << "<a href=\"00index.html\"> index </a>" ;
+
+ if(cnum<dim)
+ IChtml << " - <a href=\"" << string("IC_")+num2str(cnum+1)
+ <<".html\">next</a><p>";
+
+ IChtml << "<hr><p>" << endl;
+ }
+
+ volume4D<float> tempVol;
+
+ if(melodat.get_IC().Storage()>0)
+ {//Output raw IC map
+
+ tempVol.setmatrix(melodat.get_IC().Row(cnum),
+ melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0),
+ tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(),
+ float(-0.0));
+
+ volume<float> newvol;
+ miscpic newpic;
+
+ // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
+ // float(0.0)) * map1.max()).robustmin();
+ float map1max = map1.percentile(0.99);
+ float map2min = map2.percentile(0.01);
+ //float map2max = (map2 + binarise(tempVol[0],float(0.0),
+ // tempVol[0].max()) * map2.min()).robustmax();
+
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
+ float(0.0),
+ float(0.0),
+ float(0.01), map1max, float(-0.01), map2min,
+ 0, 0);
+
+ char instr[10000];
+
+ sprintf(instr," ");
+ strcat(instr,axials_instr.c_str());
+ strcat(instr,string(report.appendDir(prefix+
+ ".png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - raw Z transformed IC map (1 - 99 percentile)"));
+ newpic.set_cbar(string("ysb"));
+
+ newpic.slicer(newvol, instr);
+ }
+
+ IChtml << "<img BORDER=0 SRC=\""+ prefix+
+ ".png\"><p>" << endl;
+
+ {//plot time course
+ miscplot newplot;
+
+ if(opts.tr.value()>0.0)
+ newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
+ report.appendDir(string("t")+
+ num2str(cnum)+".png"),
+ string("Timecourse (in seconds); TR = ")+
+ float2str(opts.tr.value(),0,2,0)+" s",
+ opts.tr.value(),150,4,1);
+ else
+ newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
+ report.appendDir(string("t")+
+ num2str(cnum)+".png"),
+ string("Timecourse (in TRs)"));
+ write_ascii_matrix(report.appendDir(string("t")
+ +num2str(cnum)+".txt"),
+ melodat.get_Tmodes(cnum-1));
+ IChtml << "<A HREF=\"" << string("t")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
+ }//time series plot
+
+ {//plot frequency
+ miscplot newplot;
+ int fact = int(std::pow(10.0,
+ int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
+
+ if(opts.tr.value()>0.0)
+ newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
+ report.appendDir(string("f")+
+ num2str(cnum)+".png"),
+ string("FFT of timecourse (in Hz / ") +
+ num2str(fact)+")",
+ fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()),
+ 150,0,2);
+ else
+ newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
+ report.appendDir(string("f")+
+ num2str(cnum)+".png"),
+ string(string("FFT of timecourse (in cycles); ")
+ +"frequency(Hz)=cycles/("
+ +num2str(melodat.get_Tmodes(0).Nrows())
+ +"* TR); period(s)=("
+ +num2str(melodat.get_Tmodes(0).Nrows())
+ +"* TR)/cycles"));
+ write_ascii_matrix(report.appendDir(string("f")
+ +num2str(cnum)+".txt"),
+ melodat.get_Tmodes(cnum-1));
+ IChtml << "<A HREF=\"" << string("f")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
+ }//frequency plot
+
+ { //finish IC page
+ IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
+ << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "FMRIB Software Library</A>.</FONT>" << endl
+ << "</BODY></HTML>" << endl;
+ } //finish IC page
+ }
+ }
+
+ void MelodicReport::PPCA_rep(){
+
+ {//plot time course
+ report << "<p><hr><b>PCA estimates </b> <p>" << endl;
+
+ Matrix what;
+ miscplot newplot;
+
+ what = melodat.get_EV();
+ what &= melodat.get_EVP();
+
+ newplot.add_label("ordered Eigenvalues");
+ newplot.add_label("% of expl. variance");
+
+ if(melodat.get_PPCA().Storage()>0){
+ what = what.Columns(1,melodat.get_PPCA().Nrows());
+ if(opts.allPPCA.value()&&melodat.get_PPCA().Ncols()==7){
+ what &= melodat.get_PPCA().Columns(3,7).t();
+ newplot.add_label("Laplace");
+ newplot.add_label("BIC");
+ newplot.add_label("MDL");
+ newplot.add_label("RRN");
+ newplot.add_label("AIC");
+ }else{
+ what &= melodat.get_PPCA().Column(1).t();
+ newplot.add_label("dim. estimate");
+ }
+ }
+
+ newplot.set_Ylabel_fmt("%.2f");
+ newplot.add_xlabel("Number of included components");
+ newplot.set_yrange(0.0,1.02);
+ newplot.grid_swapdefault();
+ newplot.timeseries(what,
+ report.appendDir("EVplot.png"),
+ string("Eigenspectrum Analysis"),
+ 0,450,4,0);
+
+ report << "<img BORDER=0 SRC=\"EVplot.png\"><p>" << endl;
+ }//time series plot
+ }
+
+ void MelodicReport::Smode_rep(){
+ if(melodat.get_Smodes().size()>0){
+ report << "<p><hr><b>TICA Subject/Session modes </b> <br>" << endl;
+ miscplot newplot;
+ report << "Boxplots show the relative response amplitudes across the "
+ << "session/subject domain (" << melodat.get_numfiles()
+ << " input files). Components have been sorted in decreasing order of "
+ << " the median response per component. <br><br>";
+
+ outMsize("Smode.at(0)", melodat.get_Smodes().at(0));
+ Matrix allmodes = melodat.get_Smodes().at(0);
+ for(int ctr = 1; ctr < (int)melodat.get_Smodes().size();++ctr)
+ allmodes |= melodat.get_Smodes().at(ctr);
+
+ outMsize("allmodes", allmodes);
+ newplot.add_xlabel("Component No.");
+ newplot.add_ylabel("");
+ if(allmodes.Ncols()<100)
+ newplot.set_xysize(100+30*allmodes.Ncols(),300);
+ else
+ newplot.set_xysize(1200,300);
+ newplot.boxplot(allmodes,report.appendDir(string("bp_Smodes.png")),
+ string("Subject/Session modes"));
+ report << "<img BORDER=0 SRC=\"bp_Smodes.png\"><p>" << endl;
+ }
+ }
+}