From 25e40276fb8dbeaf0be2b916c497d0bdd289380c Mon Sep 17 00:00:00 2001 From: Brianna Major Date: Mon, 16 Mar 2026 10:47:38 -0400 Subject: [PATCH 1/7] =?UTF-8?q?Add=20de=20la=20Vall=C3=A9e=20Poussin=20ker?= =?UTF-8?q?nel=20on=20SO(3)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implements the radially symmetric kernel K(ω) = C · cos(ω/2)^(2κ) with normalization constant C = B(3/2, 1/2) / B(3/2, κ + 1/2) and shape parameter κ derived analytically from the half-width. - Abstract SO3Kernel base class for rotation group kernels - Concrete DeLaValleePoussinKernel with eval, eval_centered, and misorientation_angle methods Signed-off-by: Brianna Major --- hexrd/texture/__init__.py | 0 hexrd/texture/kernels.py | 233 ++++++++++++++++++++++ tests/test_delavalleepoussin_kernel.py | 255 +++++++++++++++++++++++++ 3 files changed, 488 insertions(+) create mode 100644 hexrd/texture/__init__.py create mode 100644 hexrd/texture/kernels.py create mode 100644 tests/test_delavalleepoussin_kernel.py diff --git a/hexrd/texture/__init__.py b/hexrd/texture/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/hexrd/texture/kernels.py b/hexrd/texture/kernels.py new file mode 100644 index 000000000..dcf2a5289 --- /dev/null +++ b/hexrd/texture/kernels.py @@ -0,0 +1,233 @@ +""" +SO(3) Kernel Functions for Texture Analysis + +Implements radial basis functions on the rotation group SO(3) used for +constructing smooth orientation distribution functions. +""" + +import warnings + +import numpy as np +from abc import ABC, abstractmethod +from scipy.special import beta as betafn + + +class SO3Kernel(ABC): + """ + Abstract base class for kernels on the SO(3) rotation group. + + All SO(3) kernels should inherit from this class and implement + the eval() method for kernel evaluation. + """ + + @abstractmethod + def eval( + self, R1: np.ndarray, R2: np.ndarray + ) -> np.ndarray: + """ + Evaluate kernel between two rotations. + + Parameters + ---------- + R1, R2 : array_like + Rotation matrices of shape (..., 3, 3) + + Returns + ------- + numpy.ndarray + Kernel values + """ + pass + + +class DeLaValleePoussinKernel(SO3Kernel): + """ + De la Vallée Poussin kernel on SO(3). + + A radially symmetric kernel on SO(3) defined by: + + K(ω) = C · cos(ω/2)^(2κ) + + where ω is the misorientation angle, κ is the shape parameter, and + C = B(3/2, 1/2) / B(3/2, κ + 1/2) is the normalization constant + (B denotes the Beta function). + + The halfwidth h is the angle at which the kernel drops to half its + peak value. It is related to κ analytically by: + + κ = ln(0.5) / (2 · ln(cos(h/2))) + + Parameters + ---------- + halfwidth : float + Half-width parameter in radians — the angle at which K drops + to half its maximum. Must be > 0, typically in [π/180, π/2]. + + Attributes + ---------- + halfwidth : float + Half-width parameter in radians + kappa : float + Shape parameter κ derived from half-width + norm_constant : float + Normalization constant from the Beta function + + Examples + -------- + >>> kernel = DeLaValleePoussinKernel(halfwidth=np.radians(15)) + >>> R1 = np.eye(3) + >>> R2 = np.eye(3) # Same orientation + >>> value = kernel.eval(R1, R2) # Maximum value = C + """ + + def __init__(self, halfwidth: float) -> None: + """ + Initialize de la Vallée Poussin kernel. + + Parameters + ---------- + halfwidth : float + Half-width parameter in radians + """ + if halfwidth <= 0: + raise ValueError("Half-width must be positive") + if halfwidth > np.pi / 2: + warnings.warn( + f"Large half-width " + f"{np.degrees(halfwidth):.1f}° may produce " + f"very broad distributions", + UserWarning, + ) + + self._halfwidth = float(halfwidth) + + # Shape parameter κ from the half-maximum condition K(h) = K(0)/2: + # cos(h/2)^(2κ) = 1/2 + # κ = ln(0.5) / (2·ln(cos(h/2))) + self._kappa = 0.5 * np.log(0.5) / np.log(np.cos(halfwidth / 2.0)) + + # Normalization constant: C = B(3/2, 1/2) / B(3/2, κ + 1/2) + self._C = betafn(1.5, 0.5) / betafn(1.5, self._kappa + 0.5) + + @property + def halfwidth(self) -> float: + """float: Half-width in radians (angle where K = K_max / 2).""" + return self._halfwidth + + @property + def kappa(self) -> float: + """float: Shape parameter κ.""" + return self._kappa + + @property + def norm_constant(self) -> float: + """float: Normalization constant from Beta function.""" + return self._C + + def misorientation_angle( + self, R1: np.ndarray, R2: np.ndarray + ) -> np.ndarray: + """ + Calculate misorientation angle between rotation matrices. + + Uses the formula: cos(ω) = (trace(R1^T @ R2) - 1) / 2 + + Parameters + ---------- + R1, R2 : array_like + Rotation matrices of shape (..., 3, 3) + + Returns + ------- + numpy.ndarray + Misorientation angles in radians, shape matches input broadcasting + """ + R1 = np.asarray(R1) + R2 = np.asarray(R2) + + if R1.shape[-2:] != (3, 3) or R2.shape[-2:] != (3, 3): + raise ValueError( + "Input matrices must have shape (..., 3, 3)" + ) + + # Relative rotation: R1^T @ R2 + R1_T = np.swapaxes(R1, -2, -1) + relative = np.matmul(R1_T, R2) + trace = np.trace(relative, axis1=-2, axis2=-1) + + # cos(ω) = (tr(R) - 1) / 2, clamped for numerical safety + cos_omega = np.clip((trace - 1.0) / 2.0, -1.0, 1.0) + return np.arccos(cos_omega) + + def eval( + self, R1: np.ndarray, R2: np.ndarray + ) -> np.ndarray: + """ + Evaluate de la Vallée Poussin kernel between rotations. + + K(ω) = C · cos(ω/2)^(2κ) + + Parameters + ---------- + R1, R2 : array_like + Rotation matrices of shape (..., 3, 3) + + Returns + ------- + numpy.ndarray + Kernel values, shape matches broadcasting of R1 and R2 + + Examples + -------- + >>> kernel = DeLaValleePoussinKernel( + ... halfwidth=np.radians(10) + ... ) + >>> value = kernel.eval(np.eye(3), np.eye(3)) + """ + omega = self.misorientation_angle(R1, R2) + co2 = np.cos(omega / 2.0) + return self._C * co2 ** (2.0 * self._kappa) + + def eval_centered( + self, R: np.ndarray, center: np.ndarray + ) -> np.ndarray: + """ + Evaluate kernel centered at a specific orientation. + + Equivalent to eval(center, R) with clearer semantics + for unimodal distributions. + + Parameters + ---------- + R : array_like + Rotation matrices to evaluate at, shape (..., 3, 3) + center : array_like + Central (modal) orientation, shape (3, 3) + + Returns + ------- + numpy.ndarray + Kernel values relative to center orientation + """ + return self.eval(center, R) + + def __repr__(self) -> str: + """String representation of kernel.""" + hw_deg = np.degrees(self.halfwidth) + return ( + f"DeLaValleePoussinKernel(" + f"halfwidth={self.halfwidth:.6f} rad " + f"= {hw_deg:.2f}°, " + f"kappa={self.kappa:.3f})" + ) + + def __str__(self) -> str: + """Human-readable description.""" + hw_deg = np.degrees(self.halfwidth) + return ( + f"de la Vallée Poussin kernel with " + f"{hw_deg:.1f}° half-width\n" + f"kappa = {self.kappa:.3f}, " + f"C = {self.norm_constant:.6e}\n" + f"K(ω) = C · cos(ω/2)^(2κ)" + ) diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py new file mode 100644 index 000000000..90246af4c --- /dev/null +++ b/tests/test_delavalleepoussin_kernel.py @@ -0,0 +1,255 @@ +"""Tests for DeLaValleePoussinKernel.""" + +import numpy as np +import unittest +from scipy.special import beta as betafn + +from hexrd.texture.kernels import DeLaValleePoussinKernel + + +MTEX_KERNEL_KAPPA = 90.0 +MTEX_KERNEL_OMEGA = np.array([ + 0.0000000000, + 0.5235987756, + 0.7853981634, + 1.0471975512, + 1.5707963268, + 3.1415926536, +]) +MTEX_KERNEL_VALUES = np.array([ + 1532.2892935006, + 2.9869009813, + 0.0009910670, + 0.0000000087, + 0.0000000000, + 0.0000000000, +]) +MTEX_HALFWIDTH_RAD = 0.1745329252 +MTEX_HALFWIDTH_KAPPA = 90.9031059932 + + +def _rotation_about_z(angle: float) -> np.ndarray: + cos_angle = np.cos(angle) + sin_angle = np.sin(angle) + return np.array([ + [cos_angle, -sin_angle, 0.0], + [sin_angle, cos_angle, 0.0], + [0.0, 0.0, 1.0], + ]) + + +class TestDeLaValleePoussinKernel(unittest.TestCase): + """Test DeLaValleePoussinKernel.""" + + def test_kernel_matches_mtex_ground_truth_values(self): + """Test kernel evaluation against values generated from MTEX.""" + halfwidth = 2.0 * np.arccos( + 0.5 ** (1.0 / (2.0 * MTEX_KERNEL_KAPPA)) + ) + kernel = DeLaValleePoussinKernel(halfwidth=halfwidth) + + identity = np.eye(3) + rotations = np.array([ + _rotation_about_z(angle) + for angle in MTEX_KERNEL_OMEGA + ]) + + values = kernel.eval(identity, rotations) + + np.testing.assert_allclose( + kernel.kappa, MTEX_KERNEL_KAPPA, rtol=0.0, atol=1e-10 + ) + np.testing.assert_allclose( + values, MTEX_KERNEL_VALUES, rtol=1e-8, atol=1e-10 + ) + + def test_halfwidth_conversion_matches_mtex_ground_truth(self): + """Test half-width to kappa conversion against MTEX output.""" + kernel = DeLaValleePoussinKernel(halfwidth=MTEX_HALFWIDTH_RAD) + + np.testing.assert_allclose( + kernel.kappa, MTEX_HALFWIDTH_KAPPA, rtol=1e-10, atol=1e-10 + ) + np.testing.assert_allclose( + kernel.halfwidth, MTEX_HALFWIDTH_RAD, rtol=0.0, atol=1e-10 + ) + + def test_kappa_from_halfwidth(self): + """Test that κ is derived analytically from half-width.""" + halfwidth_rad = np.radians(4.0) + kernel = DeLaValleePoussinKernel(halfwidth=halfwidth_rad) + + # Regression guard: verify against the formula directly + expected_kappa = ( + 0.5 * np.log(0.5) + / np.log(np.cos(halfwidth_rad / 2.0)) + ) + self.assertAlmostEqual( + kernel.kappa, expected_kappa, places=10 + ) + + def test_kernel_evaluation_at_identity(self): + """Test that K(0) equals the normalization constant. + + At ω = 0, cos(0/2) = 1, so K(0) = C · 1^(2κ) = C. + """ + halfwidth_rad = np.radians(4.0) + kernel = DeLaValleePoussinKernel(halfwidth=halfwidth_rad) + + identity = np.eye(3) + value = kernel.eval(identity, identity) + + # Verify against independently computed norm constant + expected_C = ( + betafn(1.5, 0.5) + / betafn(1.5, kernel.kappa + 0.5) + ) + self.assertAlmostEqual(value, expected_C, places=5) + self.assertAlmostEqual( + value, kernel.norm_constant, places=10 + ) + + def test_kernel_decreases_with_misorientation(self): + """Test that kernel value decreases as misorientation grows.""" + halfwidth_rad = np.radians(4.0) + kernel = DeLaValleePoussinKernel(halfwidth=halfwidth_rad) + + identity = np.eye(3) + + # 22.5° rotation about x-axis + rotated = np.array([ + [1.0000, 0.0000, 0.0000], + [0.0000, 0.9239, -0.3827], + [0.0000, 0.3827, 0.9239], + ]) + + value_identity = kernel.eval(identity, identity) + value_rotated = kernel.eval(identity, rotated) + + # 22.5° is well beyond the 4° halfwidth, so kernel ≈ 0 + self.assertGreater(value_identity, value_rotated) + self.assertAlmostEqual(value_rotated, 0.0, places=3) + + def test_halfwidth_is_half_maximum(self): + """Test that K(halfwidth) = K(0) / 2.""" + halfwidth_rad = np.radians(4.0) + kernel = DeLaValleePoussinKernel(halfwidth=halfwidth_rad) + + # K(0) = norm_constant + peak = kernel.norm_constant + + # Evaluate at exactly the halfwidth angle + co2 = np.cos(halfwidth_rad / 2.0) + value_at_hw = ( + kernel.norm_constant * co2 ** (2.0 * kernel.kappa) + ) + + self.assertAlmostEqual(value_at_hw, peak / 2.0, places=5) + + # --- Invalid input tests --- + + def test_negative_halfwidth_raises(self): + """Test that negative half-width raises ValueError.""" + with self.assertRaises(ValueError): + DeLaValleePoussinKernel(halfwidth=-0.1) + + def test_zero_halfwidth_raises(self): + """Test that zero half-width raises ValueError.""" + with self.assertRaises(ValueError): + DeLaValleePoussinKernel(halfwidth=0.0) + + def test_large_halfwidth_warns(self): + """Test that half-width > π/2 emits a warning.""" + with self.assertWarns(UserWarning): + DeLaValleePoussinKernel(halfwidth=np.radians(100)) + + # --- misorientation_angle tests --- + + def test_misorientation_angle_identity(self): + """Test misorientation angle between identical rotations is 0.""" + kernel = DeLaValleePoussinKernel(halfwidth=np.radians(10)) + angle = kernel.misorientation_angle(np.eye(3), np.eye(3)) + self.assertAlmostEqual(float(angle), 0.0, places=10) + + def test_misorientation_angle_known_rotation(self): + """Test misorientation angle for a known rotation.""" + kernel = DeLaValleePoussinKernel(halfwidth=np.radians(10)) + + # 90° rotation about z-axis + R90z = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1], + ], dtype=float) + + angle = kernel.misorientation_angle(np.eye(3), R90z) + self.assertAlmostEqual( + float(angle), np.radians(90), places=8 + ) + + def test_misorientation_angle_bad_shape(self): + """Test that non-(3,3) inputs raise ValueError.""" + kernel = DeLaValleePoussinKernel(halfwidth=np.radians(10)) + with self.assertRaises(ValueError): + kernel.misorientation_angle( + np.eye(2), np.eye(3) + ) + + # --- eval_centered test --- + + def test_eval_centered_matches_eval(self): + """Test that eval_centered(R, center) == eval(center, R).""" + kernel = DeLaValleePoussinKernel( + halfwidth=np.radians(10) + ) + center = np.eye(3) + R = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1], + ], dtype=float) + + val_eval = kernel.eval(center, R) + val_centered = kernel.eval_centered(R, center) + self.assertAlmostEqual( + float(val_eval), float(val_centered), places=10 + ) + + # --- Batch evaluation test --- + + def test_batch_evaluation(self): + """Test eval with a batch of rotation matrices (N, 3, 3).""" + kernel = DeLaValleePoussinKernel( + halfwidth=np.radians(10) + ) + identity = np.eye(3) + + # Build batch: identity + 90° rotations about x, y, z + R90x = np.array([ + [1, 0, 0], + [0, 0, -1], + [0, 1, 0], + ], dtype=float) + R90y = np.array([ + [ 0, 0, 1], + [ 0, 1, 0], + [-1, 0, 0], + ], dtype=float) + R90z = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1], + ], dtype=float) + batch = np.stack([identity, R90x, R90y, R90z]) + + values = kernel.eval(identity, batch) + + self.assertEqual(values.shape, (4,)) + # First entry (identity vs identity) should be the peak + self.assertAlmostEqual( + values[0], kernel.norm_constant, places=5 + ) + # Remaining entries (90° away) should be near zero + for val in values[1:]: + self.assertAlmostEqual(float(val), 0.0, places=5) + From 92e2d319191aa5bd1f9803526571c12f3838fac5 Mon Sep 17 00:00:00 2001 From: Brianna Major Date: Wed, 3 Jun 2026 10:53:35 -0400 Subject: [PATCH 2/7] Move texture kernels under phase_transition Signed-off-by: Brianna Major --- hexrd/{texture => phase_transition}/__init__.py | 0 hexrd/phase_transition/texture/__init__.py | 9 +++++++++ hexrd/{ => phase_transition}/texture/kernels.py | 0 tests/test_delavalleepoussin_kernel.py | 2 +- 4 files changed, 10 insertions(+), 1 deletion(-) rename hexrd/{texture => phase_transition}/__init__.py (100%) create mode 100644 hexrd/phase_transition/texture/__init__.py rename hexrd/{ => phase_transition}/texture/kernels.py (100%) diff --git a/hexrd/texture/__init__.py b/hexrd/phase_transition/__init__.py similarity index 100% rename from hexrd/texture/__init__.py rename to hexrd/phase_transition/__init__.py diff --git a/hexrd/phase_transition/texture/__init__.py b/hexrd/phase_transition/texture/__init__.py new file mode 100644 index 000000000..fab91ee72 --- /dev/null +++ b/hexrd/phase_transition/texture/__init__.py @@ -0,0 +1,9 @@ +from hexrd.phase_transition.texture.kernels import ( + DeLaValleePoussinKernel, + SO3Kernel, +) + +__all__ = [ + 'DeLaValleePoussinKernel', + 'SO3Kernel', +] diff --git a/hexrd/texture/kernels.py b/hexrd/phase_transition/texture/kernels.py similarity index 100% rename from hexrd/texture/kernels.py rename to hexrd/phase_transition/texture/kernels.py diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py index 90246af4c..99e0ab4fa 100644 --- a/tests/test_delavalleepoussin_kernel.py +++ b/tests/test_delavalleepoussin_kernel.py @@ -4,7 +4,7 @@ import unittest from scipy.special import beta as betafn -from hexrd.texture.kernels import DeLaValleePoussinKernel +from hexrd.phase_transition.texture.kernels import DeLaValleePoussinKernel MTEX_KERNEL_KAPPA = 90.0 From 5487980aeb4d3249ba7905a3cceb3a1104ddf4c8 Mon Sep 17 00:00:00 2001 From: Brianna Major Date: Wed, 3 Jun 2026 10:54:41 -0400 Subject: [PATCH 3/7] Add kernel string representation coverage Signed-off-by: Brianna Major --- tests/test_delavalleepoussin_kernel.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py index 99e0ab4fa..1fe4c8e95 100644 --- a/tests/test_delavalleepoussin_kernel.py +++ b/tests/test_delavalleepoussin_kernel.py @@ -108,6 +108,20 @@ def test_kernel_evaluation_at_identity(self): self.assertAlmostEqual( value, kernel.norm_constant, places=10 ) + + def test_string_representations(self): + """Test string representation methods.""" + kernel = DeLaValleePoussinKernel(halfwidth=np.radians(4.0)) + + repr_str = repr(kernel) + self.assertIn('DeLaValleePoussinKernel', repr_str) + self.assertIn('halfwidth=', repr_str) + self.assertIn('kappa=', repr_str) + + str_repr = str(kernel) + self.assertIn('de la Vallée Poussin kernel', str_repr) + self.assertIn('half-width', str_repr) + self.assertIn('K(ω)', str_repr) def test_kernel_decreases_with_misorientation(self): """Test that kernel value decreases as misorientation grows.""" From a41638f9d327e09fd76c4d3032929892c68edb15 Mon Sep 17 00:00:00 2001 From: Brianna Major Date: Wed, 3 Jun 2026 11:03:31 -0400 Subject: [PATCH 4/7] Use core rotations for symmetry-aware misorientation Add optional crystal and sample symmetry support to DeLaValleePoussinKernel. When symmetry is provided, delegate misorientation-angle calculation to hexrd.core.rotations.misorientation. Signed-off-by: Brianna Major --- hexrd/phase_transition/texture/kernels.py | 108 +++++++++++++++++++++- tests/test_delavalleepoussin_kernel.py | 52 +++++++++++ 2 files changed, 158 insertions(+), 2 deletions(-) diff --git a/hexrd/phase_transition/texture/kernels.py b/hexrd/phase_transition/texture/kernels.py index dcf2a5289..738e6662b 100644 --- a/hexrd/phase_transition/texture/kernels.py +++ b/hexrd/phase_transition/texture/kernels.py @@ -5,12 +5,40 @@ constructing smooth orientation distribution functions. """ +from abc import ABC, abstractmethod +from typing import Optional, Union import warnings import numpy as np -from abc import ABC, abstractmethod from scipy.special import beta as betafn +from hexrd.core import rotations + + +_IDENTITY_SYMMETRY = np.c_[1.0, 0.0, 0.0, 0.0].T +_SAMPLE_SYMMETRY_TO_LAUE = { + 'triclinic': 'ci', + 'monoclinic': 'c2h', + 'orthorhombic': 'd2h', +} +_Symmetry = Optional[Union[str, np.ndarray]] + + +def _symmetry_quaternions( + symmetry: _Symmetry, + *, + symtype: str, +) -> Optional[np.ndarray]: + if symmetry is None: + return None + if isinstance(symmetry, np.ndarray): + return symmetry + if not isinstance(symmetry, str): + raise TypeError(f"{symtype} symmetry must be a string or ndarray") + if symtype == 'sample': + symmetry = _SAMPLE_SYMMETRY_TO_LAUE.get(symmetry.lower(), symmetry) + return rotations.quatOfLaueGroup(symmetry) + class SO3Kernel(ABC): """ @@ -62,6 +90,14 @@ class DeLaValleePoussinKernel(SO3Kernel): halfwidth : float Half-width parameter in radians — the angle at which K drops to half its maximum. Must be > 0, typically in [π/180, π/2]. + crystal_symmetry : str or numpy.ndarray, optional + Crystal symmetry used to compute symmetry-reduced misorientation + angles. Strings are interpreted as Laue group labels. + sample_symmetry : str or numpy.ndarray, optional + Sample symmetry used to compute symmetry-reduced misorientation + angles. Strings are interpreted as Laue group labels, with + triclinic, monoclinic, and orthorhombic sample labels mapped to + their corresponding finite Laue groups. Attributes ---------- @@ -80,7 +116,12 @@ class DeLaValleePoussinKernel(SO3Kernel): >>> value = kernel.eval(R1, R2) # Maximum value = C """ - def __init__(self, halfwidth: float) -> None: + def __init__( + self, + halfwidth: float, + crystal_symmetry: _Symmetry = None, + sample_symmetry: _Symmetry = None, + ) -> None: """ Initialize de la Vallée Poussin kernel. @@ -88,6 +129,10 @@ def __init__(self, halfwidth: float) -> None: ---------- halfwidth : float Half-width parameter in radians + crystal_symmetry : str or numpy.ndarray, optional + Crystal symmetry for symmetry-reduced misorientation angles + sample_symmetry : str or numpy.ndarray, optional + Sample symmetry for symmetry-reduced misorientation angles """ if halfwidth <= 0: raise ValueError("Half-width must be positive") @@ -108,6 +153,22 @@ def __init__(self, halfwidth: float) -> None: # Normalization constant: C = B(3/2, 1/2) / B(3/2, κ + 1/2) self._C = betafn(1.5, 0.5) / betafn(1.5, self._kappa + 0.5) + self._crystal_symmetry = _symmetry_quaternions( + crystal_symmetry, + symtype='crystal', + ) + self._sample_symmetry = _symmetry_quaternions( + sample_symmetry, + symtype='sample', + ) + + @property + def has_symmetry(self) -> bool: + """bool: Whether symmetry reduction is enabled.""" + return ( + self._crystal_symmetry is not None + or self._sample_symmetry is not None + ) @property def halfwidth(self) -> float: @@ -150,6 +211,9 @@ def misorientation_angle( "Input matrices must have shape (..., 3, 3)" ) + if self.has_symmetry: + return self._symmetry_reduced_misorientation_angle(R1, R2) + # Relative rotation: R1^T @ R2 R1_T = np.swapaxes(R1, -2, -1) relative = np.matmul(R1_T, R2) @@ -159,6 +223,46 @@ def misorientation_angle( cos_omega = np.clip((trace - 1.0) / 2.0, -1.0, 1.0) return np.arccos(cos_omega) + def _symmetry_reduced_misorientation_angle( + self, + R1: np.ndarray, + R2: np.ndarray, + ) -> np.ndarray: + output_shape = np.broadcast_shapes(R1.shape[:-2], R2.shape[:-2]) + R1_flat = np.broadcast_to( + R1, + output_shape + (3, 3), + ).reshape(-1, 3, 3) + R2_flat = np.broadcast_to( + R2, + output_shape + (3, 3), + ).reshape(-1, 3, 3) + + angles = np.empty(R1_flat.shape[0], dtype=float) + symmetries = self._misorientation_symmetries() + for i, (r1, r2) in enumerate(zip(R1_flat, R2_flat)): + q1 = rotations.quatOfRotMat(r1).reshape(4, 1) + q2 = rotations.quatOfRotMat(r2).reshape(4, 1) + angle, _ = rotations.misorientation( + q1, + q2, + symmetries=symmetries, + ) + angles[i] = angle[0] + + if output_shape == (): + return angles[0] + return angles.reshape(output_shape) + + def _misorientation_symmetries(self) -> tuple[np.ndarray, ...]: + crystal_symmetry = self._crystal_symmetry + if crystal_symmetry is None and self._sample_symmetry is not None: + crystal_symmetry = _IDENTITY_SYMMETRY + symmetries = (crystal_symmetry,) if crystal_symmetry is not None else () + if self._sample_symmetry is not None: + symmetries = symmetries + (self._sample_symmetry,) + return symmetries + def eval( self, R1: np.ndarray, R2: np.ndarray ) -> np.ndarray: diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py index 1fe4c8e95..fa393df59 100644 --- a/tests/test_delavalleepoussin_kernel.py +++ b/tests/test_delavalleepoussin_kernel.py @@ -201,6 +201,58 @@ def test_misorientation_angle_known_rotation(self): float(angle), np.radians(90), places=8 ) + def test_misorientation_angle_uses_crystal_symmetry(self): + """Test symmetry-reduced misorientation angle.""" + kernel = DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry='d4h', + ) + + R90z = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1], + ], dtype=float) + + angle = kernel.misorientation_angle(np.eye(3), R90z) + self.assertAlmostEqual(float(angle), 0.0, places=8) + + def test_symmetric_kernel_value_is_peak_for_equivalent_rotations(self): + """Test symmetry-equivalent rotations evaluate to the peak value.""" + kernel = DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry='d4h', + ) + + R90z = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1], + ], dtype=float) + + value = kernel.eval(np.eye(3), R90z) + self.assertAlmostEqual(value, kernel.norm_constant, places=8) + + def test_batch_misorientation_angle_uses_crystal_symmetry(self): + """Test symmetry-reduced misorientation for batches.""" + kernel = DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry='d4h', + ) + + R45z = _rotation_about_z(np.radians(45.0)) + R90z = _rotation_about_z(np.radians(90.0)) + batch = np.stack([np.eye(3), R45z, R90z]) + + angles = kernel.misorientation_angle(np.eye(3), batch) + + self.assertEqual(angles.shape, (3,)) + np.testing.assert_allclose( + angles, + [0.0, np.radians(45.0), 0.0], + atol=1e-8, + ) + def test_misorientation_angle_bad_shape(self): """Test that non-(3,3) inputs raise ValueError.""" kernel = DeLaValleePoussinKernel(halfwidth=np.radians(10)) From 91e615d88919adc58f538704b98ff63a0fc48063 Mon Sep 17 00:00:00 2001 From: Brianna Major Date: Wed, 3 Jun 2026 11:28:00 -0400 Subject: [PATCH 5/7] Polish symmetry-aware kernel API - Update misorientation_angle docstring - Use VslueError for symmetry validation failures - Update return annotations for scalar-or-batched kernel operations - Add tests covering the symmetry error behavior. Signed-off-by: Brianna Major --- hexrd/phase_transition/texture/kernels.py | 35 +++++++++++++++-------- tests/test_delavalleepoussin_kernel.py | 24 ++++++++++++++++ 2 files changed, 47 insertions(+), 12 deletions(-) diff --git a/hexrd/phase_transition/texture/kernels.py b/hexrd/phase_transition/texture/kernels.py index 738e6662b..f489cbaed 100644 --- a/hexrd/phase_transition/texture/kernels.py +++ b/hexrd/phase_transition/texture/kernels.py @@ -21,6 +21,7 @@ 'monoclinic': 'c2h', 'orthorhombic': 'd2h', } +_AngleResult = Union[float, np.ndarray] _Symmetry = Optional[Union[str, np.ndarray]] @@ -32,12 +33,19 @@ def _symmetry_quaternions( if symmetry is None: return None if isinstance(symmetry, np.ndarray): + if len(symmetry.shape) != 2 or symmetry.shape[0] != 4: + raise ValueError(f"{symtype} symmetry must have shape (4, n)") return symmetry if not isinstance(symmetry, str): - raise TypeError(f"{symtype} symmetry must be a string or ndarray") + raise ValueError(f"{symtype} symmetry must be a string or ndarray") if symtype == 'sample': symmetry = _SAMPLE_SYMMETRY_TO_LAUE.get(symmetry.lower(), symmetry) - return rotations.quatOfLaueGroup(symmetry) + try: + return rotations.quatOfLaueGroup(symmetry) + except RuntimeError as error: + raise ValueError( + f"Invalid {symtype} symmetry: {symmetry!r}" + ) from error class SO3Kernel(ABC): @@ -51,7 +59,7 @@ class SO3Kernel(ABC): @abstractmethod def eval( self, R1: np.ndarray, R2: np.ndarray - ) -> np.ndarray: + ) -> _AngleResult: """ Evaluate kernel between two rotations. @@ -62,7 +70,7 @@ def eval( Returns ------- - numpy.ndarray + float or numpy.ndarray Kernel values """ pass @@ -187,11 +195,14 @@ def norm_constant(self) -> float: def misorientation_angle( self, R1: np.ndarray, R2: np.ndarray - ) -> np.ndarray: + ) -> _AngleResult: """ Calculate misorientation angle between rotation matrices. - Uses the formula: cos(ω) = (trace(R1^T @ R2) - 1) / 2 + Without symmetry, uses the formula: + cos(ω) = (trace(R1^T @ R2) - 1) / 2. + With crystal or sample symmetry, delegates to + hexrd.core.rotations.misorientation for symmetry-reduced angles. Parameters ---------- @@ -200,7 +211,7 @@ def misorientation_angle( Returns ------- - numpy.ndarray + float or numpy.ndarray Misorientation angles in radians, shape matches input broadcasting """ R1 = np.asarray(R1) @@ -227,7 +238,7 @@ def _symmetry_reduced_misorientation_angle( self, R1: np.ndarray, R2: np.ndarray, - ) -> np.ndarray: + ) -> _AngleResult: output_shape = np.broadcast_shapes(R1.shape[:-2], R2.shape[:-2]) R1_flat = np.broadcast_to( R1, @@ -265,7 +276,7 @@ def _misorientation_symmetries(self) -> tuple[np.ndarray, ...]: def eval( self, R1: np.ndarray, R2: np.ndarray - ) -> np.ndarray: + ) -> _AngleResult: """ Evaluate de la Vallée Poussin kernel between rotations. @@ -278,7 +289,7 @@ def eval( Returns ------- - numpy.ndarray + float or numpy.ndarray Kernel values, shape matches broadcasting of R1 and R2 Examples @@ -294,7 +305,7 @@ def eval( def eval_centered( self, R: np.ndarray, center: np.ndarray - ) -> np.ndarray: + ) -> _AngleResult: """ Evaluate kernel centered at a specific orientation. @@ -310,7 +321,7 @@ def eval_centered( Returns ------- - numpy.ndarray + float or numpy.ndarray Kernel values relative to center orientation """ return self.eval(center, R) diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py index fa393df59..3de5918b1 100644 --- a/tests/test_delavalleepoussin_kernel.py +++ b/tests/test_delavalleepoussin_kernel.py @@ -177,6 +177,30 @@ def test_large_halfwidth_warns(self): with self.assertWarns(UserWarning): DeLaValleePoussinKernel(halfwidth=np.radians(100)) + def test_invalid_symmetry_type_raises_value_error(self): + """Test that invalid symmetry types raise ValueError.""" + with self.assertRaises(ValueError): + DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry=123, + ) + + def test_invalid_symmetry_string_raises_value_error(self): + """Test that invalid symmetry strings raise ValueError.""" + with self.assertRaises(ValueError): + DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry='invalid', + ) + + def test_invalid_symmetry_array_shape_raises_value_error(self): + """Test that invalid symmetry array shapes raise ValueError.""" + with self.assertRaises(ValueError): + DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry=np.eye(3), + ) + # --- misorientation_angle tests --- def test_misorientation_angle_identity(self): From 3255be6949265814c9d23904811ccb06457630a5 Mon Sep 17 00:00:00 2001 From: Brianna Major Date: Wed, 3 Jun 2026 11:49:30 -0400 Subject: [PATCH 6/7] Remove eval_centered convenience wrapper --- hexrd/phase_transition/texture/kernels.py | 23 ----------------------- tests/test_delavalleepoussin_kernel.py | 20 -------------------- 2 files changed, 43 deletions(-) diff --git a/hexrd/phase_transition/texture/kernels.py b/hexrd/phase_transition/texture/kernels.py index f489cbaed..826504b88 100644 --- a/hexrd/phase_transition/texture/kernels.py +++ b/hexrd/phase_transition/texture/kernels.py @@ -303,29 +303,6 @@ def eval( co2 = np.cos(omega / 2.0) return self._C * co2 ** (2.0 * self._kappa) - def eval_centered( - self, R: np.ndarray, center: np.ndarray - ) -> _AngleResult: - """ - Evaluate kernel centered at a specific orientation. - - Equivalent to eval(center, R) with clearer semantics - for unimodal distributions. - - Parameters - ---------- - R : array_like - Rotation matrices to evaluate at, shape (..., 3, 3) - center : array_like - Central (modal) orientation, shape (3, 3) - - Returns - ------- - float or numpy.ndarray - Kernel values relative to center orientation - """ - return self.eval(center, R) - def __repr__(self) -> str: """String representation of kernel.""" hw_deg = np.degrees(self.halfwidth) diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py index 3de5918b1..206fa0e37 100644 --- a/tests/test_delavalleepoussin_kernel.py +++ b/tests/test_delavalleepoussin_kernel.py @@ -285,26 +285,6 @@ def test_misorientation_angle_bad_shape(self): np.eye(2), np.eye(3) ) - # --- eval_centered test --- - - def test_eval_centered_matches_eval(self): - """Test that eval_centered(R, center) == eval(center, R).""" - kernel = DeLaValleePoussinKernel( - halfwidth=np.radians(10) - ) - center = np.eye(3) - R = np.array([ - [0, -1, 0], - [1, 0, 0], - [0, 0, 1], - ], dtype=float) - - val_eval = kernel.eval(center, R) - val_centered = kernel.eval_centered(R, center) - self.assertAlmostEqual( - float(val_eval), float(val_centered), places=10 - ) - # --- Batch evaluation test --- def test_batch_evaluation(self): From 8c9ec8f2d4a9a81d343dd5c0c43bd3b6fca299e0 Mon Sep 17 00:00:00 2001 From: Brianna Major Date: Wed, 3 Jun 2026 13:37:36 -0400 Subject: [PATCH 7/7] Improve test coverage Signed-off-by: Brianna Major --- tests/test_delavalleepoussin_kernel.py | 42 ++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/tests/test_delavalleepoussin_kernel.py b/tests/test_delavalleepoussin_kernel.py index 206fa0e37..894123f74 100644 --- a/tests/test_delavalleepoussin_kernel.py +++ b/tests/test_delavalleepoussin_kernel.py @@ -4,6 +4,7 @@ import unittest from scipy.special import beta as betafn +from hexrd.core import rotations from hexrd.phase_transition.texture.kernels import DeLaValleePoussinKernel @@ -42,7 +43,7 @@ class TestDeLaValleePoussinKernel(unittest.TestCase): """Test DeLaValleePoussinKernel.""" def test_kernel_matches_mtex_ground_truth_values(self): - """Test kernel evaluation against values generated from MTEX.""" + """Test kernel evaluation against ground truth generated values.""" halfwidth = 2.0 * np.arccos( 0.5 ** (1.0 / (2.0 * MTEX_KERNEL_KAPPA)) ) @@ -64,7 +65,7 @@ def test_kernel_matches_mtex_ground_truth_values(self): ) def test_halfwidth_conversion_matches_mtex_ground_truth(self): - """Test half-width to kappa conversion against MTEX output.""" + """Test half-width to kappa conversion against ground truth values.""" kernel = DeLaValleePoussinKernel(halfwidth=MTEX_HALFWIDTH_RAD) np.testing.assert_allclose( @@ -277,6 +278,43 @@ def test_batch_misorientation_angle_uses_crystal_symmetry(self): atol=1e-8, ) + def test_misorientation_angle_accepts_quaternion_array_symmetry(self): + """Test passing precomputed quaternion symmetry matches the label.""" + R90z = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1], + ], dtype=float) + + kernel_label = DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry='d4h', + ) + kernel_array = DeLaValleePoussinKernel( + halfwidth=np.radians(10), + crystal_symmetry=rotations.quatOfLaueGroup('d4h'), + ) + + angle_label = kernel_label.misorientation_angle(np.eye(3), R90z) + angle_array = kernel_array.misorientation_angle(np.eye(3), R90z) + self.assertAlmostEqual( + float(angle_array), float(angle_label), places=10 + ) + + def test_misorientation_angle_uses_sample_symmetry(self): + """Test sample-only symmetry reduces an equivalent misorientation.""" + kernel = DeLaValleePoussinKernel( + halfwidth=np.radians(10), + sample_symmetry='orthorhombic', + ) + + # 180° rotation about z is a d2h (orthorhombic) sample operation, + # so its symmetry-reduced misorientation from identity is 0. + R180z = _rotation_about_z(np.radians(180.0)) + + angle = kernel.misorientation_angle(np.eye(3), R180z) + self.assertAlmostEqual(float(angle), 0.0, places=8) + def test_misorientation_angle_bad_shape(self): """Test that non-(3,3) inputs raise ValueError.""" kernel = DeLaValleePoussinKernel(halfwidth=np.radians(10))