@@ -38,28 +38,30 @@ namespace math {
38
38
*/
39
39
template <typename F, typename ... Args>
40
40
inline double gradient_of_f (const F &f, const double &x, const double &xc,
41
- size_t n, std::ostream *msgs, const Args&... args) {
41
+ size_t n, std::ostream *msgs,
42
+ const Args &... args) {
42
43
double gradient = 0.0 ;
43
44
44
45
// Run nested autodiff in this scope
45
46
nested_rev_autodiff nested;
46
47
47
- std::tuple<decltype (deep_copy_vars (args))...>
48
- args_tuple_local_copy ( deep_copy_vars (args)...);
48
+ std::tuple<decltype (deep_copy_vars (args))...> args_tuple_local_copy (
49
+ deep_copy_vars (args)...);
49
50
50
51
Eigen::VectorXd adjoints = Eigen::VectorXd::Zero (count_vars (args...));
51
52
52
- var fx = apply ([&f, &x, &xc, msgs]( auto &&... args) {
53
- return f (x, xc, msgs, args...);
54
- }, args_tuple_local_copy);
53
+ var fx = apply (
54
+ [&f, &x, &xc, msgs]( auto &&... args) { return f (x, xc, msgs, args...); },
55
+ args_tuple_local_copy);
55
56
56
57
fx.grad ();
57
58
58
- apply ([&](auto &&... args) {
59
- accumulate_adjoints (adjoints.data (),
60
- std::forward<decltype (args)>(args)...);
61
- },
62
- std::move (args_tuple_local_copy));
59
+ apply (
60
+ [&](auto &&... args) {
61
+ accumulate_adjoints (adjoints.data (),
62
+ std::forward<decltype (args)>(args)...);
63
+ },
64
+ std::move (args_tuple_local_copy));
63
65
64
66
gradient = adjoints.coeff (n);
65
67
if (is_nan (gradient)) {
@@ -90,13 +92,11 @@ inline double gradient_of_f(const F &f, const double &x, const double &xc,
90
92
* @param args additional arguments to pass to f
91
93
* @return numeric integral of function f
92
94
*/
93
- template <typename F, typename T_a, typename T_b,
94
- typename ... Args,
95
- require_any_st_var<T_a, T_b, Args...>* = nullptr >
95
+ template <typename F, typename T_a, typename T_b, typename ... Args,
96
+ require_any_st_var<T_a, T_b, Args...> * = nullptr >
96
97
inline return_type_t <T_a, T_b, Args...> integrate_1d_impl (
97
- const F &f, const T_a &a, const T_b &b,
98
- double relative_tolerance,
99
- std::ostream *msgs, const Args&... args) {
98
+ const F &f, const T_a &a, const T_b &b, double relative_tolerance,
99
+ std::ostream *msgs, const Args &... args) {
100
100
static const char *function = " integrate_1d" ;
101
101
check_less_or_equal (function, " lower limit" , a, b);
102
102
@@ -105,65 +105,66 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
105
105
106
106
if (a_val == b_val) {
107
107
if (is_inf (a_val)) {
108
- throw_domain_error (function, " Integration endpoints are both" ,
109
- a_val, " " , " " );
108
+ throw_domain_error (function, " Integration endpoints are both" , a_val, " " ,
109
+ " " );
110
110
}
111
111
return var (0.0 );
112
112
} else {
113
- std::tuple<decltype (value_of (args))...>
114
- args_val_tuple (value_of (args)...);
113
+ std::tuple<decltype (value_of (args))...> args_val_tuple (value_of (args)...);
115
114
116
- double integral = integrate (apply ([&f, msgs](auto &&... args) {
117
- return std::bind<double >(f, std::placeholders::_1, std::placeholders::_2,
118
- msgs, args...);
119
- }, args_val_tuple),
120
- a_val, b_val, relative_tolerance);
115
+ double integral = integrate (apply (
116
+ [&f, msgs](auto &&... args) {
117
+ return std::bind<double >(
118
+ f, std::placeholders::_1,
119
+ std::placeholders::_2, msgs, args...);
120
+ },
121
+ args_val_tuple),
122
+ a_val, b_val, relative_tolerance);
121
123
122
124
size_t num_vars_ab = count_vars (a, b);
123
125
size_t num_vars_args = count_vars (args...);
124
- vari** varis =
125
- ChainableStack::instance_->memalloc_ .alloc_array <vari*>(num_vars_ab +
126
- num_vars_args);
127
- double * partials =
128
- ChainableStack::instance_->memalloc_ .alloc_array <double >(num_vars_ab +
129
- num_vars_args);
130
- double * partials_ptr = partials;
126
+ vari **varis = ChainableStack::instance_->memalloc_ .alloc_array <vari *>(
127
+ num_vars_ab + num_vars_args);
128
+ double *partials = ChainableStack::instance_->memalloc_ .alloc_array <double >(
129
+ num_vars_ab + num_vars_args);
130
+ double *partials_ptr = partials;
131
131
132
132
save_varis (varis, a, b, args...);
133
133
134
- for (size_t i = 0 ; i < num_vars_ab + num_vars_args; ++i) {
134
+ for (size_t i = 0 ; i < num_vars_ab + num_vars_args; ++i) {
135
135
partials[i] = 0.0 ;
136
136
}
137
137
138
138
if (!is_inf (a) && is_var<T_a>::value) {
139
- *partials_ptr = apply ([&f, a_val, msgs](auto &&... args) {
140
- return -f (a_val, 0.0 , msgs, args...);
141
- }, args_val_tuple);
139
+ *partials_ptr = apply (
140
+ [&f, a_val, msgs](auto &&... args) {
141
+ return -f (a_val, 0.0 , msgs, args...);
142
+ },
143
+ args_val_tuple);
142
144
partials_ptr++;
143
145
}
144
146
145
147
if (!is_inf (b) && is_var<T_b>::value) {
146
- *partials_ptr = apply ([&f, b_val, msgs](auto &&... args) {
147
- return f (b_val, 0.0 , msgs, args...);
148
- }, args_val_tuple);
148
+ *partials_ptr
149
+ = apply ([&f, b_val, msgs](
150
+ auto &&... args) { return f (b_val, 0.0 , msgs, args...); },
151
+ args_val_tuple);
149
152
partials_ptr++;
150
153
}
151
154
152
155
for (size_t n = 0 ; n < num_vars_args; ++n) {
153
156
*partials_ptr = integrate (
154
- std::bind<double >(gradient_of_f<F, Args...>, f,
155
- std::placeholders::_1, std::placeholders::_2,
156
- n, msgs, args...),
157
- a_val, b_val, relative_tolerance);
157
+ std::bind<double >(gradient_of_f<F, Args...>, f, std::placeholders::_1,
158
+ std::placeholders::_2, n, msgs, args...),
159
+ a_val, b_val, relative_tolerance);
158
160
partials_ptr++;
159
161
}
160
162
161
- return var (new precomputed_gradients_vari (integral,
162
- num_vars_ab + num_vars_args,
163
- varis, partials));
163
+ return var (new precomputed_gradients_vari (
164
+ integral, num_vars_ab + num_vars_args, varis, partials));
164
165
}
165
166
}
166
-
167
+
167
168
/* *
168
169
* Compute the integral of the single variable function f from a to b to within
169
170
* a specified relative tolerance. a and b can be finite or infinite.
@@ -225,9 +226,8 @@ inline return_type_t<T_a, T_b, T_theta> integrate_1d(
225
226
const F &f, const T_a &a, const T_b &b, const std::vector<T_theta> &theta,
226
227
const std::vector<double > &x_r, const std::vector<int > &x_i,
227
228
std::ostream *msgs, const double relative_tolerance = std::sqrt(EPSILON)) {
228
- return integrate_1d_impl (integrate_1d_adapter<F>(f), a, b,
229
- relative_tolerance, msgs,
230
- theta, x_r, x_i);
229
+ return integrate_1d_impl (integrate_1d_adapter<F>(f), a, b, relative_tolerance,
230
+ msgs, theta, x_r, x_i);
231
231
}
232
232
233
233
} // namespace math
0 commit comments