math/optimization/differential_evolution.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_DIFFERENTIAL_EVOLUTION_HPP
0008 #define BOOST_MATH_OPTIMIZATION_DIFFERENTIAL_EVOLUTION_HPP
0009 #include <atomic>
0010 #include <boost/math/optimization/detail/common.hpp>
0011 #include <cmath>
0012 #include <limits>
0013 #include <mutex>
0014 #include <random>
0015 #include <sstream>
0016 #include <stdexcept>
0017 #include <thread>
0018 #include <utility>
0019 #include <vector>
0020
0021 namespace boost::math::optimization {
0022
0023 // Storn, R., Price, K. (1997). Differential evolution-a simple and efficient heuristic for global optimization over
0024 // continuous spaces.
0025 // Journal of global optimization, 11, 341-359.
0026 // See:
0027 // https://www.cp.eng.chula.ac.th/~prabhas//teaching/ec/ec2012/storn_price_de.pdf
0028
0029 // We provide the parameters in a struct-there are too many of them and they are too unwieldy to pass individually:
0030 template <typename ArgumentContainer> struct differential_evolution_parameters {
0031   using Real = typename ArgumentContainer::value_type;
0032   using DimensionlessReal = decltype(Real()/Real());
0033   ArgumentContainer lower_bounds;
0034   ArgumentContainer upper_bounds;
0035   // mutation factor is also called scale factor or just F in the literature:
0036   DimensionlessReal mutation_factor = static_cast<DimensionlessReal>(0.65);
0037   DimensionlessReal crossover_probability = static_cast<DimensionlessReal>(0.5);
0038   // Population in each generation:
0039   size_t NP = 500;
0040   size_t max_generations = 1000;
0041   ArgumentContainer const *initial_guess = nullptr;
0042   unsigned threads = std::thread::hardware_concurrency();
0043 };
0044
0045 template <typename ArgumentContainer>
0046 void validate_differential_evolution_parameters(differential_evolution_parameters<ArgumentContainer> const &de_params) {
0047   using std::isfinite;
0048   using std::isnan;
0049   std::ostringstream oss;
0050   detail::validate_bounds(de_params.lower_bounds, de_params.upper_bounds);
0051   if (de_params.NP < 4) {
0052     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0053     oss << ": The population size must be at least 4, but requested population size of " << de_params.NP << ".";
0054     throw std::invalid_argument(oss.str());
0055   }
0056   // From: "Differential Evolution: A Practical Approach to Global Optimization (Natural Computing Series)"
0057   // > The scale factor, F in (0,1+), is a positive real number that controls the rate at which the population evolves.
0058   // > While there is no upper limit on F, effective values are seldom greater than 1.0.
0059   // ...
0060   // Also see "Limits on F", Section 2.5.1:
0061   // > This discontinuity at F = 1 reduces the number of mutants by half and can result in erratic convergence...
0062   auto F = de_params.mutation_factor;
0063   if (isnan(F) || F >= 1 || F <= 0) {
0064     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0065     oss << ": F in (0, 1) is required, but got F=" << F << ".";
0066     throw std::domain_error(oss.str());
0067   }
0068   if (de_params.max_generations < 1) {
0069     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0070     oss << ": There must be at least one generation.";
0071     throw std::invalid_argument(oss.str());
0072   }
0073   if (de_params.initial_guess) {
0074     detail::validate_initial_guess(*de_params.initial_guess, de_params.lower_bounds, de_params.upper_bounds);
0075   }
0076   if (de_params.threads == 0) {
0077     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0078     oss << ": There must be at least one thread.";
0079     throw std::invalid_argument(oss.str());
0080   }
0081 }
0082
0083 template <typename ArgumentContainer, class Func, class URBG>
0084 ArgumentContainer differential_evolution(
0085     const Func cost_function, differential_evolution_parameters<ArgumentContainer> const &de_params, URBG &gen,
0086     std::invoke_result_t<Func, ArgumentContainer> target_value =
0087         std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),
0088     std::atomic<bool> *cancellation = nullptr,
0089     std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr,
0090     std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr) {
0091   using Real = typename ArgumentContainer::value_type;
0092   using DimensionlessReal = decltype(Real()/Real());
0093   using ResultType = std::invoke_result_t<Func, ArgumentContainer>;
0094   using std::clamp;
0095   using std::isnan;
0096   using std::round;
0097   using std::uniform_real_distribution;
0098   validate_differential_evolution_parameters(de_params);
0099   const size_t dimension = de_params.lower_bounds.size();
0100   auto NP = de_params.NP;
0101   auto population = detail::random_initial_population(de_params.lower_bounds, de_params.upper_bounds, NP, gen);
0102   if (de_params.initial_guess) {
0103     population[0] = *de_params.initial_guess;
0104   }
0105   std::vector<ResultType> cost(NP, std::numeric_limits<ResultType>::quiet_NaN());
0106   std::atomic<bool> target_attained = false;
0107   // This mutex is only used if the queries are stored:
0108   std::mutex mt;
0109
0110   std::vector<std::thread> thread_pool;
0111   auto const threads = de_params.threads;
0112   for (size_t j = 0; j < threads; ++j) {
0113     // Note that if some members of the population take way longer to compute,
0114     // then this parallelization strategy is very suboptimal.
0115     // However, we tried using std::async (which should be robust to this particular problem),
0116     // but the overhead was just totally unacceptable on ARM Macs (the only platform tested).
0117     // As the economists say "there are no solutions, only tradeoffs".
0118     thread_pool.emplace_back([&, j]() {
0119       for (size_t i = j; i < cost.size(); i += threads) {
0120         cost[i] = cost_function(population[i]);
0121         if (current_minimum_cost && cost[i] < *current_minimum_cost) {
0122           *current_minimum_cost = cost[i];
0123         }
0124         if (queries) {
0125           std::scoped_lock lock(mt);
0126           queries->push_back(std::make_pair(population[i], cost[i]));
0127         }
0128         if (!isnan(target_value) && cost[i] <= target_value) {
0129           target_attained = true;
0130         }
0131       }
0132     });
0133   }
0134   for (auto &thread : thread_pool) {
0135     thread.join();
0136   }
0137
0138   std::vector<ArgumentContainer> trial_vectors(NP);
0139   for (size_t i = 0; i < NP; ++i) {
0140     if constexpr (detail::has_resize_v<ArgumentContainer>) {
0141       trial_vectors[i].resize(dimension);
0142     }
0143   }
0144   std::vector<URBG> thread_generators(threads);
0145   for (size_t j = 0; j < threads; ++j) {
0146     thread_generators[j].seed(gen());
0147   }
0148   // std::vector<bool> isn't threadsafe!
0149   std::vector<int> updated_indices(NP, 0);
0150
0151   for (size_t generation = 0; generation < de_params.max_generations; ++generation) {
0152     if (cancellation && *cancellation) {
0153       break;
0154     }
0155     if (target_attained) {
0156       break;
0157     }
0158     thread_pool.resize(0);
0159     for (size_t j = 0; j < threads; ++j) {
0160       thread_pool.emplace_back([&, j]() {
0161         auto& tlg = thread_generators[j];
0162         uniform_real_distribution<DimensionlessReal> unif01(DimensionlessReal(0), DimensionlessReal(1));
0163         for (size_t i = j; i < cost.size(); i += threads) {
0164           if (target_attained) {
0165             return;
0166           }
0167           if (cancellation && *cancellation) {
0168             return;
0169           }
0170           size_t r1, r2, r3;
0171           do {
0172             r1 = tlg() % NP;
0173           } while (r1 == i);
0174           do {
0175             r2 = tlg() % NP;
0176           } while (r2 == i || r2 == r1);
0177           do {
0178             r3 = tlg() % NP;
0179           } while (r3 == i || r3 == r2 || r3 == r1);
0180
0181           for (size_t k = 0; k < dimension; ++k) {
0182             // See equation (4) of the reference:
0183             auto guaranteed_changed_idx = tlg() % dimension;
0184             if (unif01(tlg) < de_params.crossover_probability || k == guaranteed_changed_idx) {
0185               auto tmp = population[r1][k] + de_params.mutation_factor * (population[r2][k] - population[r3][k]);
0186               auto const &lb = de_params.lower_bounds[k];
0187               auto const &ub = de_params.upper_bounds[k];
0188               // Some others recommend regenerating the indices rather than clamping;
0189               // I dunno seems like it could get stuck regenerating . . .
0190               trial_vectors[i][k] = clamp(tmp, lb, ub);
0191             } else {
0192               trial_vectors[i][k] = population[i][k];
0193             }
0194           }
0195
0196           auto const trial_cost = cost_function(trial_vectors[i]);
0197           if (isnan(trial_cost)) {
0198             continue;
0199           }
0200           if (queries) {
0201             std::scoped_lock lock(mt);
0202             queries->push_back(std::make_pair(trial_vectors[i], trial_cost));
0203           }
0204           if (trial_cost < cost[i] || isnan(cost[i])) {
0205             cost[i] = trial_cost;
0206             if (!isnan(target_value) && cost[i] <= target_value) {
0207               target_attained = true;
0208             }
0209             if (current_minimum_cost && cost[i] < *current_minimum_cost) {
0210               *current_minimum_cost = cost[i];
0211             }
0212             // Can't do this! It's a race condition!
0213             //population[i] = trial_vectors[i];
0214             // Instead mark all the indices that need to be updated:
0215             updated_indices[i] = 1;
0216           }
0217         }
0218       });
0219     }
0220     for (auto &thread : thread_pool) {
0221       thread.join();
0222     }
0223     for (size_t i = 0; i < NP; ++i) {
0224       if (updated_indices[i]) {
0225         population[i] = trial_vectors[i];
0226         updated_indices[i] = 0;
0227       }
0228     }
0229   }
0230
0231   auto it = std::min_element(cost.begin(), cost.end());
0232   return population[std::distance(cost.begin(), it)];
0233 }
0234
0235 } // namespace boost::math::optimization
0236 #endif