defl approach modified

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 Hyvrinen'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)

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()