From 2696355470c95b601a1c11516cb27217bd9682d7 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Mon, 16 Jul 2018 11:47:12 +0300 Subject: [PATCH 01/23] Rely on ciso646 and __cplusplus macro when detecting cxx --- CMakeLists.txt | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 53b2bf60..c44805d0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -138,12 +138,14 @@ if (WITH_THREADLOCAL AND NOT DEFINED HAVE_THREAD_LOCAL_STORAGE) endif () if (NLOPT_CXX OR NLOPT_PYTHON OR NLOPT_GUILE OR NLOPT_OCTAVE) - check_cxx_symbol_exists (_LIBCPP_VERSION string SYSTEM_HAS_LIBCPP) - if (SYSTEM_HAS_LIBCPP) + check_cxx_symbol_exists (__cplusplus ciso646 SYSTEM_HAS_CXX) + if (SYSTEM_HAS_CXX) check_cxx_compiler_flag ("-std=c++11" SUPPORTS_STDCXX11) if (SUPPORTS_STDCXX11) set (CMAKE_CXX_FLAGS "-std=c++11 ${CMAKE_CXX_FLAGS}") endif () + else() + message(FATAL_ERROR "The compiler doesn't support CXX." ) endif () endif () From b5a5afceaa5e606524aeb70d1c8c729c97ddb845 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Mon, 16 Jul 2018 12:04:39 +0300 Subject: [PATCH 02/23] Add CXX11 flag to cmake --- CMakeLists.txt | 3 ++- nlopt_config.h.in | 15 +++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c44805d0..ab1fa8f9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -143,9 +143,10 @@ if (NLOPT_CXX OR NLOPT_PYTHON OR NLOPT_GUILE OR NLOPT_OCTAVE) check_cxx_compiler_flag ("-std=c++11" SUPPORTS_STDCXX11) if (SUPPORTS_STDCXX11) set (CMAKE_CXX_FLAGS "-std=c++11 ${CMAKE_CXX_FLAGS}") + set (NLOPT_CXX11 ON) endif () else() - message(FATAL_ERROR "The compiler doesn't support CXX." ) + message (FATAL_ERROR "The compiler doesn't support CXX.") endif () endif () diff --git a/nlopt_config.h.in b/nlopt_config.h.in index c5dd6858..464c60b3 100644 --- a/nlopt_config.h.in +++ b/nlopt_config.h.in @@ -1,16 +1,16 @@ /*============================================================================== # NLOPT CMake configuration file -# -# NLopt is a free/open-source library for nonlinear optimization, providing -# a common interface for a number of different free optimization routines -# available online as well as original implementations of various other +# +# NLopt is a free/open-source library for nonlinear optimization, providing +# a common interface for a number of different free optimization routines +# available online as well as original implementations of various other # algorithms -# WEBSITE: http://ab-initio.mit.edu/wiki/index.php/NLopt +# WEBSITE: http://ab-initio.mit.edu/wiki/index.php/NLopt # AUTHOR: Steven G. Johnson # # This config.cmake.h.in file was created to compile NLOPT with the CMAKE utility. # Benoit Scherrer, 2010 CRL, Harvard Medical School -# Copyright (c) 2008-2009 Children's Hospital Boston +# Copyright (c) 2008-2009 Children's Hospital Boston # # Minor changes to the source was applied to make possible the compilation with # Cmake under Linux/Win32 @@ -150,6 +150,9 @@ /* Define if compiled including C++-based routines */ #cmakedefine NLOPT_CXX +/* Define if compiled including C++11-based routines */ +#cmakedefine NLOPT_CXX11 + /* Define to empty if `const' does not conform to ANSI C. */ #undef const From 9afec621b6a3de5296d01e6e7831751f9301438b Mon Sep 17 00:00:00 2001 From: sovrasov Date: Sat, 14 Jul 2018 13:55:47 +0300 Subject: [PATCH 03/23] Add a stub for AGS algrithm --- CMakeLists.txt | 3 + src/algs/ags/ags.cc | 22 ++ src/algs/ags/ags.h | 72 ++++++ src/algs/ags/data_types.hpp | 75 ++++++ src/algs/ags/evolvent.cc | 370 +++++++++++++++++++++++++++ src/algs/ags/evolvent.hpp | 48 ++++ src/algs/ags/local_optimizer.cc | 135 ++++++++++ src/algs/ags/local_optimizer.hpp | 44 ++++ src/algs/ags/solver.cc | 415 +++++++++++++++++++++++++++++++ src/algs/ags/solver.hpp | 102 ++++++++ src/api/general.c | 13 +- src/api/nlopt.h | 38 +-- src/api/optimize.c | 15 ++ 13 files changed, 1330 insertions(+), 22 deletions(-) create mode 100644 src/algs/ags/ags.cc create mode 100644 src/algs/ags/ags.h create mode 100644 src/algs/ags/data_types.hpp create mode 100644 src/algs/ags/evolvent.cc create mode 100644 src/algs/ags/evolvent.hpp create mode 100644 src/algs/ags/local_optimizer.cc create mode 100644 src/algs/ags/local_optimizer.hpp create mode 100644 src/algs/ags/solver.cc create mode 100644 src/algs/ags/solver.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ab1fa8f9..6d3af37f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -204,6 +204,8 @@ set (NLOPT_SOURCES if (NLOPT_CXX) list (APPEND NLOPT_SOURCES src/algs/stogo/global.cc src/algs/stogo/linalg.cc src/algs/stogo/local.cc src/algs/stogo/stogo.cc src/algs/stogo/tools.cc src/algs/stogo/global.h src/algs/stogo/linalg.h src/algs/stogo/local.h src/algs/stogo/stogo_config.h src/algs/stogo/stogo.h src/algs/stogo/tools.h) + list (APPEND NLOPT_SOURCES src/algs/ags/data_types.hpp src/algs/ags/evolvent.hpp src/algs/ags/evolvent.cc src/algs/ags/solver.hpp src/algs/ags/solver.cc + src/algs/ags/local_optimizer.hpp src/algs/ags/local_optimizer.cc src/algs/ags/ags.h src/algs/ags/ags.cc) endif () install (FILES ${NLOPT_HEADERS} DESTINATION ${RELATIVE_INSTALL_INCLUDE_DIR}) @@ -222,6 +224,7 @@ target_include_directories (${nlopt_lib} PRIVATE ${PROJECT_BINARY_DIR}/src/api ${PROJECT_BINARY_DIR} src/algs/stogo + src/algs/ags src/util src/algs/direct src/algs/cdirect diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc new file mode 100644 index 00000000..e3db8b10 --- /dev/null +++ b/src/algs/ags/ags.cc @@ -0,0 +1,22 @@ +// A C-callable front-end to the AGS global-optimization library. +// -- Vladislav Sovrasov + +#include "ags.h" +#include "solver.hpp" + +/* +int ags_minimize(int n, + objective_func fgrad, void *data, + double *x, double *minf, + const double *l, const double *u, +#ifdef NLOPT_UTIL_H + nlopt_stopping *stop, +#else + long int maxeval, double maxtime, +#endif + int nrandom) + */ +int ags_minimize() +{ + return 1; +} diff --git a/src/algs/ags/ags.h b/src/algs/ags/ags.h new file mode 100644 index 00000000..0d03cfd6 --- /dev/null +++ b/src/algs/ags/ags.h @@ -0,0 +1,72 @@ +/* A C-callable front-end to the AGS global-optimization library. + -- Vladislav Sovrasov */ + +#ifndef AGS_H +#define AGS_H + +#include "nlopt-util.h" + +#ifdef __cplusplus +extern "C" { +#endif + +/* search for the global minimum of the function fgrad(n, x, grad, data) + inside a simple n-dimensional hyperrectangle. + + Input: + + n: dimension of search space (number of decision variables) + + fgrad: the objective function of the form fgrad(n, x, grad, data), + returning the objective function at x, where + n: dimension of search space + x: pointer to array of length n, point to evaluate + grad: if non-NULL, an array of length n which + should on return be the gradient d(fgrad)/dx + [ if NULL, grad should be ignored ] + data: arbitrary data pointer, whose value is the + data argument of ags_minimize + + data: arbitrary pointer to any auxiliary data needed by fgrad + + l, u: arrays of length n giving the lower and upper bounds of the + search space + + maxeval: if nonzero, a maximum number of fgrad evaluations + maxtime: if nonzero, a maximum time (in seconds) + -- REPLACED in NLopt by nlopt_stopping *stop + + nrandom: number of randomized search points to use per box, + in addition to 2*n+1 deterministic search points + (0 for a deterministic algorithm). + + Output: + + minf: the minimum value of the objective function found + + x: pointer to array of length n, giving the location of the minimum + + Return value: 0 if no minimum found, 1 otherwise. + + */ + +int ags_minimize(); +/* +int ags_minimize(int n, + objective_func fgrad, void *data, + double *x, double *minf, + const double *l, const double *u, +#ifdef NLOPT_UTIL_H + nlopt_stopping *stop, +#else + long int maxeval, double maxtime, +#endif + int nrandom); +*/ +extern int ags_verbose; /* set to nonzero for verbose output */ + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/algs/ags/data_types.hpp b/src/algs/ags/data_types.hpp new file mode 100644 index 00000000..94524284 --- /dev/null +++ b/src/algs/ags/data_types.hpp @@ -0,0 +1,75 @@ +/* +Copyright (C) 2018 Sovrasov V. - All Rights Reserved + * You may use, distribute and modify this code under the + * terms of the MIT license. + * You should have received a copy of the MIT license with + * this file. If not visit https://opensource.org/licenses/MIT +*/ +#pragma once + +#include +#include + +#define NLP_SOLVER_ERROR(msg) throw std::runtime_error(std::string(msg)) +#define NLP_SOLVER_ASSERT(expr, msg) if(!(expr)) NLP_SOLVER_ERROR(msg) + +namespace ags +{ + +const unsigned solverMaxDim = 5; +const unsigned solverMaxConstraints = 10; + +template +class IGOProblem +{ +public: + ~IGOProblem() {} + + virtual fptype Calculate(const fptype* y, int fNumber) const = 0; + virtual int GetConstraintsNumber() const = 0; + virtual int GetDimension() const = 0; + virtual void GetBounds(fptype* left, fptype* right) const = 0; + virtual int GetOptimumPoint(fptype* y) const = 0; + virtual fptype GetOptimumValue() const = 0 ; +}; + +struct Trial +{ + double x; + double y[solverMaxDim]; + double g[solverMaxConstraints + 1]; + int idx; + Trial() {} + Trial(double _x) : x(_x) {} +}; + +struct Interval +{ + Trial pl; + Trial pr; + double R; + double delta; + Interval() {} + Interval(const Trial& _pl, const Trial& _pr) : pl(_pl), pr(_pr) {} +}; + +struct CompareIntervals +{ + bool operator() (const Interval* i1, const Interval* i2) const + { + return i1->pl.x < i2->pl.x; + } +}; + +class CompareByR +{ +public: + bool operator() (const Interval* i1, const Interval* i2) const + { + return i1->R < i2->R; + } +}; + + + +} diff --git a/src/algs/ags/evolvent.cc b/src/algs/ags/evolvent.cc new file mode 100644 index 00000000..759995cc --- /dev/null +++ b/src/algs/ags/evolvent.cc @@ -0,0 +1,370 @@ +/* +Copyright (C) 2018 Sovrasov V. - All Rights Reserved + * You may use, distribute and modify this code under the + * terms of the MIT license. + * You should have received a copy of the MIT license with + * this file. If not visit https://opensource.org/licenses/MIT +*/ +#include "evolvent.hpp" +#include +#include +#include + +using namespace ags; + +void mapd(double x, int m, double* y, int n, int key = 1); /* map x to y */ +void invmad(int, double [], int, int *, double [], int, int); /* map y to x */ +void xyd(double *xx, int m, double y[], int n); /* get preimage */ + +Evolvent::Evolvent() +{ + mIsInitialized = false; +} + +Evolvent::~Evolvent() +{ +} + +Evolvent::Evolvent(int dimension, int tightness, const double* lb, const double* ub, MapType type) +{ + assert(tightness > 2); + + mDimension = dimension; + mTightness = tightness; + mMapType = type; + + p2.resize(mDimension); + mShiftScalars.resize(mDimension); + mRho.resize(mDimension); + for (int i = 0; i < mDimension; i++) + { + mRho[i] = ub[i] - lb[i]; + mShiftScalars[i] = 0.5*(lb[i] + ub[i]); + } + + switch (mMapType) + { + case Simple: + mMapKey = 1; + break; + case Linear: + mMapKey = 2; + break; + case Noninjective: + mMapKey = 3; + break; + } + + mIsInitialized = true; +} + +void Evolvent::TransformToStandardCube(const double *y, double *z) +{ + for (int i = 0; i < mDimension; i++) + z[i] = (y[i] - mShiftScalars[i]) / mRho[i]; +} + +void Evolvent::TransformToSearchDomain(const double *y, double *z) +{ + for (int i = 0; i < mDimension; i++) + z[i] = mRho[i]*y[i] + mShiftScalars[i]; +} + +void Evolvent::GetImage(double x, double y[]) +{ + if(mDimension != 1) + mapd(x, mTightness, y, mDimension, mMapKey); + else + y[0] = x - 0.5; + + TransformToSearchDomain(y, y); +} + +int Evolvent::GetAllPreimages(const double *p, double xp[]) +{ + int preimNumber = 1; + TransformToStandardCube(p, p2.data()); + if(mMapType == Noninjective) + invmad(mTightness, xp, noninjectiveMaxPreimages, &preimNumber, p2.data(), mDimension, 4); + else + xyd(xp, mTightness, p2.data(), mDimension); + + return preimNumber; +} + +void xyd(double *xx, int m, double y[], int n) +{ + /* calculate preimage x for nearest level m center to y */ + /* (x - left boundary point of level m interval) */ + int n1, nexp, l, iu[10], iv[10]; + + double x, r1; + double r; + int iw[11]; + int i, j, it, is; + void numbr(int *, const int, const int, int&, int*, int*); + + n1 = n - 1; + for (nexp = 1, i = 0; i0) ? floor(dr) : ceil(dr); + + del = 1. / (mne - dr); + d1 = del*(incr + 0.5); + for (kx = -1; kx= n) { + xyd(&x, m, y, n); + dr = x*mne; + dd = dr - fmod(dr, 1.0); + //dd = (dr>0) ? floor(dr) : ceil(dr); + dr = dd / nexp; + dd = dd - dr + fmod(dr, 1.0); + //dd = dd - ((dr>0) ? floor(dr) : ceil(dr)); + x = dd*del; + if (kx>kp) break; + k = kx++; /* label 9 */ + if (kx == 0) xp[0] = x; + else { + while (k >= 0) { + dr = fabs(x - xp[k]); /* label 11 */ + if (dr <= d1) { + for (kx--; k= 0) && (u[i] = (u[i] <= 0.0) ? 1 : -1)<0; i--); + if (i<0)break; + } + *kxx = ++kx; + +} + +void mapd(double x, int m, double* y, int n, int key) +{ + /* mapping y(x) : 1 - center, 2 - line, 3 - node */ + // use key = 1 + + int n1, nexp, l, iq, iu[10], iv[10]; + double d, mne, dd, dr;//,tmp; + double p, r; + int iw[11]; + int it, is, i, j, k; + void node(int is, int n1, int nexp, int& l, int& iq, int iu[], int iv[]); + + p = 0.0; + n1 = n - 1; + for (nexp = 1, i = 0; i 2) { + dr = mne / nexp; + dr = dr - fmod(dr, 1.0); + //dr=(dr>0)?floor(dr):ceil(dr); + dd = mne - dr; + dr = d*dd; + dd = dr - fmod(dr, 1.0); + //dd=(dr>0)?floor(dr):ceil(dr); + dr = dd + (dd - 1) / (nexp - 1); + dd = dr - fmod(dr, 1.0); + //dd=(dr>0)?floor(dr):ceil(dr); + d = dd*(1. / (mne - 1.0)); + } + + for (j = 0; j0) || ((iq == 0) && (is == 0))) k = l; + else if (iq<0) k = (it == n1) ? 0 : n1; + r = r*0.5; + it = l; + for (i = 0; i= iff) { + if ((is == iff) && (is != 1)) { l = i; iq = -1; } + is = is - iff; + k2 = 1; + } + else { + k2 = -1; + if ((is == (iff - 1)) && (is != 0)) { l = i; iq = 1; } + } + j = -k1*k2; + iv[i] = j; + iu[i] = j; + k1 = k2; + } + iv[l] = iv[l] * iq; + iv[n1] = -iv[n1]; + } +} diff --git a/src/algs/ags/evolvent.hpp b/src/algs/ags/evolvent.hpp new file mode 100644 index 00000000..f85258f2 --- /dev/null +++ b/src/algs/ags/evolvent.hpp @@ -0,0 +1,48 @@ +/* +Copyright (C) 2018 Sovrasov V. - All Rights Reserved + * You may use, distribute and modify this code under the + * terms of the MIT license. + * You should have received a copy of the MIT license with + * this file. If not visit https://opensource.org/licenses/MIT +*/ +#pragma once + +#include + +namespace ags +{ + +enum MapType { + Simple = 1, Linear = 2, Noninjective = 3 +}; + +constexpr int noninjectiveMaxPreimages = 32; + +class Evolvent +{ +protected: + int mDimension; + int mTightness; + + std::vector mRho; + std::vector mShiftScalars; + std::vector p2; + + void TransformToStandardCube(const double *y, double *z); + void TransformToSearchDomain(const double *y, double *z); + + bool mIsInitialized; +private: + MapType mMapType; + int mMapKey; + +public: + Evolvent(); + Evolvent(int dimension, int tightness, const double* lb, const double* ub, MapType type = Simple); + ~Evolvent(); + + virtual void GetImage(double x, double y[]); + virtual int GetAllPreimages(const double* p, double xp[]); +}; + +} diff --git a/src/algs/ags/local_optimizer.cc b/src/algs/ags/local_optimizer.cc new file mode 100644 index 00000000..6c5f2538 --- /dev/null +++ b/src/algs/ags/local_optimizer.cc @@ -0,0 +1,135 @@ +/* +Copyright (C) 2018 Sovrasov V. - All Rights Reserved + * You may use, distribute and modify this code under the + * terms of the MIT license. + * You should have received a copy of the MIT license with + * this file. If not visit https://opensource.org/licenses/MIT +*/ +#include "local_optimizer.hpp" + +#include +#include +#include + +using namespace ags; + +#define MAX_LOCAL_ITERATIONS_NUMBER 20 + +void HookeJeevesOptimizer::SetParameters(double eps, double step, double stepMult) +{ + NLP_SOLVER_ASSERT(eps > 0 && step > 0 && stepMult > 0, "Wrong papameters of the local optimizer"); + mEps = eps; + mStep = step; + mStepMultiplier = stepMult; +} + +Trial HookeJeevesOptimizer::Optimize(std::shared_ptr> problem, + const Trial& startPoint, std::vector& trialsCounters) +{ + mProblem = problem; + mStartPoint = startPoint; + mTrialsCounters = std::vector(mProblem->GetConstraintsNumber() + 1, 0); + + int k = 0, i=0; + bool needRestart = true; + double currentFValue, nextFValue; + + while (i < MAX_LOCAL_ITERATIONS_NUMBER) { + i++; + if (needRestart) { + k = 0; + mCurrentPoint = mStartPoint; + mCurrentResearchDirection = mStartPoint; + currentFValue = ComputeObjective(mCurrentPoint.y); + needRestart = false; + } + + std::swap(mPreviousResearchDirection, mCurrentResearchDirection); + mCurrentResearchDirection = mCurrentPoint; + nextFValue = MakeResearch(mCurrentResearchDirection.y); + + if (currentFValue > nextFValue) { + DoStep(); + k++; + currentFValue = nextFValue; + } + else if (mStep > mEps) { + if (k != 0) + std::swap(mStartPoint, mPreviousResearchDirection); + else + mStep /= mStepMultiplier; + needRestart = true; + } + else + break; + } + + mPreviousResearchDirection.idx = 0; + while (mPreviousResearchDirection.idx < mProblem->GetConstraintsNumber()) + { + mTrialsCounters[mPreviousResearchDirection.idx]++; + mPreviousResearchDirection.g[mPreviousResearchDirection.idx] = + mProblem->Calculate(mPreviousResearchDirection.y, mPreviousResearchDirection.idx); + if (mPreviousResearchDirection.g[mPreviousResearchDirection.idx] > 0) + break; + mPreviousResearchDirection.idx++; + } + + if (mPreviousResearchDirection.idx == mProblem->GetConstraintsNumber()) + { + mPreviousResearchDirection.g[mPreviousResearchDirection.idx] = + mProblem->Calculate(mPreviousResearchDirection.y, mPreviousResearchDirection.idx); + mTrialsCounters[mPreviousResearchDirection.idx]++; + } + + for(size_t i = 0; i < mTrialsCounters.size(); i++) + trialsCounters[i] += mTrialsCounters[i]; + + return mPreviousResearchDirection; +} + +void HookeJeevesOptimizer::DoStep() +{ + for (int i = 0; i < mProblem->GetDimension(); i++) + mCurrentPoint.y[i] = (1 + mStepMultiplier)*mCurrentResearchDirection.y[i] - + mStepMultiplier*mPreviousResearchDirection.y[i]; +} + +double HookeJeevesOptimizer::ComputeObjective(const double* x) const +{ + for (int i = 0; i <= mProblem->GetConstraintsNumber(); i++) + { + double value = mProblem->Calculate(x, i); + mTrialsCounters[i]++; + if (i < mProblem->GetConstraintsNumber() && value > 0) + return std::numeric_limits::max(); + else if (i == mProblem->GetConstraintsNumber()) + return value; + } + return std::numeric_limits::max(); +} + +double HookeJeevesOptimizer::MakeResearch(double* startPoint) +{ + double bestValue = ComputeObjective(startPoint); + + for (int i = 0; i < mProblem->GetDimension(); i++) + { + startPoint[i] += mStep; + double rightFvalue = ComputeObjective(startPoint); + + if (rightFvalue > bestValue) + { + startPoint[i] -= 2 * mStep; + double leftFValue = ComputeObjective(startPoint); + if (leftFValue > bestValue) + startPoint[i] += mStep; + else + bestValue = leftFValue; + } + else + bestValue = rightFvalue; + } + + return bestValue; +} diff --git a/src/algs/ags/local_optimizer.hpp b/src/algs/ags/local_optimizer.hpp new file mode 100644 index 00000000..03c45acb --- /dev/null +++ b/src/algs/ags/local_optimizer.hpp @@ -0,0 +1,44 @@ +/* +Copyright (C) 2018 Sovrasov V. - All Rights Reserved + * You may use, distribute and modify this code under the + * terms of the MIT license. + * You should have received a copy of the MIT license with + * this file. If not visit https://opensource.org/licenses/MIT +*/ +#pragma once + +#include "data_types.hpp" + +#include +#include + +namespace ags +{ + +class HookeJeevesOptimizer +{ +private: + double mEps; + double mStep; + double mStepMultiplier; + + mutable std::vector mTrialsCounters; + + std::shared_ptr> mProblem; + + Trial mCurrentPoint; + Trial mStartPoint; + Trial mCurrentResearchDirection; + Trial mPreviousResearchDirection; + + void DoStep(); + double ComputeObjective(const double* x) const; + double MakeResearch(double*); + +public: + void SetParameters(double eps, double step, double stepMult); + Trial Optimize(std::shared_ptr> problem, + const Trial& startPoint, std::vector& trialsCounters); +}; + +} diff --git a/src/algs/ags/solver.cc b/src/algs/ags/solver.cc new file mode 100644 index 00000000..35b308a3 --- /dev/null +++ b/src/algs/ags/solver.cc @@ -0,0 +1,415 @@ +/* +Copyright (C) 2018 Sovrasov V. - All Rights Reserved + * You may use, distribute and modify this code under the + * terms of the MIT license. + * You should have received a copy of the MIT license with + * this file. If not visit https://opensource.org/licenses/MIT +*/ +#include "solver.hpp" + +#include +#include +#include +#include + +using namespace ags; + +namespace +{ + const double zeroHLevel = 1e-12; + + class ProblemInternal : public IGOProblem + { + private: + std::vector mFunctions; + std::vector mLeftBound; + std::vector mRightBound; + + unsigned mDimension; + unsigned mConstraintsNumber; + + public: + ProblemInternal(const std::vector& functions, + const std::vector& leftBound, const std::vector& rightBound) + { + mFunctions = functions; + mConstraintsNumber = mFunctions.size() - 1; + mDimension = leftBound.size(); + mLeftBound = leftBound; + mRightBound = rightBound; + } + + double Calculate(const double* y, int fNumber) const + { + return mFunctions[fNumber](y); + } + int GetConstraintsNumber() const + { + return mConstraintsNumber; + } + int GetDimension() const + { + return mDimension; + } + void GetBounds(double* left, double* right) const + { + for(size_t i = 0; i < mDimension; i++) + { + left[i] = mLeftBound[i]; + right[i] = mRightBound[i]; + } + } + int GetOptimumPoint(double* y) const {return 0;} + double GetOptimumValue() const {return 0;} + }; +} + +NLPSolver::NLPSolver() {} + +void NLPSolver::SetParameters(const SolverParameters& params) +{ + mParameters = params; +} + +void NLPSolver::SetProblem(std::shared_ptr> problem) +{ + mProblem = problem; + NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= solverMaxConstraints, + "Current implementation supports up to " + std::to_string(solverMaxConstraints) + + " nonlinear inequality constraints"); + InitLocalOptimizer(); +} + +void NLPSolver::SetProblem(const std::vector& functions, + const std::vector& leftBound, const std::vector& rightBound) +{ + NLP_SOLVER_ASSERT(leftBound.size() == rightBound.size(), "Inconsistent dimensions of bounds"); + NLP_SOLVER_ASSERT(leftBound.size() > 0, "Zero problem dimension"); + mProblem = std::make_shared(functions, leftBound, rightBound); + NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= solverMaxConstraints, + "Current implementation supports up to " + std::to_string(solverMaxConstraints) + + " nonlinear inequality constraints"); + InitLocalOptimizer(); +} + +std::vector NLPSolver::GetCalculationsStatistics() const +{ + return mCalculationsCounters; +} + +std::vector NLPSolver::GetHolderConstantsEstimations() const +{ + return mHEstimations; +} + +void NLPSolver::InitLocalOptimizer() +{ + std::vector leftBound(mProblem->GetDimension()); + std::vector rightBound(mProblem->GetDimension()); + mProblem->GetBounds(leftBound.data(), rightBound.data()); + + double maxSize = 0; + for(size_t i = 0; i < leftBound.size(); i++) + maxSize = std::max(rightBound[i] - leftBound[i], maxSize); + + NLP_SOLVER_ASSERT(maxSize > 0, "Empty search domain"); + + mLocalOptimizer.SetParameters(maxSize / 1000, maxSize / 100, 2); +} + +void NLPSolver::InitDataStructures() +{ + double leftDomainBound[solverMaxDim], rightDomainBound[solverMaxDim]; + mProblem->GetBounds(leftDomainBound, rightDomainBound); + mEvolvent = Evolvent(mProblem->GetDimension(), mParameters.evolventDensity, + leftDomainBound, rightDomainBound); + + mNextPoints.resize(mParameters.numPoints); + mOptimumEstimation.idx = -1; + + mZEstimations.resize(mProblem->GetConstraintsNumber() + 1); + std::fill(mZEstimations.begin(), mZEstimations.end(), + std::numeric_limits::max()); + mNextIntervals.resize(mParameters.numPoints); + mHEstimations.resize(mProblem->GetConstraintsNumber() + 1); + std::fill(mHEstimations.begin(), mHEstimations.end(), 1.0); + mCalculationsCounters.resize(mProblem->GetConstraintsNumber() + 1); + std::fill(mCalculationsCounters.begin(), mCalculationsCounters.end(), 0); + mQueue = PriorityQueue(); + mIterationsCounter = 0; + mMinDelta = std::numeric_limits::max(); + mMaxIdx = -1; +} + +void NLPSolver::ClearDataStructures() +{ + for (const auto& ptr : mSearchInformation) + delete ptr; + mSearchInformation.clear(); + mQueue = PriorityQueue(); +} + +Trial NLPSolver::Solve() +{ + bool needStop = false; + InitDataStructures(); + FirstIteration(); + + do { + InsertIntervals(); + EstimateOptimum(); + if (mNeedRefillQueue || mQueue.size() < mParameters.numPoints) + RefillQueue(); + CalculateNextPoints(); + MakeTrials(); + needStop = mMinDelta < mParameters.eps; + mIterationsCounter++; + } while(mIterationsCounter < mParameters.itersLimit && !needStop); + + ClearDataStructures(); + + if (mParameters.refineSolution && mOptimumEstimation.idx == mProblem->GetConstraintsNumber()) { + auto localTrial = mLocalOptimizer.Optimize(mProblem, mOptimumEstimation, mCalculationsCounters); + int idx = mOptimumEstimation.idx; + if (localTrial.idx == idx && localTrial.g[idx] < mOptimumEstimation.g[idx]) + mOptimumEstimation = localTrial; + } + + return mOptimumEstimation; +} + +void NLPSolver::FirstIteration() +{ + Trial leftBound = Trial(0); + leftBound.idx = -1; + Trial rightBound = Trial(1.); + rightBound.idx = -1; + + for (size_t i = 1; i <= mParameters.numPoints; i++) + { + mNextPoints[i - 1] = Trial((double)i / (mParameters.numPoints + 1)); + mEvolvent.GetImage(mNextPoints[i - 1].x, mNextPoints[i - 1].y); + } + + MakeTrials(); + EstimateOptimum(); + + for (size_t i = 0; i <= mParameters.numPoints; i++) + { + Interval* pNewInterval; + if (i == 0) + pNewInterval = new Interval(leftBound, mNextPoints[i]); + else if (i == mParameters.numPoints) + pNewInterval = new Interval(mNextPoints[i - 1], rightBound); + else + pNewInterval = new Interval(mNextPoints[i - 1], mNextPoints[i]); + pNewInterval->delta = pow(pNewInterval->pr.x - pNewInterval->pl.x, + 1. / mProblem->GetDimension()); + mMinDelta = std::min(mMinDelta, pNewInterval->delta); + auto insRes = mSearchInformation.insert(pNewInterval); + UpdateAllH(insRes.first); + } + RefillQueue(); + CalculateNextPoints(); + MakeTrials(); +} + +void NLPSolver::MakeTrials() +{ + for (size_t i = 0; i < mNextPoints.size(); i++) + { + int idx = 0; + while(idx < mProblem->GetConstraintsNumber()) + { + mNextPoints[i].idx = idx; + double val = mProblem->Calculate(mNextPoints[i].y, idx); + mCalculationsCounters[idx]++; + mNextPoints[i].g[idx] = val; + if (val > 0) + break; + idx++; + } + + if(idx > mMaxIdx) + { + mMaxIdx = idx; + for(int i = 0; i < mMaxIdx; i++) + mZEstimations[i] = -mParameters.epsR*mHEstimations[i]; + mNeedRefillQueue = true; + } + if (idx == mProblem->GetConstraintsNumber()) + { + mCalculationsCounters[idx]++; + mNextPoints[i].idx = idx; + mNextPoints[i].g[idx] = mProblem->Calculate(mNextPoints[i].y, idx); + } + if(mNextPoints[i].idx == mMaxIdx && + mNextPoints[i].g[mMaxIdx] < mZEstimations[mMaxIdx]) + { + mZEstimations[mMaxIdx] = mNextPoints[i].g[mMaxIdx]; + mNeedRefillQueue = true; + } + } +} + +void NLPSolver::InsertIntervals() +{ + for (size_t i = 0; i < mParameters.numPoints; i++) + { + Interval* pOldInterval = mNextIntervals[i]; + Interval* pNewInterval = new Interval(mNextPoints[i], pOldInterval->pr); + pOldInterval->pr = mNextPoints[i]; + pOldInterval->delta = pow(pOldInterval->pr.x - pOldInterval->pl.x, + 1. / mProblem->GetDimension()); + pNewInterval->delta = pow(pNewInterval->pr.x - pNewInterval->pl.x, + 1. / mProblem->GetDimension()); + mMinDelta = std::min(mMinDelta, pNewInterval->delta); + mMinDelta = std::min(mMinDelta, pOldInterval->delta); + + auto insResult = mSearchInformation.insert(pNewInterval); + bool wasInserted = insResult.second; + if(!wasInserted) + throw std::runtime_error("Error during interval insertion."); + + UpdateAllH(insResult.first); + UpdateAllH(--insResult.first); + + if(!mNeedRefillQueue) + { + pNewInterval->R = CalculateR(pNewInterval); + mNextIntervals[i]->R = CalculateR(mNextIntervals[i]); + mQueue.push(pNewInterval); + mQueue.push(pOldInterval); + } + } +} + +void NLPSolver::CalculateNextPoints() +{ + for(size_t i = 0; i < mParameters.numPoints; i++) + { + mNextIntervals[i] = mQueue.top(); + mQueue.pop(); + mNextPoints[i].x = GetNextPointCoordinate(mNextIntervals[i]); + + if (mNextPoints[i].x >= mNextIntervals[i]->pr.x || mNextPoints[i].x <= mNextIntervals[i]->pl.x) + throw std::runtime_error("The next point is outside of the subdivided interval"); + + mEvolvent.GetImage(mNextPoints[i].x, mNextPoints[i].y); + } +} + +void NLPSolver::RefillQueue() +{ + mQueue = PriorityQueue(); + for (const auto& pInterval : mSearchInformation) + { + pInterval->R = CalculateR(pInterval); + mQueue.push(pInterval); + } + mNeedRefillQueue = false; +} + +void NLPSolver::EstimateOptimum() +{ + for (size_t i = 0; i < mNextPoints.size(); i++) + { + if (mOptimumEstimation.idx < mNextPoints[i].idx || + mOptimumEstimation.idx == mNextPoints[i].idx && + mOptimumEstimation.g[mOptimumEstimation.idx] > mNextPoints[i].g[mNextPoints[i].idx]) + { + mOptimumEstimation = mNextPoints[i]; + mNeedRefillQueue = true; + } + } +} + +void NLPSolver::UpdateH(double newValue, int index) +{ + if (newValue > mHEstimations[index] || mHEstimations[index] == 1.0 && newValue > zeroHLevel) + { + mHEstimations[index] = newValue; + mNeedRefillQueue = true; + } +} + +void NLPSolver::UpdateAllH(std::set::iterator iterator) +{ + Interval* pInterval = *iterator; + if (pInterval->pl.idx < 0) + return; + + if (pInterval->pl.idx == pInterval->pr.idx) + UpdateH(fabs(pInterval->pr.g[pInterval->pr.idx] - pInterval->pl.g[pInterval->pl.idx]) / + pInterval->delta, pInterval->pl.idx); + else + { + auto rightIterator = iterator; + auto leftIterator = iterator; + //right lookup + ++rightIterator; + while(rightIterator != mSearchInformation.end() && (*rightIterator)->pl.idx < pInterval->pl.idx) + ++rightIterator; + if (rightIterator != mSearchInformation.end() && (*rightIterator)->pl.idx >= pInterval->pl.idx) + { + int idx = pInterval->pl.idx; + UpdateH(fabs((*rightIterator)->pl.g[idx] - pInterval->pl.g[idx]) / + pow((*rightIterator)->pl.x - pInterval->pl.x, 1. / mProblem->GetDimension()), idx); + } + + //left lookup + --leftIterator; + while(leftIterator != mSearchInformation.begin() && (*leftIterator)->pl.idx < pInterval->pl.idx) + --leftIterator; + if (leftIterator != mSearchInformation.begin() && (*leftIterator)->pl.idx >= pInterval->pl.idx) + { + int idx = pInterval->pl.idx; + UpdateH(fabs((*leftIterator)->pl.g[idx] - pInterval->pl.g[idx]) / + pow(pInterval->pl.x - (*leftIterator)->pl.x, 1. / mProblem->GetDimension()), idx); + } + } +} + +double NLPSolver::CalculateR(Interval* i) const +{ + if(i->pl.idx == i->pr.idx) + { + const int v = i->pr.idx; + return i->delta + pow((i->pr.g[v] - i->pl.g[v]) / (mParameters.r * mHEstimations[v]), 2) / i->delta - + 2.*(i->pr.g[v] + i->pl.g[v] - 2*mZEstimations[v]) / (mParameters.r * mHEstimations[v]); + } + else if(i->pl.idx < i->pr.idx) + return 2*i->delta - 4*(i->pr.g[i->pr.idx] - mZEstimations[i->pr.idx]) / (mParameters.r * mHEstimations[i->pr.idx]); + else + return 2*i->delta - 4*(i->pl.g[i->pl.idx] - mZEstimations[i->pl.idx]) / (mParameters.r * mHEstimations[i->pl.idx]); +} + +double NLPSolver::GetNextPointCoordinate(Interval* i) const +{ + double x; + if(i->pr.idx == i->pl.idx) + { + const int v = i->pr.idx; + double dg = i->pr.g[v] - i->pl.g[v]; + x = 0.5 * (i->pr.x + i->pl.x) - + 0.5*((dg > 0.) ? 1. : -1.) * pow(fabs(dg) / mHEstimations[v], mProblem->GetDimension()) / mParameters.r; + } + else + x = 0.5 * (i->pr.x + i->pl.x); + + if (x >= i->pr.x || x <= i->pl.x) + throw std::runtime_error("Point is outside of interval\n"); + + return x; +} + +bool solver_utils::checkVectorsDiff(const double* y1, const double* y2, size_t dim, double eps) +{ + for (size_t i = 0; i < dim; i++) + { + if (fabs(y1[i] - y2[i]) > eps) + return true; + } + + return false; +} diff --git a/src/algs/ags/solver.hpp b/src/algs/ags/solver.hpp new file mode 100644 index 00000000..75ac4196 --- /dev/null +++ b/src/algs/ags/solver.hpp @@ -0,0 +1,102 @@ +/* +Copyright (C) 2018 Sovrasov V. - All Rights Reserved + * You may use, distribute and modify this code under the + * terms of the MIT license. + * You should have received a copy of the MIT license with + * this file. If not visit https://opensource.org/licenses/MIT +*/ +#pragma once + +#include "data_types.hpp" +#include "evolvent.hpp" +#include "local_optimizer.hpp" + +#include +#include +#include +#include + +namespace ags +{ + +struct SolverParameters +{ + double eps = 0.01; //method tolerance. Less value -- better search precision, less probability of early stop. + double r = 3; //reliability parameter. Higher value of r -- slower convergence, higher chance to cache the global minima. + unsigned numPoints = 1; //number of new points per iteration. > 1 is useless in current implementation. + unsigned itersLimit = 20000; // max number of iterations. + unsigned evolventDensity = 12; // density of evolvent. By default density is 2^-12 on hybercube [0,1]^N, + // which means that maximum search accuracyis 2^-12. If search hypercube is large the density can be increased accordingly to achieve better accuracy. + double epsR = 0.001; // parameter which prevents method from paying too much attention to constraints. Greater values of this parameter speed up convergence, + // but global minima can be lost. + bool refineSolution = false; //if true, the fibal solution will be refined with the HookJeves method. + + SolverParameters() {} + SolverParameters(double _eps, double _r, + double epsR_, unsigned _trialsLimit) : + eps(_eps), r(_r), itersLimit(_trialsLimit), epsR(epsR_) + {} +}; + +class NLPSolver +{ +protected: + using PriorityQueue = + std::priority_queue, CompareByR>; + + HookeJeevesOptimizer mLocalOptimizer; + + SolverParameters mParameters; + std::shared_ptr> mProblem; + Evolvent mEvolvent; + + std::vector mHEstimations; + std::vector mZEstimations; + std::vector mNextPoints; + PriorityQueue mQueue; + std::set mSearchInformation; + std::vector mNextIntervals; + Trial mOptimumEstimation; + + std::vector mCalculationsCounters; + unsigned mIterationsCounter; + bool mNeedRefillQueue; + double mMinDelta; + int mMaxIdx; + + void InitLocalOptimizer(); + void FirstIteration(); + void MakeTrials(); + void InsertIntervals(); + void CalculateNextPoints(); + void RefillQueue(); + void EstimateOptimum(); + + void InitDataStructures(); + void ClearDataStructures(); + + void UpdateAllH(std::set::iterator); + void UpdateH(double newValue, int index); + double CalculateR(Interval*) const; + double GetNextPointCoordinate(Interval*) const; + +public: + using FuncPtr = double(*)(const double*); + NLPSolver(); + + void SetParameters(const SolverParameters& params); + void SetProblem(std::shared_ptr> problem); + void SetProblem(const std::vector& functions, + const std::vector& leftBound, const std::vector& rightBound); + + Trial Solve(); + std::vector GetCalculationsStatistics() const; + std::vector GetHolderConstantsEstimations() const; +}; + +namespace solver_utils +{ + bool checkVectorsDiff(const double* y1, const double* y2, size_t dim, double eps); +} + +} diff --git a/src/api/general.c b/src/api/general.c index 410771b5..5b536ee4 100644 --- a/src/api/general.c +++ b/src/api/general.c @@ -7,17 +7,17 @@ * distribute, sublicense, and/or sell copies of the Software, and to * permit persons to whom the Software is furnished to do so, subject to * the following conditions: - * + * * The above copyright notice and this permission notice shall be * included in all copies or substantial portions of the Software. - * + * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include "nlopt-internal.h" @@ -81,7 +81,12 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { "Multi-level single-linkage (MLSL), quasi-random (global, needs sub-algorithm)", "Sequential Quadratic Programming (SQP) (local, derivative)", "CCSA (Conservative Convex Separable Approximations) with simple quadratic approximations (local, derivative)", - "ESCH evolutionary strategy" + "ESCH evolutionary strategy", +#ifdef NLOPT_CXX + "AGS (global, no-derivative)" +#else + "AGS (NOT COMPILED)" +#endif }; const char * NLOPT_STDCALL nlopt_algorithm_name(nlopt_algorithm a) diff --git a/src/api/nlopt.h b/src/api/nlopt.h index 328b6d09..1805ab01 100644 --- a/src/api/nlopt.h +++ b/src/api/nlopt.h @@ -7,17 +7,17 @@ * distribute, sublicense, and/or sell copies of the Software, and to * permit persons to whom the Software is furnished to do so, subject to * the following conditions: - * + * * The above copyright notice and this permission notice shall be * included in all copies or substantial portions of the Software. - * + * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #ifndef NLOPT_H @@ -67,7 +67,7 @@ typedef void (*nlopt_mfunc)(unsigned m, double *result, double *gradient, /* NULL if not needed */ void *func_data); -/* A preconditioner, which preconditions v at x to return vpre. +/* A preconditioner, which preconditions v at x to return vpre. (The meaning of "preconditioning" is algorithm-dependent.) */ typedef void (*nlopt_precond)(unsigned n, const double *x, const double *v, double *vpre, void *data); @@ -75,10 +75,10 @@ typedef void (*nlopt_precond)(unsigned n, const double *x, const double *v, typedef enum { /* Naming conventions: - NLOPT_{G/L}{D/N}_* - = global/local derivative/no-derivative optimization, - respectively - + NLOPT_{G/L}{D/N}_* + = global/local derivative/no-derivative optimization, + respectively + *_RAND algorithms involve some randomization. *_NOSCAL algorithms are *not* scaled to a unit hypercube @@ -151,6 +151,8 @@ typedef enum { NLOPT_GN_ESCH, + NLOPT_AGS, + NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; @@ -198,9 +200,9 @@ NLOPT_EXTERN(nlopt_opt) nlopt_copy(const nlopt_opt opt); NLOPT_EXTERN(nlopt_result) nlopt_optimize(nlopt_opt opt, double *x, double *opt_f); -NLOPT_EXTERN(nlopt_result) nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, +NLOPT_EXTERN(nlopt_result) nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, void *f_data); -NLOPT_EXTERN(nlopt_result) nlopt_set_max_objective(nlopt_opt opt, nlopt_func f, +NLOPT_EXTERN(nlopt_result) nlopt_set_max_objective(nlopt_opt opt, nlopt_func f, void *f_data); NLOPT_EXTERN(nlopt_result) nlopt_set_precond_min_objective(nlopt_opt opt, nlopt_func f, nlopt_precond pre, void *f_data); @@ -213,12 +215,12 @@ NLOPT_EXTERN(const char*) nlopt_get_errmsg(nlopt_opt opt); /* constraints: */ -NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bounds(nlopt_opt opt, +NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bounds(nlopt_opt opt, const double *lb); NLOPT_EXTERN(nlopt_result) nlopt_set_lower_bounds1(nlopt_opt opt, double lb); -NLOPT_EXTERN(nlopt_result) nlopt_get_lower_bounds(const nlopt_opt opt, +NLOPT_EXTERN(nlopt_result) nlopt_get_lower_bounds(const nlopt_opt opt, double *lb); -NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bounds(nlopt_opt opt, +NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bounds(nlopt_opt opt, const double *ub); NLOPT_EXTERN(nlopt_result) nlopt_set_upper_bounds1(nlopt_opt opt, double ub); NLOPT_EXTERN(nlopt_result) nlopt_get_upper_bounds(const nlopt_opt opt, @@ -283,7 +285,7 @@ NLOPT_EXTERN(int) nlopt_get_force_stop(const nlopt_opt opt); /* more algorithm-specific parameters */ -NLOPT_EXTERN(nlopt_result) nlopt_set_local_optimizer(nlopt_opt opt, +NLOPT_EXTERN(nlopt_result) nlopt_set_local_optimizer(nlopt_opt opt, const nlopt_opt local_opt); NLOPT_EXTERN(nlopt_result) nlopt_set_population(nlopt_opt opt, unsigned pop); @@ -292,12 +294,12 @@ NLOPT_EXTERN(unsigned) nlopt_get_population(const nlopt_opt opt); NLOPT_EXTERN(nlopt_result) nlopt_set_vector_storage(nlopt_opt opt, unsigned dim); NLOPT_EXTERN(unsigned) nlopt_get_vector_storage(const nlopt_opt opt); -NLOPT_EXTERN(nlopt_result) nlopt_set_default_initial_step(nlopt_opt opt, +NLOPT_EXTERN(nlopt_result) nlopt_set_default_initial_step(nlopt_opt opt, const double *x); -NLOPT_EXTERN(nlopt_result) nlopt_set_initial_step(nlopt_opt opt, +NLOPT_EXTERN(nlopt_result) nlopt_set_initial_step(nlopt_opt opt, const double *dx); NLOPT_EXTERN(nlopt_result) nlopt_set_initial_step1(nlopt_opt opt, double dx); -NLOPT_EXTERN(nlopt_result) nlopt_get_initial_step(const nlopt_opt opt, +NLOPT_EXTERN(nlopt_result) nlopt_get_initial_step(const nlopt_opt opt, const double *x, double *dx); /* the following are functions mainly designed to be used internally @@ -323,7 +325,7 @@ NLOPT_EXTERN(void) nlopt_munge_data(nlopt_opt opt, #if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__==3 && __GNUC_MINOR__ > 0)) # define NLOPT_DEPRECATED __attribute__((deprecated)) #else -# define NLOPT_DEPRECATED +# define NLOPT_DEPRECATED #endif typedef double (*nlopt_func_old)(int n, const double *x, diff --git a/src/api/optimize.c b/src/api/optimize.c index 9592653b..ea328df6 100644 --- a/src/api/optimize.c +++ b/src/api/optimize.c @@ -33,6 +33,7 @@ #ifdef NLOPT_CXX # include "stogo.h" +# include "ags.h" #endif #include "cdirect.h" @@ -347,6 +348,7 @@ static int elimdim_wrapcheck(nlopt_opt opt) case NLOPT_LN_SBPLX: case NLOPT_GN_ISRES: case NLOPT_GN_ESCH: + case NLOPT_AGS: case NLOPT_GD_STOGO: case NLOPT_GD_STOGO_RAND: return 1; @@ -504,6 +506,19 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) break; } + case NLOPT_AGS: +#ifdef NLOPT_CXX + if (!finite_domain(n, lb, ub)) + RETURN_ERR(NLOPT_INVALID_ARGS, opt, + "finite domain required for global algorithm"); + //if (!ags_minimize(ni, f, f_data, x, minf, lb, ub, &stop, 0)) + if (!ags_minimize()) + return NLOPT_FAILURE; + break; +#else + return NLOPT_INVALID_ARGS; +#endif + case NLOPT_GD_STOGO: case NLOPT_GD_STOGO_RAND: #ifdef NLOPT_CXX From de8fa3c5f2e2ffe0e020d299e2cc7cc61e2f38ec Mon Sep 17 00:00:00 2001 From: sovrasov Date: Sun, 15 Jul 2018 14:00:02 +0300 Subject: [PATCH 04/23] Finish basic integration of AGS --- src/algs/ags/ags.cc | 92 ++++++++++++++++++++++++++++++++++------- src/algs/ags/ags.h | 21 ++++------ src/algs/ags/solver.cc | 4 +- src/algs/ags/solver.hpp | 3 +- src/api/optimize.c | 4 +- 5 files changed, 89 insertions(+), 35 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index e3db8b10..59d40cd6 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -3,20 +3,84 @@ #include "ags.h" #include "solver.hpp" +#include -/* -int ags_minimize(int n, - objective_func fgrad, void *data, - double *x, double *minf, - const double *l, const double *u, -#ifdef NLOPT_UTIL_H - nlopt_stopping *stop, -#else - long int maxeval, double maxtime, -#endif - int nrandom) - */ -int ags_minimize() +double ags_r = 3; +double eps_res = 0.001; +unsigned evolvent_density = 11; +int ags_verbose = 0; + +int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, + double *x, double *minf, const double *l, const double *u, nlopt_stopping *stop) { - return 1; + int ret_code = NLOPT_SUCCESS; + + if (n > ags::solverMaxDim) + return NLOPT_INVALID_ARGS; + if(m != nlopt_count_constraints(m, fc) || m > ags::solverMaxConstraints) + return NLOPT_INVALID_ARGS; + + std::vector lb(l, l + n); + std::vector ub(u, u + n); + std::vector functions; + for (unsigned i = 0; i < m; i++) + { + if (fc[i].m != 1) + return NLOPT_INVALID_ARGS; + functions.push_back([fc, data, n, i](const double* x) {return fc[i].f(n, x, NULL, data);}); + } + functions.push_back([func, data, n](const double* x) {return func(n, x, NULL, data);}); + + ags::SolverParameters params; + params.r = ags_r; + params.itersLimit = stop->maxeval; + params.eps = 0; + params.evolventDensity = evolvent_density; + params.epsR = eps_res; + + ags::NLPSolver solver; + solver.SetParameters(params); + solver.SetProblem(functions, lb, ub); + + ags::Trial optPoint; + try + { + optPoint = solver.Solve(); + } + catch (const std::runtime_error& exp) + { + std::cerr << "AGS internal error: " << std::string(exp.what()) << std::endl; + return NLOPT_FAILURE; + } + + if (ags_verbose) + { + auto calcCounters = solver.GetCalculationsStatistics(); + auto holderConstEstimations = solver.GetHolderConstantsEstimations(); + + std::cout << std::string(20, '-') << "AGS statistics: " << std::string(20, '-') << std::endl; + for (size_t i = 0; i < calcCounters.size() - 1; i++) + std::cout << "Number of calculations of constraint # " << i << ": " << calcCounters[i] << "\n"; + std::cout << "Number of calculations of objective: " << calcCounters.back() << "\n";; + + for (size_t i = 0; i < holderConstEstimations.size() - 1; i++) + std::cout << "Estimation of Holder constant of function # " << i << ": " << holderConstEstimations[i] << "\n"; + std::cout << "Estimation of Holder constant of objective: " << holderConstEstimations.back() << "\n"; + if (optPoint.idx != m) + std::cout << "Feasible point not found" << "\n"; + std::cout << std::string(40, '-') << std::endl; + } + + if (m == optPoint.idx) + { + memcpy(x, optPoint.y, n*sizeof(x[0])); + *minf = optPoint.g[optPoint.idx]; + } + else //feasible point not found. + return NLOPT_FAILURE; + + if (solver.GetCalculationsStatistics()[0] >= params.itersLimit) + return NLOPT_MAXEVAL_REACHED; + + return ret_code; } diff --git a/src/algs/ags/ags.h b/src/algs/ags/ags.h index 0d03cfd6..9fc7783f 100644 --- a/src/algs/ags/ags.h +++ b/src/algs/ags/ags.h @@ -50,20 +50,13 @@ extern "C" { */ -int ags_minimize(); -/* -int ags_minimize(int n, - objective_func fgrad, void *data, - double *x, double *minf, - const double *l, const double *u, -#ifdef NLOPT_UTIL_H - nlopt_stopping *stop, -#else - long int maxeval, double maxtime, -#endif - int nrandom); -*/ -extern int ags_verbose; /* set to nonzero for verbose output */ +int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, + double *x, double *minf, const double *l, const double *u, nlopt_stopping *stop); + +extern double ags_r; +extern double eps_res; +extern unsigned evolvent_density; +extern int ags_verbose; #ifdef __cplusplus } diff --git a/src/algs/ags/solver.cc b/src/algs/ags/solver.cc index 35b308a3..83c09593 100644 --- a/src/algs/ags/solver.cc +++ b/src/algs/ags/solver.cc @@ -212,6 +212,7 @@ void NLPSolver::FirstIteration() RefillQueue(); CalculateNextPoints(); MakeTrials(); + mIterationsCounter += 2; } void NLPSolver::MakeTrials() @@ -397,9 +398,6 @@ double NLPSolver::GetNextPointCoordinate(Interval* i) const else x = 0.5 * (i->pr.x + i->pl.x); - if (x >= i->pr.x || x <= i->pl.x) - throw std::runtime_error("Point is outside of interval\n"); - return x; } diff --git a/src/algs/ags/solver.hpp b/src/algs/ags/solver.hpp index 75ac4196..6c998ef6 100644 --- a/src/algs/ags/solver.hpp +++ b/src/algs/ags/solver.hpp @@ -15,6 +15,7 @@ Copyright (C) 2018 Sovrasov V. - All Rights Reserved #include #include #include +#include namespace ags { @@ -81,7 +82,7 @@ class NLPSolver double GetNextPointCoordinate(Interval*) const; public: - using FuncPtr = double(*)(const double*); + using FuncPtr = std::function; NLPSolver(); void SetParameters(const SolverParameters& params); diff --git a/src/api/optimize.c b/src/api/optimize.c index ea328df6..03f45952 100644 --- a/src/api/optimize.c +++ b/src/api/optimize.c @@ -511,9 +511,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) if (!finite_domain(n, lb, ub)) RETURN_ERR(NLOPT_INVALID_ARGS, opt, "finite domain required for global algorithm"); - //if (!ags_minimize(ni, f, f_data, x, minf, lb, ub, &stop, 0)) - if (!ags_minimize()) - return NLOPT_FAILURE; + return ags_minimize(ni, f, f_data, opt->m, opt->fc, x, minf, lb, ub, &stop); break; #else return NLOPT_INVALID_ARGS; From e050013646f267612fba9d3051dd62311623f993 Mon Sep 17 00:00:00 2001 From: sovrasov Date: Sun, 15 Jul 2018 14:41:15 +0300 Subject: [PATCH 05/23] Clenup ags header, change cmake for test --- src/algs/ags/ags.cc | 4 ++-- src/algs/ags/ags.h | 50 +++++------------------------------------ src/algs/ags/solver.hpp | 2 +- src/api/nlopt.h | 2 +- test/CMakeLists.txt | 4 ++-- 5 files changed, 12 insertions(+), 50 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index 59d40cd6..a9ec9c1e 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -7,7 +7,7 @@ double ags_r = 3; double eps_res = 0.001; -unsigned evolvent_density = 11; +unsigned evolvent_density = 12; int ags_verbose = 0; int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, @@ -34,7 +34,7 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons ags::SolverParameters params; params.r = ags_r; params.itersLimit = stop->maxeval; - params.eps = 0; + params.eps = 1e-64; params.evolventDensity = evolvent_density; params.epsR = eps_res; diff --git a/src/algs/ags/ags.h b/src/algs/ags/ags.h index 9fc7783f..72221903 100644 --- a/src/algs/ags/ags.h +++ b/src/algs/ags/ags.h @@ -10,53 +10,15 @@ extern "C" { #endif -/* search for the global minimum of the function fgrad(n, x, grad, data) - inside a simple n-dimensional hyperrectangle. - - Input: - - n: dimension of search space (number of decision variables) - - fgrad: the objective function of the form fgrad(n, x, grad, data), - returning the objective function at x, where - n: dimension of search space - x: pointer to array of length n, point to evaluate - grad: if non-NULL, an array of length n which - should on return be the gradient d(fgrad)/dx - [ if NULL, grad should be ignored ] - data: arbitrary data pointer, whose value is the - data argument of ags_minimize - - data: arbitrary pointer to any auxiliary data needed by fgrad - - l, u: arrays of length n giving the lower and upper bounds of the - search space - - maxeval: if nonzero, a maximum number of fgrad evaluations - maxtime: if nonzero, a maximum time (in seconds) - -- REPLACED in NLopt by nlopt_stopping *stop - - nrandom: number of randomized search points to use per box, - in addition to 2*n+1 deterministic search points - (0 for a deterministic algorithm). - - Output: - - minf: the minimum value of the objective function found - - x: pointer to array of length n, giving the location of the minimum - - Return value: 0 if no minimum found, 1 otherwise. - - */ - int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, double *x, double *minf, const double *l, const double *u, nlopt_stopping *stop); -extern double ags_r; -extern double eps_res; -extern unsigned evolvent_density; -extern int ags_verbose; +extern double ags_r; //reliability parameter. Higher value of r -- slower convergence, higher chance to cache the global minima. +extern double eps_res; // parameter which prevents method from paying too much attention to constraints. Greater values of this parameter speed up convergence, +// but global minima can be lost. +extern unsigned evolvent_density; // density of evolvent. By default density is 2^-12 on hybercube [0,1]^N, +// which means that maximum search accuracyis 2^-12. If search hypercube is large the density can be increased accordingly to achieve better accuracy. +extern int ags_verbose; //print additional info #ifdef __cplusplus } diff --git a/src/algs/ags/solver.hpp b/src/algs/ags/solver.hpp index 6c998ef6..c1af38b4 100644 --- a/src/algs/ags/solver.hpp +++ b/src/algs/ags/solver.hpp @@ -22,7 +22,7 @@ namespace ags struct SolverParameters { - double eps = 0.01; //method tolerance. Less value -- better search precision, less probability of early stop. + double eps = 0.01; //method tolerance in Holder metric on 1d interval. Less value -- better search precision, less probability of early stop. double r = 3; //reliability parameter. Higher value of r -- slower convergence, higher chance to cache the global minima. unsigned numPoints = 1; //number of new points per iteration. > 1 is useless in current implementation. unsigned itersLimit = 20000; // max number of iterations. diff --git a/src/api/nlopt.h b/src/api/nlopt.h index 1805ab01..5ced305f 100644 --- a/src/api/nlopt.h +++ b/src/api/nlopt.h @@ -151,7 +151,7 @@ typedef enum { NLOPT_GN_ESCH, - NLOPT_AGS, + NLOPT_AGS, NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f1fff664..3d4f0d09 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -7,12 +7,12 @@ target_link_libraries (testopt ${nlopt_lib}) target_include_directories (testopt PRIVATE ${NLOPT_PRIVATE_INCLUDE_DIRS}) add_dependencies (tests testopt) -foreach (algo_index RANGE 29)# 42 +foreach (algo_index RANGE 29)# 43 foreach (obj_index RANGE 1)# 21 set (enable_ TRUE) # cxx stogo if (NOT NLOPT_CXX) - if (algo_index STREQUAL 8 OR algo_index STREQUAL 9) + if (algo_index STREQUAL 8 OR algo_index STREQUAL 9 OR algo_index STREQUAL 43) set (enable_ FALSE) endif () endif () From 765785132d2de262123570c28f200773d8d794ac Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Mon, 16 Jul 2018 18:24:13 +0300 Subject: [PATCH 06/23] AGS: add stop by reaching required value --- src/algs/ags/ags.cc | 2 ++ src/algs/ags/solver.cc | 10 ++++++---- src/algs/ags/solver.hpp | 3 +++ 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index a9ec9c1e..a0a07757 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -4,6 +4,7 @@ #include "ags.h" #include "solver.hpp" #include +#include double ags_r = 3; double eps_res = 0.001; @@ -37,6 +38,7 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons params.eps = 1e-64; params.evolventDensity = evolvent_density; params.epsR = eps_res; + params.stopVal = stop->minf_max; ags::NLPSolver solver; solver.SetParameters(params); diff --git a/src/algs/ags/solver.cc b/src/algs/ags/solver.cc index 83c09593..5e7935e0 100644 --- a/src/algs/ags/solver.cc +++ b/src/algs/ags/solver.cc @@ -8,7 +8,6 @@ Copyright (C) 2018 Sovrasov V. - All Rights Reserved #include "solver.hpp" #include -#include #include #include @@ -151,7 +150,7 @@ void NLPSolver::ClearDataStructures() Trial NLPSolver::Solve() { - bool needStop = false; + mNeedStop = false; InitDataStructures(); FirstIteration(); @@ -162,9 +161,9 @@ Trial NLPSolver::Solve() RefillQueue(); CalculateNextPoints(); MakeTrials(); - needStop = mMinDelta < mParameters.eps; + mNeedStop = mNeedStop || mMinDelta < mParameters.eps; mIterationsCounter++; - } while(mIterationsCounter < mParameters.itersLimit && !needStop); + } while(mIterationsCounter < mParameters.itersLimit && !mNeedStop); ClearDataStructures(); @@ -321,6 +320,9 @@ void NLPSolver::EstimateOptimum() { mOptimumEstimation = mNextPoints[i]; mNeedRefillQueue = true; + if (mOptimumEstimation.idx == mProblem->GetConstraintsNumber() && + mOptimumEstimation.g[mOptimumEstimation.idx] < mParameters.stopVal) + mNeedStop = true; } } } diff --git a/src/algs/ags/solver.hpp b/src/algs/ags/solver.hpp index c1af38b4..2c41d211 100644 --- a/src/algs/ags/solver.hpp +++ b/src/algs/ags/solver.hpp @@ -16,6 +16,7 @@ Copyright (C) 2018 Sovrasov V. - All Rights Reserved #include #include #include +#include namespace ags { @@ -23,6 +24,7 @@ namespace ags struct SolverParameters { double eps = 0.01; //method tolerance in Holder metric on 1d interval. Less value -- better search precision, less probability of early stop. + double stopVal = std::numeric_limits::lowest(); //method stops after objective becomes less than this value double r = 3; //reliability parameter. Higher value of r -- slower convergence, higher chance to cache the global minima. unsigned numPoints = 1; //number of new points per iteration. > 1 is useless in current implementation. unsigned itersLimit = 20000; // max number of iterations. @@ -62,6 +64,7 @@ class NLPSolver std::vector mCalculationsCounters; unsigned mIterationsCounter; bool mNeedRefillQueue; + bool mNeedStop; double mMinDelta; int mMaxIdx; From ef16e49b9d94ef0437a1f05521e61207a32a0434 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Mon, 16 Jul 2018 18:45:45 +0300 Subject: [PATCH 07/23] AGS: add stop by timer --- src/algs/ags/ags.cc | 2 +- src/algs/ags/solver.cc | 7 ++++++- src/algs/ags/solver.hpp | 1 + 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index a0a07757..92e07798 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -47,7 +47,7 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons ags::Trial optPoint; try { - optPoint = solver.Solve(); + optPoint = solver.Solve([stop](){ return (bool)nlopt_stop_evalstime(stop); }); } catch (const std::runtime_error& exp) { diff --git a/src/algs/ags/solver.cc b/src/algs/ags/solver.cc index 5e7935e0..1e03eb5b 100644 --- a/src/algs/ags/solver.cc +++ b/src/algs/ags/solver.cc @@ -149,6 +149,11 @@ void NLPSolver::ClearDataStructures() } Trial NLPSolver::Solve() +{ + return Solve([](){return false;}); +} + +Trial NLPSolver::Solve(std::function externalStopFunc) { mNeedStop = false; InitDataStructures(); @@ -161,7 +166,7 @@ Trial NLPSolver::Solve() RefillQueue(); CalculateNextPoints(); MakeTrials(); - mNeedStop = mNeedStop || mMinDelta < mParameters.eps; + mNeedStop = mNeedStop || mMinDelta < mParameters.eps || externalStopFunc(); mIterationsCounter++; } while(mIterationsCounter < mParameters.itersLimit && !mNeedStop); diff --git a/src/algs/ags/solver.hpp b/src/algs/ags/solver.hpp index 2c41d211..4c57df46 100644 --- a/src/algs/ags/solver.hpp +++ b/src/algs/ags/solver.hpp @@ -93,6 +93,7 @@ class NLPSolver void SetProblem(const std::vector& functions, const std::vector& leftBound, const std::vector& rightBound); + Trial Solve(std::function externalStopFunc); Trial Solve(); std::vector GetCalculationsStatistics() const; std::vector GetHolderConstantsEstimations() const; From 7f24c6f58fd0d0fce86b2b7106479b1ef1d965cf Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Mon, 16 Jul 2018 19:42:34 +0300 Subject: [PATCH 08/23] AGS: add correct return code for max_time stop --- src/algs/ags/ags.cc | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index 92e07798..fb5ae579 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -47,7 +47,14 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons ags::Trial optPoint; try { - optPoint = solver.Solve([stop](){ return (bool)nlopt_stop_evalstime(stop); }); + auto external_stop_func = [stop, &ret_code](){ + if (nlopt_stop_evalstime(stop)) { + ret_code = NLOPT_MAXTIME_REACHED; + return true; + } + else return false; + }; + optPoint = solver.Solve(external_stop_func); } catch (const std::runtime_error& exp) { From cf122f72243d69b582f6749d5af13a150fc24299 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Mon, 16 Jul 2018 20:00:19 +0300 Subject: [PATCH 09/23] AGS: stop instead of throwing an exception --- src/algs/ags/solver.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/algs/ags/solver.cc b/src/algs/ags/solver.cc index 1e03eb5b..5ad48cd6 100644 --- a/src/algs/ags/solver.cc +++ b/src/algs/ags/solver.cc @@ -298,7 +298,8 @@ void NLPSolver::CalculateNextPoints() mNextPoints[i].x = GetNextPointCoordinate(mNextIntervals[i]); if (mNextPoints[i].x >= mNextIntervals[i]->pr.x || mNextPoints[i].x <= mNextIntervals[i]->pl.x) - throw std::runtime_error("The next point is outside of the subdivided interval"); + mNeedStop = true; + //std::cout << "Warning: resolution of evolvent is not enough to continue the search"; mEvolvent.GetImage(mNextPoints[i].x, mNextPoints[i].y); } From 886d96323d1846dfb39aa0b0250e9ee4ee19d1a3 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Tue, 17 Jul 2018 16:21:06 +0300 Subject: [PATCH 10/23] Get rid of unused code --- src/algs/ags/ags.cc | 4 +- src/algs/ags/ags.h | 1 + src/algs/ags/evolvent.cc | 173 +------------------------------------- src/algs/ags/evolvent.hpp | 13 +-- 4 files changed, 7 insertions(+), 184 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index fb5ae579..62425ea1 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -9,6 +9,7 @@ double ags_r = 3; double eps_res = 0.001; unsigned evolvent_density = 12; +int ags_refine_loc = 0; int ags_verbose = 0; int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, @@ -39,6 +40,7 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons params.evolventDensity = evolvent_density; params.epsR = eps_res; params.stopVal = stop->minf_max; + params.refineSolution = (bool)ags_refine_loc; ags::NLPSolver solver; solver.SetParameters(params); @@ -48,7 +50,7 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons try { auto external_stop_func = [stop, &ret_code](){ - if (nlopt_stop_evalstime(stop)) { + if (nlopt_stop_time(stop)) { ret_code = NLOPT_MAXTIME_REACHED; return true; } diff --git a/src/algs/ags/ags.h b/src/algs/ags/ags.h index 72221903..6dabcd13 100644 --- a/src/algs/ags/ags.h +++ b/src/algs/ags/ags.h @@ -18,6 +18,7 @@ extern double eps_res; // parameter which prevents method from paying too much a // but global minima can be lost. extern unsigned evolvent_density; // density of evolvent. By default density is 2^-12 on hybercube [0,1]^N, // which means that maximum search accuracyis 2^-12. If search hypercube is large the density can be increased accordingly to achieve better accuracy. +extern int ags_refine_loc; //refine the final optimum using built-in local optimizer extern int ags_verbose; //print additional info #ifdef __cplusplus diff --git a/src/algs/ags/evolvent.cc b/src/algs/ags/evolvent.cc index 759995cc..b5d2c8b9 100644 --- a/src/algs/ags/evolvent.cc +++ b/src/algs/ags/evolvent.cc @@ -13,8 +13,6 @@ Copyright (C) 2018 Sovrasov V. - All Rights Reserved using namespace ags; void mapd(double x, int m, double* y, int n, int key = 1); /* map x to y */ -void invmad(int, double [], int, int *, double [], int, int); /* map y to x */ -void xyd(double *xx, int m, double y[], int n); /* get preimage */ Evolvent::Evolvent() { @@ -25,15 +23,13 @@ Evolvent::~Evolvent() { } -Evolvent::Evolvent(int dimension, int tightness, const double* lb, const double* ub, MapType type) +Evolvent::Evolvent(int dimension, int tightness, const double* lb, const double* ub) { assert(tightness > 2); mDimension = dimension; mTightness = tightness; - mMapType = type; - p2.resize(mDimension); mShiftScalars.resize(mDimension); mRho.resize(mDimension); for (int i = 0; i < mDimension; i++) @@ -42,19 +38,6 @@ Evolvent::Evolvent(int dimension, int tightness, const double* lb, const double* mShiftScalars[i] = 0.5*(lb[i] + ub[i]); } - switch (mMapType) - { - case Simple: - mMapKey = 1; - break; - case Linear: - mMapKey = 2; - break; - case Noninjective: - mMapKey = 3; - break; - } - mIsInitialized = true; } @@ -73,165 +56,13 @@ void Evolvent::TransformToSearchDomain(const double *y, double *z) void Evolvent::GetImage(double x, double y[]) { if(mDimension != 1) - mapd(x, mTightness, y, mDimension, mMapKey); + mapd(x, mTightness, y, mDimension); else y[0] = x - 0.5; TransformToSearchDomain(y, y); } -int Evolvent::GetAllPreimages(const double *p, double xp[]) -{ - int preimNumber = 1; - TransformToStandardCube(p, p2.data()); - if(mMapType == Noninjective) - invmad(mTightness, xp, noninjectiveMaxPreimages, &preimNumber, p2.data(), mDimension, 4); - else - xyd(xp, mTightness, p2.data(), mDimension); - - return preimNumber; -} - -void xyd(double *xx, int m, double y[], int n) -{ - /* calculate preimage x for nearest level m center to y */ - /* (x - left boundary point of level m interval) */ - int n1, nexp, l, iu[10], iv[10]; - - double x, r1; - double r; - int iw[11]; - int i, j, it, is; - void numbr(int *, const int, const int, int&, int*, int*); - - n1 = n - 1; - for (nexp = 1, i = 0; i0) ? floor(dr) : ceil(dr); - - del = 1. / (mne - dr); - d1 = del*(incr + 0.5); - for (kx = -1; kx= n) { - xyd(&x, m, y, n); - dr = x*mne; - dd = dr - fmod(dr, 1.0); - //dd = (dr>0) ? floor(dr) : ceil(dr); - dr = dd / nexp; - dd = dd - dr + fmod(dr, 1.0); - //dd = dd - ((dr>0) ? floor(dr) : ceil(dr)); - x = dd*del; - if (kx>kp) break; - k = kx++; /* label 9 */ - if (kx == 0) xp[0] = x; - else { - while (k >= 0) { - dr = fabs(x - xp[k]); /* label 11 */ - if (dr <= d1) { - for (kx--; k= 0) && (u[i] = (u[i] <= 0.0) ? 1 : -1)<0; i--); - if (i<0)break; - } - *kxx = ++kx; - -} - void mapd(double x, int m, double* y, int n, int key) { /* mapping y(x) : 1 - center, 2 - line, 3 - node */ diff --git a/src/algs/ags/evolvent.hpp b/src/algs/ags/evolvent.hpp index f85258f2..84b80c42 100644 --- a/src/algs/ags/evolvent.hpp +++ b/src/algs/ags/evolvent.hpp @@ -12,12 +12,6 @@ Copyright (C) 2018 Sovrasov V. - All Rights Reserved namespace ags { -enum MapType { - Simple = 1, Linear = 2, Noninjective = 3 -}; - -constexpr int noninjectiveMaxPreimages = 32; - class Evolvent { protected: @@ -26,23 +20,18 @@ class Evolvent std::vector mRho; std::vector mShiftScalars; - std::vector p2; void TransformToStandardCube(const double *y, double *z); void TransformToSearchDomain(const double *y, double *z); bool mIsInitialized; -private: - MapType mMapType; - int mMapKey; public: Evolvent(); - Evolvent(int dimension, int tightness, const double* lb, const double* ub, MapType type = Simple); + Evolvent(int dimension, int tightness, const double* lb, const double* ub); ~Evolvent(); virtual void GetImage(double x, double y[]); - virtual int GetAllPreimages(const double* p, double xp[]); }; } From bd981b9ce2afd899f16159e3289b9711194761aa Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Tue, 17 Jul 2018 16:30:39 +0300 Subject: [PATCH 11/23] AGS: use NLOPT_CXX11 macro --- CMakeLists.txt | 2 ++ src/api/general.c | 2 +- src/api/optimize.c | 2 +- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d3af37f..14b8954d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -204,6 +204,8 @@ set (NLOPT_SOURCES if (NLOPT_CXX) list (APPEND NLOPT_SOURCES src/algs/stogo/global.cc src/algs/stogo/linalg.cc src/algs/stogo/local.cc src/algs/stogo/stogo.cc src/algs/stogo/tools.cc src/algs/stogo/global.h src/algs/stogo/linalg.h src/algs/stogo/local.h src/algs/stogo/stogo_config.h src/algs/stogo/stogo.h src/algs/stogo/tools.h) +endif () +if (NLOPT_CXX11) list (APPEND NLOPT_SOURCES src/algs/ags/data_types.hpp src/algs/ags/evolvent.hpp src/algs/ags/evolvent.cc src/algs/ags/solver.hpp src/algs/ags/solver.cc src/algs/ags/local_optimizer.hpp src/algs/ags/local_optimizer.cc src/algs/ags/ags.h src/algs/ags/ags.cc) endif () diff --git a/src/api/general.c b/src/api/general.c index 5b536ee4..6e75074a 100644 --- a/src/api/general.c +++ b/src/api/general.c @@ -82,7 +82,7 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { "Sequential Quadratic Programming (SQP) (local, derivative)", "CCSA (Conservative Convex Separable Approximations) with simple quadratic approximations (local, derivative)", "ESCH evolutionary strategy", -#ifdef NLOPT_CXX +#ifdef NLOPT_CXX11 "AGS (global, no-derivative)" #else "AGS (NOT COMPILED)" diff --git a/src/api/optimize.c b/src/api/optimize.c index 03f45952..9a45d9d3 100644 --- a/src/api/optimize.c +++ b/src/api/optimize.c @@ -507,7 +507,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) } case NLOPT_AGS: -#ifdef NLOPT_CXX +#ifdef NLOPT_CXX11 if (!finite_domain(n, lb, ub)) RETURN_ERR(NLOPT_INVALID_ARGS, opt, "finite domain required for global algorithm"); From 1a3dc058e9a7414bcf42e052677fac88ac873faa Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Tue, 17 Jul 2018 18:52:55 +0300 Subject: [PATCH 12/23] AGS: updated documentation --- doc/docs/NLopt_Algorithms.md | 18 ++++++++++++++++++ src/algs/ags/ags.h | 2 ++ 2 files changed, 20 insertions(+) diff --git a/doc/docs/NLopt_Algorithms.md b/doc/docs/NLopt_Algorithms.md index ba5ff25d..b482f10b 100644 --- a/doc/docs/NLopt_Algorithms.md +++ b/doc/docs/NLopt_Algorithms.md @@ -123,6 +123,24 @@ Some references on StoGO are: Only bound-constrained problems are supported by this algorithm. +### AGS + +This algorithm adapted from [this repo](https://github.com/sovrasov/glob_search_nlp_solver). +AGS can handle arbitrary objectives and nonlinear inequality constraints. Also bound constraints are required for this method. To guarantee convergence, objectives and constraints should satisfy the Lipschitz condition on the specified hyperrectangle. +AGS employs the Hilbert curve to reduce the source problem to the univariate one. The algorithm divides the univariate space into intervals, generating new points by using posterior probabilities. On each trial AGS tries to evaluate the constraints consequently one by one. If some constraint is violated at this point, the next ones won't be evaluated. If all constraints are preserved, i.e. the trial point is feasible, AGS will evaluate the objective. Thus, some of constraints (except the first one) and objective can be partially undefined inside the search hyperrectangle. + +Limitation of machine arithmetic don't allow to build a tight approximation for Hilbert when the space dimension is greater than 5, so this implementation of AGS is restricted in that sense. + +AGS, like StoGO, is written in C++, but it requires C++11. If the library is built with [C++](NLopt_Installation.md) and compiler supports C++11, AGS will be built too. + +AGS is specified within NLopt by `NLOPT_AGS`. Additional parameters of AGS which are not adjustable from the common NLOpt interface are declared and described in `ags.h`. +References: +- Yaroslav D. Sergeyev, Dmitri L. Markin: An algorithm for solving global optimization problems with nonlinear constraints, Journal of Global Optimization, 7(4), pp 407–419, 1995 +- Strongin R.G., Sergeyev Ya.D., 2000. Global optimization with non-convex constraints. Sequential and parallel algorithms. Kluwer Academic +Publishers, Dordrecht. +- Gergel V. and Lebedev I.: Heterogeneous Parallel Computations for Solving Global Optimization Problems. Proc. Comput. Science 66, pp. 53–62 (2015) +- [Implementation] (https://github.com/sovrasov/multicriterial-go) of AGS for constrained multi-objective problems. + ### ISRES (Improved Stochastic Ranking Evolution Strategy) This is my implementation of the "Improved Stochastic Ranking Evolution Strategy" (ISRES) algorithm for nonlinearly-constrained global optimization (or at least semi-global; although it has heuristics to escape local optima, I'm not aware of a convergence proof), based on the method described in: diff --git a/src/algs/ags/ags.h b/src/algs/ags/ags.h index 6dabcd13..bcb4d327 100644 --- a/src/algs/ags/ags.h +++ b/src/algs/ags/ags.h @@ -10,6 +10,8 @@ extern "C" { #endif +//The algorithm supports 3 types of stop criterions: stop by execution time, stop by value and stop by exceeding limit of iterations. + int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, double *x, double *minf, const double *l, const double *u, nlopt_stopping *stop); From d7f3030aba65f57bd69e606741c26cfa1ca3e713 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Tue, 17 Jul 2018 18:57:08 +0300 Subject: [PATCH 13/23] AGS: fix wrong ifdef --- src/api/optimize.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/api/optimize.c b/src/api/optimize.c index 9a45d9d3..4dd46d08 100644 --- a/src/api/optimize.c +++ b/src/api/optimize.c @@ -33,6 +33,8 @@ #ifdef NLOPT_CXX # include "stogo.h" +#endif +#ifdef NLOPT_CXX11 # include "ags.h" #endif From fb9a8b7673cc5c899bb0aba4db458e1fe3ceca21 Mon Sep 17 00:00:00 2001 From: sovrasov Date: Tue, 17 Jul 2018 19:41:00 +0300 Subject: [PATCH 14/23] AGS: use spaces rather than tabs --- src/algs/ags/ags.cc | 94 ++++++++++++++++++++++----------------------- 1 file changed, 47 insertions(+), 47 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index 62425ea1..74065b8c 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -15,40 +15,40 @@ int ags_verbose = 0; int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, double *x, double *minf, const double *l, const double *u, nlopt_stopping *stop) { - int ret_code = NLOPT_SUCCESS; + int ret_code = NLOPT_SUCCESS; - if (n > ags::solverMaxDim) - return NLOPT_INVALID_ARGS; - if(m != nlopt_count_constraints(m, fc) || m > ags::solverMaxConstraints) - return NLOPT_INVALID_ARGS; + if (n > ags::solverMaxDim) + return NLOPT_INVALID_ARGS; + if(m != nlopt_count_constraints(m, fc) || m > ags::solverMaxConstraints) + return NLOPT_INVALID_ARGS; - std::vector lb(l, l + n); - std::vector ub(u, u + n); - std::vector functions; - for (unsigned i = 0; i < m; i++) - { - if (fc[i].m != 1) - return NLOPT_INVALID_ARGS; - functions.push_back([fc, data, n, i](const double* x) {return fc[i].f(n, x, NULL, data);}); - } - functions.push_back([func, data, n](const double* x) {return func(n, x, NULL, data);}); + std::vector lb(l, l + n); + std::vector ub(u, u + n); + std::vector functions; + for (unsigned i = 0; i < m; i++) + { + if (fc[i].m != 1) + return NLOPT_INVALID_ARGS; + functions.push_back([fc, data, n, i](const double* x) {return fc[i].f(n, x, NULL, data);}); + } + functions.push_back([func, data, n](const double* x) {return func(n, x, NULL, data);}); - ags::SolverParameters params; - params.r = ags_r; - params.itersLimit = stop->maxeval; - params.eps = 1e-64; - params.evolventDensity = evolvent_density; - params.epsR = eps_res; - params.stopVal = stop->minf_max; + ags::SolverParameters params; + params.r = ags_r; + params.itersLimit = stop->maxeval; + params.eps = 1e-64; + params.evolventDensity = evolvent_density; + params.epsR = eps_res; + params.stopVal = stop->minf_max; params.refineSolution = (bool)ags_refine_loc; - ags::NLPSolver solver; - solver.SetParameters(params); - solver.SetProblem(functions, lb, ub); + ags::NLPSolver solver; + solver.SetParameters(params); + solver.SetProblem(functions, lb, ub); - ags::Trial optPoint; - try - { + ags::Trial optPoint; + try + { auto external_stop_func = [stop, &ret_code](){ if (nlopt_stop_time(stop)) { ret_code = NLOPT_MAXTIME_REACHED; @@ -56,16 +56,16 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons } else return false; }; - optPoint = solver.Solve(external_stop_func); - } - catch (const std::runtime_error& exp) - { - std::cerr << "AGS internal error: " << std::string(exp.what()) << std::endl; - return NLOPT_FAILURE; - } + optPoint = solver.Solve(external_stop_func); + } + catch (const std::runtime_error& exp) + { + std::cerr << "AGS internal error: " << std::string(exp.what()) << std::endl; + return NLOPT_FAILURE; + } - if (ags_verbose) - { + if (ags_verbose) + { auto calcCounters = solver.GetCalculationsStatistics(); auto holderConstEstimations = solver.GetHolderConstantsEstimations(); @@ -80,18 +80,18 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons if (optPoint.idx != m) std::cout << "Feasible point not found" << "\n"; std::cout << std::string(40, '-') << std::endl; - } + } - if (m == optPoint.idx) - { - memcpy(x, optPoint.y, n*sizeof(x[0])); - *minf = optPoint.g[optPoint.idx]; - } - else //feasible point not found. - return NLOPT_FAILURE; + if (m == optPoint.idx) + { + memcpy(x, optPoint.y, n*sizeof(x[0])); + *minf = optPoint.g[optPoint.idx]; + } + else //feasible point not found. + return NLOPT_FAILURE; - if (solver.GetCalculationsStatistics()[0] >= params.itersLimit) - return NLOPT_MAXEVAL_REACHED; + if (solver.GetCalculationsStatistics()[0] >= params.itersLimit) + return NLOPT_MAXEVAL_REACHED; return ret_code; } From bdf87001eb4330378c46ddcdf1f1e246cf7fe69b Mon Sep 17 00:00:00 2001 From: sovrasov Date: Tue, 17 Jul 2018 21:11:22 +0300 Subject: [PATCH 15/23] AGS: fix enum name --- doc/docs/NLopt_Algorithms.md | 2 +- src/api/nlopt.h | 2 +- src/api/optimize.c | 4 ++-- src/api/options.c | 3 ++- 4 files changed, 6 insertions(+), 5 deletions(-) diff --git a/doc/docs/NLopt_Algorithms.md b/doc/docs/NLopt_Algorithms.md index b482f10b..8fb4da3c 100644 --- a/doc/docs/NLopt_Algorithms.md +++ b/doc/docs/NLopt_Algorithms.md @@ -133,7 +133,7 @@ Limitation of machine arithmetic don't allow to build a tight approximation for AGS, like StoGO, is written in C++, but it requires C++11. If the library is built with [C++](NLopt_Installation.md) and compiler supports C++11, AGS will be built too. -AGS is specified within NLopt by `NLOPT_AGS`. Additional parameters of AGS which are not adjustable from the common NLOpt interface are declared and described in `ags.h`. +AGS is specified within NLopt by `NLOPT_GN_AGS`. Additional parameters of AGS which are not adjustable from the common NLOpt interface are declared and described in `ags.h`. References: - Yaroslav D. Sergeyev, Dmitri L. Markin: An algorithm for solving global optimization problems with nonlinear constraints, Journal of Global Optimization, 7(4), pp 407–419, 1995 - Strongin R.G., Sergeyev Ya.D., 2000. Global optimization with non-convex constraints. Sequential and parallel algorithms. Kluwer Academic diff --git a/src/api/nlopt.h b/src/api/nlopt.h index 5ced305f..998db3cd 100644 --- a/src/api/nlopt.h +++ b/src/api/nlopt.h @@ -151,7 +151,7 @@ typedef enum { NLOPT_GN_ESCH, - NLOPT_AGS, + NLOPT_GN_AGS, NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */ } nlopt_algorithm; diff --git a/src/api/optimize.c b/src/api/optimize.c index 4dd46d08..71614d21 100644 --- a/src/api/optimize.c +++ b/src/api/optimize.c @@ -350,7 +350,7 @@ static int elimdim_wrapcheck(nlopt_opt opt) case NLOPT_LN_SBPLX: case NLOPT_GN_ISRES: case NLOPT_GN_ESCH: - case NLOPT_AGS: + case NLOPT_GN_AGS: case NLOPT_GD_STOGO: case NLOPT_GD_STOGO_RAND: return 1; @@ -508,7 +508,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) break; } - case NLOPT_AGS: + case NLOPT_GN_AGS: #ifdef NLOPT_CXX11 if (!finite_domain(n, lb, ub)) RETURN_ERR(NLOPT_INVALID_ARGS, opt, diff --git a/src/api/options.c b/src/api/options.c index 04e0bfe7..fab409bf 100644 --- a/src/api/options.c +++ b/src/api/options.c @@ -437,7 +437,8 @@ static int inequality_ok(nlopt_algorithm algorithm) { || AUGLAG_ALG(algorithm) || algorithm == NLOPT_GN_ISRES || algorithm == NLOPT_GN_ORIG_DIRECT - || algorithm == NLOPT_GN_ORIG_DIRECT_L); + || algorithm == NLOPT_GN_ORIG_DIRECT_L + || algorithm == NLOPT_GN_AGS); } nlopt_result From b0b183449e35db0f5465c99c8114af9adcf10f09 Mon Sep 17 00:00:00 2001 From: sovrasov Date: Tue, 17 Jul 2018 21:32:41 +0300 Subject: [PATCH 16/23] AGS: fix wrong calculation of constraints --- src/algs/ags/ags.cc | 9 +++++++-- src/algs/ags/ags.h | 1 + 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index 74065b8c..c68ca974 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -6,6 +6,7 @@ #include #include +double ags_eps = 0; double ags_r = 3; double eps_res = 0.001; unsigned evolvent_density = 12; @@ -29,14 +30,18 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons { if (fc[i].m != 1) return NLOPT_INVALID_ARGS; - functions.push_back([fc, data, n, i](const double* x) {return fc[i].f(n, x, NULL, data);}); + functions.push_back([fc, data, n, i](const double* x) { + double val = 0; + nlopt_eval_constraint(&val, NULL, &fc[i], n, x); + return val; + }); } functions.push_back([func, data, n](const double* x) {return func(n, x, NULL, data);}); ags::SolverParameters params; params.r = ags_r; params.itersLimit = stop->maxeval; - params.eps = 1e-64; + params.eps = ags_eps; params.evolventDensity = evolvent_density; params.epsR = eps_res; params.stopVal = stop->minf_max; diff --git a/src/algs/ags/ags.h b/src/algs/ags/ags.h index bcb4d327..809340a6 100644 --- a/src/algs/ags/ags.h +++ b/src/algs/ags/ags.h @@ -15,6 +15,7 @@ extern "C" { int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_constraint *fc, double *x, double *minf, const double *l, const double *u, nlopt_stopping *stop); +extern double ags_eps; //method tolerance in Holder metric on 1d interval. Less value -- better search precision, less probability of early stop. extern double ags_r; //reliability parameter. Higher value of r -- slower convergence, higher chance to cache the global minima. extern double eps_res; // parameter which prevents method from paying too much attention to constraints. Greater values of this parameter speed up convergence, // but global minima can be lost. From c7b67be658a1612a60117412c0abc3e189cefea6 Mon Sep 17 00:00:00 2001 From: sovrasov Date: Tue, 17 Jul 2018 21:34:18 +0300 Subject: [PATCH 17/23] AGS: add an example of problem with nonlinear constraints --- src/algs/ags/tst.cc | 58 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 src/algs/ags/tst.cc diff --git a/src/algs/ags/tst.cc b/src/algs/ags/tst.cc new file mode 100644 index 00000000..cbf6e995 --- /dev/null +++ b/src/algs/ags/tst.cc @@ -0,0 +1,58 @@ +#define _USE_MATH_DEFINES +#include +#include +#include + +#include "nlopt.hpp" + +double f_obj(const std::vector& x, std::vector&, void*) +{ + return -1.5*pow(x[0], 2) * exp(1 - pow(x[0], 2) + - 20.25*pow(x[0] - x[1], 2)) - pow(0.5 * (x[1] - 1)*(x[0]- 1), 4) + * exp(2 - pow(0.5 * (x[0] - 1), 4) - pow(x[1] - 1, 4)); +} + +double f_c0(const std::vector& x, std::vector&, void*) +{ + return 0.01*(pow(x[0] - 2.2, 2) + pow(x[1] - 1.2, 2) - 2.25); +} +double f_c1(const std::vector& x, std::vector&, void*) +{ + return 100 * (1 - pow(x[0] - 2, 2) / 1.44 - pow(0.5*x[1], 2)); +} +double f_c2(const std::vector& x, std::vector&, void*) +{ + return 10 * (x[1] - 1.5 - 1.5*sin(2*M_PI*(x[0] - 1.75))); +} + +int ags_verbose = 1; +double ags_eps = 0.001; +double eps_res = 0.1; + +int main(int argc, char **argv) +{ + nlopt::opt opt(nlopt::GN_AGS, 2); + + opt.set_lower_bounds({0, -1}); + opt.set_upper_bounds({4, 3}); + + opt.add_inequality_constraint(f_c0, NULL, 0); + opt.add_inequality_constraint(f_c1, NULL, 0); + opt.add_inequality_constraint(f_c2, NULL, 0); + opt.set_min_objective(f_obj, NULL); + opt.set_maxeval(2000); + + double minf; + std::vector x(2); + + try { + nlopt::result result = opt.optimize(x, minf); + std::cout << "found minimum at f(" << x[0] << "," << x[1] << ") = " + << std::setprecision(10) << minf << std::endl; + } + catch(std::exception &e) { + std::cout << "nlopt failed: " << e.what() << std::endl; + } + + return EXIT_SUCCESS; +} From 7743d22530c39c09d52048890206f95dd91174a7 Mon Sep 17 00:00:00 2001 From: sovrasov Date: Tue, 17 Jul 2018 23:17:53 +0300 Subject: [PATCH 18/23] AGS: update docs --- doc/docs/NLopt_Algorithms.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/docs/NLopt_Algorithms.md b/doc/docs/NLopt_Algorithms.md index 8fb4da3c..4f0a769a 100644 --- a/doc/docs/NLopt_Algorithms.md +++ b/doc/docs/NLopt_Algorithms.md @@ -31,7 +31,7 @@ Better yet, run some algorithm for a really long time until the minimum *f* Global optimization ------------------- -All of the global-optimization algorithms currently require you to specify bound constraints on all the optimization parameters. Of these algorithms, only ISRES and ORIG_DIRECT support nonlinear inequality constraints, and only ISRES supports nonlinear equality constraints. (However, any of them can be applied to nonlinearly constrained problems by combining them with the [augmented Lagrangian method](#Augmented_Lagrangian_algorithm.md) below.) +All of the global-optimization algorithms currently require you to specify bound constraints on all the optimization parameters. Of these algorithms, only ISRES, AGS and ORIG_DIRECT support nonlinear inequality constraints, and only ISRES supports nonlinear equality constraints. (However, any of them can be applied to nonlinearly constrained problems by combining them with the [augmented Lagrangian method](#Augmented_Lagrangian_algorithm.md) below.) **Something you should consider** is that, after running the global optimization, it is often worthwhile to then use the global optimum as a starting point for a local optimization to "polish" the optimum to a greater accuracy. (Many of the global optimization algorithms devote more effort to searching the global parameter space than in finding the precise position of the local optimum accurately.) @@ -127,19 +127,19 @@ Only bound-constrained problems are supported by this algorithm. This algorithm adapted from [this repo](https://github.com/sovrasov/glob_search_nlp_solver). AGS can handle arbitrary objectives and nonlinear inequality constraints. Also bound constraints are required for this method. To guarantee convergence, objectives and constraints should satisfy the Lipschitz condition on the specified hyperrectangle. -AGS employs the Hilbert curve to reduce the source problem to the univariate one. The algorithm divides the univariate space into intervals, generating new points by using posterior probabilities. On each trial AGS tries to evaluate the constraints consequently one by one. If some constraint is violated at this point, the next ones won't be evaluated. If all constraints are preserved, i.e. the trial point is feasible, AGS will evaluate the objective. Thus, some of constraints (except the first one) and objective can be partially undefined inside the search hyperrectangle. +AGS is derivative-free and employs the Hilbert curve to reduce the source problem to the univariate one. The algorithm divides the univariate space into intervals, generating new points by using posterior probabilities. On each trial AGS tries to evaluate the constraints consequently one by one. If some constraint is violated at this point, the next ones won't be evaluated. If all constraints are preserved, i.e. the trial point is feasible, AGS will evaluate the objective. Thus, some of constraints (except the first one) and objective can be partially undefined inside the search hyperrectangle. Current implementation of AGS doesn't support vector constraints. -Limitation of machine arithmetic don't allow to build a tight approximation for Hilbert when the space dimension is greater than 5, so this implementation of AGS is restricted in that sense. +Limitations of the machine arithmetic don't allow to build a tight approximation for Hilbert when the space dimension is greater than 5, so this implementation of AGS is restricted in that sense. AGS, like StoGO, is written in C++, but it requires C++11. If the library is built with [C++](NLopt_Installation.md) and compiler supports C++11, AGS will be built too. -AGS is specified within NLopt by `NLOPT_GN_AGS`. Additional parameters of AGS which are not adjustable from the common NLOpt interface are declared and described in `ags.h`. +AGS is specified within NLopt by `NLOPT_GN_AGS`. Additional parameters of AGS which are not adjustable from the common NLOpt interface are declared and described in `ags.h`. Also an example of solving a constrained problem is given in the AGS source folder. References: - Yaroslav D. Sergeyev, Dmitri L. Markin: An algorithm for solving global optimization problems with nonlinear constraints, Journal of Global Optimization, 7(4), pp 407–419, 1995 - Strongin R.G., Sergeyev Ya.D., 2000. Global optimization with non-convex constraints. Sequential and parallel algorithms. Kluwer Academic Publishers, Dordrecht. - Gergel V. and Lebedev I.: Heterogeneous Parallel Computations for Solving Global Optimization Problems. Proc. Comput. Science 66, pp. 53–62 (2015) -- [Implementation] (https://github.com/sovrasov/multicriterial-go) of AGS for constrained multi-objective problems. +- [Implementation](https://github.com/sovrasov/multicriterial-go) of AGS for constrained multi-objective problems. ### ISRES (Improved Stochastic Ranking Evolution Strategy) From 7b6aa88d4e4c4fbecabbb853c3b0d9fc4a289659 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Wed, 18 Jul 2018 10:24:14 +0300 Subject: [PATCH 19/23] Fix minor issues --- doc/docs/NLopt_Algorithms.md | 2 +- src/api/options.c | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/docs/NLopt_Algorithms.md b/doc/docs/NLopt_Algorithms.md index 4f0a769a..35597632 100644 --- a/doc/docs/NLopt_Algorithms.md +++ b/doc/docs/NLopt_Algorithms.md @@ -31,7 +31,7 @@ Better yet, run some algorithm for a really long time until the minimum *f* Global optimization ------------------- -All of the global-optimization algorithms currently require you to specify bound constraints on all the optimization parameters. Of these algorithms, only ISRES, AGS and ORIG_DIRECT support nonlinear inequality constraints, and only ISRES supports nonlinear equality constraints. (However, any of them can be applied to nonlinearly constrained problems by combining them with the [augmented Lagrangian method](#Augmented_Lagrangian_algorithm.md) below.) +All of the global-optimization algorithms currently require you to specify bound constraints on all the optimization parameters. Of these algorithms, only ISRES, AGS, and ORIG_DIRECT support nonlinear inequality constraints, and only ISRES supports nonlinear equality constraints. (However, any of them can be applied to nonlinearly constrained problems by combining them with the [augmented Lagrangian method](#Augmented_Lagrangian_algorithm.md) below.) **Something you should consider** is that, after running the global optimization, it is often worthwhile to then use the global optimum as a starting point for a local optimization to "polish" the optimum to a greater accuracy. (Many of the global optimization algorithms devote more effort to searching the global parameter space than in finding the precise position of the local optimum accurately.) diff --git a/src/api/options.c b/src/api/options.c index fab409bf..c60729b6 100644 --- a/src/api/options.c +++ b/src/api/options.c @@ -438,7 +438,7 @@ static int inequality_ok(nlopt_algorithm algorithm) { || algorithm == NLOPT_GN_ISRES || algorithm == NLOPT_GN_ORIG_DIRECT || algorithm == NLOPT_GN_ORIG_DIRECT_L - || algorithm == NLOPT_GN_AGS); + || algorithm == NLOPT_GN_AGS); } nlopt_result From f66ac2d86ba5618a6181a13272ae93d175616bb5 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Wed, 18 Jul 2018 11:13:13 +0300 Subject: [PATCH 20/23] AGS: allow up to 10 dimenstions instead of 5 --- doc/docs/NLopt_Algorithms.md | 2 +- src/algs/ags/ags.cc | 6 +++++- src/algs/ags/data_types.hpp | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/doc/docs/NLopt_Algorithms.md b/doc/docs/NLopt_Algorithms.md index 35597632..59fddffe 100644 --- a/doc/docs/NLopt_Algorithms.md +++ b/doc/docs/NLopt_Algorithms.md @@ -129,7 +129,7 @@ This algorithm adapted from [this repo](https://github.com/sovrasov/glob_search_ AGS can handle arbitrary objectives and nonlinear inequality constraints. Also bound constraints are required for this method. To guarantee convergence, objectives and constraints should satisfy the Lipschitz condition on the specified hyperrectangle. AGS is derivative-free and employs the Hilbert curve to reduce the source problem to the univariate one. The algorithm divides the univariate space into intervals, generating new points by using posterior probabilities. On each trial AGS tries to evaluate the constraints consequently one by one. If some constraint is violated at this point, the next ones won't be evaluated. If all constraints are preserved, i.e. the trial point is feasible, AGS will evaluate the objective. Thus, some of constraints (except the first one) and objective can be partially undefined inside the search hyperrectangle. Current implementation of AGS doesn't support vector constraints. -Limitations of the machine arithmetic don't allow to build a tight approximation for Hilbert when the space dimension is greater than 5, so this implementation of AGS is restricted in that sense. +Limitations of the machine arithmetic don't allow to build a tight approximation for Hilbert when the space dimension is greater than 5, so this implementation of AGS is restricted in that sense. It supports up to 10 dimensions, but the method can stop early in case of 6 and more ones. AGS, like StoGO, is written in C++, but it requires C++11. If the library is built with [C++](NLopt_Installation.md) and compiler supports C++11, AGS will be built too. diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index c68ca974..3310635d 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -5,6 +5,7 @@ #include "solver.hpp" #include #include +#include double ags_eps = 0; double ags_r = 3; @@ -23,6 +24,9 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons if(m != nlopt_count_constraints(m, fc) || m > ags::solverMaxConstraints) return NLOPT_INVALID_ARGS; + if (ags_verbose && n > 5) + std::cout << "Warning: AGS is unstable when dimension > 5" << std::endl; + std::vector lb(l, l + n); std::vector ub(u, u + n); std::vector functions; @@ -63,7 +67,7 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons }; optPoint = solver.Solve(external_stop_func); } - catch (const std::runtime_error& exp) + catch (const std::exception& exp) { std::cerr << "AGS internal error: " << std::string(exp.what()) << std::endl; return NLOPT_FAILURE; diff --git a/src/algs/ags/data_types.hpp b/src/algs/ags/data_types.hpp index 94524284..9b29d8e2 100644 --- a/src/algs/ags/data_types.hpp +++ b/src/algs/ags/data_types.hpp @@ -16,7 +16,7 @@ Copyright (C) 2018 Sovrasov V. - All Rights Reserved namespace ags { -const unsigned solverMaxDim = 5; +const unsigned solverMaxDim = 10; const unsigned solverMaxConstraints = 10; template From 0b603636d85fb4e87977872ab888c5e9cc1bd2d2 Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Wed, 18 Jul 2018 12:31:45 +0300 Subject: [PATCH 21/23] AGS: fix warnings --- src/algs/ags/ags.cc | 4 ++-- src/algs/ags/evolvent.cc | 2 +- src/algs/ags/solver.cc | 10 +++++----- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index 3310635d..644040d4 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -86,12 +86,12 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons for (size_t i = 0; i < holderConstEstimations.size() - 1; i++) std::cout << "Estimation of Holder constant of function # " << i << ": " << holderConstEstimations[i] << "\n"; std::cout << "Estimation of Holder constant of objective: " << holderConstEstimations.back() << "\n"; - if (optPoint.idx != m) + if (optPoint.idx != (int)m) std::cout << "Feasible point not found" << "\n"; std::cout << std::string(40, '-') << std::endl; } - if (m == optPoint.idx) + if ((int)m == optPoint.idx) { memcpy(x, optPoint.y, n*sizeof(x[0])); *minf = optPoint.g[optPoint.idx]; diff --git a/src/algs/ags/evolvent.cc b/src/algs/ags/evolvent.cc index b5d2c8b9..b72cf86d 100644 --- a/src/algs/ags/evolvent.cc +++ b/src/algs/ags/evolvent.cc @@ -72,7 +72,7 @@ void mapd(double x, int m, double* y, int n, int key) double d, mne, dd, dr;//,tmp; double p, r; int iw[11]; - int it, is, i, j, k; + int it, is = 0, i, j, k; void node(int is, int n1, int nexp, int& l, int& iq, int iu[], int iv[]); p = 0.0; diff --git a/src/algs/ags/solver.cc b/src/algs/ags/solver.cc index 5ad48cd6..289592de 100644 --- a/src/algs/ags/solver.cc +++ b/src/algs/ags/solver.cc @@ -73,7 +73,7 @@ void NLPSolver::SetParameters(const SolverParameters& params) void NLPSolver::SetProblem(std::shared_ptr> problem) { mProblem = problem; - NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= solverMaxConstraints, + NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= (int)solverMaxConstraints, "Current implementation supports up to " + std::to_string(solverMaxConstraints) + " nonlinear inequality constraints"); InitLocalOptimizer(); @@ -85,7 +85,7 @@ void NLPSolver::SetProblem(const std::vector& functions, NLP_SOLVER_ASSERT(leftBound.size() == rightBound.size(), "Inconsistent dimensions of bounds"); NLP_SOLVER_ASSERT(leftBound.size() > 0, "Zero problem dimension"); mProblem = std::make_shared(functions, leftBound, rightBound); - NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= solverMaxConstraints, + NLP_SOLVER_ASSERT(mProblem->GetConstraintsNumber() <= (int)solverMaxConstraints, "Current implementation supports up to " + std::to_string(solverMaxConstraints) + " nonlinear inequality constraints"); InitLocalOptimizer(); @@ -321,8 +321,8 @@ void NLPSolver::EstimateOptimum() for (size_t i = 0; i < mNextPoints.size(); i++) { if (mOptimumEstimation.idx < mNextPoints[i].idx || - mOptimumEstimation.idx == mNextPoints[i].idx && - mOptimumEstimation.g[mOptimumEstimation.idx] > mNextPoints[i].g[mNextPoints[i].idx]) + (mOptimumEstimation.idx == mNextPoints[i].idx && + mOptimumEstimation.g[mOptimumEstimation.idx] > mNextPoints[i].g[mNextPoints[i].idx])) { mOptimumEstimation = mNextPoints[i]; mNeedRefillQueue = true; @@ -335,7 +335,7 @@ void NLPSolver::EstimateOptimum() void NLPSolver::UpdateH(double newValue, int index) { - if (newValue > mHEstimations[index] || mHEstimations[index] == 1.0 && newValue > zeroHLevel) + if (newValue > mHEstimations[index] || (mHEstimations[index] == 1.0 && newValue > zeroHLevel)) { mHEstimations[index] = newValue; mNeedRefillQueue = true; From 3cceeed4c8ad5cbf9fab1f33dac6b256fc3b8aeb Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Wed, 18 Jul 2018 13:49:12 +0300 Subject: [PATCH 22/23] AGS: fix zero evaluations counter, set default maxeval --- src/algs/ags/ags.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/algs/ags/ags.cc b/src/algs/ags/ags.cc index 644040d4..63f3a3ee 100644 --- a/src/algs/ags/ags.cc +++ b/src/algs/ags/ags.cc @@ -40,11 +40,13 @@ int ags_minimize(unsigned n, nlopt_func func, void *data, unsigned m, nlopt_cons return val; }); } - functions.push_back([func, data, n](const double* x) {return func(n, x, NULL, data);}); + functions.push_back([func, data, n, stop](const double* x) { + ++ *(stop->nevals_p); + return func(n, x, NULL, data);}); ags::SolverParameters params; params.r = ags_r; - params.itersLimit = stop->maxeval; + params.itersLimit = stop->maxeval != 0 ? stop->maxeval : 5000; params.eps = ags_eps; params.evolventDensity = evolvent_density; params.epsR = eps_res; From 5e1e09f95e5483de772046992d0635264829bbde Mon Sep 17 00:00:00 2001 From: Vladislav Sovrasov Date: Thu, 19 Jul 2018 14:22:29 +0300 Subject: [PATCH 23/23] AGS: fix generation of test suite --- CMakeLists.txt | 4 +++- test/CMakeLists.txt | 6 +++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 14b8954d..6f91eb00 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -143,7 +143,9 @@ if (NLOPT_CXX OR NLOPT_PYTHON OR NLOPT_GUILE OR NLOPT_OCTAVE) check_cxx_compiler_flag ("-std=c++11" SUPPORTS_STDCXX11) if (SUPPORTS_STDCXX11) set (CMAKE_CXX_FLAGS "-std=c++11 ${CMAKE_CXX_FLAGS}") - set (NLOPT_CXX11 ON) + if (NLOPT_CXX) + set (NLOPT_CXX11 ON) + endif () endif () else() message (FATAL_ERROR "The compiler doesn't support CXX.") diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3d4f0d09..8871312d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -12,10 +12,14 @@ foreach (algo_index RANGE 29)# 43 set (enable_ TRUE) # cxx stogo if (NOT NLOPT_CXX) - if (algo_index STREQUAL 8 OR algo_index STREQUAL 9 OR algo_index STREQUAL 43) + if (algo_index STREQUAL 8 OR algo_index STREQUAL 9) set (enable_ FALSE) endif () endif () + # cxx11 ags + if (NOT NLOPT_CXX11 AND algo_index STREQUAL 43) + set (enable_ FALSE) + endif () # L-BFGS if (algo_index STREQUAL 10) set (enable_ FALSE)