diff --git a/optimization/optimizers/CMakeLists.txt b/optimization/optimizers/CMakeLists.txt
index bd410e1e04c085e98eecb25403a0d81f50545f2f..c585fe9762a67535cd384d04adbaab1bea5e7760 100644
--- a/optimization/optimizers/CMakeLists.txt
+++ b/optimization/optimizers/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory(cma_es)
add_subdirectory(conjugate_gradient_descent)
add_subdirectory(gradient_descent)
add_subdirectory(nelder_mead)
diff --git a/optimization/optimizers/cma_es/CMakeLists.txt b/optimization/optimizers/cma_es/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..94a71510805d26d469a2244e1c3df98d71dad872
--- /dev/null
+++ b/optimization/optimizers/cma_es/CMakeLists.txt
@@ -0,0 +1,8 @@
+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()
diff --git a/optimization/optimizers/cma_es/cma_es.h b/optimization/optimizers/cma_es/cma_es.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ce85d8645b2811c49cf0779470915a0b5d77162
--- /dev/null
+++ b/optimization/optimizers/cma_es/cma_es.h
@@ -0,0 +1,114 @@
+#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
diff --git a/optimization/optimizers/cma_es/main.cpp b/optimization/optimizers/cma_es/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e629fa7ae3f6ab9455d575a0f594b3ab55b0e6d5
--- /dev/null
+++ b/optimization/optimizers/cma_es/main.cpp
@@ -0,0 +1,83 @@
+#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;
+}