math/optimization/random_search.hpp

0001 /*
0002  * Copyright Nick Thompson, 2024
0003  * Use, modification and distribution are subject to the
0004  * Boost Software License, Version 1.0. (See accompanying file
0005  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006  */
0007 #ifndef BOOST_MATH_OPTIMIZATION_RANDOM_SEARCH_HPP
0008 #define BOOST_MATH_OPTIMIZATION_RANDOM_SEARCH_HPP
0009 #include <atomic>
0010 #include <cmath>
0011 #include <limits>
0012 #include <mutex>
0013 #include <random>
0014 #include <sstream>
0015 #include <stdexcept>
0016 #include <thread>
0017 #include <utility>
0018 #include <vector>
0019 #include <boost/math/optimization/detail/common.hpp>
0020
0021 namespace boost::math::optimization {
0022
0023 template <typename ArgumentContainer> struct random_search_parameters {
0024   using Real = typename ArgumentContainer::value_type;
0025   ArgumentContainer lower_bounds;
0026   ArgumentContainer upper_bounds;
0027   size_t max_function_calls = 10000*std::thread::hardware_concurrency();
0028   ArgumentContainer const *initial_guess = nullptr;
0029   unsigned threads = std::thread::hardware_concurrency();
0030 };
0031
0032 template <typename ArgumentContainer>
0033 void validate_random_search_parameters(random_search_parameters<ArgumentContainer> const &params) {
0034   using std::isfinite;
0035   using std::isnan;
0036   std::ostringstream oss;
0037   detail::validate_bounds(params.lower_bounds, params.upper_bounds);
0038   if (params.initial_guess) {
0039     detail::validate_initial_guess(*params.initial_guess, params.lower_bounds, params.upper_bounds);
0040   }
0041   if (params.threads == 0) {
0042     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0043     oss << ": There must be at least one thread.";
0044     throw std::invalid_argument(oss.str());
0045   }
0046 }
0047
0048 template <typename ArgumentContainer, class Func, class URBG>
0049 ArgumentContainer random_search(
0050     const Func cost_function,
0051     random_search_parameters<ArgumentContainer> const &params,
0052     URBG &gen,
0053     std::invoke_result_t<Func, ArgumentContainer> target_value = std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),
0054     std::atomic<bool> *cancellation = nullptr,
0055     std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,
0056     std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr)
0057  {
0058   using Real = typename ArgumentContainer::value_type;
0059   using DimensionlessReal = decltype(Real()/Real());
0060   using ResultType = std::invoke_result_t<Func, ArgumentContainer>;
0061   using std::isnan;
0062   using std::uniform_real_distribution;
0063   validate_random_search_parameters(params);
0064   const size_t dimension = params.lower_bounds.size();
0065   std::atomic<bool> target_attained = false;
0066   // Unfortunately, the "minimum_cost" variable can either be passed
0067   // (for observability) or not (if the user doesn't care).
0068   // That makes this a bit awkward . . .
0069   std::atomic<ResultType> lowest_cost = std::numeric_limits<ResultType>::infinity();
0070
0071   ArgumentContainer best_vector;
0072   if constexpr (detail::has_resize_v<ArgumentContainer>) {
0073     best_vector.resize(dimension, std::numeric_limits<Real>::quiet_NaN());
0074   }
0075   if (params.initial_guess) {
0076     auto initial_cost = cost_function(*params.initial_guess);
0077     if (!isnan(initial_cost)) {
0078       lowest_cost = initial_cost;
0079       best_vector = *params.initial_guess;
0080       if (current_minimum_cost) {
0081         *current_minimum_cost = initial_cost;
0082       }
0083     }
0084   }
0085   std::mutex mt;
0086   std::vector<std::thread> thread_pool;
0087   std::atomic<size_t> function_calls = 0;
0088   for (unsigned j = 0; j < params.threads; ++j) {
0089     auto seed = gen();
0090     thread_pool.emplace_back([&, seed]() {
0091       URBG g(seed);
0092       ArgumentContainer trial_vector;
0093       // This vector is empty unless the user requests the queries be stored:
0094       std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> local_queries;
0095       if constexpr (detail::has_resize_v<ArgumentContainer>) {
0096           trial_vector.resize(dimension, std::numeric_limits<Real>::quiet_NaN());
0097       }
0098       while (function_calls < params.max_function_calls) {
0099         if (cancellation && *cancellation) {
0100             break;
0101         }
0102         if (target_attained) {
0103             break;
0104         }
0105         // Fill trial vector:
0106         uniform_real_distribution<DimensionlessReal> unif01(DimensionlessReal(0), DimensionlessReal(1));
0107         for (size_t i = 0; i < dimension; ++i) {
0108             trial_vector[i] = params.lower_bounds[i] + (params.upper_bounds[i] - params.lower_bounds[i])*unif01(g);
0109         }
0110         ResultType trial_cost = cost_function(trial_vector);
0111         ++function_calls;
0112         if (isnan(trial_cost)) {
0113           continue;
0114         }
0115         if (trial_cost < lowest_cost) {
0116           lowest_cost = trial_cost;
0117           if (current_minimum_cost) {
0118             *current_minimum_cost = trial_cost;
0119           }
0120           // We expect to need to acquire this lock with decreasing frequency
0121           // as the computation proceeds:
0122           std::scoped_lock lock(mt);
0123           best_vector = trial_vector;
0124         }
0125         if (queries) {
0126           local_queries.push_back(std::make_pair(trial_vector, trial_cost));
0127         }
0128         if (!isnan(target_value) && trial_cost <= target_value) {
0129           target_attained = true;
0130         }
0131       }
0132       if (queries) {
0133         std::scoped_lock lock(mt);
0134         queries->insert(queries->begin(), local_queries.begin(), local_queries.end());
0135       }
0136     });
0137   }
0138   for (auto &thread : thread_pool) {
0139     thread.join();
0140   }
0141   return best_vector;
0142 }
0143
0144 } // namespace boost::math::optimization
0145 #endif