1
+ // VolEsti (volume computation and sampling library)
2
+
3
+ // Copyright (c) 2012-2024 Vissarion Fisikopoulos
4
+ // Copyright (c) 2018-2020 Apostolos Chalkis
5
+
6
+ // Licensed under GNU LGPL.3, see LICENCE file
7
+
8
+ #ifndef ODE_SOLVERS_ORACLE_AUTODIFF_FUNCTORS_HPP
9
+ #define ODE_SOLVERS_ORACLE_AUTODIFF_FUNCTORS_HPP
10
+
11
+ #include " Eigen/Eigen"
12
+ #include < autodiff/forward/real.hpp>
13
+ #include < autodiff/forward/real/eigen.hpp>
14
+ #include " cartesian_geom/cartesian_kernel.h"
15
+ #include " cartesian_geom/autopoint.h"
16
+
17
+ struct AutoDiffFunctor {
18
+
19
+ template <typename NT>
20
+ struct parameters {
21
+ unsigned int order;
22
+ NT L; // Lipschitz constant for gradient
23
+ NT m; // Strong convexity constant
24
+ NT kappa; // Condition number
25
+ Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> data;
26
+ parameters () : order(2 ), L(4 ), m(4 ), kappa(1 ){};
27
+ };
28
+
29
+ template <typename NT>
30
+ struct FunctionFunctor_internal {
31
+ using Autopoint = autopoint<NT>;
32
+ using Coeff = typename autopoint<NT>::Coeff;
33
+ using FT = typename autopoint<NT>::FT;
34
+ using Point = typename Cartesian<NT>::Point ;
35
+
36
+ static std::function<FT(const Autopoint &, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &)> pdf;
37
+
38
+ FT static result_internal (const Coeff &x, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &data){
39
+ return pdf (x, data); //
40
+ }
41
+
42
+ // external interface
43
+ Point static differentiate (Point const &x0, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &data) {
44
+ Autopoint x = Autopoint (x0.getCoefficients ()); // cast into autopoint
45
+ auto x1 = x.getCoefficients ();
46
+ Coeff y = autodiff::gradient (result_internal, autodiff::wrt (x1), autodiff::at (x1, data));
47
+ auto result = y.template cast <NT>();
48
+ return -1 * Point (result);
49
+ }
50
+
51
+ NT static result (Point const &x0, const Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> &data) {
52
+ Autopoint x = Autopoint (x0.getCoefficients ()); // cast to autopoint
53
+ auto x1 = x.getCoefficients ();
54
+ return result_internal (x1, data).val ();
55
+ }
56
+ };
57
+
58
+ template <typename Point >
59
+ struct GradientFunctor {
60
+ using NT = typename Point ::FT;
61
+
62
+ FunctionFunctor_internal<NT> F;
63
+ parameters<NT> ¶ms;
64
+
65
+ GradientFunctor (parameters<NT> ¶ms_) : params(params_){};
66
+
67
+ // The index i represents the state vector index
68
+ Point operator ()(unsigned int const &i, std::vector<Point > const &xs, NT const &t) const {
69
+ // std::cout<<"calling gradient functor"<<std::flush;
70
+ if (i == params.order - 1 ) {
71
+ return F.differentiate (xs[0 ], params.data );
72
+ }
73
+ else {
74
+ return xs[i + 1 ]; // returns derivative
75
+ }
76
+ }
77
+ };
78
+
79
+ template <typename Point >
80
+ struct FunctionFunctor {
81
+ using NT = typename Point ::FT;
82
+ parameters<NT> ¶ms;
83
+
84
+ FunctionFunctor (parameters<NT> ¶ms_) : params(params_){};
85
+ // The index i represents the state vector index
86
+ FunctionFunctor_internal<NT> F;
87
+
88
+ NT operator ()(Point const &x) const {
89
+ return F.result (x, params.data );
90
+ }
91
+ };
92
+ };
93
+
94
+ #endif // ODE_SOLVERS_ORACLE_AUTODIFF_FUNCTORS_HPP
0 commit comments