From 95a4f02ce1df9db434026fc807ea2ac748147319 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 30 Mar 2026 12:37:17 +0200 Subject: [PATCH 1/7] Add dR_drotvec and Python bindings Implement derivative of rotation matrix w.r.t. rotation vector (Rodrigues differentiation) and expose it to Python. C++: add Rotation::dR_drotvec(), dR_drotvec(const vec::fixed<3>&) implementing the Gallego & Yezzi (2015) formula with a small-angle limit and an overload that validates generic vec sizes. Returns an arma::cube (3x3x3) where slice k = dR/d(omega_k). Python wrappers: add a Rotation.dR_drotvec() method and a free dR_drotvec(rotvec) function that convert the arma::cube into a (3,3,3) NumPy array (result[k] = dR/d(omega_k)) and validate input size for the free function. --- .../Libraries/Maths/rotation.cpp | 31 ++++++++ src/Simulation/Maths/rotation.cpp | 72 +++++++++++++++++++ 2 files changed, 103 insertions(+) diff --git a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp index 33d3bdae..e7e2260f 100644 --- a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp +++ b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp @@ -139,7 +139,38 @@ void register_rotation(py::module_& m) { py::arg("B"), py::arg("active") = true, simcoon_docs::apply_stress_concentration) + .def("dR_drotvec", + [](const simcoon::Rotation& self) { + arma::cube c = self.dR_drotvec(); + // Return as (3, 3, 3) numpy array: result[k] = dR/d(omega_k) + py::array_t result({3, 3, 3}); + auto buf = result.mutable_unchecked<3>(); + for (int k = 0; k < 3; ++k) + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + buf(k, i, j) = c(i, j, k); + return result; + }, + simcoon_docs::dR_drotvec) + ; + + // Free function — takes a rotvec directly + m.def("dR_drotvec", + [](py::array_t rotvec) { + validate_vector_size(rotvec, 3, "rotvec"); + vec omega = carma::arr_to_col(rotvec); + arma::cube c = simcoon::dR_drotvec(omega); + py::array_t result({3, 3, 3}); + auto buf = result.mutable_unchecked<3>(); + for (int k = 0; k < 3; ++k) + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + buf(k, i, j) = c(i, j, k); + return result; + }, + py::arg("rotvec"), + simcoon_docs::dR_drotvec_free); } } // namespace simpy diff --git a/src/Simulation/Maths/rotation.cpp b/src/Simulation/Maths/rotation.cpp index c349efed..d17c5f01 100755 --- a/src/Simulation/Maths/rotation.cpp +++ b/src/Simulation/Maths/rotation.cpp @@ -972,4 +972,76 @@ mat rotate_stress_concentration(const mat &B, const mat &DR, const bool &active) return vs*(B*trans(ve)); } +// ============================================================================= +// dR/d(rotvec) — Rodrigues differentiation (Gallego & Yezzi 2015) +// ============================================================================= + +cube Rotation::dR_drotvec() const { + return simcoon::dR_drotvec(as_rotvec()); +} + +cube dR_drotvec(const vec::fixed<3>& omega) { + double theta = norm(omega); + cube result(3, 3, 3); + + // Skew-symmetric matrix of omega + mat::fixed<3,3> W; + W.zeros(); + W(0,1) = -omega(2); W(0,2) = omega(1); + W(1,0) = omega(2); W(1,2) = -omega(0); + W(2,0) = -omega(1); W(2,1) = omega(0); + + if (theta < simcoon::iota) { + // Small angle limit: dR/d(omega_k) = [e_k]x + result.zeros(); + result(1,2,0) = -1.0; result(2,1,0) = 1.0; // [e0]x + result(0,2,1) = 1.0; result(2,0,1) = -1.0; // [e1]x + result(0,1,2) = -1.0; result(1,0,2) = 1.0; // [e2]x + return result; + } + + double s = sin(theta), c = cos(theta); + double t2 = theta * theta; + + // Rodrigues coefficients and their derivatives w.r.t. theta + double a = s / theta; // sin(theta)/theta + double b = (1.0 - c) / t2; // (1-cos(theta))/theta^2 + double da = (c * theta - s) / t2; // d(sin(theta)/theta)/d(theta) + double db = (s * theta - 2.0 * (1.0 - c)) / (t2 * theta); // d((1-cos)/theta^2)/d(theta) + + mat::fixed<3,3> W2 = W * W; + + for (int k = 0; k < 3; ++k) { + // dW/d(omega_k) = [e_k]x + mat::fixed<3,3> dW; + dW.zeros(); + // Using direct assignment for the skew of unit vector e_k + if (k == 0) { + dW(1,2) = -1.0; dW(2,1) = 1.0; + } else if (k == 1) { + dW(0,2) = 1.0; dW(2,0) = -1.0; + } else { + dW(0,1) = -1.0; dW(1,0) = 1.0; + } + + // dW^2/d(omega_k) = dW*W + W*dW + mat::fixed<3,3> dW2 = dW * W + W * dW; + + // d(theta)/d(omega_k) = omega_k / theta + double dtk = omega(k) / theta; + + result.slice(k) = da * dtk * W + a * dW + db * dtk * W2 + b * dW2; + } + + return result; +} + +cube dR_drotvec(const vec& omega) { + if (omega.n_elem != 3) { + throw invalid_argument("dR_drotvec: rotation vector must have 3 elements"); + } + vec::fixed<3> omega_f = {omega(0), omega(1), omega(2)}; + return dR_drotvec(omega_f); +} + } //namespace simcoon From 34757bfd36e3c686499bbf45108e61178d51e112 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 30 Mar 2026 12:37:46 +0200 Subject: [PATCH 2/7] Declare dR_drotvec derivatives for rotation Add declarations and documentation for computing derivatives of the rotation matrix R(omega) with respect to rotation vector components. Introduces a member function Rotation::dR_drotvec() and two free-function overloads (fixed-size arma::vec::fixed<3> and dynamic arma::vec) that return a 3x3x3 arma::cube with slice k = dR/d(omega_k). Comments reference exact differentiation of the Rodrigues formula (Gallego & Yezzi, 2015) and clarify return shape and purpose. --- include/simcoon/Simulation/Maths/rotation.hpp | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/include/simcoon/Simulation/Maths/rotation.hpp b/include/simcoon/Simulation/Maths/rotation.hpp index 02274599..22077c1f 100755 --- a/include/simcoon/Simulation/Maths/rotation.hpp +++ b/include/simcoon/Simulation/Maths/rotation.hpp @@ -450,8 +450,41 @@ class Rotation { // Identity quaternion is [0, 0, 0, 1] or [0, 0, 0, -1] return std::abs(std::abs(_quat(3)) - 1.0) < tol; } + + /** + * @brief Compute derivatives of the rotation matrix w.r.t. rotation vector components. + * + * Uses the exact differentiation of the Rodrigues formula + * (Gallego & Yezzi, J. Math. Imaging Vis., 2015). + * + * @return 3x3x3 cube where slice k is dR/d(omega_k) + */ + arma::cube dR_drotvec() const; }; +/** + * @brief Compute derivatives of R(omega) w.r.t. rotation vector components. + * + * Given a rotation vector omega (axis * angle), computes the three 3x3 + * matrices dR/d(omega_k) for k = 0, 1, 2 using exact differentiation of + * the Rodrigues formula. + * + * @param omega Rotation vector (3 elements) + * @return 3x3x3 cube where slice k is dR/d(omega_k) + * + * @details Reference: Gallego & Yezzi, "A Compact Formula for the + * Derivative of a 3-D Rotation in Exponential Coordinates", + * J. Math. Imaging Vis., 2015. + */ +arma::cube dR_drotvec(const arma::vec::fixed<3>& omega); + +/** + * @brief Compute derivatives of R(omega) w.r.t. rotation vector components (dynamic size). + * @param omega Rotation vector (3 elements) + * @return 3x3x3 cube where slice k is dR/d(omega_k) + */ +arma::cube dR_drotvec(const arma::vec& omega); + // ========================================================================= // Convenience Free Functions // ========================================================================= From d5d1b02dff6fe07f8b0b6707945fb9f5bffe4f78 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 30 Mar 2026 12:38:00 +0200 Subject: [PATCH 3/7] Add dR_drotvec to Rotation class Introduce Rotation.dR_drotvec() to compute exact derivatives of the rotation matrix w.r.t. rotation vector components using the Rodrigues formula differentiation (Gallego & Yezzi, 2015). Returns (3,3,3) for single rotations or (N,3,3,3) for batches. Single-instance calls delegate to the C++ backend via _to_cpp(), while batched results are constructed by calling _CppRotation.from_quat for each quaternion. --- python-setup/simcoon/rotation.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/python-setup/simcoon/rotation.py b/python-setup/simcoon/rotation.py index cdc6cc7e..9f233410 100644 --- a/python-setup/simcoon/rotation.py +++ b/python-setup/simcoon/rotation.py @@ -379,6 +379,29 @@ def as_voigt_strain_rotation(self, active=True): """ return self._voigt_strain_matrices(active) + # ------------------------------------------------------------------ + # Rotation vector derivatives + # ------------------------------------------------------------------ + + def dR_drotvec(self): + """Derivatives of the rotation matrix w.r.t. rotation vector components. + + Uses the exact differentiation of the Rodrigues formula + (Gallego & Yezzi, J. Math. Imaging Vis., 2015). + + Returns + ------- + numpy.ndarray + Single: (3, 3, 3) where ``result[k]`` is dR/d(omega_k). + Batch: (N, 3, 3, 3) where ``result[n, k]`` is dR_n/d(omega_k). + """ + if not self._is_batch: + return self._to_cpp().dR_drotvec() + return np.array([ + _CppRotation.from_quat(q).dR_drotvec() + for q in self.as_quat() + ]) + # ------------------------------------------------------------------ # Compatibility helpers # ------------------------------------------------------------------ From c61dd39847473ff4b3d7ada9be398cf765e84f9b Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 30 Mar 2026 12:38:17 +0200 Subject: [PATCH 4/7] Add docs for rotation matrix derivatives Add documentation strings for dR_drotvec and dR_drotvec_free in doc_rotation.hpp. Documents computing the derivatives of the rotation matrix w.r.t. rotation vector components (returns a (3,3,3) array), cites the Rodrigues exact differentiation (Gallego & Yezzi 2015), and includes an example and notes for the free-function variant that accepts a rotation vector directly. --- .../docs/Libraries/Maths/doc_rotation.hpp | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp b/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp index a9cdc0df..5e004ccc 100644 --- a/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp +++ b/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp @@ -275,4 +275,48 @@ constexpr auto apply_stress_concentration = R"pbdoc( The rotated stress concentration tensor, same shape as input. )pbdoc"; +constexpr auto dR_drotvec = R"pbdoc( + Compute the derivatives of the rotation matrix w.r.t. rotation vector components. + + Uses the exact differentiation of the Rodrigues formula + (Gallego & Yezzi, J. Math. Imaging Vis., 2015). + + Returns + ------- + numpy.ndarray + A (3, 3, 3) array where ``result[k]`` is :math:`\partial R / \partial \omega_k`. + + Examples + -------- + .. code-block:: python + + import simcoon as smc + import numpy as np + + r = smc.Rotation.from_rotvec([0.1, 0.2, 0.3]) + dR = r.dR_drotvec() # (3, 3, 3) + # dR[0] = dR/d(omega_0), dR[1] = dR/d(omega_1), dR[2] = dR/d(omega_2) +)pbdoc"; + +constexpr auto dR_drotvec_free = R"pbdoc( + Compute the derivatives of R(omega) w.r.t. rotation vector components. + + Free function version — takes a rotation vector directly. + + Parameters + ---------- + rotvec : numpy.ndarray + A 1D array of 3 elements: rotation vector (axis * angle). + + Returns + ------- + numpy.ndarray + A (3, 3, 3) array where ``result[k]`` is :math:`\partial R / \partial \omega_k`. + + References + ---------- + Gallego & Yezzi, "A Compact Formula for the Derivative of a 3-D + Rotation in Exponential Coordinates", J. Math. Imaging Vis., 2015. +)pbdoc"; + } // namespace simcoon_docs From d0549ab26cc723f518ed9542f18a21457dffef30 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 30 Mar 2026 23:56:04 +0200 Subject: [PATCH 5/7] trigger CI From f3d980dffc1310a9752429f9458d0f16858123f4 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Fri, 17 Apr 2026 15:23:32 +0200 Subject: [PATCH 6/7] Refactor dR_drotvec with skew helpers Introduce internal skew() and skew_basis() helpers in an anonymous namespace and rewrite dR_drotvec to use them. Small-angle handling is simplified by assigning the three skew basis matrices in a loop; full-angle case now precomputes W and W2 and uses the helpers to build dW and dW2, reducing duplicated code and improving clarity/efficiency. Also remove the previous dynamic-vec overload and explicit manual skew assignments in favor of fixed-size operations. --- include/simcoon/Simulation/Maths/rotation.hpp | 7 -- python-setup/simcoon/rotation.py | 12 ++- .../Libraries/Maths/rotation.cpp | 36 ++++----- src/Simulation/Maths/rotation.cpp | 79 ++++++++----------- 4 files changed, 53 insertions(+), 81 deletions(-) diff --git a/include/simcoon/Simulation/Maths/rotation.hpp b/include/simcoon/Simulation/Maths/rotation.hpp index 22077c1f..2ffb78cb 100755 --- a/include/simcoon/Simulation/Maths/rotation.hpp +++ b/include/simcoon/Simulation/Maths/rotation.hpp @@ -478,13 +478,6 @@ class Rotation { */ arma::cube dR_drotvec(const arma::vec::fixed<3>& omega); -/** - * @brief Compute derivatives of R(omega) w.r.t. rotation vector components (dynamic size). - * @param omega Rotation vector (3 elements) - * @return 3x3x3 cube where slice k is dR/d(omega_k) - */ -arma::cube dR_drotvec(const arma::vec& omega); - // ========================================================================= // Convenience Free Functions // ========================================================================= diff --git a/python-setup/simcoon/rotation.py b/python-setup/simcoon/rotation.py index 9f233410..1131647c 100644 --- a/python-setup/simcoon/rotation.py +++ b/python-setup/simcoon/rotation.py @@ -10,7 +10,7 @@ import numpy as np from scipy.spatial.transform import Rotation as ScipyRotation -from simcoon._core import _CppRotation +from simcoon._core import _CppRotation, dR_drotvec as _dR_drotvec class Rotation(ScipyRotation): @@ -395,12 +395,10 @@ def dR_drotvec(self): Single: (3, 3, 3) where ``result[k]`` is dR/d(omega_k). Batch: (N, 3, 3, 3) where ``result[n, k]`` is dR_n/d(omega_k). """ - if not self._is_batch: - return self._to_cpp().dR_drotvec() - return np.array([ - _CppRotation.from_quat(q).dR_drotvec() - for q in self.as_quat() - ]) + rotvec = self.as_rotvec() + if rotvec.ndim == 1: + return _dR_drotvec(rotvec) + return np.array([_dR_drotvec(r) for r in rotvec]) # ------------------------------------------------------------------ # Compatibility helpers diff --git a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp index e7e2260f..7f6609a3 100644 --- a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp +++ b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp @@ -37,6 +37,19 @@ namespace { to_string(cols) + "), got (" + to_string(shape[0]) + ", " + to_string(shape[1]) + ")"); } } + + // Return a (3,3,3) numpy array where result[k,i,j] == cube(i,j,k). + // (Armadillo stores slices along the last axis; we expose slice index first + // so users can write result[k] for dR/d(omega_k).) + py::array_t cube333_to_kij_numpy(const arma::cube& c) { + py::array_t result({3, 3, 3}); + auto buf = result.mutable_unchecked<3>(); + for (int k = 0; k < 3; ++k) + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + buf(k, i, j) = c(i, j, k); + return result; + } } // anonymous namespace void register_rotation(py::module_& m) { @@ -141,33 +154,18 @@ void register_rotation(py::module_& m) { .def("dR_drotvec", [](const simcoon::Rotation& self) { - arma::cube c = self.dR_drotvec(); - // Return as (3, 3, 3) numpy array: result[k] = dR/d(omega_k) - py::array_t result({3, 3, 3}); - auto buf = result.mutable_unchecked<3>(); - for (int k = 0; k < 3; ++k) - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - buf(k, i, j) = c(i, j, k); - return result; + return cube333_to_kij_numpy(self.dR_drotvec()); }, simcoon_docs::dR_drotvec) ; - // Free function — takes a rotvec directly m.def("dR_drotvec", [](py::array_t rotvec) { validate_vector_size(rotvec, 3, "rotvec"); - vec omega = carma::arr_to_col(rotvec); - arma::cube c = simcoon::dR_drotvec(omega); - py::array_t result({3, 3, 3}); - auto buf = result.mutable_unchecked<3>(); - for (int k = 0; k < 3; ++k) - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - buf(k, i, j) = c(i, j, k); - return result; + auto r = rotvec.unchecked<1>(); + vec::fixed<3> omega = {r(0), r(1), r(2)}; + return cube333_to_kij_numpy(simcoon::dR_drotvec(omega)); }, py::arg("rotvec"), simcoon_docs::dR_drotvec_free); diff --git a/src/Simulation/Maths/rotation.cpp b/src/Simulation/Maths/rotation.cpp index d17c5f01..220a79b3 100755 --- a/src/Simulation/Maths/rotation.cpp +++ b/src/Simulation/Maths/rotation.cpp @@ -972,76 +972,59 @@ mat rotate_stress_concentration(const mat &B, const mat &DR, const bool &active) return vs*(B*trans(ve)); } -// ============================================================================= -// dR/d(rotvec) — Rodrigues differentiation (Gallego & Yezzi 2015) -// ============================================================================= - cube Rotation::dR_drotvec() const { return simcoon::dR_drotvec(as_rotvec()); } +namespace { + // Cross-product matrix [v]x such that [v]x * u = v x u. + inline mat::fixed<3,3> skew(const vec::fixed<3>& v) { + mat::fixed<3,3> S; + S.zeros(); + S(0,1) = -v(2); S(0,2) = v(1); + S(1,0) = v(2); S(1,2) = -v(0); + S(2,0) = -v(1); S(2,1) = v(0); + return S; + } + + inline mat::fixed<3,3> skew_basis(int k) { + vec::fixed<3> e; + e.zeros(); + e(k) = 1.0; + return skew(e); + } +} + cube dR_drotvec(const vec::fixed<3>& omega) { double theta = norm(omega); cube result(3, 3, 3); - // Skew-symmetric matrix of omega - mat::fixed<3,3> W; - W.zeros(); - W(0,1) = -omega(2); W(0,2) = omega(1); - W(1,0) = omega(2); W(1,2) = -omega(0); - W(2,0) = -omega(1); W(2,1) = omega(0); - if (theta < simcoon::iota) { - // Small angle limit: dR/d(omega_k) = [e_k]x - result.zeros(); - result(1,2,0) = -1.0; result(2,1,0) = 1.0; // [e0]x - result(0,2,1) = 1.0; result(2,0,1) = -1.0; // [e1]x - result(0,1,2) = -1.0; result(1,0,2) = 1.0; // [e2]x + // d/d(omega_k) [ exp([omega]x) ] at omega=0 is [e_k]x. + for (int k = 0; k < 3; ++k) { + result.slice(k) = skew_basis(k); + } return result; } + mat::fixed<3,3> W = skew(omega); + mat::fixed<3,3> W2 = W * W; + double s = sin(theta), c = cos(theta); double t2 = theta * theta; - - // Rodrigues coefficients and their derivatives w.r.t. theta - double a = s / theta; // sin(theta)/theta - double b = (1.0 - c) / t2; // (1-cos(theta))/theta^2 - double da = (c * theta - s) / t2; // d(sin(theta)/theta)/d(theta) - double db = (s * theta - 2.0 * (1.0 - c)) / (t2 * theta); // d((1-cos)/theta^2)/d(theta) - - mat::fixed<3,3> W2 = W * W; + double a = s / theta; + double b = (1.0 - c) / t2; + double da = (c * theta - s) / t2; + double db = (s * theta - 2.0 * (1.0 - c)) / (t2 * theta); for (int k = 0; k < 3; ++k) { - // dW/d(omega_k) = [e_k]x - mat::fixed<3,3> dW; - dW.zeros(); - // Using direct assignment for the skew of unit vector e_k - if (k == 0) { - dW(1,2) = -1.0; dW(2,1) = 1.0; - } else if (k == 1) { - dW(0,2) = 1.0; dW(2,0) = -1.0; - } else { - dW(0,1) = -1.0; dW(1,0) = 1.0; - } - - // dW^2/d(omega_k) = dW*W + W*dW + mat::fixed<3,3> dW = skew_basis(k); mat::fixed<3,3> dW2 = dW * W + W * dW; - - // d(theta)/d(omega_k) = omega_k / theta double dtk = omega(k) / theta; - result.slice(k) = da * dtk * W + a * dW + db * dtk * W2 + b * dW2; } return result; } -cube dR_drotvec(const vec& omega) { - if (omega.n_elem != 3) { - throw invalid_argument("dR_drotvec: rotation vector must have 3 elements"); - } - vec::fixed<3> omega_f = {omega(0), omega(1), omega(2)}; - return dR_drotvec(omega_f); -} - } //namespace simcoon From b504c4fa6ef9544ac8d92c13954e04376a78179b Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 20 Apr 2026 17:10:43 +0200 Subject: [PATCH 7/7] Expose dR_drotvec slices on last axis Change the dR_drotvec cube convention so the slice axis is last (result[:, :, k]) to match project-wide cube convention ((3,3,N), (6,6,N), ...). Updates: C++ Python wrappers (use carma::cube_to_arr), Python Rotation API (batch stacked along axis=-1), headers and docs wording, and added tests covering shapes, small-angle branch, finite-difference agreement, batch consistency, and input validation. --- include/simcoon/Simulation/Maths/rotation.hpp | 6 +- python-setup/simcoon/rotation.py | 9 +- .../docs/Libraries/Maths/doc_rotation.hpp | 9 +- .../Libraries/Maths/rotation.cpp | 17 +-- .../test/test_core/test_dR_drotvec.py | 141 ++++++++++++++++++ 5 files changed, 159 insertions(+), 23 deletions(-) create mode 100644 simcoon-python-builder/test/test_core/test_dR_drotvec.py diff --git a/include/simcoon/Simulation/Maths/rotation.hpp b/include/simcoon/Simulation/Maths/rotation.hpp index 2ffb78cb..1922a775 100755 --- a/include/simcoon/Simulation/Maths/rotation.hpp +++ b/include/simcoon/Simulation/Maths/rotation.hpp @@ -457,7 +457,8 @@ class Rotation { * Uses the exact differentiation of the Rodrigues formula * (Gallego & Yezzi, J. Math. Imaging Vis., 2015). * - * @return 3x3x3 cube where slice k is dR/d(omega_k) + * @return 3x3x3 cube where slice(k) is dR/d(omega_k) + * (Python exposure: result[:, :, k]). */ arma::cube dR_drotvec() const; }; @@ -470,7 +471,8 @@ class Rotation { * the Rodrigues formula. * * @param omega Rotation vector (3 elements) - * @return 3x3x3 cube where slice k is dR/d(omega_k) + * @return 3x3x3 cube where slice(k) is dR/d(omega_k) + * (Python exposure: result[:, :, k]). * * @details Reference: Gallego & Yezzi, "A Compact Formula for the * Derivative of a 3-D Rotation in Exponential Coordinates", diff --git a/python-setup/simcoon/rotation.py b/python-setup/simcoon/rotation.py index 1131647c..58d4a887 100644 --- a/python-setup/simcoon/rotation.py +++ b/python-setup/simcoon/rotation.py @@ -392,13 +392,16 @@ def dR_drotvec(self): Returns ------- numpy.ndarray - Single: (3, 3, 3) where ``result[k]`` is dR/d(omega_k). - Batch: (N, 3, 3, 3) where ``result[n, k]`` is dR_n/d(omega_k). + Single: (3, 3, 3) where ``result[:, :, k]`` is dR/d(omega_k). + Batch: (3, 3, 3, N) where ``result[:, :, k, n]`` is dR_n/d(omega_k). + + The slice axis is last, matching simcoon's project-wide cube + convention ((3,3,N), (6,6,N), ...). """ rotvec = self.as_rotvec() if rotvec.ndim == 1: return _dR_drotvec(rotvec) - return np.array([_dR_drotvec(r) for r in rotvec]) + return np.stack([_dR_drotvec(r) for r in rotvec], axis=-1) # ------------------------------------------------------------------ # Compatibility helpers diff --git a/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp b/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp index 5e004ccc..2886d905 100644 --- a/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp +++ b/simcoon-python-builder/include/simcoon/docs/Libraries/Maths/doc_rotation.hpp @@ -284,7 +284,9 @@ constexpr auto dR_drotvec = R"pbdoc( Returns ------- numpy.ndarray - A (3, 3, 3) array where ``result[k]`` is :math:`\partial R / \partial \omega_k`. + A (3, 3, 3) array where ``result[:, :, k]`` is + :math:`\partial R / \partial \omega_k`. The slice axis is last, + matching simcoon's project-wide cube convention. Examples -------- @@ -295,7 +297,7 @@ constexpr auto dR_drotvec = R"pbdoc( r = smc.Rotation.from_rotvec([0.1, 0.2, 0.3]) dR = r.dR_drotvec() # (3, 3, 3) - # dR[0] = dR/d(omega_0), dR[1] = dR/d(omega_1), dR[2] = dR/d(omega_2) + # dR[:, :, 0] = dR/d(omega_0), dR[:, :, 1] = dR/d(omega_1), dR[:, :, 2] = dR/d(omega_2) )pbdoc"; constexpr auto dR_drotvec_free = R"pbdoc( @@ -311,7 +313,8 @@ constexpr auto dR_drotvec_free = R"pbdoc( Returns ------- numpy.ndarray - A (3, 3, 3) array where ``result[k]`` is :math:`\partial R / \partial \omega_k`. + A (3, 3, 3) array where ``result[:, :, k]`` is + :math:`\partial R / \partial \omega_k`. References ---------- diff --git a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp index 7f6609a3..c08bf633 100644 --- a/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp +++ b/simcoon-python-builder/src/python_wrappers/Libraries/Maths/rotation.cpp @@ -37,19 +37,6 @@ namespace { to_string(cols) + "), got (" + to_string(shape[0]) + ", " + to_string(shape[1]) + ")"); } } - - // Return a (3,3,3) numpy array where result[k,i,j] == cube(i,j,k). - // (Armadillo stores slices along the last axis; we expose slice index first - // so users can write result[k] for dR/d(omega_k).) - py::array_t cube333_to_kij_numpy(const arma::cube& c) { - py::array_t result({3, 3, 3}); - auto buf = result.mutable_unchecked<3>(); - for (int k = 0; k < 3; ++k) - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - buf(k, i, j) = c(i, j, k); - return result; - } } // anonymous namespace void register_rotation(py::module_& m) { @@ -154,7 +141,7 @@ void register_rotation(py::module_& m) { .def("dR_drotvec", [](const simcoon::Rotation& self) { - return cube333_to_kij_numpy(self.dR_drotvec()); + return carma::cube_to_arr(self.dR_drotvec()); }, simcoon_docs::dR_drotvec) @@ -165,7 +152,7 @@ void register_rotation(py::module_& m) { validate_vector_size(rotvec, 3, "rotvec"); auto r = rotvec.unchecked<1>(); vec::fixed<3> omega = {r(0), r(1), r(2)}; - return cube333_to_kij_numpy(simcoon::dR_drotvec(omega)); + return carma::cube_to_arr(simcoon::dR_drotvec(omega)); }, py::arg("rotvec"), simcoon_docs::dR_drotvec_free); diff --git a/simcoon-python-builder/test/test_core/test_dR_drotvec.py b/simcoon-python-builder/test/test_core/test_dR_drotvec.py new file mode 100644 index 00000000..aeddeeb1 --- /dev/null +++ b/simcoon-python-builder/test/test_core/test_dR_drotvec.py @@ -0,0 +1,141 @@ +"""Tests for ``dR_drotvec``: derivatives of R(omega) w.r.t. rotation-vector components. + +Layout convention (matches simcoon's project-wide cube convention): + single : shape (3, 3, 3), ``result[:, :, k]`` is dR/d(omega_k) + batch : shape (3, 3, 3, N), ``result[:, :, k, n]`` is dR_n/d(omega_k) +""" + +import numpy as np +import pytest + +import simcoon as sim + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture +def rotvec_general(): + """Generic non-axis-aligned rotvec (moderate angle ~0.37 rad).""" + return np.array([0.1, -0.2, 0.3]) + + +@pytest.fixture +def rotvecs_batch(): + """A batch of diverse rotvecs, including identity.""" + return np.array([ + [0.1, -0.2, 0.3], + [0.0, 0.0, 0.5], + [0.7, 0.0, 0.0], + [0.0, 0.0, 0.0], # identity + [np.pi / 2 - 0.05, 0.0, 0.0], # near pi/2 + ]) + + +# Reference skew matrices [e_k]x for k=0,1,2 +E_SKEW = np.array([ + [[0.0, 0.0, 0.0], [0.0, 0.0, -1.0], [0.0, 1.0, 0.0]], # [e_0]x + [[0.0, 0.0, 1.0], [0.0, 0.0, 0.0], [-1.0, 0.0, 0.0]], # [e_1]x + [[0.0, -1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]], # [e_2]x +]) + + +# --------------------------------------------------------------------------- +# Shape + API symmetry +# --------------------------------------------------------------------------- + +class TestShape: + def test_single_shape(self, rotvec_general): + assert sim.dR_drotvec(rotvec_general).shape == (3, 3, 3) + + def test_member_shape(self, rotvec_general): + r = sim.Rotation.from_rotvec(rotvec_general) + assert r.dR_drotvec().shape == (3, 3, 3) + + def test_batch_shape(self, rotvecs_batch): + r = sim.Rotation.from_rotvec(rotvecs_batch) + assert r.dR_drotvec().shape == (3, 3, 3, len(rotvecs_batch)) + + def test_free_eq_member(self, rotvec_general): + r = sim.Rotation.from_rotvec(rotvec_general) + np.testing.assert_array_equal(sim.dR_drotvec(rotvec_general), r.dR_drotvec()) + + +# --------------------------------------------------------------------------- +# Small-angle / identity branch +# --------------------------------------------------------------------------- + +class TestSmallAngle: + def test_zero_rotvec_equals_basis_skews(self): + """At omega=0: dR/d(omega_k) == [e_k]x for each k.""" + result = sim.dR_drotvec(np.zeros(3)) + for k in range(3): + np.testing.assert_allclose(result[:, :, k], E_SKEW[k], atol=0.0) + + def test_identity_rotation_matches_zero_rotvec(self): + result = sim.Rotation.identity().dR_drotvec() + expected = sim.dR_drotvec(np.zeros(3)) + np.testing.assert_allclose(result, expected, atol=0.0) + + def test_below_iota_threshold_uses_small_angle_branch(self): + """Very small rotvec should give near-basis-skew slices.""" + tiny = np.array([1e-14, 0.0, 0.0]) + result = sim.dR_drotvec(tiny) + for k in range(3): + np.testing.assert_allclose(result[:, :, k], E_SKEW[k], atol=1e-12) + + +# --------------------------------------------------------------------------- +# Finite-difference agreement +# --------------------------------------------------------------------------- + +def _R(omega): + return sim.Rotation.from_rotvec(omega).as_matrix() + + +def _fd_slice(omega, k, eps=1e-6): + pert = omega.copy() + pert[k] += eps + return (_R(pert) - _R(omega - np.eye(3)[k] * eps)) / (2.0 * eps) + + +class TestFiniteDifference: + @pytest.mark.parametrize("rotvec", [ + np.array([0.1, -0.2, 0.3]), + np.array([0.0, 0.0, 0.9]), + np.array([1.2, -0.4, 0.3]), # |omega| ~ 1.3 rad + np.array([np.pi / 2 - 0.01, 0.0, 0.0]), + ]) + def test_fd_matches_analytic(self, rotvec): + analytic = sim.dR_drotvec(rotvec) + for k in range(3): + fd = _fd_slice(rotvec, k) + np.testing.assert_allclose(analytic[:, :, k], fd, atol=1e-7, rtol=1e-7) + + +# --------------------------------------------------------------------------- +# Batch consistency +# --------------------------------------------------------------------------- + +class TestBatch: + def test_batch_slice_matches_single_call(self, rotvecs_batch): + r = sim.Rotation.from_rotvec(rotvecs_batch) + batch = r.dR_drotvec() + for n, rv in enumerate(rotvecs_batch): + single = sim.dR_drotvec(rv) + np.testing.assert_allclose(batch[:, :, :, n], single, atol=0.0) + + +# --------------------------------------------------------------------------- +# Input validation +# --------------------------------------------------------------------------- + +class TestValidation: + def test_rejects_wrong_size(self): + with pytest.raises(Exception): + sim.dR_drotvec(np.zeros(2)) + + def test_rejects_2d_input(self): + with pytest.raises(Exception): + sim.dR_drotvec(np.zeros((3, 1)))