Line data Source code
1 : /*************************************************************************
2 : *
3 : * Project
4 : * _____ _____ __ __ _____
5 : * / ____| __ \| \/ | __ \
6 : * ___ _ __ ___ _ __ | | __| |__) | \ / | |__) |
7 : * / _ \| '_ \ / _ \ '_ \| | |_ | ___/| |\/| | ___/
8 : *| (_) | |_) | __/ | | | |__| | | | | | | |
9 : * \___/| .__/ \___|_| |_|\_____|_| |_| |_|_|
10 : * | |
11 : * |_|
12 : *
13 : * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
14 : *
15 : * This software is licensed as described in the file LICENSE, which
16 : * you should have received as part of this distribution. The terms
17 : * among other details are referenced in the official documentation
18 : * seen here : https://akielaries.github.io/openGPMP/ along with
19 : * important files seen in this project.
20 : *
21 : * You may opt to use, copy, modify, merge, publish, distribute
22 : * and/or sell copies of the Software, and permit persons to whom
23 : * the Software is furnished to do so, under the terms of the
24 : * LICENSE file. As this is an Open Source effort, all implementations
25 : * must be of the same methodology.
26 : *
27 : *
28 : *
29 : * This software is distributed on an AS IS basis, WITHOUT
30 : * WARRANTY OF ANY KIND, either express or implied.
31 : *
32 : ************************************************************************/
33 : #include <iostream>
34 : #include <openGPMP/optim/function.hpp>
35 : #include <vector>
36 :
37 1005 : std::vector<double> gpmp::optim::Func::generate_random_point(
38 : const std::vector<double> &lower_bounds,
39 : const std::vector<double> &upper_bounds) const {
40 1005 : if (lower_bounds.size() != upper_bounds.size()) {
41 1 : throw std::invalid_argument(
42 2 : "Lower and upper bounds must have the same dimension.");
43 : }
44 :
45 1004 : std::vector<double> point;
46 3013 : for (size_t i = 0; i < lower_bounds.size(); ++i) {
47 : double random_value =
48 2009 : lower_bounds[i] + static_cast<double>(rand()) / RAND_MAX *
49 2009 : (upper_bounds[i] - lower_bounds[i]);
50 2009 : point.push_back(random_value);
51 : }
52 :
53 1004 : return point;
54 0 : }
55 :
56 : std::vector<double>
57 5 : gpmp::optim::Func::generate_fibonacci_sequence(size_t length) const {
58 5 : std::vector<double> sequence;
59 :
60 : // double golden_ratio = (1.0 + std::sqrt(5.0)) / 2.0;
61 5 : double fibonacci_prev = 0.0;
62 5 : double fibonacci_curr = 1.0;
63 25 : for (size_t i = 0; i < length; ++i) {
64 20 : sequence.push_back(fibonacci_prev);
65 20 : double fibonacci_next = fibonacci_curr + fibonacci_prev;
66 20 : fibonacci_prev = fibonacci_curr;
67 20 : fibonacci_curr = fibonacci_next;
68 : }
69 :
70 10 : return sequence;
71 0 : }
72 :
73 : std::vector<double>
74 4 : gpmp::optim::Func::vector_addition(const std::vector<double> &a,
75 : const std::vector<double> &b) const {
76 4 : if (a.size() != b.size()) {
77 1 : throw std::invalid_argument(
78 2 : "Vector dimensions do not match for addition.");
79 : }
80 :
81 3 : std::vector<double> result;
82 12 : for (size_t i = 0; i < a.size(); ++i) {
83 9 : result.push_back(a[i] + b[i]);
84 : }
85 :
86 3 : return result;
87 0 : }
88 :
89 : std::vector<double>
90 3 : gpmp::optim::Func::vector_subtraction(const std::vector<double> &a,
91 : const std::vector<double> &b) const {
92 3 : if (a.size() != b.size()) {
93 1 : throw std::invalid_argument(
94 2 : "Vector dimensions do not match for subtraction.");
95 : }
96 :
97 2 : std::vector<double> result;
98 8 : for (size_t i = 0; i < a.size(); ++i) {
99 6 : result.push_back(a[i] - b[i]);
100 : }
101 :
102 2 : return result;
103 0 : }
104 :
105 3 : std::vector<double> gpmp::optim::Func::vector_scalar_multiply(
106 : double scalar,
107 : const std::vector<double> &vec) const {
108 3 : std::vector<double> result;
109 12 : for (double value : vec) {
110 9 : result.push_back(scalar * value);
111 : }
112 :
113 3 : return result;
114 0 : }
115 :
116 206 : double gpmp::optim::Func::calculate_midpoint(double a,
117 : double b,
118 : double fraction) const {
119 206 : return a + fraction * (b - a);
120 : }
121 :
122 0 : double gpmp::optim::Func::golden_section_search(
123 : const std::function<double(double)> &func,
124 : double a,
125 : double b,
126 : double tol) {
127 0 : if (!is_valid_interval(a, b)) {
128 0 : throw std::invalid_argument(
129 0 : "Invalid interval: lower bound must be less than upper bound.");
130 : }
131 :
132 0 : const double inv_phi = (std::sqrt(5.0) - 1.0) / 2.0;
133 0 : double x1 = b - inv_phi * (b - a);
134 0 : double x2 = a + inv_phi * (b - a);
135 :
136 0 : return golden_section_search_minimize(func, a, b, tol, x1, x2);
137 : }
138 :
139 0 : double gpmp::optim::Func::linear_interpolation(double x,
140 : double x0,
141 : double x1,
142 : double y0,
143 : double y1) {
144 0 : return y0 + (y1 - y0) / (x1 - x0) * (x - x0);
145 : }
146 :
147 0 : double gpmp::optim::Func::cubic_interpolation(double x,
148 : double x0,
149 : double x1,
150 : double y0,
151 : double y1,
152 : double y0_prime,
153 : double y1_prime) {
154 0 : double h = x1 - x0;
155 0 : double t = (x - x0) / h;
156 0 : double t2 = t * t;
157 0 : double t3 = t2 * t;
158 :
159 0 : double a = 2 * t3 - 3 * t2 + 1;
160 0 : double b = t3 - 2 * t2 + t;
161 0 : double c = t3 - t2;
162 0 : double d = -2 * t3 + 3 * t2;
163 :
164 0 : return a * y0 + b * (h * y0_prime) + c * (h * y1_prime) + d * y1;
165 : }
166 :
167 0 : double gpmp::optim::Func::golden_section_search_minimize(
168 : const std::function<double(double)> &func,
169 : double a,
170 : double b,
171 : double tol,
172 : double x1,
173 : double x2) {
174 0 : double f1 = func(x1);
175 0 : double f2 = func(x2);
176 :
177 0 : while ((b - a) > tol) {
178 0 : if (f1 < f2) {
179 0 : b = x2;
180 0 : x2 = x1;
181 0 : x1 = a + b - x2;
182 0 : f2 = f1;
183 0 : f1 = func(x1);
184 : } else {
185 0 : a = x1;
186 0 : x1 = x2;
187 0 : x2 = a + b - x1;
188 0 : f1 = f2;
189 0 : f2 = func(x2);
190 : }
191 : }
192 :
193 0 : return (f1 < f2) ? x1 : x2;
194 : }
195 :
196 : // Helper methods for multivariate golden section search
197 : // (Implementation of these methods requires more details about the multivariate
198 : // function)
199 :
200 0 : std::vector<double> gpmp::optim::Func::golden_section_search_multivariate(
201 : const std::function<double(const std::vector<double> &)> &func,
202 : const std::vector<double> &lower_bounds,
203 : const std::vector<double> &upper_bounds,
204 : double tol) {
205 0 : if (lower_bounds.size() != upper_bounds.size()) {
206 0 : throw std::invalid_argument(
207 0 : "Lower and upper bounds must have the same dimension.");
208 : }
209 :
210 : // Initialize x1 and x2 based on the golden section ratio for each dimension
211 0 : std::vector<double> x1(lower_bounds.size());
212 0 : std::vector<double> x2(lower_bounds.size());
213 :
214 0 : for (size_t i = 0; i < lower_bounds.size(); ++i) {
215 0 : if (!is_valid_interval(lower_bounds[i], upper_bounds[i])) {
216 0 : throw std::invalid_argument("Invalid interval for dimension " +
217 0 : std::to_string(i));
218 : }
219 :
220 0 : const double inv_phi = (std::sqrt(5.0) - 1.0) / 2.0;
221 0 : x1[i] = upper_bounds[i] - inv_phi * (upper_bounds[i] - lower_bounds[i]);
222 0 : x2[i] = lower_bounds[i] + inv_phi * (upper_bounds[i] - lower_bounds[i]);
223 : }
224 :
225 : return golden_section_search_minimize_multivariate(func,
226 : lower_bounds,
227 : upper_bounds,
228 : tol,
229 : x1,
230 0 : x2);
231 0 : }
232 :
233 2 : double gpmp::optim::Func::random_search(
234 : const std::function<double(const std::vector<double> &)> &func,
235 : const std::vector<double> &lower_bounds,
236 : const std::vector<double> &upper_bounds,
237 : size_t max_iterations) {
238 2 : if (lower_bounds.size() != upper_bounds.size()) {
239 1 : throw std::invalid_argument(
240 2 : "Lower and upper bounds must have the same dimension.");
241 : }
242 :
243 1 : size_t dimension = lower_bounds.size();
244 1 : std::vector<double> best_point(dimension);
245 1 : double best_value = std::numeric_limits<double>::infinity();
246 :
247 1001 : for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
248 : std::vector<double> random_point =
249 1000 : generate_random_point(lower_bounds, upper_bounds);
250 1000 : double value = func(random_point);
251 :
252 1000 : if (value < best_value) {
253 12 : best_value = value;
254 12 : best_point = random_point;
255 : }
256 1000 : }
257 :
258 1 : return best_value;
259 1 : }
260 :
261 : // Curve-fitting methods
262 :
263 : std::vector<double>
264 2 : gpmp::optim::Func::fit_linear(const std::vector<double> &x,
265 : const std::vector<double> &y) {
266 2 : size_t n = x.size();
267 :
268 2 : if (n != y.size() || n < 2) {
269 1 : throw std::invalid_argument(
270 2 : "Invalid input dimensions for linear curve fitting.");
271 : }
272 :
273 1 : double sum_x = 0.0, sum_y = 0.0, sum_xy = 0.0, sum_x_squared = 0.0;
274 :
275 6 : for (size_t i = 0; i < n; ++i) {
276 5 : sum_x += x[i];
277 5 : sum_y += y[i];
278 5 : sum_xy += x[i] * y[i];
279 5 : sum_x_squared += x[i] * x[i];
280 : }
281 :
282 1 : double a =
283 1 : (n * sum_xy - sum_x * sum_y) / (n * sum_x_squared - sum_x * sum_x);
284 1 : double b = (sum_y - a * sum_x) / n;
285 :
286 1 : return {a, b};
287 : }
288 :
289 2 : double gpmp::optim::Func::fibonacci_search(
290 : const std::function<double(const std::vector<double> &)> &func,
291 : const std::vector<double> &lower_bounds,
292 : const std::vector<double> &upper_bounds,
293 : size_t max_iterations) {
294 :
295 2 : if (lower_bounds.size() != upper_bounds.size()) {
296 1 : throw std::invalid_argument(
297 2 : "Lower and upper bounds must have the same dimension.");
298 : }
299 :
300 1 : size_t dimension = lower_bounds.size();
301 : std::vector<double> fib_sequence = generate_fibonacci_sequence(
302 1 : max_iterations + 2); // Adjusted for proper Fibonacci generation
303 :
304 : // Print Fibonacci sequence for debugging
305 1 : std::vector<double> a = lower_bounds;
306 1 : std::vector<double> b = upper_bounds;
307 :
308 3 : for (size_t k = 0; k < max_iterations; ++k) {
309 2 : double lambda = (fib_sequence[max_iterations - k + 1] /
310 2 : fib_sequence[max_iterations - k +
311 2 : 2]); // Corrected lambda calculation
312 :
313 2 : std::vector<double> x1(dimension);
314 2 : std::vector<double> x2(dimension);
315 :
316 6 : for (size_t i = 0; i < dimension; ++i) {
317 4 : x1[i] = a[i] + lambda * (b[i] - a[i]);
318 4 : x2[i] = b[i] + lambda * (a[i] - b[i]);
319 : }
320 :
321 2 : if (func(x1) < func(x2)) {
322 0 : b = x2;
323 : } else {
324 2 : a = x1;
325 : }
326 2 : }
327 :
328 2 : return func(a);
329 1 : }
330 :
331 2 : double gpmp::optim::Func::ternary_search(
332 : const std::function<double(const std::vector<double> &)> &func,
333 : const std::vector<double> &lower_bounds,
334 : const std::vector<double> &upper_bounds,
335 : size_t max_iterations) const {
336 2 : double a = lower_bounds[0];
337 2 : double b = upper_bounds[0];
338 :
339 2 : if (lower_bounds >= upper_bounds) {
340 1 : throw std::invalid_argument("Invalid bounds for ternary search method");
341 : }
342 :
343 101 : for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
344 100 : double m1 = calculate_midpoint(a, b, 1.0 / 3.0);
345 100 : double m2 = calculate_midpoint(a, b, 2.0 / 3.0);
346 :
347 100 : double f_m1 = func({m1});
348 100 : double f_m2 = func({m2});
349 :
350 100 : if (f_m1 < f_m2) {
351 50 : b = m2;
352 : } else {
353 50 : a = m1;
354 : }
355 : }
356 :
357 1 : return calculate_midpoint(a, b, 0.5);
358 : }
359 :
360 : // First-order methods
361 :
362 1 : std::vector<double> gpmp::optim::Func::bisection_method(
363 : const std::function<double(const std::vector<double> &)> &func,
364 : double lower_bound,
365 : double upper_bound,
366 : size_t max_iterations) {
367 1 : if (lower_bound >= upper_bound) {
368 1 : throw std::invalid_argument("Invalid bounds for bisection method.");
369 : }
370 :
371 0 : double a = lower_bound;
372 0 : double b = upper_bound;
373 :
374 0 : for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
375 0 : double midpoint = (a + b) / 2.0;
376 :
377 : // if (func({midpoint}) == 0.0) {
378 0 : if (fabs(func({midpoint}) - 0.0f) <
379 0 : std::numeric_limits<double>::epsilon()) {
380 0 : return {midpoint};
381 0 : } else if (func({midpoint}) * func({a}) < 0) {
382 0 : b = midpoint;
383 : } else {
384 0 : a = midpoint;
385 : }
386 : }
387 :
388 0 : return {(a + b) / 2.0};
389 : }
390 :
391 : // Curve-fitting methods
392 :
393 0 : std::vector<double> gpmp::optim::Func::newton_method(
394 : const std::function<double(const std::vector<double> &)> &func,
395 : const std::function<double(const std::vector<double> &)> &derivative,
396 : double initial_guess,
397 : size_t max_iterations) {
398 0 : double x = initial_guess;
399 :
400 0 : for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
401 0 : double fx = func({x});
402 0 : double dfx = derivative({x});
403 :
404 0 : if ((fabs(dfx) - 0.0f) < std::numeric_limits<double>::epsilon()) {
405 0 : throw std::runtime_error("Newton's method: Derivative is zero.");
406 : }
407 :
408 0 : x = x - fx / dfx;
409 : }
410 :
411 0 : return {x};
412 : }
413 :
414 0 : std::vector<double> gpmp::optim::Func::regula_falsi(
415 : const std::function<double(const std::vector<double> &)> &func,
416 : double lower_bound,
417 : double upper_bound,
418 : size_t max_iterations) {
419 0 : if (func({lower_bound}) * func({upper_bound}) >= 0.0) {
420 0 : throw std::invalid_argument("Invalid bounds for regula falsi method.");
421 : }
422 :
423 0 : double a = lower_bound;
424 0 : double b = upper_bound;
425 :
426 0 : for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
427 0 : double fa = func({a});
428 0 : double fb = func({b});
429 :
430 0 : if (fa * fb >= 0.0) {
431 0 : throw std::runtime_error("Regula falsi method: Function values at "
432 0 : "the bounds have the same sign.");
433 : }
434 :
435 0 : double c = (a * fb - b * fa) / (fb - fa);
436 :
437 0 : if (fabs(func({c}) - 0.0f) < std::numeric_limits<double>::epsilon()) {
438 0 : return {c};
439 0 : } else if (fa * func({c}) < 0.0) {
440 0 : b = c;
441 : } else {
442 0 : a = c;
443 : }
444 : }
445 :
446 0 : return {(a + b) / 2.0};
447 : }
448 :
449 0 : std::vector<double> gpmp::optim::Func::cubic_fit(const std::vector<double> &x,
450 : const std::vector<double> &y) {
451 0 : size_t n = x.size();
452 :
453 0 : if (n != y.size() || n < 4) {
454 0 : throw std::invalid_argument(
455 0 : "Invalid input dimensions for cubic curve fitting.");
456 : }
457 :
458 0 : double sum_x = 0.0, sum_y = 0.0, sum_x_squared = 0.0, sum_x_cubed = 0.0,
459 0 : sum_xy = 0.0, sum_x_squared_y = 0.0, sum_squared = 0.0;
460 :
461 0 : for (size_t i = 0; i < n; ++i) {
462 0 : double xi_squared = x[i] * x[i];
463 0 : double xi_cubed = xi_squared * x[i];
464 :
465 0 : sum_x += x[i];
466 0 : sum_y += y[i];
467 0 : sum_x_squared += xi_squared;
468 0 : sum_x_cubed += xi_cubed;
469 0 : sum_xy += x[i] * y[i];
470 0 : sum_x_squared_y += xi_squared * y[i];
471 : }
472 :
473 0 : double determinant = n * sum_x_squared * sum_x_squared * sum_x_squared -
474 0 : sum_x_squared * sum_x_squared * sum_x * sum_x_squared -
475 0 : sum_x * sum_x * sum_x_squared * sum_x_squared +
476 0 : sum_x_squared * sum_x * sum_x * sum_x_squared +
477 0 : sum_x * sum_x * sum_x * sum_x_squared -
478 0 : n * sum_x * sum_x * sum_x * sum_x;
479 :
480 0 : double a = (n * sum_x_squared * sum_x_squared_y -
481 0 : sum_x_squared * sum_x * sum_xy * sum_x_squared -
482 0 : sum_x * sum_x_squared_y * sum_x_squared +
483 0 : sum_x_squared * sum_x * sum_x * sum_xy +
484 0 : sum_x * sum_x * sum_x_squared * sum_x * sum_y -
485 0 : n * sum_x * sum_x * sum_x * sum_xy) /
486 : determinant;
487 :
488 0 : double b = (n * sum_x_squared * sum_x_squared * sum_xy -
489 0 : sum_x_squared_y * sum_x_squared * sum_x_squared +
490 0 : sum_x * sum_x * sum_x_squared_y * sum_x_squared -
491 0 : sum_x_squared * sum_x * sum_x * sum_x * sum_xy -
492 0 : sum_x_squared * sum_x_squared * sum_x * sum_y +
493 0 : n * sum_x * sum_x * sum_x * sum_x * sum_xy) /
494 : determinant;
495 :
496 0 : double c =
497 0 : (n * sum_x_squared * sum_x_squared * sum_x * sum_xy -
498 0 : sum_x_squared * sum_x_squared_y * sum_x_squared * sum_x +
499 0 : sum_x * sum_x * sum_x_squared_y * sum_x_squared * sum_x_squared -
500 0 : sum_x_squared * sum_x * sum_x * sum_x * sum_x_squared -
501 0 : sum_x * sum_x_squared * sum_x_squared * sum_x * sum_y +
502 0 : n * sum_x * sum_x * sum_x * sum_squared * sum_x * sum_y) /
503 : determinant;
504 :
505 0 : double d =
506 0 : (n * sum_x_squared * sum_x_squared * sum_x * sum_x * sum_xy -
507 0 : sum_x_squared * sum_x_squared_y * sum_x_squared * sum_x_squared *
508 0 : sum_x +
509 0 : sum_x_squared * sum_x * sum_x_squared_y * sum_x_squared *
510 0 : sum_x_squared -
511 0 : sum_x_squared * sum_x * sum_x * sum_x * sum_x_squared -
512 0 : sum_x * sum_x_squared * sum_x_squared * sum_x * sum_x_squared +
513 0 : n * sum_x * sum_x * sum_x * sum_squared * sum_x * sum_x_squared) /
514 : determinant;
515 :
516 0 : return {a, b, c, d};
517 : }
518 :
519 0 : std::vector<double> gpmp::optim::Func::nelder_mead(
520 : const std::function<double(const std::vector<double> &)> &func,
521 : std::vector<double> initial_point,
522 : double tolerance,
523 : size_t max_iterations) {
524 0 : size_t n = initial_point.size();
525 0 : std::vector<std::vector<double>> simplex(n + 1, initial_point);
526 :
527 0 : for (size_t i = 0; i < n; ++i) {
528 0 : simplex[i][i] += 1.0;
529 : }
530 :
531 0 : std::vector<double> values(n + 1);
532 0 : for (size_t i = 0; i <= n; ++i) {
533 0 : values[i] = func(simplex[i]);
534 : }
535 :
536 0 : for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
537 : size_t highest_index =
538 0 : std::distance(values.begin(),
539 0 : std::max_element(values.begin(), values.end()));
540 : size_t lowest_index =
541 0 : std::distance(values.begin(),
542 0 : std::min_element(values.begin(), values.end()));
543 :
544 : std::vector<double> centroid =
545 0 : calculate_centroid(simplex, highest_index);
546 :
547 : // Reflect
548 : std::vector<double> reflected_point =
549 0 : reflect(simplex[highest_index], centroid, 1.0);
550 0 : double reflected_value = func(reflected_point);
551 :
552 0 : if (reflected_value < values[lowest_index]) {
553 : // Expand
554 : std::vector<double> expanded_point =
555 0 : reflect(simplex[highest_index], centroid, 2.0);
556 0 : double expanded_value = func(expanded_point);
557 :
558 0 : if (expanded_value < reflected_value) {
559 0 : simplex[highest_index] = expanded_point;
560 0 : values[highest_index] = expanded_value;
561 : } else {
562 0 : simplex[highest_index] = reflected_point;
563 0 : values[highest_index] = reflected_value;
564 : }
565 0 : } else if (reflected_value >= values[lowest_index] &&
566 0 : reflected_value < values[highest_index]) {
567 0 : simplex[highest_index] = reflected_point;
568 0 : values[highest_index] = reflected_value;
569 : } else {
570 : // Contract
571 : std::vector<double> contracted_point =
572 0 : reflect(simplex[highest_index], centroid, 0.5);
573 0 : double contracted_value = func(contracted_point);
574 :
575 0 : if (contracted_value < values[highest_index]) {
576 0 : simplex[highest_index] = contracted_point;
577 0 : values[highest_index] = contracted_value;
578 : } else {
579 : // Shrink
580 0 : for (size_t i = 0; i <= n; ++i) {
581 0 : if (i != lowest_index) {
582 0 : for (size_t j = 0; j < n; ++j) {
583 0 : simplex[i][j] = 0.5 * (simplex[i][j] +
584 0 : simplex[lowest_index][j]);
585 : }
586 0 : values[i] = func(simplex[i]);
587 : }
588 : }
589 : }
590 0 : }
591 :
592 : // Check convergence
593 0 : double range = calculate_range(values);
594 0 : if (range < tolerance) {
595 0 : return simplex[lowest_index];
596 : }
597 0 : }
598 :
599 0 : return simplex[std::distance(
600 : values.begin(),
601 0 : std::min_element(values.begin(), values.end()))];
602 0 : }
603 :
604 : // Helper methods for Nelder–Mead
605 :
606 0 : std::vector<double> gpmp::optim::Func::calculate_centroid(
607 : const std::vector<std::vector<double>> &simplex,
608 : size_t exclude_index) {
609 0 : size_t n = simplex[0].size();
610 0 : std::vector<double> centroid(n, 0.0);
611 :
612 0 : for (size_t i = 0; i < simplex.size(); ++i) {
613 0 : if (i != exclude_index) {
614 0 : for (size_t j = 0; j < n; ++j) {
615 0 : centroid[j] += simplex[i][j];
616 : }
617 : }
618 : }
619 :
620 0 : for (size_t j = 0; j < n; ++j) {
621 0 : centroid[j] /= (simplex.size() - 1);
622 : }
623 :
624 0 : return centroid;
625 : }
626 :
627 : std::vector<double>
628 0 : gpmp::optim::Func::reflect(const std::vector<double> &point,
629 : const std::vector<double> ¢roid,
630 : double reflection_coefficient) {
631 0 : size_t n = point.size();
632 0 : std::vector<double> reflected_point(n);
633 :
634 0 : for (size_t i = 0; i < n; ++i) {
635 0 : reflected_point[i] =
636 0 : centroid[i] + reflection_coefficient * (centroid[i] - point[i]);
637 : }
638 :
639 0 : return reflected_point;
640 : }
641 :
642 0 : double gpmp::optim::Func::calculate_range(const std::vector<double> &values) {
643 0 : double max_value = *std::max_element(values.begin(), values.end());
644 0 : double min_value = *std::min_element(values.begin(), values.end());
645 :
646 0 : return max_value - min_value;
647 : }
648 :
649 : std::vector<double>
650 0 : gpmp::optim::Func::golden_section_search_minimize_multivariate(
651 : const std::function<double(const std::vector<double> &)> &func,
652 : const std::vector<double> &lower_bounds,
653 : const std::vector<double> &upper_bounds,
654 : double tol,
655 : const std::vector<double> &x1,
656 : const std::vector<double> &x2) {
657 :
658 0 : std::vector<double> best_point;
659 0 : double best_value = std::numeric_limits<double>::infinity();
660 :
661 : // Perform golden section search for each dimension
662 0 : for (size_t i = 0; i < lower_bounds.size(); ++i) {
663 0 : double a = lower_bounds[i];
664 0 : double b = upper_bounds[i];
665 0 : double x1_i = x1[i];
666 0 : double x2_i = x2[i];
667 :
668 0 : double result = golden_section_search_minimize(
669 0 : [&](double t) {
670 0 : std::vector<double> point(x1);
671 0 : point[i] = calculate_midpoint(x1_i, x2_i, t);
672 0 : return func(point);
673 0 : },
674 : a,
675 : b,
676 : tol,
677 : 0.382,
678 : 0.618);
679 :
680 0 : if (result < best_value) {
681 0 : best_value = result;
682 0 : best_point = x1;
683 0 : best_point[i] = calculate_midpoint(x1_i, x2_i, result);
684 : }
685 : }
686 :
687 0 : return best_point;
688 0 : }
689 :
690 0 : bool gpmp::optim::Func::is_valid_interval(double a, double b) {
691 0 : return a < b;
692 : }
|
