|
| 1 | +#include "kondakov_v_global_search/mpi/include/ops_mpi.hpp" |
| 2 | + |
| 3 | +#include <mpi.h> |
| 4 | + |
| 5 | +#include <algorithm> |
| 6 | +#include <array> |
| 7 | +#include <cmath> |
| 8 | +#include <cstddef> |
| 9 | +#include <iterator> |
| 10 | +#include <limits> |
| 11 | +#include <utility> |
| 12 | +#include <vector> |
| 13 | + |
| 14 | +#include "kondakov_v_global_search/common/include/common.hpp" |
| 15 | + |
| 16 | +namespace kondakov_v_global_search { |
| 17 | + |
| 18 | +KondakovVGlobalSearchMPI::KondakovVGlobalSearchMPI(const InType &in) { |
| 19 | + MPI_Comm_rank(MPI_COMM_WORLD, &world_rank_); |
| 20 | + MPI_Comm_size(MPI_COMM_WORLD, &world_size_); |
| 21 | + SetTypeOfTask(GetStaticTypeOfTask()); |
| 22 | + GetInput() = in; |
| 23 | + GetOutput() = {}; |
| 24 | +} |
| 25 | + |
| 26 | +double KondakovVGlobalSearchMPI::EvaluateFunction(double x) { |
| 27 | + const auto &cfg = GetInput(); |
| 28 | + switch (cfg.func_type) { |
| 29 | + case FunctionType::kQuadratic: { |
| 30 | + double t = cfg.func_param; |
| 31 | + return (x - t) * (x - t); |
| 32 | + } |
| 33 | + case FunctionType::kSine: |
| 34 | + return std::sin(x) + (0.1 * x); |
| 35 | + case FunctionType::kAbs: |
| 36 | + return std::abs(x); |
| 37 | + default: |
| 38 | + return std::numeric_limits<double>::quiet_NaN(); |
| 39 | + } |
| 40 | +} |
| 41 | + |
| 42 | +bool KondakovVGlobalSearchMPI::IsRoot() const { |
| 43 | + return world_rank_ == 0; |
| 44 | +} |
| 45 | + |
| 46 | +bool KondakovVGlobalSearchMPI::ValidationImpl() { |
| 47 | + const auto &cfg = GetInput(); |
| 48 | + bool local_valid = cfg.left < cfg.right && cfg.accuracy > 0.0 && cfg.reliability > 0.0 && cfg.max_iterations > 0; |
| 49 | + |
| 50 | + bool global_valid = false; |
| 51 | + MPI_Allreduce(&local_valid, &global_valid, 1, MPI_C_BOOL, MPI_LAND, MPI_COMM_WORLD); |
| 52 | + return global_valid; |
| 53 | +} |
| 54 | + |
| 55 | +bool KondakovVGlobalSearchMPI::PreProcessingImpl() { |
| 56 | + const auto &cfg = GetInput(); |
| 57 | + |
| 58 | + if (IsRoot()) { |
| 59 | + points_x_.clear(); |
| 60 | + values_y_.clear(); |
| 61 | + points_x_.reserve(cfg.max_iterations + (world_size_ * 10)); |
| 62 | + values_y_.reserve(cfg.max_iterations + (world_size_ * 10)); |
| 63 | + |
| 64 | + double f_a = EvaluateFunction(cfg.left); |
| 65 | + double f_b = EvaluateFunction(cfg.right); |
| 66 | + if (!std::isfinite(f_a) || !std::isfinite(f_b)) { |
| 67 | + return false; |
| 68 | + } |
| 69 | + |
| 70 | + points_x_ = {cfg.left, cfg.right}; |
| 71 | + values_y_ = {f_a, f_b}; |
| 72 | + best_point_ = (f_a < f_b) ? cfg.left : cfg.right; |
| 73 | + best_value_ = std::min(f_a, f_b); |
| 74 | + } |
| 75 | + |
| 76 | + SyncGlobalState(); |
| 77 | + return true; |
| 78 | +} |
| 79 | + |
| 80 | +void KondakovVGlobalSearchMPI::SyncGlobalState() { |
| 81 | + int n = IsRoot() ? static_cast<int>(points_x_.size()) : 0; |
| 82 | + MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD); |
| 83 | + |
| 84 | + if (!IsRoot()) { |
| 85 | + points_x_.resize(n); |
| 86 | + values_y_.resize(n); |
| 87 | + } |
| 88 | + |
| 89 | + if (n > 0) { |
| 90 | + MPI_Bcast(points_x_.data(), n, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
| 91 | + MPI_Bcast(values_y_.data(), n, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
| 92 | + } |
| 93 | + |
| 94 | + best_value_ = std::numeric_limits<double>::max(); |
| 95 | + for (std::size_t i = 0; i < points_x_.size(); ++i) { |
| 96 | + if (values_y_[i] < best_value_) { |
| 97 | + best_value_ = values_y_[i]; |
| 98 | + best_point_ = points_x_[i]; |
| 99 | + } |
| 100 | + } |
| 101 | +} |
| 102 | + |
| 103 | +double KondakovVGlobalSearchMPI::ComputeAdaptiveLipschitzEstimate(double r) const { |
| 104 | + const double min_slope = 1e-2; |
| 105 | + if (points_x_.size() < 2) { |
| 106 | + return r * min_slope; |
| 107 | + } |
| 108 | + |
| 109 | + double max_slope = min_slope; |
| 110 | + for (std::size_t i = 1; i < points_x_.size(); ++i) { |
| 111 | + double dx = points_x_[i] - points_x_[i - 1]; |
| 112 | + if (dx <= 0.0) { |
| 113 | + continue; |
| 114 | + } |
| 115 | + double dy = std::abs(values_y_[i] - values_y_[i - 1]); |
| 116 | + if (!std::isfinite(dy)) { |
| 117 | + continue; |
| 118 | + } |
| 119 | + double slope = dy / dx; |
| 120 | + if (std::isfinite(slope) && slope > max_slope) { |
| 121 | + max_slope = slope; |
| 122 | + } |
| 123 | + } |
| 124 | + return r * max_slope; |
| 125 | +} |
| 126 | + |
| 127 | +double KondakovVGlobalSearchMPI::IntervalMerit(std::size_t i, double l_est) const { |
| 128 | + double x_l = points_x_[i - 1]; |
| 129 | + double x_r = points_x_[i]; |
| 130 | + double f_l = values_y_[i - 1]; |
| 131 | + double f_r = values_y_[i]; |
| 132 | + double h = x_r - x_l; |
| 133 | + double df = f_r - f_l; |
| 134 | + return (l_est * h) - (2.0 * (f_l + f_r)) + ((df * df) / (l_est * h)); |
| 135 | +} |
| 136 | + |
| 137 | +double KondakovVGlobalSearchMPI::ProposeTrialPoint(std::size_t i, double l_est) const { |
| 138 | + double x_l = points_x_[i - 1]; |
| 139 | + double x_r = points_x_[i]; |
| 140 | + double f_l = values_y_[i - 1]; |
| 141 | + double f_r = values_y_[i]; |
| 142 | + double mid = 0.5 * (x_l + x_r); |
| 143 | + double asym = (f_r - f_l) / (2.0 * l_est); |
| 144 | + double cand = mid - asym; |
| 145 | + if (cand <= x_l || cand >= x_r) { |
| 146 | + cand = mid; |
| 147 | + } |
| 148 | + return cand; |
| 149 | +} |
| 150 | + |
| 151 | +std::size_t KondakovVGlobalSearchMPI::LocateInsertionIndex(double x) const { |
| 152 | + auto it = std::ranges::lower_bound(points_x_, x); |
| 153 | + return static_cast<std::size_t>(std::distance(points_x_.begin(), it)); |
| 154 | +} |
| 155 | + |
| 156 | +void KondakovVGlobalSearchMPI::InsertEvaluation(double x, double fx) { |
| 157 | + auto idx = LocateInsertionIndex(x); |
| 158 | + points_x_.insert(points_x_.begin() + static_cast<std::vector<double>::difference_type>(idx), x); |
| 159 | + values_y_.insert(values_y_.begin() + static_cast<std::vector<double>::difference_type>(idx), fx); |
| 160 | + if (fx < best_value_) { |
| 161 | + best_value_ = fx; |
| 162 | + best_point_ = x; |
| 163 | + } |
| 164 | +} |
| 165 | + |
| 166 | +void KondakovVGlobalSearchMPI::SelectIntervalsToRefine(double l_est, |
| 167 | + std::vector<std::pair<double, std::size_t>> &merits) { |
| 168 | + merits.clear(); |
| 169 | + for (std::size_t i = 1; i < points_x_.size(); ++i) { |
| 170 | + merits.emplace_back(IntervalMerit(i, l_est), i); |
| 171 | + } |
| 172 | + std::ranges::sort(merits, [](const auto &a, const auto &b) { return a.first > b.first; }); |
| 173 | +} |
| 174 | + |
| 175 | +bool KondakovVGlobalSearchMPI::CheckConvergence(const Params &cfg, |
| 176 | + const std::vector<std::pair<double, std::size_t>> &merits) { |
| 177 | + if (merits.empty()) { |
| 178 | + return false; |
| 179 | + } |
| 180 | + double width = points_x_[merits[0].second] - points_x_[merits[0].second - 1]; |
| 181 | + if (width <= cfg.accuracy) { |
| 182 | + has_converged_ = true; |
| 183 | + return true; |
| 184 | + } |
| 185 | + return false; |
| 186 | +} |
| 187 | + |
| 188 | +void KondakovVGlobalSearchMPI::GatherAndBroadcastTrialResults(const std::vector<std::pair<double, std::size_t>> &merits, |
| 189 | + int num_trials, double l_est) { |
| 190 | + double local_x = 0.0; |
| 191 | + double local_fx = 0.0; |
| 192 | + int local_count = 0; |
| 193 | + |
| 194 | + if (world_rank_ < num_trials && !merits.empty()) { |
| 195 | + std::size_t idx = merits[world_rank_].second; |
| 196 | + double x = ProposeTrialPoint(idx, l_est); |
| 197 | + double fx = EvaluateFunction(x); |
| 198 | + if (std::isfinite(fx)) { |
| 199 | + local_x = x; |
| 200 | + local_fx = fx; |
| 201 | + local_count = 2; |
| 202 | + } |
| 203 | + } |
| 204 | + |
| 205 | + std::vector<int> counts(world_size_); |
| 206 | + MPI_Allgather(&local_count, 1, MPI_INT, counts.data(), 1, MPI_INT, MPI_COMM_WORLD); |
| 207 | + |
| 208 | + std::vector<int> displs(world_size_); |
| 209 | + for (int i = 1; i < world_size_; ++i) { |
| 210 | + displs[i] = displs[i - 1] + counts[i - 1]; |
| 211 | + } |
| 212 | + int total = (world_size_ > 0) ? (displs.back() + counts.back()) : 0; |
| 213 | + |
| 214 | + std::vector<double> recv_buf; |
| 215 | + if (total > 0) { |
| 216 | + recv_buf.resize(total); |
| 217 | + } |
| 218 | + std::array<double, 2> send_buf = {local_x, local_fx}; |
| 219 | + const double *send_ptr = (local_count > 0) ? send_buf.data() : nullptr; |
| 220 | + |
| 221 | + MPI_Allgatherv(send_ptr, local_count, MPI_DOUBLE, recv_buf.data(), counts.data(), displs.data(), MPI_DOUBLE, |
| 222 | + MPI_COMM_WORLD); |
| 223 | + |
| 224 | + if (IsRoot()) { |
| 225 | + for (int i = 0; i < total; i += 2) { |
| 226 | + InsertEvaluation(recv_buf[i], recv_buf[i + 1]); |
| 227 | + } |
| 228 | + } |
| 229 | +} |
| 230 | + |
| 231 | +bool KondakovVGlobalSearchMPI::RunImpl() { |
| 232 | + const auto &cfg = GetInput(); |
| 233 | + std::vector<std::pair<double, std::size_t>> merits; |
| 234 | + double l_est = ComputeAdaptiveLipschitzEstimate(cfg.reliability); |
| 235 | + |
| 236 | + for (int step = 0; step < cfg.max_iterations; ++step) { |
| 237 | + if (step % 10 == 0) { |
| 238 | + l_est = ComputeAdaptiveLipschitzEstimate(cfg.reliability); |
| 239 | + } |
| 240 | + |
| 241 | + SelectIntervalsToRefine(l_est, merits); |
| 242 | + if (CheckConvergence(cfg, merits)) { |
| 243 | + break; |
| 244 | + } |
| 245 | + |
| 246 | + int num_trials = std::min(static_cast<int>(merits.size()), world_size_); |
| 247 | + GatherAndBroadcastTrialResults(merits, num_trials, l_est); |
| 248 | + SyncGlobalState(); |
| 249 | + total_evals_ += num_trials; |
| 250 | + } |
| 251 | + |
| 252 | + if (IsRoot()) { |
| 253 | + GetOutput() = |
| 254 | + Solution{.argmin = best_point_, .value = best_value_, .iterations = total_evals_, .converged = has_converged_}; |
| 255 | + } |
| 256 | + return true; |
| 257 | +} |
| 258 | + |
| 259 | +bool KondakovVGlobalSearchMPI::PostProcessingImpl() { |
| 260 | + Solution sol; |
| 261 | + if (IsRoot()) { |
| 262 | + sol = GetOutput(); |
| 263 | + } |
| 264 | + |
| 265 | + MPI_Bcast(&sol.argmin, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
| 266 | + MPI_Bcast(&sol.value, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
| 267 | + MPI_Bcast(&sol.iterations, 1, MPI_INT, 0, MPI_COMM_WORLD); |
| 268 | + int converged = sol.converged ? 1 : 0; |
| 269 | + MPI_Bcast(&converged, 1, MPI_INT, 0, MPI_COMM_WORLD); |
| 270 | + sol.converged = (converged != 0); |
| 271 | + |
| 272 | + GetOutput() = sol; |
| 273 | + return true; |
| 274 | +} |
| 275 | + |
| 276 | +} // namespace kondakov_v_global_search |
0 commit comments