diff --git a/meldata.cc b/meldata.cc
index 1cb150b20e9f13a09c6d650a4b4ebf3b577b6145..f7457755f4424cce50342c16f1ebe6bacc028b9a 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -63,7 +63,7 @@ namespace Melodic{
}
meanC = mean(tmpData,2);
- tmpData = remmean(tmpData,2);
+ // tmpData = remmean(tmpData,2);
//convert to power spectra
if(opts.pspec.value()){
@@ -86,7 +86,7 @@ namespace Melodic{
//variance - normalisation
if(opts.varnorm.value()){
message(" Normalising by voxel-wise variance ...");
- stdDev = varnorm(tmpData,tmpData.Nrows(),2.3);
+ stdDev = varnorm(tmpData,std::min(30,tmpData.Nrows()-1),2.3);
stdDevi = pow(stdDev,-1);
message(" done" << endl);
}
@@ -131,9 +131,10 @@ namespace Melodic{
void MelodicData::set_TSmode()
{
- message("Calculating T- and S-modes " << endl);
Matrix tmp, tmpT, tmpS, tmpT2, tmpS2;
tmp = expand_dimred(mixMatrix);
+ saveascii(tmp,"exp");
+ saveascii(mixMatrix,"mix");
tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
tmpS = zeros(numfiles, tmp.Ncols());
krfact(tmp,tmpT,tmpS);
@@ -175,11 +176,15 @@ namespace Melodic{
// bool tmpvarnorm = opts.varnorm.value();
// opts.varnorm.set_T(false);
alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
- for(int ctr = 1; ctr < numfiles; ctr++){
- tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
+ if(opts.debug.value())
+ save4D(alldat,string("preproc_dat") + num2str(1));
+ for(int ctr = 1; ctr < numfiles; ctr++){
+ tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
+ if(opts.debug.value())
+ save4D(tmpData /numfiles,string("preproc_dat") + num2str(ctr+1));
alldat += tmpData / numfiles;
}
-
+
// opts.varnorm.set_T(tmpvarnorm);
//variance - normalisation
// if(opts.varnorm.value()){
@@ -199,25 +204,37 @@ namespace Melodic{
if(numfiles>1)
message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);
+ if(opts.debug.value())
+ save4D(alldat,"alldat");
//estimate model order
ColumnVector tmpPPCA;
RowVector AdjEV, PercEV;
Matrix Corr, tmpE;
int order;
- order = ppca_dim(alldat, RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
+ order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
- if(opts.pca_dim.value() == 0){
+ if(opts.pca_dim.value() == 0){
opts.pca_dim.set_T(order);
PPCA=tmpPPCA;
}
order = opts.pca_dim.value();
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");
+ outMsize("PercEV",PercEV); saveascii(PercEV,"PercEV");
+ outMsize("tmpPPCA",tmpPPCA); saveascii(tmpPPCA,"tmpPPCA");
+ outMsize("whiteMatrix",whiteMatrix); saveascii(whiteMatrix,"whiteMatrix");
+ outMsize("dewhiteMatrix",dewhiteMatrix); saveascii(dewhiteMatrix,"dewhiteMatrix");
+ cerr << "Order: " << order << endl;
+ }
+
EV = AdjEV;
EVP = PercEV;
-
if(numfiles == 1){
Data = alldat;
Matrix tmp = Identity(Data.Nrows());
@@ -230,15 +247,21 @@ namespace Melodic{
// if(opts.varnorm.value())
// varnorm(tmpData,stdDev);
- // whiten (separate / joint)
+ // whiten (separate / joint)
+ Matrix newWM,newDWM;
if(!opts.joined_whiten.value()){
message(" Individual whitening in a " << order << " dimensional subspace " << endl);
std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
- calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
+ calc_white(pcaE, pcaD, order, newWM, newDWM);
+ }else{
+ std_pca(whiteMatrix*tmpData, RXweight, Corr, pcaE, pcaD);
+ calc_white(pcaE, pcaD, order, newWM, newDWM);
+ newDWM=(dewhiteMatrix*newDWM);
+ newWM=(newWM*whiteMatrix);
}
- tmpData = whiteMatrix * tmpData;
- DWM.push_back(dewhiteMatrix);
- WM.push_back(whiteMatrix);
+ DWM.push_back(newDWM);
+ WM.push_back(newWM);
+ tmpData = newWM * tmpData;
//concatenate Data
if(Data.Storage() == 0)
@@ -251,7 +274,8 @@ namespace Melodic{
message(" Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl);
meanC=mean(Data,2);
-
+ if(opts.debug.value())
+ save4D(Data,"concat_data");
//save the mean & mask
save_volume(Mask,logger.appendDir("mask"));
save_volume(Mean,logger.appendDir("mean"));
@@ -295,17 +319,17 @@ namespace Melodic{
//read in post-proc design matrices etc
if(opts.fn_Tdesign.value().length()>0)
- Tdes = read_vest(opts.fn_Tdesign.value());
+ Tdes = read_ascii_matrix(opts.fn_Tdesign.value());
if(opts.fn_Sdesign.value().length()>0)
- Sdes = read_vest(opts.fn_Sdesign.value());
+ Sdes = read_ascii_matrix(opts.fn_Sdesign.value());
if(opts.fn_Tcon.value().length()>0)
- Tcon = read_vest(opts.fn_Tcon.value());
+ Tcon = read_ascii_matrix(opts.fn_Tcon.value());
if(opts.fn_Scon.value().length()>0)
- Scon = read_vest(opts.fn_Scon.value());
+ Scon = read_ascii_matrix(opts.fn_Scon.value());
if(opts.fn_TconF.value().length()>0)
- TconF = read_vest(opts.fn_TconF.value());
+ TconF = read_ascii_matrix(opts.fn_TconF.value());
if(opts.fn_SconF.value().length()>0)
- SconF = read_vest(opts.fn_SconF.value());
+ SconF = read_ascii_matrix(opts.fn_SconF.value());
Tdes = remmean(Tdes,1);
Sdes = remmean(Sdes,1);
diff --git a/melhlprfns.cc b/melhlprfns.cc
index eaa1d0eee3968addd7a05b5b7bd09c817427c471..78780ea30dc3572a11818150d5e465d90b161d84 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -132,7 +132,7 @@ namespace Melodic{
Matrix tmpE, white, dewhite;
RowVector tmpD, tmpD2;
- std_pca(in, Corr, tmpE, tmpD);
+ std_pca(remmean(in,2), Corr, tmpE, tmpD);
calc_white(tmpE,tmpD, dim, white, dewhite);
Matrix ws = white * in;
diff --git a/melica.cc b/melica.cc
index 1379340ef85398412b09fde334bc47847316f45c..17a3c958a2da252ca37af5c07413a657303a4dda 100644
--- a/melica.cc
+++ b/melica.cc
@@ -32,7 +32,7 @@ namespace Melodic{
// see www.cis.hut.fi/projects/ica/fastica/
//initialize matrices
- Matrix redUMM_old;
+ Matrix redUMM_old, redUMM_veryold;
Matrix tmpU;
//srand((unsigned int)timer(NULL));
@@ -50,12 +50,13 @@ namespace Melodic{
symm_orth(redUMM);
- int itt_ctr,itt_ctr2=1,cum_itt=0;
+ int itt_ctr,itt_ctr2=1,cum_itt=0,newmaxitts = opts.maxNumItt.value();
double minAbsSin = 1.0;
if(opts.approach.value() == string("tica"))
opts.maxNumItt.set_T(opts.rank1interval.value());
do{// TICA loop
itt_ctr = 1;
+ redUMM_veryold = redUMM;
do{ // da loop!!!
redUMM_old = redUMM;
//calculate IC estimates
@@ -92,6 +93,9 @@ namespace Melodic{
// orthogonalize the unmixing-matrix
symm_orth(redUMM);
+ if(opts.approach.value() == string("tica")){
+ message(" TICA:");
+ }
//termination condition : angle between old and new < epsilon
minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_old)).Minimum();
message(" Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl);
@@ -113,9 +117,12 @@ namespace Melodic{
outMsize(" 2nd unmmixing matrix ", temp);
redUMM = melodat.get_white() * temp;
outMsize(" redUMM ", redUMM);
+ minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_veryold)).Minimum();
+ message(" Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl);
+ if(abs(minAbsSin) < opts.epsilon.value()){ break;}
}
} while((itt_ctr2 < 3 || itt_ctr >= opts.maxNumItt.value()) &&
- (opts.approach.value() == string("tica")) && cum_itt < 100*opts.maxNumItt.value());
+ (opts.approach.value() == string("tica")) && cum_itt < newmaxitts);
if(itt_ctr>=opts.maxNumItt.value()){
cerr << " No convergence after " << cum_itt <<" steps "<<endl;
@@ -459,6 +466,7 @@ namespace Melodic{
melodat.set_IC(temp);
melodat.set_ICstats(scales);
melodat.sort();
+ message("Calculating T- and S-modes " << endl);
melodat.set_TSmode();
}
diff --git a/meloptions.h b/meloptions.h
index 1a6913bf2c97b2fbdab271164251ab88548f8572..35f4e65306b3e835ef600382e0e01175c6d5a78d 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -251,7 +251,7 @@ class MelodicOptions {
false, requires_argument),
fn_TconF(string("--Tconf"), string(""),
string(" F-contrast matrix across time-domain"),
- false, requires_argument),
+ false, requires_argument, false),
fn_Sdesign(string("--Sdes"), string(""),
string(" design matrix across subject-domain"),
false, requires_argument),
diff --git a/melreport.cc b/melreport.cc
index 53b9fc08a33a9627d4cc717771bfe9c9971c48a1..501e8152fdd019175ccf5ba136665d24c110606d 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -692,7 +692,7 @@ namespace Melodic{
newplot.add_label("dim. estimate");
}
- //CB's latest BUG!!! newplot.grid_swapdefault();
+ newplot.grid_swapdefault();
newplot.timeseries(what,
report.appendDir("EVplot.png"),
string("Eigenspectrum Analysis"),