Random search (#1071)

doc/math.qbk

@@ -726,6 +726,8 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [endmathpart] [/mathpart roots Root Finding Algorithms]
 [mathpart optimization Optimization]
 [include optimization/differential_evolution.qbk]
+[include optimization/jso.qbk]
+[include optimization/random_search.qbk]
 [endmathpart] [/mathpart optimization Optimization]
 
 [mathpart poly Polynomials and Rational Functions]

doc/optimization/random_search.qbk

@@ -0,0 +1,91 @@
+[/
+Copyright (c) 2024 Nick Thompson
+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)
+]
+
+[section:random_search Random Search]
+
+[heading Synopsis]
+
+``
+    #include <boost/math/optimization/random_search.hpp>
+    
+    namespace boost::math::optimization {
+
+    template <typename ArgumentContainer> struct random_search_parameters {
+        using Real = typename ArgumentContainer::value_type;
+        ArgumentContainer lower_bounds;
+        ArgumentContainer upper_bounds;
+        size_t max_function_calls = 0;
+        ArgumentContainer const * initial_guess = nullptr;
+    };
+
+    template <typename ArgumentContainer, class Func, class URBG>
+    ArgumentContainer random_search(
+        const Func cost_function,
+        random_search_parameters<ArgumentContainer> const &params,
+        URBG &gen,
+        std::invoke_result_t<Func, ArgumentContainer> target_value = std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),
+        std::atomic<bool> *cancellation = nullptr,
+        std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,
+        std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr);
+
+    } // namespaces
+``
+
+The `random_search` function searches for a global minimum of a function.
+There is no special sauce to this algorithm-it merely blasts function calls over threads.
+It's existence is justified by the "No free lunch" theorem in optimization,
+which "establishes that for any algorithm, any elevated performance over one class of problems is offset by performance over another class."
+In practice, it is not clear that the conditions of the NFL theorem holds,
+and on test cases, `random_search` is slower and less accurate than (say) differential evolution and jSO.
+However, it is often the case that rapid convergence is not the goal:
+For example, we often want to spend some time exploring the cost function surface before moving to a faster converging algorithm.
+In addition, random search is embarrassingly parallel, which allows us to avoid Amdahl's law-induced performance problems.
+
+
+[heading Parameters]
+
+    `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`.
+    `max_function_calls`: Defaults to 10000*threads.
+    `initial_guess`: An optional guess for where we should start looking for solutions. This is provided for consistency with other optimization functions-it's not particularly useful for this function.
+
+[heading The function]
+
+``
+template <typename ArgumentContainer, class Func, class URBG>
+ArgumentContainer random_search(const Func cost_function,
+    random_search_parameters<ArgumentContainer> const &params,
+    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::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,
+    std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr)
+``
+
+Parameters:
+
+    `cost_function`: The cost function to be minimized.
+    `params`: The parameters to the algorithm as described above.
+    `gen`: A uniform random bit generator, like `std::mt19937_64`.
+    `value_to_reach`: An optional value that, if reached, stops the optimization. This is the most robust way to terminate the calculation, but in most cases the optimal value of the cost function is unknown. If it is, use it! Physical considerations can often be used to find optimal values for cost functions.
+    `cancellation`: An optional atomic boolean to allow the user to stop the computation and gracefully return the best result found up to that point. N.B.: Cancellation is not immediate; the in-progress generation finishes.
+    `current_minimum_cost`: An optional atomic variable to store the current minimum cost during optimization. This allows developers to (e.g.) plot the progress of the minimization over time and in conjunction with the `cancellation` argument allow halting the computation when the progress stagnates.
+    `queries`: An optional vector to store intermediate results during optimization. This is useful for debugging and perhaps volume rendering of the objective function after the calculation is complete.
+
+Returns:
+
+The argument vector corresponding to the minimum cost found by the optimization.
+
+[h4:examples Examples]
+
+An example exhibiting graceful cancellation and progress observability can be studied in [@../../example/random_search_example.cpp random_search_example.cpp].
+
+[h4:references References]
+
+* D. H. Wolpert and W. G. Macready, ['No free lunch theorems for optimization.] IEEE Transactions on Evolutionary Computation, vol. 1, no. 1, pp. 67-82, April 1997, doi: 10.1109/4235.585893.
+
+[endsect] [/section:random_search]

example/naive_monte_carlo_example.cpp

@@ -52,6 +52,7 @@ void display_progress(double progress,
 
 int main()
 {
+    using std::abs;
     double exact = 1.3932039296856768591842462603255;
     double A = 1.0 / boost::math::pow<3>(boost::math::constants::pi<double>());
     auto g = [&](std::vector<double> const & x)

example/random_search_example.cpp

@@ -0,0 +1,79 @@
+/*
+ * 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)
+ */
+#if __APPLE__ || __linux__
+#include <signal.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <future>
+#include <chrono>
+#include <iostream>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/optimization/random_search.hpp>
+
+using boost::math::optimization::random_search_parameters;
+using boost::math::optimization::random_search;
+using namespace std::chrono_literals;
+
+template <class Real> Real rastrigin(std::vector<Real> const &v) {
+  using std::cos;
+  using boost::math::constants::two_pi;
+  Real A = 10;
+  Real y = 10 * v.size();
+  for (auto x : v) {
+    y += x * x - A * cos(two_pi<Real>() * x);
+  }
+  return y;
+}
+
+std::atomic<bool> cancel = false;
+
+void ctrl_c_handler(int){
+    cancel = true;
+    std::cout << "Cancellation requested-this could take a second . . ." << std::endl;
+}
+
+int main() {
+  std::cout << "Running random search on Rastrigin function (global minimum = (0,0,...,0))\n";
+  signal(SIGINT, ctrl_c_handler);
+  using ArgType = std::vector<double>;
+  auto params = random_search_parameters<ArgType>();
+  params.lower_bounds.resize(3, -5.12);
+  params.upper_bounds.resize(3, 5.12);
+  params.max_function_calls = 100000000;
+  // Leave one thread available to respond to ctrl-C:
+  params.threads = std::thread::hardware_concurrency() - 1;
+  std::random_device rd;
+  std::mt19937_64 gen(rd());
+
+  // By definition, the value of the function which a target value is provided must be <= target_value.
+  double target_value = 1e-3;
+  std::atomic<double> current_minimum_cost;
+  std::cout << "Hit ctrl-C to gracefully terminate the optimization." << std::endl;
+  auto f = [&]() {
+    return random_search(rastrigin<double>, params, gen, target_value, &cancel, &current_minimum_cost);
+  };
+  auto future = std::async(std::launch::async, f);
+  std::future_status status = future.wait_for(3ms);
+  while (!cancel && (status != std::future_status::ready)) {
+    status = future.wait_for(3ms);
+    std::cout << "Current cost is " << current_minimum_cost << "\r";
+  }
+
+  auto local_minima = future.get();
+  std::cout << "Local minimum is {";
+  for (size_t i = 0; i < local_minima.size() - 1; ++i) {
+    std::cout << local_minima[i] << ", ";
+  }
+  std::cout << local_minima.back() << "}.\n";
+  std::cout << "Final cost: " << current_minimum_cost << "\n";
+}
+#else
+#warning "Signal handling for the random search example only works on Linux and Mac."
+int main() {
+    return 0;
+}
+#endif

include/boost/math/optimization/jso.hpp

@@ -173,10 +173,10 @@ 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, 0.3);
+  std::vector<Real> M_F(H, static_cast<Real>(0.3));
   // Algorithm 1: jSO algorithm, Line 4:
   // "Set all values in M_CR to 0.8":
-  std::vector<Real> M_CR(H, 0.8);
+  std::vector<Real> M_CR(H, static_cast<Real>(0.8));
 
   std::uniform_real_distribution<Real> unif01(Real(0), Real(1));
   bool keep_going = !target_attained;
@@ -203,17 +203,17 @@ 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 = 0.9;
-      Real mu_CR = 0.9;
+      Real mu_F = static_cast<Real>(0.9);
+      Real mu_CR = static_cast<Real>(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 = 0.0;
+      Real crossover_probability = static_cast<Real>(0);
       if (mu_CR >= 0) {
         using std::normal_distribution;
-        normal_distribution<Real> normal(mu_CR, 0.1);
+        normal_distribution<Real> normal(mu_CR, static_cast<Real>(0.1));
         crossover_probability = normal(gen);
         // Clamp comes from L-SHADE description:
         crossover_probability = clamp(crossover_probability, Real(0), Real(1));
@@ -233,7 +233,7 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
       // Algorithm 1, jSO, Line 24-27:
       // Note the adjustments to the pseudocode given in the reference
       // implementation.
-      cauchy_distribution<Real> cauchy(mu_F, 0.1);
+      cauchy_distribution<Real> cauchy(mu_F, static_cast<Real>(0.1));
       Real F;
       do {
         F = cauchy(gen);
@@ -253,13 +253,13 @@ jso(const Func cost_function, jso_parameters<ArgumentContainer> &jso_params,
       Real p = Real(0.25) * (1 - static_cast<Real>(function_evaluations) /
                                      (2 * jso_params.max_function_evaluations));
       // Equation (4) of the reference:
-      Real Fw = 1.2 * F;
+      Real Fw = static_cast<Real>(1.2) * F;
       if (10 * function_evaluations < 4 * jso_params.max_function_evaluations) {
         if (10 * function_evaluations <
             2 * jso_params.max_function_evaluations) {
-          Fw = 0.7 * F;
+          Fw = static_cast<Real>(0.7) * F;
         } else {
-          Fw = 0.8 * F;
+          Fw = static_cast<Real>(0.8) * F;
         }
       }
       // Algorithm 1, jSO, Line 28:

include/boost/math/optimization/random_search.hpp

@@ -0,0 +1,148 @@
+/*
+ * 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_RANDOM_SEARCH_HPP
+#define BOOST_MATH_OPTIMIZATION_RANDOM_SEARCH_HPP
+#include <algorithm>
+#include <array>
+#include <atomic>
+#include <cmath>
+#include <limits>
+#include <list>
+#include <mutex>
+#include <random>
+#include <sstream>
+#include <stdexcept>
+#include <thread>
+#include <type_traits>
+#include <utility>
+#include <vector>
+#include <boost/math/optimization/detail/common.hpp>
+
+namespace boost::math::optimization {
+
+template <typename ArgumentContainer> struct random_search_parameters {
+  using Real = typename ArgumentContainer::value_type;
+  ArgumentContainer lower_bounds;
+  ArgumentContainer upper_bounds;
+  size_t max_function_calls = 10000*std::thread::hardware_concurrency();
+  ArgumentContainer const *initial_guess = nullptr;
+  unsigned threads = std::thread::hardware_concurrency();
+};
+
+template <typename ArgumentContainer>
+void validate_random_search_parameters(random_search_parameters<ArgumentContainer> const &params) {
+  using std::isfinite;
+  using std::isnan;
+  std::ostringstream oss;
+  detail::validate_bounds(params.lower_bounds, params.upper_bounds);
+  if (params.initial_guess) {
+    detail::validate_initial_guess(*params.initial_guess, params.lower_bounds, params.upper_bounds);
+  }
+  if (params.threads == 0) {
+    oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
+    oss << ": There must be at least one thread.";
+    throw std::invalid_argument(oss.str());
+  }
+}
+
+template <typename ArgumentContainer, class Func, class URBG>
+ArgumentContainer random_search(
+    const Func cost_function,
+    random_search_parameters<ArgumentContainer> const &params,
+    URBG &gen,
+    std::invoke_result_t<Func, ArgumentContainer> target_value = std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),
+    std::atomic<bool> *cancellation = nullptr,
+    std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,
+    std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr)
+ {
+  using Real = typename ArgumentContainer::value_type;
+  using ResultType = std::invoke_result_t<Func, ArgumentContainer>;
+  using std::isnan;
+  using std::uniform_real_distribution;
+  validate_random_search_parameters(params);
+  const size_t dimension = params.lower_bounds.size();
+  std::atomic<bool> target_attained = false;
+  // Unfortunately, the "minimum_cost" variable can either be passed
+  // (for observability) or not (if the user doesn't care).
+  // That makes this a bit awkward . . .
+  std::atomic<ResultType> lowest_cost = std::numeric_limits<ResultType>::infinity();
+
+  ArgumentContainer best_vector;
+  if constexpr (detail::has_resize_v<ArgumentContainer>) {
+    best_vector.resize(dimension, std::numeric_limits<Real>::quiet_NaN());
+  }
+  if (params.initial_guess) {
+    auto initial_cost = cost_function(*params.initial_guess);
+    if (!isnan(initial_cost)) {
+      lowest_cost = initial_cost;
+      best_vector = *params.initial_guess;
+      if (current_minimum_cost) {
+        *current_minimum_cost = initial_cost;
+      }
+    }
+  }
+  std::mutex mt;
+  std::vector<std::thread> thread_pool;
+  std::atomic<size_t> function_calls = 0;
+  for (unsigned j = 0; j < params.threads; ++j) {
+    auto seed = gen();
+    thread_pool.emplace_back([&, seed]() {
+      URBG g(seed);
+      ArgumentContainer trial_vector;
+      // This vector is empty unless the user requests the queries be stored:
+      std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> local_queries;
+      if constexpr (detail::has_resize_v<ArgumentContainer>) {
+          trial_vector.resize(dimension, std::numeric_limits<Real>::quiet_NaN());
+      }
+      while (function_calls < params.max_function_calls) {
+        if (cancellation && *cancellation) {
+            break;
+        }
+        if (target_attained) {
+            break;
+        }
+        // Fill trial vector: 
+        uniform_real_distribution<Real> unif01(Real(0), Real(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);
+        }
+        ResultType trial_cost = cost_function(trial_vector);
+        ++function_calls;
+        if (isnan(trial_cost)) {
+          continue;
+        }
+        if (trial_cost < lowest_cost) {
+          lowest_cost = trial_cost;
+          if (current_minimum_cost) {
+            *current_minimum_cost = trial_cost;
+          }
+          // We expect to need to acquire this lock with decreasing frequency
+          // as the computation proceeds:
+          std::scoped_lock lock(mt);
+          best_vector = trial_vector;
+        }
+        if (queries) {
+          local_queries.push_back(std::make_pair(trial_vector, trial_cost));
+        }
+        if (!isnan(target_value) && trial_cost <= target_value) {
+          target_attained = true;
+        }
+      }
+      if (queries) {
+        std::scoped_lock lock(mt);
+        queries->insert(queries->begin(), local_queries.begin(), local_queries.end());
+      }
+    });
+  }
+  for (auto &thread : thread_pool) {
+    thread.join();
+  }
+  return best_vector;
+}
+
+} // namespace boost::math::optimization
+#endif

test/Jamfile.v2

@@ -1257,6 +1257,7 @@ test-suite interpolators :
    [ run test_standalone_asserts.cpp ]
    [ run differential_evolution_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
    [ run jso_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
+   [ run random_search_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
 ;
 
 test-suite quadrature :

test/cardinal_cubic_b_spline_test.cpp

@@ -5,7 +5,7 @@
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 
 #define BOOST_TEST_MODULE test_cubic_b_spline
-
+#include <cmath>
 #include <random>
 #include <functional>
 #include <boost/random/uniform_real_distribution.hpp>
@@ -21,6 +21,7 @@
 #  include <stdfloat>
 #endif
 
+using std::sin;
 using boost::multiprecision::cpp_bin_float_50;
 using boost::math::constants::third;
 using boost::math::constants::half;
@@ -101,10 +102,10 @@ void test_constant_function()
 {
     std::cout << "Testing that constants are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
     std::vector<Real> v(500);
-    Real constant = 50.2;
+    const Real constant = Real(50.2);
     for (size_t i = 0; i < v.size(); ++i)
     {
-        v[i] = Real(50.2);
+        v[i] = constant;
     }
 
     Real step = static_cast<Real>(0.02);

test/differential_evolution_test.cpp

@@ -56,7 +56,7 @@ void test_parameter_checks() {
   bool caught = false;
   try {
     validate_differential_evolution_parameters(de_params);
-  } catch(std::exception const & e) {
+  } catch(std::exception const &) {
     caught = true;
   }
   CHECK_TRUE(caught);
@@ -65,7 +65,7 @@ void test_parameter_checks() {
   de_params.NP = 1;
   try {
     validate_differential_evolution_parameters(de_params);
-  } catch(std::exception const & e) {
+  } catch(std::exception const &) {
     caught = true;
   }
   CHECK_TRUE(caught);
@@ -129,7 +129,7 @@ template <class Real> void test_rastrigin() {
   }
 
   // By definition, the value of the function which a target value is provided must be <= target_value.
-  Real target_value = 1e-3;
+  auto target_value = static_cast<Real>(1e-3);
   local_minima = differential_evolution(rastrigin<Real>, de_params, gen, target_value);
   CHECK_LE(rastrigin(local_minima), target_value);
 }

test/fourier_transform_daubechies_test.cpp

@@ -142,7 +142,6 @@ void test_ten_lectures_eq_5_1_19() {
 
            sum += tpl;
            sum += tml;
-           Real term = tpl + tml;
            ++l;
        }
        // With arg promotion, I can get this to 13 ULPS:

test/gauss_quadrature_test.cpp

@@ -293,7 +293,12 @@ void test_ca()
     Real tol = expected_error<Points>(test_ca_error_id);
     Real L1;
 
-    auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
+    auto f1 = [](const Real& x) { 
+      if (x == 0) {
+         return static_cast<Real>(1);
+      }
+      return atan(x)/(x*(x*x + 1)) ;
+    };
     Real Q = gauss<Real, Points>::integrate(f1, 0, 1, &L1);
     Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
@@ -415,6 +420,7 @@ void test_left_limit_infinite()
 template<class Complex>
 void test_complex_lambert_w()
 {
+    #ifndef BOOST_MATH_STANDALONE
     std::cout << "Testing that complex-valued integrands are integrated correctly by Gaussian quadrature on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n";
     typedef typename Complex::value_type Real;
     Real tol = 10e-9;
@@ -438,7 +444,6 @@ void test_complex_lambert_w()
     };
 
     //N[ProductLog[2+3*I], 150]
-    #ifndef BOOST_MATH_STANDALONE
     Complex Q = gauss<Real, 30>::integrate(lw, (Real) 0, pi<Real>());
     BOOST_CHECK_CLOSE_FRACTION(Q.real(), BOOST_MATH_TEST_VALUE(Real, 1.0900765344857908463017778267816696498710210863535777805644), tol);
     BOOST_CHECK_CLOSE_FRACTION(Q.imag(), BOOST_MATH_TEST_VALUE(Real, 0.5301397207748388014268602135741217419287056313827031782979), tol);

test/math_unit_test.hpp

@@ -11,6 +11,7 @@
 #include <iomanip>
 #include <cmath> // for std::isnan
 #include <string>
+#include <type_traits>
 #include <boost/math/tools/assert.hpp>
 #include <boost/math/special_functions/next.hpp>
 #include <boost/math/special_functions/trunc.hpp>
@@ -155,26 +156,28 @@ bool check_le(Real lesser, Real greater, std::string const & filename, std::stri
     using std::abs;
     using std::isnan;
 
-    if (isnan(lesser))
-    {
-        std::ios_base::fmtflags f( std::cerr.flags() );
-        std::cerr << std::setprecision(3);
-        std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n"
-                  << " \033[0m Lesser value is a nan\n";
-        std::cerr.flags(f);
-        ++detail::global_error_count;
-        return false;
-    }
+    if (std::is_floating_point<Real>::value) {
+        if (boost::math::isnan(lesser))
+        {
+            std::ios_base::fmtflags f( std::cerr.flags() );
+            std::cerr << std::setprecision(3);
+            std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n"
+                    << " \033[0m Lesser value is a nan\n";
+            std::cerr.flags(f);
+            ++detail::global_error_count;
+            return false;
+        }
 
-    if (isnan(greater))
-    {
-        std::ios_base::fmtflags f( std::cerr.flags() );
-        std::cerr << std::setprecision(3);
-        std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n"
-                  << " \033[0m Greater value is a nan\n";
-        std::cerr.flags(f);
-        ++detail::global_error_count;
-        return false;
+        if (boost::math::isnan(greater))
+        {
+            std::ios_base::fmtflags f( std::cerr.flags() );
+            std::cerr << std::setprecision(3);
+            std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n"
+                    << " \033[0m Greater value is a nan\n";
+            std::cerr.flags(f);
+            ++detail::global_error_count;
+            return false;
+        }
     }
 
     if (lesser > greater)

test/random_search_test.cpp

@@ -0,0 +1,156 @@
+/*
+ * 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)
+ */
+
+#include "math_unit_test.hpp"
+#include "test_functions_for_optimization.hpp"
+#include <boost/math/optimization/random_search.hpp>
+#include <random>
+#include <limits>
+
+using boost::math::optimization::random_search;
+using boost::math::optimization::random_search_parameters;
+
+
+template <class Real> void test_ackley() {
+  std::cout << "Testing random search on Ackley function . . .\n";
+  using ArgType = std::array<Real, 2>;
+  auto rs_params = random_search_parameters<ArgType>();
+  rs_params.lower_bounds = {-5, -5};
+  rs_params.upper_bounds = {5, 5};
+
+  std::mt19937_64 gen(12345);
+  auto local_minima = random_search(ackley<Real>, rs_params, gen);
+  CHECK_LE(std::abs(local_minima[0]), Real(0.1));
+  CHECK_LE(std::abs(local_minima[1]), Real(0.1));
+
+  // Does it work with a lambda?
+  auto ack = [](std::array<Real, 2> const &x) { return ackley<Real>(x); };
+  local_minima = random_search(ack, rs_params, gen);
+  CHECK_LE(std::abs(local_minima[0]), Real(0.1));
+  CHECK_LE(std::abs(local_minima[1]), Real(0.1));
+
+  // Test that if an intial guess is the exact solution, the returned solution is the exact solution:
+  std::array<Real, 2> initial_guess{0, 0};
+  rs_params.initial_guess = &initial_guess;
+  local_minima = random_search(ack, rs_params, gen);
+  CHECK_EQUAL(local_minima[0], Real(0));
+  CHECK_EQUAL(local_minima[1], Real(0));
+
+  std::atomic<bool> cancel = false;
+  Real target_value = 0.0;
+  std::atomic<Real> current_minimum_cost = std::numeric_limits<Real>::quiet_NaN();
+  // Test query storage:
+  std::vector<std::pair<ArgType, Real>> queries;
+  local_minima = random_search(ack, rs_params, gen, target_value, &cancel, &current_minimum_cost, &queries);
+  CHECK_EQUAL(local_minima[0], Real(0));
+  CHECK_EQUAL(local_minima[1], Real(0));
+  CHECK_LE(size_t(1), queries.size());
+  for (auto const & q : queries) {
+    auto expected = ackley<Real>(q.first);
+    CHECK_EQUAL(expected, q.second);
+  }
+}
+
+
+template <class Real> void test_rosenbrock_saddle() {
+  std::cout << "Testing random search on Rosenbrock saddle . . .\n";
+  using ArgType = std::array<Real, 2>;
+  auto rs_params = random_search_parameters<ArgType>();
+  rs_params.lower_bounds = {0.5, 0.5};
+  rs_params.upper_bounds = {2.048, 2.048};
+  rs_params.max_function_calls = 20000;
+  std::mt19937_64 gen(234568);
+  auto local_minima = random_search(rosenbrock_saddle<Real>, rs_params, gen);
+
+  CHECK_ABSOLUTE_ERROR(Real(1), local_minima[0], Real(0.05));
+  CHECK_ABSOLUTE_ERROR(Real(1), local_minima[1], Real(0.05));
+
+  // Does cancellation work?
+  std::atomic<bool> cancel = true;
+  gen.seed(12345);
+  local_minima =
+      random_search(rosenbrock_saddle<Real>, rs_params, gen, std::numeric_limits<Real>::quiet_NaN(), &cancel);
+  CHECK_GE(std::abs(local_minima[0] - Real(1)), std::sqrt(std::numeric_limits<Real>::epsilon()));
+}
+
+
+template <class Real> void test_rastrigin() {
+  std::cout << "Testing random search on Rastrigin function (global minimum = (0,0,...,0))\n";
+  using ArgType = std::vector<Real>;
+  auto rs_params = random_search_parameters<ArgType>();
+  rs_params.lower_bounds.resize(3, static_cast<Real>(-5.12));
+  rs_params.upper_bounds.resize(3, static_cast<Real>(5.12));
+  rs_params.max_function_calls = 1000000;
+  std::mt19937_64 gen(34567);
+
+  // By definition, the value of the function which a target value is provided must be <= target_value.
+  Real target_value = 2.0;
+  auto local_minima = random_search(rastrigin<Real>, rs_params, gen, target_value);
+  CHECK_LE(rastrigin(local_minima), target_value);
+}
+
+
+// Tests NaN return types and return type != input type:
+void test_sphere() {
+  std::cout << "Testing random search on sphere . . .\n";
+  using ArgType = std::vector<float>;
+  auto rs_params = random_search_parameters<ArgType>();
+  rs_params.lower_bounds.resize(4, -1);
+  rs_params.upper_bounds.resize(4, 1);
+  rs_params.max_function_calls = 100000;
+  std::mt19937_64 gen(56789);
+  auto local_minima = random_search(sphere, rs_params, gen);
+  for (auto x : local_minima) {
+    CHECK_ABSOLUTE_ERROR(0.0f, x, 0.5f);
+  }
+}
+
+
+template<typename Real>
+void test_three_hump_camel() {
+  std::cout << "Testing random search on three hump camel . . .\n";
+  using ArgType = std::array<Real, 2>;
+  auto rs_params = random_search_parameters<ArgType>();
+  rs_params.lower_bounds[0] = -5.0;
+  rs_params.lower_bounds[1] = -5.0;
+  rs_params.upper_bounds[0] = 5.0;
+  rs_params.upper_bounds[1] = 5.0;
+  std::mt19937_64 gen(56789);
+  auto local_minima = random_search(three_hump_camel<Real>, rs_params, gen);
+  for (auto x : local_minima) {
+    CHECK_ABSOLUTE_ERROR(0.0f, x, 0.2f);
+  }
+}
+
+
+template<typename Real>
+void test_beale() {
+  std::cout << "Testing random search on the Beale function . . .\n";
+  using ArgType = std::array<Real, 2>;
+  auto rs_params = random_search_parameters<ArgType>();
+  rs_params.lower_bounds[0] = -5.0;
+  rs_params.lower_bounds[1] = -5.0;
+  rs_params.upper_bounds[0]= 5.0;
+  rs_params.upper_bounds[1]= 5.0;
+  std::mt19937_64 gen(56789);
+  auto local_minima = random_search(beale<Real>, rs_params, gen);
+  CHECK_ABSOLUTE_ERROR(Real(3), local_minima[0], Real(0.1));
+  CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(0.1));
+}
+
+int main() {
+#if defined(__clang__) || defined(_MSC_VER)
+  test_ackley<float>();
+  test_ackley<double>();
+  test_rosenbrock_saddle<double>();
+  test_rastrigin<double>();
+  test_three_hump_camel<float>();
+  test_beale<double>();
+#endif
+  test_sphere();
+  return boost::math::test::report_errors();
+}

test/test_functions_for_optimization.hpp

@@ -32,8 +32,8 @@ template <typename Real> auto rosenbrock_saddle(std::array<Real, 2> const &v) {
 template <class Real> Real rastrigin(std::vector<Real> const &v) {
   using std::cos;
   using boost::math::constants::two_pi;
-  Real A = 10;
-  Real y = 10 * v.size();
+  auto A = static_cast<Real>(10);
+  auto y = static_cast<Real>(10 * v.size());
   for (auto x : v) {
     y += x * x - A * cos(two_pi<Real>() * x);
   }

tools/hyp_2f0_data.cpp

@@ -46,7 +46,6 @@ int main(int, char* [])
       "This program will generate spot tests for 2F0:\n";
 
    std::string line;
-   bool cont;
 
 #if 1
    arg1 = make_periodic_param(mp_t(-20), mp_t(-1), 19);
@@ -58,6 +57,7 @@ int main(int, char* [])
 
 #else
 
+   bool cont;
    do {
       get_user_parameter_info(arg1, "a1");
       get_user_parameter_info(arg2, "a2");

tools/hyp_2f1_data.cpp

@@ -59,7 +59,6 @@ int main(int, char* [])
       "This program will generate spot tests for 2F0:\n";
 
    std::string line;
-   bool cont;
 
    std::vector<mp_t> v;
    random_ns::mt19937 rnd;
@@ -97,4 +96,3 @@ int main(int, char* [])
    return 0;
 }
 
-

tools/igamma_temme_large_coef.cpp

@@ -134,7 +134,7 @@ void calculate_terms(double sigma, double a, unsigned bits)
    }
    cout << "Max n required:  " << max_n << endl;
 
-   unsigned max_k;
+   unsigned max_k = 0;
    for(unsigned k = 1; k < 10000; ++k)
    {
       double term = tools::real_cast<double>(Coeff(0, k) * pow(a, -((double)k)));