simplex_algorithm.h

Go to the documentation of this file.

//=============================================================================================================
#ifndef SIMPLEXALGORITHM_H
#define SIMPLEXALGORITHM_H
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "utils_global.h"
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
 
#define ALPHA 1.0
#define BETA 0.5
#define GAMMA 2.0
 
//=============================================================================================================
// DEFINE NAMESPACE MNELIB
//=============================================================================================================
 
namespace UTILSLIB
{
 
//=============================================================================================================
class UTILSSHARED_EXPORT SimplexAlgorithm
{
 
protected:
    //=========================================================================================================
    SimplexAlgorithm();
 
public:
    //=========================================================================================================
    template <typename T>
    static bool simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
                                 Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
                                 T ftol,
                                 T (*func)(const Eigen::Matrix<T,Eigen::Dynamic, 1>& x, const void *user_data),
                                 const void *user_data,
                                 int max_eval,
                                 int &neval,
                                 int report,
                                 bool (*report_func)(int loop, const Eigen::Matrix<T,Eigen::Dynamic, 1>& fitpar, double fval));
 
private:
 
    template <typename T>
    static T tryit(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &p,
                   Eigen::Matrix<T,Eigen::Dynamic, 1> &y,
                   Eigen::Matrix<T,Eigen::Dynamic, 1> &psum,
                   T (*func)(  const Eigen::Matrix<T,Eigen::Dynamic, 1> &x,const void *user_data),                     /* The function to be evaluated */
                   const void *user_data,                                      /* Data to be passed to the above function in each evaluation */
                   int   ihi,
                   int &neval,
                   T fac);
};
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
template <typename T>
bool SimplexAlgorithm::simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
                                        Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
                                        T ftol,
                                        T (*func)(const Eigen::Matrix<T,Eigen::Dynamic, 1>& x, const void *user_data),
                                        const void *user_data,
                                        int max_eval,
                                        int &neval,
                                        int report,
                                        bool (*report_func)(int loop, const Eigen::Matrix<T,Eigen::Dynamic, 1>& fitpar, double fval))
{
    int   ndim = p.cols();  /* Number of variables */
    int   i,ilo,ihi,inhi;
    int   mpts = ndim+1;
    T ytry,ysave,rtol;
    Eigen::Matrix<T,Eigen::Dynamic, 1> psum(ndim);
    bool  result = true;
    int   count = 0;
    int   loop  = 1;
 
    neval = 0;
    psum = p.colwise().sum();
 
    if (report_func != NULL && report > 0) {
        report_func(0,static_cast< Eigen::Matrix<T,Eigen::Dynamic, 1> >(p.row(0)),-1.0);
    }
 
    for (;;count++,loop++) {
        ilo = 1;
        ihi  =  y[1]>y[2] ? (inhi = 2,1) : (inhi = 1,2);
        for (i = 0; i < mpts; i++) {
            if (y[i]  <  y[ilo])
                ilo = i;
            if (y[i] > y[ihi]) {
                inhi = ihi;
                ihi = i;
            } else if (y[i] > y[inhi])
                if (i !=  ihi)
                    inhi = i;
        }
        rtol = 2.0*std::fabs(y[ihi]-y[ilo])/(std::fabs(y[ihi])+std::fabs(y[ilo]));
        /*
        * Report that we are proceeding...
        */
        if (count == report && report_func != NULL) {
            if (!report_func(loop,static_cast< Eigen::Matrix<T,Eigen::Dynamic, 1> >(p.row(ilo)),y[ilo])) {
                qCritical("Interation interrupted.");
                result = false;
                break;
            }
            count = 0;
        }
        if (rtol < ftol) break;
        if (neval >=  max_eval) {
            qCritical("Maximum number of evaluations exceeded.");
            result  =  false;
            break;
        }
        ytry = tryit<T>(p,y,psum,func,user_data,ihi,neval,-ALPHA);
        if (ytry <= y[ilo])
            tryit<T>(p,y,psum,func,user_data,ihi,neval,GAMMA);
        else if (ytry >= y[inhi]) {
            ysave = y[ihi];
            ytry = tryit<T>(p,y,psum,func,user_data,ihi,neval,BETA);
            if (ytry >= ysave) {
                for (i = 0; i < mpts; i++) {
                    if (i !=  ilo) {
                        psum = 0.5 * ( p.row(i) + p.row(ilo) );
                        p.row(i) = psum;
                        y[i] = (*func)(psum,user_data);
                    }
                }
                neval +=  ndim;
                psum = p.colwise().sum();
            }
        }
    }
 
    return result;
}
 
//=============================================================================================================
 
template <typename T>
T SimplexAlgorithm::tryit(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &p,
                          Eigen::Matrix<T,Eigen::Dynamic, 1> &y,
                          Eigen::Matrix<T,Eigen::Dynamic, 1> &psum,
                          T (*func)(  const Eigen::Matrix<T,Eigen::Dynamic, 1> &x,const void *user_data),                     /* The function to be evaluated */
                          const void *user_data,                                      /* Data to be passed to the above function in each evaluation */
                          int   ihi,
                          int &neval,
                          T fac)
{
    int ndim = p.cols();
    T fac1,fac2,ytry;
 
    Eigen::Matrix<T,Eigen::Dynamic, 1> ptry(ndim);
 
    fac1 = (1.0-fac)/ndim;
    fac2 = fac1-fac;
 
    ptry = psum * fac1 - p.row(ihi).transpose() * fac2;
 
    ytry = (*func)(ptry,user_data);
    ++neval;
 
    if (ytry < y[ihi]) {
        y[ihi] = ytry;
 
        psum += ptry - p.row(ihi).transpose();
        p.row(ihi) = ptry;
    }
 
    return ytry;
}
} //NAMESPACE
 
#undef ALPHA
#undef BETA
#undef GAMMA
 
#endif // SIMPLEXALGORITHM_H