388 lines · plain
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 at4// 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 are12// lucidly read in autodiff.cpp are forced to be split into multiple structs/functions in this file for13// C++11.14// * Use of typename RootType and SizeType is a hack to prevent Visual Studio 2015 from compiling functions15// that are never called, that would otherwise produce compiler errors. Also forces functions to be inline.16 17#ifndef BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP18#error \19 "Do not #include this file directly. This should only be #included by autodiff.hpp for C++11 compatibility."20#endif21 22#include <type_traits>23#include <boost/math/tools/mp.hpp>24 25namespace boost {26namespace math {27 28namespace mp = tools::meta_programming;29 30namespace differentiation {31inline namespace autodiff_v1 {32namespace detail {33 34template <typename RealType, size_t Order>35fvar<RealType, Order>::fvar(root_type const& ca, bool const is_variable) {36 fvar_cpp11(is_fvar<RealType>{}, ca, is_variable);37}38 39template <typename RealType, size_t Order>40template <typename RootType>41void 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 47template <typename RealType, size_t Order>48template <typename RootType>49void 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 58template <typename RealType, size_t Order>59template <typename... Orders>60get_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 66template <typename RealType, size_t Order>67template <typename... Orders>68get_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)"75template <typename RealType, size_t Order>76template <typename... Orders>77get_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 81template <typename T, typename... Ts>82constexpr T product(Ts...) {83 return static_cast<T>(1);84}85 86template <typename T, typename... Ts>87constexpr 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)"92template <typename RealType, size_t Order>93template <typename... Orders>94get_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 102template <typename RootType, typename Func>103class 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 117template <typename RealType, size_t Order>118template <typename Func, typename Fvar, typename... Fvars>119promote<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 135template <typename RealType, size_t Order>136template <typename Func, typename Fvar, typename... Fvars>137promote<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 i144 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 161template <typename RealType, size_t Order>162template <typename Func, typename Fvar, typename... Fvars>163promote<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 183template <typename RealType, size_t Order>184template <typename Func, typename Fvar, typename... Fvars>185promote<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 i192 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 210template <typename RealType, size_t Order>211template <typename SizeType>212fvar<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 227template <typename RealType, size_t Order>228template <typename SizeType>229fvar<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 247template <typename RealType, size_t Order>248fvar<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 256template <typename RealType, size_t Order>257template <typename SizeType>258fvar<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 269template <typename RealType, size_t Order>270template <typename SizeType>271fvar<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 283template <typename RealType, size_t Order>284fvar<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 290template <typename RealType, size_t Order>291template <typename RootType>292fvar<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 302template <typename RealType, size_t Order>303template <typename RootType>304fvar<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 316template <typename RealType, size_t Order>317fvar<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 322template <typename RealType, size_t Order>323template <typename RootType>324fvar<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 329template <typename RealType, size_t Order>330template <typename RootType>331fvar<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 336template <typename RealType, size_t Order>337fvar<RealType, Order>& fvar<RealType, Order>::negate() {338 return negate_cpp11(is_fvar<RealType>{}, static_cast<root_type>(*this));339}340 341template <typename RealType, size_t Order>342template <typename RootType>343fvar<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 348template <typename RealType, size_t Order>349template <typename RootType>350fvar<RealType, Order>& fvar<RealType, Order>::set_root_cpp11(std::false_type, RootType const& root) {351 v.front() = root;352 return *this;353}354 355template <typename RealType, size_t Order>356fvar<RealType, Order>& fvar<RealType, Order>::set_root(root_type const& root) {357 return set_root_cpp11(is_fvar<RealType>{}, root);358}359 360template <typename RealType, size_t Order, size_t... Is>361auto make_fvar_for_tuple(mp::index_sequence<Is...>, RealType const& ca)362 -> decltype(make_fvar<RealType, zero<Is>::value..., Order>(ca)) {363 return make_fvar<RealType, zero<Is>::value..., Order>(ca);364}365 366template <typename RealType, size_t... Orders, size_t... Is, typename... RealTypes>367auto 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>{},369 ca)...)) {370 return std::make_tuple(make_fvar_for_tuple<RealType, Orders>(mp::make_index_sequence<Is>{}, ca)...);371}372 373} // namespace detail374 375template <typename RealType, size_t... Orders, typename... RealTypes>376auto make_ftuple(RealTypes const&... ca)377 -> decltype(detail::make_ftuple_impl<RealType, Orders...>(mp::index_sequence_for<RealTypes...>{},378 ca...)) {379 static_assert(sizeof...(Orders) == sizeof...(RealTypes),380 "Number of Orders must match number of function parameters.");381 return detail::make_ftuple_impl<RealType, Orders...>(mp::index_sequence_for<RealTypes...>{}, ca...);382}383 384} // namespace autodiff_v1385} // namespace differentiation386} // namespace math387} // namespace boost388