From c8a80ec36386bdc69340be8e4d121ec13a0a6d14 Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Tue, 12 Aug 2003 17:39:13 +0000
Subject: [PATCH] defl approach modified
---
melica.cc | 101 +++++++++++++++++++++++++++++++++++++++++++++-----
meloptions.cc | 4 +-
2 files changed, 95 insertions(+), 10 deletions(-)
diff --git a/melica.cc b/melica.cc
index 981e76e..1d5f87d 100644
--- a/melica.cc
+++ b/melica.cc
@@ -121,18 +121,21 @@ namespace Melodic{
//redUMM = zeros(dim); // got to start somewhere
redUMM = melodat.get_white()*
unifrnd(melodat.get_white().Ncols(),opts.numICs.value());
+
+
+ redUMM = zeros(melodat.get_white().Nrows(),opts.numICs.value());
+ Matrix guess;
+
int guesses=0;
if(opts.guessfname.value().size()>1){
message(" Use columns in " << opts.guessfname.value() << " as initial values for the mixing matrix " <<endl);
- Matrix guess;
guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
guesses = guess.Ncols();
- redUMM.Columns(1,guesses) = guess;
}
int ctrIC = 1;
int numRestart = 0;
-
+
while(ctrIC<=opts.numICs.value()){
message(" Extracting IC " << ctrIC << " ... ");
@@ -140,15 +143,19 @@ namespace Melodic{
ColumnVector w_old;
ColumnVector tmpU;
if(ctrIC <= guesses){
- w = redUMM.Column(ctrIC);}
+ w = guess.Column(ctrIC);}
else{
w = unifrnd(dim,1);}
-
+
w = w - redUMM * redUMM.t() * w;
w = w / norm2(w);
+
+ w_old = zeros(w.Nrows(),1);
+
int itt_ctr = 1;
do{
+
w_old = w;
tmpU = Data.t() * w;
@@ -185,8 +192,10 @@ namespace Melodic{
w = w / norm2(w);
//termination condition : angle between old and new < epsilon
- if(norm2(w-w_old) < opts.epsilon.value() ||
- norm2(w+w_old) < opts.epsilon.value()){break;}
+ if((norm2(w-w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10) ||
+ (norm2(w+w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10)){
+ break;
+ }
//cout << norm2(w-w_old) << " " << norm2(w+w_old) << endl;
itt_ctr++;
} while(itt_ctr <= opts.maxNumItt.value());
@@ -223,9 +232,83 @@ namespace Melodic{
void MelodicICA::ica_maxent(const Matrix &Data)
{
- cerr << " Not fully implemented yet " << endl;
- // cout << "maxent" <<endl;
+ // based on Aapo Hyvärinen's fastica method
+ // see www.cis.hut.fi/projects/ica/fastica/
+
+ message(" MAXENT " << endl);
+
+
+ //initialize matrices
+ Matrix redUMM_old;
+ Matrix tmpU;
+ Matrix gtmpU;
+ double lambda = 0.015/std::log((double)melodat.get_white().Ncols());
+
+ //srand((unsigned int)timer(NULL));
+
+ redUMM = melodat.get_white()*
+ unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
+
+ if(opts.guessfname.value().size()>1){
+ message(" Use columns in " << opts.guessfname.value()
+ << " as initial values for the mixing matrix " <<endl);
+ Matrix guess ;
+ guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
+ redUMM.Columns(1,guess.Ncols()) = guess;
+ }
+
+ // symm_orth(redUMM);
+
+ int itt_ctr=1;
+ double minAbsSin = 1.0;
+
+ Matrix Id;
+ Id = Identity(redUMM.Ncols());
+
+ cerr << " nonlinearity : " << opts.nonlinearity.value() << endl;
+
+ do{ // da loop!!!
+
+ redUMM_old = redUMM;
+
+ //calculate IC estimates
+ tmpU = Data.t() * redUMM;
+
+ if(opts.nonlinearity.value()=="tanh"){
+ //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;
+ gtmpU = tanh(opts.nlconst1.value()*tmpU);
+ redUMM = redUMM + lambda*(Id+(1-2*gtmpU.t()*tmpU))*redUMM;
+ }
+ if(opts.nonlinearity.value()=="gauss"){
+ gtmpU = pow(1+exp(-(opts.nlconst2.value()/2) * tmpU),-1);
+ redUMM = redUMM + lambda*(Id - (gtmpU.t()-tmpU.t())*tmpU)*redUMM;
+ }
+
+ // redUMM = redUMM + ( tmpU *tmpU.t()) * redUMM;
+
+ //termination condition : angle between old and new < epsilon
+ minAbsSin = abs(1 - diag(abs(redUMM.t()*redUMM_old)).Minimum());
+ message(" Step no. " << itt_ctr << " change : " << minAbsSin << endl);
+ if(abs(minAbsSin) < opts.epsilon.value()){ break;}
+
+ itt_ctr++;
+ } while(itt_ctr < opts.maxNumItt.value());
+
+ if(itt_ctr>=opts.maxNumItt.value()){
+ cerr << " No convergence after " << itt_ctr <<" steps "<<endl;
+ } else {
+ message(" Convergence after " << itt_ctr <<" steps " << endl << endl);
+ no_convergence = false;
+ {Matrix temp(melodat.get_dewhite() * redUMM);
+ melodat.set_mix(temp);}
+ {Matrix temp(redUMM.t()*melodat.get_white());
+ melodat.set_unmix(temp);}
+ }
}
+
void MelodicICA::ica_jade(const Matrix &Data)
diff --git a/meloptions.cc b/meloptions.cc
index 743aa23..cb8719a 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -193,7 +193,7 @@ void MelodicOptions::print_usage(int argc, char *argv[])
{
print_version();
options.usage();
- /* cout << "Usage: " << argv[0] << " ";
+ /* cout << "Usage: " << argv[0] << " ";
Have own usage output here
*/
@@ -202,12 +202,14 @@ void MelodicOptions::print_usage(int argc, char *argv[])
void MelodicOptions::print_version()
{
cout << "MELODIC Version " << version;
+ cout.flush();
}
void MelodicOptions::print_copyright()
{
print_version();
cout << endl << "Copyright (C) 1999-2002 University of Oxford " << endl;
+ cout.flush();
}
void MelodicOptions::status()
--
GitLab