optimization/detail/common.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_DETAIL_COMMON_HPP
0008 #define BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP
0009 #include <algorithm> // for std::sort
0010 #include <cmath>
0011 #include <limits>
0012 #include <sstream>
0013 #include <stdexcept>
0014 #include <random>
0015 #include <type_traits>  // for std::false_type
0016
0017 namespace boost::math::optimization::detail {
0018
0019 template <typename T, typename = void> struct has_resize : std::false_type {};
0020
0021 template <typename T>
0022 struct has_resize<T, std::void_t<decltype(std::declval<T>().resize(size_t{}))>> : std::true_type {};
0023
0024 template <typename T> constexpr bool has_resize_v = has_resize<T>::value;
0025
0026 template <typename ArgumentContainer>
0027 void validate_bounds(ArgumentContainer const &lower_bounds, ArgumentContainer const &upper_bounds) {
0028   using std::isfinite;
0029   std::ostringstream oss;
0030   if (lower_bounds.size() == 0) {
0031     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0032     oss << ": The dimension of the problem cannot be zero.";
0033     throw std::domain_error(oss.str());
0034   }
0035   if (upper_bounds.size() != lower_bounds.size()) {
0036     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0037     oss << ": There must be the same number of lower bounds as upper bounds, but given ";
0038     oss << upper_bounds.size() << " upper bounds, and " << lower_bounds.size() << " lower bounds.";
0039     throw std::domain_error(oss.str());
0040   }
0041   for (size_t i = 0; i < lower_bounds.size(); ++i) {
0042     auto lb = lower_bounds[i];
0043     auto ub = upper_bounds[i];
0044     if (lb > ub) {
0045       oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0046       oss << ": The upper bound must be greater than or equal to the lower bound, but the upper bound is " << ub
0047           << " and the lower is " << lb << ".";
0048       throw std::domain_error(oss.str());
0049     }
0050     if (!isfinite(lb)) {
0051       oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0052       oss << ": The lower bound must be finite, but got " << lb << ".";
0053       oss << " For infinite bounds, emulate with std::numeric_limits<Real>::lower() or use a standard infinite->finite "
0054              "transform.";
0055       throw std::domain_error(oss.str());
0056     }
0057     if (!isfinite(ub)) {
0058       oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0059       oss << ": The upper bound must be finite, but got " << ub << ".";
0060       oss << " For infinite bounds, emulate with std::numeric_limits<Real>::max() or use a standard infinite->finite "
0061              "transform.";
0062       throw std::domain_error(oss.str());
0063     }
0064   }
0065 }
0066
0067 template <typename ArgumentContainer, class URBG>
0068 std::vector<ArgumentContainer> random_initial_population(ArgumentContainer const &lower_bounds,
0069                                                          ArgumentContainer const &upper_bounds,
0070                                                          size_t initial_population_size, URBG &&gen) {
0071   using Real = typename ArgumentContainer::value_type;
0072   using DimensionlessReal = decltype(Real()/Real());
0073   constexpr bool has_resize = detail::has_resize_v<ArgumentContainer>;
0074   std::vector<ArgumentContainer> population(initial_population_size);
0075   auto const dimension = lower_bounds.size();
0076   for (size_t i = 0; i < population.size(); ++i) {
0077     if constexpr (has_resize) {
0078       population[i].resize(dimension);
0079     } else {
0080       // Argument type must be known at compile-time; like std::array:
0081       if (population[i].size() != dimension) {
0082         std::ostringstream oss;
0083         oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0084         oss << ": For containers which do not have resize, the default size must be the same as the dimension, ";
0085         oss << "but the default container size is " << population[i].size() << " and the dimension of the problem is "
0086             << dimension << ".";
0087         oss << " The function argument container type is " << typeid(ArgumentContainer).name() << ".\n";
0088         throw std::runtime_error(oss.str());
0089       }
0090     }
0091   }
0092
0093   // Why don't we provide an option to initialize with (say) a Gaussian distribution?
0094   // > If the optimum's location is fairly well known,
0095   // > a Gaussian distribution may prove somewhat faster, although it
0096   // > may also increase the probability that the population will converge prematurely.
0097   // > In general, uniform distributions are preferred, since they best reflect
0098   // > the lack of knowledge about the optimum's location.
0099   //  - Differential Evolution: A Practical Approach to Global Optimization
0100   // That said, scipy uses Latin Hypercube sampling and says self-avoiding sequences are preferable.
0101   // So this is something that could be investigated and potentially improved.
0102   using std::uniform_real_distribution;
0103   uniform_real_distribution<DimensionlessReal> dis(DimensionlessReal(0), DimensionlessReal(1));
0104   for (size_t i = 0; i < population.size(); ++i) {
0105     for (size_t j = 0; j < dimension; ++j) {
0106       auto const &lb = lower_bounds[j];
0107       auto const &ub = upper_bounds[j];
0108       population[i][j] = lb + dis(gen) * (ub - lb);
0109     }
0110   }
0111
0112   return population;
0113 }
0114
0115 template <typename ArgumentContainer>
0116 void validate_initial_guess(ArgumentContainer const &initial_guess, ArgumentContainer const &lower_bounds,
0117                             ArgumentContainer const &upper_bounds) {
0118   using std::isfinite;
0119   std::ostringstream oss;
0120   auto const dimension = lower_bounds.size();
0121   if (initial_guess.size() != dimension) {
0122     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0123     oss << ": The initial guess must have the same dimensions as the problem,";
0124     oss << ", but the problem size is " << dimension << " and the initial guess has " << initial_guess.size()
0125         << " elements.";
0126     throw std::domain_error(oss.str());
0127   }
0128   for (size_t i = 0; i < dimension; ++i) {
0129     auto lb = lower_bounds[i];
0130     auto ub = upper_bounds[i];
0131     if (!isfinite(initial_guess[i])) {
0132       oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0133       oss << ": At index " << i << ", the initial guess is " << initial_guess[i]
0134           << ", make sure all elements of the initial guess are finite.";
0135       throw std::domain_error(oss.str());
0136     }
0137     if (initial_guess[i] < lb || initial_guess[i] > ub) {
0138       oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0139       oss << ": At index " << i << " the initial guess " << initial_guess[i] << " is not in the bounds [" << lb << ", "
0140           << ub << "].";
0141       throw std::domain_error(oss.str());
0142     }
0143   }
0144 }
0145
0146 // Return indices corresponding to the minimum function values.
0147 template <typename Real> std::vector<size_t> best_indices(std::vector<Real> const &function_values) {
0148   using std::isnan;
0149   const size_t n = function_values.size();
0150   std::vector<size_t> indices(n);
0151   for (size_t i = 0; i < n; ++i) {
0152     indices[i] = i;
0153   }
0154
0155   std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) {
0156     if (isnan(function_values[a])) {
0157       return false;
0158     }
0159     if (isnan(function_values[b])) {
0160       return true;
0161     }
0162     return function_values[a] < function_values[b];
0163   });
0164   return indices;
0165 }
0166
0167 template<typename RandomAccessContainer>
0168 auto weighted_lehmer_mean(RandomAccessContainer const & values, RandomAccessContainer const & weights) {
0169   using std::isfinite;
0170   if (values.size() != weights.size()) {
0171     std::ostringstream oss;
0172     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0173     oss << ": There must be the same number of weights as values, but got " << values.size() << " values and " << weights.size() << " weights.";
0174     throw std::logic_error(oss.str());
0175   }
0176   if (values.size() == 0) {
0177     std::ostringstream oss;
0178     oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0179     oss << ": There must at least one value provided.";
0180     throw std::logic_error(oss.str());
0181   }
0182   using Real = typename RandomAccessContainer::value_type;
0183   Real numerator = 0;
0184   Real denominator = 0;
0185   for (size_t i = 0; i < values.size(); ++i) {
0186     if (weights[i] < 0 || !isfinite(weights[i])) {
0187       std::ostringstream oss;
0188       oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0189       oss << ": All weights must be positive and finite, but got received weight " << weights[i] << " at index " << i << " of " << weights.size() << ".";
0190       throw std::domain_error(oss.str());
0191     }
0192     Real tmp = weights[i]*values[i];
0193     numerator += tmp*values[i];
0194     denominator += tmp;
0195   }
0196   return numerator/denominator;
0197 }
0198
0199 } // namespace boost::math::optimization::detail
0200 #endif