Fix race conditions in differential evolution (#1063)

Through a combination of silly mistakes, I missed a pile of race conditions in the OpenMP threading.

Switch to C++ threading. Note that this change requires serial generation of trial vectors.

Hopefully I can figure out to parallelize the generation of trial vectors to reduce the serial section a la Ahmdahl's law,
while simultaneously keeping thread sanitizer happy.

Author: NAThompson

.gitignore

@@ -27,3 +27,4 @@ cmake-build-debug/*
 .cmake/*
 build.ninja
 .ninja*
+a.out

README.md

@@ -38,11 +38,13 @@ All the implementations are fully generic and support the use of arbitrary "real
 
 These functions also provide the basis of support for the TR1 special functions.
 
-### Root Finding and Function Minimisation
+### Root Finding
 
 A comprehensive set of root-finding algorithms over the real line, both with derivatives and derivative free.
 
-Also function minimisation via Brent's Method.
+### Optimization
+
+Minimization of cost functions via Brent's method and differential evolution.
 
 ### Polynomials and Rational Functions
 

doc/math.qbk

@@ -724,6 +724,9 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [mathpart root_finding Root Finding \& Minimization Algorithms]
 [include roots/roots_overview.qbk]
 [endmathpart] [/mathpart roots Root Finding Algorithms]
+[mathpart optimization Optimization]
+[include optimization/differential_evolution.qbk]
+[endmathpart] [/mathpart optimization Optimization]
 
 [mathpart poly Polynomials and Rational Functions]
 [include internals/polynomial.qbk]

doc/roots/differential_evolution.qbk renamed to doc/optimization/differential_evolution.qbk

@@ -10,16 +10,16 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 [heading Synopsis]
 
 ``
-    #include <boost/math/tools/differential_evolution.hpp>
+    #include <boost/math/optimization/differential_evolution.hpp>
 
-    namespace boost::math::tools {
+    namespace boost::math::optimization {
 
     template <typename ArgumentContainer> struct differential_evolution_parameters {
         using Real = typename ArgumentContainer::value_type;
         ArgumentContainer lower_bounds;
         ArgumentContainer upper_bounds;
-        Real F = static_cast<Real>(0.65);
-        double crossover_ratio = 0.5;
+        Real mutation_factor = static_cast<Real>(0.65);
+        double crossover_probability = 0.5;
         // Population in each generation:
         size_t NP = 200;
         size_t max_generations = 1000;
@@ -47,8 +47,8 @@ We justify this design choice by reference to the "No free lunch" theorem of Wol
 
     `lower_bounds`: A container representing the lower bounds of the optimization space along each dimension. The `.size()` of the bounds should return the dimension of the problem.
     `upper_bounds`: A container representing the upper bounds of the optimization space along each dimension. It should have the same size of `lower_bounds`, and each element should be >= the corresponding element of `lower_bounds`.
-    `F`: The scale factor controlling the rate at which the population evolves (default is 0.65).
-    `crossover_ratio`: The crossover ratio determining the trade-off between exploration and exploitation (default is 0.5).
+    `mutation_factor`: The F or scale factor controlling the rate at which the population evolves (default is 0.65).
+    `crossover_probability`: The crossover probability determining the trade-off between exploration and exploitation (default is 0.5).
     `NP`: The population size (default is 200). Parallelization occurs over the population, so this should be "large".
     `max_generations`: The maximum number of generations for the optimization (default is 1000).
     `threads`: The number of threads to use for parallelization (default is the hardware concurrency). If the objective function is already multithreaded, then this should be set to 1 to prevent oversubscription.
@@ -64,7 +64,7 @@ The most robust way of decreasing the probability of getting stuck in a local mi
 template <typename ArgumentContainer, class Func, class URBG>
 ArgumentContainer differential_evolution(const Func cost_function,
     differential_evolution_parameters<ArgumentContainer> const &de_params,
-    URBG &g,
+    URBG &gen,
     std::invoke_result_t<Func, ArgumentContainer> value_to_reach = std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),
     std::atomic<bool> *cancellation = nullptr,
     std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr,
@@ -102,6 +102,12 @@ Supported termination criteria are explicit requests for termination, value-to-r
 Price, Storn, and Lampinen, Section 2.8 also list population statistics and lack of accepted trials over many generations as sensible termination criteria.
 These could be supported if there is demand.
 
+Parallelization with `std::thread` does have overhead-especially for very fast function calls.
+We found that the function call needs to be roughly a microsecond for unambigous parallelization wins.
+However, we have not provided conditional parallelization as computationally inexpensive cost functions are the exception; not the rule.
+If there is a clear use case for conditional parallelization (cheap cost function in very high dimensions is a good example),
+we can provide it.
+
 [h4:references References]
 
 * Price, Kenneth, Rainer M. Storn, and Jouni A. Lampinen. ['Differential evolution: a practical approach to global optimization.] Springer Science & Business Media, 2006.

doc/roots/roots_overview.qbk

@@ -21,7 +21,6 @@ There are several fully-worked __root_finding_examples, including:
 [include quartic_roots.qbk]
 [include root_finding_examples.qbk]
 [include minima.qbk]
-[include differential_evolution.qbk]
 [include root_comparison.qbk]
 
 [/ roots_overview.qbk

example/differential_evolution.cpp

@@ -5,10 +5,10 @@
  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  */
 #include <iostream>
-#include <boost/math/tools/differential_evolution.hpp>
+#include <boost/math/optimization/differential_evolution.hpp>
 
-using boost::math::tools::differential_evolution_parameters;
-using boost::math::tools::differential_evolution;
+using boost::math::optimization::differential_evolution_parameters;
+using boost::math::optimization::differential_evolution;
 
 double rosenbrock(std::vector<double> const & x) {
    double result = 0;

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

@@ -0,0 +1,165 @@
+/*
+ * Copyright Nick Thompson, 2024
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#ifndef BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP
+#define BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+#include <random>
+
+namespace boost::math::optimization::detail {
+
+template <typename T, typename = void> struct has_resize : std::false_type {};
+
+template <typename T>
+struct has_resize<T, std::void_t<decltype(std::declval<T>().resize(size_t{}))>> : std::true_type {};
+
+template <typename T> constexpr bool has_resize_v = has_resize<T>::value;
+
+template <typename ArgumentContainer>
+void validate_bounds(ArgumentContainer const &lower_bounds, ArgumentContainer const &upper_bounds) {
+  using std::isfinite;
+  std::ostringstream oss;
+  if (lower_bounds.size() == 0) {
+    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+    oss << ": The dimension of the problem cannot be zero.";
+    throw std::domain_error(oss.str());
+  }
+  if (upper_bounds.size() != lower_bounds.size()) {
+    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+    oss << ": There must be the same number of lower bounds as upper bounds, but given ";
+    oss << upper_bounds.size() << " upper bounds, and " << lower_bounds.size() << " lower bounds.";
+    throw std::domain_error(oss.str());
+  }
+  for (size_t i = 0; i < lower_bounds.size(); ++i) {
+    auto lb = lower_bounds[i];
+    auto ub = upper_bounds[i];
+    if (lb > ub) {
+      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+      oss << ": The upper bound must be greater than or equal to the lower bound, but the upper bound is " << ub
+          << " and the lower is " << lb << ".";
+      throw std::domain_error(oss.str());
+    }
+    if (!isfinite(lb)) {
+      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+      oss << ": The lower bound must be finite, but got " << lb << ".";
+      oss << " For infinite bounds, emulate with std::numeric_limits<Real>::lower() or use a standard infinite->finite "
+             "transform.";
+      throw std::domain_error(oss.str());
+    }
+    if (!isfinite(ub)) {
+      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+      oss << ": The upper bound must be finite, but got " << ub << ".";
+      oss << " For infinite bounds, emulate with std::numeric_limits<Real>::max() or use a standard infinite->finite "
+             "transform.";
+      throw std::domain_error(oss.str());
+    }
+  }
+}
+
+template <typename ArgumentContainer, class URBG>
+std::vector<ArgumentContainer> random_initial_population(ArgumentContainer const &lower_bounds,
+                                                         ArgumentContainer const &upper_bounds,
+                                                         size_t initial_population_size, URBG &&gen) {
+  using Real = typename ArgumentContainer::value_type;
+  constexpr bool has_resize = detail::has_resize_v<ArgumentContainer>;
+  std::vector<ArgumentContainer> population(initial_population_size);
+  auto const dimension = lower_bounds.size();
+  for (size_t i = 0; i < population.size(); ++i) {
+    if constexpr (has_resize) {
+      population[i].resize(dimension);
+    } else {
+      // Argument type must be known at compile-time; like std::array:
+      if (population[i].size() != dimension) {
+        std::ostringstream oss;
+        oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+        oss << ": For containers which do not have resize, the default size must be the same as the dimension, ";
+        oss << "but the default container size is " << population[i].size() << " and the dimension of the problem is "
+            << dimension << ".";
+        oss << " The function argument container type is " << typeid(ArgumentContainer).name() << ".\n";
+        throw std::runtime_error(oss.str());
+      }
+    }
+  }
+
+  // Why don't we provide an option to initialize with (say) a Gaussian distribution?
+  // > If the optimum's location is fairly well known,
+  // > a Gaussian distribution may prove somewhat faster, although it
+  // > may also increase the probability that the population will converge prematurely.
+  // > In general, uniform distributions are preferred, since they best reflect
+  // > the lack of knowledge about the optimum's location.
+  //  - Differential Evolution: A Practical Approach to Global Optimization
+  // 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));
+  for (size_t i = 0; i < population.size(); ++i) {
+    for (size_t j = 0; j < dimension; ++j) {
+      auto const &lb = lower_bounds[j];
+      auto const &ub = upper_bounds[j];
+      population[i][j] = lb + dis(gen) * (ub - lb);
+    }
+  }
+
+  return population;
+}
+
+template <typename ArgumentContainer>
+void validate_initial_guess(ArgumentContainer const &initial_guess, ArgumentContainer const &lower_bounds,
+                            ArgumentContainer const &upper_bounds) {
+  std::ostringstream oss;
+  auto const dimension = lower_bounds.size();
+  if (initial_guess.size() != dimension) {
+    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+    oss << ": The initial guess must have the same dimensions as the problem,";
+    oss << ", but the problem size is " << dimension << " and the initial guess has " << initial_guess.size()
+        << " elements.";
+    throw std::domain_error(oss.str());
+  }
+  for (size_t i = 0; i < dimension; ++i) {
+    auto lb = lower_bounds[i];
+    auto ub = upper_bounds[i];
+    if (!isfinite(initial_guess[i])) {
+      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+      oss << ": At index " << i << ", the initial guess is " << initial_guess[i]
+          << ", make sure all elements of the initial guess are finite.";
+      throw std::domain_error(oss.str());
+    }
+    if (initial_guess[i] < lb || initial_guess[i] > ub) {
+      oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+      oss << ": At index " << i << " the initial guess " << initial_guess[i] << " is not in the bounds [" << lb << ", "
+          << ub << "].";
+      throw std::domain_error(oss.str());
+    }
+  }
+}
+
+// Return indices corresponding to the minimum function values.
+template <typename Real> std::vector<size_t> best_indices(std::vector<Real> const &function_values) {
+  using std::isnan;
+  const size_t n = function_values.size();
+  std::vector<size_t> indices(n);
+  for (size_t i = 0; i < n; ++i) {
+    indices[i] = i;
+  }
+
+  std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) {
+    if (isnan(function_values[a])) {
+      return false;
+    }
+    if (isnan(function_values[b])) {
+      return true;
+    }
+    return function_values[a] < function_values[b];
+  });
+  return indices;
+}
+
+} // namespace boost::math::optimization::detail
+#endif