Add complex<dual<T,N>> Eigen support, fix rpart bug (cdc386c4) · Commits · Michael Tesch / cppduals

CMakeLists.txt

+1 −1

Original line number	Diff line number	Diff line
		@@ -3,7 +3,7 @@
		#
		cmake_minimum_required (VERSION 3.14)
		project (cppduals
		VERSION 0.8.1
		VERSION 0.8.2
		LANGUAGES C CXX
		)
		include (GNUInstallDirs)

README.md

+7 −0

Original line number	Diff line number	Diff line
		@@ -258,6 +258,13 @@ also licensed under MPL-2.0.
		ChangeLog
		=========

		v0.8.2
		------

		- add `complex<dual<T,N>>` matrix-level `rpart`/`dpart`/`dconj` to `multidual_eigen`.
		- fix `rpart(complex<dual<T,N>>)` to return `complex<T>` (was incorrectly returning `dual<T,N>`).
		- add Eigen integration tests for `dual<T,N>` and `complex<dual<T,N>>` matrices.

		v0.8.1
		------

duals/multidual

+3 −3

Original line number	Diff line number	Diff line
		@@ -367,9 +367,9 @@ public:
		template <class T, int N> T rpart(const dual<T,N> & x) { return x.rpart(); }
		template <class T, int N> T dpart(const dual<T,N> & x, int i = 0) { return x.dpart(i); }

		/// Real part of complex<dual<T,N>>
		template <class T, int N> dual<T,N> rpart(const std::complex<dual<T,N>> & x)
		{ return x.real(); }
		/// R-part of complex<dual<T,N>> is non-dual complex<T> (not to be confused with real())
		template <class T, int N> std::complex<T> rpart(const std::complex<dual<T,N>> & x)
		{ return std::complex<T>(x.real().rpart(), x.imag().rpart()); }

		/// Dual part of complex<dual<T,N>> — returns complex of the i-th partial
		template <class T, int N> std::complex<T> dpart(const std::complex<dual<T,N>> & x, int i = 0)

duals/multidual_eigen

+92 −6

Original line number	Diff line number	Diff line
		@@ -37,6 +37,11 @@

		namespace duals {

		// Trait to detect complex<dual<T,N>> scalars.
		template<typename S> struct is_complex_dual : std::false_type {};
		template<typename T, int N>
		struct is_complex_dual<std::complex<dual<T,N>>> : std::true_type {};

		/// Unary functor: extract rpart from each element of a dual<T,N> matrix.
		template<typename T, int N>
		struct CwiseMDRpartOp {
		@@ -63,11 +68,50 @@ struct CwiseMDDconjOp {
		EIGEN_STRONG_INLINE dual<T,N> operator()(const dual<T,N> & x) const { return dconj(x); }
		};

		// --- complex<dual<T,N>> functors ---

		/// Unary functor: rpart of complex<dual<T,N>> → complex<T>
		template<typename T, int N>
		struct CwiseMDCRpartOp {
		typedef std::complex<T> result_type;
		EIGEN_DEVICE_FUNC
		EIGEN_STRONG_INLINE std::complex<T> operator()(const std::complex<dual<T,N>> & x) const {
		return rpart(x);
		}
		};

		/// Unary functor: dpart(i) of complex<dual<T,N>> → complex<T>
		template<typename T, int N>
		struct CwiseMDCDpartOp {
		int idx;
		typedef std::complex<T> result_type;
		CwiseMDCDpartOp(int i = 0) : idx(i) {}
		EIGEN_DEVICE_FUNC
		EIGEN_STRONG_INLINE std::complex<T> operator()(const std::complex<dual<T,N>> & x) const {
		return dpart(x, idx);
		}
		};

		/// Unary functor: dual-conjugate each element of a complex<dual<T,N>> matrix.
		template<typename T, int N>
		struct CwiseMDCDconjOp {
		typedef std::complex<dual<T,N>> result_type;
		EIGEN_DEVICE_FUNC
		EIGEN_STRONG_INLINE std::complex<dual<T,N>> operator()(const std::complex<dual<T,N>> & x) const {
		return std::complex<dual<T,N>>(dconj(x.real()), dconj(x.imag()));
		}
		};

		// ===================================================================
		// Matrix-level rpart / dpart / dconj — dual<T,N> scalars
		// ===================================================================

		/// Extract the "real part" of a dual<T,N>-valued matrix.
		template <typename Derived>
		auto rpart(const Eigen::EigenBase<Derived> & x)
		-> decltype(x.derived().unaryExpr(CwiseMDRpartOp<
		typename Derived::Scalar::value_type, Derived::Scalar::num_vars>()))
		-> std::enable_if_t<is_dual<typename Derived::Scalar>::value,
		decltype(x.derived().unaryExpr(CwiseMDRpartOp<
		typename Derived::Scalar::value_type, Derived::Scalar::num_vars>()))>
		{
		using S = typename Derived::Scalar;
		return x.derived().unaryExpr(CwiseMDRpartOp<typename S::value_type, S::num_vars>());
		@@ -76,8 +120,9 @@ auto rpart(const Eigen::EigenBase<Derived> & x)
		/// Extract the i-th dual part of a dual<T,N>-valued matrix.
		template <typename Derived>
		auto dpart(const Eigen::EigenBase<Derived> & x, int i = 0)
		-> decltype(x.derived().unaryExpr(CwiseMDDpartOp<
		typename Derived::Scalar::value_type, Derived::Scalar::num_vars>(i)))
		-> std::enable_if_t<is_dual<typename Derived::Scalar>::value,
		decltype(x.derived().unaryExpr(CwiseMDDpartOp<
		typename Derived::Scalar::value_type, Derived::Scalar::num_vars>(i)))>
		{
		using S = typename Derived::Scalar;
		return x.derived().unaryExpr(CwiseMDDpartOp<typename S::value_type, S::num_vars>(i));
		@@ -86,13 +131,54 @@ auto dpart(const Eigen::EigenBase<Derived> & x, int i = 0)
		/// Dual-conjugate a dual<T,N>-valued matrix.
		template <typename Derived>
		auto dconj(const Eigen::EigenBase<Derived> & x)
		-> decltype(x.derived().unaryExpr(CwiseMDDconjOp<
		typename Derived::Scalar::value_type, Derived::Scalar::num_vars>()))
		-> std::enable_if_t<is_dual<typename Derived::Scalar>::value,
		decltype(x.derived().unaryExpr(CwiseMDDconjOp<
		typename Derived::Scalar::value_type, Derived::Scalar::num_vars>()))>
		{
		using S = typename Derived::Scalar;
		return x.derived().unaryExpr(CwiseMDDconjOp<typename S::value_type, S::num_vars>());
		}

		// ===================================================================
		// Matrix-level rpart / dpart / dconj — complex<dual<T,N>> scalars
		// ===================================================================

		/// Extract the non-dual part of a complex<dual<T,N>>-valued matrix → Matrix<complex<T>>.
		template <typename Derived>
		auto rpart(const Eigen::EigenBase<Derived> & x)
		-> std::enable_if_t<is_complex_dual<typename Derived::Scalar>::value,
		decltype(x.derived().unaryExpr(CwiseMDCRpartOp<
		typename Derived::Scalar::value_type::value_type,
		Derived::Scalar::value_type::num_vars>()))>
		{
		using D = typename Derived::Scalar::value_type; // dual<T,N>
		return x.derived().unaryExpr(CwiseMDCRpartOp<typename D::value_type, D::num_vars>());
		}

		/// Extract the i-th dual part of a complex<dual<T,N>>-valued matrix → Matrix<complex<T>>.
		template <typename Derived>
		auto dpart(const Eigen::EigenBase<Derived> & x, int i = 0)
		-> std::enable_if_t<is_complex_dual<typename Derived::Scalar>::value,
		decltype(x.derived().unaryExpr(CwiseMDCDpartOp<
		typename Derived::Scalar::value_type::value_type,
		Derived::Scalar::value_type::num_vars>(i)))>
		{
		using D = typename Derived::Scalar::value_type; // dual<T,N>
		return x.derived().unaryExpr(CwiseMDCDpartOp<typename D::value_type, D::num_vars>(i));
		}

		/// Dual-conjugate a complex<dual<T,N>>-valued matrix.
		template <typename Derived>
		auto dconj(const Eigen::EigenBase<Derived> & x)
		-> std::enable_if_t<is_complex_dual<typename Derived::Scalar>::value,
		decltype(x.derived().unaryExpr(CwiseMDCDconjOp<
		typename Derived::Scalar::value_type::value_type,
		Derived::Scalar::value_type::num_vars>()))>
		{
		using D = typename Derived::Scalar::value_type; // dual<T,N>
		return x.derived().unaryExpr(CwiseMDCDconjOp<typename D::value_type, D::num_vars>());
		}

		} // namespace duals

		namespace Eigen {

tests/test_multidual.cpp

+85 −0

Original line number	Diff line number	Diff line
		@@ -5,6 +5,8 @@

		#include "gtest/gtest.h"
		#include <duals/multidual>
		#include <duals/multidual_eigen>
		#include <Eigen/Dense>
		#include <complex>

		using duals::dual;
		@@ -462,3 +464,86 @@ TEST(multidual_gradient, x_sin_y) {
		EXPECT_NEAR(f.dpart(0), sin(M_PI/3), 1e-12); // df/dx = sin(y)
		EXPECT_NEAR(f.dpart(1), 2 * cos(M_PI/3), 1e-12); // df/dy = x*cos(y)
		}

		// ===================================================================
		// Eigen integration: dual<T,N> and complex<dual<T,N>> matrices
		// ===================================================================

		TEST(multidual_eigen, matrix_rpart_dpart) {
		using D3 = dual<double, 3>;
		Eigen::Matrix<D3, 2, 2> M;
		M(0,0) = D3(1.0, {2.0, 3.0, 4.0});
		M(0,1) = D3(5.0, {6.0, 7.0, 8.0});
		M(1,0) = D3(9.0, {10.0, 11.0, 12.0});
		M(1,1) = D3(13.0, {14.0, 15.0, 16.0});

		Eigen::Matrix<double, 2, 2> R = duals::rpart(M);
		EXPECT_DOUBLE_EQ(R(0,0), 1.0);
		EXPECT_DOUBLE_EQ(R(0,1), 5.0);
		EXPECT_DOUBLE_EQ(R(1,0), 9.0);
		EXPECT_DOUBLE_EQ(R(1,1), 13.0);

		Eigen::Matrix<double, 2, 2> D0 = duals::dpart(M, 0);
		EXPECT_DOUBLE_EQ(D0(0,0), 2.0);
		EXPECT_DOUBLE_EQ(D0(1,1), 14.0);

		Eigen::Matrix<double, 2, 2> D2 = duals::dpart(M, 2);
		EXPECT_DOUBLE_EQ(D2(0,0), 4.0);
		EXPECT_DOUBLE_EQ(D2(0,1), 8.0);
		}

		TEST(multidual_eigen, complex_dual_rpart_dpart) {
		using D2 = dual<double, 2>;
		using CD2 = std::complex<D2>;
		Eigen::Matrix<CD2, 2, 1> v;
		v(0) = CD2(D2(1.0, {2.0, 3.0}), D2(4.0, {5.0, 6.0}));
		v(1) = CD2(D2(7.0, {8.0, 9.0}), D2(10.0, {11.0, 12.0}));

		// rpart strips duals: complex<dual> → complex<double>
		Eigen::Matrix<complexd, 2, 1> R = duals::rpart(v);
		EXPECT_DOUBLE_EQ(R(0).real(), 1.0);
		EXPECT_DOUBLE_EQ(R(0).imag(), 4.0);
		EXPECT_DOUBLE_EQ(R(1).real(), 7.0);
		EXPECT_DOUBLE_EQ(R(1).imag(), 10.0);

		// dpart(i) extracts i-th partial as complex<double>
		Eigen::Matrix<complexd, 2, 1> D0 = duals::dpart(v, 0);
		EXPECT_DOUBLE_EQ(D0(0).real(), 2.0);
		EXPECT_DOUBLE_EQ(D0(0).imag(), 5.0);

		Eigen::Matrix<complexd, 2, 1> D1 = duals::dpart(v, 1);
		EXPECT_DOUBLE_EQ(D1(0).real(), 3.0);
		EXPECT_DOUBLE_EQ(D1(0).imag(), 6.0);
		EXPECT_DOUBLE_EQ(D1(1).real(), 9.0);
		EXPECT_DOUBLE_EQ(D1(1).imag(), 12.0);
		}

		TEST(multidual_eigen, scalar_promotion) {
		using D2 = dual<double, 2>;
		Eigen::Matrix<D2, 2, 2> M;
		M(0,0) = D2(1.0, {1.0, 0.0});
		M(0,1) = D2(0.0, {0.0, 0.0});
		M(1,0) = D2(0.0, {0.0, 0.0});
		M(1,1) = D2(2.0, {0.0, 1.0});

		// scalar * matrix should promote
		auto R = 3.0 * M;
		EXPECT_DOUBLE_EQ(R(0,0).rpart(), 3.0);
		EXPECT_DOUBLE_EQ(R(0,0).dpart(0), 3.0);
		EXPECT_DOUBLE_EQ(R(1,1).rpart(), 6.0);
		EXPECT_DOUBLE_EQ(R(1,1).dpart(1), 3.0);
		}

		TEST(multidual_eigen, matrix_multiply) {
		using D2 = dual<double, 2>;
		Eigen::Matrix<D2, 2, 2> A, B;
		A(0,0) = D2::variable(1.0, 0); A(0,1) = D2(0.0);
		A(1,0) = D2(0.0); A(1,1) = D2::variable(2.0, 1);
		B.setIdentity();

		auto C = A * B;
		EXPECT_DOUBLE_EQ(C(0,0).rpart(), 1.0);
		EXPECT_DOUBLE_EQ(C(0,0).dpart(0), 1.0);
		EXPECT_DOUBLE_EQ(C(1,1).rpart(), 2.0);
		EXPECT_DOUBLE_EQ(C(1,1).dpart(1), 1.0);
		}