]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | // Copyright Matthew Pulver 2018 - 2019. |
2 | // Distributed under the Boost Software License, Version 1.0. | |
3 | // (See accompanying file LICENSE_1_0.txt or copy at | |
4 | // https://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | // Contributors: | |
7 | // * Kedar R. Bhat - C++11 compatibility. | |
8 | ||
9 | // Notes: | |
10 | // * Any changes to this file should always be downstream from autodiff.cpp. | |
11 | // C++17 is a higher-level language and is easier to maintain. For example, a number of functions which are | |
12 | // lucidly read in autodiff.cpp are forced to be split into multiple structs/functions in this file for | |
13 | // C++11. | |
14 | // * Use of typename RootType and SizeType is a hack to prevent Visual Studio 2015 from compiling functions | |
15 | // that are never called, that would otherwise produce compiler errors. Also forces functions to be inline. | |
16 | ||
17 | #ifndef BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP | |
18 | #error \ | |
19 | "Do not #include this file directly. This should only be #included by autodiff.hpp for C++11 compatibility." | |
20 | #endif | |
21 | ||
1e59de90 TL |
22 | #include <type_traits> |
23 | #include <boost/math/tools/mp.hpp> | |
92f5a8d4 TL |
24 | |
25 | namespace boost { | |
26 | namespace math { | |
1e59de90 TL |
27 | |
28 | namespace mp = tools::meta_programming; | |
29 | ||
92f5a8d4 TL |
30 | namespace differentiation { |
31 | inline namespace autodiff_v1 { | |
32 | namespace detail { | |
33 | ||
34 | template <typename RealType, size_t Order> | |
35 | fvar<RealType, Order>::fvar(root_type const& ca, bool const is_variable) { | |
36 | fvar_cpp11(is_fvar<RealType>{}, ca, is_variable); | |
37 | } | |
38 | ||
39 | template <typename RealType, size_t Order> | |
40 | template <typename RootType> | |
41 | void fvar<RealType, Order>::fvar_cpp11(std::true_type, RootType const& ca, bool const is_variable) { | |
42 | v.front() = RealType(ca, is_variable); | |
43 | if (0 < Order) | |
44 | std::fill(v.begin() + 1, v.end(), static_cast<RealType>(0)); | |
45 | } | |
46 | ||
47 | template <typename RealType, size_t Order> | |
48 | template <typename RootType> | |
49 | void fvar<RealType, Order>::fvar_cpp11(std::false_type, RootType const& ca, bool const is_variable) { | |
50 | v.front() = ca; | |
51 | if (0 < Order) { | |
52 | v[1] = static_cast<root_type>(static_cast<int>(is_variable)); | |
53 | if (1 < Order) | |
54 | std::fill(v.begin() + 2, v.end(), static_cast<RealType>(0)); | |
55 | } | |
56 | } | |
57 | ||
58 | template <typename RealType, size_t Order> | |
59 | template <typename... Orders> | |
60 | get_type_at<RealType, sizeof...(Orders)> fvar<RealType, Order>::at_cpp11(std::true_type, | |
61 | size_t order, | |
62 | Orders...) const { | |
63 | return v.at(order); | |
64 | } | |
65 | ||
66 | template <typename RealType, size_t Order> | |
67 | template <typename... Orders> | |
68 | get_type_at<RealType, sizeof...(Orders)> fvar<RealType, Order>::at_cpp11(std::false_type, | |
69 | size_t order, | |
70 | Orders... orders) const { | |
71 | return v.at(order).at(orders...); | |
72 | } | |
73 | ||
74 | // Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)" | |
75 | template <typename RealType, size_t Order> | |
76 | template <typename... Orders> | |
77 | get_type_at<RealType, sizeof...(Orders)> fvar<RealType, Order>::at(size_t order, Orders... orders) const { | |
78 | return at_cpp11(std::integral_constant<bool, sizeof...(orders) == 0>{}, order, orders...); | |
79 | } | |
80 | ||
81 | template <typename T, typename... Ts> | |
82 | constexpr T product(Ts...) { | |
83 | return static_cast<T>(1); | |
84 | } | |
85 | ||
86 | template <typename T, typename... Ts> | |
87 | constexpr T product(T factor, Ts... factors) { | |
88 | return factor * product<T>(factors...); | |
89 | } | |
90 | ||
91 | // Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)" | |
92 | template <typename RealType, size_t Order> | |
93 | template <typename... Orders> | |
94 | get_type_at<fvar<RealType, Order>, sizeof...(Orders)> fvar<RealType, Order>::derivative( | |
95 | Orders... orders) const { | |
96 | static_assert(sizeof...(Orders) <= depth, | |
97 | "Number of parameters to derivative(...) cannot exceed fvar::depth."); | |
98 | return at(static_cast<size_t>(orders)...) * | |
99 | product(boost::math::factorial<root_type>(static_cast<unsigned>(orders))...); | |
100 | } | |
101 | ||
102 | template <typename RootType, typename Func> | |
103 | class Curry { | |
104 | Func const& f_; | |
105 | size_t const i_; | |
106 | ||
107 | public: | |
108 | template <typename SizeType> // typename SizeType to force inline constructor. | |
109 | Curry(Func const& f, SizeType i) : f_(f), i_(static_cast<std::size_t>(i)) {} | |
110 | template <typename... Indices> | |
111 | RootType operator()(Indices... indices) const { | |
112 | using unsigned_t = typename std::make_unsigned<typename std::common_type<Indices>::type...>::type; | |
113 | return f_(i_, static_cast<unsigned_t>(indices)...); | |
114 | } | |
115 | }; | |
116 | ||
117 | template <typename RealType, size_t Order> | |
118 | template <typename Func, typename Fvar, typename... Fvars> | |
119 | promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_coefficients( | |
120 | size_t const order, | |
121 | Func const& f, | |
122 | Fvar const& cr, | |
123 | Fvars&&... fvars) const { | |
124 | fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0); | |
125 | size_t i = order < order_sum ? order : order_sum; | |
126 | using return_type = promote<fvar<RealType, Order>, Fvar, Fvars...>; | |
127 | return_type accumulator = cr.apply_coefficients( | |
128 | order - i, Curry<typename return_type::root_type, Func>(f, i), std::forward<Fvars>(fvars)...); | |
129 | while (i--) | |
130 | (accumulator *= epsilon) += cr.apply_coefficients( | |
131 | order - i, Curry<typename return_type::root_type, Func>(f, i), std::forward<Fvars>(fvars)...); | |
132 | return accumulator; | |
133 | } | |
134 | ||
135 | template <typename RealType, size_t Order> | |
136 | template <typename Func, typename Fvar, typename... Fvars> | |
137 | promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_coefficients_nonhorner( | |
138 | size_t const order, | |
139 | Func const& f, | |
140 | Fvar const& cr, | |
141 | Fvars&&... fvars) const { | |
142 | fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0); | |
143 | fvar<RealType, Order> epsilon_i = fvar<RealType, Order>(1); // epsilon to the power of i | |
144 | using return_type = promote<fvar<RealType, Order>, Fvar, Fvars...>; | |
145 | return_type accumulator = cr.apply_coefficients_nonhorner( | |
146 | order, Curry<typename return_type::root_type, Func>(f, 0), std::forward<Fvars>(fvars)...); | |
147 | size_t const i_max = order < order_sum ? order : order_sum; | |
148 | for (size_t i = 1; i <= i_max; ++i) { | |
149 | epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); | |
150 | accumulator += epsilon_i.epsilon_multiply( | |
151 | i, | |
152 | 0, | |
153 | cr.apply_coefficients_nonhorner( | |
154 | order - i, Curry<typename return_type::root_type, Func>(f, i), std::forward<Fvars>(fvars)...), | |
155 | 0, | |
156 | 0); | |
157 | } | |
158 | return accumulator; | |
159 | } | |
160 | ||
161 | template <typename RealType, size_t Order> | |
162 | template <typename Func, typename Fvar, typename... Fvars> | |
163 | promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_derivatives( | |
164 | size_t const order, | |
165 | Func const& f, | |
166 | Fvar const& cr, | |
167 | Fvars&&... fvars) const { | |
168 | fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0); | |
169 | size_t i = order < order_sum ? order : order_sum; | |
170 | using return_type = promote<fvar<RealType, Order>, Fvar, Fvars...>; | |
171 | return_type accumulator = | |
172 | cr.apply_derivatives( | |
173 | order - i, Curry<typename return_type::root_type, Func>(f, i), std::forward<Fvars>(fvars)...) / | |
174 | factorial<root_type>(static_cast<unsigned>(i)); | |
175 | while (i--) | |
176 | (accumulator *= epsilon) += | |
177 | cr.apply_derivatives( | |
178 | order - i, Curry<typename return_type::root_type, Func>(f, i), std::forward<Fvars>(fvars)...) / | |
179 | factorial<root_type>(static_cast<unsigned>(i)); | |
180 | return accumulator; | |
181 | } | |
182 | ||
183 | template <typename RealType, size_t Order> | |
184 | template <typename Func, typename Fvar, typename... Fvars> | |
185 | promote<fvar<RealType, Order>, Fvar, Fvars...> fvar<RealType, Order>::apply_derivatives_nonhorner( | |
186 | size_t const order, | |
187 | Func const& f, | |
188 | Fvar const& cr, | |
189 | Fvars&&... fvars) const { | |
190 | fvar<RealType, Order> const epsilon = fvar<RealType, Order>(*this).set_root(0); | |
191 | fvar<RealType, Order> epsilon_i = fvar<RealType, Order>(1); // epsilon to the power of i | |
192 | using return_type = promote<fvar<RealType, Order>, Fvar, Fvars...>; | |
193 | return_type accumulator = cr.apply_derivatives_nonhorner( | |
194 | order, Curry<typename return_type::root_type, Func>(f, 0), std::forward<Fvars>(fvars)...); | |
195 | size_t const i_max = order < order_sum ? order : order_sum; | |
196 | for (size_t i = 1; i <= i_max; ++i) { | |
197 | epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); | |
198 | accumulator += epsilon_i.epsilon_multiply( | |
199 | i, | |
200 | 0, | |
201 | cr.apply_derivatives_nonhorner( | |
202 | order - i, Curry<typename return_type::root_type, Func>(f, i), std::forward<Fvars>(fvars)...) / | |
203 | factorial<root_type>(static_cast<unsigned>(i)), | |
204 | 0, | |
205 | 0); | |
206 | } | |
207 | return accumulator; | |
208 | } | |
209 | ||
210 | template <typename RealType, size_t Order> | |
211 | template <typename SizeType> | |
212 | fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply_cpp11(std::true_type, | |
213 | SizeType z0, | |
214 | size_t isum0, | |
215 | fvar<RealType, Order> const& cr, | |
216 | size_t z1, | |
217 | size_t isum1) const { | |
218 | size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; | |
219 | size_t const m1 = order_sum + isum1 < Order + z1 ? Order + z1 - (order_sum + isum1) : 0; | |
220 | size_t const i_max = m0 + m1 < Order ? Order - (m0 + m1) : 0; | |
221 | fvar<RealType, Order> retval = fvar<RealType, Order>(); | |
222 | for (size_t i = 0, j = Order; i <= i_max; ++i, --j) | |
223 | retval.v[j] = epsilon_inner_product(z0, isum0, m0, cr, z1, isum1, m1, j); | |
224 | return retval; | |
225 | } | |
226 | ||
227 | template <typename RealType, size_t Order> | |
228 | template <typename SizeType> | |
229 | fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply_cpp11(std::false_type, | |
230 | SizeType z0, | |
231 | size_t isum0, | |
232 | fvar<RealType, Order> const& cr, | |
233 | size_t z1, | |
234 | size_t isum1) const { | |
235 | using ssize_t = typename std::make_signed<std::size_t>::type; | |
236 | RealType const zero(0); | |
237 | size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; | |
238 | size_t const m1 = order_sum + isum1 < Order + z1 ? Order + z1 - (order_sum + isum1) : 0; | |
239 | size_t const i_max = m0 + m1 < Order ? Order - (m0 + m1) : 0; | |
240 | fvar<RealType, Order> retval = fvar<RealType, Order>(); | |
241 | for (size_t i = 0, j = Order; i <= i_max; ++i, --j) | |
242 | retval.v[j] = std::inner_product( | |
243 | v.cbegin() + ssize_t(m0), v.cend() - ssize_t(i + m1), cr.v.crbegin() + ssize_t(i + m0), zero); | |
244 | return retval; | |
245 | } | |
246 | ||
247 | template <typename RealType, size_t Order> | |
248 | fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply(size_t z0, | |
249 | size_t isum0, | |
250 | fvar<RealType, Order> const& cr, | |
251 | size_t z1, | |
252 | size_t isum1) const { | |
253 | return epsilon_multiply_cpp11(is_fvar<RealType>{}, z0, isum0, cr, z1, isum1); | |
254 | } | |
255 | ||
256 | template <typename RealType, size_t Order> | |
257 | template <typename SizeType> | |
258 | fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply_cpp11(std::true_type, | |
259 | SizeType z0, | |
260 | size_t isum0, | |
261 | root_type const& ca) const { | |
262 | fvar<RealType, Order> retval(*this); | |
263 | size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; | |
264 | for (size_t i = m0; i <= Order; ++i) | |
265 | retval.v[i] = retval.v[i].epsilon_multiply(z0, isum0 + i, ca); | |
266 | return retval; | |
267 | } | |
268 | ||
269 | template <typename RealType, size_t Order> | |
270 | template <typename SizeType> | |
271 | fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply_cpp11(std::false_type, | |
272 | SizeType z0, | |
273 | size_t isum0, | |
274 | root_type const& ca) const { | |
275 | fvar<RealType, Order> retval(*this); | |
276 | size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; | |
277 | for (size_t i = m0; i <= Order; ++i) | |
278 | if (retval.v[i] != static_cast<RealType>(0)) | |
279 | retval.v[i] *= ca; | |
280 | return retval; | |
281 | } | |
282 | ||
283 | template <typename RealType, size_t Order> | |
284 | fvar<RealType, Order> fvar<RealType, Order>::epsilon_multiply(size_t z0, | |
285 | size_t isum0, | |
286 | root_type const& ca) const { | |
287 | return epsilon_multiply_cpp11(is_fvar<RealType>{}, z0, isum0, ca); | |
288 | } | |
289 | ||
290 | template <typename RealType, size_t Order> | |
291 | template <typename RootType> | |
292 | fvar<RealType, Order>& fvar<RealType, Order>::multiply_assign_by_root_type_cpp11(std::true_type, | |
293 | bool is_root, | |
294 | RootType const& ca) { | |
295 | auto itr = v.begin(); | |
296 | itr->multiply_assign_by_root_type(is_root, ca); | |
297 | for (++itr; itr != v.end(); ++itr) | |
298 | itr->multiply_assign_by_root_type(false, ca); | |
299 | return *this; | |
300 | } | |
301 | ||
302 | template <typename RealType, size_t Order> | |
303 | template <typename RootType> | |
304 | fvar<RealType, Order>& fvar<RealType, Order>::multiply_assign_by_root_type_cpp11(std::false_type, | |
305 | bool is_root, | |
306 | RootType const& ca) { | |
307 | auto itr = v.begin(); | |
308 | if (is_root || *itr != 0) | |
309 | *itr *= ca; // Skip multiplication of 0 by ca=inf to avoid nan, except when is_root. | |
310 | for (++itr; itr != v.end(); ++itr) | |
311 | if (*itr != 0) | |
312 | *itr *= ca; | |
313 | return *this; | |
314 | } | |
315 | ||
316 | template <typename RealType, size_t Order> | |
317 | fvar<RealType, Order>& fvar<RealType, Order>::multiply_assign_by_root_type(bool is_root, | |
318 | root_type const& ca) { | |
319 | return multiply_assign_by_root_type_cpp11(is_fvar<RealType>{}, is_root, ca); | |
320 | } | |
321 | ||
322 | template <typename RealType, size_t Order> | |
323 | template <typename RootType> | |
324 | fvar<RealType, Order>& fvar<RealType, Order>::negate_cpp11(std::true_type, RootType const&) { | |
325 | std::for_each(v.begin(), v.end(), [](RealType& r) { r.negate(); }); | |
326 | return *this; | |
327 | } | |
328 | ||
329 | template <typename RealType, size_t Order> | |
330 | template <typename RootType> | |
331 | fvar<RealType, Order>& fvar<RealType, Order>::negate_cpp11(std::false_type, RootType const&) { | |
332 | std::for_each(v.begin(), v.end(), [](RealType& a) { a = -a; }); | |
333 | return *this; | |
334 | } | |
335 | ||
336 | template <typename RealType, size_t Order> | |
337 | fvar<RealType, Order>& fvar<RealType, Order>::negate() { | |
338 | return negate_cpp11(is_fvar<RealType>{}, static_cast<root_type>(*this)); | |
339 | } | |
340 | ||
341 | template <typename RealType, size_t Order> | |
342 | template <typename RootType> | |
343 | fvar<RealType, Order>& fvar<RealType, Order>::set_root_cpp11(std::true_type, RootType const& root) { | |
344 | v.front().set_root(root); | |
345 | return *this; | |
346 | } | |
347 | ||
348 | template <typename RealType, size_t Order> | |
349 | template <typename RootType> | |
350 | fvar<RealType, Order>& fvar<RealType, Order>::set_root_cpp11(std::false_type, RootType const& root) { | |
351 | v.front() = root; | |
352 | return *this; | |
353 | } | |
354 | ||
355 | template <typename RealType, size_t Order> | |
356 | fvar<RealType, Order>& fvar<RealType, Order>::set_root(root_type const& root) { | |
357 | return set_root_cpp11(is_fvar<RealType>{}, root); | |
358 | } | |
359 | ||
360 | template <typename RealType, size_t Order, size_t... Is> | |
1e59de90 | 361 | auto make_fvar_for_tuple(mp::index_sequence<Is...>, RealType const& ca) |
92f5a8d4 TL |
362 | -> decltype(make_fvar<RealType, zero<Is>::value..., Order>(ca)) { |
363 | return make_fvar<RealType, zero<Is>::value..., Order>(ca); | |
364 | } | |
365 | ||
366 | template <typename RealType, size_t... Orders, size_t... Is, typename... RealTypes> | |
1e59de90 TL |
367 | auto make_ftuple_impl(mp::index_sequence<Is...>, RealTypes const&... ca) |
368 | -> decltype(std::make_tuple(make_fvar_for_tuple<RealType, Orders>(mp::make_index_sequence<Is>{}, | |
92f5a8d4 | 369 | ca)...)) { |
1e59de90 | 370 | return std::make_tuple(make_fvar_for_tuple<RealType, Orders>(mp::make_index_sequence<Is>{}, ca)...); |
92f5a8d4 TL |
371 | } |
372 | ||
373 | } // namespace detail | |
374 | ||
375 | template <typename RealType, size_t... Orders, typename... RealTypes> | |
376 | auto make_ftuple(RealTypes const&... ca) | |
1e59de90 | 377 | -> decltype(detail::make_ftuple_impl<RealType, Orders...>(mp::index_sequence_for<RealTypes...>{}, |
92f5a8d4 TL |
378 | ca...)) { |
379 | static_assert(sizeof...(Orders) == sizeof...(RealTypes), | |
380 | "Number of Orders must match number of function parameters."); | |
1e59de90 | 381 | return detail::make_ftuple_impl<RealType, Orders...>(mp::index_sequence_for<RealTypes...>{}, ca...); |
92f5a8d4 TL |
382 | } |
383 | ||
384 | } // namespace autodiff_v1 | |
385 | } // namespace differentiation | |
386 | } // namespace math | |
387 | } // namespace boost |