Texture unimodal odf

hexrd/phase_transition/texture/__init__.py

@@ -3,9 +3,23 @@
     SO3Kernel,
 )
 from hexrd.phase_transition.texture.uniform_odf import UniformODF
+from hexrd.phase_transition.texture.unimodal_odf import UnimodalODF
+from hexrd.phase_transition.texture.evaluation import (
+    validate_orientations,
+    eval_odf,
+    eval_odf_batch,
+    eval_at_identity,
+    eval_random_orientations,
+)
 
 __all__ = [
     'DeLaValleePoussinKernel',
     'SO3Kernel',
     'UniformODF',
+    'UnimodalODF',
+    'validate_orientations',
+    'eval_odf',
+    'eval_odf_batch',
+    'eval_at_identity',
+    'eval_random_orientations',
 ]

hexrd/phase_transition/texture/evaluation.py

@@ -0,0 +1,212 @@
+"""
+ODF Evaluation Functions
+
+Provides common evaluation functionality for orientation distribution functions.
+This module serves as a foundation for more complex ODF implementations.
+"""
+
+import numpy as np
+from scipy.spatial.transform import Rotation
+
+
+def validate_orientations(orientations):
+    """
+    Validate and normalize orientation input.
+
+    Parameters
+    ----------
+    orientations : array_like
+        Orientation matrices, expected shape (..., 3, 3)
+
+    Returns
+    -------
+    numpy.ndarray
+        Validated orientation matrices
+
+    Raises
+    ------
+    ValueError
+        If orientations don't have proper shape or aren't valid rotation matrices
+    """
+    orientations = np.asarray(orientations)
+
+    # Check basic shape requirements
+    if orientations.ndim < 2:
+        raise ValueError("Orientations must have at least 2 dimensions")
+
+    if orientations.shape[-2:] != (3, 3):
+        raise ValueError(
+            f"Orientations must have shape (..., 3, 3), got {orientations.shape}"
+        )
+
+    # Optional: Check if matrices are approximately orthogonal
+    # For performance, we'll skip this validation by default
+    # but could add a 'validate_rotation' parameter in the future
+
+    return orientations
+
+
+def eval_odf(odf, orientations, validate_input=True):
+    """
+    Generic ODF evaluation function.
+
+    Provides a unified interface for evaluating orientation distribution functions
+    at specified orientations. Supports both single orientations and batch evaluation.
+
+    Parameters
+    ----------
+    odf : ODF object
+        Any ODF object with an eval() method (e.g., UniformODF, UnimodalODF)
+    orientations : array_like
+        Orientation matrices of shape (..., 3, 3). Can be:
+        - Single 3x3 matrix
+        - Array of shape (N, 3, 3) for N orientations
+        - Higher dimensional arrays with shape (..., 3, 3)
+    validate_input : bool, optional
+        Whether to validate input orientations, default True
+
+    Returns
+    -------
+    numpy.ndarray
+        ODF values evaluated at the given orientations.
+        Shape matches the leading dimensions of input orientations.
+
+    Examples
+    --------
+    >>> from hexrd.phase_transition.texture import UniformODF, eval_odf
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>>
+    >>> # Single orientation
+    >>> R = np.eye(3)
+    >>> value = eval_odf(odf, R)  # scalar result
+    >>>
+    >>> # Multiple orientations
+    >>> Rs = np.array([np.eye(3), np.eye(3)])  # shape (2, 3, 3)
+    >>> values = eval_odf(odf, Rs)  # shape (2,)
+    """
+    # Validate ODF object
+    if not hasattr(odf, 'eval'):
+        raise TypeError(
+            f"ODF object {type(odf).__name__} must have an eval() method"
+        )
+
+    # Validate and process orientations
+    if validate_input:
+        orientations = validate_orientations(orientations)
+
+    # Delegate to ODF's eval method
+    return odf.eval(orientations)
+
+
+def eval_odf_batch(odf, orientations, chunk_size=10000, validate_input=True):
+    """
+    Evaluate ODF on large batches of orientations with memory management.
+
+    For very large orientation datasets, this function processes orientations
+    in chunks to avoid memory issues while maintaining vectorization benefits.
+
+    Parameters
+    ----------
+    odf : ODF object
+        ODF object with eval() method
+    orientations : array_like
+        Large array of orientation matrices, shape (N, 3, 3)
+    chunk_size : int, optional
+        Number of orientations to process per chunk, default 10000
+    validate_input : bool, optional
+        Whether to validate input orientations, default True
+
+    Returns
+    -------
+    numpy.ndarray
+        ODF values, shape (N,)
+
+    Examples
+    --------
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>> # Large batch of orientations
+    >>> Rs = np.random.normal(size=(50000, 3, 3))  # Not actual rotations!
+    >>> # In practice, use proper rotation matrices
+    >>> values = eval_odf_batch(odf, Rs, chunk_size=5000)
+    """
+    if validate_input:
+        orientations = validate_orientations(orientations)
+
+    # Handle single orientation or small batches normally
+    if orientations.ndim == 2 or orientations.shape[0] <= chunk_size:
+        return eval_odf(odf, orientations, validate_input=False)
+
+    # Process large batches in chunks
+    n_orientations = orientations.shape[0]
+    results = []
+
+    for start_idx in range(0, n_orientations, chunk_size):
+        end_idx = min(start_idx + chunk_size, n_orientations)
+        chunk = orientations[start_idx:end_idx]
+        chunk_results = eval_odf(odf, chunk, validate_input=False)
+        results.append(chunk_results)
+
+    return np.concatenate(results)
+
+
+def eval_at_identity(odf):
+    """
+    Evaluate ODF at the identity orientation.
+
+    Convenience function for getting ODF value at the identity rotation.
+    Useful for normalization checks and texture strength calculations.
+
+    Parameters
+    ----------
+    odf : ODF object
+        ODF to evaluate
+
+    Returns
+    -------
+    float
+        ODF value at identity orientation
+
+    Examples
+    --------
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>> value = eval_at_identity(odf)  # 1.0 for a uniform ODF (MRD)
+    """
+    identity = np.eye(3)
+    return float(eval_odf(odf, identity, validate_input=False))
+
+
+def eval_random_orientations(odf, n_orientations=1000, seed=None):
+    """
+    Evaluate ODF at random orientations for statistical analysis.
+
+    Generates random rotation matrices and evaluates the ODF at these
+    orientations. Useful for normalization checks and Monte Carlo integration.
+
+    Parameters
+    ----------
+    odf : ODF object
+        ODF to evaluate
+    n_orientations : int, optional
+        Number of random orientations to generate, default 1000
+    seed : int, optional
+        Random seed for reproducibility
+
+    Returns
+    -------
+    tuple of (numpy.ndarray, numpy.ndarray)
+        orientations : Random rotation matrices, shape (n_orientations, 3, 3)
+        values : ODF values at these orientations, shape (n_orientations,)
+
+    Examples
+    --------
+    >>> odf = UniformODF('oh', 'triclinic')
+    >>> orientations, values = eval_random_orientations(odf, n_orientations=100)
+    >>> print(f"Mean ODF value: {np.mean(values)}")  # Should be ~1.0 (MRD)
+    """
+    # Haar-uniform random rotations on SO(3).
+    orientations = Rotation.random(n_orientations, random_state=seed).as_matrix()
+
+    # Evaluate ODF at these orientations
+    values = eval_odf(odf, orientations, validate_input=False)
+
+    return orientations, values

hexrd/phase_transition/texture/kernels.py

@@ -1,7 +1,7 @@
 """
 SO(3) Kernel Functions for Texture Analysis
 
-Implements radial basis functions on the rotation group SO(3) used for 
+Implements radial basis functions on the rotation group SO(3) used for
 constructing smooth orientation distribution functions.
 """
 
@@ -48,6 +48,22 @@ def _symmetry_quaternions(
         ) from error
 
 
+def _is_trivial_symmetry(quats: Optional[np.ndarray]) -> bool:
+    """
+    Return True if a symmetry group performs no reduction.
+
+    A group is trivial when it is absent (None) or contains only the
+    identity rotation (e.g. triclinic). The identity quaternion has a
+    scalar part w = ±1 with a zero vector part, while every proper,
+    non-identity rotation has |w| < 1. A group whose elements are all
+    ±identity therefore leaves misorientation angles unchanged, so it can
+    safely take the fast (no-symmetry) evaluation path.
+    """
+    if quats is None:
+        return True
+    return bool(np.allclose(np.abs(quats[0]), 1.0))
+
+
 class SO3Kernel(ABC):
     """
     Abstract base class for kernels on the SO(3) rotation group.
@@ -161,21 +177,53 @@ def __init__(
 
         # 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(
+        self._crystal_symmetry_quats = _symmetry_quaternions(
             crystal_symmetry,
             symtype='crystal',
         )
-        self._sample_symmetry = _symmetry_quaternions(
+        self._sample_symmetry_quats = _symmetry_quaternions(
             sample_symmetry,
             symtype='sample',
         )
 
+        # Retain the original symmetry labels (when given as strings) so the
+        # kernel can act as the single source of truth for symmetry. When
+        # symmetry is supplied as a quaternion array, no label is available.
+        self._crystal_symmetry_label = (
+            crystal_symmetry if isinstance(crystal_symmetry, str) else None
+        )
+        self._sample_symmetry_label = (
+            sample_symmetry if isinstance(sample_symmetry, str) else None
+        )
+
+    @property
+    def crystal_symmetry(self) -> Optional[str]:
+        """
+        str or None: Crystal symmetry label, if built from a label.
+
+        Returns None when the symmetry was supplied directly as a
+        quaternion array, even if symmetry reduction is active. Use
+        ``has_symmetry`` to test whether reduction is enabled.
+        """
+        return self._crystal_symmetry_label
+
+    @property
+    def sample_symmetry(self) -> Optional[str]:
+        """
+        str or None: Sample symmetry label, if built from a label.
+
+        Returns None when the symmetry was supplied directly as a
+        quaternion array, even if symmetry reduction is active. Use
+        ``has_symmetry`` to test whether reduction is enabled.
+        """
+        return self._sample_symmetry_label
+
     @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
+        """bool: Whether non-trivial symmetry reduction is enabled."""
+        return not (
+            _is_trivial_symmetry(self._crystal_symmetry_quats)
+            and _is_trivial_symmetry(self._sample_symmetry_quats)
         )
 
     @property
@@ -251,6 +299,12 @@ def _symmetry_reduced_misorientation_angle(
 
         angles = np.empty(R1_flat.shape[0], dtype=float)
         symmetries = self._misorientation_symmetries()
+        # TODO(perf): this per-pair Python loop calls rotations.misorientation
+        # once per orientation and is ~1400x slower than the no-symmetry path
+        # for large grids (e.g. ~2 s for 3000 cubic-symmetry evaluations).
+        # Vectorize by applying the symmetry quaternion group with broadcasting
+        # and taking the minimum angle. Revisit on the texture-normalization
+        # branch before evaluating ODFs on dense (1e5-1e6) grids.
         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)
@@ -266,12 +320,12 @@ def _symmetry_reduced_misorientation_angle(
         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 = self._crystal_symmetry_quats
+        if crystal_symmetry is None and self._sample_symmetry_quats 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,)
+        if self._sample_symmetry_quats is not None:
+            symmetries = symmetries + (self._sample_symmetry_quats,)
         return symmetries
 
     def eval(