From ab9734f00d7ee0c4dd0e85b02fa448063262056e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Fri, 1 Nov 2024 16:00:09 +0100 Subject: [PATCH 1/3] [advalue] Introduce the ADValue class for templated traceless forward mode This was initially developed within the Dune module dune-fufem (https://gitlab.dune-project.org/fufem/dune-fufem). The class `adolc::ADValue` implementes a traceless forward mode templated with respect to the underlying type `T`, the maximal tracked derivative order `maxOrder` and the domain dimension `dim`. Currently only `maxOrder<=2` is implemented. Using the special value `dynamicDim` allows to use a runtime dimension. In the latter case the default is to use `std::vector` as storage. If support for `boost::pool` is activated, this is used instead. Notable differences to `adtl::adouble` are: * Being able to use other types that `T=double` e.g. allows for low/high precision or interval arithmetics. * Support for second order derivatives. The interface is also flexible for later implementation of higher order. * Fixing the dimension at compile time allows the compiler to do better optimization and generate faster code. Furthermore this guarantees thread-safety (x) by construction and allows to use different dimensions at the same time within a single program. * In the dynamic dimension case, the dimension is set per variable, which allows to have different dimensions at the same time. This is implemented by using one `boost:pool` per dimension instead of a single one for a globally fixed dimension. This also avoids the memory leaks of `adtl::adouble` when changing the global dimension. Since the pools are `thread_local` instead of global, this again provides thread-safety (X). The last point could probably also be implemented for `adtl::adouble`. (x) "Thread-safety" hear does _not_ mean that concurrent access to a single `ADValue` object is safe. Instead it means that concurrent accesss to different objects in different threads is safe and not prone to race-conditions due to internally used global variables. While higher order derivatives could also be computed using `adtl_hov::adouble`, this requires several calls to compute mixed derivatives which is significantly slower compared to computing all at once. Furthermore the latter cannot be used in concurrent threads and with different dimensions at once. --- ADOL-C/include/adolc/CMakeLists.txt | 1 + ADOL-C/include/adolc/advalue.h | 1233 +++++++++++++++++++++++++++ 2 files changed, 1234 insertions(+) create mode 100644 ADOL-C/include/adolc/advalue.h diff --git a/ADOL-C/include/adolc/CMakeLists.txt b/ADOL-C/include/adolc/CMakeLists.txt index ff7fab57..24906c57 100644 --- a/ADOL-C/include/adolc/CMakeLists.txt +++ b/ADOL-C/include/adolc/CMakeLists.txt @@ -11,6 +11,7 @@ install(FILES adtl_indo.h adutilsc.h adutils.h + advalue.h advector.h checkpointing.h convolut.h diff --git a/ADOL-C/include/adolc/advalue.h b/ADOL-C/include/adolc/advalue.h new file mode 100644 index 00000000..def13cf4 --- /dev/null +++ b/ADOL-C/include/adolc/advalue.h @@ -0,0 +1,1233 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef ADOLC_ADVALUE_HH +#define ADOLC_ADVALUE_HH + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#if USE_BOOST_POOL +#include +#endif + +namespace adolc { + +// Utility class providing a consecutive range +// of integral values. +template class range { + T begin_; + T end_; + + struct iterator { + T pos; + T operator*() const { return pos; } + T &operator++() { return ++pos; } + bool operator!=(const iterator &other) const { return pos != other.pos; } + }; + +public: + range(const T &begin, const T &end) : begin_(begin), end_(end) {} + range(const T &size) : range(T(0), size) {} + iterator begin() const { return iterator{begin_}; } + iterator end() const { return iterator{end_}; } +}; + +/** + * \brief Special value to indicate a dynamically selected dimension + * \ingroup TapeLess + * + * This value can be passed as dimension template parameter to `ADValue` + * to indicate that the dimension should be selected dynamically. + */ +inline constexpr std::size_t dynamicDim = + std::numeric_limits::max(); + +template class ADValue; + +/** + * \brief Traits class for checking if F is an `ADValue` + * \ingroup TapeLess + */ +template struct IsADValue; + +template struct IsADValue : public std::false_type {}; + +template +struct IsADValue> : public std::true_type {}; + +/** + * \brief Short cut to `IsADValue::value` + * \ingroup TapeLess + */ +template constexpr bool IsADValue_t = IsADValue::value; + +namespace Impl { + +// Helper function computing the total number of partial derivatives +static constexpr std::size_t fullDerivativeCount(std::size_t dim, + std::size_t maxOrder) { + if (maxOrder == 1) + return dim; + if (maxOrder == 2) + return dim + dim * dim; + else + return 0; +} + +// Helper function computing the reduced number of partial derivatives +// exploiting symmetry +static std::size_t reducedDerivativeCount(std::size_t dim, + std::size_t maxOrder) { + if (maxOrder == 1) + return dim; + if (maxOrder == 2) + return dim + (dim + 1) * dim / 2; + else + return 0; +} + +// Base class storing derivatives up to given maxOrder. +// If the dimension is given statically, we use an +// std::array to store all partial derivatives. +// This implementation stores the full Hessian, +// including duplicates of the mixed derivatives. +// +// For statically known size the compiler generates +// faster code for this, probably by using auto-vectorization +// with SIMD instructions. +template +class DerivativeStorage + : public std::array { + using Base = std::array; + +public: + using Base::operator[]; + using Base::size; + + // Export the dimension. + static constexpr std::size_t dimension() { return dim; } + + void updateDimension(const DerivativeStorage &other) {} + + // Given an index i, this exports the number n, such that + // all derivatives (i,j) with j static constexpr std::size_t derivativeIndex(I... i) { + return storedDerivativeIndex(i...); + } +}; + +// If the dimension is dynamic, we either use a buffer +// managed by boost::pool<> or a std::vector as storage. +#if USE_BOOST_POOL + +// This class uses a simple array-like storage of dynamic size +// with a special allocation pattern: Allocations are done using +// a boost::pool<> which provides a pool allocation strategy for +// fixed block size. To allow for dynamic sizes, this uses a static +// thread_local container storing a pool for each requested size. +// In order to avoid costly lookup for the pool, a pointer to the +// pool is stored in the object on copy and assignment. +template class PoolVector { + class Pool : public boost::pool<> { + public: + Pool(std::size_t size) : boost::pool<>(size * sizeof(T)), size_(size) {} + + std::size_t size() const { return size_; } + + T *alloc() { return reinterpret_cast(this->malloc()); } + + void dealloc(T *p) { this->free(p); } + + private: + std::size_t size_; + }; + + static Pool &threadLocalPool(std::size_t size) { + thread_local std::vector> pools; + if (pools.size() < size + 1) { + pools.resize(size + 1); + pools[size] = std::make_unique(size); + } + return *(pools[size]); + } + + bool nonEmpty() const { return pool_; } + +public: + PoolVector() = default; + + PoolVector(std::size_t size, const T &defaultValue) + : pool_(&threadLocalPool(size)) { + data_ = pool_->alloc(); + for (auto i : range(size)) + data_[i] = defaultValue; + } + + PoolVector(Pool *pool) : pool_(pool) { + if (pool_) + data_ = pool_->alloc(); + } + + PoolVector(Pool *pool, const T &defaultValue) : pool_(pool) { + if (pool_) { + data_ = pool_->alloc(); + for (auto i : range(size())) + data_[i] = defaultValue; + } + } + + PoolVector(const PoolVector &other) : pool_(other.pool_) { + if (other.nonEmpty()) { + data_ = pool_->alloc(); + for (auto i : range(size())) + data_[i] = other.data_[i]; + } + } + + PoolVector(PoolVector &&other) : pool_(other.pool_), data_(other.data_) { + other.pool_ = nullptr; + other.data_ = nullptr; + } + + ~PoolVector() { + if (nonEmpty()) + pool_->dealloc(data_); + } + + PoolVector &operator=(const PoolVector &other) { + if (nonEmpty()) { + if (pool_ == other.pool_) { + for (auto i : range(size())) + data_[i] = other.data_[i]; + return *this; + } else + pool_->dealloc(data_); + } + pool_ = other.pool_; + if (other.nonEmpty()) { + data_ = pool_->alloc(); + for (auto i : range(size())) + data_[i] = other.data_[i]; + } else + data_ = nullptr; + return *this; + } + + PoolVector &operator=(PoolVector &&other) { + if (nonEmpty()) + pool_->dealloc(data_); + pool_ = other.pool_; + data_ = other.data_; + other.pool_ = nullptr; + other.data_ = nullptr; + return *this; + } + + void resize(std::size_t newSize, const T &defaultValue) { + auto *newPool = &threadLocalPool(newSize); + auto *newDerivatives = newPool->alloc(); + if (nonEmpty()) { + for (auto i : range(size())) + newDerivatives[i] = data_[i]; + pool_->dealloc(data_); + } + for (auto i : range(size(), newSize)) + newDerivatives[i] = defaultValue; + pool_ = newPool; + data_ = newDerivatives; + } + + T &operator[](std::size_t k) { return data_[k]; } + + const T &operator[](auto k) const { return data_[k]; } + + std::size_t size() const { + if (pool_) + return pool_->size(); + else + return 0; + } + + Pool *pool() const { return pool_; } + +protected: + Pool *pool_ = nullptr; + T *data_ = nullptr; +}; + +template using DynamicArray = PoolVector; + +template +auto derivativeStorageWithSameDimension( + const DerivativeStorage &other) { + return DerivativeStorage(other.dimension(), other.pool()); +} + +#else + +template using DynamicArray = std::vector; + +template +auto derivativeStorageWithSameDimension( + const DerivativeStorage &other) { + return DerivativeStorage(other.dimension()); +} + +#endif + +// If the dimension is given dynamically, we use an +// std::vector or related boost::pool<> based container +// to store the partial derivatives. +// This implementation only stores the lower triangle +// of the Hessian, avoiding duplicates of the mixed derivatives. +template +class DerivativeStorage : public DynamicArray { + using Base = DynamicArray; + using Base::resize; + +public: + using Base::operator[]; + using Base::size; + + DerivativeStorage() = default; + + DerivativeStorage(const DerivativeStorage &other) = default; + + DerivativeStorage(DerivativeStorage &&other) + : Base(std::move(other)), dimension_(other.dimension_) { + other.dimension_ = 0; + } + + explicit DerivativeStorage(std::size_t dim) + : Base(reducedDerivativeCount(dim, maxOrder), 0), dimension_(dim) {} + + template 0), int> = 0> + DerivativeStorage(std::size_t dim, Args &&...args) + : Base(std::forward(args)...), dimension_(dim) {} + + DerivativeStorage &operator=(const DerivativeStorage &other) = default; + + DerivativeStorage &operator=(DerivativeStorage &&other) { + Base::operator=(std::move(other)); + dimension_ = other.dimension_; + other.dimension_ = 0; + return *this; + } + + std::size_t dimension() const { return dimension_; } + + void updateDimension(const DerivativeStorage &other) { + if (dimension() < other.dimension()) { + auto copy = DerivativeStorage(other.dimension()); + if constexpr (maxOrder >= 1) + for (auto i : range(dimension())) + copy[copy.storedDerivativeIndex(i)] = + (*this)[storedDerivativeIndex(i)]; + if constexpr (maxOrder >= 2) + for (auto i : range(dimension())) + for (auto j : range(storedDerivatives(i))) + copy[copy.storedDerivativeIndex(i, j)] = + (*this)[storedDerivativeIndex(i, j)]; + *this = copy; + } + } + + // Given an index i, this exports the number n, such that + // all derivatives (i,j) with j= j) + return storedDerivativeIndex(i, j); + else + return storedDerivativeIndex(j, i); + } + +protected: + std::size_t dimension_ = 0; +}; + +// Generate a convenience function for accessing the +// derivatives stored in an `ADValue` from indices for +// the deriavative directions. In contranst to adValue.partial(i...) +// this only supports to pass indices i... that are +// actually stored. +template auto storedDerivativeAccess(ADV &adValue) { + return [&](auto... i) -> decltype(auto) { + if constexpr (sizeof...(i) == 0) + return adValue.value(); + else + return adValue.derivativeStorage()[adValue.derivativeStorage() + .storedDerivativeIndex(i...)]; + }; +} + +template auto paddedDerivativeAccess(ADV &adValue) { + using Value = typename ADV::Value; + return [&](auto... i) -> Value { + if (((i < adValue.derivativeStorage().dimension()) && ...)) + return storedDerivativeAccess(adValue)(i...); + else + return Value(0); + }; +} + +template +auto createWithMaxDimension(const ADValue &x, + const ADValue &y) { + if constexpr (dim == dynamicDim) { + if (x.dimension() >= y.dimension()) + return ADValue( + 0, derivativeStorageWithSameDimension(x.derivativeStorage())); + else + return ADValue( + 0, derivativeStorageWithSameDimension(y.derivativeStorage())); + } else + return ADValue(); +} + +} // end namespace Impl + +// Forward declaration +template X inv(const X &x); + +/** + * \brief A value class for traceless automated differentiation + * \ingroup TapeLess + * + * \tparam V Type of the stored value + * \tparam m Maximal derivative order to be computed + * \tparam n Dimension of the domain space (defaults to `dynamicDim` for a +dynamically selected value) + * + * This class behaves like a scalar within some function + * evaluation but simultaneously computes all derivatives + * up to a given maximal order. + * Currently only `m<=2` is supported. + * + * The following provides a usage example, where this class is used + * to evaluate second order derivatives of a trivariate function. + * \code{.cpp} +// Use double as underlying type +using T = double; + +// Raw input vector for function evaluation +auto x_raw = std::array{ 23., 42., 237. }; + +// Function to evaluate +auto f = [](auto x) { + return log(1. + x[0]*x[0]) / (2. + cos(x[1])*sin(x[2])); +}; + +// Typedef to AD-aware version of T. This will track up to +// second order derivatives with respect to three directions. +// If the number of directions is omitted, it will be determined +// dynamically, which is in general significantly slower. +using T_AD = ADValue; + +// Translate x to AD-aware version +auto x = std::array(); +for(std::size_t i=0; i<3; ++i) + x[i] = T_AD(x_raw[i], i); + +// Evaluate the function +auto y = f(x); + +// Print plain function value. Alternatively to y.value() +// we could use y.value() or the cast T(y) +std::cout << y.value() << std::endl; + +// Print all first order partial derivatives +for(std::size_t i=0; i<3; ++i) + std::cout << y.partial(i) << std::endl; + +// Print all second order partial derivatives +for(std::size_t i=0; i<3; ++i) + for(std::size_t j=0; j<3; ++j) + std::cout << y.partial(i, j) << std::endl; + * \endcode + */ +template class ADValue { + + static_assert(m <= 2, "ADValue only supports maxOrder<=2."); + +public: + //! Maximal order of derivatives that will be tracked + static constexpr std::size_t maxOrder = m; + + //! Underlying value type + using Value = V; + + //! Interbal type for flat storage of partial derivatives + using DerivativeStorage = typename Impl::DerivativeStorage; + + // Constructors + + /** + * \brief Defaulted default constructor + * + * This will default initialize the stored value + * and partial derivatives. If `dim==dynamicDim`, + * then dynamic dimension value will be initialized to 0. + */ + ADValue() = default; + + //! Defaulted copy constructor + ADValue(const ADValue &) = default; + + //! Defaulted move constructor + ADValue(ADValue &&) = default; + + /** + * \brief Initialize as constant + * \param value Initialize `ADValue` with this value + * + * Initialize from the given value and consider this + * `ADValue` to be a constant which is invariant under + * the function inputs. + * This assigns all derivatives to zero. + */ + template , int> = 0> + ADValue(const T &value) : derivatives_() { + (*this) = value; + } + + /** + * \brief Initialize as function input + * \param value Initialize `ADValue` with this value + * \param k Consider `ADValue` as `k`-th input variable + * + * Initialize from the given value and consider this + * `ADValue` as `k`-th input variable of the function. + * This assigns the k-th first order derivative to 1 + * and all other ones to zero. + * + * This constructor is only available, if the dimension + * is statically know, i.e. if `dim!=dynamicDim` and if `maxOrder>0`. + */ + template 0), int> = 0> + ADValue(const T &value, std::size_t k) : derivatives_() { + value_ = value; + if constexpr (maxOrder >= 1) + derivatives_[derivatives_.storedDerivativeIndex(k)] = 1; + } + + /** + * \brief Initialize as function input + * \param value Initialize `ADValue` with this value + * \param k Consider `ADValue` as `k`-th input variable + * + * Initialize from the given value and consider this + * `ADValue` as `k`-th input variable of the function. + * This assigns the `k`-th first order derivative to 1 + * and all other ones to zero. + * + * This constructor is only available, if the dimension + * is dynamic, i.e. if `dim==dynamicDim` and if `maxOrder>0`. + * It will initialize the dynamic dimension value by `k+1`. + */ + template 0), int> = 0> + ADValue(const T &value, std::size_t k) : ADValue(value, k, k + 1) {} + + /** + * \brief Initialize as function input + * \param value Initialize `ADValue` with this value + * \param k Consider `ADValue` as `k`-th input variable + * \param d The dynamic domain dimension + * + * Initialize from the given value and consider this + * `ADValue` as `k`-th input variable of the function. + * This assigns the `k`-th first order derivative to 1 + * and all other ones to zero. + * + * This constructor is only available, if the dimension + * is dynamic know, i.e. if `dim==dynamicDim` and if `maxOrder>0`. + */ + template 0), int> = 0> + ADValue(const T &value, std::size_t k, std::size_t d) : derivatives_(d) { + value_ = value; + if constexpr (maxOrder >= 1) + derivatives_[derivatives_.storedDerivativeIndex(k)] = 1; + } + + /** + * \brief Initialize from raw value and derivatives data + * \param value Initializer for value + * \param derivatives Initializer for partial derivative storage + */ + template + ADValue(const T &value, const DerivativeStorage &derivatives) + : value_(value), derivatives_(derivatives) {} + + /** + * \brief Initialize from raw value and derivatives data + * \param value Initializer for value + * \param derivatives Initializer for partial derivative storage + */ + template + ADValue(const T &value, DerivativeStorage &&derivatives) + : value_(value), derivatives_(std::move(derivatives)) {} + + /** + * \brief Assignment from raw value + * + * This will treat this as a constant afterwards. + */ + template , int> = 0> + ADValue &operator=(const T &value) { + value_ = value; + for (auto k : range(derivatives_.size())) + derivatives_[k] = 0; + return *this; + } + + // Assignment operators + + //! Defaulted copy assignment + ADValue &operator=(const ADValue &) = default; + + //! Defaulted move assignment + ADValue &operator=(ADValue &&) = default; + + // Comparison operators + + //! Three way comparison based on the stored value + auto operator<=>(const ADValue &rhs) const { return value_ <=> rhs.value_; } + + //! Three way comparison based on the stored value + template auto operator<=>(const Other &rhs) const { + return value_ <=> rhs; + } + + //! Equality comparison based on the stored value + bool operator==(const ADValue &rhs) const { return value_ == rhs.value_; } + + //! Equality comparison based on the stored value + template bool operator==(const Other &rhs) const { + return value_ == rhs; + } + + // Custom user interface + + //! Dimension of the domain space + auto dimension() const { return derivatives_.dimension(); } + + //! Const access to value + const Value &value() const { return value_; } + + //! Mutable access to value + Value &value() { return value_; } + + /** + * \brief Mutable access to stored derivatives + * \returns Container with stored derivatives + * + * The partial derivatives are stored in an implementation + * defined flat container. + * From the returned storage object `s = adValue.derivativeStorage()` + * the derivative in directions `i, j, ...` can be obtained + * using `s[s.derivativeIndex(i, j, ...)]`. Notice that, + * depending on the internal storage, this may return a + * reference to the same storage location for potentially + * duplicate derivatives, i.e. `s[s.derivativeIndex(1,0)]` + * may internally flip the indices and return the same + * location as `s[s.derivativeIndex(0,1)]`. In contrast + * `s[s.storedDerivativeIndex(i, j, ...)]` will not do + * any index flipping, but only supports to pass the + * indices of actually stored derivatives. + */ + auto &derivativeStorage() { return derivatives_; } + + /** + * \brief Const access to stored derivatives + * \returns Container with stored derivatives + * + * This has the same semantics as the non-const + * method with the same name but only provides + * const access to the stored derivatives. + */ + const auto &derivativeStorage() const { return derivatives_; } + + /** + * \brief Obtain `i...` partial derivative + * + * `adValue.partial(i,j, ...)` returns the partial + * derivative in directions `i,j,...`. Passing no + * arguments returns the value. + */ + template const auto &partial(I... i) const { + if constexpr (sizeof...(i) == 0) + return value_; + else + return derivatives_[derivatives_.derivativeIndex(i...)]; + } + + /** + * \brief Obtain `i...` partial derivative + * + * `adValue.partial(i,j, ...)` returns the partial + * derivative in directions `i,j,...`. Passing no + * arguments returns the value. + */ + template auto &partial(I... i) { + if constexpr (sizeof...(i) == 0) + return value_; + else + return derivatives_[derivatives_.derivativeIndex(i...)]; + } + + /** + * \brief Const cast to `Value` + * + * This is equivalent to `adValue.value()`. + */ + explicit operator Value &() const { return value_; } + + /** + * \brief Mutable cast to `Value` + * + * This is equivalent to `adValue.value()`. + */ + explicit operator const Value &() const { return value_; } + + // Arithmetic operators to make this behave like a raw Value object + + // Sign operators + + ADValue operator+() const { return *this; } + + ADValue operator-() const { return -1 * (*this); } + + // Binary operators for an ADValue and another type + + template + friend ADValue operator+(const ADValue &x, const Other &y) { + auto z = x; + z += y; + return z; + } + + template + friend ADValue operator+(const Other &x, const ADValue &y) { + return y + x; + } + + template + friend ADValue operator-(const ADValue &x, const Other &y) { + auto z = x; + z -= y; + return z; + } + + template + friend ADValue operator-(const Other &x, const ADValue &y) { + auto z = y; + z *= -1; + z += x; + return z; + } + + template + friend ADValue operator*(const ADValue &x, const Other &y) { + auto z = x; + z *= y; + return z; + } + + template + friend ADValue operator*(const Other &x, const ADValue &y) { + return y * x; + } + + template + friend ADValue operator/(const ADValue &x, const Other &y) { + auto z = x; + z /= y; + return z; + } + + template + friend ADValue operator/(const Other &x, const ADValue &y) { + return x * inv(y); + } + + // Compound assignment from another type + + template ADValue &operator+=(const Other &y) { + value_ += y; + return *this; + } + + template ADValue &operator-=(const Other &y) { + value_ -= y; + return *this; + } + + template ADValue &operator*=(const Other &y) { + value_ *= y; + for (auto k : range(derivatives_.size())) + derivatives_[k] *= y; + return *this; + } + + template ADValue &operator/=(const Other &y) { + value_ /= y; + for (auto k : range(derivatives_.size())) + derivatives_[k] /= y; + return *this; + } + + // Binary operators with two ADValues + + friend ADValue operator+(const ADValue &x, const ADValue &y) { + if constexpr (dim != dynamicDim) { + auto z = x; + z += y; + return z; + } else { + auto z = Impl::createWithMaxDimension(x, y); + auto X = Impl::paddedDerivativeAccess(x); + auto Y = Impl::paddedDerivativeAccess(y); + auto Z = Impl::storedDerivativeAccess(z); + Z() = X() + Y(); + if constexpr (maxOrder >= 1) + for (auto i : range(z.derivatives_.dimension())) { + Z(i) = X(i) + Y(i); + if constexpr (maxOrder >= 2) + for (auto j : range(z.derivatives_.storedDerivatives(i))) + Z(i, j) = X(i, j) + Y(i, j); + } + return z; + } + } + + friend ADValue operator-(const ADValue &x, const ADValue &y) { + if constexpr (dim != dynamicDim) { + auto z = x; + z -= y; + return z; + } else { + auto z = Impl::createWithMaxDimension(x, y); + auto X = Impl::paddedDerivativeAccess(x); + auto Y = Impl::paddedDerivativeAccess(y); + auto Z = Impl::storedDerivativeAccess(z); + Z() = X() - Y(); + if constexpr (maxOrder >= 1) + for (auto i : range(z.derivatives_.dimension())) { + Z(i) = X(i) - Y(i); + if constexpr (maxOrder >= 2) + for (auto j : range(z.derivatives_.storedDerivatives(i))) + Z(i, j) = X(i, j) - Y(i, j); + } + return z; + } + } + + friend ADValue operator*(const ADValue &x, const ADValue &y) { + auto z = Impl::createWithMaxDimension(x, y); + auto X = Impl::paddedDerivativeAccess(x); + auto Y = Impl::paddedDerivativeAccess(y); + auto Z = Impl::storedDerivativeAccess(z); + Z() = X() * Y(); + if constexpr (maxOrder >= 1) + for (auto i : range(z.derivatives_.dimension())) { + Z(i) = X(i) * Y() + X() * Y(i); + if constexpr (maxOrder >= 2) + for (auto j : range(z.derivatives_.storedDerivatives(i))) + Z(i, j) = X() * Y(i, j) + Y() * X(i, j) + X(i) * Y(j) + X(j) * Y(i); + } + return z; + } + + friend ADValue operator/(const ADValue &x, const ADValue &y) { + auto z = Impl::createWithMaxDimension(x, y); + auto X = Impl::paddedDerivativeAccess(x); + auto Y = Impl::paddedDerivativeAccess(y); + auto Z = Impl::storedDerivativeAccess(z); + Z() = X() / Y(); + if constexpr (maxOrder >= 1) { + auto Y_squared = Y() * Y(); + for (auto i : range(z.derivatives_.dimension())) { + Z(i) = (X(i) * Y() - X() * Y(i)) / Y_squared; + if constexpr (maxOrder >= 2) { + auto Y_cubed = Y_squared * Y(); + for (auto j : range(z.derivatives_.storedDerivatives(i))) + Z(i, j) = + (X(i, j) * Y_squared - X() * Y(i, j) * Y() + + 2 * X() * Y(i) * Y(j) - (X(i) * Y(j) + X(j) * Y(i)) * Y()) / + Y_cubed; + } + } + } + return z; + } + + // Compound assignment from an ADValue + + ADValue &operator+=(const ADValue &y) { + derivatives_.updateDimension(y.derivatives_); + auto X = Impl::storedDerivativeAccess(*this); + auto Y = Impl::storedDerivativeAccess(y); + X() += Y(); + if constexpr (maxOrder >= 1) + for (auto i : range(y.derivatives_.dimension())) { + X(i) += Y(i); + if constexpr (maxOrder >= 2) + for (auto j : range(y.derivatives_.storedDerivatives(i))) + X(i, j) += Y(i, j); + } + return *this; + } + + ADValue &operator-=(const ADValue &y) { + derivatives_.updateDimension(y.derivatives_); + auto X = Impl::storedDerivativeAccess(*this); + auto Y = Impl::storedDerivativeAccess(y); + X() -= Y(); + if constexpr (maxOrder >= 1) + for (auto i : range(y.derivatives_.dimension())) { + X(i) -= Y(i); + if constexpr (maxOrder >= 2) + for (auto j : range(y.derivatives_.storedDerivatives(i))) + X(i, j) -= Y(i, j); + } + return *this; + } + + ADValue &operator*=(const ADValue &y) { + derivatives_.updateDimension(y.derivatives_); + auto X = Impl::storedDerivativeAccess(*this); + auto Y = Impl::storedDerivativeAccess(y); + if constexpr (maxOrder >= 2) { + for (auto i : range(y.derivatives_.dimension())) { + for (auto j : range(y.derivatives_.storedDerivatives(i))) + X(i, j) = X() * Y(i, j) + Y() * X(i, j) + X(i) * Y(j) + X(j) * Y(i); + for (auto j : range(y.derivatives_.storedDerivatives(i), + derivatives_.storedDerivatives(i))) + X(i, j) = Y() * X(i, j) + X(j) * Y(i); + } + for (auto i : + range(y.derivatives_.dimension(), derivatives_.dimension())) { + for (auto j : range(y.derivatives_.dimension())) + X(i, j) = Y() * X(i, j) + X(i) * Y(j); + for (auto j : range(y.derivatives_.dimension(), + derivatives_.storedDerivatives(i))) + X(i, j) = Y() * X(i, j); + } + } + if constexpr (maxOrder >= 1) { + for (auto i : range(y.derivatives_.dimension())) + X(i) = X(i) * Y() + X() * Y(i); + for (auto i : range(y.derivatives_.dimension(), derivatives_.dimension())) + X(i) = X(i) * Y(); + } + X() = X() * Y(); + return *this; + } + + ADValue &operator/=(const ADValue &y) { + (*this) *= inv(y); + return *this; + } + +private: + Value value_; + DerivativeStorage derivatives_; +}; + +/** + * \brief Helper function for implementing `ADValue`-aware nonlinear functions. + * \ingroup TapeLess + * \param x Inner value, this can be an `ADValue` or a raw value + * \param f Outer function + * + * This computes the composition with custom nonlinear function. + * The function `f` should implement derivatives of order `k` + * using `f(std::integral_constant(), x)`. + * If `x` is a raw value, this simply returns + * `f(std::integral_constant(), x)`. + * If `x` is an `ADValue`, then the return value is obtained + * by composition of the inner function represented by `x` + * with `f`, i.e. the partial derivatives are set using the chain rule. + * Currently only derivatives up to order two are supported. + */ +template +auto adCompose(const Value &x, const Derivatives &f) { + auto _0 = std::integral_constant(); + auto _1 = std::integral_constant(); + auto _2 = std::integral_constant(); + if constexpr (IsADValue_t) { + // (f o x)' = f'(x)*x' + // (f o x)'' = f''(x)*x'*x' + f'(x)*x'' + auto y = x; + auto X = Impl::storedDerivativeAccess(x); + auto Y = Impl::storedDerivativeAccess(y); + Y() = f(_0, X()); + if constexpr (Value::maxOrder >= 1) { + auto df_x = f(_1, X()); + for (auto i : range(y.derivativeStorage().dimension())) + Y(i) *= df_x; + if constexpr (Value::maxOrder >= 2) { + auto ddf_x = f(_2, X()); + for (auto i : range(y.derivativeStorage().dimension())) + for (auto j : range(y.derivativeStorage().storedDerivatives(i))) + Y(i, j) = X(i, j) * df_x + ddf_x * X(i) * X(j); + } + } + return y; + } else + return f(_0, x); +} + +// Some scalar mathematical functions + +/** + * \brief `ADValue`-aware nonlinear reciprocal function + * \ingroup TapeLess + */ +template X inv(const X &x) { + return adCompose(x, [](auto k, auto x) { + constexpr auto order = decltype(k)::value; + if constexpr (order == 0) + return 1 / x; + if constexpr (order == 1) + return -1 / (x * x); + if constexpr (order == 2) + return 2 / (x * x * x); + static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); + }); +} + +/** + * \brief `ADValue`-aware fabs-function + * \ingroup TapeLess + */ +template auto fabs(const X &x) { + using std::fabs; + return adCompose(x, [](auto k, auto x) { + constexpr auto order = decltype(k)::value; + if constexpr (order == 0) + return fabs(x); + else if constexpr (order == 1) + return x < 0 ? -1. : 1.; + else + return 0; + }); +} + +/** + * \brief `ADValue`-aware abs-function + * \ingroup TapeLess + */ +template auto abs(const X &x) { + using std::abs; + return adCompose(x, [](auto k, auto x) { + constexpr auto order = decltype(k)::value; + if constexpr (order == 0) + return abs(x); + else if constexpr (order == 1) + return x < 0 ? -1. : 1.; + else + return 0; + }); +} + +/** + * \brief `ADValue`-aware sin-function + * \ingroup TapeLess + */ +template auto sin(const X &x) { + using std::cos; + using std::sin; + return adCompose(x, [](auto order, auto x) { + if ((order % 4) == 0) + return sin(x); + if ((order % 4) == 1) + return cos(x); + if ((order % 4) == 2) + return -sin(x); + return -cos(x); + }); +} + +/** + * \brief `ADValue`-aware cos-function + * \ingroup TapeLess + */ +template auto cos(const X &x) { + using std::cos; + using std::sin; + return adCompose(x, [](auto order, auto x) { + if ((order % 4) == 0) + return cos(x); + if ((order % 4) == 1) + return -sin(x); + if ((order % 4) == 2) + return -cos(x); + return sin(x); + }); +} + +/** + * \brief `ADValue`-aware exp-function + * \ingroup TapeLess + */ +template auto exp(const X &x) { + using std::exp; + return adCompose(x, [](auto order, auto x) { return exp(x); }); +} + +/** + * \brief `ADValue`-aware log-function + * \ingroup TapeLess + */ +template auto log(const X &x) { + using std::log; + using std::pow; + return adCompose(x, [](auto k, auto x) { + constexpr auto order = decltype(k)::value; + if constexpr (order == 0) + return log(x); + if constexpr (order == 1) + return 1. / x; + if constexpr (order == 2) + return -1. / (x * x); + static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); + }); +} + +/** + * \brief `ADValue`-aware sqrt-function + * \ingroup TapeLess + */ +template auto sqrt(const X &x) { + using std::sqrt; + return adCompose(x, [](auto k, auto x) { + constexpr auto order = decltype(k)::value; + if constexpr (order == 0) + return sqrt(x); + if constexpr (order == 1) + return 1. / (2. * sqrt(x)); + if constexpr (order == 2) + return -1. / (4. * x * sqrt(x)); + static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); + }); +} + +/** + * \brief `ADValue`-aware pow-function with constant exponent + * \ingroup TapeLess + */ +template , int> = 0> +auto pow(const X &x, const Y &y) { + using std::pow; + return adCompose(x, [y](auto k, auto x) { + constexpr auto order = decltype(k)::value; + if constexpr (order == 0) + return pow(x, y); + if constexpr (order == 1) + return y * pow(x, y - 1); + if constexpr (order == 2) + return y * (y - 1.) * pow(x, y - 2.); + static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); + }); +} + +/** + * \brief `ADValue`-aware pow-function with constant base + * \ingroup TapeLess + */ +template , int> = 0> +auto pow(const X &x, const Y &y) { + using std::log; + using std::pow; + return adCompose(y, [x](auto k, auto y) { + constexpr auto order = decltype(k)::value; + if constexpr (order == 0) + return pow(x, y); + if constexpr (order == 1) + return log(x) * pow(x, y); + if constexpr (order == 2) + return log(x) * log(x) * pow(x, y); + static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); + }); +} + +/** + * \brief `ADValue`-aware pow-function variable base and exponent + * \ingroup TapeLess + */ +template , int> = 0> +auto pow(const V &x, const V &y) { + using std::log; + using std::pow; + auto z = Impl::createWithMaxDimension(x, y); + auto X = Impl::paddedDerivativeAccess(x); + auto Y = Impl::paddedDerivativeAccess(y); + auto Z = Impl::storedDerivativeAccess(z); + // Here we use: + // z = pow(x, y) = exp(y*log(x)) + // z' = exp(y*log(x)) * (y'*log(x) + y/x*x') = y*x^(y-1)*x' + log(x)*z*y' + Z() = pow(X(), Y()); + if constexpr (V::maxOrder >= 1) { + auto pow_X_Y1 = pow(X(), Y() - 1); + auto log_X = log(X()); + auto tmp1 = Y() * pow_X_Y1; + auto tmp2 = Z() * log_X; + for (auto i : range(z.derivativeStorage().dimension())) + Z(i) = tmp1 * X(i) + tmp2 * Y(i); + if constexpr (V::maxOrder >= 2) { + auto tmp3 = Y() * (Y() - 1) * pow(X(), Y() - 2); + auto tmp4 = tmp2 * log_X; + auto tmp5 = (1 + log_X * Y()) * pow_X_Y1; + for (auto i : range(z.derivativeStorage().dimension())) + for (auto j : range(z.derivativeStorage().storedDerivatives(i))) + Z(i, j) = tmp1 * X(i, j) + tmp2 * Y(i, j) + tmp3 * X(i) * X(j) + + tmp4 * Y(i) * Y(j) + tmp5 * (X(i) * Y(j) + Y(i) * X(j)); + } + } + return z; +} + +} // end namespace adolc + +#endif // DUNE_FUFEM_FUNCTIONS_ADVALUE_HH From bfa05459c081f54fb22f92429095da6af19bb10c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Fri, 15 Nov 2024 23:29:49 +0100 Subject: [PATCH 2/3] [test][advalue] Add rudimentary test of ADValue The test is still missing checks for many features. --- ADOL-C/boost-test/CMakeLists.txt | 1 + ADOL-C/boost-test/tracelessADValue.cpp | 392 +++++++++++++++++++++++++ 2 files changed, 393 insertions(+) create mode 100644 ADOL-C/boost-test/tracelessADValue.cpp diff --git a/ADOL-C/boost-test/CMakeLists.txt b/ADOL-C/boost-test/CMakeLists.txt index 137d40de..17bd37e9 100644 --- a/ADOL-C/boost-test/CMakeLists.txt +++ b/ADOL-C/boost-test/CMakeLists.txt @@ -17,6 +17,7 @@ set(SOURCE_FILES adouble.cpp main.cpp traceCompositeTests.cpp + tracelessADValue.cpp tracelessCompositeTests.cpp tracelessOperatorScalar.cpp tracelessOperatorVector.cpp diff --git a/ADOL-C/boost-test/tracelessADValue.cpp b/ADOL-C/boost-test/tracelessADValue.cpp new file mode 100644 index 00000000..8afbbe08 --- /dev/null +++ b/ADOL-C/boost-test/tracelessADValue.cpp @@ -0,0 +1,392 @@ +#define BOOST_TEST_DYN_LINK + +#include +#include + +#include + +//************************************************** +//* Test for the traceless forward mode +//* based on adolc::ADValue +//* +//* Author: Carsten Graeser +//************************************************** + +//************************************************** +//* Using ADValue requires C++20 or later +//************************************************** +#if __cplusplus >= 202002L + +#include + +//************************************************** +//* Some utilities for testing +//************************************************** + +/** + * \brief Helper function for checking an ADValue + * + * \param checkDim Correct dimension of argument space + * \param value Correct value + */ +template +void check(const adolc::ADValue &adValue, std::size_t checkDim, + T value) { + BOOST_TEST(adValue.dimension() == checkDim); + // Check value + BOOST_TEST(adValue.partial() == value); +} + +/** + * \brief Helper function for checking an ADValue + * + * \param checkDim Correct dimension of argument space + * \param value Correct value + * \param d Unary callback computing correct first order partial derivatives + */ +template +void check(const adolc::ADValue &adValue, std::size_t checkDim, + T value, D_Callback &&d) { + check(adValue, checkDim, value); + // Check 1st order derivatives + if constexpr (order >= 1) + for (std::size_t i = 0; i < adValue.dimension(); ++i) + BOOST_TEST(adValue.partial(i) == d(i)); +} + +/** + * \brief Helper function for checking an ADValue + * + * \param checkDim Correct dimension of argument space + * \param value Correct value + * \param d Unary callback computing correct first order partial derivatives + * \param dd Binary callback computing correct second order partial derivatives + */ +template +void check(const adolc::ADValue &adValue, std::size_t checkDim, + T value, D_Callback &&d, DD_Callback &&dd) { + check(adValue, checkDim, value, d); + // Check 2nd order derivatives + if constexpr (order >= 2) + for (std::size_t i = 0; i < adValue.dimension(); ++i) + for (std::size_t j = 0; j < adValue.dimension(); ++j) + BOOST_TEST(adValue.partial(i, j) == dd(i, j)); +} + +/** + * \brief Helper function object returning zero for any input arguments + */ +constexpr auto zero = [](auto... i) { return 0; }; + +/** + * \brief Create example value of static size + * + * The resulting ADValue's value is set to seed, + * while the (i0,...in)-th partial derivative + * is seed*i0*...*in. + */ +template +auto exampleValue(T seed) { + auto x = adolc::ADValue(); + x.partial() = seed; + if constexpr (order >= 1) { + for (std::size_t i = 0; i < dim; ++i) + x.partial(i) = seed * i; + } + if constexpr (order >= 2) { + for (std::size_t i = 0; i < dim; ++i) + for (std::size_t j = 0; j < dim; ++j) + x.partial(i, j) = seed * i * j; + } + return x; +} + +/** + * \brief Create example value of dynamic size + * + * The resulting ADValue's value is set to seed, + * while the (i0,...in)-th partial derivative + * is seed*i0*...*in. + */ +template +auto exampleValue(std::size_t dim, T seed) { + if constexpr (order == 0) + return adolc::ADValue(seed); + else { + if (dim == 0) + return adolc::ADValue(seed); + auto x = adolc::ADValue(seed, 0, dim); + if constexpr (order >= 1) { + for (std::size_t i = 0; i < dim; ++i) + x.partial(i) = seed * i; + } + if constexpr (order >= 2) { + for (std::size_t i = 0; i < dim; ++i) + for (std::size_t j = 0; j < dim; ++j) + x.partial(i, j) = seed * i * j; + } + return x; + } +} + +/** + * \brief Call check with a few combinations of base type order and dimension + */ +template void defaultTestSuite(Check &&checkCallBack) { + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); + checkCallBack.template operator()(); +} + +//************************************************** +//* Test of individual feature of ADValue +//************************************************** + +BOOST_AUTO_TEST_SUITE(traceless_advalue) + +BOOST_AUTO_TEST_CASE(ADValueConstructor) { + +#if USE_BOOST_POOL + std::cout << "Testing ADValue with boost-pool support." << std::endl; +#else + std::cout << "Testing ADValue without boost-pool support." << std::endl; +#endif + + defaultTestSuite([]() { + // Check default dimension value + using ADValue_dynamic = adolc::ADValue; + using ADValue_default = adolc::ADValue; + BOOST_TEST((std::is_same_v)); + + // Construct as constant + T x = 42; + check(adolc::ADValue(x), dim, x, zero, zero); + check(adolc::ADValue(x), 0, x, zero, zero); + + if constexpr (order >= 1) { + // Construct as k-th input argument + for (std::size_t k = 0; k < dim; ++k) { + auto D = [&](auto i) { return i == k; }; + check(adolc::ADValue(x, k), dim, x, D, zero); + check(adolc::ADValue(x, k, dim), dim, x, D, + zero); + } + } + }); +} + +BOOST_AUTO_TEST_CASE(ADValueAssign) { + defaultTestSuite([]() { + T aValue = 42; + T bValue = 23; + T cValue = 237; + + { + auto a = adolc::ADValue(aValue); + auto b = adolc::ADValue(bValue); + check(a, dim, aValue, zero, zero); + a = b; + check(a, dim, bValue, zero, zero); + a = cValue; + check(a, dim, cValue, zero, zero); + } + { + auto a = adolc::ADValue(aValue); + auto b = adolc::ADValue(bValue); + check(a, 0, aValue, zero, zero); + a = b; + check(a, 0, bValue, zero, zero); + a = cValue; + check(a, 0, cValue, zero, zero); + } + + if constexpr ((dim > 0) and (order > 0)) { + { + auto a = adolc::ADValue(aValue, 0); + auto b = adolc::ADValue(bValue, dim - 1); + check(a, dim, aValue, [](auto i) { return i == 0; }, zero); + a = b; + check(a, dim, bValue, [](auto i) { return i == (dim - 1); }, zero); + a = cValue; + check(a, dim, cValue, zero, zero); + } + { + auto a = adolc::ADValue(aValue, 0, 1); + auto b = + adolc::ADValue(bValue, dim - 1, dim); + check(a, 1, aValue, [](auto i) { return i == 0; }, zero); + a = b; + check(a, dim, bValue, [](auto i) { return i == (dim - 1); }, zero); + a = cValue; + check(a, dim, cValue, zero, zero); + } + } + }); +} + +BOOST_AUTO_TEST_CASE(ADValueSum) { + defaultTestSuite([]() { + T aValue = 42; + T bValue = 23; + + { + auto a = exampleValue(aValue); + auto b = exampleValue(bValue); + + auto c = a + b; + check( + c, dim, aValue + bValue, + [&](auto i) { return (aValue + bValue) * i; }, + [&](auto i, auto j) { return (aValue + bValue) * i * j; }); + + auto d = a; + d += b; + check( + d, dim, aValue + bValue, + [&](auto i) { return (aValue + bValue) * i; }, + [&](auto i, auto j) { return (aValue + bValue) * i * j; }); + } + + { + auto a = exampleValue(dim, aValue); + auto b = exampleValue(dim, bValue); + + auto c = a + b; + check( + c, dim, aValue + bValue, + [&](auto i) { return (aValue + bValue) * i; }, + [&](auto i, auto j) { return (aValue + bValue) * i * j; }); + + auto d = a; + d += b; + check( + d, dim, aValue + bValue, + [&](auto i) { return (aValue + bValue) * i; }, + [&](auto i, auto j) { return (aValue + bValue) * i * j; }); + } + }); +} + +BOOST_AUTO_TEST_CASE(ADValueDiff) { + defaultTestSuite([]() { + T aValue = 42; + T bValue = 23; + + { + auto a = exampleValue(aValue); + auto b = exampleValue(bValue); + + auto c = a - b; + check( + c, dim, aValue - bValue, + [&](auto i) { return (aValue - bValue) * i; }, + [&](auto i, auto j) { return (aValue - bValue) * i * j; }); + + auto d = a; + d -= b; + check( + d, dim, aValue - bValue, + [&](auto i) { return (aValue - bValue) * i; }, + [&](auto i, auto j) { return (aValue - bValue) * i * j; }); + } + + { + auto a = exampleValue(dim, aValue); + auto b = exampleValue(dim, bValue); + + auto c = a - b; + check( + c, dim, aValue - bValue, + [&](auto i) { return (aValue - bValue) * i; }, + [&](auto i, auto j) { return (aValue - bValue) * i * j; }); + + auto d = a; + d -= b; + check( + d, dim, aValue - bValue, + [&](auto i) { return (aValue - bValue) * i; }, + [&](auto i, auto j) { return (aValue - bValue) * i * j; }); + } + }); +} + +BOOST_AUTO_TEST_CASE(ADValueMult) { + defaultTestSuite([]() { + T aValue = 42; + T bValue = 23; + + { + auto a = exampleValue(aValue); + auto b = exampleValue(bValue); + + auto c = a * b; + check( + c, dim, aValue * bValue, + [&](auto i) { return 2 * aValue * bValue * i; }, + [&](auto i, auto j) { return 4 * aValue * bValue * i * j; }); + + auto d = a; + d *= b; + check( + d, dim, aValue * bValue, + [&](auto i) { return 2 * aValue * bValue * i; }, + [&](auto i, auto j) { return 4 * aValue * bValue * i * j; }); + } + + { + auto a = exampleValue(dim, aValue); + auto b = exampleValue(dim, bValue); + + auto c = a * b; + check( + c, dim, aValue * bValue, + [&](auto i) { return 2 * aValue * bValue * i; }, + [&](auto i, auto j) { return 4 * aValue * bValue * i * j; }); + + auto d = a; + d *= b; + check( + d, dim, aValue * bValue, + [&](auto i) { return 2 * aValue * bValue * i; }, + [&](auto i, auto j) { return 4 * aValue * bValue * i * j; }); + } + }); +} + +//************************************************** +//* ToDo: Add checks for the following features +//* - Division of ADValue and ADValue +//* - Arithmetic operations of ADValue and raw value +//* - Nonlinear functions +//************************************************** + +BOOST_AUTO_TEST_SUITE_END() + +#else //__cplusplus >= 202002L + +BOOST_AUTO_TEST_SUITE(traceless_advalue) +BOOST_AUTO_TEST_CASE(ADValueNotTested) { + std::cout << "ADOL-C Warning: ADValue is not tested since test is not " + "compiled with C++20 support." + << std::endl; +} +BOOST_AUTO_TEST_SUITE_END() + +#endif //__cplusplus >= 202002L From 7c200771baf4e0c08d2990f898231979d0a6b867 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Thu, 23 Jan 2025 17:04:05 +0100 Subject: [PATCH 3/3] [advalue] Convert constants if necessary Unfortunately we cannot rely on the fact that number-like type `T` support arithmetic operations with `int` or `double`. E.g. `T = boost::numeric::interval` only supports arithmetic operations between `T` and `S`. Hence we cannot use `t-1`. While `t-1.` works for `S=double` it still fails for `S=float`. As a remedy we use check if constants of the raw type are compatible and convert the constant to `T` otherwise. Doing the conversion always might be costly for other types `T` that support cheap operations with elementary type (e.g. `T=ADValue` itself). --- ADOL-C/include/adolc/advalue.h | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/ADOL-C/include/adolc/advalue.h b/ADOL-C/include/adolc/advalue.h index def13cf4..68e8a7a3 100644 --- a/ADOL-C/include/adolc/advalue.h +++ b/ADOL-C/include/adolc/advalue.h @@ -1027,12 +1027,14 @@ auto adCompose(const Value &x, const Derivatives &f) { template X inv(const X &x) { return adCompose(x, [](auto k, auto x) { constexpr auto order = decltype(k)::value; + constexpr auto canUseInt = requires() { 1 / x; }; + using Constant = std::conditional_t; if constexpr (order == 0) - return 1 / x; + return Constant(1) / x; if constexpr (order == 1) - return -1 / (x * x); + return Constant(-1) / (x * x); if constexpr (order == 2) - return 2 / (x * x * x); + return Constant(2) / (x * x * x); static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); }); } @@ -1125,12 +1127,14 @@ template auto log(const X &x) { using std::pow; return adCompose(x, [](auto k, auto x) { constexpr auto order = decltype(k)::value; + constexpr auto canUseInt = requires() { 1 / x; }; + using Constant = std::conditional_t; if constexpr (order == 0) return log(x); if constexpr (order == 1) - return 1. / x; + return Constant(1) / x; if constexpr (order == 2) - return -1. / (x * x); + return Constant(-1) / (x * x); static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); }); } @@ -1143,12 +1147,14 @@ template auto sqrt(const X &x) { using std::sqrt; return adCompose(x, [](auto k, auto x) { constexpr auto order = decltype(k)::value; + constexpr auto canUseDouble = requires() { 1. / x; }; + using Constant = std::conditional_t; if constexpr (order == 0) return sqrt(x); if constexpr (order == 1) - return 1. / (2. * sqrt(x)); + return Constant(1. / 2.) / sqrt(x); if constexpr (order == 2) - return -1. / (4. * x * sqrt(x)); + return Constant(-1. / 4.) / (x * sqrt(x)); static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); }); } @@ -1162,12 +1168,14 @@ auto pow(const X &x, const Y &y) { using std::pow; return adCompose(x, [y](auto k, auto x) { constexpr auto order = decltype(k)::value; + constexpr auto canUseInt = requires() { y - 1; }; + using Constant = std::conditional_t; if constexpr (order == 0) return pow(x, y); if constexpr (order == 1) - return y * pow(x, y - 1); + return y * pow(x, y - Constant(1)); if constexpr (order == 2) - return y * (y - 1.) * pow(x, y - 2.); + return y * (y - Constant(1)) * pow(x, y - Constant(2)); static_assert(order <= 2, "Only derivatives up to order 2 are implemented"); }); }