source: Sophya/trunk/SophyaLib/NTools/simplex.h@ 3638

#ifndef SIMPLEX_SEEN
#define SIMPLEX_SEEN

#include "machdefs.h"
#include "tvector.h"
#include "generalfit.h"
#include <string>
#include <vector>

// ---------------------------------------------------------------
//             (C) UPS-LAL  R. Ansari , Aout 2004 
// ---------------------------------------------------------------
// Classe de minimisation (optimisation) non lineaire dans
// un espace multidimensionnel suivant la methode des simplex.
// Simplex : Figure geometrique formee de N+1 sommets (points) 
//  (points) dans un espace a N dimension 
// La methode d'optimisation codee ici dans la classe MinZSimplex 
// est basee sur celle decrite dans Numerical Recipes ,
// Chapitre X, Minimization or Maximization of Functions 
// L'algorithme est ameliore sur deux aspects :
//  - Test d'arret plus complexe 
//  - Ajout d'une transformation du simplex  (ExpandHigh)
//    Reflection,ReflecExpand,ContractHigh,ContractLow+ExpandHigh 
// ----------------------------------------------------------------

using namespace std;

namespace SOPHYA {

//! Interface definition for multivariable function (used by SimplexMinmizer)
class MinZFunction {
public:
  //! Constructor with nVar defining the number of variables
  MinZFunction(unsigned int nVar); 
  virtual ~MinZFunction();
  //! Returns the value of the function for the point defined by xp[]
  virtual double Value(double const xp[])=0;
  //! Returns the number of variables (dimension of xp[]) in method Value()
  inline int     NVar() const {return mNVar;}  
protected:
  const int mNVar;  //!< nombre de variables f(x,y,z,...)  
};

//!  Wrapper class implementing MinZFunction for a GeneralXi2 associated with a GeneralFitData
class MinZFuncXi2 : public MinZFunction {
public: 
  MinZFuncXi2(GeneralXi2* gxi2, GeneralFitData* gd);
  virtual ~MinZFuncXi2(); 
  virtual double Value(double const xp[]);
protected:
  GeneralXi2 *mGXi2;
  GeneralFitData *mGData;
};

//! Class implementing the Simplex minimisation method 
class MinZSimplex {
public:
  static int     AutoTest(int tsel=-1, int prtlev=0);

                 MinZSimplex(MinZFunction *mzf);
  virtual        ~MinZSimplex();

  //! Return the parameter space dimension 
  inline int     NDim() { return mZF->NVar(); }
  // Simplex initial 
  //! Defines the initial point (center of simplex figure)
  inline void    SetInitialPoint(Vector& point) { mPoint0 = point; }
  //! Defines the step along each dimension to construct the simplex from initial point
  inline void    SetInitialStep(Vector& step) { mStep0 = step; }

  //! Set the info/debug print level  
  inline void    SetPrtLevel(int lev=0) { mPrt = lev; }
  //! Return the current print level
  inline int     PrtLevel() { return mPrt; }

  //! Set the maximum number of iteration
  inline void    SetMaxIter(int max = 100000) { mMaxIter = max; }
  //! Return the current  max iter 
  inline int     MaxIter() { return mMaxIter; }
  //! Return the number of iterations performed
  inline int     NbIter() { return mIter; }
  //! Return the stop reason
  inline int     StopReason() { return mStop; }
  //! Return the stop reason and a description string (\b s)
         int     StopReason(string& s);

  // On minimise f(x) f=mZF->Value() , 
  //             f_max = max(f) sur simplex , f_min = min(f) sur simplex 
  //             fm = (abs(f_max)+abs(f_min))
  //             [Delta f] = abs(f_max-f_min)
  //             [Delta f/f]simplex = 2.*Delta f / fm
  //             fm2 = (abs(f_max)+abs(f_max(iter-1)))
  //             [Delta f_max/f_max]iter = [f_max(iter-1)-f_max]/fm2 
  // Test d'arret :  
  //  fm < mTol0                                                       OU 
  //  [Delta f/f]simplex              < mTol1   mRep1 fois de suite    OU 
  //  [Delta f_max/f_max]iter         < mTol2   mRep2 fois de suite 

  //! Define the tolerances for the various convergence tests
  inline void    SetStopTolerance(double tol0=1.e-39, double tol1 = 1.e-3, int rep1=5, 
                                  double tol2=1.e-4, int rep2=5)   
                 { mTol0 = tol0; mTol1 = tol1; mRep1 = rep1; mTol2 = tol2, mRep2 = rep2; }

  // Controle des facteurs de deformation du simplex:
  // Alpha  = Facteur d'homothetie pour la reflexion simple (Reflection)
  // Gamma  = Facteur d'homothetie pour la reflexion+expansion (ReflecExpand)
  // Beta   = Facteur d'homothetie pour la contraction pour le sommet haut f_max (ContractHigh)
  // Beta2  = Facteur d'homothetie pour la contraction vers le sommet bas f_min (ContractLow)
  // Gamma2 = Facteur d'homothetie pour la l'extension pour le sommet haut ExpandHigh

  //! Define the similarity (homothetic) factors for the different simplex transformations
  inline void    SetControls(double alpha=1., double beta=0.5, double beta2=0.5, 
                             double gamma=2.0, double gamma2=2.0)
                 { mAlpha = alpha; mBeta = beta; mBeta2 = beta2;  mGamma = gamma; mGamma2 = gamma2;}  

  //! Return the the homothetic factor for Reflection
  inline double  Alpha() { return mAlpha; }
  //! Return the the homothetic factor for ContractHigh (contraction away from high point)
  inline double  Beta() { return mBeta; }
  //! Return the the homothetic factor for ContractLow (contraction toward the low point)
  inline double  Beta2() { return mBeta2; }
  //! Return the the homothetic factor for ReflecExpand (reflection+expansion)
  inline double  Gamma() { return mGamma; }
  //! Return the the homothetic factor for ExpandHigh (expansion along high point)
  inline double  Gamma2() { return mGamma2; }

  // Fonction de minimisation - Rc=0, OK, convergence , Rc>0, Non convergence 
  // Voir StopReason(string& s)
  virtual int    Minimize(Vector& fpoint);
  
protected:
  int            FindMinMax12(Vector& fval, int& ilo, int& ihi, int& inhi);
  inline double  Value(Vector& point) { return mZF->Value(point.Data()); }
  
  MinZFunction* mZF;
  int mMaxIter, mIter;
  int mStop;
  int mPrt;

  Vector mPoint0;
  Vector mStep0;
  //  vector< Vector > simplex0;
  //  vector< Vector > simplex_cur;
  double mAlpha, mBeta, mBeta2, mGamma, mGamma2;
  double mTol0, mTol1, mTol2;
  int mRep1,mRep2;
};

} // namespace SOPHYA

#endif