[CI SKIP][ci skip] Units compatibility for optimization

include/boost/math/optimization/detail/common.hpp

@@ -68,6 +68,7 @@ std::vector<ArgumentContainer> random_initial_population(ArgumentContainer const
                                                          ArgumentContainer const &upper_bounds,
                                                          size_t initial_population_size, URBG &&gen) {
   using Real = typename ArgumentContainer::value_type;
+  using DimensionlessReal = decltype(Real()/Real());
   constexpr bool has_resize = detail::has_resize_v<ArgumentContainer>;
   std::vector<ArgumentContainer> population(initial_population_size);
   auto const dimension = lower_bounds.size();
@@ -98,7 +99,7 @@ std::vector<ArgumentContainer> random_initial_population(ArgumentContainer const
   // That said, scipy uses Latin Hypercube sampling and says self-avoiding sequences are preferable.
   // So this is something that could be investigated and potentially improved.
   using std::uniform_real_distribution;
-  uniform_real_distribution<Real> dis(Real(0), Real(1));
+  uniform_real_distribution<DimensionlessReal> dis(DimensionlessReal(0), DimensionlessReal(1));
   for (size_t i = 0; i < population.size(); ++i) {
     for (size_t j = 0; j < dimension; ++j) {
       auto const &lb = lower_bounds[j];
@@ -163,6 +164,7 @@ template <typename Real> std::vector<size_t> best_indices(std::vector<Real> cons
 
 template<typename RandomAccessContainer>
 auto weighted_lehmer_mean(RandomAccessContainer const & values, RandomAccessContainer const & weights) {
+  using std::isfinite;
   if (values.size() != weights.size()) {
     std::ostringstream oss;
     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
@@ -179,7 +181,7 @@ auto weighted_lehmer_mean(RandomAccessContainer const & values, RandomAccessCont
   Real numerator = 0;
   Real denominator = 0;
   for (size_t i = 0; i < values.size(); ++i) {
-    if (weights[i] < 0 || !std::isfinite(weights[i])) {
+    if (weights[i] < 0 || !isfinite(weights[i])) {
       std::ostringstream oss;
       oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
       oss << ": All weights must be positive and finite, but got received weight " << weights[i] << " at index " << i << " of " << weights.size() << ".";

include/boost/math/optimization/differential_evolution.hpp

@@ -35,11 +35,12 @@ namespace boost::math::optimization {
 // We provide the parameters in a struct-there are too many of them and they are too unwieldy to pass individually:
 template <typename ArgumentContainer> struct differential_evolution_parameters {
   using Real = typename ArgumentContainer::value_type;
+  using DimensionlessReal = decltype(Real()/Real());
   ArgumentContainer lower_bounds;
   ArgumentContainer upper_bounds;
   // mutation factor is also called scale factor or just F in the literature:
-  Real mutation_factor = static_cast<Real>(0.65);
-  double crossover_probability = 0.5;
+  DimensionlessReal mutation_factor = static_cast<DimensionlessReal>(0.65);
+  DimensionlessReal crossover_probability = static_cast<DimensionlessReal>(0.5);
   // Population in each generation:
   size_t NP = 500;
   size_t max_generations = 1000;
@@ -94,6 +95,7 @@ ArgumentContainer differential_evolution(
     std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr,
     std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr) {
   using Real = typename ArgumentContainer::value_type;
+  using DimensionlessReal = decltype(Real()/Real());
   using ResultType = std::invoke_result_t<Func, ArgumentContainer>;
   using std::clamp;
   using std::isnan;
@@ -159,7 +161,7 @@ ArgumentContainer differential_evolution(
     // I first tried seeding thread-local random number generators with the global generator,
     // but ThreadSanitizer didn't like it. I *suspect* there's a way around this, but
     // even if it's formally threadsafe, there's still a lot of effort to get it computationally reproducible.
-    uniform_real_distribution<Real> unif01(Real(0), Real(1));
+    uniform_real_distribution<DimensionlessReal> unif01(DimensionlessReal(0), DimensionlessReal(1));
     for (size_t i = 0; i < NP; ++i) {
       size_t r1, r2, r3;
       do {

include/boost/math/optimization/jso.hpp

@@ -36,6 +36,7 @@ namespace boost::math::optimization {
 // IEEE Transactions on evolutionary computation, 13(5), 945-958."
 template <typename ArgumentContainer> struct jso_parameters {
   using Real = typename ArgumentContainer::value_type;
+  using DimensionlessReal = decltype(Real()/Real());
   ArgumentContainer lower_bounds;
   ArgumentContainer upper_bounds;
   // Population in the first generation.
@@ -107,6 +108,7 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
                           std::invoke_result_t<Func, ArgumentContainer>>>
         *queries = nullptr) {
   using Real = typename ArgumentContainer::value_type;
+  using DimensionlessReal = decltype(Real()/Real());
   validate_jso_parameters(jso_params);
 
   using ResultType = std::invoke_result_t<Func, ArgumentContainer>;
@@ -116,6 +118,7 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
   using std::isnan;
   using std::max;
   using std::round;
+  using std::isfinite;
   using std::uniform_real_distribution;
 
   // Refer to the referenced paper, pseudocode on page 1313:
@@ -173,12 +176,12 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
   // last bullet, which claims this should be set to 0.3. The reference
   // implementation also does 0.3:
   size_t H = 5;
-  std::vector<Real> M_F(H, static_cast<Real>(0.3));
+  std::vector<DimensionlessReal> M_F(H, static_cast<DimensionlessReal>(0.3));
   // Algorithm 1: jSO algorithm, Line 4:
   // "Set all values in M_CR to 0.8":
-  std::vector<Real> M_CR(H, static_cast<Real>(0.8));
+  std::vector<DimensionlessReal> M_CR(H, static_cast<DimensionlessReal>(0.8));
 
-  std::uniform_real_distribution<Real> unif01(Real(0), Real(1));
+  std::uniform_real_distribution<DimensionlessReal> unif01(DimensionlessReal(0), DimensionlessReal(1));
   bool keep_going = !target_attained;
   if (cancellation && *cancellation) {
     keep_going = false;
@@ -190,8 +193,8 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
   size_t minimum_population_size = (max)(size_t(4), size_t(jso_params.threads));
   while (keep_going) {
     // Algorithm 1, jSO, Line 7:
-    std::vector<Real> S_CR;
-    std::vector<Real> S_F;
+    std::vector<DimensionlessReal> S_CR;
+    std::vector<DimensionlessReal> S_F;
     // Equation 9 of L-SHADE:
     std::vector<ResultType> delta_f;
     for (size_t i = 0; i < population.size(); ++i) {
@@ -203,45 +206,45 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
       // I confess I find it weird to store the historical memory if we're just
       // gonna ignore it, but that's what the paper and the reference
       // implementation says!
-      Real mu_F = static_cast<Real>(0.9);
-      Real mu_CR = static_cast<Real>(0.9);
+      DimensionlessReal mu_F = static_cast<DimensionlessReal>(0.9);
+      DimensionlessReal mu_CR = static_cast<DimensionlessReal>(0.9);
       if (ri != H - 1) {
         mu_F = M_F[ri];
         mu_CR = M_CR[ri];
       }
       // Algorithm 1, jSO, Line 14-18:
-      Real crossover_probability = static_cast<Real>(0);
+      DimensionlessReal crossover_probability = static_cast<DimensionlessReal>(0);
       if (mu_CR >= 0) {
         using std::normal_distribution;
-        normal_distribution<Real> normal(mu_CR, static_cast<Real>(0.1));
+        normal_distribution<DimensionlessReal> normal(mu_CR, static_cast<DimensionlessReal>(0.1));
         crossover_probability = normal(gen);
         // Clamp comes from L-SHADE description:
-        crossover_probability = clamp(crossover_probability, Real(0), Real(1));
+        crossover_probability = clamp(crossover_probability, DimensionlessReal(0), DimensionlessReal(1));
       }
       // Algorithm 1, jSO, Line 19-23:
       // Note that the pseudocode uses a "max_generations parameter",
       // but the reference implementation does not.
       // Since we already require specification of max_function_evaluations,
       // the pseudocode adds an unnecessary parameter.
       if (4 * function_evaluations < jso_params.max_function_evaluations) {
-        crossover_probability = (max)(crossover_probability, Real(0.7));
+        crossover_probability = (max)(crossover_probability, DimensionlessReal(0.7));
       } else if (2 * function_evaluations <
                  jso_params.max_function_evaluations) {
-        crossover_probability = (max)(crossover_probability, Real(0.6));
+        crossover_probability = (max)(crossover_probability, DimensionlessReal(0.6));
       }
 
       // Algorithm 1, jSO, Line 24-27:
       // Note the adjustments to the pseudocode given in the reference
       // implementation.
-      cauchy_distribution<Real> cauchy(mu_F, static_cast<Real>(0.1));
-      Real F;
+      cauchy_distribution<DimensionlessReal> cauchy(mu_F, static_cast<DimensionlessReal>(0.1));
+      DimensionlessReal F;
       do {
         F = cauchy(gen);
         if (F > 1) {
           F = 1;
         }
       } while (F <= 0);
-      Real threshold = static_cast<Real>(7) / static_cast<Real>(10);
+      DimensionlessReal threshold = static_cast<DimensionlessReal>(7) / static_cast<DimensionlessReal>(10);
       if ((10 * function_evaluations <
            6 * jso_params.max_function_evaluations) &&
           (F > threshold)) {
@@ -250,16 +253,16 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
       // > p value for mutation strategy linearly decreases from pmax to pmin
       // during the evolutionary process, > where pmax = 0.25 in jSO and pmin =
       // pmax/2.
-      Real p = Real(0.25) * (1 - static_cast<Real>(function_evaluations) /
+      DimensionlessReal p = DimensionlessReal(0.25) * (1 - static_cast<DimensionlessReal>(function_evaluations) /
                                      (2 * jso_params.max_function_evaluations));
       // Equation (4) of the reference:
-      Real Fw = static_cast<Real>(1.2) * F;
+      DimensionlessReal Fw = static_cast<DimensionlessReal>(1.2) * F;
       if (10 * function_evaluations < 4 * jso_params.max_function_evaluations) {
         if (10 * function_evaluations <
             2 * jso_params.max_function_evaluations) {
-          Fw = static_cast<Real>(0.7) * F;
+          Fw = static_cast<DimensionlessReal>(0.7) * F;
         } else {
-          Fw = static_cast<Real>(0.8) * F;
+          Fw = static_cast<DimensionlessReal>(0.8) * F;
         }
       }
       // Algorithm 1, jSO, Line 28:
@@ -367,13 +370,13 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
     // If there are no successful updates this generation, we do not update the
     // historical memory:
     if (S_CR.size() > 0) {
-      std::vector<Real> weights(S_CR.size(),
-                                std::numeric_limits<Real>::quiet_NaN());
+      std::vector<DimensionlessReal> weights(S_CR.size(),
+                                std::numeric_limits<DimensionlessReal>::quiet_NaN());
       ResultType delta_sum = static_cast<ResultType>(0);
       for (auto const &delta : delta_f) {
         delta_sum += delta;
       }
-      if (delta_sum <= 0 || !std::isfinite(delta_sum)) {
+      if (delta_sum <= 0 || !isfinite(delta_sum)) {
         std::ostringstream oss;
         oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
         oss << "\n\tYou've hit a bug: The sum of improvements must be strictly "

include/boost/math/optimization/random_search.hpp

@@ -60,6 +60,7 @@ ArgumentContainer random_search(
     std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr)
  {
   using Real = typename ArgumentContainer::value_type;
+  using DimensionlessReal = decltype(Real()/Real());
   using ResultType = std::invoke_result_t<Func, ArgumentContainer>;
   using std::isnan;
   using std::uniform_real_distribution;
@@ -106,7 +107,7 @@ ArgumentContainer random_search(
             break;
         }
         // Fill trial vector: 
-        uniform_real_distribution<Real> unif01(Real(0), Real(1));
+        uniform_real_distribution<DimensionlessReal> unif01(DimensionlessReal(0), DimensionlessReal(1));
         for (size_t i = 0; i < dimension; ++i) {
             trial_vector[i] = params.lower_bounds[i] + (params.upper_bounds[i] - params.lower_bounds[i])*unif01(g);
         }

test/differential_evolution_test.cpp

@@ -188,6 +188,19 @@ void test_beale() {
   CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(2e-4));
 }
 
+#if BOOST_MATH_TEST_UNITS_COMPATIBILITY
+void test_dimensioned_sphere() {
+  std::cout << "Testing differential evolution on dimensioned sphere . . .\n";
+  using ArgType = std::vector<quantity<length>>;
+  auto params = differential_evolution_parameters<ArgType>();
+  params.lower_bounds.resize(4, -1.0*meter);
+  params.upper_bounds.resize(4, 1*meter);
+  params.threads = 2;
+  std::mt19937_64 gen(56789);
+  auto local_minima = differential_evolution(dimensioned_sphere, params, gen);
+}
+#endif
+
 int main() {
 
 #if defined(__clang__) || defined(_MSC_VER)
@@ -197,6 +210,9 @@ int main() {
   test_rastrigin<float>();
   test_three_hump_camel<float>();
   test_beale<double>();
+#endif
+#if BOOST_MATH_TEST_UNITS_COMPATIBILITY
+  test_dimensioned_sphere();
 #endif
   test_sphere();
   test_parameter_checks();

test/jso_test.cpp

@@ -143,6 +143,19 @@ void test_beale() {
   CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(2e-4));
 }
 
+#if BOOST_MATH_TEST_UNITS_COMPATIBILITY
+void test_dimensioned_sphere() {
+  std::cout << "Testing jso on dimensioned sphere . . .\n";
+  using ArgType = std::vector<quantity<length>>;
+  auto params = jso_parameters<ArgType>();
+  params.lower_bounds.resize(4, -1.0*meter);
+  params.upper_bounds.resize(4, 1*meter);
+  params.threads = 2;
+  std::mt19937_64 gen(56789);
+  auto local_minima = jso(dimensioned_sphere, params, gen);
+}
+#endif
+
 int main() {
 #if defined(__clang__) || defined(_MSC_VER)
   test_ackley<float>();
@@ -151,6 +164,9 @@ int main() {
   test_rastrigin<double>();
   test_three_hump_camel<float>();
   test_beale<double>();
+#endif
+#if BOOST_MATH_TEST_UNITS_COMPATIBILITY
+  test_dimensioned_sphere();
 #endif
   test_sphere();
   test_weighted_lehmer_mean();

test/random_search_test.cpp

@@ -10,7 +10,6 @@
 #include <boost/math/optimization/random_search.hpp>
 #include <random>
 #include <limits>
-
 using boost::math::optimization::random_search;
 using boost::math::optimization::random_search_parameters;
 
@@ -149,6 +148,21 @@ void test_beale() {
   CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(0.1));
 }
 
+#if BOOST_MATH_TEST_UNITS_COMPATIBILITY
+void test_dimensioned_sphere() {
+  std::cout << "Testing random search on dimensioned sphere . . .\n";
+  using ArgType = std::vector<quantity<length>>;
+  auto rs_params = random_search_parameters<ArgType>();
+  rs_params.lower_bounds.resize(4, -1.0*meter);
+  rs_params.upper_bounds.resize(4, 1*meter);
+  rs_params.max_function_calls = 100000;
+  rs_params.threads = 2;
+  std::mt19937_64 gen(56789);
+  auto local_minima = random_search(dimensioned_sphere, rs_params, gen);
+}
+
+#endif
+
 int main() {
 #if defined(__clang__) || defined(_MSC_VER)
   test_ackley<float>();
@@ -157,6 +171,9 @@ int main() {
   test_rastrigin<double>();
   test_three_hump_camel<float>();
   test_beale<double>();
+#endif
+#if BOOST_MATH_TEST_UNITS_COMPATIBILITY
+  test_dimensioned_sphere();
 #endif
   test_sphere();
   return boost::math::test::report_errors();

test/test_functions_for_optimization.hpp

@@ -6,7 +6,32 @@
  */
 #ifndef TEST_FUNCTIONS_FOR_OPTIMIZATION_HPP
 #define TEST_FUNCTIONS_FOR_OPTIMIZATION_HPP
+#include <array>
+#include <vector>
 #include <boost/math/constants/constants.hpp>
+#if __has_include(<boost/units/systems/si/length.hpp>)
+#define BOOST_MATH_TEST_UNITS_COMPATIBILITY 1
+#include <boost/units/systems/si/length.hpp>
+#include <boost/units/systems/si/area.hpp>
+#include <boost/units/cmath.hpp>
+#include <boost/units/quantity.hpp>
+#include <boost/units/systems/si/io.hpp>
+using namespace boost::units;
+using namespace boost::units::si;
+
+// This *should* return an area, but see: https://github.com/boostorg/units/issues/58
+// This sadly prevents std::atomic<quantity<area>>.
+// Nonetheless, we *do* get some information making the argument type dimensioned,
+// even if it would be better to get the full information:
+double dimensioned_sphere(std::vector<quantity<length>> const & v) {
+  quantity<area> r(0.0*meters*meters);
+  for (auto const & x : v) {
+    r += (x * x);
+  }
+  quantity<area> scale(1.0*meters*meters);
+  return static_cast<double>(r/scale);
+}
+#endif
 
 // Taken from: https://en.wikipedia.org/wiki/Test_functions_for_optimization
 template <typename Real> Real ackley(std::array<Real, 2> const &v) {