ceph/src/boost/boost/math/differentiation/autodiff_cpp11.hpp

   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
  22 #include <type_traits>
  23 #include <boost/math/tools/mp.hpp>
  24
  25 namespace boost {
  26 namespace math {
  27
  28 namespace mp = tools::meta_programming;
  29
  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>
 361 auto 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
 366 template <typename RealType, size_t... Orders, size_t... Is, typename... RealTypes>
 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>{},
 369                                                                       ca)...)) {
 370   return std::make_tuple(make_fvar_for_tuple<RealType, Orders>(mp::make_index_sequence<Is>{}, ca)...);
 371 }
 372
 373 }  // namespace detail
 374
 375 template <typename RealType, size_t... Orders, typename... RealTypes>
 376 auto 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_v1
 385 }  // namespace differentiation
 386 }  // namespace math
 387 }  // namespace boost