From 1ec3187856c69622b7a3c8050e608355aa177cab Mon Sep 17 00:00:00 2001
From: Erik Strand <erik.strand@cba.mit.edu>
Date: Thu, 16 Apr 2020 17:09:09 -0400
Subject: [PATCH] Don't template Nelder-Mead on the objective

---
 apps/compare_convergence.cpp                  |  2 +-
 optimization/optimizers/nelder_mead/main.cpp  |  4 +-
 .../optimizers/nelder_mead/nelder_mead.h      | 32 +++++++++------
 .../optimizers/nelder_mead/nelder_mead_log.h  | 40 ++++++++++++-------
 .../optimizers/nelder_mead/nelder_mead_vis.h  | 10 ++---
 5 files changed, 53 insertions(+), 35 deletions(-)

diff --git a/apps/compare_convergence.cpp b/apps/compare_convergence.cpp
index ef25945..b154ac4 100644
--- a/apps/compare_convergence.cpp
+++ b/apps/compare_convergence.cpp
@@ -25,7 +25,7 @@ int main() {
     Scalar value_threshold = 1e-8;
     std::vector<uint32_t> dims = {2, 4, 8, 16, 32, 64, 128, 256, 512, 1024};
 
-    NelderMead<Objective, -1> nm(max_evaluations, value_threshold);
+    NelderMead<-1> nm(max_evaluations, value_threshold);
 
     // We'll set the learning rate in the loop.
     GradientDescent<-1> gd(0, max_evaluations, gradient_threshold);
diff --git a/optimization/optimizers/nelder_mead/main.cpp b/optimization/optimizers/nelder_mead/main.cpp
index b2819e2..b70e989 100644
--- a/optimization/optimizers/nelder_mead/main.cpp
+++ b/optimization/optimizers/nelder_mead/main.cpp
@@ -48,7 +48,7 @@ int main(int const argc, char const** argv) {
     Objective objective;
     objective.dim() = dim;
 
-    NelderMead<Objective, -1> optimizer(max_evaluations, relative_y_tolerance);
+    NelderMead<-1> optimizer(max_evaluations, relative_y_tolerance);
     VectorXs minimum = optimizer.optimize(objective, simplex);
     std::cout << "n evaluations: " << optimizer.n_evaluations() << '\n';
     std::cout << "final point: " << minimum << '\n';
@@ -61,7 +61,7 @@ int main(int const argc, char const** argv) {
     }
 
     if (!vis_file_path.empty()) {
-        json data = NelderMeadVis<Objective, -1>{optimizer};
+        json data = NelderMeadVis<-1>{optimizer};
         std::ofstream vis_file(vis_file_path);
         vis_file << data.dump(4) << '\n';
     }
diff --git a/optimization/optimizers/nelder_mead/nelder_mead.h b/optimization/optimizers/nelder_mead/nelder_mead.h
index 36826b1..058f34c 100644
--- a/optimization/optimizers/nelder_mead/nelder_mead.h
+++ b/optimization/optimizers/nelder_mead/nelder_mead.h
@@ -9,15 +9,15 @@ namespace optimization {
 
 //--------------------------------------------------------------------------------------------------
 // TODO record termination condition
-template <typename Objective, int32_t D = -1>
-class NelderMead : public NelderMeadLog<Objective, D> {
+template <int32_t D = -1>
+class NelderMead : public NelderMeadLog<D> {
 public:
     // vector in the space we're searching
-    using VectorD = typename NelderMeadLog<Objective, D>::VectorD;
+    using VectorD = typename NelderMeadLog<D>::VectorD;
     // vector in one higher dimension (e.g. to store a real value for each simplex vertex)
-    using VectorN = typename NelderMeadLog<Objective, D>::VectorN;
+    using VectorN = typename NelderMeadLog<D>::VectorN;
     // D x N matrix (each column is a VectorD holding a vertex, each row is a VectorN)
-    using MatrixDN = typename NelderMeadLog<Objective, D>::MatrixDN;
+    using MatrixDN = typename NelderMeadLog<D>::MatrixDN;
 
     NelderMead(uint32_t me = 1000, Scalar ry = 1e-8)
         : max_evaluations_(me), relative_y_tolerance_(ry)
@@ -29,14 +29,17 @@ public:
     VectorN const& simplex_values() const { return simplex_values_; }
 
     // Constructs a simplex from an initial point by offset all coordinate values.
+    template <typename Objective>
     VectorD optimize(Objective& objective, VectorD const& initial_point, Scalar offset = 1);
     // Each row of simplex is a vertex.
+    template <typename Objective>
     VectorD optimize(Objective& objective, MatrixDN const& simplex);
 
     VectorD const& min_point() const { return simplex_vertices_[i_lowest_]; }
     Scalar min_value() const { return simplex_values_[i_lowest_]; }
 
 private:
+    template <typename Objective>
     Scalar try_new_point(Objective& objective, Scalar factor);
 
     // For fixed size D (and thus N), dim_ and n_vertices_ are redundant.
@@ -63,8 +66,9 @@ private:
 };
 
 //..................................................................................................
-template <typename Objective, int32_t D>
-auto NelderMead<Objective, D>::optimize(
+template <int32_t D>
+template <typename Objective>
+auto NelderMead<D>::optimize(
     Objective& objective,
     VectorD const& initial_point,
     Scalar offset
@@ -80,8 +84,9 @@ auto NelderMead<Objective, D>::optimize(
 }
 
 //..................................................................................................
-template <typename Objective, int32_t D>
-auto NelderMead<Objective, D>::optimize(Objective& objective, MatrixDN const& simplex) -> VectorD {
+template <int32_t D>
+template <typename Objective>
+auto NelderMead<D>::optimize(Objective& objective, MatrixDN const& simplex) -> VectorD {
     simplex_vertices_ = simplex;
     n_vertices_ = simplex_vertices_.cols();
     dim_ = simplex_vertices_.rows();
@@ -96,6 +101,8 @@ auto NelderMead<Objective, D>::optimize(Objective& objective, MatrixDN const& si
     }
     n_evaluations_ = n_vertices_;
 
+    NelderMeadLog<D>::initialize(objective);
+
     while (true) {
         // Find lowest, highest, and second highest points.
         i_lowest_ = (simplex_values_[0] < simplex_values_[1]) ? 0 : 1;
@@ -116,7 +123,7 @@ auto NelderMead<Objective, D>::optimize(Objective& objective, MatrixDN const& si
         }
 
         // Log the simplex.
-        NelderMeadLog<Objective, D>::push_back(simplex_vertices_, simplex_values_);
+        NelderMeadLog<D>::push_back(simplex_vertices_, simplex_values_);
 
         // Evaluate termination conditions.
         Scalar const difference = std::abs(simplex_values_[i_highest_] - simplex_values_[i_lowest_]);
@@ -170,8 +177,9 @@ auto NelderMead<Objective, D>::optimize(Objective& objective, MatrixDN const& si
 }
 
 //..................................................................................................
-template <typename Objective, int32_t D>
-Scalar NelderMead<Objective, D>::try_new_point(Objective& objective, Scalar factor) {
+template <int32_t D>
+template <typename Objective>
+Scalar NelderMead<D>::try_new_point(Objective& objective, Scalar factor) {
     // Generate a new point by reflecting/expanding/contracting the worst point.
     Scalar const t1 = (Scalar(1) - factor) / dim_;
     Scalar const t2 = factor - t1;
diff --git a/optimization/optimizers/nelder_mead/nelder_mead_log.h b/optimization/optimizers/nelder_mead/nelder_mead_log.h
index 0e39c7d..4229ad2 100644
--- a/optimization/optimizers/nelder_mead/nelder_mead_log.h
+++ b/optimization/optimizers/nelder_mead/nelder_mead_log.h
@@ -14,7 +14,7 @@
 namespace optimization {
 
 //--------------------------------------------------------------------------------------------------
-template <typename Objective, int32_t D>
+template <int32_t D>
 struct NelderMeadState {
     // If D != -1, N == D + 1. If D == -1, N == -1.
     static constexpr int32_t N = [D_ = D](){
@@ -48,32 +48,42 @@ struct NelderMeadState {
 //--------------------------------------------------------------------------------------------------
 // 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).
-template <typename Objective, int32_t D>
+template <int32_t D>
 struct NelderMeadLog {
-    using VectorD = typename NelderMeadState<Objective, D>::VectorD;
-    using VectorN = typename NelderMeadState<Objective, D>::VectorN;
-    using MatrixDN = typename NelderMeadState<Objective, D>::MatrixDN;
+    using VectorD = typename NelderMeadState<D>::VectorD;
+    using VectorN = typename NelderMeadState<D>::VectorN;
+    using MatrixDN = typename NelderMeadState<D>::MatrixDN;
 
     void reserve(uint32_t n) VIS_ONLY_METHOD;
+    template <typename Objective>
+    void initialize(Objective const&) VIS_ONLY_METHOD;
     void push_back(MatrixDN const& simplex, VectorN const& values) VIS_ONLY_METHOD;
     void clear() VIS_ONLY_METHOD;
 
     #ifdef VISUALIZE
-    std::vector<NelderMeadState<Objective, D>> states;
+    std::string objective_name;
+    std::vector<NelderMeadState<D>> states;
     #endif
 };
 
 #ifdef VISUALIZE
 
 //..................................................................................................
-template <typename Objective, int32_t D>
-void NelderMeadLog<Objective, D>::reserve(uint32_t n) {
+template <int32_t D>
+void NelderMeadLog<D>::reserve(uint32_t n) {
     states.reserve(n);
 }
 
 //..................................................................................................
-template <typename Objective, int32_t D>
-void NelderMeadLog<Objective, D>::push_back(
+template <int32_t D>
+template <typename Objective>
+void NelderMeadLog<D>::initialize(Objective const&) {
+    objective_name = Objective::name;
+}
+
+//..................................................................................................
+template <int32_t D>
+void NelderMeadLog<D>::push_back(
     MatrixDN const& simplex,
     VectorN const& values
 ) {
@@ -82,8 +92,8 @@ void NelderMeadLog<Objective, D>::push_back(
 
 
 //--------------------------------------------------------------------------------------------------
-template <typename Objective, int32_t D>
-void to_json(nlohmann::json& j, NelderMeadState<Objective, D> const& state) {
+template <int32_t D>
+void to_json(nlohmann::json& j, NelderMeadState<D> const& state) {
     j = nlohmann::json{
         {"simplex", transpose_for_json(state.simplex)},
         {"values", state.values}
@@ -91,11 +101,11 @@ void to_json(nlohmann::json& j, NelderMeadState<Objective, D> const& state) {
 }
 
 //--------------------------------------------------------------------------------------------------
-template <typename Objective, int32_t D>
-void to_json(nlohmann::json& j, NelderMeadLog<Objective, D> const& log) {
+template <int32_t D>
+void to_json(nlohmann::json& j, NelderMeadLog<D> const& log) {
     j = nlohmann::json{
         {"algorithm", "nelder-mead"},
-        {"objective", Objective::name},
+        {"objective", log.objective_name},
         {"data", log.states}
     };
 }
diff --git a/optimization/optimizers/nelder_mead/nelder_mead_vis.h b/optimization/optimizers/nelder_mead/nelder_mead_vis.h
index a81d474..648b1bc 100644
--- a/optimization/optimizers/nelder_mead/nelder_mead_vis.h
+++ b/optimization/optimizers/nelder_mead/nelder_mead_vis.h
@@ -8,17 +8,17 @@
 namespace optimization {
 
 //--------------------------------------------------------------------------------------------------
-template <typename Objective, int32_t D>
+template <int32_t D>
 struct NelderMeadVis {
-    NelderMeadLog<Objective, D> const& log;
+    NelderMeadLog<D> const& log;
 };
 
 //..................................................................................................
-template <typename Objective, int32_t D>
-void to_json(nlohmann::json& j, NelderMeadVis<Objective, D> const& vis) {
+template <int32_t D>
+void to_json(nlohmann::json& j, NelderMeadVis<D> const& vis) {
     j = nlohmann::json{
         {"algorithm", "nelder-mead"},
-        {"objective", Objective::name},
+        {"objective", vis.log.objective_name},
         {"data", {
             {"polygons", nlohmann::json::array()}
         }}
-- 
GitLab