simplex_algorithm.h

Go to the documentation of this file.

//=============================================================================================================
 
#ifndef SIMPLEXALGORITHM_H
#define SIMPLEXALGORITHM_H
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "math_global.h"
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
 
//=============================================================================================================
// DEFINE NAMESPACE UTILSLIB
//=============================================================================================================
 
namespace UTILSLIB
{
 
//=============================================================================================================
class MATHSHARED_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 stol,
                                 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_lo,
                                                     double fval_hi,
                                                     double par_diff));
 
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>
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;
}
 
//=============================================================================================================
 
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 stol,
                                        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_lo,
                                                            double fval_hi,
                                                            double par_diff))
{
    constexpr int MIN_STOL_LOOP = 5;
    int   ndim = p.cols();
    int   i,ilo,ihi,inhi;
    int   mpts = ndim+1;
    T ytry,ysave,rtol;
    double dsum;
    Eigen::Matrix<T,Eigen::Dynamic, 1> psum(ndim);
    bool  result = true;
    int   count = 0;
    int   loop  = 1;
 
    neval = 0;
    psum = p.colwise().sum();
 
    constexpr T kAlpha = static_cast<T>(1.0);
    constexpr T kBeta  = static_cast<T>(0.5);
    constexpr T kGamma = static_cast<T>(2.0);
 
    if (report_func != nullptr && report > 0)
        report_func(0, static_cast<Eigen::Matrix<T,Eigen::Dynamic, 1>>(p.row(0)), -1.0, -1.0, 0.0);
 
    dsum = 0.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 != nullptr) {
            if (!report_func(loop, static_cast<Eigen::Matrix<T,Eigen::Dynamic, 1>>(p.row(ilo)),
                             y[ilo], y[ihi], std::sqrt(dsum))) {
                qWarning("Iteration interrupted.");
                result = false;
                break;
            }
            count = 0;
        }
        if (rtol < ftol) break;
        if (neval >=  max_eval) {
            qWarning("Maximum number of evaluations exceeded.");
            result  =  false;
            break;
        }
        if (stol > 0) {  /* Has the simplex collapsed? */
            dsum = (p.row(ilo) - p.row(ihi)).squaredNorm();
            if (loop > MIN_STOL_LOOP && std::sqrt(dsum) < stol)
                break;
        }
        ytry = tryit<T>(p,y,psum,func,user_data,ihi,neval,-kAlpha);
        if (ytry <= y[ilo])
            tryit<T>(p,y,psum,func,user_data,ihi,neval,kGamma);
        else if (ytry >= y[inhi]) {
            ysave = y[ihi];
            ytry = tryit<T>(p,y,psum,func,user_data,ihi,neval,kBeta);
            if (ytry >= ysave) {
                for (i = 0; i < mpts; i++) {
                    if (i !=  ilo) {
                        psum = static_cast<T>(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;
}
 
} //NAMESPACE
 
#endif // SIMPLEXALGORITHM_H