reference Data for ODF sampling (hybridIA) + skeleton of a test

python/damask/_rotation.py

@@ -1,6 +1,8 @@
 import numpy as np
 
 from . import mechanics
+from . import util
+from . import grid_filters
 
 _P = -1
 
@@ -647,13 +649,56 @@ class Rotation:
 
         return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize()
 
-    # for compatibility (old names do not follow convention)
-    fromEulers     = from_Eulers
-    fromQuaternion = from_quaternion
-    asAxisAngle    = as_axis_angle
+    # for compatibility
     __mul__        = __matmul__
 
 
+    @staticmethod
+    def from_ODF(weights,Eulers,N=500,degrees=True,fractions=True,seed=None):
+        """
+        Sample discrete orientations from a binned ODF.
+
+        References
+        ----------
+        P. Eisenlohr, F. Roters, Computational Materials Science 42(4), 670-678, 2008
+        https://doi.org/10.1016/j.commatsci.2007.09.015
+
+        Parameters
+        ----------
+        fname : file, str, or pathlib.Path
+            ODF file containing normalized probability (labeled 'intensity')
+            on a grid in Euler space (labeled 'euler').
+        N : integer, optional
+            Number of discrete orientations to be sampled from the given ODF.
+            Defaults to 500.
+        degrees : boolean, optional
+            Euler grid values are in degrees. Defaults to True.
+        fractions : boolean, optional
+            ODF values correspond to volume fractions, not probability density.
+            Defaults to True.
+        seed: {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
+            A seed to initialize the BitGenerator. Defaults to None, i.e. unpredictable entropy
+            will be pulled from the OS.
+
+        Notes
+        -----
+        Explain here the different things that need to be considered
+
+        """
+        def _dg(eu,deg):
+            """Return infinitesimal Euler space volume of bin(s)."""
+            Eulers_sorted = eu[np.lexsort((eu[:,0],eu[:,1],eu[:,2]))]
+            steps,size,_ = grid_filters.cell_coord0_gridSizeOrigin(Eulers_sorted)
+            delta = np.radians(size/steps) if deg else size/steps
+            return delta[0]*2.0*np.sin(delta[1]/2.0)*delta[2] / 8.0 / np.pi**2 * np.sin(np.radians(eu[:,1]) if deg else eu[:,1])
+
+        dg = 1.0 if fractions else _dg(Eulers,degrees)
+        dV_V = dg * np.maximum(0.0,weights.squeeze())
+        orientations = Rotation.from_Eulers(Eulers[util.hybrid_IA(dV_V,N,seed)],degrees)
+
+        return orientations
+
+
     @staticmethod
     def from_spherical_component(center,sigma,N=500,degrees=True,seed=None):
         """

python/tests/test_Rotation.py

@@ -1,11 +1,11 @@
-import os
-
 import pytest
 import numpy as np
 from scipy import stats
 
 from damask import Rotation
+from damask import Table
 from damask import _rotation
+from damask import grid_filters
 
 n = 1000
 atol=1.e-4
@@ -13,7 +13,7 @@ atol=1.e-4
 @pytest.fixture
 def reference_dir(reference_dir_base):
     """Directory containing reference results."""
-    return os.path.join(reference_dir_base,'Rotation')
+    return reference_dir_base/'Rotation'
 
 @pytest.fixture
 def set_of_rotations(set_of_quaternions):
@@ -943,3 +943,39 @@ class TestRotation:
         sigma_out = np.degrees(np.std(dist))
         p = np.average(p)
         assert (.9 < sigma/sigma_out < 1.1) and p > 1e-2, f'{sigma/sigma_out},{p}'
+
+
+    @pytest.mark.parametrize('fractions',[True,False])
+    @pytest.mark.parametrize('degrees',[True,False])
+    @pytest.mark.parametrize('N',[2**13,2**14,2**15])
+    def test_ODF_cell(self,reference_dir,fractions,degrees,N):
+        steps = np.array([144,36,36])
+        limits = np.array([360.,90.,90.])
+        rng = tuple(zip(np.zeros(3),limits))
+
+        weights = Table.load(reference_dir/'ODF_experimental_cell.txt').get('intensity').flatten()
+        Eulers = grid_filters.cell_coord0(steps,limits)
+        Eulers = np.radians(Eulers) if not degrees else Eulers
+
+        Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees,fractions).as_Eulers(True)
+        weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
+
+        if fractions: assert np.sqrt(((weights_r - weights) ** 2).mean()) < 4
+
+    @pytest.mark.parametrize('degrees',[True,False])
+    @pytest.mark.parametrize('N',[2**13,2**14,2**15])
+    def test_ODF_node(self,reference_dir,degrees,N):
+        steps = np.array([144,36,36])
+        limits = np.array([360.,90.,90.])
+        rng = tuple(zip(np.zeros(3)-limits/steps*.5,limits-limits/steps*.5))
+
+        weights = Table.load(reference_dir/'ODF_experimental.txt').get('intensity')
+        weights = weights.reshape(steps+1,order='F')[:-1,:-1,:-1].reshape(-1,order='F')
+
+        Eulers = grid_filters.node_coord0(steps,limits)[:-1,:-1,:-1]
+        Eulers = np.radians(Eulers) if not degrees else Eulers
+
+        Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees).as_Eulers(True)
+        weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
+
+        assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5