diff --git a/optimization/optimizers/cma_es/CMakeLists.txt b/optimization/optimizers/cma_es/CMakeLists.txt index 1eefff0d475009deecd255459725582ed4d99b8e..4e4b2650bb313760299289a6e3c82c59a5b0742f 100644 --- a/optimization/optimizers/cma_es/CMakeLists.txt +++ b/optimization/optimizers/cma_es/CMakeLists.txt @@ -3,7 +3,5 @@ add_executable(cma_es ) target_link_libraries(cma_es optimization_lib cma-es clara) -if (VISUALIZE) - make_plot_target(cma_es 2d ARGS -d 2) - make_plot_target(cma_es 10d ARGS -d 10 -n 10000) -endif() +make_plot_target(cma_es 2d ARGS -d 2) +make_plot_target(cma_es 10d ARGS -d 10 -n 10000) diff --git a/optimization/optimizers/cma_es/cma_es.h b/optimization/optimizers/cma_es/cma_es.h index 8fb5e6fc0a4924124d84394f22d6e59caba1549c..dbf8f123e84e00ee4a203991a6543d05579cc9f9 100644 --- a/optimization/optimizers/cma_es/cma_es.h +++ b/optimization/optimizers/cma_es/cma_es.h @@ -1,50 +1,70 @@ #ifndef OPTIMIZATION_CMA_ES_H #define OPTIMIZATION_CMA_ES_H -#include "cma_es_log.h" -#include "utils/vector.h" #include "cmaes_interface.h" +#include "logs/nothing.h" +#include "utils/vector.h" #include <iostream> -#include <iomanip> namespace optimization { //-------------------------------------------------------------------------------------------------- -class CmaEs : public CmaEsLog { +class CmaEs { public: - CmaEs() {} + CmaEs(): points_(nullptr), values_(nullptr) {} 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) + rel_y_tol_(ryt), + points_(nullptr), + values_(nullptr) {} - - uint32_t n_evaluations() const { return n_evaluations_; } - uint32_t n_iterations() const { return n_iterations_; } + ~CmaEs() { if (points_ != nullptr) { cmaes_exit(&cma_es_); } } uint32_t dim() const { return dim_; } uint32_t& dim() { return dim_; } uint32_t pop_size() const { return pop_size_; } uint32_t& pop_size() { return pop_size_; } + uint32_t n_evaluations() { return cmaes_Get(&cma_es_, "eval"); } + uint32_t n_iterations() const { return n_iterations_; } + // TODO See if I can make these functions const by modifying some CMA-ES guts. + inline Eigen::Map<const Eigen::MatrixXd> point(); + Scalar value() { return cmaes_Get(&cma_es_, "fbestever"); } + template <typename Objective> - VectorXs optimize( + Eigen::Map<const Eigen::MatrixXd> optimize( Objective& objective, VectorXs const& initial_point, Scalar initial_std_dev = 0.3 ); + template <typename Objective, typename Log> + Eigen::Map<const Eigen::MatrixXd> optimize( + Objective& objective, + VectorXs const& initial_point, + Scalar initial_std_dev, + Log& log + ); template <typename Objective> - VectorXs optimize( + Eigen::Map<const Eigen::MatrixXd> optimize( Objective& objective, VectorXs const& initial_point, VectorXs const& initial_std_dev ); + template <typename Objective, typename Log> + Eigen::Map<const Eigen::MatrixXd> optimize( + Objective& objective, + VectorXs const& initial_point, + VectorXs const& initial_std_dev, + Log& log + ); private: + // hyperparameters uint32_t dim_; uint32_t pop_size_; int32_t seed_; @@ -52,90 +72,119 @@ private: uint32_t max_iterations_; Scalar rel_y_tol_; - uint32_t n_evaluations_; + // algorithm state uint32_t n_iterations_; + cmaes_t cma_es_; + // These are used to communicate with the CMA-ES object. They also help with memory management, + // in that cma_es_ has memory that eventually needs to be freed if and only if points_ and + // values_ are nullptr. + double* const* points_; // points to cma_es_'s population points when we're optimizing + double* values_; // points to cma_es_'s population values when we're optimizing }; +//.................................................................................................. +Eigen::Map<const Eigen::MatrixXd> CmaEs::point() { + double const* best_point = cmaes_GetPtr(&cma_es_, "xbestever"); + return Eigen::Map<const Eigen::MatrixXd>(best_point, dim_, 1); +} + //.................................................................................................. template <typename Objective> -VectorXs CmaEs::optimize( +Eigen::Map<const Eigen::MatrixXd> CmaEs::optimize( Objective& objective, VectorXs const& initial_point, Scalar initial_std_dev +) { + CmaEsLogNothing log; + return optimize(objective, initial_point, initial_std_dev, log); +} + +//.................................................................................................. +template <typename Objective, typename Log> +Eigen::Map<const Eigen::MatrixXd> CmaEs::optimize( + Objective& objective, + VectorXs const& initial_point, + Scalar initial_std_dev, + Log& log ) { VectorXs std_dev_vector; std_dev_vector.resize(initial_point.size()); std_dev_vector.fill(initial_std_dev); - return optimize(objective, initial_point, std_dev_vector); + return optimize(objective, initial_point, std_dev_vector, log); } //.................................................................................................. template <typename Objective> -VectorXs CmaEs::optimize( +Eigen::Map<const Eigen::MatrixXd> CmaEs::optimize( Objective& objective, VectorXs const& initial_point, VectorXs const& initial_std_dev ) { - n_evaluations_ = 0; + CmaEsLogNothing log; + return optimize(objective, initial_point, initial_std_dev, log); +} + +//.................................................................................................. +template <typename Objective, typename Log> +Eigen::Map<const Eigen::MatrixXd> CmaEs::optimize( + Objective& objective, + VectorXs const& initial_point, + VectorXs const& initial_std_dev, + Log& log +) { n_iterations_ = 0; + if (points_ != nullptr) { + cmaes_exit(&cma_es_); + } - // 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); + // This is used for various copy operations. Shouldn't be necessary. + VectorXs point_vec = initial_point; + { + // 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 std_dev = initial_std_dev; + + cmaes_init_para( + &cma_es_, + dim_, + point_vec.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_); - CmaEsLog::initialize(objective); - std::cout << cmaes_SayHello(&cma_es) << '\n'; + log.initialize(objective); + std::cout << cmaes_SayHello(&cma_es_) << '\n'; - while(!cmaes_TestForTermination(&cma_es)) { + while (!cmaes_TestForTermination(&cma_es_)) { // Get a generation of points. - points = cmaes_SamplePopulation(&cma_es); + 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); + 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]); + point_vec = point_map; + objective.eval(point_vec, values_[i]); } - CmaEsLog::push_back(points, values, dim_, pop_size_); + log.push_back(points_, values_, dim_, pop_size_); // Update the search distribution used for cmaes_SamplePopulation(). - cmaes_UpdateDistribution(&cma_es, values); + cmaes_UpdateDistribution(&cma_es_, values_); ++n_iterations_; } - std::cout << cmaes_TestForTermination(&cma_es) << '\n'; - //cmaes_WriteToFile(&cma_es, "all", "allcmaes.dat"); // write final results + std::cout << cmaes_TestForTermination(&cma_es_) << '\n'; // 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; + return point(); } } diff --git a/optimization/optimizers/cma_es/cma_es_log.h b/optimization/optimizers/cma_es/logs/everything.h similarity index 64% rename from optimization/optimizers/cma_es/cma_es_log.h rename to optimization/optimizers/cma_es/logs/everything.h index 00aff539f390fe47689ff08675512c71095065de..de81575163a64c9d3ad2fe872c789466f4e79c1c 100644 --- a/optimization/optimizers/cma_es/cma_es_log.h +++ b/optimization/optimizers/cma_es/logs/everything.h @@ -1,15 +1,10 @@ -#ifndef OPTIMIZATION_CMA_ES_LOG_H -#define OPTIMIZATION_CMA_ES_LOG_H +#ifndef OPTIMIZATION_CMA_ES_LOGS_EVERYTHING_H +#define OPTIMIZATION_CMA_ES_LOGS_EVERYTHING_H -#include "utils/vector.h" #include "utils/matrix.h" -#include "utils/vis_only.h" - -#ifdef VISUALIZE -#include <vector> -#include "json.hpp" +#include "utils/vector.h" #include "utils/eigen_json.h" -#endif +#include <vector> namespace optimization { @@ -22,41 +17,24 @@ struct CmaEsState { }; //-------------------------------------------------------------------------------------------------- -// This is used as a base class rather than a member so that the empty base class optimization can -// be applied (the member would take up space even if it is an empty class). -struct CmaEsLog { - inline void reserve(uint32_t n) VIS_ONLY_METHOD; +struct CmaEsLogEverything { + inline void reserve(uint32_t n) { states.reserve(n); } template <typename Objective> - inline void initialize(Objective const&) VIS_ONLY_METHOD; + inline void initialize(Objective const&) { objective_name = Objective::name; } inline void push_back( double const* const* points, double const* values, uint32_t dim, uint32_t pop_size - ) VIS_ONLY_METHOD; - inline void clear() VIS_ONLY_METHOD; + ); + inline void clear(); - #ifdef VISUALIZE std::string objective_name; std::vector<CmaEsState> states; - #endif }; -#ifdef VISUALIZE - //.................................................................................................. -void CmaEsLog::reserve(uint32_t n) { - states.reserve(n); -} - -//.................................................................................................. -template <typename Objective> -void CmaEsLog::initialize(Objective const&) { - objective_name = Objective::name; -} - -//.................................................................................................. -void CmaEsLog::push_back( +void CmaEsLogEverything::push_back( double const* const* points, double const* values, uint32_t dim, @@ -74,7 +52,8 @@ void CmaEsLog::push_back( } //.................................................................................................. -void CmaEsLog::clear() { +void CmaEsLogEverything::clear() { + objective_name.clear(); states.clear(); } @@ -87,7 +66,7 @@ inline void to_json(nlohmann::json& j, CmaEsState const& state) { } //-------------------------------------------------------------------------------------------------- -inline void to_json(nlohmann::json& j, CmaEsLog const& log) { +inline void to_json(nlohmann::json& j, CmaEsLogEverything const& log) { j = nlohmann::json{ {"algorithm", "gradient descent"}, {"objective", log.objective_name}, @@ -95,8 +74,6 @@ inline void to_json(nlohmann::json& j, CmaEsLog const& log) { }; } -#endif - } #endif diff --git a/optimization/optimizers/cma_es/cma_es_vis.h b/optimization/optimizers/cma_es/logs/everything_vis.h similarity index 88% rename from optimization/optimizers/cma_es/cma_es_vis.h rename to optimization/optimizers/cma_es/logs/everything_vis.h index 7dee87080796c035ce425d04386b46ac496fafbf..84fc55958d5b4f3696f6f45a3b0de0b7fa40e315 100644 --- a/optimization/optimizers/cma_es/cma_es_vis.h +++ b/optimization/optimizers/cma_es/logs/everything_vis.h @@ -1,15 +1,13 @@ -#ifndef OPTIMIZATION_CMA_ES_VIS_H -#define OPTIMIZATION_CMA_ES_VIS_H +#ifndef OPTIMIZATION_CMA_ES_LOGS_EVERYTHING_VIS_H +#define OPTIMIZATION_CMA_ES_LOGS_EVERYTHING_VIS_H -#include <iostream> -#include "cma_es_log.h" -#include "utils/eigen_json.h" +#include "everything.h" namespace optimization { //-------------------------------------------------------------------------------------------------- struct CmaEsVis { - CmaEsLog const& log; + CmaEsLogEverything const& log; }; //.................................................................................................. diff --git a/optimization/optimizers/cma_es/logs/nothing.h b/optimization/optimizers/cma_es/logs/nothing.h new file mode 100644 index 0000000000000000000000000000000000000000..cb6d1600646fd377b85bfb61963bfef7da860e61 --- /dev/null +++ b/optimization/optimizers/cma_es/logs/nothing.h @@ -0,0 +1,19 @@ +#ifndef OPTIMIZATION_CMA_ES_LOGS_NOTHING_H +#define OPTIMIZATION_CMA_ES_LOGS_NOTHING_H + +#include <cstdint> + +namespace optimization { + +//-------------------------------------------------------------------------------------------------- +struct CmaEsLogNothing { + inline void reserve(uint32_t) {} + template <typename Objective> + inline void initialize(Objective const&) {} + inline void push_back(double const* const*, double const*, uint32_t, uint32_t) {} + inline void clear() {} +}; + +} + +#endif diff --git a/optimization/optimizers/cma_es/main.cpp b/optimization/optimizers/cma_es/main.cpp index 8148c9f4700b536a0e49a9bce23367458fc4bcaf..eb159394daf3b5f522a7db99e22409d01cdf7cd7 100644 --- a/optimization/optimizers/cma_es/main.cpp +++ b/optimization/optimizers/cma_es/main.cpp @@ -2,17 +2,12 @@ #include "cma_es.h" #include "objectives/paraboloid.h" #include "objectives/rosenbrock.h" -#include <iostream> - -#ifdef VISUALIZE -#include "cma_es_vis.h" -#include "utils/eigen_json.h" +#include "logs/everything_vis.h" #include <fstream> - -using json = nlohmann::json; -#endif +#include <iostream> using namespace optimization; +using json = nlohmann::json; //-------------------------------------------------------------------------------------------------- int main(int const argc, char const** argv) { @@ -65,24 +60,29 @@ int main(int const argc, char const** argv) { 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); + CmaEsLogEverything log; + + if (log_file_path.empty() && vis_file_path.empty()) { + optimizer.optimize(objective, initial_point, initial_std_dev); + } else { + optimizer.optimize(objective, initial_point, initial_std_dev, log); + } + std::cout << "n evaluations: " << optimizer.n_evaluations() << '\n'; std::cout << "n generations: " << optimizer.n_iterations() << '\n'; - std::cout << "final point: " << minimum << '\n'; + std::cout << "final point: " << optimizer.point() << '\n'; - #ifdef VISUALIZE if (!log_file_path.empty()) { - json data = optimizer; + json data = log; std::ofstream log_file(log_file_path); log_file << data.dump(4) << '\n'; } if (!vis_file_path.empty()) { - json data = CmaEsVis{optimizer}; + json data = CmaEsVis{log}; std::ofstream vis_file(vis_file_path); vis_file << data.dump(4) << '\n'; } - #endif return 0; }