-
Notifications
You must be signed in to change notification settings - Fork 226
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Fix race conditions in differential evolution
Through a combination of silly mistakes, I missed a pile of race conditions in the OpenMP threading. Switch to C++ threading. Note that this change requires serial generation of trial vectors. Hopefully I can figure out to parallelize the generation of trial vectors to reduce the serial section a la Ahmdahl's law, while simultaneously keeping thread sanitizer happy.
- Loading branch information
1 parent
79b4015
commit e655c8a
Showing
12 changed files
with
557 additions
and
395 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -27,3 +27,4 @@ cmake-build-debug/* | |
| .cmake/* | ||
| build.ninja | ||
| .ninja* | ||
| a.out | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,6 @@ | ||
| #include <iostream> | ||
| #include <boost/math/optimization/jso.hpp> | ||
|
|
||
| int main() { | ||
| std::cout << "Yo!\n"; | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,165 @@ | ||
| /* | ||
| * Copyright Nick Thompson, 2024 | ||
| * Use, modification and distribution are subject to the | ||
| * Boost Software License, Version 1.0. (See accompanying file | ||
| * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | ||
| */ | ||
| #ifndef BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP | ||
| #define BOOST_MATH_OPTIMIZATION_DETAIL_COMMON_HPP | ||
| #include <algorithm> | ||
| #include <cmath> | ||
| #include <limits> | ||
| #include <sstream> | ||
| #include <stdexcept> | ||
| #include <random> | ||
|
|
||
| namespace boost::math::optimization::detail { | ||
|
|
||
| template <typename T, typename = void> struct has_resize : std::false_type {}; | ||
|
|
||
| template <typename T> | ||
| struct has_resize<T, std::void_t<decltype(std::declval<T>().resize(size_t{}))>> : std::true_type {}; | ||
|
|
||
| template <typename T> constexpr bool has_resize_v = has_resize<T>::value; | ||
|
|
||
| template <typename ArgumentContainer> | ||
| void validate_bounds(ArgumentContainer const &lower_bounds, ArgumentContainer const &upper_bounds) { | ||
| using std::isfinite; | ||
| std::ostringstream oss; | ||
| if (lower_bounds.size() == 0) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": The dimension of the problem cannot be zero."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| if (upper_bounds.size() != lower_bounds.size()) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": There must be the same number of lower bounds as upper bounds, but given "; | ||
| oss << upper_bounds.size() << " upper bounds, and " << lower_bounds.size() << " lower bounds."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| for (size_t i = 0; i < lower_bounds.size(); ++i) { | ||
| auto lb = lower_bounds[i]; | ||
| auto ub = upper_bounds[i]; | ||
| if (lb > ub) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": The upper bound must be greater than or equal to the lower bound, but the upper bound is " << ub | ||
| << " and the lower is " << lb << "."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| if (!isfinite(lb)) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": The lower bound must be finite, but got " << lb << "."; | ||
| oss << " For infinite bounds, emulate with std::numeric_limits<Real>::lower() or use a standard infinite->finite " | ||
| "transform."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| if (!isfinite(ub)) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": The upper bound must be finite, but got " << ub << "."; | ||
| oss << " For infinite bounds, emulate with std::numeric_limits<Real>::max() or use a standard infinite->finite " | ||
| "transform."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| template <typename ArgumentContainer, class URBG> | ||
| std::vector<ArgumentContainer> random_initial_population(ArgumentContainer const &lower_bounds, | ||
| ArgumentContainer const &upper_bounds, | ||
| size_t initial_population_size, URBG &&gen) { | ||
| using Real = typename ArgumentContainer::value_type; | ||
| constexpr bool has_resize = detail::has_resize_v<ArgumentContainer>; | ||
| std::vector<ArgumentContainer> population(initial_population_size); | ||
| auto const dimension = lower_bounds.size(); | ||
| for (size_t i = 0; i < population.size(); ++i) { | ||
| if constexpr (has_resize) { | ||
| population[i].resize(dimension); | ||
| } else { | ||
| // Argument type must be known at compile-time; like std::array: | ||
| if (population[i].size() != dimension) { | ||
| std::ostringstream oss; | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": For containers which do not have resize, the default size must be the same as the dimension, "; | ||
| oss << "but the default container size is " << population[i].size() << " and the dimension of the problem is " | ||
| << dimension << "."; | ||
| oss << " The function argument container type is " << typeid(ArgumentContainer).name() << ".\n"; | ||
| throw std::runtime_error(oss.str()); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| // Why don't we provide an option to initialize with (say) a Gaussian distribution? | ||
| // > If the optimum's location is fairly well known, | ||
| // > a Gaussian distribution may prove somewhat faster, although it | ||
| // > may also increase the probability that the population will converge prematurely. | ||
| // > In general, uniform distributions are preferred, since they best reflect | ||
| // > the lack of knowledge about the optimum's location. | ||
| // - Differential Evolution: A Practical Approach to Global Optimization | ||
| // 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<Real> dis(Real(0), Real(1)); | ||
| for (size_t i = 0; i < population.size(); ++i) { | ||
| for (size_t j = 0; j < dimension; ++j) { | ||
| auto const &lb = lower_bounds[j]; | ||
| auto const &ub = upper_bounds[j]; | ||
| population[i][j] = lb + dis(gen) * (ub - lb); | ||
| } | ||
| } | ||
|
|
||
| return population; | ||
| } | ||
|
|
||
| template <typename ArgumentContainer> | ||
| void validate_initial_guess(ArgumentContainer const &initial_guess, ArgumentContainer const &lower_bounds, | ||
| ArgumentContainer const &upper_bounds) { | ||
| std::ostringstream oss; | ||
| auto const dimension = lower_bounds.size(); | ||
| if (initial_guess.size() != dimension) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": The initial guess must have the same dimensions as the problem,"; | ||
| oss << ", but the problem size is " << dimension << " and the initial guess has " << initial_guess.size() | ||
| << " elements."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| for (size_t i = 0; i < dimension; ++i) { | ||
| auto lb = lower_bounds[i]; | ||
| auto ub = upper_bounds[i]; | ||
| if (!isfinite(initial_guess[i])) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": At index " << i << ", the initial guess is " << initial_guess[i] | ||
| << ", make sure all elements of the initial guess are finite."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| if (initial_guess[i] < lb || initial_guess[i] > ub) { | ||
| oss << __FILE__ << ":" << __LINE__ << ":" << __func__; | ||
| oss << ": At index " << i << " the initial guess " << initial_guess[i] << " is not in the bounds [" << lb << ", " | ||
| << ub << "]."; | ||
| throw std::domain_error(oss.str()); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| // Return indices corresponding to the minimum function values. | ||
| template <typename Real> std::vector<size_t> best_indices(std::vector<Real> const &function_values) { | ||
| using std::isnan; | ||
| const size_t n = function_values.size(); | ||
| std::vector<size_t> indices(n); | ||
| for (size_t i = 0; i < n; ++i) { | ||
| indices[i] = i; | ||
| } | ||
|
|
||
| std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) { | ||
| if (isnan(function_values[a])) { | ||
| return false; | ||
| } | ||
| if (isnan(function_values[b])) { | ||
| return true; | ||
| } | ||
| return function_values[a] < function_values[b]; | ||
| }); | ||
| return indices; | ||
| } | ||
|
|
||
| } // namespace boost::math::optimization::detail | ||
| #endif |
Oops, something went wrong.