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