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, typename CostFunc, typename ReportFunc>
    static bool simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
                                 Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
                                 T ftol,
                                 T stol,
                                 CostFunc&& func,
                                 int max_eval,
                                 int &neval,
                                 int report,
                                 ReportFunc&& report_func);
 
    //=========================================================================================================
    template <typename T, typename CostFunc>
    static bool simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
                                 Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
                                 T ftol,
                                 T stol,
                                 CostFunc&& func,
                                 int max_eval,
                                 int &neval);
 
private:
 
    template <typename T, typename CostFunc>
    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,
                   CostFunc&& func,
                   int   ihi,
                   int &neval,
                   T fac);
};
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
template <typename T, typename CostFunc>
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,
                          CostFunc&& func,
                          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);
    ++neval;
 
    if (ytry < y[ihi]) {
        y[ihi] = ytry;
 
        psum += ptry - p.row(ihi).transpose();
        p.row(ihi) = ptry;
    }
 
    return ytry;
}
 
//=============================================================================================================
 
template <typename T, typename CostFunc>
bool SimplexAlgorithm::simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
                                        Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
                                        T ftol,
                                        T stol,
                                        CostFunc&& func,
                                        int max_eval,
                                        int &neval)
{
    auto no_report = [](int, const Eigen::Matrix<T,Eigen::Dynamic,1>&, double, double, double) { return true; };
    return simplex_minimize<T>(p, y, ftol, stol, std::forward<CostFunc>(func), max_eval, neval, -1, no_report);
}
 
//=============================================================================================================
 
template <typename T, typename CostFunc, typename ReportFunc>
bool SimplexAlgorithm::simplex_minimize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& p,
                                        Eigen::Matrix<T,Eigen::Dynamic, 1>& y,
                                        T ftol,
                                        T stol,
                                        CostFunc&& func,
                                        int max_eval,
                                        int &neval,
                                        int report,
                                        ReportFunc&& report_func)
{
    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 > 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) {
            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,ihi,neval,-kAlpha);
        if (ytry <= y[ilo])
            tryit<T>(p,y,psum,func,ihi,neval,kGamma);
        else if (ytry >= y[inhi]) {
            ysave = y[ihi];
            ytry = tryit<T>(p,y,psum,func,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);
                    }
                }
                neval +=  ndim;
                psum = p.colwise().sum();
            }
        }
    }
    return result;
}
 
} //NAMESPACE
 
#endif // SIMPLEXALGORITHM_H