Add C++ implementation of COBYLA

see #136

Author: jschueller

.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
+
+      - name: Test
+        run: ctest --test-dir build --output-on-failure --timeout 100 --schedule-random -V -j3

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)

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)

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`.

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)

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 <iostream>
+#include <cmath>
+#include <Eigen/Core>
+#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<double>::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;
+    }
+}

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 <iostream>
+#include <cmath>
+#include <iomanip>
+#include <Eigen/Core>
+#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<double>::infinity(),
+                           std::numeric_limits<double>::infinity(),
+                           std::numeric_limits<double>::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<double>::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<double>::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<double>::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<double>::infinity(),
+                 -std::numeric_limits<double>::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<double>::infinity()),
+        Vector3d(0.5, 0.25, std::numeric_limits<double>::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;
+}

cpp/src/CMakeLists.txt

@@ -0,0 +1,48 @@
+add_library (primacpp INTERFACE)
+target_include_directories (primacpp INTERFACE
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+target_link_libraries (primacpp INTERFACE Eigen3::Eigen)
+target_compile_options (primacpp INTERFACE $<$<CXX_COMPILER_ID:MSVC>:/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
+)