fixed parameter problem in Downhill Simplex

DownhillSimplexOptimizer.cpp

@@ -16,13 +16,15 @@ using namespace optimization;
 
 DownhillSimplexOptimizer::DownhillSimplexOptimizer(OptLogBase *loger): SuperClass(loger)
 {
-	m_abort			= false;
-	m_simplexInitialized	= false;
-	m_alpha			= 1.0;
-	m_beta			= 0.5;
-	m_gamma			= 1.0;
-	m_rankdeficiencythresh= 0.01;
-	m_rankcheckenabled= false;
+  m_abort = false;
+  m_simplexInitialized = false;
+  //note that we use the negative value of alpha in amoeba
+  m_alpha = 1.0;
+  m_beta = 0.5;
+  // m_gamma was previously 1.0, but this actually deactivates extrapolation. We use 2.0 as suggested in Num. Recipes
+  m_gamma = 2.0;
+  m_rankdeficiencythresh = 0.01;
+  m_rankcheckenabled = false;
 }
 
 DownhillSimplexOptimizer::DownhillSimplexOptimizer(const DownhillSimplexOptimizer &opt) : SuperClass(opt)
@@ -164,7 +166,12 @@ int DownhillSimplexOptimizer::amoeba()
                                                 // successful steps 
     //int count_not_improve=0;
     //double last_best_value=1e99;	
- 
+
+  if (m_verbose)
+  {
+    std::cerr << "ndim: " << ndim << std::endl;
+  }
+
 	const int spts=ndim+1;						// number of vertices
     int i,j, max_val;
     
@@ -184,7 +191,6 @@ int DownhillSimplexOptimizer::amoeba()
 	*/
 	m_startTime = clock();
 
-    
     for (j=0;j<ndim;j++)
 	{ 											// get sum of vertex-coordinates
 		double sum=0.0;
@@ -192,7 +198,19 @@ int DownhillSimplexOptimizer::amoeba()
 		for (i=0;i<spts;i++)
 			sum += m_vertices[j][i];
 		psum[0][j]=sum;
+    std::cerr << sum << " " ;
+    }
+    std::cerr << std::endl;
+    
+  if (m_verbose)
+  {
+    std::cerr << "initial psum: ";
+    for (j=0;j<ndim;j++)
+    { 
+      std::cerr << psum[0][j] << " " ;
     }
+    std::cerr << std::endl;    
+  }
     
     for (;;)
 	{											// loop until terminating
@@ -202,11 +220,11 @@ int DownhillSimplexOptimizer::amoeba()
 			{
 				for(int v = 0; v < ndim ; v++)
 				{
-					std::cout << m_vertices[v][u] << " ";
+					std::cerr << m_vertices[v][u] << " ";
 				}
-				std::cout<< " " << m_y[u][0] << std::endl;
+				std::cerr<< " " << m_y[u][0] << std::endl;
 			}
-			std::cout << std::endl;
+			std::cerr << std::endl;
 		}
 
 
@@ -361,7 +379,12 @@ int DownhillSimplexOptimizer::amoeba()
 				break;
 			}
 		}
-        
+    
+    //informative output
+    if (m_verbose)
+    {    
+      std::cerr << "start new iteration with amotry -alpha, i.e., reflect worst point through simplex" << std::endl;
+    }
 		
 													// Begin a new iteration. 
 													// First extrapolate by the 
@@ -373,6 +396,11 @@ int DownhillSimplexOptimizer::amoeba()
 		ytry=amotry(psum,ihi,-m_alpha);
 		if (ytry < m_y[ilo][0])
 		{
+      //informative output
+      if (m_verbose)
+      {
+        std::cerr << "reflected point is better than best point, perform further extrapolation with gamma" << std::endl;
+      }
 			
 							 						// result is better than best
 													// point, try additional 
@@ -388,7 +416,12 @@ int DownhillSimplexOptimizer::amoeba()
 		{
 		
 			if (ytry >= m_y[inhi][0]) 
-			{ 			
+			{
+        //informative output
+        if (m_verbose)
+        {
+          std::cerr << "reflected point is worse then second worst, looking for intermediate point with beta" << std::endl;
+        }
 				
 														// The reflected point is worse
 														// than the second
@@ -405,8 +438,11 @@ int DownhillSimplexOptimizer::amoeba()
 					#ifdef OPT_DEBUG        
 						std::cout<<"Case three .. contract around lowest point" << std::endl;
 					#endif 
-
-
+          //informative output
+          if (m_verbose)
+          {
+            std::cerr <<  "Intermediate point is also worse, contract around current best point with factor 0.5." << std::endl;
+          }
 					for (i=0;i<spts;i++) 
 					{						 			//  contract around lowest 
 						// (best) point 
@@ -461,9 +497,25 @@ double DownhillSimplexOptimizer::amotry(matrix_type & psum, int ihi, double fac)
     matrix_type ptry(1,ndim);
     fac1=(1.0-fac)/ndim;
     fac2=fac1-fac;
-    
+      
     for (int j=0;j<ndim;j++) 
-    ptry[0][j]=psum[0][j]*fac1-m_vertices[j][ihi]*fac2;
+    {
+      ptry[0][j]=psum[0][j]*fac1-m_vertices[j][ihi]*fac2;
+    }
+
+    //informative output
+    if (m_verbose)
+    {
+      std::cerr << "amotry fac: " <<  fac << std::endl;
+      std::cerr << "fac1: " << fac1 << " fac2: " << fac2 << std::endl;
+      std::cerr << "ptry: ";
+      for (int j=0;j<ndim;j++) 
+      {
+        std::cerr << ptry[0][j] << " ";
+      }
+      std::cerr << std::endl;
+      
+    }
     
     if (checkParameters(!ptry)) 
 	{ 
@@ -490,12 +542,11 @@ double DownhillSimplexOptimizer::amotry(matrix_type & psum, int ihi, double fac)
 }
 
 
-void DownhillSimplexOptimizer::setDownhillParams(double alpha, double beta, double gamma)
+void DownhillSimplexOptimizer::setDownhillParams(const double alpha, const double beta, const double gamma)
 {
-
-	m_alpha = alpha;
-	m_beta = beta;
-	m_gamma = gamma;
+  m_alpha = alpha;
+  m_beta = beta;
+  m_gamma = gamma;
 }
 
 
@@ -510,6 +561,21 @@ void DownhillSimplexOptimizer::setRankCheckStatus(bool status)
     m_rankcheckenabled= status;
 }
 
+double DownhillSimplexOptimizer::getDownhillParameterAlpha() const
+{
+  return m_alpha;
+}
+
+double DownhillSimplexOptimizer::getDownhillParameterBeta() const
+{
+  return m_beta;
+}
+
+double DownhillSimplexOptimizer::getDownhillParameterGamma() const
+{
+  return m_gamma;
+}
+
 
 bool DownhillSimplexOptimizer::getRankCheckStatus()
 {

DownhillSimplexOptimizer.h

@@ -73,7 +73,7 @@ class DownhillSimplexOptimizer : public SimpleOptimizer
 		///	@param gamma
 		///	
 		///
-		void setDownhillParams(double alpha, double beta, double gamma);
+		void setDownhillParams(const double alpha, const double beta, const double gamma);
         
         ///
         /// Sets the rank deficiency threshold
@@ -87,6 +87,27 @@ class DownhillSimplexOptimizer : public SimpleOptimizer
         ///
         void setRankCheckStatus(bool status);
         
+        /**
+        * @brief get the downhill simplex parameter alpha, needed for the reflection
+        * @date 27-09-2012
+        * @author Alexander Freytag
+        */        
+        double getDownhillParameterAlpha() const;
+        
+        /**
+        * @brief get the downhill simplex parameter alpha, needed for the contraction
+        * @date 27-09-2012
+        * @author Alexander Freytag
+        */        
+        double getDownhillParameterBeta() const;
+        
+        /**
+        * @brief get the downhill simplex parameter alpha, needed for the expansion
+        * @date 27-09-2012
+        * @author Alexander Freytag
+        */        
+        double getDownhillParameterGamma() const;
+        
         ///
         /// Returns the status of rankcheck
         /// @retval true, if rankcheck is enabled
@@ -178,6 +199,7 @@ class DownhillSimplexOptimizer : public SimpleOptimizer
         /// Rank check status: if false, a rank check is disabled (default)
         ///
         bool m_rankcheckenabled;
+        
 
 };
 
@@ -192,5 +214,7 @@ class DownhillSimplexOptimizer : public SimpleOptimizer
 		///	@return (numerical) rank of A
 		///
 		unsigned int getRank(const optimization::matrix_type& A, double numZero);
+    
 
+    
 #endif