From 04ab822bb3bd9a683550ba07aaceb2e59f4a6036 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 10 Feb 2024 12:36:12 -0800 Subject: [PATCH] Units compatibility for optimization --- include/boost/math/optimization/cma_es.hpp | 105 ++++++++++-------- .../boost/math/optimization/detail/common.hpp | 6 +- .../optimization/differential_evolution.hpp | 8 +- include/boost/math/optimization/jso.hpp | 47 ++++---- .../boost/math/optimization/random_search.hpp | 3 +- test/cma_es_test.cpp | 15 +++ test/differential_evolution_test.cpp | 16 +++ test/jso_test.cpp | 16 +++ test/random_search_test.cpp | 19 +++- test/test_functions_for_optimization.hpp | 30 +++++ 10 files changed, 187 insertions(+), 78 deletions(-) diff --git a/include/boost/math/optimization/cma_es.hpp b/include/boost/math/optimization/cma_es.hpp index 910df46916..9470c60a71 100644 --- a/include/boost/math/optimization/cma_es.hpp +++ b/include/boost/math/optimization/cma_es.hpp @@ -42,6 +42,7 @@ namespace boost::math::optimization { template struct cma_es_parameters { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); ArgumentContainer lower_bounds; ArgumentContainer upper_bounds; size_t max_generations = 1000; @@ -51,12 +52,13 @@ template struct cma_es_parameters { // and rounded up to the nearest multiple of threads: size_t population_size = 0; // In the reference, learning_rate = c_m: - Real learning_rate = 1; + DimensionlessReal learning_rate = 1; }; template void validate_cma_es_parameters(cma_es_parameters ¶ms) { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); using std::isfinite; using std::isnan; using std::log; @@ -77,7 +79,7 @@ void validate_cma_es_parameters(cma_es_parameters ¶ms) { //params.population_size = k*params.threads; params.population_size = static_cast(4 + floor(3*log(n))); } - if (params.learning_rate <= Real(0) || !std::isfinite(params.learning_rate)) { + if (params.learning_rate <= DimensionlessReal(0) || !isfinite(params.learning_rate)) { oss << __FILE__ << ":" << __LINE__ << ":" << __func__; oss << ": The learning rate must be > 0, but got " << params.learning_rate << "."; throw std::invalid_argument(oss.str()); @@ -94,8 +96,8 @@ ArgumentContainer cma_es( std::atomic> *current_minimum_cost = nullptr, std::vector>> *queries = nullptr) { - using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); using ResultType = std::invoke_result_t; using std::abs; using std::log; @@ -115,81 +117,81 @@ ArgumentContainer cma_es( std::atomic lowest_cost = std::numeric_limits::infinity(); ArgumentContainer best_vector; // p_{c} := evolution path, equation (24) of the arxiv review: - Eigen::Vector p_c(n); + Eigen::Vector p_c(n); // p_{\sigma} := conjugate evolution path, equation (31) of the arxiv review: - Eigen::Vector p_sigma(n); + Eigen::Vector p_sigma(n); if constexpr (detail::has_resize_v) { best_vector.resize(n, std::numeric_limits::quiet_NaN()); } for (size_t i = 0; i < n; ++i) { - p_c[i] = Real(0); - p_sigma[i] = Real(0); + p_c[i] = DimensionlessReal(0); + p_sigma[i] = DimensionlessReal(0); } // Table 1, \mu = floor(\lambda/2): size_t mu = params.population_size/2; - std::vector w_prime(params.population_size, std::numeric_limits::quiet_NaN()); + std::vector w_prime(params.population_size, std::numeric_limits::quiet_NaN()); for (size_t i = 0; i < params.population_size; ++i) { // Equation (49), but 0-indexed: - w_prime[i] = log(static_cast(params.population_size + 1)/(2*(i+1))); + w_prime[i] = log(static_cast(params.population_size + 1)/(2*(i+1))); } // Table 1, notes at top: - Real positive_weight_sum = 0; - Real sq_weight_sum = 0; + DimensionlessReal positive_weight_sum = 0; + DimensionlessReal sq_weight_sum = 0; for (size_t i = 0; i < mu; ++i) { BOOST_MATH_ASSERT(w_prime[i] > 0); positive_weight_sum += w_prime[i]; sq_weight_sum += w_prime[i]*w_prime[i]; } - Real mu_eff = positive_weight_sum*positive_weight_sum/sq_weight_sum; + DimensionlessReal mu_eff = positive_weight_sum*positive_weight_sum/sq_weight_sum; BOOST_MATH_ASSERT(1 <= mu_eff); BOOST_MATH_ASSERT(mu_eff <= mu); - Real negative_weight_sum = 0; + DimensionlessReal negative_weight_sum = 0; sq_weight_sum = 0; for (size_t i = mu; i < params.population_size; ++i) { BOOST_MATH_ASSERT(w_prime[i] <= 0); negative_weight_sum += w_prime[i]; sq_weight_sum += w_prime[i]*w_prime[i]; } - Real mu_eff_m = negative_weight_sum*negative_weight_sum/sq_weight_sum; + DimensionlessReal mu_eff_m = negative_weight_sum*negative_weight_sum/sq_weight_sum; // Equation (54): - Real c_m = params.learning_rate; + DimensionlessReal c_m = params.learning_rate; // Equation (55): - Real c_sigma = (mu_eff + 2)/(n + mu_eff + 5); + DimensionlessReal c_sigma = (mu_eff + 2)/(n + mu_eff + 5); BOOST_MATH_ASSERT(c_sigma < 1); - Real d_sigma = 1 + 2*(max)(Real(0), sqrt((mu_eff - 1)/(n + 1)) - 1) + c_sigma; + DimensionlessReal d_sigma = 1 + 2*(max)(DimensionlessReal(0), sqrt(DimensionlessReal((mu_eff - 1)/(n + 1))) - DimensionlessReal(1)) + c_sigma; // Equation (56): - Real c_c = (4 + mu_eff/n)/(n + 4 + 2*mu_eff/n); + DimensionlessReal c_c = (4 + mu_eff/n)/(n + 4 + 2*mu_eff/n); BOOST_MATH_ASSERT(c_c <= 1); // Equation (57): - Real c_1 = Real(2)/(pow(n + 1.3, 2) + mu_eff); + DimensionlessReal c_1 = DimensionlessReal(2)/(pow(n + 1.3, 2) + mu_eff); // Equation (58) - Real c_mu = (min)(1 - c_1, 2*(Real(0.25) + mu_eff + 1/mu_eff - 2)/((n+2)*(n+2) + mu_eff)); - BOOST_MATH_ASSERT(c_1 + c_mu <= Real(1)); + DimensionlessReal c_mu = (min)(1 - c_1, 2*(DimensionlessReal(0.25) + mu_eff + 1/mu_eff - 2)/((n+2)*(n+2) + mu_eff)); + BOOST_MATH_ASSERT(c_1 + c_mu <= DimensionlessReal(1)); // Equation (50): - Real alpha_mu_m = 1 + c_1/c_mu; + DimensionlessReal alpha_mu_m = 1 + c_1/c_mu; // Equation (51): - Real alpha_mu_eff_m = 1 + 2*mu_eff_m/(mu_eff + 2); + DimensionlessReal alpha_mu_eff_m = 1 + 2*mu_eff_m/(mu_eff + 2); // Equation (52): - Real alpha_m_pos_def = (1- c_1 - c_mu)/(n*c_mu); + DimensionlessReal alpha_m_pos_def = (1- c_1 - c_mu)/(n*c_mu); // Equation (53): - std::vector weights(params.population_size, std::numeric_limits::quiet_NaN()); + std::vector weights(params.population_size, std::numeric_limits::quiet_NaN()); for (size_t i = 0; i < mu; ++i) { weights[i] = w_prime[i]/positive_weight_sum; } - Real min_alpha = (min)(alpha_mu_m, (min)(alpha_mu_eff_m, alpha_m_pos_def)); + DimensionlessReal min_alpha = (min)(alpha_mu_m, (min)(alpha_mu_eff_m, alpha_m_pos_def)); for (size_t i = mu; i < params.population_size; ++i) { weights[i] = min_alpha*w_prime[i]/abs(negative_weight_sum); } // mu:= number of parents, lambda := number of offspring. - Eigen::Matrix C = Eigen::Matrix::Identity(n, n); + Eigen::Matrix C = Eigen::Matrix::Identity(n, n); ArgumentContainer mean_vector; - Real sigma = 1; + // See the footnote in Figure 6 of the arxiv review: + // We should consider the more robust initialization described there. . . + Real sigma = DimensionlessReal(0.3)*(params.upper_bounds[0] - params.lower_bounds[0]);; if (params.initial_guess) { mean_vector = *params.initial_guess; } else { - // See the footnote in Table 1 of the arxiv review: - sigma = Real(0.3)*(params.upper_bounds[0] - params.lower_bounds[0]); mean_vector = detail::random_initial_population(params.lower_bounds, params.upper_bounds, 1, gen)[0]; } auto initial_cost = cost_function(mean_vector); @@ -221,11 +223,11 @@ ArgumentContainer cma_es( #endif size_t generation = 0; - std::vector> ys(params.population_size); + std::vector> ys(params.population_size); std::vector xs(params.population_size); std::vector costs(params.population_size, std::numeric_limits::quiet_NaN()); - Eigen::Vector weighted_avg_y(n); - Eigen::Vector z(n); + Eigen::Vector weighted_avg_y(n); + Eigen::Vector z(n); if constexpr (detail::has_resize_v) { for (auto & x : xs) { x.resize(n, std::numeric_limits::quiet_NaN()); @@ -234,7 +236,7 @@ ArgumentContainer cma_es( for (auto & y : ys) { y.resize(n); } - normal_distribution dis(Real(0), Real(1)); + normal_distribution dis(DimensionlessReal(0), DimensionlessReal(1)); do { if (cancellation && *cancellation) { break; @@ -243,14 +245,16 @@ ArgumentContainer cma_es( // Section B.2 "Strategy internal numerical effort": // "In practice, the re-calculation of B and D needs to be done not until about // max(1, floor(1/(10n(c_1+c_mu)))) generations." - Eigen::SelfAdjointEigenSolver> eigensolver(C); + // Note that sigma can be dimensionless, in which case C carries the units. + // This is a weird decision-we're not gonna do that! + Eigen::SelfAdjointEigenSolver> eigensolver(C); if (eigensolver.info() != Eigen::Success) { std::ostringstream oss; oss << __FILE__ << ":" << __LINE__ << ":" << __func__; oss << ": Could not decompose the covariance matrix as BDB^{T}."; throw std::logic_error(oss.str()); } - Eigen::Matrix B = eigensolver.eigenvectors(); + Eigen::Matrix B = eigensolver.eigenvectors(); // Eigen returns D^2, in the notation of the survey: auto D = eigensolver.eigenvalues(); // So make it better: @@ -278,7 +282,10 @@ ArgumentContainer cma_es( for (size_t i = 0; i < n; ++i) { z[i] = dis(gen); } - Eigen::Vector Dz = D.array()*z.array(); + Eigen::Vector Dz(n); + for (size_t i = 0; i < n; ++i) { + Dz[i] = D[i]*z[i]; + } y = B*Dz; for (size_t i = 0; i < n; ++i) { BOOST_MATH_ASSERT(!isnan(mean_vector[i])); @@ -331,37 +338,37 @@ ArgumentContainer cma_es( mean_vector[j] = mean_vector[j] + c_m*sigma*weighted_avg_y[j]; } // Equation (43), Figure 6: Start with C^{-1/2}_{w} - Eigen::Vector inv_D_B_transpose_y = B.transpose()*weighted_avg_y; + Eigen::Vector inv_D_B_transpose_y = B.transpose()*weighted_avg_y; for (size_t j = 0; j < inv_D_B_transpose_y.size(); ++j) { inv_D_B_transpose_y[j] /= D[j]; } - Eigen::Vector C_inv_sqrt_y_avg = B*inv_D_B_transpose_y; + Eigen::Vector C_inv_sqrt_y_avg = B*inv_D_B_transpose_y; // Equation (43), Figure 6: - Real p_sigma_norm = 0; + DimensionlessReal p_sigma_norm = 0; for (size_t j = 0; j < n; ++j) { p_sigma[j] = (1-c_sigma)*p_sigma[j] + sqrt(c_sigma*(2-c_sigma)*mu_eff)*C_inv_sqrt_y_avg[j]; p_sigma_norm += p_sigma[j]*p_sigma[j]; } p_sigma_norm = sqrt(p_sigma_norm); // A: Algorithm Summary: E[||N(0,1)||]: - const Real expectation_norm_0I = sqrt(static_cast(n))*(Real(1) - Real(1)/(4*n) + Real(1)/(21*n*n)); + const DimensionlessReal expectation_norm_0I = sqrt(static_cast(n))*(DimensionlessReal(1) - DimensionlessReal(1)/(4*n) + DimensionlessReal(1)/(21*n*n)); // Equation (44), Figure 6: sigma = sigma*exp(c_sigma*(p_sigma_norm/expectation_norm_0I -1)/d_sigma); // A: Algorithm Summary: - Real h_sigma = 0; - Real rhs = (Real(1.4) + Real(2)/(n+1))*expectation_norm_0I*sqrt(1 - pow(1-c_sigma, 2*(generation+1))); + DimensionlessReal h_sigma = 0; + DimensionlessReal rhs = (DimensionlessReal(1.4) + DimensionlessReal(2)/(n+1))*expectation_norm_0I*sqrt(1 - pow(1-c_sigma, 2*(generation+1))); if (p_sigma_norm < rhs) { h_sigma = 1; } // Equation (45), Figure 6: p_c = (1-c_c)*p_c + h_sigma*sqrt(c_c*(2-c_c)*mu_eff)*weighted_avg_y; - Real delta_h_sigma = (1-h_sigma)*c_c*(2-c_c); - Real weight_sum = 0; + DimensionlessReal delta_h_sigma = (1-h_sigma)*c_c*(2-c_c); + DimensionlessReal weight_sum = 0; for (auto & w : weights) { weight_sum += w; } // Equation (47), Figure 6: - Real K = (1 + c_1*delta_h_sigma - c_1 - c_mu*weight_sum); + DimensionlessReal K = (1 + c_1*delta_h_sigma - c_1 - c_mu*weight_sum); // Can these operations be sped up using `.selfadjointView`? // Maybe: A.selfadjointView().rankUpdate(p_c, c_1);? C = K*C + c_1*p_c*p_c.transpose(); @@ -370,12 +377,12 @@ ArgumentContainer cma_es( C += c_mu*weights[i]*ys[indices[i]]*ys[indices[i]].transpose(); } for (size_t i = params.population_size/2; i < params.population_size; ++i) { - Eigen::Vector D_inv_BTy = B.transpose()*ys[indices[i]]; + Eigen::Vector D_inv_BTy = B.transpose()*ys[indices[i]]; for (size_t j = 0; j < n; ++j) { D_inv_BTy[j] /= D[j]; } - Real squared_norm = D_inv_BTy.squaredNorm(); - Real K2 = c_mu*weights[i]/squared_norm; + DimensionlessReal squared_norm = D_inv_BTy.squaredNorm(); + DimensionlessReal K2 = c_mu*weights[i]/squared_norm; C += K2*ys[indices[i]]*ys[indices[i]].transpose(); } } while (generation++ < params.max_generations); diff --git a/include/boost/math/optimization/detail/common.hpp b/include/boost/math/optimization/detail/common.hpp index 5853c26f8f..f98fbd62bd 100644 --- a/include/boost/math/optimization/detail/common.hpp +++ b/include/boost/math/optimization/detail/common.hpp @@ -68,6 +68,7 @@ std::vector random_initial_population(ArgumentContainer const ArgumentContainer const &upper_bounds, size_t initial_population_size, URBG &&gen) { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); constexpr bool has_resize = detail::has_resize_v; std::vector population(initial_population_size); auto const dimension = lower_bounds.size(); @@ -98,7 +99,7 @@ std::vector random_initial_population(ArgumentContainer const // That said, scipy uses Latin Hypercube sampling and says self-avoiding sequences are preferable. // So this is something that could be investigated and potentially improved. using std::uniform_real_distribution; - uniform_real_distribution dis(Real(0), Real(1)); + uniform_real_distribution dis(DimensionlessReal(0), DimensionlessReal(1)); for (size_t i = 0; i < population.size(); ++i) { for (size_t j = 0; j < dimension; ++j) { auto const &lb = lower_bounds[j]; @@ -163,6 +164,7 @@ template std::vector best_indices(std::vector cons template auto weighted_lehmer_mean(RandomAccessContainer const & values, RandomAccessContainer const & weights) { + using std::isfinite; if (values.size() != weights.size()) { std::ostringstream oss; oss << __FILE__ << ":" << __LINE__ << ":" << __func__; @@ -179,7 +181,7 @@ auto weighted_lehmer_mean(RandomAccessContainer const & values, RandomAccessCont Real numerator = 0; Real denominator = 0; for (size_t i = 0; i < values.size(); ++i) { - if (weights[i] < 0 || !std::isfinite(weights[i])) { + if (weights[i] < 0 || !isfinite(weights[i])) { std::ostringstream oss; oss << __FILE__ << ":" << __LINE__ << ":" << __func__; oss << ": All weights must be positive and finite, but got received weight " << weights[i] << " at index " << i << " of " << weights.size() << "."; diff --git a/include/boost/math/optimization/differential_evolution.hpp b/include/boost/math/optimization/differential_evolution.hpp index 36ccfa781d..6298f3bf97 100644 --- a/include/boost/math/optimization/differential_evolution.hpp +++ b/include/boost/math/optimization/differential_evolution.hpp @@ -35,11 +35,12 @@ namespace boost::math::optimization { // We provide the parameters in a struct-there are too many of them and they are too unwieldy to pass individually: template struct differential_evolution_parameters { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); ArgumentContainer lower_bounds; ArgumentContainer upper_bounds; // mutation factor is also called scale factor or just F in the literature: - Real mutation_factor = static_cast(0.65); - double crossover_probability = 0.5; + DimensionlessReal mutation_factor = static_cast(0.65); + DimensionlessReal crossover_probability = static_cast(0.5); // Population in each generation: size_t NP = 500; size_t max_generations = 1000; @@ -94,6 +95,7 @@ ArgumentContainer differential_evolution( std::vector>> *queries = nullptr, std::atomic> *current_minimum_cost = nullptr) { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); using ResultType = std::invoke_result_t; using std::clamp; using std::isnan; @@ -159,7 +161,7 @@ ArgumentContainer differential_evolution( // I first tried seeding thread-local random number generators with the global generator, // but ThreadSanitizer didn't like it. I *suspect* there's a way around this, but // even if it's formally threadsafe, there's still a lot of effort to get it computationally reproducible. - uniform_real_distribution unif01(Real(0), Real(1)); + uniform_real_distribution unif01(DimensionlessReal(0), DimensionlessReal(1)); for (size_t i = 0; i < NP; ++i) { size_t r1, r2, r3; do { diff --git a/include/boost/math/optimization/jso.hpp b/include/boost/math/optimization/jso.hpp index ed31a5235b..e48d93836d 100644 --- a/include/boost/math/optimization/jso.hpp +++ b/include/boost/math/optimization/jso.hpp @@ -36,6 +36,7 @@ namespace boost::math::optimization { // IEEE Transactions on evolutionary computation, 13(5), 945-958." template struct jso_parameters { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); ArgumentContainer lower_bounds; ArgumentContainer upper_bounds; // Population in the first generation. @@ -107,6 +108,7 @@ jso(const Func cost_function, jso_parameters &jso_params, std::invoke_result_t>> *queries = nullptr) { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); validate_jso_parameters(jso_params); using ResultType = std::invoke_result_t; @@ -116,6 +118,7 @@ jso(const Func cost_function, jso_parameters &jso_params, using std::isnan; using std::max; using std::round; + using std::isfinite; using std::uniform_real_distribution; // Refer to the referenced paper, pseudocode on page 1313: @@ -173,12 +176,12 @@ jso(const Func cost_function, jso_parameters &jso_params, // last bullet, which claims this should be set to 0.3. The reference // implementation also does 0.3: size_t H = 5; - std::vector M_F(H, static_cast(0.3)); + std::vector M_F(H, static_cast(0.3)); // Algorithm 1: jSO algorithm, Line 4: // "Set all values in M_CR to 0.8": - std::vector M_CR(H, static_cast(0.8)); + std::vector M_CR(H, static_cast(0.8)); - std::uniform_real_distribution unif01(Real(0), Real(1)); + std::uniform_real_distribution unif01(DimensionlessReal(0), DimensionlessReal(1)); bool keep_going = !target_attained; if (cancellation && *cancellation) { keep_going = false; @@ -190,8 +193,8 @@ jso(const Func cost_function, jso_parameters &jso_params, size_t minimum_population_size = (max)(size_t(4), size_t(jso_params.threads)); while (keep_going) { // Algorithm 1, jSO, Line 7: - std::vector S_CR; - std::vector S_F; + std::vector S_CR; + std::vector S_F; // Equation 9 of L-SHADE: std::vector delta_f; for (size_t i = 0; i < population.size(); ++i) { @@ -203,20 +206,20 @@ jso(const Func cost_function, jso_parameters &jso_params, // I confess I find it weird to store the historical memory if we're just // gonna ignore it, but that's what the paper and the reference // implementation says! - Real mu_F = static_cast(0.9); - Real mu_CR = static_cast(0.9); + DimensionlessReal mu_F = static_cast(0.9); + DimensionlessReal mu_CR = static_cast(0.9); if (ri != H - 1) { mu_F = M_F[ri]; mu_CR = M_CR[ri]; } // Algorithm 1, jSO, Line 14-18: - Real crossover_probability = static_cast(0); + DimensionlessReal crossover_probability = static_cast(0); if (mu_CR >= 0) { using std::normal_distribution; - normal_distribution normal(mu_CR, static_cast(0.1)); + normal_distribution normal(mu_CR, static_cast(0.1)); crossover_probability = normal(gen); // Clamp comes from L-SHADE description: - crossover_probability = clamp(crossover_probability, Real(0), Real(1)); + crossover_probability = clamp(crossover_probability, DimensionlessReal(0), DimensionlessReal(1)); } // Algorithm 1, jSO, Line 19-23: // Note that the pseudocode uses a "max_generations parameter", @@ -224,24 +227,24 @@ jso(const Func cost_function, jso_parameters &jso_params, // Since we already require specification of max_function_evaluations, // the pseudocode adds an unnecessary parameter. if (4 * function_evaluations < jso_params.max_function_evaluations) { - crossover_probability = (max)(crossover_probability, Real(0.7)); + crossover_probability = (max)(crossover_probability, DimensionlessReal(0.7)); } else if (2 * function_evaluations < jso_params.max_function_evaluations) { - crossover_probability = (max)(crossover_probability, Real(0.6)); + crossover_probability = (max)(crossover_probability, DimensionlessReal(0.6)); } // Algorithm 1, jSO, Line 24-27: // Note the adjustments to the pseudocode given in the reference // implementation. - cauchy_distribution cauchy(mu_F, static_cast(0.1)); - Real F; + cauchy_distribution cauchy(mu_F, static_cast(0.1)); + DimensionlessReal F; do { F = cauchy(gen); if (F > 1) { F = 1; } } while (F <= 0); - Real threshold = static_cast(7) / static_cast(10); + DimensionlessReal threshold = static_cast(7) / static_cast(10); if ((10 * function_evaluations < 6 * jso_params.max_function_evaluations) && (F > threshold)) { @@ -250,16 +253,16 @@ jso(const Func cost_function, jso_parameters &jso_params, // > p value for mutation strategy linearly decreases from pmax to pmin // during the evolutionary process, > where pmax = 0.25 in jSO and pmin = // pmax/2. - Real p = Real(0.25) * (1 - static_cast(function_evaluations) / + DimensionlessReal p = DimensionlessReal(0.25) * (1 - static_cast(function_evaluations) / (2 * jso_params.max_function_evaluations)); // Equation (4) of the reference: - Real Fw = static_cast(1.2) * F; + DimensionlessReal Fw = static_cast(1.2) * F; if (10 * function_evaluations < 4 * jso_params.max_function_evaluations) { if (10 * function_evaluations < 2 * jso_params.max_function_evaluations) { - Fw = static_cast(0.7) * F; + Fw = static_cast(0.7) * F; } else { - Fw = static_cast(0.8) * F; + Fw = static_cast(0.8) * F; } } // Algorithm 1, jSO, Line 28: @@ -367,13 +370,13 @@ jso(const Func cost_function, jso_parameters &jso_params, // If there are no successful updates this generation, we do not update the // historical memory: if (S_CR.size() > 0) { - std::vector weights(S_CR.size(), - std::numeric_limits::quiet_NaN()); + std::vector weights(S_CR.size(), + std::numeric_limits::quiet_NaN()); ResultType delta_sum = static_cast(0); for (auto const &delta : delta_f) { delta_sum += delta; } - if (delta_sum <= 0 || !std::isfinite(delta_sum)) { + if (delta_sum <= 0 || !isfinite(delta_sum)) { std::ostringstream oss; oss << __FILE__ << ":" << __LINE__ << ":" << __func__; oss << "\n\tYou've hit a bug: The sum of improvements must be strictly " diff --git a/include/boost/math/optimization/random_search.hpp b/include/boost/math/optimization/random_search.hpp index b1637b72a4..b419c24064 100644 --- a/include/boost/math/optimization/random_search.hpp +++ b/include/boost/math/optimization/random_search.hpp @@ -60,6 +60,7 @@ ArgumentContainer random_search( std::vector>> *queries = nullptr) { using Real = typename ArgumentContainer::value_type; + using DimensionlessReal = decltype(Real()/Real()); using ResultType = std::invoke_result_t; using std::isnan; using std::uniform_real_distribution; @@ -106,7 +107,7 @@ ArgumentContainer random_search( break; } // Fill trial vector: - uniform_real_distribution unif01(Real(0), Real(1)); + uniform_real_distribution unif01(DimensionlessReal(0), DimensionlessReal(1)); for (size_t i = 0; i < dimension; ++i) { trial_vector[i] = params.lower_bounds[i] + (params.upper_bounds[i] - params.lower_bounds[i])*unif01(g); } diff --git a/test/cma_es_test.cpp b/test/cma_es_test.cpp index f221094223..d2016997e6 100644 --- a/test/cma_es_test.cpp +++ b/test/cma_es_test.cpp @@ -146,6 +146,18 @@ void test_beale() { CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(0.1)); } +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY +void test_dimensioned_sphere() { + std::cout << "Testing CMA-ES on dimensioned sphere . . .\n"; + using ArgType = std::vector>; + auto params = cma_es_parameters(); + params.lower_bounds.resize(4, -1.0*meter); + params.upper_bounds.resize(4, 1*meter); + std::mt19937_64 gen(56789); + auto local_minima = cma_es(dimensioned_sphere, params, gen); +} +#endif + int main() { #if (defined(__clang__) || defined(_MSC_VER)) test_ackley(); @@ -154,6 +166,9 @@ int main() { test_rastrigin(); test_three_hump_camel(); test_beale(); +#endif +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY + test_dimensioned_sphere(); #endif test_sphere(); return boost::math::test::report_errors(); diff --git a/test/differential_evolution_test.cpp b/test/differential_evolution_test.cpp index fdbb97dfad..b3d70c901f 100644 --- a/test/differential_evolution_test.cpp +++ b/test/differential_evolution_test.cpp @@ -188,6 +188,19 @@ void test_beale() { CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(2e-4)); } +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY +void test_dimensioned_sphere() { + std::cout << "Testing differential evolution on dimensioned sphere . . .\n"; + using ArgType = std::vector>; + auto params = differential_evolution_parameters(); + params.lower_bounds.resize(4, -1.0*meter); + params.upper_bounds.resize(4, 1*meter); + params.threads = 2; + std::mt19937_64 gen(56789); + auto local_minima = differential_evolution(dimensioned_sphere, params, gen); +} +#endif + int main() { #if defined(__clang__) || defined(_MSC_VER) @@ -197,6 +210,9 @@ int main() { test_rastrigin(); test_three_hump_camel(); test_beale(); +#endif +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY + test_dimensioned_sphere(); #endif test_sphere(); test_parameter_checks(); diff --git a/test/jso_test.cpp b/test/jso_test.cpp index 8085ec95cf..443257a6c0 100644 --- a/test/jso_test.cpp +++ b/test/jso_test.cpp @@ -143,6 +143,19 @@ void test_beale() { CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(2e-4)); } +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY +void test_dimensioned_sphere() { + std::cout << "Testing jso on dimensioned sphere . . .\n"; + using ArgType = std::vector>; + auto params = jso_parameters(); + params.lower_bounds.resize(4, -1.0*meter); + params.upper_bounds.resize(4, 1*meter); + params.threads = 2; + std::mt19937_64 gen(56789); + auto local_minima = jso(dimensioned_sphere, params, gen); +} +#endif + int main() { #if defined(__clang__) || defined(_MSC_VER) test_ackley(); @@ -151,6 +164,9 @@ int main() { test_rastrigin(); test_three_hump_camel(); test_beale(); +#endif +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY + test_dimensioned_sphere(); #endif test_sphere(); test_weighted_lehmer_mean(); diff --git a/test/random_search_test.cpp b/test/random_search_test.cpp index f671b93707..acb690abbd 100644 --- a/test/random_search_test.cpp +++ b/test/random_search_test.cpp @@ -10,7 +10,6 @@ #include #include #include - using boost::math::optimization::random_search; using boost::math::optimization::random_search_parameters; @@ -149,6 +148,21 @@ void test_beale() { CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(0.1)); } +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY +void test_dimensioned_sphere() { + std::cout << "Testing random search on dimensioned sphere . . .\n"; + using ArgType = std::vector>; + auto rs_params = random_search_parameters(); + rs_params.lower_bounds.resize(4, -1.0*meter); + rs_params.upper_bounds.resize(4, 1*meter); + rs_params.max_function_calls = 100000; + rs_params.threads = 2; + std::mt19937_64 gen(56789); + auto local_minima = random_search(dimensioned_sphere, rs_params, gen); +} + +#endif + int main() { #if defined(__clang__) || defined(_MSC_VER) test_ackley(); @@ -157,6 +171,9 @@ int main() { test_rastrigin(); test_three_hump_camel(); test_beale(); +#endif +#if BOOST_MATH_TEST_UNITS_COMPATIBILITY + test_dimensioned_sphere(); #endif test_sphere(); return boost::math::test::report_errors(); diff --git a/test/test_functions_for_optimization.hpp b/test/test_functions_for_optimization.hpp index 4b8567b25d..513b97d8d0 100644 --- a/test/test_functions_for_optimization.hpp +++ b/test/test_functions_for_optimization.hpp @@ -6,7 +6,37 @@ */ #ifndef TEST_FUNCTIONS_FOR_OPTIMIZATION_HPP #define TEST_FUNCTIONS_FOR_OPTIMIZATION_HPP +#include +#include #include +#if __has_include() +// This is the only system boost.units still works on. +// I imagine this will start to fail at some point, +// and we'll have to remove this test as well. +#if defined(__APPLE__) +#define BOOST_MATH_TEST_UNITS_COMPATIBILITY 1 +#include +#include +#include +#include +#include +using namespace boost::units; +using namespace boost::units::si; + +// This *should* return an area, but see: https://github.com/boostorg/units/issues/58 +// This sadly prevents std::atomic>. +// Nonetheless, we *do* get some information making the argument type dimensioned, +// even if it would be better to get the full information: +double dimensioned_sphere(std::vector> const & v) { + quantity r(0.0*meters*meters); + for (auto const & x : v) { + r += (x * x); + } + quantity scale(1.0*meters*meters); + return static_cast(r/scale); +} +#endif +#endif // Taken from: https://en.wikipedia.org/wiki/Test_functions_for_optimization template Real ackley(std::array const &v) {