146 lines · plain
1/*2 * Copyright Nick Thompson, 20243 * Use, modification and distribution are subject to the4 * Boost Software License, Version 1.0. (See accompanying file5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 */7#ifndef BOOST_MATH_OPTIMIZATION_RANDOM_SEARCH_HPP8#define BOOST_MATH_OPTIMIZATION_RANDOM_SEARCH_HPP9#include <atomic>10#include <cmath>11#include <limits>12#include <mutex>13#include <random>14#include <sstream>15#include <stdexcept>16#include <thread>17#include <utility>18#include <vector>19#include <boost/math/optimization/detail/common.hpp>20 21namespace boost::math::optimization {22 23template <typename ArgumentContainer> struct random_search_parameters {24 using Real = typename ArgumentContainer::value_type;25 ArgumentContainer lower_bounds;26 ArgumentContainer upper_bounds;27 size_t max_function_calls = 10000*std::thread::hardware_concurrency();28 ArgumentContainer const *initial_guess = nullptr;29 unsigned threads = std::thread::hardware_concurrency();30};31 32template <typename ArgumentContainer>33void validate_random_search_parameters(random_search_parameters<ArgumentContainer> const ¶ms) {34 using std::isfinite;35 using std::isnan;36 std::ostringstream oss;37 detail::validate_bounds(params.lower_bounds, params.upper_bounds);38 if (params.initial_guess) {39 detail::validate_initial_guess(*params.initial_guess, params.lower_bounds, params.upper_bounds);40 }41 if (params.threads == 0) {42 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;43 oss << ": There must be at least one thread.";44 throw std::invalid_argument(oss.str());45 }46}47 48template <typename ArgumentContainer, class Func, class URBG>49ArgumentContainer random_search(50 const Func cost_function,51 random_search_parameters<ArgumentContainer> const ¶ms,52 URBG &gen,53 std::invoke_result_t<Func, ArgumentContainer> target_value = std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),54 std::atomic<bool> *cancellation = nullptr,55 std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,56 std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr)57 {58 using Real = typename ArgumentContainer::value_type;59 using DimensionlessReal = decltype(Real()/Real());60 using ResultType = std::invoke_result_t<Func, ArgumentContainer>;61 using std::isnan;62 using std::uniform_real_distribution;63 validate_random_search_parameters(params);64 const size_t dimension = params.lower_bounds.size();65 std::atomic<bool> target_attained = false;66 // Unfortunately, the "minimum_cost" variable can either be passed67 // (for observability) or not (if the user doesn't care).68 // That makes this a bit awkward . . .69 std::atomic<ResultType> lowest_cost = std::numeric_limits<ResultType>::infinity();70 71 ArgumentContainer best_vector;72 if constexpr (detail::has_resize_v<ArgumentContainer>) {73 best_vector.resize(dimension, std::numeric_limits<Real>::quiet_NaN());74 }75 if (params.initial_guess) {76 auto initial_cost = cost_function(*params.initial_guess);77 if (!isnan(initial_cost)) {78 lowest_cost = initial_cost;79 best_vector = *params.initial_guess;80 if (current_minimum_cost) {81 *current_minimum_cost = initial_cost;82 }83 }84 }85 std::mutex mt;86 std::vector<std::thread> thread_pool;87 std::atomic<size_t> function_calls = 0;88 for (unsigned j = 0; j < params.threads; ++j) {89 auto seed = gen();90 thread_pool.emplace_back([&, seed]() {91 URBG g(seed);92 ArgumentContainer trial_vector;93 // This vector is empty unless the user requests the queries be stored:94 std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> local_queries;95 if constexpr (detail::has_resize_v<ArgumentContainer>) {96 trial_vector.resize(dimension, std::numeric_limits<Real>::quiet_NaN());97 }98 while (function_calls < params.max_function_calls) {99 if (cancellation && *cancellation) {100 break;101 }102 if (target_attained) {103 break;104 }105 // Fill trial vector: 106 uniform_real_distribution<DimensionlessReal> unif01(DimensionlessReal(0), DimensionlessReal(1));107 for (size_t i = 0; i < dimension; ++i) {108 trial_vector[i] = params.lower_bounds[i] + (params.upper_bounds[i] - params.lower_bounds[i])*unif01(g);109 }110 ResultType trial_cost = cost_function(trial_vector);111 ++function_calls;112 if (isnan(trial_cost)) {113 continue;114 }115 if (trial_cost < lowest_cost) {116 lowest_cost = trial_cost;117 if (current_minimum_cost) {118 *current_minimum_cost = trial_cost;119 }120 // We expect to need to acquire this lock with decreasing frequency121 // as the computation proceeds:122 std::scoped_lock lock(mt);123 best_vector = trial_vector;124 }125 if (queries) {126 local_queries.push_back(std::make_pair(trial_vector, trial_cost));127 }128 if (!isnan(target_value) && trial_cost <= target_value) {129 target_attained = true;130 }131 }132 if (queries) {133 std::scoped_lock lock(mt);134 queries->insert(queries->begin(), local_queries.begin(), local_queries.end());135 }136 });137 }138 for (auto &thread : thread_pool) {139 thread.join();140 }141 return best_vector;142}143 144} // namespace boost::math::optimization145#endif146