237 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_DIFFERENTIAL_EVOLUTION_HPP8#define BOOST_MATH_OPTIMIZATION_DIFFERENTIAL_EVOLUTION_HPP9#include <atomic>10#include <boost/math/optimization/detail/common.hpp>11#include <cmath>12#include <limits>13#include <mutex>14#include <random>15#include <sstream>16#include <stdexcept>17#include <thread>18#include <utility>19#include <vector>20 21namespace boost::math::optimization {22 23// Storn, R., Price, K. (1997). Differential evolution-a simple and efficient heuristic for global optimization over24// continuous spaces.25// Journal of global optimization, 11, 341-359.26// See:27// https://www.cp.eng.chula.ac.th/~prabhas//teaching/ec/ec2012/storn_price_de.pdf28 29// We provide the parameters in a struct-there are too many of them and they are too unwieldy to pass individually:30template <typename ArgumentContainer> struct differential_evolution_parameters {31 using Real = typename ArgumentContainer::value_type;32 using DimensionlessReal = decltype(Real()/Real());33 ArgumentContainer lower_bounds;34 ArgumentContainer upper_bounds;35 // mutation factor is also called scale factor or just F in the literature:36 DimensionlessReal mutation_factor = static_cast<DimensionlessReal>(0.65);37 DimensionlessReal crossover_probability = static_cast<DimensionlessReal>(0.5);38 // Population in each generation:39 size_t NP = 500;40 size_t max_generations = 1000;41 ArgumentContainer const *initial_guess = nullptr;42 unsigned threads = std::thread::hardware_concurrency();43};44 45template <typename ArgumentContainer>46void validate_differential_evolution_parameters(differential_evolution_parameters<ArgumentContainer> const &de_params) {47 using std::isfinite;48 using std::isnan;49 std::ostringstream oss;50 detail::validate_bounds(de_params.lower_bounds, de_params.upper_bounds);51 if (de_params.NP < 4) {52 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;53 oss << ": The population size must be at least 4, but requested population size of " << de_params.NP << ".";54 throw std::invalid_argument(oss.str());55 }56 // From: "Differential Evolution: A Practical Approach to Global Optimization (Natural Computing Series)"57 // > The scale factor, F in (0,1+), is a positive real number that controls the rate at which the population evolves.58 // > While there is no upper limit on F, effective values are seldom greater than 1.0.59 // ...60 // Also see "Limits on F", Section 2.5.1:61 // > This discontinuity at F = 1 reduces the number of mutants by half and can result in erratic convergence...62 auto F = de_params.mutation_factor;63 if (isnan(F) || F >= 1 || F <= 0) {64 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;65 oss << ": F in (0, 1) is required, but got F=" << F << ".";66 throw std::domain_error(oss.str());67 }68 if (de_params.max_generations < 1) {69 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;70 oss << ": There must be at least one generation.";71 throw std::invalid_argument(oss.str());72 }73 if (de_params.initial_guess) {74 detail::validate_initial_guess(*de_params.initial_guess, de_params.lower_bounds, de_params.upper_bounds);75 }76 if (de_params.threads == 0) {77 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;78 oss << ": There must be at least one thread.";79 throw std::invalid_argument(oss.str());80 }81}82 83template <typename ArgumentContainer, class Func, class URBG>84ArgumentContainer differential_evolution(85 const Func cost_function, differential_evolution_parameters<ArgumentContainer> const &de_params, URBG &gen,86 std::invoke_result_t<Func, ArgumentContainer> target_value =87 std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),88 std::atomic<bool> *cancellation = nullptr,89 std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,90 std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr) {91 using Real = typename ArgumentContainer::value_type;92 using DimensionlessReal = decltype(Real()/Real());93 using ResultType = std::invoke_result_t<Func, ArgumentContainer>;94 using std::clamp;95 using std::isnan;96 using std::round;97 using std::uniform_real_distribution;98 validate_differential_evolution_parameters(de_params);99 const size_t dimension = de_params.lower_bounds.size();100 auto NP = de_params.NP;101 auto population = detail::random_initial_population(de_params.lower_bounds, de_params.upper_bounds, NP, gen);102 if (de_params.initial_guess) {103 population[0] = *de_params.initial_guess;104 }105 std::vector<ResultType> cost(NP, std::numeric_limits<ResultType>::quiet_NaN());106 std::atomic<bool> target_attained = false;107 // This mutex is only used if the queries are stored:108 std::mutex mt;109 110 std::vector<std::thread> thread_pool;111 auto const threads = de_params.threads;112 for (size_t j = 0; j < threads; ++j) {113 // Note that if some members of the population take way longer to compute,114 // then this parallelization strategy is very suboptimal.115 // However, we tried using std::async (which should be robust to this particular problem),116 // but the overhead was just totally unacceptable on ARM Macs (the only platform tested).117 // As the economists say "there are no solutions, only tradeoffs".118 thread_pool.emplace_back([&, j]() {119 for (size_t i = j; i < cost.size(); i += threads) {120 cost[i] = cost_function(population[i]);121 if (current_minimum_cost && cost[i] < *current_minimum_cost) {122 *current_minimum_cost = cost[i];123 }124 if (queries) {125 std::scoped_lock lock(mt);126 queries->push_back(std::make_pair(population[i], cost[i]));127 }128 if (!isnan(target_value) && cost[i] <= target_value) {129 target_attained = true;130 }131 }132 });133 }134 for (auto &thread : thread_pool) {135 thread.join();136 }137 138 std::vector<ArgumentContainer> trial_vectors(NP);139 for (size_t i = 0; i < NP; ++i) {140 if constexpr (detail::has_resize_v<ArgumentContainer>) {141 trial_vectors[i].resize(dimension);142 }143 }144 std::vector<URBG> thread_generators(threads);145 for (size_t j = 0; j < threads; ++j) {146 thread_generators[j].seed(gen());147 }148 // std::vector<bool> isn't threadsafe!149 std::vector<int> updated_indices(NP, 0);150 151 for (size_t generation = 0; generation < de_params.max_generations; ++generation) {152 if (cancellation && *cancellation) {153 break;154 }155 if (target_attained) {156 break;157 }158 thread_pool.resize(0);159 for (size_t j = 0; j < threads; ++j) {160 thread_pool.emplace_back([&, j]() {161 auto& tlg = thread_generators[j];162 uniform_real_distribution<DimensionlessReal> unif01(DimensionlessReal(0), DimensionlessReal(1));163 for (size_t i = j; i < cost.size(); i += threads) {164 if (target_attained) {165 return;166 }167 if (cancellation && *cancellation) {168 return;169 }170 size_t r1, r2, r3;171 do {172 r1 = tlg() % NP;173 } while (r1 == i);174 do {175 r2 = tlg() % NP;176 } while (r2 == i || r2 == r1);177 do {178 r3 = tlg() % NP;179 } while (r3 == i || r3 == r2 || r3 == r1);180 181 for (size_t k = 0; k < dimension; ++k) {182 // See equation (4) of the reference:183 auto guaranteed_changed_idx = tlg() % dimension;184 if (unif01(tlg) < de_params.crossover_probability || k == guaranteed_changed_idx) {185 auto tmp = population[r1][k] + de_params.mutation_factor * (population[r2][k] - population[r3][k]);186 auto const &lb = de_params.lower_bounds[k];187 auto const &ub = de_params.upper_bounds[k];188 // Some others recommend regenerating the indices rather than clamping;189 // I dunno seems like it could get stuck regenerating . . .190 trial_vectors[i][k] = clamp(tmp, lb, ub);191 } else {192 trial_vectors[i][k] = population[i][k];193 }194 }195 196 auto const trial_cost = cost_function(trial_vectors[i]);197 if (isnan(trial_cost)) {198 continue;199 }200 if (queries) {201 std::scoped_lock lock(mt);202 queries->push_back(std::make_pair(trial_vectors[i], trial_cost));203 }204 if (trial_cost < cost[i] || isnan(cost[i])) {205 cost[i] = trial_cost;206 if (!isnan(target_value) && cost[i] <= target_value) {207 target_attained = true;208 }209 if (current_minimum_cost && cost[i] < *current_minimum_cost) {210 *current_minimum_cost = cost[i];211 }212 // Can't do this! It's a race condition!213 //population[i] = trial_vectors[i];214 // Instead mark all the indices that need to be updated:215 updated_indices[i] = 1;216 }217 }218 });219 }220 for (auto &thread : thread_pool) {221 thread.join();222 }223 for (size_t i = 0; i < NP; ++i) {224 if (updated_indices[i]) {225 population[i] = trial_vectors[i];226 updated_indices[i] = 0;227 }228 }229 }230 231 auto it = std::min_element(cost.begin(), cost.end());232 return population[std::distance(cost.begin(), it)];233}234 235} // namespace boost::math::optimization236#endif237