Write a wrapper for CMA-ES and test in 2d

optimization/optimizers/CMakeLists.txt

+add_subdirectory(cma_es)
 add_subdirectory(conjugate_gradient_descent)
 add_subdirectory(gradient_descent)
 add_subdirectory(nelder_mead)

optimization/optimizers/cma_es/CMakeLists.txt

+if (VISUALIZE)
+    add_executable(cma_es
+        main.cpp
+    )
+    target_link_libraries(cma_es optimization_lib cma-es clara)
+    #make_plot_target(cma_es 2d ARGS -d 2)
+    #make_plot_target(cma_es 10d ARGS -d 10)
+endif()

optimization/optimizers/cma_es/cma_es.h

+#ifndef OPTIMIZATION_CMA_ES_H
+#define OPTIMIZATION_CMA_ES_H
+
+#include "utils/vector.h"
+#include "cmaes_interface.h"
+#include <iostream>
+#include <iomanip>
+
+namespace optimization {
+
+//--------------------------------------------------------------------------------------------------
+class CmaES {
+public:
+    CmaES() {}
+    CmaES(uint32_t d, uint32_t ps, int32_t s, uint32_t me, uint32_t mi, Scalar ryt) :
+        dim_(d),
+        pop_size_(ps),
+        seed_(s),
+        max_evaulations_(me),
+        max_iterations_(mi),
+        rel_y_tol_(ryt)
+    {}
+
+    uint32_t n_evaluations() const { return n_evaluations_; }
+    uint32_t n_iterations() const { return n_iterations_; }
+
+    template <typename Objective>
+    VectorXs optimize(
+        Objective& objective,
+        VectorXs const& initial_point,
+        VectorXs const& initial_std_dev
+    );
+
+private:
+    uint32_t dim_;
+    uint32_t pop_size_;
+    int32_t seed_;
+    uint32_t max_evaulations_;
+    uint32_t max_iterations_;
+    Scalar rel_y_tol_;
+
+    uint32_t n_evaluations_;
+    uint32_t n_iterations_;
+};
+
+//..................................................................................................
+template <typename Objective>
+VectorXs CmaES::optimize(
+    Objective& objective,
+    VectorXs const& initial_point,
+    VectorXs const& initial_std_dev
+) {
+    n_evaluations_ = 0;
+    n_iterations_ = 0;
+
+    // Decalare the CMA-ES object, and the arrays we will use to communicate with it.
+    cmaes_t cma_es;
+    double* const* points;
+    double* values;
+
+    // For some reason cmaes_init_para takes double*, not double const*, so I have to copy these.
+    // TODO: Looking at Hansen's code, there's no reason they aren't const. I should submit a PR.
+    VectorXs point = initial_point;
+    VectorXs std_dev = initial_std_dev;
+
+    cmaes_init_para(
+        &cma_es,
+        dim_,
+        point.data(),
+        std_dev.data(),
+        seed_,
+        pop_size_,
+        "none"
+    );
+    cma_es.sp.stopMaxFunEvals = max_evaulations_;
+    cma_es.sp.stopMaxIter = max_iterations_;
+    cma_es.sp.stopTolFun = rel_y_tol_;
+    values = cmaes_init_final(&cma_es);
+
+    std::cout << cmaes_SayHello(&cma_es) << '\n';
+
+    while(!cmaes_TestForTermination(&cma_es)) {
+        // Get a generation of points.
+        points = cmaes_SamplePopulation(&cma_es);
+
+        // TODO: This wants to be parallelized.
+        for (uint32_t i = 0; i < pop_size_; ++i) {
+            auto point_map = Eigen::Map<const Eigen::MatrixXd>(points[i], dim_, 1);
+            // TODO: make objectives handle Maps so we don't have to copy the data.
+            point = point_map;
+            objective.eval(point, values[i]);
+        }
+
+        // Update the search distribution used for cmaes_SamplePopulation().
+        cmaes_UpdateDistribution(&cma_es, values);
+        ++n_iterations_;
+    }
+
+    std::cout << cmaes_TestForTermination(&cma_es) << '\n';
+    //cmaes_WriteToFile(&cma_es, "all", "allcmaes.dat");  // write final results
+
+    // Get best estimator for the optimum, xmean.
+    double const* final_point_data = cmaes_GetPtr(&cma_es, "xmean");
+    auto point_map = Eigen::Map<const Eigen::MatrixXd>(final_point_data, dim_, 1);
+    point = point_map;
+    n_evaluations_ = cmaes_Get(&cma_es, "eval");
+    cmaes_exit(&cma_es);
+
+    return point;
+}
+
+}
+
+#endif

optimization/optimizers/cma_es/main.cpp

+#include "clara.hpp"
+#include "cma_es.h"
+#include "objectives/paraboloid.h"
+#include "objectives/rosenbrock.h"
+#include "utils/eigen_json.h"
+#include <iostream>
+#include <fstream>
+
+using namespace optimization;
+using json = nlohmann::json;
+
+//--------------------------------------------------------------------------------------------------
+int main(int const argc, char const** argv) {
+    std::string log_file_path;
+    std::string vis_file_path;
+    uint32_t dim = 2;
+    uint32_t pop_size = 10;
+    int32_t seed = 0xdefceed9;
+    uint32_t max_evaluations = 1000;
+    uint32_t max_iterations = 1000;
+    Scalar rel_y_tol = 1e-8;
+    Scalar x0 = -1;
+    Scalar y0 = 2;
+    Scalar std_dev = 1;
+
+    auto const cli =
+        clara::Arg(log_file_path, "log_file_path")("Location of the optimization log") |
+        clara::Arg(vis_file_path, "vis_file_path")("Location of the visualization file") |
+        clara::Opt(dim, "dim")["-d"]["--dimension"]("Dimensionality of the objective") |
+        clara::Opt(pop_size, "pop_size")["-p"]["--pop-size"]("Population size") |
+        clara::Opt(seed, "seed")["-s"]["--seed"]("Seed for CMA-ES random number generator") |
+        clara::Opt(max_evaluations, "max_evaluations")["-n"]["--max-evaluations"]("Max number of function evaluations") |
+        clara::Opt(max_iterations, "max_iterations")["-n"]["--max-iterations"]("Max number of generations") |
+        clara::Opt(rel_y_tol, "rel_y_tol")["-t"]["--rel-y-tol"]("Termination condition") |
+        clara::Opt(x0, "x0")["-x"]["--x0"]("X coordinate of initial point") |
+        clara::Opt(y0, "y0")["-y"]["--y0"]("Y coordinate of initial point");
+        clara::Opt(std_dev, "std_dev")["-d"]["--std-dev"]("Initial standard deviation");
+    auto const result = cli.parse(clara::Args(argc, argv));
+    if (!result) {
+        std::cerr << "Error in command line: " << result.errorMessage() << std::endl;
+        exit(1);
+    }
+
+    VectorXs initial_point;
+    initial_point.resize(dim);
+    initial_point[0] = x0;
+    initial_point[1] = y0;
+    for (uint32_t i = 2; i < dim; ++i) {
+        initial_point[i] = -1;
+    }
+
+    VectorXs initial_std_dev;
+    initial_std_dev.resize(dim);
+    initial_std_dev.fill(std_dev);
+
+    //using Objective = Paraboloid<Vector2<Scalar>>;
+    //Objective objective(dim);
+    using Objective = Rosenbrock<-1>;
+    Objective objective;
+    objective.dim() = dim;
+
+    CmaES optimizer(dim, pop_size, seed, max_evaluations, max_iterations, rel_y_tol);
+    VectorXs minimum = optimizer.optimize(objective, initial_point, initial_std_dev);
+    std::cout << "n evaluations: " << optimizer.n_evaluations() << '\n';
+    std::cout << "n generations: " << optimizer.n_iterations() << '\n';
+    std::cout << "final point: " << minimum << '\n';
+
+    /*
+    if (!log_file_path.empty()) {
+        json data = optimizer;
+        std::ofstream log_file(log_file_path);
+        log_file << data.dump(4) << '\n';
+    }
+
+    if (!vis_file_path.empty()) {
+        json data = GradientDescentVis<-1>{optimizer};
+        std::ofstream vis_file(vis_file_path);
+        vis_file << data.dump(4) << '\n';
+    }
+    */
+
+    return 0;
+}