201 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_DETAIL_COMMON_HPP8#define BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP9#include <algorithm> // for std::sort10#include <cmath>11#include <limits>12#include <sstream>13#include <stdexcept>14#include <random>15#include <type_traits> // for std::false_type16 17namespace boost::math::optimization::detail {18 19template <typename T, typename = void> struct has_resize : std::false_type {};20 21template <typename T>22struct has_resize<T, std::void_t<decltype(std::declval<T>().resize(size_t{}))>> : std::true_type {};23 24template <typename T> constexpr bool has_resize_v = has_resize<T>::value;25 26template <typename ArgumentContainer>27void validate_bounds(ArgumentContainer const &lower_bounds, ArgumentContainer const &upper_bounds) {28 using std::isfinite;29 std::ostringstream oss;30 if (lower_bounds.size() == 0) {31 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;32 oss << ": The dimension of the problem cannot be zero.";33 throw std::domain_error(oss.str());34 }35 if (upper_bounds.size() != lower_bounds.size()) {36 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;37 oss << ": There must be the same number of lower bounds as upper bounds, but given ";38 oss << upper_bounds.size() << " upper bounds, and " << lower_bounds.size() << " lower bounds.";39 throw std::domain_error(oss.str());40 }41 for (size_t i = 0; i < lower_bounds.size(); ++i) {42 auto lb = lower_bounds[i];43 auto ub = upper_bounds[i];44 if (lb > ub) {45 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;46 oss << ": The upper bound must be greater than or equal to the lower bound, but the upper bound is " << ub47 << " and the lower is " << lb << ".";48 throw std::domain_error(oss.str());49 }50 if (!isfinite(lb)) {51 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;52 oss << ": The lower bound must be finite, but got " << lb << ".";53 oss << " For infinite bounds, emulate with std::numeric_limits<Real>::lower() or use a standard infinite->finite "54 "transform.";55 throw std::domain_error(oss.str());56 }57 if (!isfinite(ub)) {58 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;59 oss << ": The upper bound must be finite, but got " << ub << ".";60 oss << " For infinite bounds, emulate with std::numeric_limits<Real>::max() or use a standard infinite->finite "61 "transform.";62 throw std::domain_error(oss.str());63 }64 }65}66 67template <typename ArgumentContainer, class URBG>68std::vector<ArgumentContainer> random_initial_population(ArgumentContainer const &lower_bounds,69 ArgumentContainer const &upper_bounds,70 size_t initial_population_size, URBG &&gen) {71 using Real = typename ArgumentContainer::value_type;72 using DimensionlessReal = decltype(Real()/Real());73 constexpr bool has_resize = detail::has_resize_v<ArgumentContainer>;74 std::vector<ArgumentContainer> population(initial_population_size);75 auto const dimension = lower_bounds.size();76 for (size_t i = 0; i < population.size(); ++i) {77 if constexpr (has_resize) {78 population[i].resize(dimension);79 } else {80 // Argument type must be known at compile-time; like std::array:81 if (population[i].size() != dimension) {82 std::ostringstream oss;83 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;84 oss << ": For containers which do not have resize, the default size must be the same as the dimension, ";85 oss << "but the default container size is " << population[i].size() << " and the dimension of the problem is "86 << dimension << ".";87 oss << " The function argument container type is " << typeid(ArgumentContainer).name() << ".\n";88 throw std::runtime_error(oss.str());89 }90 }91 }92 93 // Why don't we provide an option to initialize with (say) a Gaussian distribution?94 // > If the optimum's location is fairly well known,95 // > a Gaussian distribution may prove somewhat faster, although it96 // > may also increase the probability that the population will converge prematurely.97 // > In general, uniform distributions are preferred, since they best reflect98 // > the lack of knowledge about the optimum's location.99 // - Differential Evolution: A Practical Approach to Global Optimization100 // That said, scipy uses Latin Hypercube sampling and says self-avoiding sequences are preferable.101 // So this is something that could be investigated and potentially improved.102 using std::uniform_real_distribution;103 uniform_real_distribution<DimensionlessReal> dis(DimensionlessReal(0), DimensionlessReal(1));104 for (size_t i = 0; i < population.size(); ++i) {105 for (size_t j = 0; j < dimension; ++j) {106 auto const &lb = lower_bounds[j];107 auto const &ub = upper_bounds[j];108 population[i][j] = lb + dis(gen) * (ub - lb);109 }110 }111 112 return population;113}114 115template <typename ArgumentContainer>116void validate_initial_guess(ArgumentContainer const &initial_guess, ArgumentContainer const &lower_bounds,117 ArgumentContainer const &upper_bounds) {118 using std::isfinite;119 std::ostringstream oss;120 auto const dimension = lower_bounds.size();121 if (initial_guess.size() != dimension) {122 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;123 oss << ": The initial guess must have the same dimensions as the problem,";124 oss << ", but the problem size is " << dimension << " and the initial guess has " << initial_guess.size()125 << " elements.";126 throw std::domain_error(oss.str());127 }128 for (size_t i = 0; i < dimension; ++i) {129 auto lb = lower_bounds[i];130 auto ub = upper_bounds[i];131 if (!isfinite(initial_guess[i])) {132 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;133 oss << ": At index " << i << ", the initial guess is " << initial_guess[i]134 << ", make sure all elements of the initial guess are finite.";135 throw std::domain_error(oss.str());136 }137 if (initial_guess[i] < lb || initial_guess[i] > ub) {138 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;139 oss << ": At index " << i << " the initial guess " << initial_guess[i] << " is not in the bounds [" << lb << ", "140 << ub << "].";141 throw std::domain_error(oss.str());142 }143 }144}145 146// Return indices corresponding to the minimum function values.147template <typename Real> std::vector<size_t> best_indices(std::vector<Real> const &function_values) {148 using std::isnan;149 const size_t n = function_values.size();150 std::vector<size_t> indices(n);151 for (size_t i = 0; i < n; ++i) {152 indices[i] = i;153 }154 155 std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) {156 if (isnan(function_values[a])) {157 return false;158 }159 if (isnan(function_values[b])) {160 return true;161 }162 return function_values[a] < function_values[b];163 });164 return indices;165}166 167template<typename RandomAccessContainer>168auto weighted_lehmer_mean(RandomAccessContainer const & values, RandomAccessContainer const & weights) {169 using std::isfinite;170 if (values.size() != weights.size()) {171 std::ostringstream oss;172 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;173 oss << ": There must be the same number of weights as values, but got " << values.size() << " values and " << weights.size() << " weights.";174 throw std::logic_error(oss.str());175 }176 if (values.size() == 0) {177 std::ostringstream oss;178 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;179 oss << ": There must at least one value provided.";180 throw std::logic_error(oss.str());181 }182 using Real = typename RandomAccessContainer::value_type;183 Real numerator = 0;184 Real denominator = 0;185 for (size_t i = 0; i < values.size(); ++i) {186 if (weights[i] < 0 || !isfinite(weights[i])) {187 std::ostringstream oss;188 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;189 oss << ": All weights must be positive and finite, but got received weight " << weights[i] << " at index " << i << " of " << weights.size() << ".";190 throw std::domain_error(oss.str());191 }192 Real tmp = weights[i]*values[i];193 numerator += tmp*values[i];194 denominator += tmp;195 }196 return numerator/denominator;197}198 199} // namespace boost::math::optimization::detail200#endif201