From 392ba7d9aa956d47c96372972d8fbc6efa8c9f89 Mon Sep 17 00:00:00 2001 From: Julien Schueller Date: Thu, 4 Jun 2026 16:35:58 +0200 Subject: [PATCH] Add C++ implementation of COBYLA see #136 --- .github/workflows/test_cpp.yml | 43 ++ CMakeLists.txt | 6 + cpp/CMakeLists.txt | 15 + cpp/README.md | 46 ++ cpp/examples/CMakeLists.txt | 5 + cpp/examples/cobyla_example.cpp | 64 +++ cpp/examples/rosenbrock.cpp | 148 +++++++ cpp/src/CMakeLists.txt | 48 ++ cpp/src/prima/bounds.hpp | 88 ++++ cpp/src/prima/checkbreak.hpp | 61 +++ cpp/src/prima/cobyla/cobyla.hpp | 536 +++++++++++++++++++++++ cpp/src/prima/cobyla/cobylb.hpp | 553 ++++++++++++++++++++++++ cpp/src/prima/cobyla/geometry.hpp | 225 ++++++++++ cpp/src/prima/cobyla/initialize.hpp | 149 +++++++ cpp/src/prima/cobyla/trustregion.hpp | 546 +++++++++++++++++++++++ cpp/src/prima/cobyla/update.hpp | 236 ++++++++++ cpp/src/prima/consts.hpp | 57 +++ cpp/src/prima/evaluate.hpp | 106 +++++ cpp/src/prima/history.hpp | 42 ++ cpp/src/prima/infos.hpp | 36 ++ cpp/src/prima/linalg.hpp | 404 +++++++++++++++++ cpp/src/prima/linear_constraints.hpp | 161 +++++++ cpp/src/prima/message.hpp | 126 ++++++ cpp/src/prima/nonlinear_constraints.hpp | 137 ++++++ cpp/src/prima/powalg.hpp | 162 +++++++ cpp/src/prima/preproc.hpp | 240 ++++++++++ cpp/src/prima/present.hpp | 22 + cpp/src/prima/prima.hpp | 317 ++++++++++++++ cpp/src/prima/project.hpp | 138 ++++++ cpp/src/prima/ratio.hpp | 44 ++ cpp/src/prima/redrho.hpp | 30 ++ cpp/src/prima/selectx.hpp | 237 ++++++++++ cpp/src/primacpp-config.cmake.in | 6 + cpp/tests/CMakeLists.txt | 20 + cpp/tests/test_bounds.cpp | 137 ++++++ cpp/tests/test_common.hpp | 34 ++ cpp/tests/test_end_to_end.cpp | 99 +++++ cpp/tests/test_get_lincon.cpp | 103 +++++ cpp/tests/test_miscellaneous.cpp | 168 +++++++ 39 files changed, 5595 insertions(+) create mode 100644 .github/workflows/test_cpp.yml create mode 100644 cpp/CMakeLists.txt create mode 100644 cpp/README.md create mode 100644 cpp/examples/CMakeLists.txt create mode 100644 cpp/examples/cobyla_example.cpp create mode 100644 cpp/examples/rosenbrock.cpp create mode 100644 cpp/src/CMakeLists.txt create mode 100644 cpp/src/prima/bounds.hpp create mode 100644 cpp/src/prima/checkbreak.hpp create mode 100644 cpp/src/prima/cobyla/cobyla.hpp create mode 100644 cpp/src/prima/cobyla/cobylb.hpp create mode 100644 cpp/src/prima/cobyla/geometry.hpp create mode 100644 cpp/src/prima/cobyla/initialize.hpp create mode 100644 cpp/src/prima/cobyla/trustregion.hpp create mode 100644 cpp/src/prima/cobyla/update.hpp create mode 100644 cpp/src/prima/consts.hpp create mode 100644 cpp/src/prima/evaluate.hpp create mode 100644 cpp/src/prima/history.hpp create mode 100644 cpp/src/prima/infos.hpp create mode 100644 cpp/src/prima/linalg.hpp create mode 100644 cpp/src/prima/linear_constraints.hpp create mode 100644 cpp/src/prima/message.hpp create mode 100644 cpp/src/prima/nonlinear_constraints.hpp create mode 100644 cpp/src/prima/powalg.hpp create mode 100644 cpp/src/prima/preproc.hpp create mode 100644 cpp/src/prima/present.hpp create mode 100644 cpp/src/prima/prima.hpp create mode 100644 cpp/src/prima/project.hpp create mode 100644 cpp/src/prima/ratio.hpp create mode 100644 cpp/src/prima/redrho.hpp create mode 100644 cpp/src/prima/selectx.hpp create mode 100644 cpp/src/primacpp-config.cmake.in create mode 100644 cpp/tests/CMakeLists.txt create mode 100644 cpp/tests/test_bounds.cpp create mode 100644 cpp/tests/test_common.hpp create mode 100644 cpp/tests/test_end_to_end.cpp create mode 100644 cpp/tests/test_get_lincon.cpp create mode 100644 cpp/tests/test_miscellaneous.cpp diff --git a/.github/workflows/test_cpp.yml b/.github/workflows/test_cpp.yml new file mode 100644 index 0000000000..f81324b6f4 --- /dev/null +++ b/.github/workflows/test_cpp.yml @@ -0,0 +1,43 @@ +name: Test CPP + +on: + pull_request: + schedule: + - cron: '0 16 4-31/4 * *' + workflow_dispatch: + +jobs: + test: + name: Run CPP tests + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, ubuntu-24.04-arm, windows-latest, macos-latest] + steps: + - uses: actions/checkout@v6 + + - name: Install dependencies (Ubuntu) + if: runner.os == 'Linux' + run: sudo apt-get update && sudo apt-get install -y libeigen3-dev + + - name: Install dependencies (macOS) + if: runner.os == 'macOS' + run: brew install eigen + + - name: Install dependencies (Windows) + if: runner.os == 'Windows' + run: vcpkg install eigen3 --triplet x64-windows + + - name: Set vcpkg toolchain (Windows) + if: runner.os == 'Windows' + shell: pwsh + run: echo "CMAKE_TOOLCHAIN_FILE=$env:VCPKG_INSTALLATION_ROOT/scripts/buildsystems/vcpkg.cmake" >> $env:GITHUB_ENV + + - name: Build + run: | + cmake -DPRIMA_ENABLE_TESTING=ON -DCMAKE_CXX_STANDARD=17 -DCMAKE_CXX_STANDARD_REQUIRED=TRUE -DCMAKE_BUILD_TYPE=Release -B build -S cpp + cmake --build build -j3 --config Release + + - name: Test + run: ctest --test-dir build --output-on-failure --timeout 100 --schedule-random -V -j3 -C Release diff --git a/CMakeLists.txt b/CMakeLists.txt index 87d008e3bf..0e874b3848 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -129,6 +129,12 @@ if (PRIMA_ENABLE_C) set(primac_target "primac") endif () +option (PRIMA_ENABLE_CPP "C++ binding" OFF) +if (PRIMA_ENABLE_CPP) + enable_language(CXX) + add_subdirectory(cpp) +endif () + # Get the version number find_package(Git) set(IS_REPO FALSE) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt new file mode 100644 index 0000000000..5f13f26884 --- /dev/null +++ b/cpp/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required (VERSION 3.18) + +project (primacpp VERSION 0.1.0 LANGUAGES CXX) + +option (BUILD_SHARED_LIBS "shared/static" ON) + +include (GNUInstallDirs) + +find_package (Eigen3 REQUIRED) +message (STATUS "Found Eigen3: ${Eigen3_DIR} (found version ${Eigen3_VERSION})") + +add_subdirectory (src) +enable_testing() +add_subdirectory (tests) +add_subdirectory (examples) diff --git a/cpp/README.md b/cpp/README.md new file mode 100644 index 0000000000..072bdfd2e7 --- /dev/null +++ b/cpp/README.md @@ -0,0 +1,46 @@ +## About + +This is a C++ translation of [Zaikun Zhang](https://www.zhangzk.net)'s [modern-Fortran reference implementation](https://github.com/libprima/prima/tree/main/fortran) +for Powell's derivative-free optimization solvers, which is available at `fortran/` under the root directory. +It is a faithful translation of the [Python translation](https://github.com/libprima/prima/tree/main/pyprima), +following the same structure, variable names, and algorithm logic to keep maintenance across languages tractable. + +Due to [bug-fixes](https://github.com/libprima/prima#bug-fixes) and [improvements](https://github.com/libprima/prima#improvements), +the modern-Fortran reference implementation by [Zaikun Zhang](https://www.zhangzk.net) +behaves differently from the original Fortran 77 implementation by [M. J. D. Powell](https://www.zhangzk.net/powell.html), +even though the algorithms are essentially the same. Therefore, it is important to point out that you are using +PRIMA rather than the original solvers if you want your results to be reproducible. + +As of June 2026, only the COBYLA solver is available in this C++ translation. +The other solvers will be translated from the Python/Fortran reference implementations in the future. + +## Building + +This is a header-only library requiring only Eigen3. To build the tests: + +```bash +cmake -S cpp -B build -DEigen3_DIR=/path/to/eigen3/cmake +cmake --build build --target test_minimize_cpp_exe +``` + +To install: + +```bash +cmake --install build --prefix /usr/local +``` + +After installation, use from another project: + +```cmake +find_package(primacpp REQUIRED) +target_link_libraries(myapp PRIVATE prima::primacpp) +``` + +## Development notes + +- Function names, variable names, and file layout follow the Fortran and Python implementations. + Keep them in sync when making changes. +- Comments are kept minimal compared to the Python/Fortran sources. When the intent is unclear, + refer to the Python or Fortran reference. +- The library is header-only. All implementation is in `.hpp` files under `src/prima/`. +- The namespace is `prima`; internals go in `prima::detail`. diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt new file mode 100644 index 0000000000..466ee407b0 --- /dev/null +++ b/cpp/examples/CMakeLists.txt @@ -0,0 +1,5 @@ +add_executable (cobyla_example_exe EXCLUDE_FROM_ALL cobyla_example.cpp) +add_executable (rosenbrock_example_exe EXCLUDE_FROM_ALL rosenbrock.cpp) + +target_link_libraries (cobyla_example_exe PRIVATE primacpp) +target_link_libraries (rosenbrock_example_exe PRIVATE primacpp) diff --git a/cpp/examples/cobyla_example.cpp b/cpp/examples/cobyla_example.cpp new file mode 100644 index 0000000000..2f35f5cd9b --- /dev/null +++ b/cpp/examples/cobyla_example.cpp @@ -0,0 +1,64 @@ +// This is an example to illustrate the usage of the COBYLA solver. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include +#include "prima/prima.hpp" + +using namespace prima; +using namespace Eigen; + +// Objective: f(x) = (x1-5)^2 + (x2-4)^2 +double objective(const VectorXd& x) { + return std::pow(x(0) - 5.0, 2) + std::pow(x(1) - 4.0, 2); +} + +int main() { + std::cout << "=== COBYLA Example ===" << std::endl; + + // Simple constrained optimization: + // min (x1-5)^2 + (x2-4)^2 + // s.t. x1^2 - 9 <= 0 (i.e., |x1| <= 3) + + VectorXd x0(2); + x0 << 0.0, 0.0; + + // Nonlinear constraint: x1^2 - 9 <= 0 + auto cons_fun = [](const VectorXd& x) -> VectorXd { + return VectorXd::Constant(1, x(0) * x(0) - 9.0); + }; + VectorXd nlc_lb(1); + nlc_lb << -std::numeric_limits::infinity(); + VectorXd nlc_ub(1); + nlc_ub << 0.0; + NonlinearConstraint nlc(cons_fun, nlc_lb, nlc_ub); + auto nlc_func = transform_constraint_function(nlc); + + MinimizeOptions opts; + opts.quiet = false; // Print progress + opts.rhoend = 1e-6; + opts.maxfun = 5000; + + auto result = minimize(objective, x0, "cobyla", nullptr, nullptr, &nlc_func, opts); + + std::cout << "\nResult:" << std::endl; + std::cout << " x = [" << result.x(0) << ", " << result.x(1) << "]" << std::endl; + std::cout << " f = " << result.fun << std::endl; + std::cout << " constraint x1^2-9 = " << (result.x(0) * result.x(0) - 9.0) << std::endl; + std::cout << " nfev = " << result.nfev << std::endl; + + // Check: x1 should be near 3, x2 near 4, f near 4 + if (std::abs(result.x(0) - 3.0) < 1e-2 && + std::abs(result.x(1) - 4.0) < 1e-2 && + std::abs(result.fun - 4.0) < 1e-2) { + std::cout << "\nExample PASSED." << std::endl; + return 0; + } else { + std::cerr << "\nExample FAILED." << std::endl; + return 1; + } +} diff --git a/cpp/examples/rosenbrock.cpp b/cpp/examples/rosenbrock.cpp new file mode 100644 index 0000000000..ed68979f1c --- /dev/null +++ b/cpp/examples/rosenbrock.cpp @@ -0,0 +1,148 @@ +// Illustration of how to use prima with COBYLA. +// +// Minimize the chained Rosenbrock function subject to various constraints. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include +#include +#include "prima/prima.hpp" + +using namespace prima; +using namespace Eigen; + +// Chained Rosenbrock function +double chrosen(const VectorXd& x) { + int n = x.size(); + double f = 0; + for (int i = 0; i < n - 1; ++i) { + f += std::pow(1 - x(i), 2) + 4 * std::pow(x(i + 1) - x(i) * x(i), 2); + } + return f; +} + +// Nonlinear inequality constraint: x(i)^2 >= x(i+1) +VectorXd nlc_ineq(const VectorXd& x) { + int n = x.size(); + VectorXd c(n - 1); + for (int i = 0; i < n - 1; ++i) { + c(i) = x(i) * x(i) - x(i + 1); + } + return c; +} + +// Nonlinear equality constraint: ||x||^2 = 1 +VectorXd nlc_eq(const VectorXd& x) { + VectorXd c(1); + c(0) = x.squaredNorm() - 1; + return c; +} + +int main() { + std::cout << std::setprecision(4); + std::cout << "Minimize the chained Rosenbrock function with three variables " + << "subject to various constraints using COBYLA.\n" << std::endl; + + VectorXd x0(3); + x0 << 0, 0, 0; + + // ---------------------------------------------------------------- // + // 1. Nonlinear constraints + // ||x||_2^2 = 1, x(i)^2 >= x(i+1) >= 0.5*x(i) >= 0 for i = 1, 2 + // ---------------------------------------------------------------- // + std::cout << "1. Nonlinear constraints --- ||x||_2^2 = 1, " + << "x(i)^2 >= x(i+1) >= 0.5*x(i) >= 0 for i = 1, 2:\n" << std::endl; + + VectorXd lb(3); lb << 0, 0, 0; + VectorXd ub(3); ub << std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::numeric_limits::infinity(); + Bounds bounds(lb, ub); + + // Linear constraints: 0.5*x(i) - x(i+1) <= 0 + MatrixXd A(2, 3); + A << 0.5, -1, 0, + 0, 0.5, -1; + LinearConstraint lin_con(A, + VectorXd::Constant(2, -std::numeric_limits::infinity()), + VectorXd::Constant(2, 0.0)); + + // Nonlinear constraints + NonlinearConstraint nlc_ineq_obj(nlc_ineq, + VectorXd::Constant(2, 0.0), + VectorXd::Constant(2, std::numeric_limits::infinity())); + NonlinearConstraint nlc_eq_obj(nlc_eq, + VectorXd::Constant(1, 0.0), + VectorXd::Constant(1, 0.0)); + + auto nlc_ineq_t = transform_constraint_function(nlc_ineq_obj); + auto nlc_eq_t = transform_constraint_function(nlc_eq_obj); + + MinimizeOptions opts; + opts.quiet = true; + + // COBYLA only handles one NonlinearConstraintFunction, so combine them + NonlinearConstraintFunction combined_nlc = [nlc_ineq_t, nlc_eq_t](const VectorXd& x) -> VectorXd { + VectorXd v1 = nlc_ineq_t(x); + VectorXd v2 = nlc_eq_t(x); + VectorXd r(v1.size() + v2.size()); + r.head(v1.size()) = v1; + r.tail(v2.size()) = v2; + return r; + }; + + auto result = minimize(chrosen, x0, "cobyla", &bounds, &lin_con, &combined_nlc, opts); + std::cout << " x = [" << result.x(0) << ", " << result.x(1) << ", " << result.x(2) << "]" << std::endl; + std::cout << " f = " << result.fun << " nfev = " << result.nfev << std::endl; + std::cout << " ||x||^2 = " << result.x.squaredNorm() << std::endl; + + // ---------------------------------------------------------------- // + // 2. Linear constraints + // sum(x) = 1, x(i+1) <= x(i) <= 1 for i = 1, 2 + // ---------------------------------------------------------------- // + std::cout << "\n2. Linear constraints --- sum(x) = 1, x(i+1) <= x(i) <= 1 for i = 1, 2:\n" << std::endl; + + Bounds bounds2( + VectorXd::Constant(3, -std::numeric_limits::infinity()), + VectorXd::Constant(3, 1.0)); + MatrixXd A2(3, 3); + A2 << -1, 1, 0, + 0, -1, 1, + 1, 1, 1; + LinearConstraint lin_con2(A2, + Vector3d(-std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), 1.0), + Vector3d(0.0, 0.0, 1.0)); + + auto result2 = minimize(chrosen, x0, "cobyla", &bounds2, &lin_con2, nullptr, opts); + std::cout << " x = [" << result2.x(0) << ", " << result2.x(1) << ", " << result2.x(2) << "]" << std::endl; + std::cout << " f = " << result2.fun << " nfev = " << result2.nfev << std::endl; + std::cout << " sum(x) = " << result2.x.sum() << std::endl; + + // ---------------------------------------------------------------- // + // 3. Bound constraints: -0.5 <= x(1) <= 0.5, 0 <= x(2) <= 0.25 + // ---------------------------------------------------------------- // + std::cout << "\n3. Bound constraints --- -0.5 <= x(1) <= 0.5, 0 <= x(2) <= 0.25:\n" << std::endl; + + Bounds bounds3( + Vector3d(-0.5, 0.0, -std::numeric_limits::infinity()), + Vector3d(0.5, 0.25, std::numeric_limits::infinity())); + + auto result3 = minimize(chrosen, x0, "cobyla", &bounds3, nullptr, nullptr, opts); + std::cout << " x = [" << result3.x(0) << ", " << result3.x(1) << ", " << result3.x(2) << "]" << std::endl; + std::cout << " f = " << result3.fun << " nfev = " << result3.nfev << std::endl; + + // ---------------------------------------------------------------- // + // 4. No constraints + // ---------------------------------------------------------------- // + std::cout << "\n4. No constraints:\n" << std::endl; + auto result4 = minimize(chrosen, x0, "cobyla", nullptr, nullptr, nullptr, opts); + std::cout << " x = [" << result4.x(0) << ", " << result4.x(1) << ", " << result4.x(2) << "]" << std::endl; + std::cout << " f = " << result4.fun << " nfev = " << result4.nfev << std::endl; + + return 0; +} diff --git a/cpp/src/CMakeLists.txt b/cpp/src/CMakeLists.txt new file mode 100644 index 0000000000..760855923a --- /dev/null +++ b/cpp/src/CMakeLists.txt @@ -0,0 +1,48 @@ +add_library (primacpp INTERFACE) +target_include_directories (primacpp INTERFACE + $ + $ +) +target_link_libraries (primacpp INTERFACE Eigen3::Eigen) +target_compile_options (primacpp INTERFACE $<$:/bigobj>) +target_compile_features (primacpp INTERFACE cxx_std_17) + +include (GNUInstallDirs) +include (CMakePackageConfigHelpers) + +install ( + TARGETS primacpp + EXPORT primacpp-targets +) + +install ( + DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/prima + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} + FILES_MATCHING PATTERN "*.hpp" +) + +install ( + EXPORT primacpp-targets + FILE primacpp-targets.cmake + NAMESPACE prima:: + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/primacpp +) + +configure_package_config_file ( + ${CMAKE_CURRENT_SOURCE_DIR}/primacpp-config.cmake.in + ${CMAKE_BINARY_DIR}/primacpp-config.cmake + INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/primacpp +) + +write_basic_package_version_file ( + ${CMAKE_BINARY_DIR}/primacpp-config-version.cmake + VERSION ${PROJECT_VERSION} + COMPATIBILITY AnyNewerVersion +) + +install ( + FILES + ${CMAKE_BINARY_DIR}/primacpp-config.cmake + ${CMAKE_BINARY_DIR}/primacpp-config-version.cmake + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/primacpp +) diff --git a/cpp/src/prima/bounds.hpp b/cpp/src/prima/bounds.hpp new file mode 100644 index 0000000000..5b7d6f48c6 --- /dev/null +++ b/cpp/src/prima/bounds.hpp @@ -0,0 +1,88 @@ +#ifndef PRIMA_CPP_BOUNDS_HPP +#define PRIMA_CPP_BOUNDS_HPP + +/* +This module handles bound constraints for optimization variables. + +bounds can either be an object with the properties lb and ub, or a list of tuples +indicating a lower bound and an upper bound for each variable. If the list contains +fewer entries than the length of x0, the remaining entries will be generated as -/+ infinity. + +Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. + +Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +*/ + +#include +#include +#include +#include + +namespace prima { + +// Container for lower and upper bounds +struct Bounds { + Eigen::VectorXd lb; + Eigen::VectorXd ub; + + Bounds() = default; + + Bounds(const Eigen::VectorXd& lb_, const Eigen::VectorXd& ub_) + : lb(lb_), ub(ub_) {} +}; + +/* +Process bounds for optimization: +- If bounds is nullptr, returns [-inf, inf] for all variables. +- Truncates bounds if longer than lenx0. +- Pads bounds with -/+inf if shorter than lenx0. +- Checks for infeasible bounds (lb > ub). + +Returns a pair of (lb, ub) vectors of length lenx0. +*/ +inline std::pair +process_bounds(const Bounds* bounds, int lenx0) { + using namespace Eigen; + + if (bounds == nullptr) { + VectorXd lb = VectorXd::Constant(lenx0, -std::numeric_limits::infinity()); + VectorXd ub = VectorXd::Constant(lenx0, std::numeric_limits::infinity()); + return {lb, ub}; + } + + VectorXd lb = bounds->lb; + VectorXd ub = bounds->ub; + + // Truncate if bounds are longer than lenx0 + if (lb.size() > lenx0) lb = lb.head(lenx0); + if (ub.size() > lenx0) ub = ub.head(lenx0); + + // Pad if bounds are shorter than lenx0 + if (lb.size() < lenx0) { + VectorXd pad = VectorXd::Constant(lenx0 - lb.size(), -std::numeric_limits::infinity()); + VectorXd tmp(lb.size() + pad.size()); + tmp << lb, pad; + lb = tmp; + } + if (ub.size() < lenx0) { + VectorXd pad = VectorXd::Constant(lenx0 - ub.size(), std::numeric_limits::infinity()); + VectorXd tmp(ub.size() + pad.size()); + tmp << ub, pad; + ub = tmp; + } + + // Check the infeasibility of the bounds + for (int i = 0; i < lenx0; ++i) { + if (lb(i) > ub(i)) { + throw std::invalid_argument( + "Some of the provided bounds are infeasible." + ); + } + } + + return {lb, ub}; +} + +} // namespace prima + +#endif // PRIMA_CPP_BOUNDS_HPP diff --git a/cpp/src/prima/checkbreak.hpp b/cpp/src/prima/checkbreak.hpp new file mode 100644 index 0000000000..1d23219d28 --- /dev/null +++ b/cpp/src/prima/checkbreak.hpp @@ -0,0 +1,61 @@ +#ifndef PRIMA_CPP_CHECKBREAK_HPP +#define PRIMA_CPP_CHECKBREAK_HPP + +// This module checks whether to break out of the solver loop. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include + +#include "infos.hpp" +#include "consts.hpp" + +namespace prima { + +// This function checks whether to break out of the solver loop in the unconstrained case. +inline int checkbreak_unc(int maxfun, int nf, double f, double ftarget, const Eigen::VectorXd& x) { + int info = INFO_DEFAULT; + + // Although X should not contain NaN unless there is a bug, we include the following for security. + // X can be Inf, as finite + finite can be Inf numerically. + for (Eigen::Index i = 0; i < x.size(); ++i) { + if (std::isnan(x[i]) || std::isinf(x[i])) { info = NAN_INF_X; break; } + } + + // Although NAN_INF_F should not happen unless there is a bug, we include the following for security. + if (std::isnan(f) || std::isinf(f)) info = NAN_INF_F; + + if (f <= ftarget) info = FTARGET_ACHIEVED; + if (nf >= maxfun) info = MAXFUN_REACHED; + + return info; +} + +// This function checks whether to break out of the solver loop in the constrained case. +inline int checkbreak_con(int maxfun, int nf, double cstrv, double ctol, + double f, double ftarget, const Eigen::VectorXd& x) { + int info = INFO_DEFAULT; + + // Although X should not contain NaN unless there is a bug, we include the following for security. + // X can be Inf, as finite + finite can be Inf numerically. + for (Eigen::Index i = 0; i < x.size(); ++i) { + if (std::isnan(x[i]) || std::isinf(x[i])) { info = NAN_INF_X; break; } + } + + // Although NAN_INF_F should not happen unless there is a bug, we include the following for security. + if (std::isnan(f) || std::isinf(f) || std::isnan(cstrv) || std::isinf(cstrv)) + info = NAN_INF_F; + + if (cstrv <= ctol && f <= ftarget) info = FTARGET_ACHIEVED; + if (nf >= maxfun) info = MAXFUN_REACHED; + + return info; +} + +} // namespace prima + +#endif // PRIMA_CPP_CHECKBREAK_HPP diff --git a/cpp/src/prima/cobyla/cobyla.hpp b/cpp/src/prima/cobyla/cobyla.hpp new file mode 100644 index 0000000000..ceb181b4dd --- /dev/null +++ b/cpp/src/prima/cobyla/cobyla.hpp @@ -0,0 +1,536 @@ +#ifndef PRIMA_CPP_COBYLA_COBYLA_HPP +#define PRIMA_CPP_COBYLA_COBYLA_HPP + +// This module provides Powell's COBYLA algorithm. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +// +// N.B.: +// +// 1. The modern-Fortran reference implementation in PRIMA contains bug fixes and improvements over the +// original Fortran 77 implementation by Powell. Consequently, the PRIMA implementation behaves differently +// from the original Fortran 77 implementation by Powell. Therefore, it is important to point out that +// you are using PRIMA rather than the original solvers if you want your results to be reproducible. +// +// 2. Compared to Powell's Fortran 77 implementation, the modern-Fortran implementation and hence any +// faithful translation like this one generally produce better solutions with fewer function evaluations, +// making them preferable for applications with expensive function evaluations. However, if function +// evaluations are not the dominant cost in your application, the Fortran 77 solvers are likely to be +// faster, as they are more efficient in terms of memory usage and flops thanks to the careful and +// ingenious (but unmaintained and unmaintainable) implementation by Powell. +// +// See the PRIMA documentation (www.libprima.net) for more information. + +#include +#include +#include +#include +#include +#include +#include + +#include "../consts.hpp" +#include "../evaluate.hpp" +#include "../linalg.hpp" +#include "../preproc.hpp" +#include "../present.hpp" + +#include "cobylb.hpp" + +namespace prima { + +struct COBYLAResult { + Eigen::VectorXd x; + double f; + Eigen::VectorXd constr; + double cstrv; + int nf; + std::vector xhist; + std::vector fhist; + std::vector chist; + std::vector conhist; + int info; +}; + +// This subroutine wraps the linear and bound constraints into a single constraint: +// AMAT*X <= BVEC. +// +// N.B.: +// +// LINCOA normalizes the linear constraints so that each constraint has a gradient +// of norm 1. However, COBYLA does not do this. +inline std::tuple, std::optional> +get_lincon(const std::optional& Aeq, + const std::optional& Aineq, + const std::optional& beq, + const std::optional& bineq, + const std::optional& xl, + const std::optional& xu) { + // Sizes + int num_vars = 0; + if (Aeq.has_value()) num_vars = Aeq->cols(); + else if (Aineq.has_value()) num_vars = Aineq->cols(); + else if (xl.has_value()) num_vars = xl->size(); + else if (xu.has_value()) num_vars = xu->size(); + else return {std::nullopt, std::nullopt}; + + //====================// + // Calculation starts // + //====================// + + // Define the indices of the nontrivial bound constraints. + std::vector ixl; + std::vector ixu; + + if (xl.has_value()) { + for (int i = 0; i < xl->size(); ++i) + if ((*xl)[i] > -BOUNDMAX) ixl.push_back(i); + } + if (xu.has_value()) { + for (int i = 0; i < xu->size(); ++i) + if ((*xu)[i] < BOUNDMAX) ixu.push_back(i); + } + + // Wrap the linear constraints. + // The bound constraint XL <= X <= XU is handled as two constraints: + // -X <= -XL, X <= XU. + // The equality constraint Aeq*X = Beq is handled as two constraints: + // -Aeq*X <= -Beq, Aeq*X <= Beq. + // N.B.: + // 1. The treatment of the equality constraints is naive. One may choose to + // eliminate them instead. + Eigen::MatrixXd idmat = Eigen::MatrixXd::Identity(num_vars, num_vars); + + int mxl = ixl.size(); + int mxu = ixu.size(); + int meq = Aeq.has_value() ? Aeq->rows() : 0; + int mineq = Aineq.has_value() ? Aineq->rows() : 0; + int total_rows = mxl + mxu + 2 * meq + mineq; + + if (total_rows == 0) return {std::nullopt, std::nullopt}; + + Eigen::MatrixXd amat(total_rows, num_vars); + Eigen::VectorXd bvec(total_rows); + + int row = 0; + for (int i : ixl) { + amat.row(row) = -idmat.row(i); + bvec[row] = -(*xl)[i]; + row++; + } + for (int i : ixu) { + amat.row(row) = idmat.row(i); + bvec[row] = (*xu)[i]; + row++; + } + if (Aeq.has_value()) { + for (int i = 0; i < Aeq->rows(); ++i) { + amat.row(row) = -Aeq->row(i); + bvec[row] = -(*beq)[i]; + row++; + } + for (int i = 0; i < Aeq->rows(); ++i) { + amat.row(row) = Aeq->row(i); + bvec[row] = (*beq)[i]; + row++; + } + } + if (Aineq.has_value()) { + for (int i = 0; i < Aineq->rows(); ++i) { + amat.row(row) = Aineq->row(i); + bvec[row] = (*bineq)[i]; + row++; + } + } + + //==================// + // Calculation ends // + //==================// + + return {amat, bvec}; +} + +// Among all the arguments, only CALCFC, M_NLCON, and X are obligatory. The others are +// OPTIONAL and you can neglect them unless you are familiar with the algorithm. Any +// unspecified optional input will take the default value detailed below. For +// instance, we may invoke the solver as follows. +// +// // First define CALCFC, M_NLCON, and X, and then do the following. +// result = cobyla(calcfc, m_nlcon, x); +// +// or +// +// // First define CALCFC, M_NLCON, X, Aineq, and Bineq, and then do the following. +// result = cobyla(calcfc, m_nlcon, x, Aineq, bineq, rhobeg=1.0e0, +// rhoend=1.0e-6); +// +// #################################################################################### +// # IMPORTANT NOTICE: The user must set M_NLCON correctly to the number of nonlinear +// # constraints, namely the size of NLCONSTR introduced below. Set it to 0 if there +// # is no nonlinear constraint. +// #################################################################################### +// +// See examples/cobyla/cobyla_example.py for a concrete example. +// +// A detailed introduction to the arguments is as follows. +// +// #################################################################################### +// # INPUTS +// #################################################################################### +// +// CALCFC +// Input, function. +// f, nlconstr = CALCFC(X) should evaluate the objective function and nonlinear +// constraints at the given vector X; it should return a tuple consisting of the +// objective function value and the nonlinear constraint value. It must be provided +// by the user, and its definition must conform to the following interface: +// //-------------------------------------------------------------------------// +// auto [f, nlconstr] = calcfc(x); +// //-------------------------------------------------------------------------// +// +// M_NLCON +// Input, scalar. +// M_NLCON must be set to the number of nonlinear constraints, namely the size of +// NLCONSTR(X). +// N.B.: +// 1. Why don't we define M_NLCON as optional and default it to 0 when it is absent? +// This is because we need to allocate memory for CONSTR_LOC using M_NLCON. To +// ensure that the size of CONSTR_LOC is correct, we require the user to specify +// M_NLCON explicitly. +// +// X +// Input, vector. +// As an input, X should be an N-dimensional vector that contains the starting +// point, N being the dimension of the problem. +// +// Aineq, Bineq +// Input, matrix of size [Mineq, N] and vector of size Mineq unless they are both +// empty, default: None and None. +// Aineq and Bineq represent the linear inequality constraints: Aineq*X <= Bineq. +// +// Aeq, Beq +// Input, matrix of size [Meq, N] and vector of size Meq unless they are both +// empty, default: None and None. +// Aeq and Beq represent the linear equality constraints: Aeq*X = Beq. +// +// XL, XU +// Input, vectors of size N unless they are both None, default: None and None. +// XL is the lower bound for X. If XL is None, X has no +// lower bound. Any entry of XL that is NaN or below -BOUNDMAX will be taken as +// -BOUNDMAX, which effectively means there is no lower bound for the corresponding +// entry of X. The value of BOUNDMAX is 0.25*HUGE(X), which is about 8.6E37 for +// single precision and 4.5E307 for double precision. XU is similar. +// +// F0 +// Input, scalar. +// F0, if present, should be set to the objective function value of the starting X. +// +// NLCONSTR0 +// Input, vector. +// NLCONSTR0, if present, should be set to the nonlinear constraint value at the +// starting X; in addition, SIZE(NLCONSTR0) must be M_NLCON, or the solver will +// abort. +// +// RHOBEG, RHOEND +// Inputs, scalars, default: RHOBEG = 1, RHOEND = 10^-6. RHOBEG and RHOEND must be +// set to the initial and final values of a trust-region radius, both being positive +// and RHOEND <= RHOBEG. Typically RHOBEG should be about one tenth of the greatest +// expected change to a variable, and RHOEND should indicate the accuracy that is +// required in the final values of the variables. +// +// FTARGET +// Input, scalar, default: -Inf. +// FTARGET is the target function value. The algorithm will terminate when a +// feasible point with a function value <= FTARGET is found. +// +// CTOL +// Input, scalar, default: sqrt(machine epsilon). +// CTOL is the tolerance of constraint violation. X is considered feasible if +// CSTRV(X) <= CTOL. +// N.B.: +// 1. CTOL is absolute, not relative. +// 2. CTOL is used only when selecting the returned X. It does not affect the +// iterations of the algorithm. +// +// CWEIGHT +// Input, scalar, default: CWEIGHT_DFT defined in common/consts.hpp. +// CWEIGHT is the weight that the constraint violation takes in the selection of the +// returned X. +// +// MAXFUN +// Input, integer scalar, default: MAXFUN_DIM_DFT*N with MAXFUN_DIM_DFT defined in +// common/consts.hpp. MAXFUN is the maximal number of calls of CALCFC. +// +// IPRINT +// Input, integer scalar, default: 0. +// The value of IPRINT should be set to 0, 1, -1, 2, -2, 3, or -3, which controls +// how much information will be printed during the computation: +// 0: there will be no printing; +// 1: a message will be printed to the screen at the return, showing the best vector +// of variables found and its objective function value; +// 2: in addition to 1, each new value of RHO is printed to the screen, with the +// best vector of variables so far and its objective function value; each new +// value of CPEN is also printed; +// 3: in addition to 2, each function evaluation with its variables will be printed +// to the screen; -1, -2, -3: the same information as 1, 2, 3 will be printed, +// not to the screen but to a file named COBYLA_output.txt; the file will be +// created if it does not exist; the new output will be appended to the end of +// this file if it already exists. +// Note that IPRINT = +/-3 can be costly in terms of time and/or space. +// +// ETA1, ETA2, GAMMA1, GAMMA2 +// Input, scalars, default: ETA1 = 0.1, ETA2 = 0.7, GAMMA1 = 0.5, and GAMMA2 = 2. +// ETA1, ETA2, GAMMA1, and GAMMA2 are parameters in the updating scheme of the +// trust-region radius detailed in the subroutine TRRAD in trustregion.hpp. Roughly +// speaking, the trust-region radius is contracted by a factor of GAMMA1 when the +// reduction ratio is below ETA1, and enlarged by a factor of GAMMA2 when the +// reduction ratio is above ETA2. It is required that 0 < ETA1 <= ETA2 < 1 and +// 0 < GAMMA1 < 1 < GAMMA2. Normally, ETA1 <= 0.25. It is NOT advised to set +// ETA1 >= 0.5. +// +// MAXFILT +// Input, scalar. +// MAXFILT is a nonnegative integer indicating the maximal length of the filter used +// for selecting the returned solution; default: MAXFILT_DFT (a value lower than +// MIN_MAXFILT is not recommended); +// see common/consts.hpp for the definitions of MAXFILT_DFT and MIN_MAXFILT. +// +// CALLBACK +// Input, function to report progress and optionally request termination. +// +// +// #################################################################################### +// # OUTPUTS +// #################################################################################### +// +// The output is a single data structure, COBYLAResult, with the following fields: +// +// X +// Output, vector. +// As an output, X will be set to an approximate minimizer. +// +// F +// Output, scalar. +// F will be set to the objective function value of X at exit. +// +// CONSTR +// Output, vector. +// CONSTR will be set to the constraint value of X at exit. +// +// CSTRV +// Output, scalar. +// CSTRV will be set to the constraint violation of X at exit, i.e., +// max([0, XL - X, X - XU, Aineq*X - Bineq, ABS(Aeq*X - Beq), NLCONSTR(X)]). +// +// NF +// Output, scalar. +// NF will be set to the number of calls of CALCFC at exit. +// +// XHIST, FHIST, CHIST, CONHIST, MAXHIST +// XHIST: Output, rank 2 array; +// FHIST: Output, rank 1 array; +// CHIST: Output, rank 1 array; +// CONHIST: Output, rank 2 array; +// MAXHIST: Input, scalar, default: MAXFUN +// XHIST, if present, will output the history of iterates; FHIST, if present, will +// output the history function values; CHIST, if present, will output the history of +// constraint violations; CONHIST, if present, will output the history of constraint +// values; MAXHIST should be a nonnegative integer, and XHIST/FHIST/CHIST/CONHIST +// will output only the history of the last MAXHIST iterations. +// Therefore, MAXHIST= 0 means XHIST/FHIST/CONHIST/CHIST will output +// nothing, while setting MAXHIST = MAXFUN requests XHIST/FHIST/CHIST/CONHIST to +// output all the history. +// +// IMPORTANT NOTICE: +// Setting MAXHIST to a large value can be costly in terms of memory for large +// problems. +// MAXHIST will be reset to a smaller value if the memory needed exceeds MAXHISTMEM +// defined in common/consts.hpp. +// Use *HIST with caution!!! (N.B.: the algorithm is NOT designed for large +// problems). +// +// INFO +// Output, scalar. +// INFO is the exit flag. It will be set to one of the following values defined in +// common/infos.hpp: +// SMALL_TR_RADIUS: the lower bound for the trust region radius is reached; +// FTARGET_ACHIEVED: the target function value is reached; +// MAXFUN_REACHED: the objective function has been evaluated MAXFUN times; +// MAXTR_REACHED: the trust region iteration has been performed MAXTR times (MAXTR = 2*MAXFUN); +// NAN_INF_X: NaN or Inf occurs in X; +// DAMAGING_ROUNDING: rounding errors are becoming damaging. +// //--------------------------------------------------------------------------// +// The following case(s) should NEVER occur unless there is a bug. +// NAN_INF_F: the objective function returns NaN or +Inf; +// NAN_INF_MODEL: NaN or Inf occurs in the model; +// TRSUBP_FAILED: a trust region step failed to reduce the model +// //--------------------------------------------------------------------------// +inline COBYLAResult +cobyla(const Calcfc& calcfc, int m_nlcon, const Eigen::VectorXd& x, + std::optional Aineq = std::nullopt, + std::optional bineq = std::nullopt, + std::optional Aeq = std::nullopt, + std::optional beq = std::nullopt, + std::optional xl = std::nullopt, + std::optional xu = std::nullopt, + std::optional f0 = std::nullopt, + std::optional nlconstr0 = std::nullopt, + std::optional rhobeg = std::nullopt, + std::optional rhoend = std::nullopt, + double ftarget = FTARGET_DEFAULT, + double ctol = CTOL_DEFAULT, + double cweight = CWEIGHT_DEFAULT, + std::optional maxfun = std::nullopt, + int iprint = IPRINT_DEFAULT, + std::optional eta1 = std::nullopt, + std::optional eta2 = std::nullopt, + double gamma1 = GAMMA1_DEFAULT, + double gamma2 = GAMMA2_DEFAULT, + std::optional maxhist = std::nullopt, + int maxfilt = 2000, + std::function callback = nullptr) { + + // Local variables + std::string solver = "COBYLA"; + int num_vars = x.size(); + + // Sizes + int mineq = Aineq.has_value() ? Aineq->rows() : 0; + int meq = Aeq.has_value() ? Aeq->rows() : 0; + int mxl = 0, mxu = 0; + if (xl.has_value()) + for (int i = 0; i < xl->size(); ++i) + if ((*xl)[i] > -BOUNDMAX) mxl++; + if (xu.has_value()) + for (int i = 0; i < xu->size(); ++i) + if ((*xu)[i] < BOUNDMAX) mxu++; + int mmm = mxl + mxu + 2 * meq + mineq + m_nlcon; + + // Wrap the linear and bound constraints into a single constraint: AMAT@X <= BVEC. + auto [amat_opt, bvec_opt] = get_lincon(Aeq, Aineq, beq, bineq, xl, xu); + + // Create constraint vector + Eigen::VectorXd constr = Eigen::VectorXd::Zero(mmm); + + // Set [F_LOC, CONSTR_LOC] to [F(X0), CONSTR(X0)] after evaluating the latter if + // needed. In this way, COBYLB only needs one interface. + // N.B.: Due to the preconditions above, there are two possibilities for F0 and + // NLCONSTR0. + // If NLCONSTR0 is present, then F0 must be present, and we assume that F(X0) = F0 + // even if F0 is NaN. + // If NLCONSTR0 is absent, then F0 must be either absent or NaN, both of which will + // be interpreted as F(X0) is not provided and we have to evaluate F(X0) and + // NLCONSTR(X0) now. + Eigen::VectorXd constr_loc; + double f_loc; + if (f0.has_value() && nlconstr0.has_value()) { + bool x_finite = true; + for (int i = 0; i < x.size(); ++i) + if (!std::isfinite(x[i])) { x_finite = false; break; } + if (x_finite) { + f_loc = moderatef(f0.value()); + if (amat_opt.has_value()) { + constr.head(mmm - m_nlcon) = moderatec(matprod(amat_opt.value(), x) - bvec_opt.value()); + } + constr.tail(m_nlcon) = moderatec(nlconstr0.value()); + } else { + Eigen::VectorXd xm = moderatex(x); + std::tie(f_loc, constr_loc) = evaluate(calcfc, xm, m_nlcon, + amat_opt.has_value() ? &amat_opt.value() : nullptr, + bvec_opt.has_value() ? &bvec_opt.value() : nullptr); + constr = constr_loc; + } + } else { + Eigen::VectorXd xm = moderatex(x); + std::tie(f_loc, constr_loc) = evaluate(calcfc, xm, m_nlcon, + amat_opt.has_value() ? &amat_opt.value() : nullptr, + bvec_opt.has_value() ? &bvec_opt.value() : nullptr); + constr = constr_loc; + } + // N.B.: Do NOT call FMSG, SAVEHIST, or SAVEFILT for the function/constraint evaluation at X0. + // They will be called during the initialization, which will read the function/constraint at X0. + double cstrv = 0; + for (int i = 0; i < constr.size(); ++i) + cstrv = std::max(cstrv, constr[i]); + + // If RHOBEG is present, use it; otherwise, RHOBEG takes the default value for + // RHOBEG, taking the value of RHOEND into account. Note that RHOEND is considered + // only if it is present and it is VALID (i.e., finite and positive). The other + // inputs are read similarly. + double rhobeg_loc; + if (rhobeg.has_value()) { + rhobeg_loc = rhobeg.value(); + } else if (rhoend.has_value() && std::isfinite(rhoend.value()) && rhoend.value() > 0) { + rhobeg_loc = std::max(10.0 * rhoend.value(), RHOBEG_DEFAULT); + } else { + rhobeg_loc = RHOBEG_DEFAULT; + } + + double rhoend_loc; + if (rhoend.has_value()) { + rhoend_loc = rhoend.value(); + } else if (rhobeg_loc > 0) { + rhoend_loc = std::max(EPS, std::min(RHOEND_DEFAULT / RHOBEG_DEFAULT * rhobeg_loc, RHOEND_DEFAULT)); + } else { + rhoend_loc = RHOEND_DEFAULT; + } + + int maxfun_loc = maxfun.has_value() ? maxfun.value() : MAXFUN_DIM_DEFAULT * num_vars; + + double eta1_loc; + if (eta1.has_value()) { + eta1_loc = eta1.value(); + } else if (eta2.has_value() && 0 < eta2.value() && eta2.value() < 1) { + eta1_loc = std::max(EPS, eta2.value() / 7.0); + } else { + eta1_loc = ETA1_DEFAULT; + } + + double eta2_loc; + if (eta2.has_value()) { + eta2_loc = eta2.value(); + } else if (0 < eta1_loc && eta1_loc < 1) { + eta2_loc = (eta1_loc + 2.0) / 3.0; + } else { + eta2_loc = ETA2_DEFAULT; + } + + int maxhist_loc = maxhist.has_value() ? maxhist.value() + : std::max({maxfun_loc, num_vars + 2, MAXFUN_DIM_DEFAULT * num_vars}); + + // Preprocess the inputs in case some of them are invalid. It does nothing if all + // inputs are valid. + PreprocResult pp = preproc(solver, num_vars, iprint, maxfun_loc, maxhist_loc, + ftarget, rhobeg_loc, rhoend_loc, + mmm, std::nullopt, maxfilt, ctol, cweight, + eta1_loc, eta2_loc, gamma1, gamma2, (mmm > 0)); + // npt is unused in COBYLA + // _x0 is unused in COBYLA + + // Further revise MAXHIST according to MAXHISTMEM, and allocate memory for the history. + // In MATLAB/Python/Julia/R implementation, we should simply set MAXHIST = MAXFUN and initialize + // CHIST = NaN(1, MAXFUN), CONHIST = NaN(M, MAXFUN), FHIST = NaN(1, MAXFUN), XHIST = NaN(N, MAXFUN) + // if they are requested; replace MAXFUN with 0 for the history that is not requested. + + // call cobylb, which performs the real calculations + auto [x_out, f_out, constr_out, cstrv_out, nf_out, + xhist_out, fhist_out, chist_out, conhist_out, info_out] = + cobylb(calcfc, pp.iprint, pp.maxfilt, pp.maxfun, + amat_opt.has_value() ? &amat_opt.value() : nullptr, + bvec_opt.has_value() ? &bvec_opt.value() : nullptr, + pp.ctol, pp.cweight, pp.eta1, pp.eta2, + pp.ftarget, pp.gamma1, pp.gamma2, + pp.rhobeg, pp.rhoend, + constr, f_loc, x, pp.maxhist, callback); + + Eigen::VectorXd final_constr = constr_out.tail(m_nlcon); + + return {x_out, f_out, final_constr, cstrv_out, nf_out, + xhist_out, fhist_out, chist_out, conhist_out, info_out}; +} + +} // namespace prima + +#endif // PRIMA_CPP_COBYLA_COBYLA_HPP diff --git a/cpp/src/prima/cobyla/cobylb.hpp b/cpp/src/prima/cobyla/cobylb.hpp new file mode 100644 index 0000000000..d4ef4536cb --- /dev/null +++ b/cpp/src/prima/cobyla/cobylb.hpp @@ -0,0 +1,553 @@ +#ifndef PRIMA_CPP_COBYLA_COBYLB_HPP +#define PRIMA_CPP_COBYLA_COBYLB_HPP + +// This module performs the major calculations of COBYLA. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +// +// N.B. (Zaikun 20220131): Powell's implementation of COBYLA uses RHO rather than +// DELTA as the trust-region radius, and RHO is never increased. DELTA does not exist +// in Powell's COBYLA code. Following the idea in Powell's other solvers (UOBYQA, ..., +// LINCOA), our code uses DELTA as the trust-region radius, while RHO works as a lower +// bound of DELTA and indicates the current resolution of the algorithm. DELTA is +// updated in a classical way subject to DELTA >= RHO, whereas RHO is updated as in +// Powell's COBYLA code and is never increased. The new implementation improves the +// performance of COBYLA. + +#include +#include +#include +#include +#include +#include + +#include "../checkbreak.hpp" +#include "../consts.hpp" +#include "../evaluate.hpp" +#include "../history.hpp" +#include "../infos.hpp" +#include "../linalg.hpp" +#include "../message.hpp" +#include "../ratio.hpp" +#include "../redrho.hpp" +#include "../selectx.hpp" +#include "geometry.hpp" +#include "initialize.hpp" +#include "trustregion.hpp" +#include "update.hpp" + +namespace prima { + +// FC RATIO: calculate the ratio between the "typical change" of F and that of CONSTR. +// See equations (12)-(13) in Section 3 of the COBYLA paper for the definition. +inline double fcratio(const Eigen::MatrixXd& conmat, const Eigen::VectorXd& fval) { + int m = conmat.rows(); + if (m == 0) return 0; + Eigen::VectorXd cmin = -conmat.rowwise().maxCoeff(); + Eigen::VectorXd cmax = -conmat.rowwise().minCoeff(); + double fmin = fval.minCoeff(); + double fmax = fval.maxCoeff(); + double r = 0; + if ((cmin.array() < 0.5 * cmax.array()).any() && fmin < fmax) { + Eigen::VectorXd mask = (cmin.array() < 0.5 * cmax.array()).select(Eigen::VectorXd::Ones(m), Eigen::VectorXd::Zero(m)); + // Powell mentioned the following alternative in section 4 of his COBYLA paper. + // According to a test on 20230610, it does not make much difference. + // denom = np.max(max(cmax, 0) - cmin, mask=(cmin < 0.5 * cmax)) + double denom = std::numeric_limits::infinity(); + for (int i = 0; i < m; ++i) { + if (mask[i] > 0) { + denom = std::min(denom, std::max(cmax[i], 0.0) - cmin[i]); + } + } + if (std::isfinite(denom) && denom > 0) + r = (fmax - fmin) / denom; + } + return r; +} + +// GET CPEN: get the penalty parameter CPEN so that PREREM = PREREF + CPEN * PREREC > 0. +// See the discussions around equation (9) of the COBYLA paper. +// +// Increase CPEN if necessary to ensure PREREM > 0. Branch back for the next loop +// if this change alters the optimal vertex of the current simplex. +// +// Note the following: +// 1. In each loop, CPEN is changed only if PREREC > 0 > PREREF, in which case +// PREREM is guaranteed positive after the update. Note that PREREC >= 0 and +// max(PREREC, PREREF) > 0 in theory. If this holds numerically as well then CPEN +// is not changed only if PREREC = 0 or PREREF >= 0, in which case PREREM is +// currently positive, explaining why CPEN needs no update. +// 2. Even without an upper bound for the loop counter, the loop can occur at most +// NUM_VARS+1 times. This is because the update of CPEN does not decrease CPEN, +// and hence it can make vertex J (J <= NUM_VARS) become the new optimal vertex +// only if CVAL[J] is less than CVAL[NUM_VARS], which can happen at most NUM_VARS +// times. See the paragraph below (9) in the COBYLA paper. After the "correct" +// optimal vertex is found, one more loop is needed to calculate CPEN, and hence +// the loop can occur at most NUM_VARS+1 times. +// Powell's code defines BARMU = -PREREF / PREREC, and CPEN is increased to +// 2*BARMU if and only if it is currently less than 1.5*BARMU, a very +// "Powellful" scheme. In our implementation, however, we set CPEN directly to +// the maximum between its current value and 2*BARMU while handling possible +// overflow. This simplifies the scheme without worsening the performance. +inline double getcpen(const Eigen::MatrixXd* amat, const Eigen::VectorXd* bvec, + Eigen::MatrixXd& conmat, double cpen, + Eigen::VectorXd& cval, double delta, + Eigen::VectorXd& fval, double rho, + Eigen::MatrixXd& sim, Eigen::MatrixXd& simi) { + int m_lcon = (bvec != nullptr) ? bvec->size() : 0; + int num_constraints = conmat.rows(); + int num_vars = sim.rows(); + int info = INFO_DEFAULT; + + // Work on local copies to avoid modifying caller state (see pyprima bug HS102) + Eigen::MatrixXd conmat_loc = conmat; + Eigen::VectorXd cval_loc = cval; + Eigen::VectorXd fval_loc = fval; + Eigen::MatrixXd sim_loc = sim; + Eigen::MatrixXd simi_loc = simi; + + Eigen::MatrixXd A(num_vars, num_constraints); + + for (int iter = 0; iter <= num_vars; ++iter) { + // Switch the best vertex of the current simplex to SIM[:, NUM_VARS] + std::tie(conmat_loc, cval_loc, fval_loc, sim_loc, simi_loc, info) = + updatepole(cpen, conmat_loc, cval_loc, fval_loc, sim_loc, simi_loc); + if (info == DAMAGING_ROUNDING) break; + + // Calculate the linear approximations to the objective and constraint functions + Eigen::VectorXd fdiff1 = (fval_loc.head(num_vars).array() - fval_loc[num_vars]).matrix(); + Eigen::VectorXd g = matprod(fdiff1, simi_loc); + A.setZero(); + if (amat != nullptr && m_lcon > 0) { + Eigen::MatrixXd amat_t = amat->transpose(); + A.leftCols(m_lcon) = amat_t; + } + int m_nlcon = num_constraints - m_lcon; + if (m_nlcon > 0) { + Eigen::VectorXd col_nlcon = conmat_loc.col(num_vars).segment(m_lcon, m_nlcon); + Eigen::MatrixXd diff = conmat_loc.block(m_lcon, 0, m_nlcon, num_vars).colwise() + - col_nlcon; + A.rightCols(m_nlcon) = matprod(diff, simi_loc).transpose(); + } + + // Calculate the trust-region trial step D. Note that D does NOT depend on CPEN. + Eigen::VectorXd d = trstlp(A, -conmat_loc.col(num_vars), delta, g); + + // Predict the change to F (PREREF) and to the constraint violation (PREREC) due to D + double preref = -inprod(d, g); // Can be negative + double ad_max = 0; + for (int i = 0; i < num_constraints; ++i) + ad_max = std::max(ad_max, conmat_loc(i, num_vars) + inprod(d, A.col(i))); + double prerec = cval_loc[num_vars] - std::max(0.0, ad_max); + + // PREREC <= 0 or PREREF >= 0 or either is NaN + if (!(prerec > 0 && preref < 0)) break; + + cpen = std::max(cpen, std::min(-2 * preref / prerec, REALMAX)); + + if (findpole(cpen, cval_loc, fval_loc) == num_vars) break; + } + + return cpen; +} + +inline std::tuple, std::vector, + std::vector, std::vector, int> +cobylb(const Calcfc& calcfc, int iprint, int maxfilt, int maxfun, + const Eigen::MatrixXd* amat, const Eigen::VectorXd* bvec, + double ctol, double cweight, double eta1, double eta2, + double ftarget, double gamma1, double gamma2, + double rhobeg, double rhoend, + const Eigen::VectorXd& constr_in, double f_in, + const Eigen::VectorXd& x_in, int maxhist, + const std::function& callback) { + + // Outputs + std::vector xhist, conhist; + std::vector fhist, chist; + + // Sizes + int m_lcon = (bvec != nullptr) ? bvec->size() : 0; + int num_constraints = constr_in.size(); + int m_nlcon = num_constraints - m_lcon; + int num_vars = x_in.size(); + + // Initialize SIM, FVAL, CONMAT, and CVAL, together with the history. + // After the initialization, SIM[:, NUM_VARS] holds the vertex of the initial + // simplex with the smallest function value (regardless of the constraint + // violation), and SIM[:, :NUM_VARS] holds the displacements from the other + // vertices to SIM[:, NUM_VARS]. FVAL, CONMAT, and CVAL hold the function values, + // constraint values, and constraint violations on the vertices in the order + // corresponding to SIM. + auto [evaluated, conmat, cval, sim, simi, fval, nf, subinfo] = + initxfc(calcfc, iprint, maxfun, constr_in, amat, bvec, ctol, f_in, + ftarget, rhobeg, x_in, xhist, fhist, chist, conhist, maxhist); + + // Initialize the filter, including xfilt, ffilt, confilt, cfilt, and nfilt. + // N.B.: The filter is used only when selecting which iterate to return. It does + // not interfere with the iterations. COBYLA is NOT a filter method but a + // trust-region method based on an L-infinity merit function. Powell's + // implementation does not use a filter to select the iterate, possibly returning + // a suboptimal iterate. + int maxfilt_actual = std::min(std::max(maxfilt, 1), maxfun); + Eigen::VectorXd cfilt = Eigen::VectorXd::Zero(maxfilt_actual); + Eigen::MatrixXd confilt = Eigen::MatrixXd::Zero(num_constraints, maxfilt_actual); + Eigen::VectorXd ffilt = Eigen::VectorXd::Zero(maxfilt_actual); + Eigen::MatrixXd xfilt = Eigen::MatrixXd::Zero(num_vars, maxfilt_actual); + + int nfilt = initfilt(conmat, ctol, cweight, cval, fval, sim, evaluated, + cfilt, confilt, ffilt, xfilt); + + // Check whether to return due to abnormal cases that may occur during init + if (subinfo != INFO_DEFAULT) { + // N.B.: selectx and findpole choose X by different standards, one cannot + // replace the other + int kopt = selectx(ffilt.head(nfilt), cfilt.head(nfilt), cweight, ctol); + Eigen::VectorXd x = xfilt.col(kopt); + double f = ffilt[kopt]; + Eigen::VectorXd constr = confilt.col(kopt); + double cstrv = cfilt[kopt]; + retmsg("COBYLA", subinfo, iprint, nf, f, x, cstrv, &constr); + return {x, f, constr, cstrv, nf, xhist, fhist, chist, conhist, subinfo}; + } + + // Set some more initial values. + // We must initialize shortd, ratio, and jdrop_tr because these get defined on + // branches that are not guaranteed to be executed, but their values are used later. + // Our initialization of CPEN differs from Powell's in two ways. First, we use the + // ratio defined in (13) of Powell's COBYLA paper to initialize CPEN. Second, we + // impose CPEN >= CPENMIN > 0. Powell's code simply initializes CPEN to 0. + double rho = rhobeg; + double delta = rhobeg; + double cpenmin = EPS; + double cpen = std::max(cpenmin, std::min(1.0e3, fcratio(conmat, fval))); + bool shortd = false; + double ratio = -1; + int jdrop_tr = 0; + Eigen::VectorXd d_last = Eigen::VectorXd::Zero(num_vars); + + // If DELTA <= GAMMA3*RHO after an update, we set DELTA to RHO. GAMMA3 must be + // less than GAMMA2. The reason: + // Imagine a very successful step with DNORM = the un-updated DELTA = RHO. TRRAD + // will update DELTA to GAMMA2*RHO. If GAMMA3 >= GAMMA2, then DELTA will be reset + // to RHO, which is not reasonable as D is very successful. See paragraph two of + // Sec 5.2.5 in T. M. Ragonneau's thesis: "Model-Based Derivative-Free Optimization + // Methods and Software." + double gamma3 = std::max(1.0, std::min(0.75 * gamma2, 1.5)); + + // MAXTR is the maximal number of trust-region iterations. Each trust-region + // iteration takes 1 or 2 function evaluations unless the trust-region step is + // short or the trust-region subproblem solver fails but the geometry step is not + // invoked. Thus the following MAXTR is unlikely to be reached. + int maxtr = 10 * maxfun; + int info = MAXTR_REACHED; + + Eigen::MatrixXd A(num_vars, num_constraints); + Eigen::VectorXd distsq(num_vars + 1); + + // Begin the iterative procedure + // After solving a trust-region subproblem, we use three boolean variables to + // control the workflow. + // SHORTD - Is the trust-region trial step too short to invoke a function eval? + // IMPROVE_GEO - Will we improve the model after the trust-region iteration? If + // yes, a geometry step will be taken, corresponding to "Branch + // (Delta)" in the COBYLA paper. + // REDUCE_RHO - Will we reduce rho after the trust-region iteration? + // COBYLA never sets IMPROVE_GEO and REDUCE_RHO to True simultaneously. + for (int tr = 0; tr < maxtr; ++tr) { + // Increase the penalty parameter CPEN, if needed, so that + // PREREM = PREREF + CPEN * PREREC > 0. + // This is the first (out of two) update of CPEN, where CPEN increases or + // remains the same. + // N.B.: CPEN and the merit function PHI = FVAL + CPEN*CVAL are used in three + // places only. + // 1. In FINDPOLE/UPDATEPOLE, deciding the optimal vertex of the simplex. + // 2. After the trust-region trial step, calculating the reduction ratio. + // 3. In GEOSTEP, deciding the direction of the geometry step. + // They do not appear explicitly in the trust-region subproblem, though the + // trust-region center (i.e. the current optimal vertex) is defined by them. + cpen = getcpen(amat, bvec, conmat, cpen, cval, delta, fval, rho, sim, simi); + + // Switch the best vertex of the current simplex to SIM[:, NUM_VARS] + std::tie(conmat, cval, fval, sim, simi, subinfo) = + updatepole(cpen, conmat, cval, fval, sim, simi); + if (subinfo == DAMAGING_ROUNDING) { info = subinfo; break; } + + // Does the interpolation set have adequate geometry? It affects improve_geo + // and reduce_rho. + bool adequate_geo = true; + for (int i = 0; i < num_vars; ++i) { + double nsq = sim.col(i).squaredNorm(); + if (nsq > 4 * delta * delta) { adequate_geo = false; break; } + } + + // Calculate the linear approximations to the objective and constraint functions. + // N.B.: TRSTLP accesses A mostly by columns, so it is more reasonable to + // store A instead of A^T. + Eigen::VectorXd fdiff2 = (fval.head(num_vars).array() - fval[num_vars]).matrix(); + Eigen::VectorXd g = matprod(fdiff2, simi); + A.setZero(); + if (amat != nullptr && m_lcon > 0) { + Eigen::MatrixXd amat_t = amat->transpose(); + A.leftCols(m_lcon) = amat_t; + } + int mnlcon = num_constraints - m_lcon; + if (mnlcon > 0) { + Eigen::VectorXd col_nlcon2 = conmat.col(num_vars).segment(m_lcon, mnlcon); + Eigen::MatrixXd diff = conmat.block(m_lcon, 0, mnlcon, num_vars).colwise() + - col_nlcon2; + A.rightCols(mnlcon) = matprod(diff, simi).transpose(); + } + + // Calculate the trust-region trial step d. Note that d does NOT depend on cpen. + d_last = trstlp(A, -conmat.col(num_vars), delta, g); + Eigen::VectorXd d = d_last; + double dnorm = std::min(delta, norm(d)); + + // Is the trust-region trial step short? N.B.: we compare DNORM with RHO, not + // DELTA. Powell's code especially defines SHORTD by SHORTD = (DNORM < 0.5 * + // RHO). In our tests 1/10 seems to work better than 1/2 or 1/4, especially + // for linearly constrained problems. + shortd = (dnorm <= 0.1 * rho); + + // Predict the change to F (PREREF) and to the constraint violation (PREREC) + // due to D. In theory: + // 1. B[:NUM_CONSTRAINTS] = -CONMAT[:, NUM_VARS] and hence + // max(B[:NUM_CONSTRAINTS] - D@A[:, :NUM_CONSTRAINTS], 0) is the L-infinity + // violation of the linearized constraints. When D=0, the violation is + // CVAL[NUM_VARS]. PREREC is the reduction of this violation achieved by D, + // which is nonnegative in theory; PREREC = 0 iff B[:NUM_CONSTRAINTS] <= 0, + // i.e. the trust-region center satisfies the linearized constraints. + // 2. PREREF may be negative or 0, but it is positive when PREREC = 0 and + // shortd is False. + // 3. Due to 2, in theory, max(PREREC, PREREF) > 0 if shortd is False. + double preref = -inprod(d, g); // Can be negative + double ad_max = 0; + for (int i = 0; i < num_constraints; ++i) + ad_max = std::max(ad_max, conmat(i, num_vars) + inprod(d, A.col(i))); + double prerec = cval[num_vars] - std::max(0.0, ad_max); + + // Evaluate PREREM, the predicted reduction in the merit function. + // In theory, PREREM >= 0 and it is 0 iff CPEN = 0 = PREREF. + double prerem = preref + cpen * prerec; + bool trfail = !(prerem > 1.0e-6 * std::min(cpen, 1.0) * rho); + + if (shortd || trfail) { + // Reduce DELTA if D is short or fails to render PREREM > 0. The latter + // can only happen due to rounding errors. + delta *= 0.1; + if (delta <= gamma3 * rho) delta = rho; + } else { + // Calculate the next value of the objective and constraint functions. + // If X is close to one of the points in the interpolation set, then we + // do not evaluate, assuming the values at the closest point. + // N.B.: If this happens, do NOT include X into the filter, as F and + // CONSTR are inaccurate. + Eigen::VectorXd x = sim.col(num_vars) + d; + distsq[num_vars] = (x - sim.col(num_vars)).squaredNorm(); + for (int i = 0; i < num_vars; ++i) + distsq[i] = (x - (sim.col(num_vars) + sim.col(i))).squaredNorm(); + + int j; + distsq.minCoeff(&j); + double f; Eigen::VectorXd constr; double cstrv; + double rhoend_sq = 1e-4 * rhoend; + rhoend_sq = rhoend_sq * rhoend_sq; + + if (distsq[j] <= rhoend_sq) { + f = fval[j]; constr = conmat.col(j); cstrv = cval[j]; + } else { + std::tie(f, constr) = evaluate(calcfc, x, m_nlcon, amat, bvec); + cstrv = 0; + for (int i = 0; i < constr.size(); ++i) cstrv = std::max(cstrv, constr[i]); + ++nf; + savehist(maxhist, x, xhist, f, fhist, cstrv, chist, constr, conhist); + nfilt = savefilt(cstrv, ctol, cweight, f, x, nfilt, cfilt, ffilt, xfilt, &constr, &confilt); + } + + fmsg("COBYLA", "Trust region", iprint, nf, delta, f, x, cstrv, &constr); + + // Evaluate ACTREM, the actual reduction in the merit function + double actrem = (fval[num_vars] + cpen * cval[num_vars]) - (f + cpen * cstrv); + + // Calculate the reduction ratio by redrat, which handles inf/nan + ratio = redrat(actrem, prerem, eta1); + + // Update DELTA. After this, DELTA < DNORM may hold. + // N.B.: + // 1. Powell's code uses RHO as the trust-region radius and updates it as: + // Reduce RHO to GAMMA1*RHO if ADEQUATE_GEO is TRUE and either SHORTD + // is TRUE or RATIO < ETA1, then revise RHO to RHOEND if its new value + // is not more than GAMMA3*RHOEND; RHO is never increased. + // 2. Our implementation uses DELTA as the trust-region radius, with RHO + // as a lower bound for DELTA. DELTA is updated in a typical trust- + // region way, and revised to RHO if its new value is not more than + // GAMMA3*RHO. RHO reflects the current resolution of the algorithm. + // 3. The same as Powell's code, we do not reduce RHO unless ADEQUATE_GEO + // is TRUE. This is also how Powell updated RHO in UOBYQA/NEWUOA/ + // BOBYQA/LINCOA. + delta = trrad(delta, dnorm, eta1, eta2, gamma1, gamma2, ratio); + if (delta <= gamma3 * rho) delta = rho; + + // Is the newly generated X better than the current best point? + bool ximproved = actrem > 0; + + // Set JDROP_TR to the index of the vertex to be replaced with X. + // JDROP_TR = -1 means there is no good point to replace, and X will not + // be included into the simplex; in this case, the geometry likely needs + // improvement, which will be handled below. + jdrop_tr = setdrop_tr(ximproved, d, delta, rho, sim, simi); + + // Update SIM, SIMI, FVAL, CONMAT, and CVAL so that SIM[:, JDROP_TR] is + // replaced with D. UPDATEXFC does nothing if JDROP_TR is -1. + std::tie(sim, simi, fval, conmat, cval, subinfo) = + updatexfc(jdrop_tr, constr, cpen, cstrv, d, f, conmat, cval, fval, sim, simi); + if (subinfo == DAMAGING_ROUNDING) { info = subinfo; break; } + + // Check whether to break due to maxfun, ftarget, etc. + subinfo = checkbreak_con(maxfun, nf, cstrv, ctol, f, ftarget, x); + if (subinfo != INFO_DEFAULT) { info = subinfo; break; } + } + // End of if SHORTD or TRFAIL. The normal trust-region calculation ends. + + // Before the next trust-region iteration, we possibly improve the geometry + // of the simplex or reduce RHO according to IMPROVE_GEO and REDUCE_RHO. + // N.B.: We must ensure that the algorithm does not set IMPROVE_GEO = True + // at infinitely many consecutive iterations without moving SIM[:, NUM_VARS] + // or reducing RHO. This is ensured by: + // 1. If an iteration sets IMPROVE_GEO to True, it must also reduce DELTA or + // set DELTA to RHO. + // 2. If SIM[:, NUM_VARS] and RHO remain unchanged, then ADEQUATE_GEO will + // become True after at most NUM_VARS invocations of GEOSTEP. + + // BAD_TRSTEP: Is the last trust-region step bad? + bool bad_trstep = shortd || trfail || ratio <= 0 || jdrop_tr == -1; + // IMPROVE_GEO: Should we take a geometry step to improve geometry? + bool improve_geo = bad_trstep && !adequate_geo; + // REDUCE_RHO: Should we enhance the resolution by reducing rho? + bool reduce_rho = bad_trstep && adequate_geo && std::max(delta, dnorm) <= rho; + + // COBYLA never sets IMPROVE_GEO and REDUCE_RHO to True simultaneously. + + // Improve the geometry of the simplex by removing a point and adding a new one. + // If the current interpolation set has acceptable geometry, skip the step. + if (improve_geo) { + bool geo_ok = true; + for (int i = 0; i < num_vars; ++i) + if (sim.col(i).squaredNorm() > 4 * delta * delta) { geo_ok = false; break; } + if (!geo_ok) { + // Decide a vertex to drop from the simplex. It will be replaced with + // SIM[:, NUM_VARS] + D to improve geometry. + // N.B.: + // 1. COBYLA never sets JDROP_GEO = num_vars. + // 2. The following JDROP_GEO comes from UOBYQA/NEWUOA/BOBYQA/LINCOA. + double max_dist = 0; + int jdrop_geo = 0; + for (int i = 0; i < num_vars; ++i) { + double nd = sim.col(i).squaredNorm(); + if (nd > max_dist) { max_dist = nd; jdrop_geo = i; } + } + + // Calculate the geometry step D + double delbar = delta / 2; + d = geostep(jdrop_geo, amat, bvec, conmat, cpen, cval, delbar, fval, simi); + Eigen::VectorXd x = sim.col(num_vars) + d; + + distsq[num_vars] = (x - sim.col(num_vars)).squaredNorm(); + for (int i = 0; i < num_vars; ++i) + distsq[i] = (x - (sim.col(num_vars) + sim.col(i))).squaredNorm(); + int j; distsq.minCoeff(&j); + double f; Eigen::VectorXd constr; double cstrv; + double rhoend_sq = 1e-4 * rhoend; + rhoend_sq = rhoend_sq * rhoend_sq; + + if (distsq[j] <= rhoend_sq) { + f = fval[j]; constr = conmat.col(j); cstrv = cval[j]; + } else { + std::tie(f, constr) = evaluate(calcfc, x, m_nlcon, amat, bvec); + cstrv = 0; + for (int i = 0; i < constr.size(); ++i) cstrv = std::max(cstrv, constr[i]); + ++nf; + savehist(maxhist, x, xhist, f, fhist, cstrv, chist, constr, conhist); + nfilt = savefilt(cstrv, ctol, cweight, f, x, nfilt, cfilt, ffilt, xfilt, &constr, &confilt); + } + + fmsg("COBYLA", "Geometry", iprint, nf, delta, f, x, cstrv, &constr); + std::tie(sim, simi, fval, conmat, cval, subinfo) = + updatexfc(jdrop_geo, constr, cpen, cstrv, d, f, conmat, cval, fval, sim, simi); + if (subinfo == DAMAGING_ROUNDING) { info = subinfo; break; } + + subinfo = checkbreak_con(maxfun, nf, cstrv, ctol, f, ftarget, x); + if (subinfo != INFO_DEFAULT) { info = subinfo; break; } + } + } + + // The calculations with the current RHO are complete. Enhance the resolution + // by reducing RHO; update DELTA and CPEN at the same time. + if (reduce_rho) { + if (rho <= rhoend) { info = SMALL_TR_RADIUS; break; } + delta = std::max(0.5 * rho, redrho(rho, rhoend)); + rho = redrho(rho, rhoend); + // The second (out of two) update of CPEN, where CPEN decreases or + // remains the same. + // Powell's code: cpen = min(cpen, fcratio(fval, conmat)), may set to 0. + cpen = std::max(cpenmin, std::min(cpen, fcratio(conmat, fval))); + Eigen::VectorXd pole_constr = conmat.col(num_vars); + rhomsg("COBYLA", iprint, nf, fval[num_vars], rho, sim.col(num_vars), + cval[num_vars], &pole_constr, cpen); + std::tie(conmat, cval, fval, sim, simi, subinfo) = + updatepole(cpen, conmat, cval, fval, sim, simi); + if (subinfo == DAMAGING_ROUNDING) { info = subinfo; break; } + } + + // Report the current best value, and check if user asks for early termination + if (callback) { + bool terminate = callback(sim.col(num_vars), fval[num_vars], nf, tr, + cval[num_vars], conmat.col(num_vars)); + if (terminate) { info = CALLBACK_TERMINATE; break; } + } + } + // End of for loop. The iterative procedure ends. + + // Return from the calculation, after trying the last trust-region step if it + // has not been tried yet. + // Ensure that D has not been updated after SHORTD == TRUE occurred, or the + // code below is incorrect. + Eigen::VectorXd x_try = sim.col(num_vars) + d_last; + if (info == SMALL_TR_RADIUS && shortd && + (x_try - sim.col(num_vars)).norm() > 1.0e-3 * rhoend && nf < maxfun) { + // UPDATEXFC or UPDATEPOLE is not called since the last trust-region step. + // Hence SIM[:, NUM_VARS] remains unchanged. Otherwise SIM[:, NUM_VARS] + D + // would not make sense. + auto [f_try, constr_try] = evaluate(calcfc, x_try, m_nlcon, amat, bvec); + double cstrv_try = 0; + for (int i = 0; i < constr_try.size(); ++i) cstrv_try = std::max(cstrv_try, constr_try[i]); + ++nf; + savehist(maxhist, x_try, xhist, f_try, fhist, cstrv_try, chist, constr_try, conhist); + nfilt = savefilt(cstrv_try, ctol, cweight, f_try, x_try, nfilt, cfilt, ffilt, xfilt, &constr_try, &confilt); + // DELTA has been updated. RHO is only indicative here. + fmsg("COBYLA", "Trust region", iprint, nf, rho, f_try, x_try, cstrv_try, &constr_try); + } + + // Return the best calculated values of the variables + // N.B.: SELECTX and FINDPOLE choose X by different standards, one cannot replace + // the other. + int kopt = selectx(ffilt.head(nfilt), cfilt.head(nfilt), std::max(cpen, cweight), ctol); + Eigen::VectorXd x_ret = xfilt.col(kopt); + double f_ret = ffilt[kopt]; + Eigen::VectorXd constr_ret = confilt.col(kopt); + double cstrv_ret = cfilt[kopt]; + + retmsg("COBYLA", info, iprint, nf, f_ret, x_ret, cstrv_ret, &constr_ret); + + return {x_ret, f_ret, constr_ret, cstrv_ret, nf, + xhist, fhist, chist, conhist, info}; +} + +} // namespace prima + +#endif // PRIMA_CPP_COBYLA_COBYLB_HPP diff --git a/cpp/src/prima/cobyla/geometry.hpp b/cpp/src/prima/cobyla/geometry.hpp new file mode 100644 index 0000000000..4d372bf66b --- /dev/null +++ b/cpp/src/prima/cobyla/geometry.hpp @@ -0,0 +1,225 @@ +#ifndef PRIMA_CPP_COBYLA_GEOMETRY_HPP +#define PRIMA_CPP_COBYLA_GEOMETRY_HPP + +// This module contains subroutines concerning the geometry-improving of the interpolation set. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include +#include + +#include "../consts.hpp" +#include "../linalg.hpp" + +namespace prima { + +inline int setdrop_tr(bool ximproved, const Eigen::VectorXd& d, double delta, + double rho, const Eigen::MatrixXd& sim, const Eigen::MatrixXd& simi) { + // This function finds (the index) of a current interpolation point to be replaced with + // the trust-region trial point. See (19)--(22) of the COBYLA paper. + // N.B.: + // 1. If XIMPROVED == True, then JDROP > 0 so that D is included into XPT. Otherwise, + // it is a bug. + // 2. COBYLA never sets JDROP = NUM_VARS + // TODO: Check whether it improves the performance if JDROP = NUM_VARS is allowed when + // XIMPROVED is True. Note that UPDATEXFC should be revised accordingly. + + int num_vars = sim.rows(); + int jdrop = -1; + + // DISTSQ[j] is the square of the distance from the jth vertex of the simplex to get "best" point so + // far, taking the trial point SIM[:, NUM_VARS] + D into account. + Eigen::VectorXd distsq(sim.cols()); + if (ximproved) { + for (int i = 0; i < num_vars; ++i) + distsq[i] = (sim.col(i) - d).squaredNorm(); + distsq[num_vars] = d.squaredNorm(); + } else { + for (int i = 0; i < num_vars; ++i) + distsq[i] = sim.col(i).squaredNorm(); + distsq[num_vars] = 0; + } + + // -------------------------------------------------------------------------------------------------- // + // The following code is Powell's scheme for defining JDROP. + // -------------------------------------------------------------------------------------------------- // + // ! JDROP = 0 by default. It cannot be removed, as JDROP may not be set below in some cases (e.g., + // ! when XIMPROVED == FALSE, MAXVAL(ABS(SIMID)) <= 1, and MAXVAL(VETA) <= EDGMAX). + // jdrop = 0 + // + // ! SIMID(J) is the value of the J-th Lagrange function at D. It is the counterpart of VLAG in UOBYQA + // ! and DEN in NEWUOA/BOBYQA/LINCOA, but it excludes the value of the (N+1)-th Lagrange function. + // simid = matprod(simi, d) + // if (any(abs(simid) > 1) .or. (ximproved .and. any(.not. is_nan(simid)))) then + // jdrop = int(maxloc(abs(simid), mask=(.not. is_nan(simid)), dim=1), kind(jdrop)) + // !!MATLAB: [~, jdrop] = max(simid, [], 'omitnan'); + // end if + // + // ! VETA(J) is the distance from the J-th vertex of the simplex to the best vertex, taking the trial + // ! point SIM(:, N+1) + D into account. + // if (ximproved) then + // veta = sqrt(sum((sim(:, 1:n) - spread(d, dim=2, ncopies=n))**2, dim=1)) + // !!MATLAB: veta = sqrt(sum((sim(:, 1:n) - d).^2)); % d should be a column! Implicit expansion + // else + // veta = sqrt(sum(sim(:, 1:n)**2, dim=1)) + // end if + // + // ! VSIG(J) (J=1, .., N) is the Euclidean distance from vertex J to the opposite face of the simplex. + // vsig = ONE / sqrt(sum(simi**2, dim=2)) + // sigbar = abs(simid) * vsig + // + // ! The following JDROP will overwrite the previous one if its premise holds. + // mask = (veta > factor_delta * delta .and. (sigbar >= factor_alpha * delta .or. sigbar >= vsig)) + // if (any(mask)) then + // jdrop = int(maxloc(veta, mask=mask, dim=1), kind(jdrop)) + // !!MATLAB: etamax = max(veta(mask)); jdrop = find(mask & ~(veta < etamax), 1, 'first'); + // end if + // + // ! Powell's code does not include the following instructions. With Powell's code, if SIMID consists + // ! of only NaN, then JDROP can be 0 even when XIMPROVED == TRUE (i.e., D reduces the merit function). + // ! With the following code, JDROP cannot be 0 when XIMPROVED == TRUE, unless VETA is all NaN, which + // ! should not happen if X0 does not contain NaN, the trust-region/geometry steps never contain NaN, + // ! and we exit once encountering an iterate containing Inf (due to overflow). + // if (ximproved .and. jdrop <= 0) then ! Write JDROP <= 0 instead of JDROP == 0 for robustness. + // jdrop = int(maxloc(veta, mask=(.not. is_nan(veta)), dim=1), kind(jdrop)) + // !!MATLAB: [~, jdrop] = max(veta, [], 'omitnan'); + // end if + // -------------------------------------------------------------------------------------------------- // + // Powell's scheme ends here. + // -------------------------------------------------------------------------------------------------- // + + // The following definition of JDROP is inspired by SETDROP_TR in UOBYQA/NEWUOA/BOBYQA/LINCOA. + // It is simpler and works better than Powell's scheme. Note that we allow JDROP to be NUM_VARS+1 if + // XIMPROVED is True, whereas Powell's code does not. + // See also (4.1) of Scheinberg-Toint-2010: Self-Correcting Geometry in Model-Based Algorithms for + // Derivative-Free Unconstrained Optimization, which refers to the strategy here as the "combined + // distance/poisedness criteria". + + double denom = std::max(rho, delta / 10.0); + double denom_sq = denom * denom; + Eigen::VectorXd weight(distsq.size()); + for (int i = 0; i < weight.size(); ++i) + weight[i] = std::max(1.0, distsq[i] / denom_sq); + // Similar to Powell's NEWUOA code. + + // Other possible definitions of weight. They work almost the same as the one above. + // weight = distsq // Similar to Powell's LINCOA code, but WRONG. See comments in LINCOA/geometry.f90. + // weight = max(1, max(25 * distsq / delta**2)) // Similar to Powell's BOBYQA code, works well. + // weight = max(1, max(10 * distsq / delta**2)) + // weight = max(1, max(1e2 * distsq / delta**2)) + // weight = max(1, max(distsq / rho**2)) ! Similar to Powell's UOBYQA + + // If 0 <= j < NUM_VARS, SIMID[j] is the value of the jth Lagrange function at D; the value of the + // (NUM_VARS+1)th Lagrange function is 1 - sum(SIMID). [SIMID, 1 - sum(SIMID)] is the counterpart of + // VLAG in UOBYQA and DEN in NEWUOA/BOBYQA/LINCOA. + Eigen::VectorXd simid = matprod(simi, d); + double simid_sum = simid.sum(); + Eigen::VectorXd score(sim.cols()); + for (int i = 0; i < num_vars; ++i) + score[i] = weight[i] * std::abs(simid[i]); + score[num_vars] = weight[num_vars] * std::abs(1.0 - simid_sum); + + // If XIMPROVED = False (D does not render a better X), set SCORE[NUM_VARS] = -1 to avoid JDROP = NUM_VARS. + if (!ximproved) + score[num_vars] = -1; + + // score[j] is NaN implies SIMID[j] is NaN, but we want abs(SIMID) to be big. So we + // exclude such j. + for (int i = 0; i < score.size(); ++i) { + if (std::isnan(score[i])) + score[i] = -1; + } + + // The following if statement works a bit better than + // `if any(score > 1) or (any(score > 0) and ximproved)` from Powell's UOBYQA and + // NEWUOA code. + bool any_positive = false; + for (int i = 0; i < score.size(); ++i) { + if (score[i] > 0) { any_positive = true; break; } + } + + if (any_positive) { // Powell's BOBYQA and LINCOA code. + double max_score = -1; + for (int i = 0; i < score.size(); ++i) { + if (score[i] > max_score) { + max_score = score[i]; + jdrop = i; + } + } + } + + // JDROP >= 1 when XIMPROVED = TRUE unless NaN occurs in DISTSQ, which should not happen if the + // starting point does not contain NaN and the trust-region/geometry steps never contain NaN. + if (ximproved && jdrop < 0) { + double max_distsq = -1; + for (int i = 0; i < distsq.size(); ++i) { + if (distsq[i] > max_distsq) { + max_distsq = distsq[i]; + jdrop = i; + } + } + } + + return jdrop; +} + +inline Eigen::VectorXd geostep(int jdrop, const Eigen::MatrixXd* amat, + const Eigen::VectorXd* bvec, + const Eigen::MatrixXd& conmat, double cpen, + const Eigen::VectorXd& cval, double delbar, + const Eigen::VectorXd& fval, const Eigen::MatrixXd& simi) { + // This function calculates a geometry step so that the geometry of the interpolation set is improved + // when SIM[: JDROP_GEO] is replaced with SIM[:, NUM_VARS] + D. See (15)--(17) of the COBYLA paper. + + int m_lcon = (bvec != nullptr) ? bvec->size() : 0; + int num_constraints = conmat.rows(); + int num_vars = simi.rows(); + + // SIMI[JDROP, :] is a vector perpendicular to the face of the simplex to the opposite of vertex + // JDROP. Set D to the vector in this direction and with length DELBAR. + Eigen::VectorXd d = simi.row(jdrop).transpose(); + d = delbar * (d / norm(d)); + + // The code below chooses the direction of D according to an approximation of the merit function. + // See (17) of the COBYLA paper and line 225 of Powell's cobylb.f. + + // Calculate the coefficients of the linear approximations to the objective and constraint functions. + // N.B.: CONMAT and SIMI have been updated after the last trust-region step, but G and A have not. + // So we cannot pass G and A from outside. + Eigen::VectorXd fdiff = (fval.head(num_vars).array() - fval[num_vars]).matrix(); + Eigen::VectorXd g = matprod(fdiff, simi); + + Eigen::MatrixXd A(num_vars, num_constraints); + A.setZero(); + if (amat != nullptr && m_lcon > 0) { + A.leftCols(m_lcon) = amat->transpose(); + } + int m_nlcon = num_constraints - m_lcon; + if (m_nlcon > 0) { + Eigen::MatrixXd diff = conmat.block(m_lcon, 0, m_nlcon, num_vars).colwise() + - conmat.col(num_vars).segment(m_lcon, m_nlcon); + A.rightCols(m_nlcon) = matprod(diff, simi).transpose(); + } + + // CVPD and CVND are the predicted constraint violation of D and -D by the linear models. + double cvpd = 0; + double cvnd = 0; + for (int i = 0; i < num_constraints; ++i) { + double adi = inprod(d, A.col(i)); + cvpd = std::max(cvpd, conmat(i, num_vars) + adi); + cvnd = std::max(cvnd, conmat(i, num_vars) - adi); + } + + if (-inprod(d, g) + cpen * cvnd < inprod(d, g) + cpen * cvpd) + d *= -1; + + return d; +} + +} // namespace prima + +#endif // PRIMA_CPP_COBYLA_GEOMETRY_HPP diff --git a/cpp/src/prima/cobyla/initialize.hpp b/cpp/src/prima/cobyla/initialize.hpp new file mode 100644 index 0000000000..86a3706ab3 --- /dev/null +++ b/cpp/src/prima/cobyla/initialize.hpp @@ -0,0 +1,149 @@ +#ifndef PRIMA_CPP_COBYLA_INITIALIZE_HPP +#define PRIMA_CPP_COBYLA_INITIALIZE_HPP + +#include +#include +#include +#include +#include + +#include "../checkbreak.hpp" +#include "../consts.hpp" +#include "../evaluate.hpp" +#include "../infos.hpp" +#include "../history.hpp" +#include "../linalg.hpp" +#include "../message.hpp" +#include "../selectx.hpp" + +namespace prima { + +inline std::tuple +initxfc(const Calcfc& calcfc, int iprint, int maxfun, + const Eigen::VectorXd& constr0, const Eigen::MatrixXd* amat, + const Eigen::VectorXd* bvec, double ctol, double f0, double ftarget, + double rhobeg, const Eigen::VectorXd& x0, + std::vector& xhist, std::vector& fhist, + std::vector& chist, std::vector& conhist, int maxhist) { + // This subroutine does the initialization concerning X, function values, and + // constraints. + + // Local variables + int num_constraints = constr0.size(); + int m_lcon = bvec ? bvec->size() : 0; + int m_nlcon = num_constraints - m_lcon; + int num_vars = x0.size(); + + // Initialize info to the default value. At return, a value different from this + // value will indicate an abnormal return. + int info = INFO_DEFAULT; + + // Initialize the simplex. It will be revised during the initialization. + Eigen::MatrixXd sim = Eigen::MatrixXd::Zero(num_vars, num_vars + 1); + sim.rightCols(1).col(0) = x0; + for (int i = 0; i < num_vars; ++i) sim(i, i) = rhobeg; + + // Initialize the matrix simi. In most cases simi is overwritten, but not always. + Eigen::MatrixXd simi = Eigen::MatrixXd::Identity(num_vars, num_vars) / rhobeg; + + // evaluated[j] = true iff the function/constraint of sim.col(j) has been evaluated. + Eigen::VectorXi evaluated = Eigen::VectorXi::Zero(num_vars + 1); + + // Initialize fval + Eigen::VectorXd fval = Eigen::VectorXd::Constant(num_vars + 1, REALMAX); + Eigen::VectorXd cval = Eigen::VectorXd::Constant(num_vars + 1, REALMAX); + Eigen::MatrixXd conmat = Eigen::MatrixXd::Constant(num_constraints, num_vars + 1, REALMAX); + + for (int k = 0; k < num_vars + 1; ++k) { + // We will evaluate F corresponding to sim.col(j). + Eigen::VectorXd x = sim.col(num_vars); + int j; + double f; + Eigen::VectorXd constr; + if (k == 0) { + j = num_vars; + f = f0; + constr = constr0; + } else { + j = k - 1; + x[j] += rhobeg; + auto [fv, cv] = evaluate(calcfc, x, m_nlcon, amat, bvec); + f = fv; constr = cv; + } + double cstrv = 0; + for (int i = 0; i < constr.size(); ++i) cstrv = std::max(cstrv, constr[i]); + + // Print a message about the function/constraint evaluation according to IPRINT. + fmsg("COBYLA", "Initialization", iprint, k, rhobeg, f, x, cstrv, &constr); + + // Save X, F, CONSTR, CSTRV into the history. + savehist(maxhist, x, xhist, f, fhist, cstrv, chist, constr, conhist); + + // Save F, CONSTR, and CSTRV to FVAL, CONMAT, and CVAL respectively. + evaluated[j] = 1; + fval[j] = f; + conmat.col(j) = constr; + cval[j] = cstrv; + + // Check whether to exit. + int subinfo = checkbreak_con(maxfun, k, cstrv, ctol, f, ftarget, x); + if (subinfo != INFO_DEFAULT) { info = subinfo; break; } + + // Exchange the new vertex of the initial simplex with the optimal vertex if necessary. + // This is the ONLY part that is essentially non-parallel. + if (j < num_vars && fval[j] < fval[num_vars]) { + std::swap(fval[j], fval[num_vars]); + std::swap(cval[j], cval[num_vars]); + conmat.col(j).swap(conmat.col(num_vars)); + sim.col(num_vars) = x; + // SIM(:, :j+1) is lower triangular + for (int ii = 0; ii <= j; ++ii) sim(j, ii) = -rhobeg; + } + } + + int nf = evaluated.sum(); + + if (evaluated.minCoeff() == 1) { + // Initialize SIMI to the inverse of SIM[:, :num_vars] + simi = inv(sim.leftCols(num_vars)); + } + + return {evaluated, conmat, cval, sim, simi, fval, nf, info}; +} + +inline int initfilt(const Eigen::MatrixXd& conmat, double ctol, double cweight, + const Eigen::VectorXd& cval, const Eigen::VectorXd& fval, + const Eigen::MatrixXd& sim, const Eigen::VectorXi& evaluated, + Eigen::VectorXd& cfilt, Eigen::MatrixXd& confilt, + Eigen::VectorXd& ffilt, Eigen::MatrixXd& xfilt) { + // This function initializes the filter (XFILT, etc) that will be used when selecting + // X at the end of the solver. + // N.B.: + // 1. Why not initialize the filters using XHIST, etc? Because the history is empty if + // the user chooses not to output it. + // 2. We decouple INITXFC and INITFILT so that it is easier to parallelize the former + // if needed. + + int num_constraints = conmat.rows(); + int num_vars = sim.rows(); + int maxfilt = ffilt.size(); + int nfilt = 0; + + for (int i = 0; i < num_vars + 1; ++i) { + if (evaluated[i]) { + Eigen::VectorXd x; + if (i < num_vars) x = sim.col(i) + sim.col(num_vars); + else x = sim.col(i); // i == num_vars, i.e. the last column + Eigen::VectorXd ci = conmat.col(i); + nfilt = savefilt(cval[i], ctol, cweight, fval[i], x, nfilt, cfilt, ffilt, xfilt, + &ci, &confilt); + } + } + + return nfilt; +} + +} // namespace prima + +#endif // PRIMA_CPP_COBYLA_INITIALIZE_HPP diff --git a/cpp/src/prima/cobyla/trustregion.hpp b/cpp/src/prima/cobyla/trustregion.hpp new file mode 100644 index 0000000000..c597fe0df0 --- /dev/null +++ b/cpp/src/prima/cobyla/trustregion.hpp @@ -0,0 +1,546 @@ +#ifndef PRIMA_CPP_COBYLA_TRUSTREGION_HPP +#define PRIMA_CPP_COBYLA_TRUSTREGION_HPP + +// This module provides subroutines concerning the trust-region calculations of COBYLA. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include +#include + +#include "../consts.hpp" +#include "../linalg.hpp" +#include "../powalg.hpp" + +namespace prima { + +namespace detail { + +// TRSTLP_SUB: the real calculations for trstlp, both stage 1 and stage 2. +// Major differences between stage 1 and stage 2: +// 1. Initialization. Stage 2 inherits the values of some variables from stage 1, so +// they are initialized in stage 1 but not in stage 2. +// 2. cviol. Updated after each iteration in stage 1, constant in stage 2. +// 3. sdirn. See the definition in the code for details. +// 4. optnew. The two stages have different objectives. +// 5. step <= cviol in stage 1. +// +// Powell's code can encounter infinite cycling, which did happen when testing CUTEst +// problems: DANWOODLS, GAUSS1LS, GAUSS2LS, GAUSS3LS, KOEBHELB, TAX13322, TAXR13322. +// Indeed, in all these cases, Inf/NaN appear in d due to extremely large values in +// A (up to 10^219). To resolve this, we set the maximal number of iterations to +// maxiter, and terminate if Inf/NaN occurs in d. +inline std::tuple +trstlp_sub(Eigen::VectorXi iact, int nact, int stage, + const Eigen::MatrixXd& A, const Eigen::VectorXd& b, double delta, + const Eigen::VectorXd& d_in, const Eigen::VectorXd& vmultc_in, + const Eigen::MatrixXd& z_in) { + int mcon = A.cols(); + int num_vars = A.rows(); + + Eigen::VectorXd zdasav = Eigen::VectorXd::Zero(z_in.cols()); + Eigen::VectorXd vmultd = Eigen::VectorXd::Zero(vmultc_in.size()); + Eigen::VectorXd zdota = Eigen::VectorXd::Zero(z_in.cols()); + + Eigen::MatrixXd z = z_in; + Eigen::VectorXd vmultc = vmultc_in; + Eigen::VectorXd d = d_in; + + int num_constraints = 0; + int icon = 0; + double cviol = 0; + + // Initialize according to stage + if (stage == 1) { + for (int i = 0; i < mcon; ++i) iact[i] = i; + nact = 0; + d.setZero(); + if (mcon == 0) { + // Quick return: no constraints in stage 1 + return {iact, nact, d, vmultc, z}; + } + cviol = std::max(0.0, -b.minCoeff()); + vmultc = Eigen::VectorXd::Constant(mcon, cviol) + b; + z = Eigen::MatrixXd::Identity(num_vars, num_vars); + + // Check whether a quick return is possible + if (cviol <= 0) { + return {iact, nact, d, vmultc, z}; + } + + bool all_b_nan = true; + for (int i = 0; i < b.size(); ++i) + if (!std::isnan(b[i])) { all_b_nan = false; break; } + if (all_b_nan) { + return {iact, nact, d, vmultc, z}; + } + + // icon = index of most violated constraint (max(-b)) + icon = 0; + double max_neg_b = -std::numeric_limits::infinity(); + for (int i = 0; i < mcon; ++i) { + if (!std::isnan(b[i]) && -b[i] > max_neg_b) { + max_neg_b = -b[i]; + icon = i; + } + } + num_constraints = mcon; + } else { + // Stage 2: check whether a quick return is possible + if (d.squaredNorm() >= delta * delta) { + return {iact, nact, d, vmultc, z}; + } + + iact[mcon - 1] = mcon - 1; + vmultc[mcon - 1] = 0; + num_constraints = mcon - 1; + icon = mcon - 1; + + // Recalculate cviol so it need not be passed from stage 1 to 2 + cviol = 0; + for (int i = 0; i < num_constraints; ++i) { + double adi = inprod(d, A.col(i)); + cviol = std::max(cviol, adi - b[i]); + } + cviol = std::max(0.0, cviol); + } + + // zdota: scalar products of Z columns with active constraint gradients + for (int k = 0; k < nact; ++k) + zdota[k] = inprod(z.col(k), A.col(iact[k])); + + // More initialization + double optold = REALMAX; + int nactold = nact; + int nfail = 0; + int maxiter = std::min(10000, 100 * std::max(num_constraints, num_vars)); + + // vmultd is computed from scratch at each iteration, but vmultc is inherited + Eigen::VectorXd sdirn(num_vars); + sdirn.setZero(); + + for (int iter = 0; iter < maxiter; ++iter) { + double optnew; + if (stage == 1) + optnew = cviol; + else + optnew = inprod(d, A.col(mcon - 1)); + + // End the current stage if 3 consecutive iterations have failed to reduce + // the best calculated objective value or to increase the number of active + // constraints since the best value was calculated. This prevents cycling, + // but there is a remote possibility of premature termination. + if (optnew < optold || nact > nactold) { + nactold = nact; + nfail = 0; + } else { + nfail += 1; + } + optold = std::min(optold, optnew); + if (nfail == 3) break; + + // If icon exceeds nact, add the constraint with index iact[icon] to the + // active set. + if (icon >= nact) { + zdasav.head(nact) = zdota.head(nact); + int nactsav = nact; + + Eigen::VectorXd new_col = A.col(iact[icon]); + auto [z_new, zdota_new, nact_new] = qradd_Rdiag(new_col, z, zdota, nact); + z = z_new; + zdota = zdota_new; + nact = nact_new; + + if (nact == nactsav + 1) { + // qradd_Rdiag succeeded: the constraint was added to the active set + if (nact != icon + 1) { + std::swap(vmultc[icon], vmultc[nact - 1]); + std::swap(iact[icon], iact[nact - 1]); + } else { + vmultc[nact - 1] = 0; + } + } else { + // qradd_Rdiag failed to add the constraint. + // If the active set is empty, there is nothing to do. + if (nact == 0) break; + + // VMULTD is calculated from scratch for the first (out of 2) time + // in this iteration. + // Note that IACT has not been updated to replace IACT[NACT] with + // IACT[ICON]. Thus A[:, IACT[:NACT]] is the UNUPDATED version + // before QRADD (note Z[:, :NACT] remains the same before and after + // QRADD). Therefore if we supply ZDOTA to LSQR (as Rdiag) as Powell + // did, we should use the UNUPDATED version, namely ZDASAV. + Eigen::MatrixXd A_sub(num_vars, nact); + for (int k = 0; k < nact; ++k) A_sub.col(k) = A.col(iact[k]); + Eigen::VectorXd v = lsqr(A_sub, A.col(iact[icon]), + z.leftCols(nact), zdasav.head(nact)); + vmultd.head(nact) = v; + + // N.B.: This can be triggered by NACT == 0 (among other + // possibilities)! This is important, because NACT will be used as + // an index in the sequel. However, the nact == 0 case is handled + // above. + bool any_positive = false; + for (int k = 0; k < nact; ++k) { + if (vmultd[k] > 0 && iact[k] <= num_constraints) { + any_positive = true; + break; + } + } + if (!any_positive) + break; + + // vmultd[NACT+1:mcon] is not used, but initialize to avoid issues + vmultd.tail(mcon - nact).setConstant(-1); + + // Revise the Lagrange multipliers. + // Only the places with vmultd > 0 and iact <= num_constraints matter + Eigen::VectorXd fracmult(nact); + for (int k = 0; k < nact; ++k) { + if (vmultd[k] > 0 && iact[k] <= num_constraints) + fracmult[k] = vmultc[k] / vmultd[k]; + else + fracmult[k] = REALMAX; + } + double frac = fracmult.minCoeff(); + + for (int k = 0; k < nact; ++k) + vmultc[k] = std::max(0.0, vmultc[k] - frac * vmultd[k]); + + // Exit if the new value of zdota[nact] is not acceptable. + // Powell's condition: not abs(zdota[nact]) > 0. + // Note that it is different from 'abs(zdota[nact]) <= 0)' as + // zdota[nact] can be NaN. + // N.B.: We cannot arrive here with nact == 0, which should have + // triggered a break above. + if (std::isnan(zdota[nact - 1]) || std::abs(zdota[nact - 1]) <= EPS * EPS) + break; + + // Reorder the active constraints so that the one to be replaced is + // at the end of the list. + vmultc[icon] = 0; + vmultc[nact - 1] = frac; + std::swap(iact[icon], iact[nact - 1]); + } + + // In stage 2, ensure that the objective continues to be treated as the + // last active constraint. + if (stage == 2 && iact[nact - 1] != (mcon - 1)) { + if (nact <= 1) { + // We must exit, as nact-2 is used as an index below. + break; + } + Eigen::MatrixXd A_sub_full(num_vars, nact); + for (int k = 0; k < nact; ++k) A_sub_full.col(k) = A.col(iact[k]); + auto [z_exc, zdota_exc] = qrexc_Rdiag(A_sub_full, z, zdota.head(nact), nact - 2); + z = z_exc; + zdota.head(nact) = zdota_exc; + std::swap(iact[nact - 2], iact[nact - 1]); + std::swap(vmultc[nact - 2], vmultc[nact - 1]); + } + + // Powell's code does not have the following. It avoids subsequent + // floating point exceptions. + if (std::isnan(zdota[nact - 1]) || std::abs(zdota[nact - 1]) <= EPS * EPS) + break; + + // Set sdirn to the direction of the next change to the current vector. + // During stage 1, sdirn gives a search direction that reduces all the + // active constraint violations by one simultaneously. + if (stage == 1) { + sdirn -= ((inprod(sdirn, A.col(iact[nact - 1])) + 1) / zdota[nact - 1]) * z.col(nact - 1); + } else { + sdirn = (-1.0 / zdota[nact - 1]) * z.col(nact - 1); + } + } else { + // icon < nact: delete the constraint with index iact[icon] from the + // active set. Reorder iact[icon:nact] into [iact[icon+1:nact], iact[icon]] + // then reduce nact to nact-1. + Eigen::MatrixXd A_sub_full(num_vars, nact); + for (int k = 0; k < nact; ++k) A_sub_full.col(k) = A.col(iact[k]); + auto [z_exc, zdota_exc] = qrexc_Rdiag(A_sub_full, z, zdota.head(nact), icon); + z = z_exc; + zdota.head(nact) = zdota_exc; + + int saved_icon_val = iact[icon]; + int saved_vmultc_val = vmultc[icon]; + for (int k = icon; k < nact - 1; ++k) { + iact[k] = iact[k + 1]; + vmultc[k] = vmultc[k + 1]; + } + iact[nact - 1] = saved_icon_val; + vmultc[nact - 1] = saved_vmultc_val; + nact -= 1; + + // In theory, nact > 0 in stage 2, as the objective function should + // always be considered as an "active constraint". However, looking at + // the code, it's not guaranteed. It did happen in stage 1 that nact + // became 0 after the reduction after almost one year of random tests. + if (stage == 2 && nact < 0) break; + // Powell's code does not have the following. Avoids exceptions. + if (nact > 0) { + if (std::isnan(zdota[nact - 1]) || std::abs(zdota[nact - 1]) <= EPS * EPS) + break; + } + + // Set sdirn to the direction of the next change to the current vector. + if (stage == 1) { + if (nact < num_vars) + sdirn -= inprod(sdirn, z.col(nact)) * z.col(nact); + // sdirn is orthogonal to z[:, nact+1] + } else { + sdirn = (-1.0 / zdota[nact - 1]) * z.col(nact - 1); + } + } + + // Calculate the step to the trust-region boundary or the step that reduces + // cviol to 0. + // The following calculation of step is adopted from NEWUOA/BOBYQA/LINCOA. + // It seems to improve the performance of COBYLA. We also found that removing + // the precaution about underflows is beneficial --- the underflows are + // harmless anyway. + double dd = delta * delta - d.squaredNorm(); + double ss = sdirn.squaredNorm(); + double sd = inprod(sdirn, d); + if (dd <= 0 || ss <= EPS * delta * delta || std::isnan(sd)) break; + + // sqrtd: square root of a discriminant. The max avoids sqrtd < abs(sd) due + // to underflow. + double sqrtd = std::max(std::sqrt(ss * dd + sd * sd), std::abs(sd)); + sqrtd = std::max(sqrtd, std::sqrt(ss * dd)); + + double step; + if (sd > 0) + step = dd / (sqrtd + sd); + else + step = (sqrtd - sd) / ss; + + // step < 0 should not happen. Step can be 0 or NaN when, e.g., sd or ss + // becomes inf. + if (step <= 0 || !std::isfinite(step)) break; + + if (stage == 1) { + if (isminor(cviol, step)) break; + step = std::min(step, cviol); + } + + // Set dnew to the new variables if step is the steplength, and reduce cviol + // to the corresponding maximum residual if stage 1 is being done. + Eigen::VectorXd dnew = d + step * sdirn; + + if (stage == 1) { + double max_res = 0; + for (int k = 0; k < nact; ++k) { + double res = inprod(dnew, A.col(iact[k])) - b[iact[k]]; + max_res = std::max(max_res, res); + } + cviol = std::max(0.0, max_res); + } + + // vmultd is computed from scratch for the second (out of 2) time in one + // iteration. vmultd[:nact] and vmultd[nact:mcon] are calculated separately + // with no coupling. vmultd will be calculated from scratch again in the + // next iteration. + // Set vmultd to the vmultc vector that would occur if d became dnew. + // A device is included to force vmultd[k] = 0 if deviations can be + // attributed to computer rounding errors. + Eigen::VectorXd vmultd_new(mcon); + { + Eigen::MatrixXd A_sub_nact(num_vars, nact); + for (int k = 0; k < nact; ++k) A_sub_nact.col(k) = A.col(iact[k]); + Eigen::VectorXd v = lsqr(A_sub_nact, dnew, z.leftCols(nact), zdota.head(nact)); + for (int k = 0; k < nact; ++k) vmultd_new[k] = -v[k]; + } + if (stage == 2) + vmultd_new[nact - 1] = std::max(0.0, vmultd_new[nact - 1]); + + // Complete vmultd by finding the new constraint residuals. + Eigen::VectorXd cvshift(mcon); + for (int k = 0; k < mcon; ++k) + cvshift[k] = cviol - (inprod(dnew, A.col(iact[k])) - b[iact[k]]); + Eigen::VectorXd cvsabs(mcon); + for (int k = 0; k < mcon; ++k) + cvsabs[k] = inprod(dnew.cwiseAbs(), A.col(iact[k]).cwiseAbs()) + std::abs(b[iact[k]]) + cviol; + apply_isminor(cvshift, cvsabs); + for (int k = nact; k < mcon; ++k) + vmultd_new[k] = cvshift[k]; + + // Calculate the fraction of the step from d to dnew that will be taken. + double frac = 1.0; + icon = -1; + for (int k = 0; k < mcon; ++k) { + if (vmultd_new[k] < 0) { + double fk = vmultc[k] / (vmultc[k] - vmultd_new[k]); + if (fk < frac) { + frac = fk; + icon = k; + } + } + } + + // Update d, vmultc, and cviol + Eigen::VectorXd dold = d; + d = (1.0 - frac) * d + frac * dnew; + for (int k = 0; k < mcon; ++k) + vmultc[k] = std::max(0.0, (1.0 - frac) * vmultc[k] + frac * vmultd_new[k]); + + // Break in the case of inf/nan in d or vmultc + if (!std::isfinite(d.cwiseAbs().sum()) || !std::isfinite(vmultc.cwiseAbs().sum())) { + d = dold; + break; + } + + if (stage == 1) { + double max_res2 = 0; + for (int k = 0; k < mcon; ++k) { + double res = inprod(d, A.col(k)) - b[k]; + max_res2 = std::max(max_res2, res); + } + // In theory, cviol = max(d@A - b, 0), yet the cviol updated as above + // can be quite different from this value if A has huge entries (> 1e20). + cviol = std::max(0.0, max_res2); + } + + // In Powell's code, the condition is icon == 0. Indeed, icon < 0 cannot + // hold unless fracmult contains only nan, which should not happen; icon >= + // mcon should never occur. + if (icon < 0 || icon >= mcon) break; + } + + return {iact, nact, d, vmultc, z}; +} + +} // namespace detail + +// TRSTLP: calculate an n-component vector d by the following two stages. +// Stage 1: d is set to the shortest vector that minimizes the greatest violation +// of the constraints A^T D <= B, subject to the Euclidean length of d being at +// most delta. If its length is strictly less than delta, then stage 2 uses the +// resultant freedom in d to minimize G^T D subject to no increase in any greatest +// constraint violation. +// +// cviol is the largest constraint violation of the current d: max(max(A^T D - b), 0). +// icon is the index of a most violated constraint if cviol is positive. +// +// nact is the number of constraints in the active set and iact[0..nact-1] are +// their indices, while the remainder of iact contains a permutation of the +// remaining constraint indices. +// N.B.: nact <= min(num_constraints, num_vars). +// +// Z is an orthogonal matrix whose first nact columns can be regarded as the +// result of Gram-Schmidt applied to the active constraint gradients. For +// j = 0..nact-1, zdota[j] is the scalar product of the jth column of Z with +// the gradient of the jth active constraint. d is the current vector of variables +// and here the residuals of the active constraints should be zero. +// +// N.B.: +// 0. In Powell's implementation, the constraints are A^T D >= B. In other words, +// the A and B in our implementation are the negative of those in Powell's. +// 1. The algorithm was NOT documented in the COBYLA paper. +// 2. As a major part of the algorithm (see trstlp_sub), the code maintains and +// updates the QR factorization of A[iact[:nact]], i.e. the gradients of all +// active constraints. Z is indeed Q, and zdota is the diagonal of R. +// 3. There are probably better algorithms available for this problem. +inline Eigen::VectorXd trstlp(const Eigen::MatrixXd& A, const Eigen::VectorXd& b, + double delta, const Eigen::VectorXd& g) { + int num_constraints = A.cols(); + int num_vars = A.rows(); + + Eigen::VectorXd vmultc; + Eigen::VectorXi iact; + int nact = 0; + Eigen::VectorXd d = Eigen::VectorXd::Zero(num_vars); + Eigen::MatrixXd z = Eigen::MatrixXd::Identity(num_vars, num_vars); + + if (num_constraints <= 0) { + // No constraints: steepest descent direction, scaled to trust region + double gn = norm(g); + if (gn > 0) d = -(delta / gn) * g; + return d; + } + + vmultc = Eigen::VectorXd::Zero(num_constraints + 1); + iact = Eigen::VectorXi::Zero(num_constraints + 1); + + // Form A_aug and b_aug. This allows the gradient of the objective to be + // regarded as the gradient of a constraint in the second stage. + Eigen::MatrixXd A_aug(num_vars, num_constraints + 1); + A_aug.leftCols(num_constraints) = A; + A_aug.rightCols(1).col(0) = g; + Eigen::VectorXd b_aug(num_constraints + 1); + b_aug.head(num_constraints) = b; + b_aug[num_constraints] = 0; + + // Scale the problem if A contains large values. Otherwise floating point + // exceptions may occur. Note that the trust-region step is scale invariant. + for (int i = 0; i < num_constraints + 1; ++i) { + double maxval = A_aug.col(i).cwiseAbs().maxCoeff(); + if (maxval > 1e12) { + double modscal = std::max(2.0 * REALMIN, 1.0 / maxval); + A_aug.col(i) *= modscal; + b_aug[i] *= modscal; + } + } + + // Stage 1: minimize the L-infinity constraint violation of the linearized constraints + { + Eigen::VectorXi iact_sub = iact.head(num_constraints); + Eigen::VectorXd vmultc_sub = vmultc.head(num_constraints); + auto [iact_out, nact_out, d_out, vmultc_out, z_out] = + detail::trstlp_sub(iact_sub, nact, 1, + A_aug.leftCols(num_constraints), + b_aug.head(num_constraints), + delta, d, vmultc_sub, z); + iact.head(num_constraints) = iact_out; + nact = nact_out; + d = d_out; + vmultc.head(num_constraints) = vmultc_out; + z = z_out; + } + + // Stage 2: minimize the linearized objective without increasing the L-infinity + // constraint violation + { + auto [iact_out, nact_out, d_out, vmultc_out, z_out] = + detail::trstlp_sub(iact, nact, 2, A_aug, b_aug, delta, d, vmultc, z); + iact = iact_out; + nact = nact_out; + d = d_out; + vmultc = vmultc_out; + z = z_out; + } + + return d; +} + +// TRRAD: update the trust-region radius according to RATIO and DNORM. +inline double trrad(double delta_in, double dnorm, double eta1, double eta2, + double gamma1, double gamma2, double ratio) { + double delta; + if (ratio <= eta1) { + // Powell's UOBYQA/NEWUOA + // Powell's COBYLA/LINCOA: delta = gamma1 * delta_in + // Powell's BOBYQA: delta = min(gamma1 * delta_in, dnorm) + delta = gamma1 * dnorm; + } else if (ratio <= eta2) { + // Powell's UOBYQA/NEWUOA/BOBYQA/LINCOA + delta = std::max(gamma1 * delta_in, dnorm); + } else { + // Powell's NEWUOA/BOBYQA + // Modified version: delta = max(delta_in, gamma2 * dnorm). Works well for + // UOBYQA. For noise-free CUTEst problems of <= 100 variables, Powell's + // version works slightly better. + // Powell's LINCOA: delta = min(max(gamma1*delta_in, gamma2*dnorm), gamma3*delta_in) + delta = std::max(gamma1 * delta_in, gamma2 * dnorm); + } + return delta; +} + +} // namespace prima + +#endif // PRIMA_CPP_COBYLA_TRUSTREGION_HPP diff --git a/cpp/src/prima/cobyla/update.hpp b/cpp/src/prima/cobyla/update.hpp new file mode 100644 index 0000000000..8fb08bb4cb --- /dev/null +++ b/cpp/src/prima/cobyla/update.hpp @@ -0,0 +1,236 @@ +//----------------------------------------------------------------------------------------------------! +// This module contains subroutines concerning the update of the interpolation set. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +// +// Python translation by Nickolai Belakovski. +// C++ translation by ... +//----------------------------------------------------------------------------------------------------! + +#ifndef PRIMA_CPP_COBYLA_UPDATE_HPP +#define PRIMA_CPP_COBYLA_UPDATE_HPP + +#include +#include +#include +#include + +#include "../consts.hpp" +#include "../infos.hpp" +#include "../linalg.hpp" + +namespace prima { + +inline int findpole(double cpen, const Eigen::VectorXd& cval, const Eigen::VectorXd& fval) { + // This subroutine identifies the best vertex of the current simplex with respect to the merit + // function PHI = F + CPEN * CSTRV. + int num_vars = fval.size() - 1; + int jopt = num_vars; + Eigen::VectorXd phi = fval + cpen * cval; + double phimin = phi.minCoeff(); + + // Identify the optimal vertex of the current simplex + // Essentially jopt = argmin(phi). However, we keep jopt = num_vars unless there + // is a strictly better choice. When there are multiple choices, we choose the jopt + // with the smallest value of cval. + if (phimin < phi[jopt] || ((cval.array() < cval[jopt]).any() && (phi.array() <= phi[jopt]).any())) { + // While we could use argmin(phi), there may be two places where phi achieves + // phimin, and in that case we should choose the one with the smallest cval. + double best_cv = std::numeric_limits::infinity(); + for (int i = 0; i <= num_vars; ++i) { + if (phi[i] <= phimin && cval[i] < best_cv) { + best_cv = cval[i]; + jopt = i; + } + } + } + return jopt; +} + +inline std::tuple +updatepole(double cpen, const Eigen::MatrixXd& conmat, const Eigen::VectorXd& cval, + const Eigen::VectorXd& fval, const Eigen::MatrixXd& sim, + const Eigen::MatrixXd& simi) { + //----------------------------------------------------------------------------------------------! + // This subroutine identifies the best vertex of the current simplex with respect to the merit + // function PHI = F + CPEN * CSTRV, and then switch this vertex to SIM[:, NUM_VARS], which Powell + // called the "pole position" in his comments. CONMAT, CVAL, FVAL, and SIMI are updated + // accordingly. + // + // N.B. 1: In precise arithmetic, the following two procedures produce the same results: + // 1) apply UPDATEPOLE to SIM twice, first with CPEN = CPEN1 and then with CPEN = CPEN2; + // 2) apply UPDATEPOLE to SIM with CPEN = CPEN2. + // In finite-precision arithmetic, however, they may produce different results unless + // CPEN1 = CPEN2. + // + // N.B. 2: When JOPT == N+1, the best vertex is already at the pole position, so there is nothing + // to switch. However, as in Powell's code, the code below will check whether SIMI is good enough + // to work as the inverse of SIM(:, 1:N) or not. If not, Powell's code would invoke an error + // return of COBYLB; our implementation, however, will try calculating SIMI from scratch; if the + // recalculated SIMI is still of poor quality, then UPDATEPOLE will return with INFO = + // DAMAGING_ROUNDING, informing COBYLB that SIMI is poor due to damaging rounding errors. + // + // N.B. 3: UPDATEPOLE should be called when and only when FINDPOLE can potentially returns a value + // other than N+1. The value of FINDPOLE is determined by CPEN, CVAL, and FVAL, the latter two + // being decided by SIM. Thus UPDATEPOLE should be called after CPEN or SIM changes. COBYLA + // updates CPEN at only two places: the beginning of each trust-region iteration, and when REDRHO + // is called; SIM is updated only by UPDATEXFC, which itself calls UPDATEPOLE internally. + // Therefore, we only need to call UPDATEPOLE after updating CPEN at the beginning of each + // trust-region iteration and after each invocation of REDRHO. + //----------------------------------------------------------------------------------------------! + + int num_constraints = conmat.rows(); + int num_vars = sim.rows(); + // INFO must be set, as it is an output. + int info = INFO_DEFAULT; + + // Identify the optimal vertex of the current simplex. + int jopt = findpole(cpen, cval, fval); + + Eigen::MatrixXd sim_new = sim; + Eigen::MatrixXd simi_new = simi; + Eigen::VectorXd fval_new = fval; + Eigen::MatrixXd conmat_new = conmat; + Eigen::VectorXd cval_new = cval; + + // Switch the best vertex to the pole position SIM[:, NUM_VARS] if it is not there already and + // update SIMI. Before the update, save a copy of SIM and SIMI. If the update is unsuccessful due + // to damaging rounding errors, we restore them and return with INFO = DAMAGING_ROUNDING. + Eigen::MatrixXd sim_old = sim; + Eigen::MatrixXd simi_old = simi; + + if (jopt >= 0 && jopt < num_vars) { + // Unless there is a bug in FINDPOLE it is guaranteed that JOPT >= 0 + // When JOPT == NUM_VARS, there is nothing to switch; in addition SIMI[JOPT, :] will be + // illegal. + sim_new.col(num_vars) += sim_new.col(jopt); + Eigen::VectorXd sim_jopt = sim_new.col(jopt); + sim_new.col(jopt).setZero(); + for (int i = 0; i < num_vars; ++i) + sim_new.col(i) -= sim_jopt; + + // The above update is equivalent to multiplying SIM[:, :NUM_VARS] from the right side by a + // matrix whose JOPT-th row is [-1, -1, ..., -1], while all the other rows are the same as + // those of the identity matrix. It is easy to check that the inverse of this matrix is + // itself. Therefore, SIMI should be updated by a multiplication with this matrix (i.e. its + // inverse) from the left side, as is done in the following line. The JOPT-th row of the + // updated SIMI is minus the sum of all rows of the original SIMI, whereas all the other rows + // remain unchanged. + simi_new.row(jopt) = -simi.colwise().sum(); + } + + // Check whether SIMI is a poor approximation to the inverse of SIM[:, :NUM_VARS] + // Calculate SIMI from scratch if the current one is damaged by rounding errors. + Eigen::MatrixXd sim_left = sim_new.leftCols(num_vars); + double erri = (matprod(simi_new, sim_left) - + Eigen::MatrixXd::Identity(num_vars, num_vars)).cwiseAbs().maxCoeff(); + double itol = 1; + if (erri > 0.1 * itol || std::isnan(erri)) { + Eigen::MatrixXd simi_test = inv(sim_left); + double erri_test = (matprod(simi_test, sim_left) - + Eigen::MatrixXd::Identity(num_vars, num_vars)).cwiseAbs().maxCoeff(); + if (erri_test < erri || (std::isnan(erri) && !std::isnan(erri_test))) { + simi_new = simi_test; + erri = erri_test; + } + } + + // If SIMI is satisfactory, then update FVAL, CONMAT, and CVAL. Otherwise restore SIM and SIMI, + // and return with INFO = DAMAGING_ROUNDING. + if (erri <= itol) { + if (jopt >= 0 && jopt < num_vars) { + std::swap(fval_new[jopt], fval_new[num_vars]); + conmat_new.col(jopt).swap(conmat_new.col(num_vars)); + std::swap(cval_new[jopt], cval_new[num_vars]); + } + } else { + info = DAMAGING_ROUNDING; + sim_new = sim_old; + simi_new = simi_old; + } + + return {conmat_new, cval_new, fval_new, sim_new, simi_new, info}; +} + +inline std::tuple +updatexfc(int jdrop, const Eigen::VectorXd& constr, double cpen, double cstrv, + const Eigen::VectorXd& d, double f, const Eigen::MatrixXd& conmat, + const Eigen::VectorXd& cval, const Eigen::VectorXd& fval, + const Eigen::MatrixXd& sim, const Eigen::MatrixXd& simi) { + // This function revises the simplex by updating the elements of SIM, SIMI, FVAL, CONMAT, and + // CVAL. + + int num_constraints = conmat.rows(); + int num_vars = sim.rows(); + int info = INFO_DEFAULT; + + // Do nothing when JDROP < 0. This can only happen after a trust-region step. + // JDROP < 0 is impossible if the input is correct. + if (jdrop < 0) { + return {sim, simi, fval, conmat, cval, info}; + } + + Eigen::MatrixXd sim_new = sim; + Eigen::MatrixXd simi_new = simi; + + if (jdrop < num_vars) { + sim_new.col(jdrop) = d; + double inprod_val = inprod(simi_new.row(jdrop), d); + Eigen::VectorXd simi_jdrop = simi_new.row(jdrop).transpose() / inprod_val; + simi_new -= outprod(matprod(simi_new, d), simi_jdrop); + simi_new.row(jdrop) = simi_jdrop.transpose(); + } else { + sim_new.col(num_vars) += d; + for (int i = 0; i < num_vars; ++i) + sim_new.col(i) -= d; + + Eigen::VectorXd simid = matprod(simi_new, d); + double sum_simid = simid.sum(); + Eigen::VectorXd sum_simi = primasum(simi_new, 0); + simi_new += outprod(simid, sum_simi / (1.0 - sum_simid)); + } + + // Check whether SIMI is a poor approximation to the inverse of SIM[:, :NUM_VARS] + // Calculate SIMI from scratch if the current one is damaged by rounding errors. + double itol = 1; + Eigen::MatrixXd left = sim_new.leftCols(num_vars); + double erri = (matprod(simi_new, left) - + Eigen::MatrixXd::Identity(num_vars, num_vars)).cwiseAbs().maxCoeff(); + if (erri > 0.1 * itol || std::isnan(erri)) { + Eigen::MatrixXd simi_test = inv(left); + double erri_test = (matprod(simi_test, left) - + Eigen::MatrixXd::Identity(num_vars, num_vars)).cwiseAbs().maxCoeff(); + if (erri_test < erri || (std::isnan(erri) && !std::isnan(erri_test))) { + simi_new = simi_test; + erri = erri_test; + } + } + + // If SIMI is satisfactory, then update FVAL, CONMAT, CVAL, and the pole position. Otherwise + // restore SIM and SIMI, and return with INFO = DAMAGING_ROUNDING. + if (erri <= itol) { + Eigen::VectorXd fval_new = fval; + fval_new[jdrop] = f; + Eigen::MatrixXd conmat_new = conmat; + conmat_new.col(jdrop) = constr; + Eigen::VectorXd cval_new = cval; + cval_new[jdrop] = cstrv; + + // Switch the best vertex to the pole position SIM[:, NUM_VARS] if it is not there already + std::tie(conmat_new, cval_new, fval_new, sim_new, simi_new, info) = + updatepole(cpen, conmat_new, cval_new, fval_new, sim_new, simi_new); + + return {sim_new, simi_new, fval_new, conmat_new, cval_new, info}; + } else { + info = DAMAGING_ROUNDING; + return {sim, simi, fval, conmat, cval, info}; + } +} + +} // namespace prima + +#endif // PRIMA_CPP_COBYLA_UPDATE_HPP diff --git a/cpp/src/prima/consts.hpp b/cpp/src/prima/consts.hpp new file mode 100644 index 0000000000..f29345573f --- /dev/null +++ b/cpp/src/prima/consts.hpp @@ -0,0 +1,57 @@ +#ifndef PRIMA_CPP_CONSTS_HPP +#define PRIMA_CPP_CONSTS_HPP + +/* +This is a module defining some constants. + +Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. + +Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +*/ + +#include +#include + +namespace prima { + +inline constexpr bool DEBUGGING = false; + +// Smallest positive normalized double (denorm_min gives smallest subnormal) +inline constexpr double REALMIN = std::numeric_limits::denorm_min(); +inline constexpr double REALMAX = std::numeric_limits::max(); +inline constexpr double FUNCMAX = 1.0e30; +inline constexpr double CONSTRMAX = FUNCMAX; +inline double EPS = std::numeric_limits::epsilon(); + +// Any bound with an absolute value at least BOUNDMAX is considered as no bound. +inline constexpr double BOUNDMAX = REALMAX / 4.0; + +// Some default values +inline constexpr double RHOBEG_DEFAULT = 1.0; +inline constexpr double RHOEND_DEFAULT = 1e-6; +inline constexpr double FTARGET_DEFAULT = -REALMAX; +inline const double CTOL_DEFAULT = std::sqrt(std::numeric_limits::epsilon()); +inline constexpr double CWEIGHT_DEFAULT = 1e8; +inline constexpr double ETA1_DEFAULT = 0.1; +inline constexpr double ETA2_DEFAULT = 0.7; +inline constexpr double GAMMA1_DEFAULT = 0.5; +inline constexpr double GAMMA2_DEFAULT = 2.0; +inline constexpr int IPRINT_DEFAULT = 0; +inline constexpr int MAXFUN_DIM_DEFAULT = 500; + +// 1MB > 10^5 * REAL64. 100 can be too small. +inline constexpr int PRIMA_MAX_HIST_MEM_MB = 300; + +// Maximal amount of memory (Byte) allowed for XHIST, FHIST, CONHIST, CHIST, and the filters. +inline constexpr double MHM = PRIMA_MAX_HIST_MEM_MB * 1.0e6; + +// Make sure that MAXHISTMEM does not exceed HUGE(0) to avoid overflow and memory errors. +inline constexpr int MAXHISTMEM = static_cast(std::min(MHM, static_cast(std::numeric_limits::max()))); + +// Maximal length of the filter used in constrained solvers. Should be positive; < 200 is not recommended. +inline constexpr int MIN_MAXFILT = 200; +inline constexpr int MAXFILT_DEFAULT = 10 * MIN_MAXFILT; + +} // namespace prima + +#endif // PRIMA_CPP_CONSTS_HPP diff --git a/cpp/src/prima/evaluate.hpp b/cpp/src/prima/evaluate.hpp new file mode 100644 index 0000000000..7706eac308 --- /dev/null +++ b/cpp/src/prima/evaluate.hpp @@ -0,0 +1,106 @@ +#ifndef PRIMA_CPP_EVALUATE_HPP +#define PRIMA_CPP_EVALUATE_HPP + +#include +#include +#include +#include +#include + +#include "consts.hpp" +#include "linalg.hpp" + +// This is a module evaluating the objective/constraint function with Nan/Inf handling. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +namespace prima { + +using Calcfc = std::function(const Eigen::VectorXd&)>; + +// This function moderates a decision variable. It replaces NaN by 0 and Inf/-Inf by +// REALMAX/-REALMAX. +inline Eigen::VectorXd moderatex(const Eigen::VectorXd& x) { + Eigen::VectorXd y = x; + for (Eigen::Index i = 0; i < y.size(); ++i) { + if (std::isnan(y[i])) y[i] = 0; + y[i] = std::max(-REALMAX, std::min(REALMAX, y[i])); + } + return y; +} + +// This function moderates the function value of a MINIMIZATION problem. It replaces +// NaN and any value above FUNCMAX by FUNCMAX. +// We may moderate huge negative function values as follows, but we decide not to. +// f = std::max(-FUNCMAX, std::min(FUNCMAX, f)); +inline double moderatef(double f) { + if (std::isnan(f)) f = FUNCMAX; + return std::max(-REALMAX, std::min(FUNCMAX, f)); +} + +// This function moderates the constraint value, the constraint demanding this value +// to be NONNEGATIVE. It replaces any value below -CONSTRMAX by -CONSTRMAX, and any +// NaN or value above CONSTRMAX by CONSTRMAX. +inline double moderatec(double c) { + if (std::isnan(c) || c > CONSTRMAX) c = CONSTRMAX; + return std::max(-CONSTRMAX, std::min(CONSTRMAX, c)); +} + +// This function moderates the constraint value, the constraint demanding this value +// to be NONNEGATIVE. It replaces any value below -CONSTRMAX by -CONSTRMAX, and any +// NaN or value above CONSTRMAX by CONSTRMAX. +inline Eigen::VectorXd moderatec(const Eigen::VectorXd& c) { + Eigen::VectorXd y = c; + for (Eigen::Index i = 0; i < y.size(); ++i) { + if (std::isnan(y[i]) || y[i] > CONSTRMAX) y[i] = CONSTRMAX; + y[i] = std::max(-CONSTRMAX, std::min(CONSTRMAX, y[i])); + } + return y; +} + +// This function evaluates CALCFC at X, returning the objective function value and the +// constraint value. Nan/Inf are handled by a moderated extreme barrier. +inline std::pair +evaluate(const Calcfc& calcfc, const Eigen::VectorXd& x, int m_nlcon, + const Eigen::MatrixXd* amat, const Eigen::VectorXd* bvec) { + + // Sizes + int m_lcon = (bvec != nullptr) ? bvec->size() : 0; + int total_con = m_lcon + m_nlcon; + Eigen::VectorXd constr(total_con); + + if (amat != nullptr) { + constr.head(m_lcon) = matprod(x, amat->transpose()) - *bvec; + } + + double f; + bool has_nan = false; + for (Eigen::Index i = 0; i < x.size(); ++i) { + if (std::isnan(x[i])) { has_nan = true; break; } + } + + if (has_nan) { + // Although this should not happen unless there is a bug, we include this case + // for robustness. + f = primasum(x); + constr.tail(m_nlcon).setConstant(f); + } else { + auto [fc, cc] = calcfc(moderatex(x)); + f = fc; + if (m_nlcon > 0) constr.tail(m_nlcon) = cc; + + // Moderated extreme barrier: replace NaN/huge objective or constraint values + // with a large but finite value. This is naive, and better approaches surely + // exist. + f = moderatef(f); + if (m_nlcon > 0) constr.tail(m_nlcon) = moderatec(constr.tail(m_nlcon)); + } + + return {f, constr}; +} + +} // namespace prima + +#endif // PRIMA_CPP_EVALUATE_HPP diff --git a/cpp/src/prima/history.hpp b/cpp/src/prima/history.hpp new file mode 100644 index 0000000000..98c15518de --- /dev/null +++ b/cpp/src/prima/history.hpp @@ -0,0 +1,42 @@ +#ifndef PRIMA_CPP_HISTORY_HPP +#define PRIMA_CPP_HISTORY_HPP + +// This module provides subroutines that handle the X/F/C histories of the solver, taking into +// account that MAXHIST may be smaller than NF. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include + +namespace prima { + +// Save the data values to the history lists. +inline void savehist(int maxhist, const Eigen::VectorXd& x, + std::vector& xhist, + double f, std::vector& fhist, + double cstrv, std::vector& chist, + const Eigen::VectorXd& constr, std::vector& conhist) { + if (static_cast(xhist.size()) < maxhist) { + xhist.push_back(x); + fhist.push_back(f); + chist.push_back(cstrv); + conhist.push_back(constr); + } else { + // This effectively accomplishes what rangehist does in the Fortran implementation. + xhist.erase(xhist.begin()); + fhist.erase(fhist.begin()); + chist.erase(chist.begin()); + conhist.erase(conhist.begin()); + xhist.push_back(x); + fhist.push_back(f); + chist.push_back(cstrv); + conhist.push_back(constr); + } +} + +} // namespace prima + +#endif // PRIMA_CPP_HISTORY_HPP diff --git a/cpp/src/prima/infos.hpp b/cpp/src/prima/infos.hpp new file mode 100644 index 0000000000..471f114ad2 --- /dev/null +++ b/cpp/src/prima/infos.hpp @@ -0,0 +1,36 @@ +#ifndef PRIMA_CPP_INFOS_HPP +#define PRIMA_CPP_INFOS_HPP + +// This is a module defining exit flags. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +namespace prima { + +constexpr int INFO_DEFAULT = 0; +constexpr int SMALL_TR_RADIUS = 0; +constexpr int FTARGET_ACHIEVED = 1; +constexpr int TRSUBP_FAILED = 2; +constexpr int MAXFUN_REACHED = 3; +constexpr int FIXED_SUCCESS = 13; +constexpr int MAXTR_REACHED = 20; +constexpr int NAN_INF_X = -1; +constexpr int NAN_INF_F = -2; +constexpr int NAN_INF_MODEL = -3; +constexpr int NO_SPACE_BETWEEN_BOUNDS = 6; +constexpr int DAMAGING_ROUNDING = 7; +constexpr int ZERO_LINEAR_CONSTRAINT = 8; +constexpr int CALLBACK_TERMINATE = 30; + +// Stop-codes. +// The following codes are used by ERROR STOP as stop-codes, which should be default integers. +constexpr int INVALID_INPUT = 100; +constexpr int ASSERTION_FAILS = 101; +constexpr int VALIDATION_FAILS = 102; +constexpr int MEMORY_ALLOCATION_FAILS = 103; + +} // namespace prima + +#endif // PRIMA_CPP_INFOS_HPP diff --git a/cpp/src/prima/linalg.hpp b/cpp/src/prima/linalg.hpp new file mode 100644 index 0000000000..e02c7d40e6 --- /dev/null +++ b/cpp/src/prima/linalg.hpp @@ -0,0 +1,404 @@ +#ifndef PRIMA_CPP_LINALG_HPP +#define PRIMA_CPP_LINALG_HPP + +/* +This module provides some basic linear algebra procedures. + +Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. + +Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include "consts.hpp" + +namespace prima { + +/* +We use naive implementations of matrix multiplication and other routines for two +reasons: +1. When Fortran is compiled in debug mode, and C++ is using these routines, we + can get bit for bit identical results as compared to Fortran. This is helpful + for comparing the two implementations. It will be particularly helpful when porting + the other implementations like LINCOA, etc. +2. On some problems this algorithm is very sensitive to errors in finite precision + arithmetic. Switching to naive implementation will slow down the algorithm, but + may be more stable. +*/ +inline bool USE_NAIVE_MATH = false; + +/* +This function tests whether x is minor compared to ref. It is used by Powell, e.g., in COBYLA. +In precise arithmetic, isminor(x, ref) is true if and only if x == 0; in floating point +arithmetic, isminor(x, ref) is true if x is 0 or its nonzero value can be attributed to +computer rounding errors according to ref. +Larger sensitivity means the function is more strict/precise, the value 0.1 being due to Powell. + +For example: +isminor(1e-20, 1e300) -> true, because in floating point arithmetic 1e-20 cannot be added to +1e300 without being rounded to 1e300. +isminor(1e300, 1e-20) -> false, because in floating point arithmetic adding 1e300 to 1e-20 +dominates the latter number. +isminor(3, 4) -> false, because 3 can be added to 4 without being rounded off +*/ +inline bool isminor(double x, double ref) { + double sensitivity = 0.1; + double refa = std::abs(ref) + sensitivity * std::abs(x); + double refb = std::abs(ref) + 2 * sensitivity * std::abs(x); + return std::abs(ref) >= refa || refa >= refb; +} + +// Element-wise isminor for vectors +inline Eigen::VectorXi isminor(const Eigen::VectorXd& x, const Eigen::VectorXd& ref) { + Eigen::VectorXi r(x.size()); + for (Eigen::Index i = 0; i < x.size(); ++i) r[i] = isminor(x[i], ref[i]) ? 1 : 0; + return r; +} + +// Dot product; uses Eigen's dot or naive implementation depending on USE_NAIVE_MATH +inline double inprod(const Eigen::VectorXd& x, const Eigen::VectorXd& y) { + if (!USE_NAIVE_MATH) return x.dot(y); + double r = 0; + for (Eigen::Index i = 0; i < x.size(); ++i) r += x[i] * y[i]; + return r; +} + +// Matrix-vector product for 1x2 shapes (row vector * matrix) +inline Eigen::VectorXd matprod12(const Eigen::VectorXd& x, const Eigen::MatrixXd& y) { + Eigen::VectorXd r(y.cols()); + for (Eigen::Index i = 0; i < y.cols(); ++i) r[i] = inprod(x, y.col(i)); + return r; +} + +// Matrix-vector product for 2x1 shapes (matrix * column vector) +inline Eigen::VectorXd matprod21(const Eigen::MatrixXd& x, const Eigen::VectorXd& y) { + Eigen::VectorXd r(x.rows()); + r.setZero(); + for (Eigen::Index i = 0; i < x.cols(); ++i) r += x.col(i) * y[i]; + return r; +} + +// Matrix-matrix product (naive triple-loop) +inline Eigen::MatrixXd matprod22(const Eigen::MatrixXd& x, const Eigen::MatrixXd& y) { + Eigen::MatrixXd r(x.rows(), y.cols()); + r.setZero(); + for (Eigen::Index i = 0; i < y.cols(); ++i) + for (Eigen::Index j = 0; j < x.cols(); ++j) + r.col(j) += x.col(i) * y(i, j); + return r; +} + +// Dispatches matprod12 or uses Eigen's multiply based on USE_NAIVE_MATH +inline Eigen::VectorXd matprod(const Eigen::VectorXd& x, const Eigen::MatrixXd& y) { + if (!USE_NAIVE_MATH) return x.transpose() * y; + return matprod12(x, y); +} + +// Dispatches matprod21 or uses Eigen's multiply based on USE_NAIVE_MATH +inline Eigen::VectorXd matprod(const Eigen::MatrixXd& x, const Eigen::VectorXd& y) { + if (!USE_NAIVE_MATH) return x * y; + return matprod21(x, y); +} + +// Dispatches matprod22 or uses Eigen's multiply based on USE_NAIVE_MATH +inline Eigen::MatrixXd matprod(const Eigen::MatrixXd& x, const Eigen::MatrixXd& y) { + if (!USE_NAIVE_MATH) return x * y; + return matprod22(x, y); +} + +// Outer product; uses Eigen or naive implementation +inline Eigen::MatrixXd outprod(const Eigen::VectorXd& x, const Eigen::VectorXd& y) { + if (!USE_NAIVE_MATH) return x * y.transpose(); + Eigen::MatrixXd r(x.size(), y.size()); + for (Eigen::Index i = 0; i < y.size(); ++i) r.col(i) = x * y[i]; + return r; +} + +/* +Solve the linear least squares problem using Powell's method. +If USE_NAIVE_MATH is false, uses Eigen's SVD-based solver. +Otherwise, uses the modified Gram-Schmidt / QR approach as in Powell's implementation, +where columns are processed in reverse order with isminor checks to detect +rank-deficient systems. +*/ +inline Eigen::VectorXd lsqr(const Eigen::MatrixXd& A, const Eigen::VectorXd& b, + const Eigen::MatrixXd& Q, const Eigen::VectorXd& Rdiag) { + int n = A.cols(); + if (n == 0) return Eigen::VectorXd(); + + // Use the Powell QR-based solver with precomputed Q and Rdiag (the diagonal of R). + // This is the same algorithm as the Fortran lsqr_Rdiag, which is the original + // COBYLA approach. See PRIMA's fortran/common/linalg.f90. + // In precise arithmetic, the algorithm performs back-substitution on the + // upper-triangular system R * x = Q^T * b. It forces x[j] = 0 when the + // computed value can be attributed to computer rounding errors (isminor). + int rank = std::min(static_cast(A.rows()), n); + Eigen::VectorXd x = Eigen::VectorXd::Zero(n); + Eigen::VectorXd y = b; + for (int i = rank - 1; i >= 0; --i) { + double yq = inprod(y, Q.col(i)); + double yqa = inprod(y.cwiseAbs(), Q.col(i).cwiseAbs()); + if (isminor(yq, yqa)) { x[i] = 0; } + else { x[i] = yq / Rdiag[i]; y -= x[i] * A.col(i); } + } + return x; +} + +// Compute sqrt(x1^2 + x2^2) with overflow/underflow protection +inline double hypot(double x1, double x2) { + if (!USE_NAIVE_MATH) return std::hypot(x1, x2); + if (!std::isfinite(x1)) return std::abs(x1); + if (!std::isfinite(x2)) return std::abs(x2); + double a = std::min(std::abs(x1), std::abs(x2)); + double b = std::max(std::abs(x1), std::abs(x2)); + if (a > std::sqrt(REALMIN) && b < std::sqrt(REALMAX / 2.1)) { + return std::sqrt(a * a + b * b); + } else if (b > 0) { + return b * std::sqrt((a / b) * (a / b) + 1); + } + return 0; +} + +// Euclidean norm; uses Eigen or naive implementation +inline double norm(const Eigen::VectorXd& x) { + if (!USE_NAIVE_MATH) return x.norm(); + // NOTE: Avoid std::pow! And exponentiation in general! + double r = 0; + for (Eigen::Index i = 0; i < x.size(); ++i) r += x[i] * x[i]; + return std::sqrt(r); +} + +/* +According to its documentation, Eigen's sum may sometimes do partial pairwise summation. +For our purposes, when comparing, we want don't want to do anything fancy, and we +just want to add things up one at a time. +*/ +inline double primasum(const Eigen::VectorXd& x) { return x.sum(); } + +inline double primasum(const Eigen::MatrixXd& x) { return x.sum(); } + +// Sum along specified axis +inline Eigen::VectorXd primasum(const Eigen::MatrixXd& x, int axis) { + if (axis == 0) return x.colwise().sum(); + return x.rowwise().sum(); +} + +// Element-wise square +inline Eigen::MatrixXd primapow2(const Eigen::MatrixXd& x) { return x.cwiseProduct(x); } + +/* +Believe it or not, x**2 is not always the same as x*x in some languages. +In Fortran they appear to be identical, but in C++ we always use x*x. +*/ +inline Eigen::VectorXd primapow2(const Eigen::VectorXd& x) { return x.cwiseProduct(x); } + +/* +As in MATLAB, planerot(x) returns a 2x2 Givens matrix G for x in R2 so that Y = G @ x has Y[1] = 0. +Roughly speaking, G = [[x[0]/R, x[1]/R], [-x[1]/R, x[0]/R]], where R = norm(x). +0. We need to take care of the possibilities of R=0, Inf, NaN, and over/underflow. +1. The G defined above is continuous with respect to X except at 0. Following this definition, + G = [[sign(x[0]), 0], [0, sign(x[0])]] if x[1] == 0, + G = [[0, sign(x[1])], [sign(x[1]), 0]] if x[0] == 0 + Yet some implementations ignore the signs, leading to discontinuity and numerical instability. +2. Difference from MATLAB: if x contains NaN or consists of only Inf, MATLAB returns a NaN matrix, + but we return an identity matrix or a matrix of +/-sqrt(2). We intend to keep G always orthogonal. +*/ +inline Eigen::MatrixXd planerot(const Eigen::Vector2d& x) { + // Define C = X(1) / R and S = X(2) / R with R = HYPOT(X(1), X(2)). Handle Inf/NaN, over/underflow. + double c, s; + if (std::isnan(x[0]) || std::isnan(x[1])) { + // In this case, MATLAB sets G to NaN(2, 2). We refrain from doing so to keep G orthogonal. + c = 1; s = 0; + } + else if (!std::isfinite(x[0]) && !std::isfinite(x[1])) { + // In this case, MATLAB sets G to NaN(2, 2). We refrain from doing so to keep G orthogonal. + c = 1.0 / std::sqrt(2) * (x[0] >= 0 ? 1 : -1); + s = 1.0 / std::sqrt(2) * (x[1] >= 0 ? 1 : -1); + } + // X(1) == 0 == X(2) + else if (std::abs(x[0]) <= 0 && std::abs(x[1]) <= 0) { c = 1; s = 0; } + /* + N.B.: + 0. With <= instead of <, this case covers X(1) == 0 == X(2), which is treated above separately + to avoid the confusing SIGN(., 0). + 1. SIGN(A, 0) = ABS(A) in Fortran but sign(0) = 0 in MATLAB, Python, Julia, and R. + 2. Taking SIGN(X(1)) into account ensures the continuity of G with respect to X except at 0. + */ + else if (std::abs(x[1]) <= EPS * std::abs(x[0])) { c = x[0] >= 0 ? 1 : -1; s = 0; } + /* + N.B.: SIGN(A, X) = ABS(A) * sign of X != A * sign of X. Therefore, it is WRONG to define G + as SIGN(RESHAPE([ZERO, -ONE, ONE, ZERO], [2, 2]), X(2)). This mistake was committed on + 20211206 and took a whole day to debug! NEVER use SIGN on arrays unless you are really sure. + */ + else if (std::abs(x[0]) <= EPS * std::abs(x[1])) { c = 0; s = x[1] >= 0 ? 1 : -1; } + else { + /* + Here is the normal case. It implements the Givens rotation in a stable & continuous way as in: + Bindel, D., Demmel, J., Kahan, W., and Marques, O. (2002). On computing Givens rotations + reliably and efficiently. ACM Transactions on Mathematical Software (TOMS), 28(2), 206-238. + N.B.: 1. Modern compilers compute SQRT(REALMIN) and SQRT(REALMAX/2.1) at compilation time. + 2. The direct calculation without involving T and U seems to work better; use it if possible. + */ + if ((std::sqrt(REALMIN) < std::abs(x[0]) && std::abs(x[0]) < std::sqrt(REALMAX / 2.1)) && + (std::sqrt(REALMIN) < std::abs(x[1]) && std::abs(x[1]) < std::sqrt(REALMAX / 2.1))) { + // Do NOT use HYPOTENUSE here; the best implementation for one may be suboptimal for the other + double r = norm(x); + c = x[0] / r; + s = x[1] / r; + } else if (std::abs(x[0]) > std::abs(x[1])) { + double t = x[1] / x[0]; + double u = std::max(1.0, std::max(std::abs(t), std::sqrt(1 + t * t))); // MAXVAL: precaution against rounding error. + u *= (x[0] >= 0 ? 1 : -1); + c = 1 / u; + s = t / u; + } else { + double t = x[0] / x[1]; + double u = std::max(1.0, std::max(std::abs(t), std::sqrt(1 + t * t))); // MAXVAL: precaution against rounding error. + u *= (x[1] >= 0 ? 1 : -1); + c = t / u; + s = 1 / u; + } + } + Eigen::Matrix2d G; + G << c, s, -s, c; // G = [c, s; -s, c] + return G; +} + +// Set entries of cvshift to zero where they are minor compared to cvsabs +inline void apply_isminor(Eigen::VectorXd& cvshift, const Eigen::VectorXd& cvsabs) { + for (Eigen::Index i = 0; i < cvshift.size(); ++i) { + if (isminor(cvshift[i], cvsabs[i])) cvshift[i] = 0; + } +} + +// Test if A is lower triangular within tolerance +inline bool istril(const Eigen::MatrixXd& A, double tol = 0) { + return (A.cwiseAbs() - A.triangularView().toDenseMatrix().cwiseAbs()).lpNorm() <= tol; +} + +// Test if A is upper triangular within tolerance +inline bool istriu(const Eigen::MatrixXd& A, double tol = 0) { + return (A.cwiseAbs() - A.triangularView().toDenseMatrix().cwiseAbs()).lpNorm() <= tol; +} + +// Compute A^T * B +inline Eigen::MatrixXd matprod_transpose(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B) { + return A.transpose() * B; +} + +/* +Compute the inverse of A. Uses Eigen's inverse by default. +If USE_NAIVE_MATH is true, exploits triangular structure (used in COBYLA) +or falls back to QR-based inversion. +*/ +inline Eigen::MatrixXd inv(const Eigen::MatrixXd& A) { + if (!USE_NAIVE_MATH) return A.inverse(); + int n = A.rows(); + if (istril(A)) { + // This case is invoked in COBYLA. + Eigen::MatrixXd R = A.transpose(); + Eigen::MatrixXd B = Eigen::MatrixXd::Zero(n, n); + for (int i = 0; i < n; ++i) { + B(i, i) = 1.0 / R(i, i); + if (i > 0) { + Eigen::MatrixXd Bsub = B.block(0, 0, i, i); + B.block(0, i, i, 1) = -matprod(Bsub, Eigen::MatrixXd(R.block(0, i, i, 1))) / R(i, i); + } + } + return B.transpose(); + } else if (istriu(A)) { + Eigen::MatrixXd B = Eigen::MatrixXd::Zero(n, n); + for (int i = 0; i < n; ++i) { + B(i, i) = 1.0 / A(i, i); + if (i > 0) { + Eigen::MatrixXd Bsub = B.block(0, 0, i, i); + B.block(0, i, i, 1) = -matprod(Bsub, Eigen::MatrixXd(A.block(0, i, i, 1))) / A(i, i); + } + } + return B; + } else { + // This is NOT the best algorithm for the inverse, but since the QR subroutine is available ... + Eigen::HouseholderQR qr(A); + Eigen::MatrixXd R = qr.matrixQR().triangularView().toDenseMatrix(); + Eigen::MatrixXd Q = qr.householderQ(); + Eigen::MatrixXd B = Eigen::MatrixXd::Zero(n, n); + for (int i = n - 1; i >= 0; --i) { + Eigen::MatrixXd Bsub = B.rightCols(n - 1 - i); + Eigen::MatrixXd Rsub = R.block(i + 1, i, n - 1 - i, 1); + B.col(i) = (Q.col(i) - matprod(Bsub, Rsub)) / R(i, i); + } + return B; + } +} + +// QR decomposition with column pivoting +inline std::tuple qr(const Eigen::MatrixXd& A) { + Eigen::ColPivHouseholderQR qr(A); + Eigen::MatrixXd Q = qr.householderQ(); + Eigen::MatrixXd R = qr.matrixQR().triangularView().toDenseMatrix(); + Eigen::VectorXi P = qr.colsPermutation().indices(); + return {Q, R, P}; +} + +// Test whether A = B^{-1} up to tolerance tol +inline bool isinv(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B, double tol = 0) { + int n = A.rows(); + double t = (tol > 0) ? tol : std::min(1e-3, 100.0 * EPS * std::max(n, 1)); + t = std::max({t, t * A.cwiseAbs().maxCoeff(), t * B.cwiseAbs().maxCoeff()}); + return ((matprod(A, B) - Eigen::MatrixXd::Identity(n, n)).cwiseAbs().maxCoeff() <= t) || + ((matprod(B, A) - Eigen::MatrixXd::Identity(n, n)).cwiseAbs().maxCoeff() <= t); +} + +// Test whether matrix A has orthonormal columns up to tolerance tol +inline bool isorth(const Eigen::MatrixXd& A, double tol = 0) { + int n = A.cols(); + if (n > A.rows()) return false; + if (std::isnan(A.cwiseAbs().sum())) return false; + Eigen::MatrixXd At = A.transpose(); + Eigen::MatrixXd M = matprod(At, A) - Eigen::MatrixXd::Identity(n, n); + if (tol > 0) return M.cwiseAbs().maxCoeff() <= std::max(tol, tol * A.cwiseAbs().maxCoeff()); + return M.cwiseAbs().maxCoeff() <= 0; +} + +/* +Get a relative tolerance for a set of arrays. Borrowed from COBYQA. + +Parameters: + arrays: vector of vectors to get the tolerance for. + +Returns: + Relative tolerance for the set of arrays. + +Throws: + std::invalid_argument if no array is provided. +*/ +inline double get_arrays_tol(const std::vector& arrays) { + if (arrays.empty()) throw std::invalid_argument("At least one array must be provided."); + Eigen::Index size = 0; + double weight = 1.0; + for (const auto& arr : arrays) { + size = std::max(size, arr.size()); + double m = 0; + for (Eigen::Index i = 0; i < arr.size(); ++i) + if (std::isfinite(arr[i])) m = std::max(m, std::abs(arr[i])); + weight = std::max(weight, m); + } + if (size == 0) size = 1; + return 10.0 * EPS * static_cast(size) * weight; +} + +// Convenience overload for two vectors +inline double get_arrays_tol(const Eigen::VectorXd& a, const Eigen::VectorXd& b) { + return get_arrays_tol({a, b}); +} + +} // namespace prima + +#endif // PRIMA_CPP_LINALG_HPP diff --git a/cpp/src/prima/linear_constraints.hpp b/cpp/src/prima/linear_constraints.hpp new file mode 100644 index 0000000000..4220e6f46d --- /dev/null +++ b/cpp/src/prima/linear_constraints.hpp @@ -0,0 +1,161 @@ +#ifndef PRIMA_CPP_LINEAR_CONSTRAINTS_HPP +#define PRIMA_CPP_LINEAR_CONSTRAINTS_HPP + +/* +This module handles linear constraints of the form lb <= A*x <= ub. + +Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. + +Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +*/ + +#include +#include +#include +#include +#include + +namespace prima { + +// Represents a set of linear constraints: lb <= A*x <= ub +struct LinearConstraint { + Eigen::MatrixXd A; + Eigen::VectorXd lb; + Eigen::VectorXd ub; + + LinearConstraint() = default; + + LinearConstraint(const Eigen::MatrixXd& A_, const Eigen::VectorXd& lb_, const Eigen::VectorXd& ub_) + : A(A_), lb(lb_), ub(ub_) {} +}; + +// Combine multiple LinearConstraint objects by vertically stacking +// their A matrices and corresponding lb/ub vectors +inline LinearConstraint combine_multiple_linear_constraints( + const std::vector& constraints +) { + if (constraints.empty()) { + return LinearConstraint(); + } + + Eigen::MatrixXd full_A = constraints[0].A; + Eigen::VectorXd full_lb = constraints[0].lb; + Eigen::VectorXd full_ub = constraints[0].ub; + + for (std::size_t i = 1; i < constraints.size(); ++i) { + Eigen::MatrixXd new_A(full_A.rows() + constraints[i].A.rows(), full_A.cols()); + new_A << full_A, constraints[i].A; + full_A = new_A; + + Eigen::VectorXd new_lb(full_lb.size() + constraints[i].lb.size()); + new_lb << full_lb, constraints[i].lb; + full_lb = new_lb; + + Eigen::VectorXd new_ub(full_ub.size() + constraints[i].ub.size()); + new_ub << full_ub, constraints[i].ub; + full_ub = new_ub; + } + + return LinearConstraint(full_A, full_lb, full_ub); +} + +// Result of separating linear constraints into equality and inequality parts +struct SeparatedLinearConstraints { + Eigen::MatrixXd A_eq; + Eigen::VectorXd b_eq; + Eigen::MatrixXd A_ineq; + Eigen::VectorXd b_ineq; +}; + +/* +Separate a LinearConstraint (lb <= A*x <= ub) into: +- Equality constraints: A_eq * x = b_eq (where lb == ub within tolerance) +- Inequality constraints: A_ineq * x <= b_ineq (where lb > -inf or ub < inf) + +Lower bounds are converted to the form -A[i,:] * x <= -lb[i]. +Upper bounds are converted to the form A[i,:] * x <= ub[i]. +*/ +inline SeparatedLinearConstraints separate_LC_into_eq_and_ineq( + const LinearConstraint& linear_constraint +) { + const double eps = std::numeric_limits::epsilon(); + const double tol = 2.0 * eps; + + const Eigen::MatrixXd& A = linear_constraint.A; + const Eigen::VectorXd& lb = linear_constraint.lb; + const Eigen::VectorXd& ub = linear_constraint.ub; + + std::vector eq_indices; + std::vector ineq_lb_indices; + std::vector ineq_ub_indices; + + for (Eigen::Index i = 0; i < lb.size(); ++i) { + if (ub(i) <= lb(i) + tol) { + eq_indices.push_back(i); + } else { + if (lb(i) > -std::numeric_limits::infinity()) { + ineq_lb_indices.push_back(i); + } + if (ub(i) < std::numeric_limits::infinity()) { + ineq_ub_indices.push_back(i); + } + } + } + + SeparatedLinearConstraints result; + + if (!eq_indices.empty()) { + result.A_eq = Eigen::MatrixXd(eq_indices.size(), A.cols()); + result.b_eq = Eigen::VectorXd(eq_indices.size()); + for (std::size_t j = 0; j < eq_indices.size(); ++j) { + result.A_eq.row(j) = A.row(eq_indices[j]); + result.b_eq(j) = 0.5 * (lb(eq_indices[j]) + ub(eq_indices[j])); + } + } + + Eigen::MatrixXd A_ineq_lb, A_ineq_ub; + Eigen::VectorXd b_ineq_lb, b_ineq_ub; + bool has_lb = false, has_ub = false; + + if (!ineq_lb_indices.empty()) { + A_ineq_lb = Eigen::MatrixXd(ineq_lb_indices.size(), A.cols()); + b_ineq_lb = Eigen::VectorXd(ineq_lb_indices.size()); + for (std::size_t j = 0; j < ineq_lb_indices.size(); ++j) { + A_ineq_lb.row(j) = -A.row(ineq_lb_indices[j]); + b_ineq_lb(j) = -lb(ineq_lb_indices[j]); + } + has_lb = true; + } + + if (!ineq_ub_indices.empty()) { + A_ineq_ub = Eigen::MatrixXd(ineq_ub_indices.size(), A.cols()); + b_ineq_ub = Eigen::VectorXd(ineq_ub_indices.size()); + for (std::size_t j = 0; j < ineq_ub_indices.size(); ++j) { + A_ineq_ub.row(j) = A.row(ineq_ub_indices[j]); + b_ineq_ub(j) = ub(ineq_ub_indices[j]); + } + has_ub = true; + } + + int total_ineq = (has_lb ? A_ineq_lb.rows() : 0) + (has_ub ? A_ineq_ub.rows() : 0); + if (total_ineq > 0) { + result.A_ineq = Eigen::MatrixXd(total_ineq, A.cols()); + result.b_ineq = Eigen::VectorXd(total_ineq); + int row = 0; + if (has_lb) { + result.A_ineq.block(0, 0, A_ineq_lb.rows(), A.cols()) = A_ineq_lb; + result.b_ineq.segment(0, A_ineq_lb.rows()) = b_ineq_lb; + row += A_ineq_lb.rows(); + } + if (has_ub) { + result.A_ineq.block(row, 0, A_ineq_ub.rows(), A.cols()) = A_ineq_ub; + result.b_ineq.segment(row, A_ineq_ub.rows()) = b_ineq_ub; + } + } + + return result; +} + +} // namespace prima + +#endif // PRIMA_CPP_LINEAR_CONSTRAINTS_HPP diff --git a/cpp/src/prima/message.hpp b/cpp/src/prima/message.hpp new file mode 100644 index 0000000000..8aae51f4f4 --- /dev/null +++ b/cpp/src/prima/message.hpp @@ -0,0 +1,126 @@ +#ifndef PRIMA_CPP_MESSAGE_HPP +#define PRIMA_CPP_MESSAGE_HPP + +/* +This module provides some functions that print messages to terminal/files. + +Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. + +Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +N.B.: +1. In case parallelism is desirable (especially during initialization), the functions may + have to be modified or disabled due to the IO operations. +2. IPRINT indicates the level of verbosity, which increases with the absolute value of IPRINT. + IPRINT = +/-3 can be expensive due to high IO operations. +*/ + +#include +#include +#include +#include +#include +#include + +#include "infos.hpp" + +namespace prima { + +// Translate an info code into a human-readable reason string +inline std::string get_info_string(const std::string& solver, int info) { + std::string reason; + if (info == FTARGET_ACHIEVED) reason = "the target function value is achieved."; + else if (info == MAXFUN_REACHED) reason = "the objective function has been evaluated MAXFUN times."; + else if (info == MAXTR_REACHED) reason = "the maximal number of trust region iterations has been reached."; + else if (info == SMALL_TR_RADIUS) reason = "the trust region radius reaches its lower bound."; + else if (info == TRSUBP_FAILED) reason = "a trust region step has failed to reduce the quadratic model."; + else if (info == NAN_INF_X) reason = "NaN or Inf occurs in x."; + else if (info == NAN_INF_F) reason = "the objective function returns NaN/+Inf."; + else if (info == NAN_INF_MODEL) reason = "NaN or Inf occurs in the models."; + else if (info == DAMAGING_ROUNDING) reason = "rounding errors are becoming damaging."; + else if (info == NO_SPACE_BETWEEN_BOUNDS) reason = "there is no space between the lower and upper bounds of variable."; + else if (info == ZERO_LINEAR_CONSTRAINT) reason = "one of the linear constraints has a zero gradient"; + else if (info == CALLBACK_TERMINATE) reason = "the callback function requested termination"; + else reason = "UNKNOWN EXIT FLAG"; + return "Return from " + solver + " because " + reason; +} + +// This function prints messages at return +inline void retmsg(const std::string& solver, int info, int iprint, int nf, + double f, const Eigen::VectorXd& x, std::optional cstrv = std::nullopt, + std::optional constr = std::nullopt) { + if (std::abs(iprint) < 1) return; + + // Decide whether the problem is truly constrained + bool is_constrained = constr.has_value() || cstrv.has_value(); + // Decide the constraint violation (N.B.: We assume that the constraint is CONSTR >= 0.) + double cstrv_loc = cstrv.value_or(0); + if (!cstrv.has_value() && constr.has_value() && constr.value()) { + double m = 0; + for (Eigen::Index i = 0; i < constr.value()->size(); ++i) + m = std::max(m, -(*constr.value())[i]); + cstrv_loc = m; + } + std::string msg = get_info_string(solver, info); + + // Print X in compact form for 2 or fewer elements, otherwise as column + char buf[1024]; + if (x.size() <= 2) { + std::snprintf(buf, sizeof(buf), "\nThe corresponding X is: [%g", x[0]); + for (Eigen::Index i = 1; i < x.size(); ++i) + std::snprintf(buf + strlen(buf), sizeof(buf) - strlen(buf), ", %g", x[i]); + std::snprintf(buf + strlen(buf), sizeof(buf) - strlen(buf), "]"); + } else { + std::snprintf(buf, sizeof(buf), "\nThe corresponding X is:\n"); + for (Eigen::Index i = 0; i < x.size(); ++i) + std::snprintf(buf + strlen(buf), sizeof(buf) - strlen(buf), "%g ", x[i]); + } + msg += buf; + + if (is_constrained) + std::printf("\nNumber of function values = %d Least value of F = %g Constraint violation = %g", + nf, f, cstrv_loc); + else + std::printf("\nNumber of function values = %d Least value of F = %g", nf, f); + + std::printf("\n%s\n", msg.c_str()); +} + +// This subroutine prints messages for each evaluation of the objective function +inline void fmsg(const std::string& solver, const std::string& state, int iprint, int nf, + double delta, double f, const Eigen::VectorXd& x, + std::optional cstrv = std::nullopt, + std::optional constr = std::nullopt) { + if (std::abs(iprint) < 2) return; + bool is_constrained = constr.has_value() || cstrv.has_value(); + double cstrv_loc = cstrv.value_or(0); + std::printf("\n%s step with radius = %g", state.c_str(), delta); + if (is_constrained) + std::printf("\nNumber of function values = %d Least value of F = %g Constraint violation = %g", + nf, f, cstrv_loc); + else + std::printf("\nNumber of function values = %d Least value of F = %g", nf, f); +} + +// This function prints messages when RHO is updated +inline void rhomsg(const std::string& solver, int iprint, int nf, double f, double rho, + const Eigen::VectorXd& x, std::optional cstrv = std::nullopt, + std::optional constr = std::nullopt, + std::optional cpen = std::nullopt) { + if (std::abs(iprint) < 2) return; + bool is_constrained = constr.has_value() || cstrv.has_value(); + double cstrv_loc = cstrv.value_or(0); + if (cpen.has_value()) + std::printf("\nNew RHO = %g CPEN = %g", rho, cpen.value()); + else + std::printf("\nNew RHO = %g", rho); + if (is_constrained) + std::printf("\nNumber of function values = %d Least value of F = %g Constraint violation = %g", + nf, f, cstrv_loc); + else + std::printf("\nNumber of function values = %d Least value of F = %g", nf, f); +} + +} // namespace prima + +#endif // PRIMA_CPP_MESSAGE_HPP diff --git a/cpp/src/prima/nonlinear_constraints.hpp b/cpp/src/prima/nonlinear_constraints.hpp new file mode 100644 index 0000000000..81e3664743 --- /dev/null +++ b/cpp/src/prima/nonlinear_constraints.hpp @@ -0,0 +1,137 @@ +#ifndef PRIMA_CPP_NONLINEAR_CONSTRAINTS_HPP +#define PRIMA_CPP_NONLINEAR_CONSTRAINTS_HPP + +/* +This module handles nonlinear constraints of the form lb <= constraint_fun(x) <= ub. + +Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. + +Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +*/ + +#include +#include +#include +#include +#include +#include +#include + +namespace prima { + +// Represents a nonlinear constraint: lb <= fun(x) <= ub +struct NonlinearConstraint { + std::function fun; + Eigen::VectorXd lb; + Eigen::VectorXd ub; + + NonlinearConstraint() = default; + + NonlinearConstraint( + std::function fun_, + const Eigen::VectorXd& lb_, + const Eigen::VectorXd& ub_ + ) : fun(std::move(fun_)), lb(lb_), ub(ub_) {} +}; + +using NonlinearConstraintFunction = std::function; + +/* +Transform a NonlinearConstraint into a function that returns constraint violations. + +For each component i: + - If lb(i) > -inf: returns lb(i) - values(i) (negative if satisfied) + - If ub(i) < inf: returns values(i) - ub(i) (negative if satisfied) + +The output is a concatenated vector of all active constraint violations. +If lb or ub has size 1 but values has size > 1, the scalar bound is broadcast. +*/ +inline NonlinearConstraintFunction transform_constraint_function( + const NonlinearConstraint& nlc +) { + return [nlc](const Eigen::VectorXd& x) -> Eigen::VectorXd { + Eigen::VectorXd values = nlc.fun(x); + + Eigen::VectorXd lb = nlc.lb; + Eigen::VectorXd ub = nlc.ub; + + // Broadcast scalar bounds to match the output size + if (lb.size() == 1 && values.size() > 1) { + lb = Eigen::VectorXd::Constant(values.size(), nlc.lb(0)); + } + if (ub.size() == 1 && values.size() > 1) { + ub = Eigen::VectorXd::Constant(values.size(), nlc.ub(0)); + } + + if (values.size() != lb.size()) { + throw std::invalid_argument( + "The number of elements in the constraint function's output " + "does not match the number of elements in the lower bound." + ); + } + if (values.size() != ub.size()) { + throw std::invalid_argument( + "The number of elements in the constraint function's output " + "does not match the number of elements in the upper bound." + ); + } + + // Count result size: one entry per finite bound + int n_out = 0; + for (Eigen::Index i = 0; i < lb.size(); ++i) { + if (lb(i) > -std::numeric_limits::infinity()) ++n_out; + if (ub(i) < std::numeric_limits::infinity()) ++n_out; + } + + Eigen::VectorXd result(n_out); + int idx = 0; + for (Eigen::Index i = 0; i < lb.size(); ++i) { + if (lb(i) > -std::numeric_limits::infinity()) { + result(idx++) = lb(i) - values(i); + } + } + for (Eigen::Index i = 0; i < ub.size(); ++i) { + if (ub(i) < std::numeric_limits::infinity()) { + result(idx++) = values(i) - ub(i); + } + } + + return result; + }; +} + +/* +Combine multiple NonlinearConstraint functions into a single function. +Each constraint is transformed into a violation function, and the results +are concatenated into a single output vector. +*/ +inline NonlinearConstraintFunction process_nl_constraints( + const std::vector& nlcs +) { + std::vector functions; + functions.reserve(nlcs.size()); + for (const auto& nlc : nlcs) { + functions.push_back(transform_constraint_function(nlc)); + } + + return [functions](const Eigen::VectorXd& x) -> Eigen::VectorXd { + std::vector all_values; + int total_size = 0; + for (const auto& fun : functions) { + Eigen::VectorXd v = fun(x); + total_size += v.size(); + all_values.push_back(std::move(v)); + } + Eigen::VectorXd result(total_size); + int idx = 0; + for (const auto& v : all_values) { + result.segment(idx, v.size()) = v; + idx += v.size(); + } + return result; + }; +} + +} // namespace prima + +#endif // PRIMA_CPP_NONLINEAR_CONSTRAINTS_HPP diff --git a/cpp/src/prima/powalg.hpp b/cpp/src/prima/powalg.hpp new file mode 100644 index 0000000000..1d2fe67402 --- /dev/null +++ b/cpp/src/prima/powalg.hpp @@ -0,0 +1,162 @@ +#ifndef PRIMA_CPP_POWALG_HPP +#define PRIMA_CPP_POWALG_HPP + +// This module provides some Powell-style linear algebra procedures. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include +#include +#include + +#include "consts.hpp" +#include "linalg.hpp" + +namespace prima { + +// This function updates the QR factorization of an MxN matrix A of full column rank, attempting to +// add a new column C to this matrix as the LAST column while maintaining the full-rankness. +// Case 1. If C is not in range(A) (theoretically, it implies N < M), then the new matrix is +// [A, C] (hstack). +// Case 2. If C is in range(A), then the new matrix is [A(:, :n-1), C] (hstack). +// N.B.: +// 0. Instead of R, this subroutine updates Rdiag, which is diag(R), with a size at most M and at +// least min(m, n+1). The number is min(m, n+1) rather than min(m, n) as n may be augmented by 1 +// in the function. +// 1. With the two cases specified as above, this function does not need A as an input. +// 2. The function changes only Q[:, nsave:m] (nsave is the original value of n) and +// R[:, n-1] (n takes the updated value). +// 3. Indeed, when C is in range(A), Powell wrote in comments that "set iOUT to the index of the +// constraint (here, column of A --- Zaikun) to be deleted, but branch if no suitable index can be +// found". The idea is to replace a column of A by C so that the new matrix still has full rank +// (such a column must exist unless C = 0). But his code essentially sets iout=n always. Maybe he +// found this worked well enough in practice. Meanwhile, Powell's code includes a snippet that can +// never be reached, which was probably intended to deal with the case that IOUT != n. +inline std::tuple +qradd_Rdiag(const Eigen::VectorXd& c, Eigen::MatrixXd& Q, + const Eigen::VectorXd& Rdiag_in, int n) { + int m = Q.cols(); + int nsave = n; // Needed for debugging (only) + Eigen::VectorXd Rdiag = Rdiag_in; + + // As in Powell's COBYLA, CQ is set to 0 at the positions with CQ being negligible as per ISMINOR. + // This may not be the best choice if the subroutine is used in other contexts, e.g. LINCOA. + Eigen::VectorXd cq = matprod(c, Q); + Eigen::VectorXd ca = c.cwiseAbs(); + Eigen::MatrixXd Qa = Q.cwiseAbs(); + Eigen::VectorXd cqa = matprod(ca, Qa); + // The line below basically makes an element of cq 0 if adding it to the corresponding element of + // cqa does not change the latter. + for (int i = 0; i < cq.size(); ++i) { + if (isminor(cq[i], cqa[i])) cq[i] = 0; + } + + // Update Q so that the columns of Q[:, n+1:m] are orthogonal to C. This is done by applying a 2D + // Givens rotation to Q[:, [k, k+1]] from the right to zero C' @ Q[:, k+1] out for K=n+1, ... m-1. + // Nothing will be done if n >= m-1. + for (int k = m - 2; k >= std::max(0, n - 1); --k) { + // Powell wrote cq[k+1] != 0 instead of abs. The two differ if cq[k+1] is NaN. + // If we apply the rotation below when cq[k+1] = 0, then cq[k] will get updated to |cq[k]|. + if (std::abs(cq[k + 1]) > 0) { + Eigen::MatrixXd G = planerot(Eigen::Vector2d(cq[k], cq[k + 1])); + Eigen::MatrixXd Qk = Q.block(0, k, Q.rows(), 2) * G.transpose(); + Q.block(0, k, Q.rows(), 2) = Qk; + cq[k] = hypot(cq[k], cq[k + 1]); + } + } + + // Augment n by 1 if C is not in range(A) + if (n < m) { + // Powell's condition for the following if: cq[n+1] != 0 + if (std::abs(cq[n]) > EPS * EPS && !isminor(cq[n], cqa[n])) { + n += 1; + } + } + + // Update Rdiag so that Rdiag[n] = cq[n] = dot(c, q[:, n]). Note that N may have been augmented. + if (n - 1 >= 0 && n - 1 < m) { // n >= m should not happen unless the input is wrong + Rdiag[n - 1] = cq[n - 1]; + } + + return {Q, Rdiag, n}; +} + +// This function updates the QR factorization for an MxN matrix A=Q*R so that the updated Q and +// R form a QR factorization of [A_0, ..., A_{I-1}, A_{I+1}, ..., A_{N-1}, A_I] which is the matrix +// obtained by rearranging columns [I, I+1, ... N-1] of A to [I+1, ..., N-1, I]. Here A is ASSUMED +// TO BE OF FULL COLUMN RANK, Q is a matrix whose columns are orthogonal, and R, which is not +// present, is an upper triangular matrix whose diagonal entries are nonzero. Q and R need not be +// square. +// N.B.: +// 0. Instead of R, this function updates Rdiag, which is diag(R), the size being n. +// 1. With L = Q.cols() = R.rows(), we have M >= L >= N. Most often L = M or N. +// 2. This function changes only Q[:, i:] and Rdiag[i:]. +// 3. (NDB 20230919) In Python, i is either icon or nact - 2, whereas in FORTRAN it is either +// icon or nact - 1. +// Used in COBYLA. +inline std::tuple +qrexc_Rdiag(const Eigen::MatrixXd& A, const Eigen::MatrixXd& Q_in, + const Eigen::VectorXd& Rdiag_in, int i) { + // Sizes + int m = A.rows(), n = A.cols(); + + // Preconditions + // assert(n >= 1 && n <= m); + // assert(i >= 0 && i < n); + // assert(Rdiag.size() == n); + // assert(Q.rows() == m && Q.cols() >= n && Q.cols() <= m); + // double tol = std::max(1.0E-8, std::min(1.0E-1, 1.0E8 * EPS * m + 1)); + // assert(isorth(Q, tol)); // Costly! + + Eigen::MatrixXd Q = Q_in; + Eigen::VectorXd Rdiag = Rdiag_in; + + if (i < 0 || i >= n) return {Q, Rdiag}; + + // Let R be the upper triangular matrix in the QR factorization, namely R = Q^T * A. + // For each k, find the Givens rotation G with G * (R[k:k+2, :]) = [hypt, 0], and update + // Q[:, k:k+2] to Q[:, k:k+2] * G^T. Then R = Q^T * A is an upper triangular matrix as long + // as A[:, [k, k+1]] is updated to A[:, [k+1, k]]. Indeed, this new upper triangular matrix can + // be obtained by first updating R[[k, k+1], :] to G * (R[[k, k+1], :]) and then exchanging its + // columns K and K+1; at the same time, entries k and k+1 of R's diagonal Rdiag become + // [hypt, -(Rdiag[k+1]/hypt) * Rdiag[k]]. + // After this is done for each k = 0, ..., n-2, we obtain the QR factorization of the matrix + // that rearranges columns [i, i+1, ... n-1] of A as [i+1, ..., n-1, i]. + // Powell's code, however, is slightly different: before everything, he first exchanged columns + // k and k+1 of Q (as well as rows k and k+1 of R). This makes sure that the entries of the + // updated Rdiag are all positive if it is the case for the original Rdiag. + for (int k = i; k < n - 1; ++k) { + Eigen::Vector2d gvec(Rdiag[k + 1], inprod(Q.col(k), A.col(k + 1))); + Eigen::MatrixXd G = planerot(gvec); + Eigen::MatrixXd cols = Q.middleCols(k + 1, 1); + Eigen::MatrixXd Qpair(Q.rows(), 2); + Qpair.col(0) = Q.col(k + 1); + Qpair.col(1) = Q.col(k); + Qpair = Qpair * G.transpose(); + Q.col(k) = Qpair.col(0); + Q.col(k + 1) = Qpair.col(1); + // Powell's code updates Rdiag in the following way: + // hypt = sqrt(Rdiag[k+1]^2 + dot(Q[:, k], A[:, k+1])^2) + // Rdiag[[k, k+1]] = [hypt, (Rdiag[k+1]/hypt) * Rdiag[k]] + // Note that Rdiag[n-1] inherits all rounding in Rdiag[i:n-1] and Q[:, i:n-1] and hence + // contains significant errors. Thus we may modify Powell's code to set only + // Rdiag[k] = hypt here and then calculate Rdiag[n] by an inner product after the loop. + // Nevertheless, we simply calculate Rdiag from scratch below. + } + + // Calculate Rdiag(i:n) from scratch + for (int k = i; k < n - 1; ++k) { + Rdiag[k] = inprod(Q.col(k), A.col(k + 1)); + } + Rdiag[n - 1] = inprod(Q.col(n - 1), A.col(i)); + + return {Q, Rdiag}; +} + +} // namespace prima + +#endif // PRIMA_CPP_POWALG_HPP diff --git a/cpp/src/prima/preproc.hpp b/cpp/src/prima/preproc.hpp new file mode 100644 index 0000000000..1baebece4a --- /dev/null +++ b/cpp/src/prima/preproc.hpp @@ -0,0 +1,240 @@ +#ifndef PRIMA_CPP_PREPROC_HPP +#define PRIMA_CPP_PREPROC_HPP + +/* +This is a module that preprocesses the inputs. + +Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. + +Dedicated to late Professor M. J. D. Powell FRS (1936--2015). +*/ + +#include +#include +#include +#include +#include + +#include "consts.hpp" +#include "present.hpp" + +namespace prima { + +struct PreprocResult { + int iprint; + int maxfun; + int maxhist; + double ftarget; + double rhobeg; + double rhoend; + int npt; + int maxfilt; + double ctol; + double cweight; + double eta1; + double eta2; + double gamma1; + double gamma2; + Eigen::VectorXd x0; +}; + +/* +This subroutine preprocesses the inputs. It does nothing to the inputs that are valid. +*/ +inline PreprocResult preproc( + const std::string& solver, int num_vars, int iprint, int maxfun, int maxhist, + double ftarget, double rhobeg, double rhoend, + std::optional num_constraints = std::nullopt, + std::optional npt = std::nullopt, + std::optional maxfilt = std::nullopt, + std::optional ctol = std::nullopt, + std::optional cweight = std::nullopt, + std::optional eta1 = std::nullopt, + std::optional eta2 = std::nullopt, + std::optional gamma1 = std::nullopt, + std::optional gamma2 = std::nullopt, + std::optional is_constrained = std::nullopt, + std::optional honour_x0 = std::nullopt, + std::optional xl = std::nullopt, + std::optional xu = std::nullopt, + std::optional x0 = std::nullopt +) { + // Read num_constraints, if necessary + int nlc = (num_constraints.has_value() && solver == "cobyla") ? num_constraints.value() : 0; + + // Decide whether the problem is truly constrained + bool is_con = is_constrained.value_or(nlc > 0); + + // Validate IPRINT + if (std::abs(iprint) > 3) iprint = IPRINT_DEFAULT; + + // Validate MAXFUN + double min_maxfun; + if (solver == "uobyqa") min_maxfun = (num_vars + 1) * (num_vars + 2) / 2 + 1; + else if (solver == "cobyla") min_maxfun = num_vars + 2; + else min_maxfun = num_vars + 3; + if (maxfun < min_maxfun) maxfun = static_cast(min_maxfun); + + // Validate MAXHIST: MAXHIST > MAXFUN is never needed + if (maxhist < 0) maxhist = maxfun; + maxhist = std::min(maxhist, maxfun); + + // Validate FTARGET + if (std::isnan(ftarget)) ftarget = FTARGET_DEFAULT; + + // Validate NPT + auto npt_val = npt; + if (npt_val.has_value() && (solver == "newuoa" || solver == "bobyqa" || solver == "lincoa")) { + if (npt_val.value() < num_vars + 2 || npt_val.value() > std::min(maxfun - 1, ((num_vars + 2) * (num_vars + 1)) / 2)) { + npt_val = std::min(maxfun - 1, 2 * num_vars + 1); + } + } + + // Validate MAXFILT + auto maxfilt_val = maxfilt; + if (maxfilt_val.has_value() && (solver == "lincoa" || solver == "cobyla")) { + int mf_in = maxfilt_val.value(); + if (maxfilt_val.value() < 1) maxfilt_val = MAXFILT_DEFAULT; + else maxfilt_val = std::max(MIN_MAXFILT, maxfilt_val.value()); + // Further revise MAXFILT according to MAXHISTMEM + double unit_memo; + if (solver == "lincoa") unit_memo = (num_vars + 2) * sizeof(double); + else if (solver == "cobyla") unit_memo = (nlc + num_vars + 2) * sizeof(double); + else unit_memo = 1; + // We cannot simply set MAXFILT = MIN(MAXFILT, MAXHISTMEM/...) without type considerations + if (maxfilt_val.value() > MAXHISTMEM / unit_memo) + maxfilt_val = static_cast(MAXHISTMEM / unit_memo); + maxfilt_val = std::min(maxfun, std::max(MIN_MAXFILT, maxfilt_val.value())); + } + + // Validate ETA1 and ETA2 + double eta1_local = eta1.value_or(ETA1_DEFAULT); + double eta2_local = eta2.value_or(ETA2_DEFAULT); + + // When the difference between ETA1 and ETA2 is tiny, force them to equal + if (eta1.has_value() && eta2.has_value()) { + if (std::abs(eta1.value() - eta2.value()) < 100 * EPS * std::max(std::abs(eta1.value()), 1.0)) + eta2_local = eta1.value(); + } + + if (eta1.has_value()) { + if (std::isnan(eta1.value())) eta1_local = ETA1_DEFAULT; + else if (eta1.value() < 0 || eta1.value() >= 1) { + // Take ETA2 into account if it has a valid value + if (eta2.has_value() && eta2.value() > 0 && eta2.value() <= 1) + eta1_local = std::max(EPS, eta2.value() / 7.0); + else eta1_local = ETA1_DEFAULT; + } else eta1_local = eta1.value(); + } + + if (eta2.has_value()) { + if (std::isnan(eta2.value())) eta2_local = ETA2_DEFAULT; + else if (eta2.value() < eta1_local || eta2.value() > 1) + eta2_local = (eta1_local + 2) / 3.0; + else eta2_local = eta2.value(); + } + + // Validate GAMMA1 and GAMMA2 + double gamma1_local = gamma1.value_or(GAMMA1_DEFAULT); + double gamma2_local = gamma2.value_or(GAMMA2_DEFAULT); + + if (gamma1.has_value()) { + if (std::isnan(gamma1.value())) gamma1_local = GAMMA1_DEFAULT; + else if (gamma1.value() <= 0 || gamma1.value() >= 1) gamma1_local = GAMMA1_DEFAULT; + else gamma1_local = gamma1.value(); + } + + if (gamma2.has_value()) { + if (std::isnan(gamma2.value())) gamma2_local = GAMMA2_DEFAULT; + else if (gamma2.value() < 1 || std::isinf(gamma2.value())) gamma2_local = GAMMA2_DEFAULT; + else gamma2_local = gamma2.value(); + } + + // Validate RHOBEG and RHOEND + + // When the data is passed from interfaces (e.g., MEX) to Fortran code, RHOBEG and RHOEND + // may change a bit. If we set RHOEND = RHOBEG in the interfaces, it may happen that + // RHOEND > RHOBEG, which is invalid. Force them to equal when the difference is tiny. + if (std::abs(rhobeg - rhoend) < 100 * EPS * std::max(std::abs(rhobeg), 1.0)) + rhoend = rhobeg; + + // Revise the default values for RHOBEG/RHOEND according to the solver + double rhobeg_default = (solver == "bobyqa") + ? std::max(EPS, std::min(RHOBEG_DEFAULT, (xu.value() - xl.value()).minCoeff() / 4.0)) + : RHOBEG_DEFAULT; + double rhoend_default = (solver == "bobyqa") + ? std::max(EPS, std::min(0.1 * rhobeg_default, RHOEND_DEFAULT)) + : RHOEND_DEFAULT; + + if (solver == "bobyqa") { + // Do NOT merge this IF into the one below! Otherwise, XU and XL may be + // accessed even if the solver is not BOBYQA. + // Do NOT make this revision if RHOBEG is not positive or not finite, + // because otherwise RHOBEG will get a huge value when XU or XL contains + // huge values that indicate unbounded variables. + if (rhobeg > (xu.value() - xl.value()).minCoeff() / 2) { + rhobeg = (xu.value() - xl.value()).minCoeff() / 4.0; + } + } + if (rhobeg <= 0 || std::isnan(rhobeg) || std::isinf(rhobeg)) { + // Take RHOEND into account if it has a valid value. We do not do + // this if the solver is BOBYQA, which requires that RHOBEG <= (XU-XL)/2. + if (std::isfinite(rhoend) && rhoend > 0 && solver != "bobyqa") + rhobeg = std::max(10 * rhoend, rhobeg_default); + else rhobeg = rhobeg_default; + } + + if (rhoend <= 0 || rhobeg < rhoend || std::isnan(rhoend) || std::isinf(rhoend)) + rhoend = std::max(EPS, std::min(0.1 * rhobeg, rhoend_default)); + + // For BOBYQA, revise X0 or RHOBEG so that the distance between X0 and + // the inactive bounds is at least RHOBEG. If HONOUR_X0 == true, revise + // RHOBEG if needed; otherwise, revise X0 if needed. + Eigen::VectorXd x0_local = x0.value_or(Eigen::VectorXd()); + if (honour_x0.has_value() && solver == "bobyqa") { + if (honour_x0.value()) { + for (int i = 0; i < num_vars; ++i) { + if (std::isfinite(xl.value()[i]) && std::abs(x0_local[i] - xl.value()[i]) <= EPS * std::max(1.0, std::abs(xl.value()[i]))) + x0_local[i] = xl.value()[i]; + if (std::isfinite(xu.value()[i]) && std::abs(x0_local[i] - xu.value()[i]) <= EPS * std::max(1.0, std::abs(xu.value()[i]))) + x0_local[i] = xu.value()[i]; + } + double min_dist = rhobeg; + for (int i = 0; i < num_vars; ++i) { + if (std::isfinite(xl.value()[i]) && std::abs(x0_local[i] - xl.value()[i]) > EPS) + min_dist = std::min(min_dist, x0_local[i] - xl.value()[i]); + if (std::isfinite(xu.value()[i]) && std::abs(x0_local[i] - xu.value()[i]) > EPS) + min_dist = std::min(min_dist, xu.value()[i] - x0_local[i]); + } + rhobeg = std::max(EPS, min_dist); + } + } + + // Validate CTOL (it can be 0) + if (ctol.has_value()) { + double c = ctol.value(); + if (std::isnan(c) || c < 0) ctol = CTOL_DEFAULT; + else ctol = c; + } + + // Validate CWEIGHT (it can be +Inf) + if (cweight.has_value()) { + double cw = cweight.value(); + if (std::isnan(cw) || cw < 0) cweight = CWEIGHT_DEFAULT; + else cweight = cw; + } + + return { + iprint, maxfun, maxhist, ftarget, rhobeg, rhoend, + npt_val.value_or(0), + maxfilt_val.value_or(0), + ctol.value_or(CTOL_DEFAULT), + cweight.value_or(CWEIGHT_DEFAULT), + eta1_local, eta2_local, gamma1_local, gamma2_local, + x0_local + }; +} + +} // namespace prima + +#endif // PRIMA_CPP_PREPROC_HPP diff --git a/cpp/src/prima/present.hpp b/cpp/src/prima/present.hpp new file mode 100644 index 0000000000..6167530815 --- /dev/null +++ b/cpp/src/prima/present.hpp @@ -0,0 +1,22 @@ +#ifndef PRIMA_CPP_PRESENT_HPP +#define PRIMA_CPP_PRESENT_HPP + +/* +This is a C++ equivalent of the Fortran 'present' function for optional arguments. +*/ + +#include + +namespace prima { + +// present check for raw pointers (nullptr check) +template +inline bool present(const T* x) { return x != nullptr; } + +// present check for std::optional +template +inline bool present(const std::optional& x) { return x.has_value(); } + +} // namespace prima + +#endif // PRIMA_CPP_PRESENT_HPP diff --git a/cpp/src/prima/prima.hpp b/cpp/src/prima/prima.hpp new file mode 100644 index 0000000000..2a6524c761 --- /dev/null +++ b/cpp/src/prima/prima.hpp @@ -0,0 +1,317 @@ +#ifndef PRIMA_CPP_PRIMA_HPP +#define PRIMA_CPP_PRIMA_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "infos.hpp" +#include "bounds.hpp" +#include "project.hpp" +#include "linalg.hpp" +#include "linear_constraints.hpp" +#include "nonlinear_constraints.hpp" +#include "evaluate.hpp" +#include "message.hpp" +#include "cobyla/cobyla.hpp" + +namespace prima { + +struct MinimizeOptions { + bool quiet = true; + double rhobeg = std::numeric_limits::quiet_NaN(); + double rhoend = std::numeric_limits::quiet_NaN(); + int maxfun = 0; + int iprint = 0; + double ftarget = -std::numeric_limits::infinity(); + double ctol = std::numeric_limits::quiet_NaN(); + double cweight = CWEIGHT_DEFAULT; + double eta1 = std::numeric_limits::quiet_NaN(); + double eta2 = std::numeric_limits::quiet_NaN(); + double gamma1 = std::numeric_limits::quiet_NaN(); + double gamma2 = std::numeric_limits::quiet_NaN(); + int maxfilt = 2000; + int maxhist = 0; + std::function callback = nullptr; +}; + +struct MinimizeResult { + Eigen::VectorXd x; + bool success = false; + int status = 0; + std::string message; + double fun = 0.0; + int nfev = 0; + double maxcv = 0.0; + Eigen::VectorXd nlconstr; + std::string method; +}; + +namespace detail { + +inline LinearConstraint adjust_linear_constraint_for_fixed( + const LinearConstraint& lc, const Eigen::VectorXi& fixed_mask, + const Eigen::VectorXd& fixed_vals, int n_total +) { + int n_free = n_total - fixed_mask.sum(); + int m = lc.lb.size(); + Eigen::VectorXd contrib = Eigen::VectorXd::Zero(m); + int f_idx = 0; + for (int i = 0; i < n_total; ++i) + if (fixed_mask(i)) { contrib += lc.A.col(i) * fixed_vals(f_idx); ++f_idx; } + Eigen::MatrixXd new_A(m, n_free); + int col = 0; + for (int i = 0; i < n_total; ++i) + if (!fixed_mask(i)) { new_A.col(col) = lc.A.col(i); ++col; } + return LinearConstraint(new_A, lc.lb - contrib, lc.ub - contrib); +} + +inline NonlinearConstraintFunction wrap_nlc_for_fixed( + const NonlinearConstraintFunction& nlc, const Eigen::VectorXi& fm, + const Eigen::VectorXd& fv, int n_total +) { + return [nlc, fm, fv, n_total](const Eigen::VectorXd& x) -> Eigen::VectorXd { + Eigen::VectorXd newx(n_total); + int fi = 0, fri = 0; + for (int i = 0; i < n_total; ++i) + newx(i) = fm(i) ? fv(fi++) : x(fri++); + return nlc(newx); + }; +} + +inline Eigen::VectorXi compute_fixed_mask( + const Eigen::VectorXd& lb, const Eigen::VectorXd& ub, double tol +) { + int n = lb.size(); + Eigen::VectorXi m(n); + for (int i = 0; i < n; ++i) + m(i) = (lb(i) <= ub(i) && ub(i) <= lb(i) + tol) ? 1 : 0; + return m; +} + +template +std::optional ptr_to_opt(const T* p) { + if (p) return *p; + return std::nullopt; +} + +inline double compute_maxcv( + const Eigen::VectorXd& x, + const Eigen::MatrixXd* A_ineq, const Eigen::VectorXd* b_ineq, + const Eigen::MatrixXd* A_eq, const Eigen::VectorXd* b_eq, + const Eigen::VectorXd* nlconstr +) { + double m = 0; + if (A_ineq && A_ineq->size() > 0) { double v = (*A_ineq * x - *b_ineq).maxCoeff(); if (v > m) m = v; } + if (A_eq && A_eq->size() > 0) { double v = (*A_eq * x - *b_eq).cwiseAbs().maxCoeff(); if (v > m) m = v; } + if (nlconstr && nlconstr->size() > 0) { double v = nlconstr->maxCoeff(); if (v > m) m = v; } + return m; +} + +inline void extract_free_variables( + const Eigen::VectorXd& full, const Eigen::VectorXi& mask, Eigen::VectorXd& free +) { + int j = 0; + for (int i = 0; i < full.size(); ++i) + if (!mask(i)) free(j++) = full(i); +} + +} // namespace detail + +inline MinimizeResult minimize( + std::function fun, + const Eigen::VectorXd& x0, + const std::string& method = "", + const Bounds* bounds = nullptr, + const LinearConstraint* linear_constraint = nullptr, + const NonlinearConstraintFunction* nonlinear_constraint_function = nullptr, + MinimizeOptions options = MinimizeOptions{} +) { + using namespace Eigen; + using namespace detail; + + MinimizeResult result; + int lenx0 = x0.size(); + std::string algo = method; + bool quiet = options.quiet; + + if (algo.empty()) algo = "cobyla"; + else { + for (auto& c : algo) c = std::tolower(c); + if (algo != "cobyla") + throw std::invalid_argument("Only COBYLA is implemented"); + } + + auto [lb, ub] = process_bounds(bounds, lenx0); + double arr_tol = get_arrays_tol(lb, ub); + VectorXi fixed_mask = compute_fixed_mask(lb, ub, arr_tol); + int n_fixed = fixed_mask.sum(); + bool any_fixed = n_fixed > 0; + bool all_fixed = (n_fixed == lenx0); + + std::function active_fun(std::move(fun)); + + if (all_fixed) { + VectorXd fv(n_fixed); + for (int i = 0; i < lenx0; ++i) + fv(i) = std::max(lb(i), std::min(ub(i), 0.5 * (lb(i) + ub(i)))); + double fval = active_fun(fv); + VectorXd nlconstr; + if (nonlinear_constraint_function) nlconstr = (*nonlinear_constraint_function)(fv); + MatrixXd A_eq, A_ineq; VectorXd b_eq, b_ineq; + if (linear_constraint) { + auto sep = separate_LC_into_eq_and_ineq(*linear_constraint); + if (sep.A_eq.size() > 0) { A_eq = sep.A_eq; b_eq = sep.b_eq; } + if (sep.A_ineq.size() > 0) { A_ineq = sep.A_ineq; b_ineq = sep.b_ineq; } + } + double maxcv = compute_maxcv(fv, + A_ineq.size() > 0 ? &A_ineq : nullptr, b_ineq.size() > 0 ? &b_ineq : nullptr, + A_eq.size() > 0 ? &A_eq : nullptr, b_eq.size() > 0 ? &b_eq : nullptr, + nlconstr.size() > 0 ? &nlconstr : nullptr); + result.x = fv; result.success = true; result.status = FIXED_SUCCESS; + result.message = "All variables were fixed by the provided bounds."; + result.fun = fval; result.nfev = 1; result.maxcv = maxcv; + result.nlconstr = nlconstr; result.method = "cobyla"; + return result; + } + + VectorXd x0_adj = x0, lb_adj = lb, ub_adj = ub, fixed_vals; + const LinearConstraint* active_lc = linear_constraint; + LinearConstraint adjusted_lc; + NonlinearConstraintFunction active_nlc_fun; + const NonlinearConstraintFunction* active_nlc = nonlinear_constraint_function; + + if (any_fixed) { + fixed_vals = VectorXd(n_fixed); + { + int j = 0; + for (int i = 0; i < lenx0; ++i) + if (fixed_mask(i)) + fixed_vals(j++) = std::max(lb(i), std::min(ub(i), 0.5 * (lb(i) + ub(i)))); + } + int n_free = lenx0 - n_fixed; + x0_adj = VectorXd(n_free); lb_adj = VectorXd(n_free); ub_adj = VectorXd(n_free); + extract_free_variables(x0, fixed_mask, x0_adj); + extract_free_variables(lb, fixed_mask, lb_adj); + extract_free_variables(ub, fixed_mask, ub_adj); + + VectorXi fm = fixed_mask; VectorXd fv = fixed_vals; + int ntot = lenx0; + auto orig = std::move(active_fun); + active_fun = [orig, fm, fv, ntot](const VectorXd& x) -> double { + VectorXd newx(ntot); int fi = 0, fri = 0; + for (int i = 0; i < ntot; ++i) + newx(i) = fm(i) ? fv(fi++) : x(fri++); + return orig(newx); + }; + if (linear_constraint) { + adjusted_lc = adjust_linear_constraint_for_fixed(*linear_constraint, fixed_mask, fixed_vals, lenx0); + active_lc = &adjusted_lc; + } + if (nonlinear_constraint_function) { + active_nlc_fun = wrap_nlc_for_fixed(*nonlinear_constraint_function, fixed_mask, fixed_vals, lenx0); + active_nlc = &active_nlc_fun; + } + } + + if (active_nlc == nullptr) { + auto proj = project(x0_adj, lb_adj, ub_adj, active_lc); + x0_adj = proj.x; + } + + // Separate linear constraints + MatrixXd A_eq, A_ineq; VectorXd b_eq, b_ineq; + if (active_lc) { + auto sep = separate_LC_into_eq_and_ineq(*active_lc); + if (sep.A_eq.rows() > 0) { A_eq = sep.A_eq; b_eq = sep.b_eq; } + if (sep.A_ineq.rows() > 0) { A_ineq = sep.A_ineq; b_ineq = sep.b_ineq; } + } + + // Evaluate f0 and nlconstr0 + double f0_val = active_fun(x0_adj); + VectorXd nlconstr0; + int m_nlcon = 0; + if (active_nlc) { + nlconstr0 = (*active_nlc)(x0_adj); + m_nlcon = nlconstr0.size(); + } + + // Build calcfc like pyprima does: returns (f, nlconstr) + // Capture active_fun and active_nlc by value to avoid dangling references + auto fun_copy = active_fun; + auto nlc_copy = active_nlc ? std::make_shared(*active_nlc) : nullptr; + Calcfc calcfc_inner; + if (nlc_copy) { + calcfc_inner = [fun_copy, nlc_copy](const VectorXd& x) -> std::pair { + return {fun_copy(x), (*nlc_copy)(x)}; + }; + } else { + calcfc_inner = [fun_copy](const VectorXd& x) -> std::pair { + return {fun_copy(x), VectorXd(0)}; + }; + } + + // Wrap bounds + linear constraints into AMAT/BVEC via get_lincon + auto [amat_opt, bvec_opt] = get_lincon( + ptr_to_opt(A_eq.size() > 0 ? &A_eq : nullptr), + ptr_to_opt(A_ineq.size() > 0 ? &A_ineq : nullptr), + ptr_to_opt(b_eq.size() > 0 ? &b_eq : nullptr), + ptr_to_opt(b_ineq.size() > 0 ? &b_ineq : nullptr), + ptr_to_opt(&lb_adj), + ptr_to_opt(&ub_adj) + ); + + auto cobyla_result = cobyla( + calcfc_inner, m_nlcon, x0_adj, + amat_opt, bvec_opt, std::nullopt, std::nullopt, + std::nullopt, std::nullopt, + f0_val, nlconstr0.size() > 0 ? std::optional(nlconstr0) : std::nullopt, + std::isnan(options.rhobeg) ? std::nullopt : std::optional(options.rhobeg), + std::isnan(options.rhoend) ? std::nullopt : std::optional(options.rhoend), + options.ftarget, + std::isnan(options.ctol) ? CTOL_DEFAULT : options.ctol, + options.cweight, + options.maxfun > 0 ? std::optional(options.maxfun) : std::nullopt, + options.iprint, + std::isnan(options.eta1) ? std::nullopt : std::optional(options.eta1), + std::isnan(options.eta2) ? std::nullopt : std::optional(options.eta2), + GAMMA1_DEFAULT, GAMMA2_DEFAULT, + options.maxhist > 0 ? std::optional(options.maxhist) : std::nullopt, + options.maxfilt, + options.callback + ); + + result.x = cobyla_result.x; + result.fun = cobyla_result.f; + result.nlconstr = cobyla_result.constr; + result.maxcv = cobyla_result.cstrv; + result.nfev = cobyla_result.nf; + result.status = cobyla_result.info; + result.success = (cobyla_result.info == SMALL_TR_RADIUS || cobyla_result.info == FTARGET_ACHIEVED); + result.message = get_info_string("COBYLA", cobyla_result.info); + result.method = "cobyla"; + + if (any_fixed) { + VectorXd full_x(lenx0); + int fi = 0, fri = 0; + for (int i = 0; i < lenx0; ++i) + full_x(i) = fixed_mask(i) ? fixed_vals(fi++) : result.x(fri++); + result.x = full_x; + } + + return result; +} + +} // namespace prima + +#endif // PRIMA_CPP_PRIMA_HPP diff --git a/cpp/src/prima/project.hpp b/cpp/src/prima/project.hpp new file mode 100644 index 0000000000..4963f8fdf6 --- /dev/null +++ b/cpp/src/prima/project.hpp @@ -0,0 +1,138 @@ +#ifndef PRIMA_CPP_COMMON_HPP +#define PRIMA_CPP_COMMON_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include "bounds.hpp" +#include "linear_constraints.hpp" +#include "infos.hpp" + +namespace prima { + +struct OptimizeResult { + Eigen::VectorXd x; + bool success = false; + int status = 0; + std::string message; + double fun = 0.0; + int nfev = 0; + double maxcv = 0.0; + Eigen::VectorXd nlconstr; + std::string method; +}; + +inline OptimizeResult project( + const Eigen::VectorXd& x0, + const Eigen::VectorXd& lb, + const Eigen::VectorXd& ub, + const LinearConstraint* linear_constraint +) { + const double max_con = 1e20; + const double eps = std::numeric_limits::epsilon(); + + OptimizeResult result; + + if (linear_constraint == nullptr) { + Eigen::VectorXd x_proj = x0.cwiseMax(lb).cwiseMin(ub); + result.x = x_proj; + return result; + } + + const Eigen::MatrixXd& A = linear_constraint->A; + const Eigen::VectorXd& lc_lb = linear_constraint->lb; + const Eigen::VectorXd& lc_ub = linear_constraint->ub; + + bool all_eq = true; + for (Eigen::Index i = 0; i < lc_lb.size(); ++i) { + if (std::abs(lc_ub(i) - lc_lb(i)) > eps) { + all_eq = false; + break; + } + } + bool bounds_trivial = (lb.maxCoeff() <= -max_con) && (ub.minCoeff() >= max_con); + + if (all_eq && bounds_trivial) { + Eigen::VectorXd b = 0.5 * (lc_lb + lc_ub); + Eigen::VectorXd rhs = b - A * x0; +#if EIGEN_VERSION_AT_LEAST(5,0,0) + Eigen::VectorXd xi = A.bdcSvd().solve(rhs); +#else + Eigen::VectorXd xi = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(rhs); +#endif + Eigen::VectorXd x_proj = (x0 + xi).cwiseMax(lb).cwiseMin(ub); + result.x = x_proj; + return result; + } + + if (linear_constraint != nullptr) { + Eigen::VectorXd ax = A * x0; + bool feasible = true; + + for (Eigen::Index i = 0; i < lc_lb.size(); ++i) { + if (ax(i) < lc_lb(i) - eps || ax(i) > lc_ub(i) + eps) { + feasible = false; + break; + } + } + for (Eigen::Index i = 0; i < x0.size(); ++i) { + if (x0(i) < lb(i) - eps || x0(i) > ub(i) + eps) { + feasible = false; + break; + } + } + + if (feasible) { + result.x = x0; + return result; + } + + Eigen::VectorXd x_proj = x0.cwiseMax(lb).cwiseMin(ub); + + Eigen::MatrixXd Aeq; + Eigen::VectorXd beq; + bool has_eq = false; + { + std::vector eq_idx; + for (Eigen::Index i = 0; i < lc_lb.size(); ++i) { + if (std::abs(lc_ub(i) - lc_lb(i)) <= eps) { + eq_idx.push_back(i); + } + } + if (!eq_idx.empty()) { + Aeq = Eigen::MatrixXd(eq_idx.size(), A.cols()); + beq = Eigen::VectorXd(eq_idx.size()); + for (std::size_t j = 0; j < eq_idx.size(); ++j) { + Aeq.row(j) = A.row(eq_idx[j]); + beq(j) = 0.5 * (lc_lb(eq_idx[j]) + lc_ub(eq_idx[j])); + } + has_eq = true; + } + } + + if (has_eq) { + Eigen::VectorXd rhs = beq - Aeq * x_proj; +#if EIGEN_VERSION_AT_LEAST(5,0,0) + Eigen::VectorXd step = Aeq.bdcSvd().solve(rhs); +#else + Eigen::VectorXd step = Aeq.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(rhs); +#endif + x_proj = (x_proj + step).cwiseMax(lb).cwiseMin(ub); + } + + result.x = x_proj; + return result; + } + + result.x = x0; + return result; +} + +} // namespace prima + +#endif // PRIMA_CPP_COMMON_HPP diff --git a/cpp/src/prima/ratio.hpp b/cpp/src/prima/ratio.hpp new file mode 100644 index 0000000000..6364259380 --- /dev/null +++ b/cpp/src/prima/ratio.hpp @@ -0,0 +1,44 @@ +#ifndef PRIMA_CPP_RATIO_HPP +#define PRIMA_CPP_RATIO_HPP + +// This module calculates the reduction ratio for trust-region methods. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include + +#include "consts.hpp" + +namespace prima { + +// This function evaluates the reduction ratio of a trust-region step, handling inf/nan properly. +inline double redrat(double ared, double pred, double rshrink) { + if (std::isnan(ared)) { + // This should not happen in unconstrained problems due to the moderated extreme barrier. + return -REALMAX; + } + if (std::isnan(pred) || pred <= 0) { + // The trust-region subproblem solver fails in this rare case. Instead of terminating as + // Powell's original code does, we set ratio as follows so that the solver may continue + // to progress. + if (ared > 0) { + // The trial point will be accepted, but the trust-region radius will be shrunk if + // rshrink > 0. + return rshrink / 2; + } + // Set the ratio to a large negative number to signify a bad trust-region step, so that the + // solver will check whether to take a geometry step or reduce rho. + return -REALMAX; + } + if (std::isinf(pred) && std::isinf(ared)) return 1; // ared/pred = NaN if calculated directly + if (std::isinf(pred) && std::isinf(-ared)) return -REALMAX; // ared/pred = NaN if calculated directly + return ared / pred; +} + +} // namespace prima + +#endif // PRIMA_CPP_RATIO_HPP diff --git a/cpp/src/prima/redrho.hpp b/cpp/src/prima/redrho.hpp new file mode 100644 index 0000000000..87249217d8 --- /dev/null +++ b/cpp/src/prima/redrho.hpp @@ -0,0 +1,30 @@ +#ifndef PRIMA_CPP_REDRHO_HPP +#define PRIMA_CPP_REDRHO_HPP + +// This module provides a function that calculates RHO when it needs to be reduced. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include + +namespace prima { + +// This function calculates RHO when it needs to be reduced. +// The scheme is shared by UOBYQA, NEWUOA, BOBYQA, LINCOA. For COBYLA, Powell's code reduces RHO by +// 'RHO *= 0.5; if RHO <= 1.5 * RHOEND: RHO = RHOEND' as specified in (11) of the COBYLA +// paper. However, this scheme seems to work better, especially after we introduce DELTA. +inline double redrho(double rho_in, double rhoend) { + double rho_ratio = rho_in / rhoend; + double rho; + if (rho_ratio > 250) rho = 0.1 * rho_in; + else if (rho_ratio <= 16) rho = rhoend; + else rho = std::sqrt(rho_ratio) * rhoend; // rho = sqrt(rho * rhoend) + return rho; +} + +} // namespace prima + +#endif // PRIMA_CPP_REDRHO_HPP diff --git a/cpp/src/prima/selectx.hpp b/cpp/src/prima/selectx.hpp new file mode 100644 index 0000000000..4d7cb66ada --- /dev/null +++ b/cpp/src/prima/selectx.hpp @@ -0,0 +1,237 @@ +#ifndef PRIMA_CPP_SELECTX_HPP +#define PRIMA_CPP_SELECTX_HPP + +// This module provides subroutines that ensure the returned X is optimal among all the calculated +// points in the sense that no other point achieves both lower function value and lower constraint +// violation at the same time. This module is needed only in the constrained case. +// +// Translated from Zaikun Zhang's modern-Fortran reference implementation in PRIMA. +// +// Dedicated to late Professor M. J. D. Powell FRS (1936--2015). + +#include +#include +#include +#include +#include + +#include "consts.hpp" +#include "present.hpp" +#include "linalg.hpp" + +namespace prima { + +// This function compares whether FC1 = (F1, C1) is (strictly) better than FC2 = (F2, C2), which +// basically means that (F1 < F2 and C1 <= C2) or (F1 <= F2 and C1 < C2). +// It takes care of the cases where some of these values are NaN or Inf, even though some cases +// should never happen due to the moderated extreme barrier. +// At return, returns TRUE if and only if (F1, C1) is better than (F2, C2). +// Here, C means constraint violation, which is a nonnegative number. +inline bool isbetter(double f1, double c1, double f2, double c2, double ctol) { + if (std::isnan(f1) || std::isinf(f1) || std::isnan(c1) || std::isinf(c1)) return false; + if (std::isnan(f2) || std::isinf(f2) || std::isnan(c2) || std::isinf(c2)) return true; + + // Even though NaN/+Inf should not occur in FC1 or FC2 due to the moderated extreme barrier, for + // security and robustness, the code below does not make this assumption. + bool ib = false; + ib = ib || (f1 < f2 && c1 <= c2); + ib = ib || (f1 <= f2 && c1 < c2); + + // If C1 <= CTOL and C2 is significantly larger/worse than CTOL, i.e., C2 > MAX(CTOL,CREF), + // then FC1 is better than FC2 as long as F1 < REALMAX. Normally CREF >= CTOL so MAX(CTOL, CREF) + // is indeed CREF. However, this may not be true if CTOL > 1E-1*CONSTRMAX. + double cref = 10 * std::max(EPS, std::min(ctol, 1.0E-2 * CONSTRMAX)); // The MIN avoids overflow. + ib = ib || (f1 < REALMAX && c1 <= ctol && (c2 > std::max(ctol, cref) || std::isnan(c2))); + return ib; +} + +// This subroutine saves X, F, and CSTRV in XFILT, FFILT, and CFILT (and CONSTR in CONFILT +// if they are present), unless a vector in XFILT[:, :NFILT] is better than X. +// If X is better than some vectors in XFILT[:, :NFILT] then these vectors will be +// removed. If X is not better than any of XFILT[:, :NFILT], but NFILT == MAXFILT, +// then we remove a column from XFILT according to the merit function +// PHI = FFILT + CWEIGHT * max(CFILT - CTOL, 0) +// N.B.: +// 1. Only XFILT[:, :NFILT] and FFILT[:, :NFILT] etc contains valid information, +// while XFILT[:, NFILT+1:MAXFILT] and FFILT[:, NFILT+1:MAXFILT] etc are not +// initialized yet. +// 2. We decide whether X is better than another by the ISBETTER function. +inline int savefilt(double cstrv, double ctol, double cweight, double f, + const Eigen::VectorXd& x, int nfilt, + Eigen::VectorXd& cfilt, Eigen::VectorXd& ffilt, + Eigen::MatrixXd& xfilt, + const Eigen::VectorXd* constr, Eigen::MatrixXd* confilt) { + int maxfilt = ffilt.size(); + + // Return immediately if any column of XFILT is better than X. + for (int i = 0; i < nfilt; ++i) { + if (isbetter(ffilt[i], cfilt[i], f, cstrv, ctol)) return nfilt; + if (ffilt[i] <= f && cfilt[i] <= cstrv) return nfilt; + } + + // Decide which columns of XFILT to keep. + Eigen::VectorXi keep = Eigen::VectorXi::Ones(nfilt); + for (int i = 0; i < nfilt; ++i) { + if (isbetter(f, cstrv, ffilt[i], cfilt[i], ctol)) keep[i] = 0; + } + + // If NFILT == MAXFILT and X is not better than any column of XFILT, then we remove the worst + // column of XFILT according to the merit function PHI = FFILT + CWEIGHT * MAX(CFILT - CTOL, 0). + int n_keep = keep.sum(); + if (n_keep == maxfilt) { // In this case, nfilt == keep.size() == n_keep == maxfilt > 0. + Eigen::VectorXd cfilt_shifted = (cfilt.array() - ctol).cwiseMax(0); + Eigen::VectorXd phi; + if (cweight <= 0) phi = ffilt; + else if (std::isinf(cweight)) { + phi = cfilt_shifted; + // We should not use CFILT here; if MAX(CFILT_SHIFTED) is attained at multiple indices, + // then we will check FFILT to exhaust the remaining degree of freedom. + } else { + phi = ffilt.array().max(-REALMAX).matrix(); + for (int i = 0; i < phi.size(); ++i) + if (std::isnan(phi[i])) phi[i] = -REALMAX; // Replace NaN with -REALMAX. + phi.array() += cweight * cfilt_shifted.array(); + } + // We select X to maximize PHI. In case there are multiple maximizers, we take the one with + // the largest CSTRV_SHIFTED; if there are more than one choices, we take the one with the + // largest F; if there are several candidates, we take the one with the largest CSTRV; if + // the last comparison still leads to more than one possibilities, then they are equally bad + // and we choose the first. + // N.B.: + // 1. This process is the opposite of selecting KOPT in SELECTX. + // 2. In finite-precision arithmetic, PHI_1 == PHI_2 and CSTRV_SHIFTED_1 == + // CSTRV_SHIFTED_2 do not ensure that F_1 == F_2! + double phimax = phi.maxCoeff(); + double cref = 0; + for (int i = 0; i < nfilt; ++i) + if (phi[i] >= phimax) cref = std::max(cref, cfilt_shifted[i]); + double fref = -std::numeric_limits::infinity(); + for (int i = 0; i < nfilt; ++i) + if (cfilt_shifted[i] >= cref) fref = std::max(fref, ffilt[i]); + int kworst = 0; + for (int i = 0; i < nfilt; ++i) + if (ffilt[i] <= fref && cfilt[i] > cfilt[kworst]) kworst = i; + if (kworst >= 0 && kworst < nfilt) keep[kworst] = 0; + n_keep = keep.sum(); + } + + // Keep the good xfilt values and remove all the ones that are strictly worse than the new x. + int idx = 0; + for (int i = 0; i < nfilt; ++i) { + if (keep[i]) { + xfilt.col(idx) = xfilt.col(i); + ffilt[idx] = ffilt[i]; + cfilt[idx] = cfilt[i]; + if (confilt && constr) confilt->col(idx) = confilt->col(i); + ++idx; + } + } + nfilt = idx; + + // Once we have removed all the vectors that are strictly worse than x, + // we add x to the filter. + xfilt.col(nfilt) = x; + ffilt[nfilt] = f; + cfilt[nfilt] = cstrv; + if (confilt && constr) confilt->col(nfilt) = *constr; + ++nfilt; + + return nfilt; +} + +// This subroutine selects X according to FHIST and CHIST, which represents (a part of) history +// of F and CSTRV. Normally, FHIST and CHIST are not the full history but only a filter, e.g. ffilt +// and CFILT generated by SAVEFILT. However, we name them as FHIST and CHIST because the [F, CSTRV] +// in a filter should not dominate each other, but this subroutine does NOT assume such a property. +// N.B.: CTOL is the tolerance of the constraint violation (CSTRV). A point is considered feasible if +// its constraint violation is at most CTOL. Note that CTOL is absolute, not relative. +inline int selectx(const Eigen::VectorXd& fhist, const Eigen::VectorXd& chist, + double cweight, double ctol) { + int nhist = fhist.size(); + + double fref, cref; + // We select X among the points with F < FREF and CSTRV < CREF. + // Do NOT use F <= FREF, because F == FREF (FUNCMAX or REALMAX) may mean F == INF in practice! + bool has_fc = false; + for (int i = 0; i < nhist; ++i) { + if (fhist[i] < FUNCMAX && chist[i] < CONSTRMAX) { has_fc = true; break; } + } + if (has_fc) { fref = FUNCMAX; cref = CONSTRMAX; } + else { + bool has_fr = false; + for (int i = 0; i < nhist; ++i) + if (fhist[i] < REALMAX && chist[i] < CONSTRMAX) { has_fr = true; break; } + if (has_fr) { fref = REALMAX; cref = CONSTRMAX; } + else { + bool has_fc2 = false; + for (int i = 0; i < nhist; ++i) + if (fhist[i] < FUNCMAX && chist[i] < REALMAX) { has_fc2 = true; break; } + if (has_fc2) { fref = FUNCMAX; cref = REALMAX; } + else { fref = REALMAX; cref = REALMAX; } + } + } + + int kopt = nhist - 1; + bool any_valid = false; + for (int i = 0; i < nhist; ++i) + if (fhist[i] < fref && chist[i] < cref) { any_valid = true; break; } + + if (any_valid) { + // Shift the constraint violations by ctol, so that cstrv <= ctol is regarded as no violation. + Eigen::VectorXd chist_shifted = (chist.array() - ctol).cwiseMax(0); + // cmin is the minimal shift constraint violation attained in the history. + double cmin = std::numeric_limits::infinity(); + for (int i = 0; i < nhist; ++i) + if (fhist[i] < fref) cmin = std::min(cmin, chist_shifted[i]); + // We consider only the points whose shifted constraint violations are at most the cref below. + // N.B.: Without taking std::max(EPS, .), cref would be 0 if cmin = 0. In that case, asking for + // cstrv_shift < cref would be WRONG! + double cref2 = std::max(EPS, 2 * cmin); + + // We use the following phi as our merit function to select X. + Eigen::VectorXd phi(nhist); + if (cweight <= 0) phi = fhist; + else if (std::isinf(cweight)) { + phi = chist_shifted; + // We should not use chist here; if MIN(chist_shifted) is attained at multiple indices, + // then we will check fhist to exhaust the remaining degree of freedom. + } else { + phi = fhist.array().max(-REALMAX) + cweight * chist_shifted.array(); + // std::max(fhist, -REALMAX) makes sure that phi will not contain NaN (unless there is a bug). + } + + // We select X to minimize phi subject to f < fref and cstrv_shift <= cref (see the comments + // above for the reason of taking "<" and "<=" in these two constraints). In case there are + // multiple minimizers, we take the one with the least cstrv_shift; if there is more than one + // choice, we take the one with the least f; if there are several candidates, we take the one + // with the least cstrv; if the last comparison still leads to more than one possibility, then + // they are equally good and we choose the first. + // N.B.: + // 1. This process is the opposite of selecting kworst in savefilt. + // 2. In finite-precision arithmetic, phi_1 == phi_2 and cstrv_shift_1 == cstrv_shift_2 do + // not ensure that f_1 == f_2! + double phimin = std::numeric_limits::infinity(); + for (int i = 0; i < nhist; ++i) + if (fhist[i] < fref && chist_shifted[i] <= cref2) + phimin = std::min(phimin, phi[i]); + + double cmin2 = std::numeric_limits::infinity(); + for (int i = 0; i < nhist; ++i) + if (fhist[i] < fref && phi[i] <= phimin) + cmin2 = std::min(cmin2, chist_shifted[i]); + + double fmin2 = std::numeric_limits::infinity(); + for (int i = 0; i < nhist; ++i) + if (chist_shifted[i] <= cmin2) + fmin2 = std::min(fmin2, fhist[i]); + + for (int i = 0; i < nhist; ++i) + if (fhist[i] <= fmin2 && chist[i] < chist[kopt]) kopt = i; + } + + return kopt; +} + +} // namespace prima + +#endif // PRIMA_CPP_SELECTX_HPP diff --git a/cpp/src/primacpp-config.cmake.in b/cpp/src/primacpp-config.cmake.in new file mode 100644 index 0000000000..8e2b3dfff4 --- /dev/null +++ b/cpp/src/primacpp-config.cmake.in @@ -0,0 +1,6 @@ +@PACKAGE_INIT@ + +include(CMakeFindDependencyMacro) +find_dependency(Eigen3) + +include("${CMAKE_CURRENT_LIST_DIR}/primacpp-targets.cmake") diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt new file mode 100644 index 0000000000..120db2ab8d --- /dev/null +++ b/cpp/tests/CMakeLists.txt @@ -0,0 +1,20 @@ +option (PRIMA_ENABLE_TESTING "Build tests" OFF) + +if (NOT TARGET tests) + add_custom_target (tests) +endif () + +macro (prima_add_cpp_test name) + add_executable(${name}_cpp_exe EXCLUDE_FROM_ALL ${name}.cpp) + if (PRIMA_ENABLE_TESTING) + set_target_properties (${name}_cpp_exe PROPERTIES EXCLUDE_FROM_ALL FALSE) + endif () + target_link_libraries(${name}_cpp_exe PRIVATE primacpp) + add_dependencies(tests ${name}_cpp_exe) + add_test(NAME ${name}_cpp COMMAND ${name}_cpp_exe) +endmacro () + +prima_add_cpp_test(test_bounds) +prima_add_cpp_test(test_get_lincon) +prima_add_cpp_test(test_end_to_end) +prima_add_cpp_test(test_miscellaneous) diff --git a/cpp/tests/test_bounds.cpp b/cpp/tests/test_bounds.cpp new file mode 100644 index 0000000000..7925867012 --- /dev/null +++ b/cpp/tests/test_bounds.cpp @@ -0,0 +1,137 @@ +#include +#include +#include +#include + +#include "prima/infos.hpp" +#include "prima/bounds.hpp" +#include "prima/project.hpp" +#include "prima/linear_constraints.hpp" +#include "prima/nonlinear_constraints.hpp" +#include "prima/linalg.hpp" +#include "prima/prima.hpp" +#include "test_common.hpp" + +using namespace prima; +using namespace Eigen; +using prima::test::TestRunner; + +void test_process_bounds(TestRunner& t) { + auto [lb, ub] = process_bounds(nullptr, 3); + t.check(lb.size() == 3, "process_bounds(nullptr, 3) size"); + t.check(std::isinf(-lb(0)) && lb(0) < 0, "process_bounds nullptr lb = -inf"); + t.check(std::isinf(ub(0)) && ub(0) > 0, "process_bounds nullptr ub = inf"); + + VectorXd lb_in(2); lb_in << 1.0, 2.0; + VectorXd ub_in(2); ub_in << 3.0, 4.0; + Bounds b(lb_in, ub_in); + auto [lb2, ub2] = process_bounds(&b, 3); + t.check(lb2.size() == 3, "process_bounds bounds size"); + t.check(std::abs(lb2(0) - 1.0) < 1e-15, "process_bounds lb[0]"); + t.check(std::abs(lb2(1) - 2.0) < 1e-15, "process_bounds lb[1]"); + t.check(std::isinf(-lb2(2)) && lb2(2) < 0, "process_bounds lb[2] padded -inf"); + t.check(std::abs(ub2(0) - 3.0) < 1e-15, "process_bounds ub[0]"); + t.check(std::abs(ub2(1) - 4.0) < 1e-15, "process_bounds ub[1]"); + t.check(std::isinf(ub2(2)) && ub2(2) > 0, "process_bounds ub[2] padded inf"); +} + +void test_eliminate_fixed_bounds(TestRunner& t) { + auto obj = [](const VectorXd& x) -> double { return x.squaredNorm(); }; + + VectorXd lb(5), ub(5); + lb << -1, -std::numeric_limits::infinity(), 1, + -std::numeric_limits::infinity(), -0.5; + ub << -0.5, -0.5, std::numeric_limits::infinity(), + std::numeric_limits::infinity(), -0.5; + Bounds bounds(lb, ub); + VectorXd x0(5); x0 << 1, 2, 3, 4, 5; + + MinimizeOptions opts; + opts.quiet = true; + auto result = minimize(obj, x0, "cobyla", &bounds, nullptr, nullptr, opts); + t.check(std::abs(result.x(0) + 0.5) < 1e-3, "x[0] == -0.5"); + t.check(std::abs(result.x(1) + 0.5) < 1e-3, "x[1] == -0.5"); + t.check(std::abs(result.x(2) - 1.0) < 1e-3, "x[2] == 1"); + t.check(std::abs(result.x(3)) < 5e-3, "x[3] == 0"); + t.check(std::abs(result.x(4) + 0.5) < 1e-3, "x[4] == -0.5"); + t.check(std::abs(result.fun - 1.75) < 1e-3, "f == 1.75"); +} + +void test_eliminate_fixed_bounds_with_linear_constraints(TestRunner& t) { + auto obj = [](const VectorXd& x) -> double { return x.squaredNorm(); }; + + VectorXd lb(3), ub(3); + lb << -1, -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(); + ub << -1, std::numeric_limits::infinity(), + std::numeric_limits::infinity(); + Bounds bounds(lb, ub); + MatrixXd A(1, 3); A << 1, 1, 1; + LinearConstraint lc(A, VectorXd::Constant(1, 9), VectorXd::Constant(1, 15)); + VectorXd x0(3); x0 << 1, 1, 1; + + MinimizeOptions opts; + opts.quiet = true; + auto result = minimize(obj, x0, "cobyla", &bounds, &lc, nullptr, opts); + t.check(std::abs(result.x(0) + 1.0) < 1e-6, "x[0] == -1"); + t.check(std::abs(result.x(1) - 5.0) < 5e-5, "x[1] == 5"); + t.check(std::abs(result.x(2) - 5.0) < 5e-5, "x[2] == 5"); + t.check(std::abs(result.fun - 51.0) < 1e-6, "f == 51"); +} + +void test_eliminate_fixed_bounds_with_nonlinear_constraints(TestRunner& t) { + auto obj = [](const VectorXd& x) -> double { return x.squaredNorm(); }; + + VectorXd lb(3), ub(3); + lb << -1, -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(); + ub << -1, std::numeric_limits::infinity(), + std::numeric_limits::infinity(); + Bounds bounds(lb, ub); + + int n = 3; + auto nlc_fun = [n](const VectorXd& x) -> VectorXd { + return VectorXd::Constant(1, x(n - 1) * x(n - 1)); + }; + NonlinearConstraint nlc(nlc_fun, VectorXd::Constant(1, 9), VectorXd::Constant(1, 15)); + auto nlc_func = transform_constraint_function(nlc); + VectorXd x0(3); x0 << 1, 1, 1; + + MinimizeOptions opts; + opts.quiet = true; + auto result = minimize(obj, x0, "cobyla", &bounds, nullptr, &nlc_func, opts); + t.check(std::abs(result.x(0) + 1.0) < 1e-6, "x[0] == -1"); + t.check(std::abs(result.x(1)) < 1e-6, "x[1] == 0"); + t.check(std::abs(result.x(2) - 3.0) < 1e-6, "x[2] == 3"); + t.check(std::abs(result.fun - 10.0) < 1e-6, "f == 10"); +} + +void test_all_fixed(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + Bounds bounds(Vector2d(1.0, 2.0), Vector2d(1.0, 2.0)); + auto obj = [](const VectorXd& x) -> double { + return std::pow(x(0) - 1.0, 2) + std::pow(x(1) - 2.0, 2); + }; + + MinimizeOptions opts; + opts.quiet = true; + auto result = minimize(obj, x0, "cobyla", &bounds, nullptr, nullptr, opts); + t.check(result.success, "all fixed success"); + t.check(result.status == FIXED_SUCCESS, "status FIXED_SUCCESS"); + t.check(std::abs(result.x(0) - 1.0) < 1e-15, "x[0] == 1"); + t.check(std::abs(result.x(1) - 2.0) < 1e-15, "x[1] == 2"); + t.check(result.nfev == 1, "nfev == 1"); +} + +int main() { + std::cout << "=== test_bounds: Bounds and Fixed Variables ===" << std::endl; + TestRunner t; + + test_process_bounds(t); + test_eliminate_fixed_bounds(t); + test_eliminate_fixed_bounds_with_linear_constraints(t); + test_eliminate_fixed_bounds_with_nonlinear_constraints(t); + test_all_fixed(t); + + return t.done(); +} diff --git a/cpp/tests/test_common.hpp b/cpp/tests/test_common.hpp new file mode 100644 index 0000000000..666eab7013 --- /dev/null +++ b/cpp/tests/test_common.hpp @@ -0,0 +1,34 @@ +#ifndef PRIMA_CPP_TEST_COMMON_HPP +#define PRIMA_CPP_TEST_COMMON_HPP + +#include +#include +#include + +namespace prima::test { + +struct TestRunner { + int failures = 0; + + void check(bool cond, const std::string& msg) { + if (!cond) { + std::cerr << "FAIL: " << msg << std::endl; + ++failures; + } else { + std::cout << "PASS: " << msg << std::endl; + } + } + + int done() { + if (failures > 0) { + std::cerr << "\n" << failures << " test(s) FAILED!" << std::endl; + return 1; + } + std::cout << "\nAll tests PASSED." << std::endl; + return 0; + } +}; + +} // namespace prima::test + +#endif // PRIMA_CPP_TEST_COMMON_HPP diff --git a/cpp/tests/test_end_to_end.cpp b/cpp/tests/test_end_to_end.cpp new file mode 100644 index 0000000000..d51713df07 --- /dev/null +++ b/cpp/tests/test_end_to_end.cpp @@ -0,0 +1,99 @@ +#include +#include +#include + +#include "prima/infos.hpp" +#include "prima/bounds.hpp" +#include "prima/linear_constraints.hpp" +#include "prima/nonlinear_constraints.hpp" +#include "prima/prima.hpp" +#include "test_common.hpp" + +using namespace prima; +using namespace Eigen; +using prima::test::TestRunner; + +static double obj(const VectorXd& x) { + return std::pow(x(0) - 1.0, 2) + std::pow(x(1) - 2.5, 2); +} + +void test_no_constraints(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + MinimizeOptions opts; + opts.quiet = true; + opts.rhoend = 1e-4; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "no constraints success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 2.5) < 1e-3, "x[1] close to 2.5"); +} + +void test_bounds(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + Bounds bounds(Vector2d(-5.0, 10.0), Vector2d(5.0, 10.0)); + MinimizeOptions opts; + opts.quiet = true; + opts.rhoend = 1e-4; + auto result = minimize(obj, x0, "cobyla", &bounds, nullptr, nullptr, opts); + t.check(result.success, "bounds success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 10.0) < 1e-3, "x[1] close to 10"); +} + +void test_linear_constraints(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + // x0 + x1 = 5 + MatrixXd A(1, 2); A << 1.0, 1.0; + LinearConstraint lc(A, VectorXd::Constant(1, 5.0), VectorXd::Constant(1, 5.0)); + MinimizeOptions opts; + opts.quiet = true; + opts.rhoend = 1e-4; + auto result = minimize(obj, x0, "cobyla", nullptr, &lc, nullptr, opts); + t.check(result.success, "linear constraints success"); + t.check(std::abs(result.x(0) - 1.75) < 1.0, "x[0] close to 1.75"); + t.check(std::abs(result.x(1) - 3.25) < 1.0, "x[1] close to 3.25"); + t.check(std::abs(result.x(0) + result.x(1) - 5.0) < 1e-3, "eq satisfied"); +} + +void test_nonlinear_constraint(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + + auto cons_fun = [](const VectorXd& x) -> VectorXd { + return VectorXd::Constant(1, x(0) * x(0) - 0.25); + }; + NonlinearConstraint nlc(cons_fun, VectorXd::Constant(1, 0.0), VectorXd::Constant(1, 0.0)); + auto nlc_func = transform_constraint_function(nlc); + + MinimizeOptions opts; + opts.quiet = true; + opts.rhoend = 1e-4; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, &nlc_func, opts); + t.check(result.success, "nonlinear constraint success"); + t.check(std::abs(result.x(0) - 0.5) < 1e-3, "x[0] close to 0.5"); + t.check(std::abs(result.x(1) - 2.5) < 0.05, "x[1] close to 2.5"); + t.check(std::abs(result.x(0) * result.x(0) - 0.25) < 1e-8, "cons satisfied"); +} + +void test_auto_detect(TestRunner& t) { + VectorXd x0(2); x0 << 3.0, 3.0; + MinimizeOptions opts; + opts.quiet = true; + opts.rhoend = 1e-4; + opts.maxfun = 100; + auto result = minimize(obj, x0, "", nullptr, nullptr, nullptr, opts); + t.check(result.success, "auto-detect success"); + t.check(result.method == "cobyla", "method == cobyla"); +} + +int main() { + std::cout << "=== test_end_to_end: Solver Integration ===" << std::endl; + TestRunner t; + + test_no_constraints(t); + test_bounds(t); + test_linear_constraints(t); + test_nonlinear_constraint(t); + test_auto_detect(t); + + return t.done(); +} diff --git a/cpp/tests/test_get_lincon.cpp b/cpp/tests/test_get_lincon.cpp new file mode 100644 index 0000000000..8a00b74d16 --- /dev/null +++ b/cpp/tests/test_get_lincon.cpp @@ -0,0 +1,103 @@ +#include +#include + +#include "prima/bounds.hpp" +#include "prima/linear_constraints.hpp" +#include "prima/cobyla/cobyla.hpp" +#include "test_common.hpp" + +using namespace prima; +using namespace Eigen; +using prima::test::TestRunner; + +void test_get_lincon(TestRunner& t) { + MatrixXd Aeq(2, 2); Aeq << 1, 2, 3, 4; + MatrixXd Aineq(2, 2); Aineq << 5, 6, 7, 8; + VectorXd beq(2); beq << 9, 10; + VectorXd bineq(2); bineq << 11, 12; + VectorXd xl(2); xl << 0, -1; + VectorXd xu(2); xu << 13, 14; + + auto [amat_opt, bvec_opt] = get_lincon(Aeq, Aineq, beq, bineq, xl, xu); + t.check(amat_opt.has_value(), "amat has value"); + t.check(bvec_opt.has_value(), "bvec has value"); + + MatrixXd amat_expected(10, 2); + amat_expected << -1, 0, 0, -1, 1, 0, 0, 1, -1, -2, -3, -4, 1, 2, 3, 4, 5, 6, 7, 8; + VectorXd bvec_expected(10); + bvec_expected << 0, 1, 13, 14, -9, -10, 9, 10, 11, 12; + + t.check(amat_opt->isApprox(amat_expected), "amat matches"); + t.check(bvec_opt->isApprox(bvec_expected), "bvec matches"); +} + +void test_get_lincon_boundmax(TestRunner& t) { + MatrixXd Aeq(2, 2); Aeq << 1, 2, 3, 4; + MatrixXd Aineq(2, 2); Aineq << 5, 6, 7, 8; + VectorXd beq(2); beq << 9, 10; + VectorXd bineq(2); bineq << 11, 12; + VectorXd xl(2); xl << -BOUNDMAX - 1, -1; + VectorXd xu(2); xu << 13, 14; + + auto [amat_opt, bvec_opt] = get_lincon(Aeq, Aineq, beq, bineq, xl, xu); + t.check(amat_opt.has_value(), "boundmax amat has value"); + t.check(bvec_opt.has_value(), "boundmax bvec has value"); + + MatrixXd amat_expected(9, 2); + amat_expected << 0, -1, 1, 0, 0, 1, -1, -2, -3, -4, 1, 2, 3, 4, 5, 6, 7, 8; + VectorXd bvec_expected(9); + bvec_expected << 1, 13, 14, -9, -10, 9, 10, 11, 12; + + t.check(amat_opt->isApprox(amat_expected), "boundmax amat matches"); + t.check(bvec_opt->isApprox(bvec_expected), "boundmax bvec matches"); +} + +void test_get_lincon_none(TestRunner& t) { + auto [amat, bvec] = get_lincon(std::nullopt, std::nullopt, + std::nullopt, std::nullopt, + std::nullopt, std::nullopt); + t.check(!amat.has_value(), "none amat is nullopt"); + t.check(!bvec.has_value(), "none bvec is nullopt"); +} + +void test_combine_linear_constraints(TestRunner& t) { + MatrixXd A1(1, 2); A1 << 1.0, 0.0; + MatrixXd A2(1, 2); A2 << 0.0, 1.0; + LinearConstraint lc1(A1, VectorXd::Constant(1, -1.0), VectorXd::Constant(1, 1.0)); + LinearConstraint lc2(A2, VectorXd::Constant(1, -2.0), VectorXd::Constant(1, 2.0)); + + auto combined = combine_multiple_linear_constraints({lc1, lc2}); + t.check(combined.A.rows() == 2, "combined A rows == 2"); + t.check(combined.A.cols() == 2, "combined A cols == 2"); + t.check(std::abs(combined.A(0, 0) - 1.0) < 1e-15, "A(0,0) == 1"); + t.check(std::abs(combined.A(1, 1) - 1.0) < 1e-15, "A(1,1) == 1"); + t.check(std::abs(combined.lb(1) + 2.0) < 1e-15, "lb(1) == -2"); + t.check(std::abs(combined.ub(1) - 2.0) < 1e-15, "ub(1) == 2"); +} + +void test_separate_LC(TestRunner& t) { + MatrixXd A(2, 2); A << 1.0, 0.0, 0.0, 1.0; + LinearConstraint lc(A, Vector2d(2.0, -1.0), Vector2d(2.0, 5.0)); + auto sep = separate_LC_into_eq_and_ineq(lc); + t.check(sep.A_eq.rows() == 1, "eq rows == 1"); + t.check(sep.A_ineq.rows() == 2, "ineq rows == 2"); + t.check(std::abs(sep.b_eq(0) - 2.0) < 1e-15, "b_eq == 2"); + + auto sep2 = separate_LC_into_eq_and_ineq( + LinearConstraint(A, Vector2d(-1.0, -2.0), Vector2d(3.0, 4.0))); + t.check(sep2.A_eq.rows() == 0, "no eq rows"); + t.check(sep2.A_ineq.rows() == 4, "4 ineq rows"); +} + +int main() { + std::cout << "=== test_get_lincon: Linear Constraint Helpers ===" << std::endl; + TestRunner t; + + test_get_lincon(t); + test_get_lincon_boundmax(t); + test_get_lincon_none(t); + test_combine_linear_constraints(t); + test_separate_LC(t); + + return t.done(); +} diff --git a/cpp/tests/test_miscellaneous.cpp b/cpp/tests/test_miscellaneous.cpp new file mode 100644 index 0000000000..f432531c32 --- /dev/null +++ b/cpp/tests/test_miscellaneous.cpp @@ -0,0 +1,168 @@ +#include +#include +#include +#include + +#include "prima/infos.hpp" +#include "prima/nonlinear_constraints.hpp" +#include "prima/prima.hpp" +#include "test_common.hpp" + +using namespace prima; +using namespace Eigen; +using prima::test::TestRunner; + +static double obj(const VectorXd& x) { + return std::pow(x(0) - 1.0, 2) + std::pow(x(1) - 2.5, 2); +} + +void test_callback_terminate(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + int call_count = 0; + auto cb = [&](const VectorXd&, double, int, int, double, const VectorXd&) -> bool { + ++call_count; return true; + }; + MinimizeOptions opts; + opts.quiet = true; + opts.callback = cb; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(call_count >= 1, "callback was called"); + t.check(result.fun < 100, "non-catastrophic result"); +} + +void test_callback_no_terminate(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + int call_count = 0; + auto cb = [&](const VectorXd&, double, int, int, double, const VectorXd&) -> bool { + ++call_count; return false; + }; + MinimizeOptions opts; + opts.quiet = true; + opts.rhoend = 1e-4; + opts.callback = cb; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "no-terminate callback success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 2.5) < 1e-3, "x[1] close to 2.5"); + t.check(call_count > 1, "callback called multiple times"); +} + +void test_rhoend_without_rhobeg(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + MinimizeOptions opts; + opts.quiet = true; + opts.rhoend = 4e-4; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "rhoend only success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 2.5) < 1e-3, "x[1] close to 2.5"); +} + +void test_rhobeg_without_rhoend(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + MinimizeOptions opts; + opts.quiet = true; + opts.rhobeg = -1; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "rhobeg only success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 2.5) < 1e-3, "x[1] close to 2.5"); +} + +void test_eta2_without_eta1(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + MinimizeOptions opts; + opts.quiet = true; + opts.eta2 = 0.7; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "eta2 only success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 2.5) < 1e-3, "x[1] close to 2.5"); +} + +void test_eta2_out_of_range(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + MinimizeOptions opts; + opts.quiet = true; + opts.eta2 = 1.7; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "eta2 out of range success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 2.5) < 1e-3, "x[1] close to 2.5"); +} + +void test_eta1_out_of_range(TestRunner& t) { + VectorXd x0(2); x0 << 5.0, 5.0; + MinimizeOptions opts; + opts.quiet = true; + opts.eta1 = 1.7; + auto result = minimize(obj, x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "eta1 out of range success"); + t.check(std::abs(result.x(0) - 1.0) < 1e-3, "x[0] close to 1"); + t.check(std::abs(result.x(1) - 2.5) < 1e-3, "x[1] close to 2.5"); +} + +void test_method_validation(TestRunner& t) { + VectorXd x0(2); x0 << 0.0, 0.0; + + bool threw = false; + try { minimize(obj, x0, "invalid"); } + catch (const std::invalid_argument&) { threw = true; } + t.check(threw, "invalid method throws"); + + threw = false; + try { minimize(obj, x0, "newuoa"); } + catch (const std::invalid_argument&) { threw = true; } + t.check(threw, "newuoa throws (only cobyla implemented)"); +} + +void test_minimize_constraint_violation(TestRunner& t) { + auto combined_nlc = [](const VectorXd& x) -> VectorXd { + VectorXd r(2); + r(0) = x(0) - 4; + r(1) = 5 - x(0); + return r; + }; + VectorXd lb2(2); + lb2 << -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(); + VectorXd ub2(2); ub2 << 0, 0; + NonlinearConstraint nlc(combined_nlc, lb2, ub2); + auto nlc_func = transform_constraint_function(nlc); + + VectorXd x0(1); x0 << 0; + MinimizeOptions opts; + opts.quiet = true; + auto result = minimize([](const VectorXd& x) { return x(0); }, + x0, "cobyla", nullptr, nullptr, &nlc_func, opts); + t.check(result.maxcv > 0.1, "constraint violation > 0.1"); +} + +void test_scalar(TestRunner& t) { + VectorXd x0(1); x0 << 5; + MinimizeOptions opts; + opts.quiet = true; + auto result = minimize([](const VectorXd& x) { return x(0) * x(0); }, + x0, "cobyla", nullptr, nullptr, nullptr, opts); + t.check(result.success, "scalar success"); + t.check(std::abs(result.x(0)) < 1e-3, "x close to 0"); + t.check(std::abs(result.fun) < 1e-3, "f close to 0"); +} + +int main() { + std::cout << "=== test_miscellaneous: Options and Edge Cases ===" << std::endl; + TestRunner t; + + test_callback_terminate(t); + test_callback_no_terminate(t); + test_rhoend_without_rhobeg(t); + test_rhobeg_without_rhoend(t); + test_eta2_without_eta1(t); + test_eta2_out_of_range(t); + test_eta1_out_of_range(t); + test_method_validation(t); + test_minimize_constraint_violation(t); + test_scalar(t); + + return t.done(); +}