From 5e1f5f83a961c34045661495b7694971e9fa971c Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Thu, 15 Nov 2012 21:00:55 +0000
Subject: [PATCH] stable changes
---
meldata.cc | 30 +++++++++++++++++++++++-------
melhlprfns.cc | 9 ++++++++-
melhlprfns.h | 1 +
melica.cc | 12 ++++++------
meloptions.h | 10 +++++-----
5 files changed, 43 insertions(+), 19 deletions(-)
diff --git a/meldata.cc b/meldata.cc
index bc4f70a..2287756 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -327,8 +327,7 @@ namespace Melodic{
PPCA=tmpPPCA;
}
order = opts.pca_dim.value();
- if(opts.debug.value())
- message(endl << "Model order : "<<order<<endl<<endl);
+ dbgmsg(endl << "Model order : "<<order<<endl);
if (opts.paradigmfname.value().length()>0){
Matrix tmpPscales;
@@ -347,7 +346,6 @@ namespace Melodic{
outMsize("tmpPPCA",tmpPPCA); saveascii(tmpPPCA,"tmpPPCA");
outMsize("whiteMatrix",whiteMatrix); saveascii(whiteMatrix,"whiteMatrix");
outMsize("dewhiteMatrix",dewhiteMatrix); saveascii(dewhiteMatrix,"dewhiteMatrix");
- cerr << "Order: " << order << endl;
}
EV = AdjEV;
@@ -360,11 +358,16 @@ namespace Melodic{
WM.push_back(tmp);
}
else {
- for(int ctr = 0; ctr < numfiles; ctr++){
+
+ dbgmsg("Multi-Subject ICA");
+ for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
if(opts.joined_vn.value() && tmpvarnorm){
- tmpData=SP(tmpData,pow(ones(tmpData.Nrows(),1)*stdDev,-1));
+ dbgmsg("tmpData normalisation"<< endl);
+ dbgmsg("stdDev " << stdDev(1,2)<< endl);
+ dbgmsg("tmpData " << tmpData.SubMatrix(1,1,1,2)<< endl);
+ SP3(tmpData,pow(stdDev,-1));
}
// whiten (separate / joint)
Matrix newWM,newDWM;
@@ -381,10 +384,11 @@ namespace Melodic{
}
else{
if(opts.debug.value())
- message(" --mod_pca ");
+ dbgmsg(" --mod_pca ");
Matrix tmp1, tmp2;
tmp1 = whiteMatrix * alldat;
- tmp1 = remmean(tmp1,2) * tmpData.t();
+ remmean(tmp1,2);
+ tmp1 *= tmpData.t();
tmp2 = pinv(tmp1.t()).t();
std_pca(tmp1 * tmpData, RXweight, Corr, pcaE, pcaD);
calc_white(pcaE, pcaD, order, newWM, newDWM);
@@ -417,6 +421,7 @@ namespace Melodic{
}
} // void setup()
+
void MelodicData::setup_misc()
{
@@ -498,6 +503,7 @@ namespace Melodic{
remmean(Tdes);
}
+
void MelodicData::save()
{
@@ -617,6 +623,7 @@ namespace Melodic{
message("...done" << endl);
} //void save()
+
int MelodicData::remove_components()
{
message("Reading " << opts.filtermix.value() << endl)
@@ -691,6 +698,7 @@ namespace Melodic{
return 0;
} // int remove_components()
+
void MelodicData::create_RXweight()
{
message("Reading the weights for the covariance R_X from file "<< opts.segment.value() << endl);
@@ -700,6 +708,7 @@ namespace Melodic{
RXweight = tmpRX.matrix(Mask);
}
+
void MelodicData::est_smoothness()
{
if(Resels == 0){
@@ -736,6 +745,7 @@ namespace Melodic{
}
}
+
unsigned long MelodicData::standardise(volume<float>& mask, volume4D<float>& R)
{
unsigned long count = 0;
@@ -777,6 +787,7 @@ namespace Melodic{
return count;
}
+
float MelodicData::est_resels(volume4D<float> R, volume<float> mask)
{
message(" Estimating data smoothness ... ");
@@ -852,6 +863,7 @@ namespace Melodic{
return resels;
}
+
void MelodicData::create_mask(volume<float>& theMask)
{
if(opts.use_mask.value() && opts.maskfname.value().size()>0){ // mask provided
@@ -926,6 +938,7 @@ namespace Melodic{
}
} //void create_mask()
+
void MelodicData::sort()
{
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
@@ -1032,3 +1045,6 @@ namespace Melodic{
}
+
+
+
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 0758bc3..3692b7d 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -183,7 +183,14 @@ namespace Melodic{
Res = SP(in,Res);
}
return Res;
- } //Matrix SP
+ } //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;
diff --git a/melhlprfns.h b/melhlprfns.h
index 8e16fd1..5903ed2 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -49,6 +49,7 @@ namespace Melodic{
RowVector varnorm(Matrix& in, Matrix& Corr, int dim = 30, float level = 1.6);
Matrix SP2(const Matrix& in, const Matrix& weights, bool econ = 0);
+ void SP3(Matrix& in, const Matrix& weights);
RowVector Feta(int n1,int n2);
RowVector cumsum(const RowVector& Inp);
diff --git a/melica.cc b/melica.cc
index 72cb01f..e6a1458 100644
--- a/melica.cc
+++ b/melica.cc
@@ -57,9 +57,9 @@ namespace Melodic {
if(opts.approach.value() == string("tica"))
opts.maxNumItt.set_T(opts.rank1interval.value());
- rank1_old = melodat.get_dewhite()*redUMM;
- rank1_old = melodat.expand_dimred(rank1_old);
- rank1_old = krapprox(rank1_old,int(rank1_old.Nrows()/melodat.get_numfiles()));
+ rank1_old = melodat.get_dewhite()*redUMM;
+ rank1_old = melodat.expand_dimred(rank1_old);
+ rank1_old = krapprox(rank1_old,int(rank1_old.Nrows()/melodat.get_numfiles()));
do{// TICA loop
itt_ctr = 1;
@@ -67,7 +67,7 @@ namespace Melodic {
redUMM_old = redUMM;
//calculate IC estimates
tmpU = Data.t() * redUMM;
-
+
//update redUMM depending on nonlinearity
if(opts.nonlinearity.value()=="pow4"){
redUMM = (Data * pow(tmpU,3.0)) / samples - 3 * redUMM;
@@ -81,7 +81,7 @@ namespace Melodic {
Matrix hyptanh;
hyptanh = tanh(opts.nlconst1.value()*tmpU);
redUMM = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
- sum(1-pow(hyptanh,2),1),redUMM))/samples;
+ sum(1-pow(hyptanh,2),1),redUMM))/samples;
}
if(opts.nonlinearity.value()=="gauss"){
Matrix tmpUsq;
@@ -98,10 +98,10 @@ namespace Melodic {
// orthogonalize the unmixing-matrix
symm_orth(redUMM);
-
if(opts.approach.value() == string("tica")){
message("");
}
+
//termination condition : angle between old and new < epsilon
minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_old)).Minimum();
diff --git a/meloptions.h b/meloptions.h
index 7e10aeb..9b06b25 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -185,12 +185,12 @@ class MelodicOptions {
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"), true,
+ joined_whiten(string("--sep_whiten"), false,
string("switch on separate whitening"),
- false, no_argument),
+ false, no_argument, false),
joined_vn(string("--sep_vn"), true,
- string("switch off joined variance nomalisation"),
- false, no_argument),
+ string("switch on separate variance nomalisation"),
+ false, no_argument, false),
dr_pca(string("--mod_pca"), true,
string("switch off modified PCA for concat ICA"),
false, no_argument, false),
@@ -203,7 +203,7 @@ class MelodicOptions {
approach(string("-a,--approach"), string("symm"),
string("approach for decomposition, 2D: defl, symm (default), 3D: tica (default), concat"),
false, requires_argument),
- nonlinearity(string("--nl"), string("pow3"),
+ nonlinearity(string("--nl"), string("tanh"),
string("\tnonlinearity: gauss, tanh, pow3, pow4"),
false, requires_argument),
varnorm(string("--vn,--varnorm"), true,
--
GitLab