holder.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUN_HOLDER_HPP
#define STAN_MATH_PRIM_FUN_HOLDER_HPP
 
#include <stan/math/prim/meta/is_eigen.hpp>
#include <stan/math/prim/meta/is_plain_type.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <memory>
#include <type_traits>
#include <utility>
 
// This was implenmented following the tutorial on edding new expressions to
// Eigen: https://eigen.tuxfamily.org/dox/TopicNewExpressionType.html
 
namespace stan {
namespace math {
 
template <class ArgType, typename... Ptrs>
class Holder;
 
}  // namespace math
}  // namespace stan
 
namespace Eigen {
namespace internal {
 
template <class ArgType, typename... Ptrs>
struct traits<stan::math::Holder<ArgType, Ptrs...>> {
  typedef typename ArgType::StorageKind StorageKind;
  typedef typename traits<ArgType>::XprKind XprKind;
  typedef typename ArgType::StorageIndex StorageIndex;
  typedef typename ArgType::Scalar Scalar;
  enum {
    // Possible flags are documented here:
    // https://eigen.tuxfamily.org/dox/group__flags.html
    Flags = (ArgType::Flags
             & (RowMajorBit | LvalueBit | LinearAccessBit | DirectAccessBit
                | PacketAccessBit | NoPreferredStorageOrderBit))
            | NestByRefBit,
    RowsAtCompileTime = ArgType::RowsAtCompileTime,
    ColsAtCompileTime = ArgType::ColsAtCompileTime,
    MaxRowsAtCompileTime = ArgType::MaxRowsAtCompileTime,
    MaxColsAtCompileTime = ArgType::MaxColsAtCompileTime,
    InnerStrideAtCompileTime = ArgType::InnerStrideAtCompileTime,
    OuterStrideAtCompileTime = ArgType::OuterStrideAtCompileTime
  };
};
 
}  // namespace internal
}  // namespace Eigen
 
namespace stan {
namespace math {
 
template <typename ArgType, typename... Ptrs>
class Holder
    : public Eigen::internal::dense_xpr_base<Holder<ArgType, Ptrs...>>::type {
 public:
  typedef typename Eigen::internal::ref_selector<Holder<ArgType, Ptrs...>>::type
      Nested;
  typename Eigen::internal::ref_selector<ArgType>::non_const_type m_arg;
  std::tuple<std::unique_ptr<Ptrs>...> m_unique_ptrs;
 
  explicit Holder(ArgType&& arg, Ptrs*... pointers)
      : m_arg(arg), m_unique_ptrs(std::unique_ptr<Ptrs>(pointers)...) {}
 
  // we need to explicitely default copy and move constructors as we are
  // defining copy and move assignment operators
  Holder(const Holder<ArgType, Ptrs...>&) = default;
  Holder(Holder<ArgType, Ptrs...>&&) = default;
 
  // all these functions just call the same on the argument
  Eigen::Index rows() const { return m_arg.rows(); }
  Eigen::Index cols() const { return m_arg.cols(); }
  Eigen::Index innerStride() const { return m_arg.innerStride(); }
  Eigen::Index outerStride() const { return m_arg.outerStride(); }
  auto* data() { return m_arg.data(); }
 
  template <typename T, require_eigen_t<T>* = nullptr>
  inline Holder<ArgType, Ptrs...>& operator=(const T& other) {
    m_arg = other;
    return *this;
  }
 
  // copy and move assignment operators need to be separately overloaded,
  // otherwise defaults will be used.
  inline Holder<ArgType, Ptrs...>& operator=(
      const Holder<ArgType, Ptrs...>& other) {
    m_arg = other;
    return *this;
  }
  inline Holder<ArgType, Ptrs...>& operator=(Holder<ArgType, Ptrs...>&& other) {
    m_arg = std::move(other.m_arg);
    return *this;
  }
};
 
}  // namespace math
}  // namespace stan
 
namespace Eigen {
namespace internal {
 
template <typename ArgType, typename... Ptrs>
struct evaluator<stan::math::Holder<ArgType, Ptrs...>>
    : evaluator_base<stan::math::Holder<ArgType, Ptrs...>> {
  typedef stan::math::Holder<ArgType, Ptrs...> XprType;
  typedef typename remove_all<ArgType>::type ArgTypeNestedCleaned;
  typedef typename XprType::CoeffReturnType CoeffReturnType;
  typedef typename XprType::Scalar Scalar;
  enum {
    CoeffReadCost = evaluator<ArgTypeNestedCleaned>::CoeffReadCost,
    // Possible flags are documented here:
    // https://eigen.tuxfamily.org/dox/group__flags.html
    Flags = evaluator<ArgTypeNestedCleaned>::Flags,
    Alignment = evaluator<ArgTypeNestedCleaned>::Alignment,
  };
 
  evaluator<ArgTypeNestedCleaned> m_argImpl;
 
  explicit evaluator(const XprType& xpr) : m_argImpl(xpr.m_arg) {}
 
  // all these functions just call the same on the argument
  EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const {
    return m_argImpl.coeff(row, col);
  }
  EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const {
    return m_argImpl.coeff(index);
  }
 
  EIGEN_STRONG_INLINE Scalar& coeffRef(Index row, Index col) {
    return m_argImpl.coeffRef(row, col);
  }
  EIGEN_STRONG_INLINE Scalar& coeffRef(Index index) {
    return m_argImpl.coeffRef(index);
  }
 
  template <int LoadMode, typename PacketType>
  EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const {
    return m_argImpl.template packet<LoadMode, PacketType>(row, col);
  }
  template <int LoadMode, typename PacketType>
  EIGEN_STRONG_INLINE PacketType packet(Index index) const {
    return m_argImpl.template packet<LoadMode, PacketType>(index);
  }
 
  template <int StoreMode, typename PacketType>
  EIGEN_STRONG_INLINE void writePacket(Index row, Index col,
                                       const PacketType& x) {
    return m_argImpl.template writePacket<StoreMode, PacketType>(row, col, x);
  }
  template <int StoreMode, typename PacketType>
  EIGEN_STRONG_INLINE void writePacket(Index index, const PacketType& x) {
    return m_argImpl.template writePacket<StoreMode, PacketType>(index, x);
  }
};
 
}  // namespace internal
}  // namespace Eigen
 
namespace stan {
namespace math {
 
template <typename T, typename... Ptrs,
          std::enable_if_t<sizeof...(Ptrs) >= 1>* = nullptr>
Holder<T, Ptrs...> holder(T&& arg, Ptrs*... pointers) {
  return Holder<T, Ptrs...>(std::forward<T>(arg), pointers...);
}
// trivial case with no pointers constructs no holder object
template <typename T>
T holder(T&& arg) {
  return std::forward<T>(arg);
}
 
namespace internal {
// the function holder_handle_element is also used in holder_cl
template <typename T>
auto holder_handle_element(T& a, T*& res) {
  res = &a;
  return std::make_tuple();
}
template <typename T,
          std::enable_if_t<!(Eigen::internal::traits<std::decay_t<T>>::Flags
                             & Eigen::NestByRefBit)>* = nullptr>
auto holder_handle_element(T&& a, std::remove_reference_t<T>*& res) {
  res = &a;
  return std::make_tuple();
}
 
template <typename T, require_t<std::is_rvalue_reference<T&&>>* = nullptr,
          std::enable_if_t<
              static_cast<bool>(Eigen::internal::traits<std::decay_t<T>>::Flags&
                                    Eigen::NestByRefBit)>* = nullptr>
auto holder_handle_element(T&& a, T*& res) {
  res = new T(std::move(a));
  return std::make_tuple(res);
}
template <typename T, require_t<std::is_rvalue_reference<T&&>>* = nullptr,
          require_not_eigen_t<T>* = nullptr>
auto holder_handle_element(T&& a, T*& res) {
  res = new T(std::move(a));
  return std::make_tuple(res);
}
 
template <typename T, std::size_t... Is, typename... Args>
auto make_holder_impl_construct_object(T&& expr, std::index_sequence<Is...>,
                                       const std::tuple<Args*...>& ptrs) {
  return holder(std::forward<T>(expr), std::get<Is>(ptrs)...);
}
 
template <typename F, std::size_t... Is, typename... Args>
auto make_holder_impl(const F& func, std::index_sequence<Is...>,
                      Args&&... args) {
  std::tuple<std::remove_reference_t<Args>*...> res;
  auto ptrs = std::tuple_cat(
      holder_handle_element(std::forward<Args>(args), std::get<Is>(res))...);
  return make_holder_impl_construct_object(
      func(*std::get<Is>(res)...),
      std::make_index_sequence<std::tuple_size<decltype(ptrs)>::value>(), ptrs);
}
 
}  // namespace internal
 
template <typename F, typename... Args,
          require_not_plain_type_t<
              decltype(std::declval<F>()(std::declval<Args&>()...))>* = nullptr>
auto make_holder(const F& func, Args&&... args) {
  return internal::make_holder_impl(func,
                                    std::make_index_sequence<sizeof...(Args)>(),
                                    std::forward<Args>(args)...);
}
 
template <typename F, typename... Args,
          require_plain_type_t<
              decltype(std::declval<F>()(std::declval<Args&>()...))>* = nullptr>
auto make_holder(const F& func, Args&&... args) {
  return func(std::forward<Args>(args)...);
}
 
}  // namespace math
}  // namespace stan
 
#endif