diff --git a/unsupported/Eigen/AutoDiff b/unsupported/Eigen/AutoDiff index 45078bcd0..0480c69ee 100644 --- a/unsupported/Eigen/AutoDiff +++ b/unsupported/Eigen/AutoDiff @@ -33,6 +33,7 @@ namespace Eigen { #include "../../Eigen/src/Core/util/DisableStupidWarnings.h" // IWYU pragma: begin_exports +#include "src/AutoDiff/CoherentPadOp.h" #include "src/AutoDiff/AutoDiffScalar.h" // #include "src/AutoDiff/AutoDiffVector.h" #include "src/AutoDiff/AutoDiffJacobian.h" diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h index 70e74c2e5..c6ffa0d5d 100644 --- a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h @@ -17,21 +17,52 @@ namespace Eigen { namespace internal { -template -struct make_coherent_impl { - static void run(A&, B&) {} -}; - -// resize a to match b is a.size()==0, and conversely. -template -void make_coherent(const A& a, const B& b) { - make_coherent_impl::run(a.const_cast_derived(), b.const_cast_derived()); -} - template struct auto_diff_special_op; -} // end namespace internal +template +struct maybe_coherent_pad_helper { + static constexpr int SizeAtCompileTime = DerivativeType::SizeAtCompileTime == Dynamic || + OtherDerivativeType::SizeAtCompileTime == Dynamic + ? Dynamic + : int(DerivativeType::SizeAtCompileTime) > + int(OtherDerivativeType::SizeAtCompileTime) + ? DerivativeType::SizeAtCompileTime + : OtherDerivativeType::SizeAtCompileTime; + using type = CoherentPadOp; + static type pad(const DerivativeType& x, const OtherDerivativeType& y) { + // CoherentPadOp uses variable_if_dynamic. In this case, `SizeAtCompileTime` might + // by Dynamic, so we need to take the runtime maximum of x, y. + return CoherentPadOp(x, numext::maxi(x.size(), y.size())); + } +}; + +// Both are fixed-sized and equal, don't need to pad. +// Both are fixed-size and this is larger than other, don't need to pad. +template +struct maybe_coherent_pad_helper< + DerivativeType, OtherDerivativeType, + std::enable_if_t= OtherDerivativeType::SizeAtCompileTime && + DerivativeType::SizeAtCompileTime != Dynamic && + OtherDerivativeType::SizeAtCompileTime != Dynamic>> { + using type = const DerivativeType&; + static const DerivativeType& pad(const DerivativeType& x, const OtherDerivativeType& /*y*/) { return x; } +}; + +template +typename maybe_coherent_pad_helper::type MaybeCoherentPad( + const DerivativeType& x, const OtherDerivativeType& y) { + return maybe_coherent_pad_helper::pad(x, y); +} + +template +auto MakeCoherentCwiseBinaryOp(const LhsDerivativeType& x, const RhsDerivativeType& y, Op op = Op()) { + const auto& lhs = MaybeCoherentPad(x, y); + const auto& rhs = MaybeCoherentPad(y, x); + return CwiseBinaryOp, remove_all_t>(lhs, rhs, op); +} + +} // namespace internal template class AutoDiffScalar; @@ -214,13 +245,10 @@ class AutoDiffScalar } template - inline AutoDiffScalar< - CwiseBinaryOp, const DerType, const internal::remove_all_t>> - operator+(const AutoDiffScalar& other) const { - internal::make_coherent(m_derivatives, other.derivatives()); - return AutoDiffScalar< - CwiseBinaryOp, const DerType, const internal::remove_all_t>>( - m_value + other.value(), m_derivatives + other.derivatives()); + inline auto operator+(const AutoDiffScalar& other) const { + return MakeAutoDiffScalar( + m_value + other.value(), + internal::MakeCoherentCwiseBinaryOp>(m_derivatives, other.derivatives())); } template @@ -245,13 +273,10 @@ class AutoDiffScalar } template - inline AutoDiffScalar< - CwiseBinaryOp, const DerType, const internal::remove_all_t>> - operator-(const AutoDiffScalar& other) const { - internal::make_coherent(m_derivatives, other.derivatives()); - return AutoDiffScalar, const DerType, - const internal::remove_all_t>>( - m_value - other.value(), m_derivatives - other.derivatives()); + inline auto operator-(const AutoDiffScalar& other) const { + return MakeAutoDiffScalar(m_value - other.value(), + internal::MakeCoherentCwiseBinaryOp>( + m_derivatives, other.derivatives())); } template @@ -264,13 +289,11 @@ class AutoDiffScalar return AutoDiffScalar, const DerType>>(-m_value, -m_derivatives); } - inline AutoDiffScalar operator*( - const Scalar& other) const { + inline auto operator*(const Scalar& other) const { return MakeAutoDiffScalar(m_value * other, m_derivatives * other); } - friend inline AutoDiffScalar operator*( - const Scalar& other, const AutoDiffScalar& a) { + friend inline auto operator*(const Scalar& other, const AutoDiffScalar& a) { return MakeAutoDiffScalar(a.value() * other, a.derivatives() * other); } @@ -290,13 +313,11 @@ class AutoDiffScalar // a.derivatives() * other); // } - inline AutoDiffScalar operator/( - const Scalar& other) const { + inline auto operator/(const Scalar& other) const { return MakeAutoDiffScalar(m_value / other, (m_derivatives * (Scalar(1) / other))); } - friend inline AutoDiffScalar operator/( - const Scalar& other, const AutoDiffScalar& a) { + friend inline auto operator/(const Scalar& other, const AutoDiffScalar& a) { return MakeAutoDiffScalar(other / a.value(), a.derivatives() * (Scalar(-other) / (a.value() * a.value()))); } @@ -317,26 +338,18 @@ class AutoDiffScalar // } template - inline AutoDiffScalar EIGEN_COMMA const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE( - DerType, Scalar, product) EIGEN_COMMA const - EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(internal::remove_all_t, Scalar, product)>, - Scalar, product)> - operator/(const AutoDiffScalar& other) const { - internal::make_coherent(m_derivatives, other.derivatives()); + inline auto operator/(const AutoDiffScalar& other) const { return MakeAutoDiffScalar(m_value / other.value(), - ((m_derivatives * other.value()) - (other.derivatives() * m_value)) * + internal::MakeCoherentCwiseBinaryOp>( + m_derivatives * other.value(), (other.derivatives() * m_value)) * (Scalar(1) / (other.value() * other.value()))); } template - inline AutoDiffScalar, const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(DerType, Scalar, product), - const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(internal::remove_all_t, Scalar, product)>> - operator*(const AutoDiffScalar& other) const { - internal::make_coherent(m_derivatives, other.derivatives()); + inline auto operator*(const AutoDiffScalar& other) const { return MakeAutoDiffScalar(m_value * other.value(), - (m_derivatives * other.value()) + (other.derivatives() * m_value)); + internal::MakeCoherentCwiseBinaryOp>( + m_derivatives * other.value(), other.derivatives() * m_value)); } inline AutoDiffScalar& operator*=(const Scalar& other) { @@ -430,64 +443,6 @@ struct auto_diff_special_op { void operator+() const; }; -template -void make_coherent_expression(CwiseBinaryOp xpr, const RefType& ref) { - make_coherent(xpr.const_cast_derived().lhs(), ref); - make_coherent(xpr.const_cast_derived().rhs(), ref); -} - -template -void make_coherent_expression(const CwiseUnaryOp& xpr, const RefType& ref) { - make_coherent(xpr.nestedExpression().const_cast_derived(), ref); -} - -// needed for compilation only -template -void make_coherent_expression(const CwiseNullaryOp&, const RefType&) {} - -template -struct make_coherent_impl, B> { - typedef Matrix A; - static void run(A& a, B& b) { - if ((A_Rows == Dynamic || A_Cols == Dynamic) && (a.size() == 0)) { - a.resize(b.size()); - a.setZero(); - } else if (B::SizeAtCompileTime == Dynamic && a.size() != 0 && b.size() == 0) { - make_coherent_expression(b, a); - } - } -}; - -template -struct make_coherent_impl> { - typedef Matrix B; - static void run(A& a, B& b) { - if ((B_Rows == Dynamic || B_Cols == Dynamic) && (b.size() == 0)) { - b.resize(a.size()); - b.setZero(); - } else if (A::SizeAtCompileTime == Dynamic && b.size() != 0 && a.size() == 0) { - make_coherent_expression(a, b); - } - } -}; - -template -struct make_coherent_impl, - Matrix> { - typedef Matrix A; - typedef Matrix B; - static void run(A& a, B& b) { - if ((A_Rows == Dynamic || A_Cols == Dynamic) && (a.size() == 0)) { - a.resize(b.size()); - a.setZero(); - } else if ((B_Rows == Dynamic || B_Cols == Dynamic) && (b.size() == 0)) { - b.resize(a.size()); - b.setZero(); - } - } -}; - } // end namespace internal template @@ -518,10 +473,7 @@ struct ScalarBinaryOpTraits, B #define EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(FUNC, CODE) \ template \ - inline Eigen::AutoDiffScalar, \ - typename Eigen::internal::traits>::Scalar, product)> \ - FUNC(const Eigen::AutoDiffScalar& x) { \ + inline auto FUNC(const Eigen::AutoDiffScalar& x) { \ using namespace Eigen; \ typedef typename Eigen::internal::traits>::Scalar Scalar; \ EIGEN_UNUSED_VARIABLE(sizeof(Scalar)); \ @@ -602,10 +554,8 @@ EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(log, using std::log; x.derivatives() * (Scalar(1) / x.value()));) template -inline Eigen::AutoDiffScalar, typename internal::traits>::Scalar, product)> -pow(const Eigen::AutoDiffScalar& x, - const typename internal::traits>::Scalar& y) { +inline auto pow(const Eigen::AutoDiffScalar& x, + const typename internal::traits>::Scalar& y) { using namespace Eigen; using std::pow; return Eigen::MakeAutoDiffScalar(pow(x.value(), y), x.derivatives() * (y * pow(x.value(), y - 1))); diff --git a/unsupported/Eigen/src/AutoDiff/CoherentPadOp.h b/unsupported/Eigen/src/AutoDiff/CoherentPadOp.h new file mode 100644 index 000000000..0fdd696c8 --- /dev/null +++ b/unsupported/Eigen/src/AutoDiff/CoherentPadOp.h @@ -0,0 +1,152 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2020 The Eigen Team. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_COHERENT_PAD_OP_H +#define EIGEN_COHERENT_PAD_OP_H + +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +namespace internal { + +// Pads a vector with zeros to a given size. +template +struct CoherentPadOp; + +template +struct traits> : public traits { + typedef typename internal::remove_all::type PlainXprType; + typedef typename internal::ref_selector::type XprNested; + typedef typename std::remove_reference_t XprNested_; + enum : int { + IsRowMajor = traits::Flags & RowMajorBit, + SizeAtCompileTime = SizeAtCompileTime_, + RowsAtCompileTime = IsRowMajor ? 1 : SizeAtCompileTime, + ColsAtCompileTime = IsRowMajor ? SizeAtCompileTime : 1, + MaxRowsAtCompileTime = RowsAtCompileTime, + MaxColsAtCompileTime = ColsAtCompileTime, + Flags = traits::Flags & ~NestByRefBit, + }; +}; + +// Pads a vector with zeros to a given size. +template +struct CoherentPadOp : public dense_xpr_base>::type { + typedef typename internal::generic_xpr_base>::type Base; + EIGEN_GENERIC_PUBLIC_INTERFACE(CoherentPadOp) + + using XprNested = typename traits::XprNested; + using XprNested_ = typename traits::XprNested_; + using NestedExpression = XprNested_; + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoherentPadOp() = delete; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoherentPadOp(const CoherentPadOp&) = default; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoherentPadOp(CoherentPadOp&& other) = default; + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoherentPadOp(const XprType& xpr, Index size) : xpr_(xpr), size_(size) { + static_assert(XprNested_::IsVectorAtCompileTime && "input type must be a vector"); + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const XprNested_& nestedExpression() const { return xpr_; } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index size() const { return size_.value(); } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index rows() const { + return traits::IsRowMajor ? Index(1) : size(); + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index cols() const { + return traits::IsRowMajor ? size() : Index(1); + } + + private: + XprNested xpr_; + const internal::variable_if_dynamic size_; +}; + +// Adapted from the Replicate evaluator. +template +struct unary_evaluator> + : evaluator_base> { + typedef CoherentPadOp XprType; + typedef typename internal::remove_all_t CoeffReturnType; + typedef typename internal::nested_eval::type ArgTypeNested; + typedef internal::remove_all_t ArgTypeNestedCleaned; + + enum { + CoeffReadCost = evaluator::CoeffReadCost, + LinearAccessMask = XprType::IsVectorAtCompileTime ? LinearAccessBit : 0, + Flags = evaluator::Flags & (HereditaryBits | LinearAccessMask | RowMajorBit), + Alignment = evaluator::Alignment + }; + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit unary_evaluator(const XprType& pad) + : m_arg(pad.nestedExpression()), m_argImpl(m_arg), m_size(pad.nestedExpression().size()) {} + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const { + EIGEN_IF_CONSTEXPR(XprType::IsRowMajor) { + if (col < m_size.value()) { + return m_argImpl.coeff(1, col); + } + } + else { + if (row < m_size.value()) { + return m_argImpl.coeff(row, 1); + } + } + return CoeffReturnType(0); + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const { + if (index < m_size.value()) { + return m_argImpl.coeff(index); + } + return CoeffReturnType(0); + } + + template + EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const { + // AutoDiff scalar's derivative must be a vector, which is enforced by static assert. + // Defer to linear access for simplicity. + EIGEN_IF_CONSTEXPR(XprType::IsRowMajor) { return packet(col); } + return packet(row); + } + + template + EIGEN_STRONG_INLINE PacketType packet(Index index) const { + constexpr int kPacketSize = unpacket_traits::size; + if (index + kPacketSize <= m_size.value()) { + return m_argImpl.template packet(index); + } else if (index < m_size.value()) { + // Partial packet. + EIGEN_ALIGN_MAX std::remove_const_t values[kPacketSize]; + const int partial = m_size.value() - index; + for (int i = 0; i < partial && i < kPacketSize; ++i) { + values[i] = m_argImpl.coeff(index + i); + } + for (int i = partial; i < kPacketSize; ++i) { + values[i] = CoeffReturnType(0); + } + return pload(values); + } + return pset1(CoeffReturnType(0)); + } + + protected: + const ArgTypeNested m_arg; + evaluator m_argImpl; + const variable_if_dynamic m_size; +}; + +} // namespace internal + +} // namespace Eigen + +#endif // EIGEN_CWISE_BINARY_OP_H diff --git a/unsupported/test/autodiff.cpp b/unsupported/test/autodiff.cpp index 9ea4e7dd4..a885a1e37 100644 --- a/unsupported/test/autodiff.cpp +++ b/unsupported/test/autodiff.cpp @@ -175,15 +175,11 @@ void forward_jacobian(const Func& f) { jref.setZero(); yref.setZero(); f(x, &yref, &jref); - // std::cerr << y.transpose() << "\n\n";; - // std::cerr << j << "\n\n";; j.setZero(); y.setZero(); AutoDiffJacobian autoj(f); autoj(x, &y, &j); - // std::cerr << y.transpose() << "\n\n";; - // std::cerr << j << "\n\n";; VERIFY_IS_APPROX(y, yref); VERIFY_IS_APPROX(j, jref); @@ -277,8 +273,6 @@ double bug_1222() { return denom.value(); } -#ifdef EIGEN_TEST_PART_5 - double bug_1223() { using std::min; typedef Eigen::AutoDiffScalar AD; @@ -338,8 +332,6 @@ double bug_1281() { return (y1 + y2 + y3).value(); } -#endif - EIGEN_DECLARE_TEST(autodiff) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(test_autodiff_scalar<1>());