From 81c3c734aac872e4ba3911d2624582b2a1d76f61 Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Thu, 9 Aug 2007 11:10:38 +0000
Subject: [PATCH] MM report option added
---
meldata.cc | 9 +++++++--
melica.cc | 9 ++++++---
melodic.cc | 9 ++++-----
3 files changed, 17 insertions(+), 10 deletions(-)
diff --git a/meldata.cc b/meldata.cc
index 0dc130d..02cfe74 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -413,16 +413,21 @@ namespace Melodic{
// ICadjust = IC;
}
else{
+
Matrix resids = stdev(Data - mixMatrix * IC);
for(int ctr=1;ctr<=resids.Ncols();ctr++)
if(resids(1,ctr) < 0.05)
resids(1,ctr)=1;
- stdNoisei = pow(stdev(Data - mixMatrix * IC)*
+ // stdNoisei = pow(stdev(Data - mixMatrix * IC)*
+ // std::sqrt((float)(Data.Nrows()-1))/
+ // std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
+ stdNoisei = pow(resids*
std::sqrt((float)(Data.Nrows()-1))/
std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
-
+
ColumnVector diagvals;
diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5);
+
save4D(SP(IC,diagvals*stdNoisei),opts.outputfname.value() + "_IC");
}
diff --git a/melica.cc b/melica.cc
index b5d9b7c..9956837 100644
--- a/melica.cc
+++ b/melica.cc
@@ -101,9 +101,11 @@ namespace Melodic{
}
//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);
- if((abs(minAbsSin) < opts.epsilon.value())&&
- (opts.approach.value()!=string("tica"))){ break;}
+ // if((abs(minAbsSin) < opts.epsilon.value())&&
+ // (opts.approach.value()!=string("tica"))){ break;}
+ if((abs(minAbsSin) < opts.epsilon.value())){ break;}
itt_ctr++;
} while(itt_ctr < opts.maxNumItt.value());
@@ -119,7 +121,7 @@ namespace Melodic{
temp = melodat.reduce_dimred(temp);
redUMM = melodat.get_white() * temp;
- message(" change : " << abs(minAbsSin2) << endl);
+ message(" change : " << minAbsSin2 << endl);
if(abs(minAbsSin2) < opts.epsilonS.value() && abs(minAbsSin) < opts.epsilon.value()){ break;}
}
} while(
@@ -467,6 +469,7 @@ namespace Melodic{
melodat.set_mix(tmp);
melodat.set_IC(temp);
+
melodat.set_ICstats(scales);
melodat.sort();
diff --git a/melodic.cc b/melodic.cc
index d6d95d8..0d04f18 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -114,7 +114,7 @@ int main(int argc, char *argv[]){
opts.pca_dim.set_T(opts.pca_dim.value()-retry*opts.retrystep);
}
else{
- if(opts.pca_dim.value()-retry*opts.retrystep < melodat.data_dim()){
+ if(opts.pca_dim.value()+retry*opts.retrystep < melodat.data_dim()){
opts.pca_dim.set_T(opts.pca_dim.value()+retry*opts.retrystep);
}else{
leaveloop = true; //stupid, but break does not compile
@@ -234,7 +234,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
message(" done" << endl);
}
- message("Reading mixing matrix " << opts.filtermix.value());
+ message("Reading mixing matrix " << opts.filtermix.value() << " ... ");
mixMatrix = read_ascii_matrix(opts.filtermix.value());
if (mixMatrix.Storage()<=0) {
cerr <<" Please specify the mixing matrix correctly" << endl;
@@ -257,13 +257,12 @@ void mmonly(Log& logger, MelodicOptions& opts,
}
}
-
melodat.tempInfo = ICvolInfo;
melodat.set_mask(Mask);
melodat.set_mean(Mean);
melodat.set_IC(ICs);
melodat.set_mix(mixMatrix);
- fmixMatrix = calc_FFT(melodat.expand_mix(), opts.logPower.value());
+ fmixMatrix = calc_FFT(mixMatrix, opts.logPower.value());
melodat.set_fmix(fmixMatrix);
fmixMatrix = pinv(mixMatrix);
melodat.set_unmix(fmixMatrix);
@@ -272,7 +271,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
Matrix mmres;
Matrix pmaps;//(ICs);
- if(opts.perf_mm.value())
+ if(opts.perf_mm.value())
mmres = mmall(logger,opts,melodat,report,pmaps);
}
--
GitLab