diff --git a/optimization/optimizers/nelder_mead/CMakeLists.txt b/optimization/optimizers/nelder_mead/CMakeLists.txt
index dd74919f1b861a1d9e441ba618927db7911a25b8..5f2bcb508ab80d2209812321c876ad5c2f320dc1 100644
--- a/optimization/optimizers/nelder_mead/CMakeLists.txt
+++ b/optimization/optimizers/nelder_mead/CMakeLists.txt
@@ -3,7 +3,5 @@ add_executable(nelder_mead
)
target_link_libraries(nelder_mead optimization_lib clara)
-if (VISUALIZE)
- make_plot_target(nelder_mead 2d ARGS -d 2)
- make_plot_target(nelder_mead 10d ARGS -d 10 -n 10000)
-endif()
+make_plot_target(nelder_mead 2d ARGS -d 2)
+make_plot_target(nelder_mead 10d ARGS -d 10 -n 10000)
diff --git a/optimization/optimizers/nelder_mead/logs/everything.h b/optimization/optimizers/nelder_mead/logs/everything.h
new file mode 100644
index 0000000000000000000000000000000000000000..f8bf8f03f4885bdc80eafcc279b5c48f0ae13f18
--- /dev/null
+++ b/optimization/optimizers/nelder_mead/logs/everything.h
@@ -0,0 +1,66 @@
+#ifndef OPTIMIZATION_NELDER_MEAD_LOGS_EVERYTHING_H
+#define OPTIMIZATION_NELDER_MEAD_LOGS_EVERYTHING_H
+
+#include "json.hpp"
+#include "../nelder_mead.h"
+#include "utils/eigen_json.h"
+#include <string>
+#include <vector>
+
+namespace optimization {
+
+//--------------------------------------------------------------------------------------------------
+template <int32_t D>
+struct NelderMeadState {
+ typename NelderMead<D>::MatrixDN simplex;
+ typename NelderMead<D>::VectorN values;
+};
+
+//--------------------------------------------------------------------------------------------------
+template <int32_t D>
+struct NelderMeadLogEverything {
+ template <typename Objective>
+ void initialize(Objective const&) { objective_name = Objective::name; }
+
+ void push_back(
+ typename NelderMead<D>::MatrixDN const& simplex,
+ typename NelderMead<D>::VectorN const& values
+ ) {
+ states.push_back({simplex, values});
+ }
+
+ void clear();
+
+ std::string objective_name;
+ std::vector<NelderMeadState<D>> states;
+};
+
+//..................................................................................................
+template <int32_t D>
+void NelderMeadLogEverything<D>::clear() {
+ objective_name.clear();
+ states.clear();
+}
+
+//--------------------------------------------------------------------------------------------------
+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}
+ };
+}
+
+//--------------------------------------------------------------------------------------------------
+template <int32_t D>
+void to_json(nlohmann::json& j, NelderMeadLogEverything<D> const& log) {
+ j = nlohmann::json{
+ {"algorithm", "gradient descent"},
+ {"objective", log.objective_name},
+ {"data", log.states}
+ };
+}
+
+}
+
+#endif
diff --git a/optimization/optimizers/nelder_mead/nelder_mead_vis.h b/optimization/optimizers/nelder_mead/logs/everything_vis.h
similarity index 87%
rename from optimization/optimizers/nelder_mead/nelder_mead_vis.h
rename to optimization/optimizers/nelder_mead/logs/everything_vis.h
index 648b1bca4333edf6c65ca77ebecc59d041d16af3..d7f07c4a02bf59571e9c570ad5b49223076006a2 100644
--- a/optimization/optimizers/nelder_mead/nelder_mead_vis.h
+++ b/optimization/optimizers/nelder_mead/logs/everything_vis.h
@@ -1,16 +1,14 @@
#ifndef OPTIMIZATION_NELDER_MEAD_VIS_H
#define OPTIMIZATION_NELDER_MEAD_VIS_H
-#include <iostream>
-#include "nelder_mead_log.h"
-#include "utils/eigen_json.h"
+#include "everything.h"
namespace optimization {
//--------------------------------------------------------------------------------------------------
template <int32_t D>
struct NelderMeadVis {
- NelderMeadLog<D> const& log;
+ NelderMeadLogEverything<D> const& log;
};
//..................................................................................................
diff --git a/optimization/optimizers/nelder_mead/logs/nothing.h b/optimization/optimizers/nelder_mead/logs/nothing.h
new file mode 100644
index 0000000000000000000000000000000000000000..bb7407076866d30ee62e35cea41c4de1cfe61d22
--- /dev/null
+++ b/optimization/optimizers/nelder_mead/logs/nothing.h
@@ -0,0 +1,20 @@
+#ifndef OPTIMIZATION_NELDER_MEAD_LOGS_NOTHING_H
+#define OPTIMIZATION_NELDER_MEAD_LOGS_NOTHING_H
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+
+namespace optimization {
+
+//--------------------------------------------------------------------------------------------------
+struct NelderMeadLogNothing {
+ template <typename Objective>
+ void initialize(Objective const&) {}
+
+ template <int32_t D, int32_t N>
+ void push_back(MatrixNs<D, N> const&, VectorNs<N> const&) {}
+};
+
+}
+
+#endif
diff --git a/optimization/optimizers/nelder_mead/main.cpp b/optimization/optimizers/nelder_mead/main.cpp
index b70e9897423b1ef498e35c9c5a7342fdad1b79d2..8d73373a5b9b55d8f06a0a81388deaa3cc0aee5e 100644
--- a/optimization/optimizers/nelder_mead/main.cpp
+++ b/optimization/optimizers/nelder_mead/main.cpp
@@ -1,18 +1,13 @@
#include "clara.hpp"
+#include "logs/everything_vis.h"
#include "nelder_mead.h"
#include "objectives/paraboloid.h"
#include "objectives/rosenbrock.h"
#include <iostream>
-
-#ifdef VISUALIZE
-#include "nelder_mead_vis.h"
-#include "utils/eigen_json.h"
#include <fstream>
-using json = nlohmann::json;
-#endif
-
using namespace optimization;
+using json = nlohmann::json;
//--------------------------------------------------------------------------------------------------
int main(int const argc, char const** argv) {
@@ -49,23 +44,29 @@ int main(int const argc, char const** argv) {
objective.dim() = dim;
NelderMead<-1> optimizer(max_evaluations, relative_y_tolerance);
- VectorXs minimum = optimizer.optimize(objective, simplex);
+ NelderMeadLogEverything<-1> log;
+
+ // Only log stuff if we're going to use it.
+ if (log_file_path.empty() && vis_file_path.empty()) {
+ optimizer.optimize(objective, simplex);
+ } else {
+ optimizer.optimize(objective, simplex, log);
+ }
+
std::cout << "n evaluations: " << optimizer.n_evaluations() << '\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 = NelderMeadVis<-1>{optimizer};
+ json data = NelderMeadVis<-1>{log};
std::ofstream vis_file(vis_file_path);
vis_file << data.dump(4) << '\n';
}
- #endif
return 0;
}
diff --git a/optimization/optimizers/nelder_mead/nelder_mead.h b/optimization/optimizers/nelder_mead/nelder_mead.h
index 058f34c37a396bc8d8b21e4b0c41b3fafbd3b04e..e03399971444c703643dbc6bcc44cc585dad924b 100644
--- a/optimization/optimizers/nelder_mead/nelder_mead.h
+++ b/optimization/optimizers/nelder_mead/nelder_mead.h
@@ -1,7 +1,7 @@
#ifndef OPTIMIZATION_NELDER_MEAD_H
#define OPTIMIZATION_NELDER_MEAD_H
-#include "nelder_mead_log.h"
+#include "logs/nothing.h"
#include <cmath>
#include <iostream>
@@ -10,54 +10,90 @@ namespace optimization {
//--------------------------------------------------------------------------------------------------
// TODO record termination condition
template <int32_t D = -1>
-class NelderMead : public NelderMeadLog<D> {
+class NelderMead {
public:
+ // If D != -1, N == D + 1. If D == -1, N == -1.
+ static constexpr int32_t N = [D_ = D](){
+ if (D_ < 0) {
+ return D_;
+ } else {
+ return D_ + 1;
+ }
+ }();
+
// vector in the space we're searching
- using VectorD = typename NelderMeadLog<D>::VectorD;
+ using VectorD = VectorNs<D>;
// vector in one higher dimension (e.g. to store a real value for each simplex vertex)
- using VectorN = typename NelderMeadLog<D>::VectorN;
+ using VectorN = VectorNs<N>;
// D x N matrix (each column is a VectorD holding a vertex, each row is a VectorN)
- using MatrixDN = typename NelderMeadLog<D>::MatrixDN;
+ // There are two reasons to store vertices as columns. First, it avoids the use of
+ // transposition. Eigen will happily tranpose for you if you assign a row vector to a column
+ // vector, but it will complain if you do something like `c_vec = c_vec2 + r_vec`. Second, Eigen
+ // defaults to column major storage, so this keeps our vertex values contiguous. For logging we
+ // transpose the matrix, so each row is a vertex, since this is more natural to deal with in
+ // numpy/matplotlib.
+ using MatrixDN = MatrixNs<D, N>;
+ // Constructor
NelderMead(uint32_t me = 1000, Scalar ry = 1e-8)
: max_evaluations_(me), relative_y_tolerance_(ry)
{}
uint32_t n_evaluations() const { return n_evaluations_; }
-
MatrixDN const& simplex_vertices() const { return simplex_vertices_; }
VectorN const& simplex_values() const { return simplex_values_; }
+ //VectorD const& min_point() const { return simplex_vertices_[i_lowest_]; }
+ decltype(auto) point() const {
+ return simplex_vertices_.col(i_lowest_);
+ }
+ Scalar value() const { return simplex_values_[i_lowest_]; }
- // Constructs a simplex from an initial point by offset all coordinate values.
+ // Constructs a simplex from an initial point by offsetting all coordinate values.
template <typename Objective>
- VectorD optimize(Objective& objective, VectorD const& initial_point, Scalar offset = 1);
- // Each row of simplex is a vertex.
+ VectorD optimize(
+ Objective& objective,
+ VectorD const& initial_point,
+ Scalar offset = 1
+ );
+ template <typename Objective, typename Log>
+ VectorD optimize(
+ Objective& objective,
+ VectorD const& initial_point,
+ Scalar offset,
+ Log& log
+ );
+
+ // Each column 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_]; }
+ template <typename Objective, typename Log>
+ VectorD optimize(Objective& objective, MatrixDN const& simplex, Log& log);
private:
template <typename Objective>
Scalar try_new_point(Objective& objective, Scalar factor);
+ // hyperparameters
+ uint32_t max_evaluations_;
+ //Scalar relative_x_tolerance_;
+ Scalar relative_y_tolerance_;
+
+ // algorithm state
// For fixed size D (and thus N), dim_ and n_vertices_ are redundant.
- uint32_t dim_;
- uint32_t n_vertices_; // == dim_ + 1
- MatrixDN simplex_vertices_;
- VectorN simplex_values_;
- VectorD vertex_sum_;
- // Invariant: simplex_values_[i_lowest] <= simplex_values_[i_second_highest] <= simplex_values_[i_highest]
+ uint32_t dim_; // we are searching in R^dim_
+ uint32_t n_vertices_; // the simplex has n_vertices_ == dim_ + 1 different vertices
+ uint32_t n_evaluations_; // number of times we've evaluated the objective
+ MatrixDN simplex_vertices_; // the simplex itself (each column is a vertex)
+ VectorN simplex_values_; // the function values at the vertices
+ VectorD vertex_sum_; // the sum of all the vertices
+ VectorD x_new_; // used for trial points
+ // Invariant:
+ // simplex_values_[i_lowest] <= simplex_values_[i_second_highest] <= simplex_values_[i_highest]
uint32_t i_lowest_;
uint32_t i_second_highest_;
uint32_t i_highest_;
- uint32_t n_evaluations_;
- uint32_t max_evaluations_;
- //Scalar relative_x_tolerance_;
- Scalar relative_y_tolerance_;
-
+ // algorithm constants
static constexpr Scalar reflection_coefficient_ = -1.0;
static constexpr Scalar expansion_coefficient_ = 2.0;
static constexpr Scalar contraction_coefficient_ = 0.5;
@@ -72,6 +108,19 @@ auto NelderMead<D>::optimize(
Objective& objective,
VectorD const& initial_point,
Scalar offset
+) -> VectorD {
+ NelderMeadLogNothing log;
+ return optimize(objective, initial_point, offset, log);
+}
+
+//..................................................................................................
+template <int32_t D>
+template <typename Objective, typename Log>
+auto NelderMead<D>::optimize(
+ Objective& objective,
+ VectorD const& initial_point,
+ Scalar offset,
+ Log& log
) -> VectorD {
MatrixDN simplex;
simplex.resize(initial_point.size(), initial_point.size() + 1);
@@ -80,20 +129,29 @@ auto NelderMead<D>::optimize(
simplex.col(i + 1) = initial_point;
simplex(i, i + 1) += offset;
}
- return optimize(objective, simplex);
+ return optimize(objective, simplex, log);
}
//..................................................................................................
template <int32_t D>
template <typename Objective>
auto NelderMead<D>::optimize(Objective& objective, MatrixDN const& simplex) -> VectorD {
+ NelderMeadLogNothing log;
+ return optimize(objective, simplex, log);
+}
+
+//..................................................................................................
+template <int32_t D>
+template <typename Objective, typename Log>
+auto NelderMead<D>::optimize(Objective& objective, MatrixDN const& simplex, Log& log) -> VectorD {
+ // initialize state
simplex_vertices_ = simplex;
n_vertices_ = simplex_vertices_.cols();
dim_ = simplex_vertices_.rows();
assert(n_vertices_ == dim_ + 1);
simplex_values_.resize(n_vertices_);
- vertex_sum_.resize(dim_);
vertex_sum_ = simplex_vertices_.rowwise().sum();
+ log.initialize(objective);
// Evaluate the objective at all simplex vertices.
for (uint32_t i = 0u; i < n_vertices_; ++i) {
@@ -101,8 +159,6 @@ auto NelderMead<D>::optimize(Objective& objective, MatrixDN const& simplex) -> V
}
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;
@@ -123,7 +179,7 @@ auto NelderMead<D>::optimize(Objective& objective, MatrixDN const& simplex) -> V
}
// Log the simplex.
- NelderMeadLog<D>::push_back(simplex_vertices_, simplex_values_);
+ log.push_back(simplex_vertices_, simplex_values_);
// Evaluate termination conditions.
Scalar const difference = std::abs(simplex_values_[i_highest_] - simplex_values_[i_lowest_]);
diff --git a/optimization/optimizers/nelder_mead/nelder_mead_log.h b/optimization/optimizers/nelder_mead/nelder_mead_log.h
deleted file mode 100644
index 4229ad2392ad406ac1fee10376c54f6971dee8fe..0000000000000000000000000000000000000000
--- a/optimization/optimizers/nelder_mead/nelder_mead_log.h
+++ /dev/null
@@ -1,117 +0,0 @@
-#ifndef OPTIMIZATION_NELDER_MEAD_LOG_H
-#define OPTIMIZATION_NELDER_MEAD_LOG_H
-
-#include "utils/vector.h"
-#include "utils/matrix.h"
-#include "utils/vis_only.h"
-
-#ifdef VISUALIZE
-#include <vector>
-#include "json.hpp"
-#include "utils/eigen_json.h"
-#endif
-
-namespace optimization {
-
-//--------------------------------------------------------------------------------------------------
-template <int32_t D>
-struct NelderMeadState {
- // If D != -1, N == D + 1. If D == -1, N == -1.
- static constexpr int32_t N = [D_ = D](){
- if (D_ < 0) {
- return D_;
- } else {
- return D_ + 1;
- }
- }();
-
- // vector in the space we're searching
- using VectorD = VectorNs<D>;
- // vector in one higher dimension (e.g. to store a real value for each simplex vertex)
- using VectorN = VectorNs<N>;
- // D x N matrix (each column is a VectorD holding a vertex, each row is a VectorN)
- // There are two reasons to store vertices as columns. First, it avoids the use of
- // transposition. Eigen will happily tranpose for you if you assign a row vector to a column
- // vector, but it will complain if you do something like `c_vec = c_vec2 + r_vec`. Second, Eigen
- // defaults to column major storage, so this keeps our vertex values contiguous. For logging we
- // transpose the matrix, so each row is a vertex, since this is more natural to deal with in
- // numpy/matplotlib.
- using MatrixDN = MatrixNs<D, N>;
-
- NelderMeadState(MatrixDN const& s, VectorN const& v)
- : simplex(s), values(v) {}
-
- MatrixDN simplex;
- VectorN values;
-};
-
-//--------------------------------------------------------------------------------------------------
-// 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 <int32_t D>
-struct NelderMeadLog {
- 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::string objective_name;
- std::vector<NelderMeadState<D>> states;
- #endif
-};
-
-#ifdef VISUALIZE
-
-//..................................................................................................
-template <int32_t D>
-void NelderMeadLog<D>::reserve(uint32_t n) {
- states.reserve(n);
-}
-
-//..................................................................................................
-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
-) {
- states.emplace_back(simplex, values);
-}
-
-
-//--------------------------------------------------------------------------------------------------
-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}
- };
-}
-
-//--------------------------------------------------------------------------------------------------
-template <int32_t D>
-void to_json(nlohmann::json& j, NelderMeadLog<D> const& log) {
- j = nlohmann::json{
- {"algorithm", "nelder-mead"},
- {"objective", log.objective_name},
- {"data", log.states}
- };
-}
-
-#endif
-
-}
-
-#endif