Optimist  0.0.0
A C++ library for optimization
Loading...
Searching...
No Matches
FiniteDifferences.hh
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi. *
3 * *
4 * The Optimist project is distributed under the BSD 2-Clause License. *
5 * *
6 * Davide Stocco Enrico Bertolazzi *
7 * University of Trento University of Trento *
8 * davide.stocco@unitn.it enrico.bertolazzi@unitn.it *
9\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10
11#pragma once
12
13#ifndef OPTIMIST_FINITE_DIFFERENCES_HH
14#define OPTIMIST_FINITE_DIFFERENCES_HH
15
16#include "Optimist.hh"
17
18namespace Optimist {
19
20 namespace FiniteDifferences {
21
26 template <typename Scalar>
27 requires TypeTrait<Scalar>::IsScalar
28 class Epsilon {
29 OPTIMIST_BASIC_CONSTANTS(Scalar)
30
31 Scalar m_epsilon_1{SQRT_EPSILON};
33 Scalar m_epsilon_2{
34 std::pow(EPSILON, 0.75)};
35 Scalar m_epsilon_3{
36 std::pow(EPSILON, 0.25)};
37 public:
42 Scalar epsilon_1(const Scalar v) const {
43 return (std::abs(v) + 1.0) * this->m_epsilon_1;
44 }
45
50 Scalar epsilon_2(const Scalar v) const {
51 return (std::abs(v) + 1.0) * this->m_epsilon_2;
52 }
53
58 Scalar epsilon_3(const Scalar v) const {
59 return (std::abs(v) + 1.0) * this->m_epsilon_3;
60 }
61 };
62
74 template <typename Vector, typename Scalar>
75 requires TypeTrait<Vector>::IsEigen && TypeTrait<Scalar>::IsScalar
76 inline void SideFiniteDifferences(const Scalar f_0,
77 const Scalar f_1,
78 const Scalar f_2,
79 const Scalar h_1,
80 const Scalar h_2,
81 const Integer i,
82 Vector &out) {
83 const Scalar d_f_1{f_1 - f_0}, d_f_2{f_2 - f_0}, d_h{h_1 - h_2};
84 if constexpr (TypeTrait<Vector>::IsFixed ||
85 TypeTrait<Vector>::IsDynamic) {
86 out(i) = ((h_1 / h_2) * d_f_2 - (h_2 / h_1) * d_f_1) / d_h;
87 } else if constexpr (TypeTrait<Vector>::IsSparse) {
88 out.coeffRef(i) = ((h_1 / h_2) * d_f_2 - (h_2 / h_1) * d_f_1) / d_h;
89 }
90 }
91
102 template <typename Vector, typename Scalar>
103 requires TypeTrait<Vector>::IsEigen && TypeTrait<Scalar>::IsScalar
104 inline void CenteredFiniteDifferences(const Scalar f_l,
105 const Scalar f_c,
106 const Scalar f_r,
107 const Scalar h,
108 const Integer i,
109 Vector &out) {
110 constexpr Scalar EPSILON{std::numeric_limits<Scalar>::epsilon()};
111 const Scalar diff_r{(f_r - f_c) / h}, diff_l{(f_c - f_l) / h};
112 Scalar weight_r{std::abs(diff_r) + EPSILON},
113 weight_l{std::abs(diff_l) + EPSILON};
114 const Scalar weight_max{std::max(weight_r, weight_l)};
115 weight_r = std::sqrt(weight_r / weight_max);
116 weight_l = std::sqrt(weight_l / weight_max);
117 if constexpr (TypeTrait<Vector>::IsFixed ||
118 TypeTrait<Vector>::IsDynamic) {
119 out(i) =
120 (diff_r * weight_l + diff_l * weight_r) / (weight_r + weight_l);
121 } else if constexpr (TypeTrait<Vector>::IsSparse) {
122 out.coeffRef(i) =
123 (diff_r * weight_l + diff_l * weight_r) / (weight_r + weight_l);
124 }
125 }
126
140 template <typename Vector, typename Matrix, typename Scalar>
141 requires TypeTrait<Vector>::IsEigen && TypeTrait<Scalar>::IsScalar
142 inline void SideFiniteDifferences(const Vector &f_0,
143 const Vector &f_1,
144 const Vector &f_2,
145 const Scalar h_1,
146 const Scalar h_2,
147 const Integer i,
148 Matrix &out) {
149 const Scalar d_h{h_1 - h_2};
150 const Scalar t_1{-(h_2 / h_1) / d_h}, t_2{(h_1 / h_2) / d_h};
151 if constexpr (TypeTrait<Vector>::IsFixed ||
152 TypeTrait<Vector>::IsDynamic) {
153 out.col(i) = t_1 * (f_1 - f_0) + t_2 * (f_2 - f_0);
154 } else if constexpr (TypeTrait<Vector>::IsSparse) {
155 const Vector tmp(t_1 * (f_1 - f_0) + t_2 * (f_2 - f_0));
156 for (Eigen::Index j{0}; j < tmp.size(); ++j) {
157 out.coeffRef(j, i) = tmp.coeff(j);
158 }
159 }
160 }
161
174 template <typename Vector, typename Matrix, typename Scalar>
175 requires TypeTrait<Vector>::IsEigen && TypeTrait<Scalar>::IsScalar
176 inline void CenteredFiniteDifferences(const Vector &f_l,
177 const Vector &f_c,
178 const Vector &f_r,
179 const Scalar h,
180 const Integer i,
181 Matrix &&out) {
182 constexpr Scalar EPSILON{std::numeric_limits<Scalar>::epsilon()};
183 const Vector diff_r((f_r - f_c) / h), diff_l((f_c - f_l) / h);
184 if constexpr (TypeTrait<Vector>::IsFixed ||
185 TypeTrait<Vector>::IsDynamic) {
186 Vector weight_r(diff_r.array().abs() + EPSILON),
187 weight_l(diff_l.array().abs() + EPSILON);
188 const Vector weight_max(weight_r.array().max(weight_l.array()));
189 weight_r = (weight_r.array() / weight_max.array()).sqrt();
190 weight_l = (weight_l.array() / weight_max.array()).sqrt();
191 out.col(i) = (diff_r.array() * weight_l.array() +
192 diff_l.array() * weight_r.array()) /
193 (weight_r.array() + weight_l.array());
194 } else if constexpr (TypeTrait<Vector>::IsSparse) {
195 Scalar weight_r, weight_l, weight_max;
196 for (Eigen::Index j{0}; j < diff_r.size(); ++j) {
197 weight_r = std::abs(diff_r.coeff(j)) + EPSILON;
198 weight_l = std::abs(diff_l.coeff(j)) + EPSILON;
199 weight_max = std::max(weight_r, weight_l);
200 weight_r = std::sqrt(weight_r / weight_max);
201 weight_l = std::sqrt(weight_l / weight_max);
202 out.coeffRef(j, i) =
203 (diff_r.coeff(j) * weight_l + diff_l.coeff(j) * weight_r) /
204 (weight_r + weight_l);
205 }
206 }
207 }
208
220 template <typename Function,
221 typename Vector,
222 typename Scalar = typename Vector::Scalar>
223 requires std::
224 is_invocable_r_v<bool, Function, const Vector &, Scalar &> &&
225 TypeTrait<Scalar>::IsScalar && TypeTrait<Vector>::IsEigen
226 inline bool Gradient(Function &&function, const Vector &x, Vector &out) {
227 Epsilon<Scalar> eps;
228 Scalar v_c{0.0};
229 if (!function(x, v_c) || !std::isfinite(v_c)) {
230 return false;
231 }
232 Eigen::Index dim_x{x.size()};
233 if constexpr (TypeTrait<Vector>::IsDynamic ||
234 TypeTrait<Vector>::IsSparse) {
235 out.resize(dim_x);
236 }
237 if constexpr (TypeTrait<Vector>::IsSparse) {
238 out.reserve(dim_x);
239 }
240 out.setZero();
241 for (Eigen::Index i{0}; i < dim_x; ++i) {
242 Vector v_x(x);
243 Scalar tmp;
244 if constexpr (TypeTrait<Vector>::IsFixed ||
245 TypeTrait<Vector>::IsDynamic) {
246 tmp = x(i);
247 } else if constexpr (TypeTrait<Vector>::IsSparse) {
248 tmp = x.coeff(i);
249 }
250 Scalar h_1{eps.epsilon_1(tmp)}, h_2{eps.epsilon_2(tmp)};
251 if constexpr (TypeTrait<Vector>::IsFixed ||
252 TypeTrait<Vector>::IsDynamic) {
253 v_x(i) = tmp + h_1;
254 } else if constexpr (TypeTrait<Vector>::IsSparse) {
255 v_x.coeffRef(i) = tmp + h_1;
256 }
257 Scalar v_r{0.0};
258 bool is_finite_r{function(v_x, v_r) && std::isfinite(v_r)};
259 if constexpr (TypeTrait<Vector>::IsFixed ||
260 TypeTrait<Vector>::IsDynamic) {
261 v_x(i) = tmp - h_1;
262 } else if constexpr (TypeTrait<Vector>::IsSparse) {
263 v_x.coeffRef(i) = tmp - h_1;
264 }
265 Scalar v_l{0.0};
266 bool is_finite_l{function(v_x, v_l) && std::isfinite(v_l)};
267 Eigen::Index ic{(is_finite_r && is_finite_l)
268 ? 0
269 : (is_finite_r ? 1 : (is_finite_l ? -1 : -2))};
270 switch (ic) {
271 case 0: {
272 CenteredFiniteDifferences(v_l, v_c, v_r, h_1, i, out);
273 break;
274 }
275 case 1: {
276 if constexpr (TypeTrait<Vector>::IsFixed ||
277 TypeTrait<Vector>::IsDynamic) {
278 v_x(i) = tmp + h_2;
279 } else if constexpr (TypeTrait<Vector>::IsSparse) {
280 v_x.coeffRef(i) = tmp + h_2;
281 }
282 Scalar v_rr{0.0};
283 bool is_finite_rr{function(v_x, v_rr) && std::isfinite(v_rr)};
284 if (is_finite_rr) {
285 SideFiniteDifferences(v_c, v_r, v_rr, h_1, h_2, i, out);
286 } else {
287 if constexpr (TypeTrait<Vector>::IsFixed ||
288 TypeTrait<Vector>::IsDynamic) {
289 out(i) = (v_r - v_c) / h_1;
290 } else if constexpr (TypeTrait<Vector>::IsSparse) {
291 out.coeffRef(i) = (v_r - v_c) / h_1;
292 }
293 }
294 break;
295 }
296 case -1: {
297 if constexpr (TypeTrait<Vector>::IsFixed ||
298 TypeTrait<Vector>::IsDynamic) {
299 v_x(i) = tmp - h_2;
300 } else if constexpr (TypeTrait<Vector>::IsSparse) {
301 v_x.coeffRef(i) = tmp - h_2;
302 }
303 Scalar v_ll{0.0};
304 bool is_finite_ll{function(v_x, v_ll) && std::isfinite(v_ll)};
305 if (is_finite_ll) {
306 SideFiniteDifferences(v_c, v_l, v_ll, -h_1, -h_2, i, out);
307 } else {
308 if constexpr (TypeTrait<Vector>::IsFixed ||
309 TypeTrait<Vector>::IsDynamic) {
310 out(i) = (v_c - v_l) / h_1;
311 } else if constexpr (TypeTrait<Vector>::IsSparse) {
312 out.coeffRef(i) = (v_c - v_l) / h_1;
313 }
314 }
315 break;
316 }
317 case -2: {
318 if constexpr (TypeTrait<Vector>::IsFixed ||
319 TypeTrait<Vector>::IsDynamic) {
320 out(i) = 0.0;
321 }
322 return false;
323 }
324 }
325 }
326 if constexpr (TypeTrait<Vector>::IsFixed ||
327 TypeTrait<Vector>::IsDynamic) {
328 return out.allFinite();
329 } else if constexpr (TypeTrait<Vector>::IsSparse) {
330 for (Eigen::Index j{0}; j < out.size(); ++j) {
331 if (!std::isfinite(out.coeff(j))) {
332 return false;
333 }
334 }
335 return true;
336 }
337 }
338
351 template <typename Function,
352 typename Vector,
353 typename Matrix,
354 typename Scalar = typename Vector::Scalar>
355 requires std::
356 is_invocable_r_v<bool, Function, const Vector &, Vector &> &&
357 TypeTrait<Scalar>::IsScalar && TypeTrait<Vector>::IsEigen &&
358 TypeTrait<Matrix>::IsEigen
359 inline bool Jacobian(Function &&function, const Vector &x, Matrix &out) {
360 Epsilon<Scalar> eps;
361 Vector v_c;
362 if (!function(x, v_c)) {
363 return false;
364 }
365 if constexpr (TypeTrait<Vector>::IsFixed ||
366 TypeTrait<Vector>::IsDynamic) {
367 if (!v_c.allFinite()) {
368 return false;
369 }
370 } else if constexpr (TypeTrait<Vector>::IsSparse) {
371 for (Eigen::Index r{0}; r < v_c.rows(); ++r) {
372 if (!std::isfinite(v_c.coeff(r))) {
373 return false;
374 }
375 }
376 }
377 Eigen::Index dim_x{x.size()};
378 if constexpr (TypeTrait<Matrix>::IsDynamic ||
379 TypeTrait<Matrix>::IsSparse) {
380 out.resize(dim_x, dim_x);
381 }
382 if constexpr (TypeTrait<Matrix>::IsSparse) {
383 out.reserve(dim_x * dim_x);
384 }
385 out.setZero();
386 for (Eigen::Index j{0}; j < dim_x; ++j) {
387 Vector v_x(x);
388 Scalar tmp;
389 if constexpr (TypeTrait<Vector>::IsFixed ||
390 TypeTrait<Vector>::IsDynamic) {
391 tmp = x(j);
392 } else if constexpr (TypeTrait<Vector>::IsSparse) {
393 tmp = x.coeff(j);
394 }
395 Scalar h_1{eps.epsilon_1(tmp)}, h_2{eps.epsilon_2(tmp)};
396 if constexpr (TypeTrait<Vector>::IsFixed ||
397 TypeTrait<Vector>::IsDynamic) {
398 v_x(j) = tmp + h_1;
399 } else if constexpr (TypeTrait<Vector>::IsSparse) {
400 v_x.coeffRef(j) = tmp + h_1;
401 }
402 Vector v_r;
403 bool is_finite_r{function(v_x, v_r)};
404 if constexpr (TypeTrait<Vector>::IsFixed ||
405 TypeTrait<Vector>::IsDynamic) {
406 is_finite_r = is_finite_r && v_r.allFinite();
407 } else if constexpr (TypeTrait<Vector>::IsSparse) {
408 for (Eigen::Index r{0}; r < v_r.rows(); ++r) {
409 is_finite_r = is_finite_r && std::isfinite(v_r.coeff(r));
410 }
411 }
412 if constexpr (TypeTrait<Vector>::IsFixed ||
413 TypeTrait<Vector>::IsDynamic) {
414 v_x(j) = tmp - h_1;
415 } else if constexpr (TypeTrait<Vector>::IsSparse) {
416 v_x.coeffRef(j) = tmp - h_1;
417 }
418 Vector v_l;
419 bool is_finite_l{function(v_x, v_l)};
420 if constexpr (TypeTrait<Vector>::IsFixed ||
421 TypeTrait<Vector>::IsDynamic) {
422 is_finite_l = is_finite_l && v_l.allFinite();
423 } else if constexpr (TypeTrait<Vector>::IsSparse) {
424 for (Eigen::Index r{0}; r < v_l.rows(); ++r) {
425 is_finite_l = is_finite_l && std::isfinite(v_l.coeff(r));
426 }
427 }
428 Eigen::Index ic{(is_finite_r && is_finite_l)
429 ? 0
430 : (is_finite_r ? 1 : (is_finite_l ? -1 : -2))};
431 switch (ic) {
432 case 0:
433 CenteredFiniteDifferences(v_l, v_c, v_r, h_1, j, out);
434 break;
435 case 1: {
436 if constexpr (TypeTrait<Vector>::IsFixed ||
437 TypeTrait<Vector>::IsDynamic) {
438 v_x(j) = tmp + h_2;
439 } else if constexpr (TypeTrait<Vector>::IsSparse) {
440 v_x.coeffRef(j) = tmp + h_2;
441 }
442 Vector v_rr;
443 bool is_finite_rr{function(v_x, v_rr)};
444 if constexpr (TypeTrait<Vector>::IsFixed ||
445 TypeTrait<Vector>::IsDynamic) {
446 is_finite_rr = is_finite_rr && v_rr.allFinite();
447 } else if constexpr (TypeTrait<Vector>::IsSparse) {
448 for (Eigen::Index r{0}; r < v_rr.rows(); ++r) {
449 is_finite_rr = is_finite_rr && std::isfinite(v_rr.coeff(r));
450 }
451 }
452 if (is_finite_rr) {
453 SideFiniteDifferences(v_c, v_r, v_rr, h_1, h_2, j, out);
454 } else {
455 if constexpr (TypeTrait<Vector>::IsFixed ||
456 TypeTrait<Vector>::IsDynamic) {
457 out.col(j) = (v_r - v_c) / h_1;
458 } else if constexpr (TypeTrait<Vector>::IsSparse) {
459 for (Eigen::Index r{0}; r < v_r.rows(); ++r) {
460 out.coeffRef(r, j) = (v_r.coeff(r) - v_c.coeff(r)) / h_1;
461 }
462 }
463 }
464 break;
465 }
466 case -1: {
467 if constexpr (TypeTrait<Vector>::IsFixed ||
468 TypeTrait<Vector>::IsDynamic) {
469 v_x(j) = tmp - h_2;
470 } else if constexpr (TypeTrait<Vector>::IsSparse) {
471 v_x.coeffRef(j) = tmp - h_2;
472 }
473 Vector v_ll;
474 bool is_finite_ll{function(v_x, v_ll)};
475 if constexpr (TypeTrait<Vector>::IsFixed ||
476 TypeTrait<Vector>::IsDynamic) {
477 is_finite_ll = is_finite_ll && v_ll.allFinite();
478 } else if constexpr (TypeTrait<Vector>::IsSparse) {
479 for (Eigen::Index r{0}; r < v_ll.rows(); ++r) {
480 is_finite_ll = is_finite_ll && std::isfinite(v_ll.coeff(r));
481 }
482 }
483 if (is_finite_ll) {
484 SideFiniteDifferences(v_c, v_l, v_ll, -h_1, -h_2, j, out);
485 } else {
486 if constexpr (TypeTrait<Vector>::IsFixed ||
487 TypeTrait<Vector>::IsDynamic) {
488 out.col(j) = (v_c - v_l) / h_1;
489 } else if constexpr (TypeTrait<Vector>::IsSparse) {
490 for (Eigen::Index r{0}; r < v_l.rows(); ++r) {
491 out.coeffRef(r, j) = (v_c.coeff(r) - v_l.coeff(r)) / h_1;
492 }
493 }
494 }
495 break;
496 }
497 case -2: {
498 if constexpr (TypeTrait<Vector>::IsFixed ||
499 TypeTrait<Vector>::IsDynamic) {
500 out.col(j).setZero();
501 } else if constexpr (TypeTrait<Vector>::IsSparse) {
502 for (typename Matrix::InnerIterator it(out, j); it; ++it) {
503 it.valueRef() = 0.0;
504 }
505 }
506 return false;
507 }
508 }
509 }
510 if constexpr (TypeTrait<Vector>::IsFixed ||
511 TypeTrait<Vector>::IsDynamic) {
512 return out.allFinite();
513 } else if constexpr (TypeTrait<Vector>::IsSparse) {
514 for (Eigen::Index r{0}; r < out.rows(); ++r) {
515 for (Eigen::Index c{0}; c < out.cols(); ++c) {
516 if (!std::isfinite(out.coeff(r, c))) {
517 return false;
518 }
519 }
520 }
521 return true;
522 }
523 }
524
537 template <typename Function,
538 typename Vector,
539 typename Matrix,
540 typename Scalar = typename Vector::Scalar>
541 requires std::
542 is_invocable_r_v<bool, Function, const Vector &, Scalar &> &&
543 TypeTrait<Scalar>::IsScalar && TypeTrait<Vector>::IsEigen &&
544 TypeTrait<Matrix>::IsEigen
545 inline bool Hessian(Function &&function, const Vector &x, Matrix &out) {
546 Epsilon<Scalar> eps;
547 Eigen::Index dim_x{x.size()};
548 if constexpr (TypeTrait<Matrix>::IsDynamic ||
549 TypeTrait<Matrix>::IsSparse) {
550 out.resize(dim_x, dim_x);
551 }
552 if constexpr (TypeTrait<Matrix>::IsSparse) {
553 out.reserve(dim_x * dim_x);
554 }
555 out.setZero();
556 Scalar fc{0.0};
557 if (!function(x, fc) || !std::isfinite(fc)) {
558 return false;
559 }
560 for (Eigen::Index j{0}; j < dim_x; ++j) {
561 Scalar tmp_j;
562 if constexpr (TypeTrait<Vector>::IsFixed ||
563 TypeTrait<Vector>::IsDynamic) {
564 tmp_j = x(j);
565 } else if constexpr (TypeTrait<Vector>::IsSparse) {
566 tmp_j = x.coeff(j);
567 }
568 Scalar h_j{eps.epsilon_3(tmp_j)};
569 Vector v_x(x);
570 if constexpr (TypeTrait<Vector>::IsFixed ||
571 TypeTrait<Vector>::IsDynamic) {
572 v_x(j) = tmp_j + h_j;
573 } else if constexpr (TypeTrait<Vector>::IsSparse) {
574 v_x.coeffRef(j) = tmp_j + h_j;
575 }
576 Scalar fp;
577 if (!function(v_x, fp) || !std::isfinite(fp)) {
578 return false;
579 }
580 if constexpr (TypeTrait<Vector>::IsFixed ||
581 TypeTrait<Vector>::IsDynamic) {
582 v_x(j) = tmp_j - h_j;
583 } else if constexpr (TypeTrait<Vector>::IsSparse) {
584 v_x.coeffRef(j) = tmp_j - h_j;
585 }
586 Scalar fm;
587 if (!function(v_x, fm) || !std::isfinite(fm)) {
588 return false;
589 }
590 if constexpr (TypeTrait<Vector>::IsFixed ||
591 TypeTrait<Vector>::IsDynamic) {
592 out(j, j) = ((fp + fm) - 2.0 * fc) / (h_j * h_j);
593 } else if constexpr (TypeTrait<Vector>::IsSparse) {
594 out.coeffRef(j, j) = ((fp + fm) - 2.0 * fc) / (h_j * h_j);
595 }
596 for (Eigen::Index i{j + 1}; i < dim_x; ++i) {
597 Scalar tmp_i;
598 if constexpr (TypeTrait<Vector>::IsFixed ||
599 TypeTrait<Vector>::IsDynamic) {
600 tmp_i = x(i);
601 } else if constexpr (TypeTrait<Vector>::IsSparse) {
602 tmp_i = x.coeff(i);
603 }
604 Scalar h_i{eps.epsilon_3(tmp_i)};
605 if constexpr (TypeTrait<Vector>::IsFixed ||
606 TypeTrait<Vector>::IsDynamic) {
607 v_x(i) = tmp_i + h_i;
608 v_x(j) = tmp_j + h_j;
609 } else if constexpr (TypeTrait<Vector>::IsSparse) {
610 v_x.coeffRef(i) = tmp_i + h_i;
611 v_x.coeffRef(j) = tmp_j + h_j;
612 }
613 Scalar fpp;
614 if (!function(v_x, fpp) || !std::isfinite(fpp)) {
615 return false;
616 }
617 if constexpr (TypeTrait<Vector>::IsFixed ||
618 TypeTrait<Vector>::IsDynamic) {
619 v_x(i) = tmp_i - h_i;
620 } else if constexpr (TypeTrait<Vector>::IsSparse) {
621 v_x.coeffRef(i) = tmp_i - h_i;
622 }
623 Scalar fmp;
624 if (!function(v_x, fmp) || !std::isfinite(fmp)) {
625 return false;
626 }
627 if constexpr (TypeTrait<Vector>::IsFixed ||
628 TypeTrait<Vector>::IsDynamic) {
629 v_x(j) = tmp_j - h_j;
630 } else if constexpr (TypeTrait<Vector>::IsSparse) {
631 v_x.coeffRef(j) = tmp_j - h_j;
632 }
633 Scalar fmm;
634 if (!function(v_x, fmm) || !std::isfinite(fmm)) {
635 return false;
636 }
637 if constexpr (TypeTrait<Vector>::IsFixed ||
638 TypeTrait<Vector>::IsDynamic) {
639 v_x(i) = tmp_i + h_i;
640 } else if constexpr (TypeTrait<Vector>::IsSparse) {
641 v_x.coeffRef(i) = tmp_i + h_i;
642 }
643 Scalar fpm;
644 if (!function(v_x, fpm) || !std::isfinite(fpm)) {
645 return false;
646 }
647 Scalar h_ij{4.0 * h_i * h_j},
648 value{((fpp + fmm) - (fpm + fmp)) / h_ij};
649 if constexpr (TypeTrait<Vector>::IsFixed ||
650 TypeTrait<Vector>::IsDynamic) {
651 out(j, i) = out(i, j) = value;
652 v_x(i) = tmp_i;
653 } else if constexpr (TypeTrait<Vector>::IsSparse) {
654 out.coeffRef(j, i) = out.coeffRef(i, j) = value;
655 v_x.coeffRef(i) = tmp_i;
656 }
657 }
658 if constexpr (TypeTrait<Vector>::IsFixed ||
659 TypeTrait<Vector>::IsDynamic) {
660 v_x(j) = tmp_j;
661 } else if constexpr (TypeTrait<Vector>::IsSparse) {
662 v_x.coeffRef(j) = tmp_j;
663 }
664 }
665 if constexpr (TypeTrait<Vector>::IsFixed ||
666 TypeTrait<Vector>::IsDynamic) {
667 return out.allFinite();
668 } else if constexpr (TypeTrait<Vector>::IsSparse) {
669 for (Eigen::Index r{0}; r < out.rows(); ++r) {
670 for (Eigen::Index c{0}; c < out.cols(); ++c) {
671 if (!std::isfinite(out.coeff(r, c))) {
672 return false;
673 }
674 }
675 }
676 return true;
677 }
678 }
679
680 } // namespace FiniteDifferences
681
682} // namespace Optimist
683
684#endif // OPTIMIST_FINITE_DIFFERENCES_HH
OPTIMIST_BASIC_CONSTANTS
#define OPTIMIST_BASIC_CONSTANTS(Scalar)
Definition Optimist.hh:70
Optimist::FiniteDifferences::Epsilon
Definition FiniteDifferences.hh:28
Optimist::FiniteDifferences::Epsilon::m_epsilon_3
Scalar m_epsilon_3
Definition FiniteDifferences.hh:35
Optimist::FiniteDifferences::Epsilon::epsilon_2
Scalar epsilon_2(const Scalar v) const
Definition FiniteDifferences.hh:50
Optimist::FiniteDifferences::Epsilon::m_epsilon_2
Scalar m_epsilon_2
Definition FiniteDifferences.hh:33
Optimist::FiniteDifferences::Epsilon::epsilon_1
Scalar epsilon_1(const Scalar v) const
Definition FiniteDifferences.hh:42
Optimist::FiniteDifferences::Epsilon::epsilon_3
Scalar epsilon_3(const Scalar v) const
Definition FiniteDifferences.hh:58
Optimist::FiniteDifferences::Epsilon::m_epsilon_1
Scalar m_epsilon_1
Definition FiniteDifferences.hh:31
Optimist::Function
Class container for the vector-valued function (both input and output are vectors).
Definition Function.hh:221
Optimist::FiniteDifferences
Definition FiniteDifferences.hh:20
Optimist::FiniteDifferences::Hessian
bool Hessian(Function &&function, const Vector &x, Matrix &out)
Definition FiniteDifferences.hh:545
Optimist::FiniteDifferences::Jacobian
bool Jacobian(Function &&function, const Vector &x, Matrix &out)
Definition FiniteDifferences.hh:359
Optimist::FiniteDifferences::CenteredFiniteDifferences
void CenteredFiniteDifferences(const Scalar f_l, const Scalar f_c, const Scalar f_r, const Scalar h, const Integer i, Vector &out)
Definition FiniteDifferences.hh:104
Optimist::FiniteDifferences::Gradient
bool Gradient(Function &&function, const Vector &x, Vector &out)
Definition FiniteDifferences.hh:226
Optimist::FiniteDifferences::SideFiniteDifferences
void SideFiniteDifferences(const Scalar f_0, const Scalar f_1, const Scalar f_2, const Scalar h_1, const Scalar h_2, const Integer i, Vector &out)
Definition FiniteDifferences.hh:76
Optimist
Namespace for the Optimist library.
Definition Optimist.hh:90
Optimist::Integer
OPTIMIST_DEFAULT_INTEGER_TYPE Integer
The Integer type as used for the API.
Definition Optimist.hh:98
Optimist::TypeTrait
Definition Optimist.hh:114