tests.test_diagnostics

  1"""> PyDRex: tests for texture diagnostics."""
  2
  3import numpy as np
  4import pytest
  5from numpy import random as rn
  6from scipy.spatial.transform import Rotation
  7
  8from pydrex import diagnostics as _diagnostics
  9from pydrex import geometry as _geo
 10from pydrex import stats as _stats
 11
 12
 13class TestElasticityComponents:
 14    """Test symmetry decomposition of elastic tensors."""
 15
 16    def test_olivine_Browaeys2004(self):
 17        C = np.array(
 18            [
 19                [192, 66, 60, 0, 0, 0],
 20                [66, 160, 56, 0, 0, 0],
 21                [60, 56, 272, 0, 0, 0],
 22                [0, 0, 0, 60, 0, 0],
 23                [0, 0, 0, 0, 62, 0],
 24                [0, 0, 0, 0, 0, 49],
 25            ]
 26        )
 27        out = _diagnostics.elasticity_components([C])
 28        # FIXME: How do they get 15.2% for hexagonal? It isn't the ratio of the norms
 29        # nor the ratio of the squared norms...
 30        expected = {
 31            "bulk_modulus": 109.8,
 32            "shear_modulus": 63.7,
 33            "percent_anisotropy": 20.7,
 34            # "percent_hexagonal": 15.2,
 35            # "percent_tetragonal": 0.625,
 36            # "percent_orthorhombic": 4.875,
 37            "percent_monoclinic": 0,
 38            "percent_triclinic": 0,
 39        }
 40        for k, v in expected.items():
 41            assert np.isclose(out[k][0], v, atol=0.1, rtol=0), f"{k}: {out[k]} != {v}"
 42        np.testing.assert_allclose(out["hexagonal_axis"][0], [0, 0, 1])
 43
 44    def test_enstatite_Browaeys2004(self):
 45        C = np.array(
 46            [
 47                [225, 54, 72, 0, 0, 0],
 48                [54, 214, 53, 0, 0, 0],
 49                [72, 53, 178, 0, 0, 0],
 50                [0, 0, 0, 78, 0, 0],
 51                [0, 0, 0, 0, 82, 0],
 52                [0, 0, 0, 0, 0, 76],
 53            ]
 54        )
 55        out = _diagnostics.elasticity_components([C])
 56        # FIXME: Test remaining percentages when I figure out how they get the values.
 57        expected = {
 58            "bulk_modulus": 108.3,
 59            "shear_modulus": 76.4,
 60            "percent_anisotropy": 9.2,
 61            # "percent_hexagonal": 4.3,
 62            # "percent_tetragonal": ?,  # + ortho = 4.9
 63            # "percent_orthorhombic": ?,  # + tetra = 4.9
 64            "percent_monoclinic": 0,
 65            "percent_triclinic": 0,
 66        }
 67        for k, v in expected.items():
 68            assert np.isclose(out[k][0], v, atol=0.1, rtol=0), f"{k}: {out[k]} != {v}"
 69        np.testing.assert_allclose(out["hexagonal_axis"][0], [0, 0, 1])
 70
 71
 72class TestSymmetryPGR:
 73    """Test Point-Girdle-Random (eigenvalue) symmetry diagnostics."""
 74
 75    def test_pointX(self):
 76        """Test diagnostics of point symmetry aligned to the X axis."""
 77        # Initial orientations within 10°.
 78        orientations = (
 79            Rotation.from_rotvec(
 80                np.stack(
 81                    [
 82                        [0, x * np.pi / 18 - np.pi / 36, x * np.pi / 18 - np.pi / 36]
 83                        for x in rn.default_rng().random(100)
 84                    ]
 85                )
 86            )
 87            .inv()
 88            .as_matrix()
 89        )
 90        np.testing.assert_allclose(
 91            _diagnostics.symmetry_pgr(orientations, axis="a"), (1, 0, 0), atol=0.05
 92        )
 93
 94    def test_random(self):
 95        """Test diagnostics of random grain orientations."""
 96        orientations = Rotation.random(1000).as_matrix()
 97        np.testing.assert_allclose(
 98            _diagnostics.symmetry_pgr(orientations, axis="a"), (0, 0, 1), atol=0.15
 99        )
100
101    def test_girdle(self):
102        """Test diagnostics of girdled orientations."""
103        rng = rn.default_rng()
104        a = np.zeros(1000)
105        b = np.zeros(1000)
106        c = rng.normal(0, 1.0, size=1000)
107        d = rng.normal(0, 1.0, size=1000)
108        orientations = Rotation.from_quat(np.column_stack([a, b, c, d])).as_matrix()
109        np.testing.assert_allclose(
110            _diagnostics.symmetry_pgr(orientations, axis="a"), (0, 1, 0), atol=0.1
111        )
112
113
114class TestVolumeWeighting:
115    """Tests for volumetric resampling of orientation data."""
116
117    def test_output_shape(self):
118        """Test that we get the correct output shape."""
119        orientations = [
120            Rotation.random(1000).as_matrix(),
121            Rotation.random(1000).as_matrix(),
122        ]
123        fractions = [np.full(1000, 1 / 1000), np.full(1000, 1 / 1000)]
124        new_orientations, new_fractions = _stats.resample_orientations(
125            orientations, fractions
126        )
127        np.testing.assert_array_equal(
128            np.asarray(orientations).shape, new_orientations.shape
129        )
130        np.testing.assert_array_equal(np.asarray(fractions).shape, new_fractions.shape)
131        new_orientations, new_fractions = _stats.resample_orientations(
132            orientations, fractions, n_samples=500
133        )
134        assert new_orientations.shape[1] == 500
135        assert new_fractions.shape[1] == 500
136
137    def test_upsample(self):
138        """Test upsampling of the raw orientation data."""
139        orientations = (
140            Rotation.from_rotvec(
141                [
142                    [0, 0, 0],
143                    [0, 0, np.pi / 6],
144                    [np.pi / 6, 0, 0],
145                ]
146            )
147            .inv()
148            .as_matrix()
149        )
150        fractions = np.array([0.25, 0.6, 0.15])
151        new_orientations, _ = _stats.resample_orientations(
152            [orientations],
153            [fractions],
154            25,
155        )
156        assert np.all(a in orientations for a in new_orientations[0])
157
158    def test_downsample(self):
159        """Test downsampling of orientation data."""
160        orientations = (
161            Rotation.from_rotvec(
162                [
163                    [0, 0, 0],
164                    [0, 0, np.pi / 6],
165                    [np.pi / 6, 0, 0],
166                ]
167            )
168            .inv()
169            .as_matrix()
170        )
171        fractions = np.array([0.25, 0.6, 0.15])
172        new_orientations = _stats.resample_orientations(
173            [orientations],
174            [fractions],
175            2,
176        )
177        assert np.all(a in orientations for a in new_orientations[0])
178
179    def test_common_input_errors(self):
180        """Test that exceptions are raised for bad input data."""
181        orientations = Rotation.from_rotvec(
182            [
183                [0, 0, 0],
184                [0, 0, np.pi / 6],
185                [np.pi / 6, 0, 0],
186            ]
187        )
188        fractions = np.array([0.25, 0.6, 0.15])
189        with pytest.raises(ValueError):
190            # SciPy Rotation instance is not valid input.
191            _stats.resample_orientations(orientations, fractions)
192            # Input must be a stack of orientations.
193            _stats.resample_orientations(orientations.as_matrix(), fractions)
194            # Input must be a stack of fractions.
195            _stats.resample_orientations([orientations.as_matrix()], fractions)
196            # First two dimensions of inputs must match.
197            _stats.resample_orientations([orientations.as_matrix()], [fractions[:-1]])
198
199        # This is the proper way to do it:
200        _stats.resample_orientations([orientations.as_matrix()], [fractions])
201
202
203class TestBinghamStats:
204    """Tests for antipodally symmetric (bingham) statistics."""
205
206    def test_average_0(self):
207        """Test bingham average of vectors aligned to the reference frame."""
208        orientations = Rotation.from_rotvec([[0, 0, 0]] * 10).inv().as_matrix()
209        a_mean = _diagnostics.bingham_average(orientations, axis="a")
210        b_mean = _diagnostics.bingham_average(orientations, axis="b")
211        c_mean = _diagnostics.bingham_average(orientations, axis="c")
212        np.testing.assert_array_equal(a_mean, [1, 0, 0])
213        np.testing.assert_array_equal(b_mean, [0, 1, 0])
214        np.testing.assert_array_equal(c_mean, [0, 0, 1])
215
216    def test_average_twopoles90Z(self):
217        """Test bingham average of vectors rotated by ±90° around Z."""
218        orientations = (
219            Rotation.from_rotvec(
220                [
221                    [0, 0, -np.pi / 2],
222                    [0, 0, np.pi / 2],
223                ]
224            )
225            .inv()
226            .as_matrix()
227        )
228        a_mean = _diagnostics.bingham_average(orientations, axis="a")
229        b_mean = _diagnostics.bingham_average(orientations, axis="b")
230        c_mean = _diagnostics.bingham_average(orientations, axis="c")
231        np.testing.assert_array_equal(a_mean, [0, 1, 0])
232        np.testing.assert_array_equal(b_mean, [1, 0, 0])
233        np.testing.assert_array_equal(c_mean, [0, 0, 1])
234
235    def test_average_spread10X(self):
236        """Test bingham average of vectors spread within 10° of the ±X-axis."""
237        orientations = (
238            Rotation.from_rotvec(
239                np.stack(
240                    [
241                        [0, x * np.pi / 18 - np.pi / 36, x * np.pi / 18 - np.pi / 36]
242                        for x in rn.default_rng().random(100)
243                    ]
244                )
245            )
246            .inv()
247            .as_matrix()
248        )
249        a_mean = _diagnostics.bingham_average(orientations, axis="a")
250        b_mean = _diagnostics.bingham_average(orientations, axis="b")
251        c_mean = _diagnostics.bingham_average(orientations, axis="c")
252        np.testing.assert_allclose(np.abs(a_mean), [1.0, 0, 0], atol=np.sin(np.pi / 18))
253        np.testing.assert_allclose(np.abs(b_mean), [0, 1.0, 0], atol=np.sin(np.pi / 18))
254        np.testing.assert_allclose(np.abs(c_mean), [0, 0, 1.0], atol=np.sin(np.pi / 18))
255
256
257class TestMIndex:
258    """Tests for the M-index texture strength diagnostic."""
259
260    # TODO: Add tests for other (non-orthorhombic) lattice symmetries.
261
262    def test_texture_uniform_ortho(self, seed):
263        """Test with random (uniform distribution) orthorhombic grain orientations."""
264        orientations = Rotation.random(1000, random_state=seed).as_matrix()
265        assert np.isclose(
266            _diagnostics.misorientation_index(
267                orientations, _geo.LatticeSystem.orthorhombic
268            ),
269            0.05,
270            atol=1e-2,
271            rtol=0,
272        )
273
274    def test_texture_spread10X_ortho(self, seed):
275        """Test for orthorhombic grains spread within 10° of the ±X axis."""
276        orientations = (
277            Rotation.from_rotvec(
278                np.stack(
279                    [
280                        [0, x * np.pi / 18 - np.pi / 36, x * np.pi / 18 - np.pi / 36]
281                        for x in rn.default_rng(seed=seed).random(1000)
282                    ]
283                )
284            )
285            .inv()
286            .as_matrix()
287        )
288        assert np.isclose(
289            _diagnostics.misorientation_index(
290                orientations, _geo.LatticeSystem.orthorhombic
291            ),
292            0.99,
293            atol=1e-2,
294            rtol=0,
295        )
296
297    def test_texture_spread45X_ortho(self, seed):
298        """Test for orthorhombic grains spread within 45° of the ±X axis."""
299        orientations = Rotation.from_rotvec(
300            np.stack(
301                [
302                    [0, x * np.pi / 2 - np.pi / 4, x * np.pi / 2 - np.pi / 4]
303                    for x in rn.default_rng(seed=seed).random(1000)
304                ]
305            )
306        ).as_matrix()
307        assert np.isclose(
308            _diagnostics.misorientation_index(
309                orientations, _geo.LatticeSystem.orthorhombic
310            ),
311            0.81,
312            atol=1e-2,
313            rtol=0,
314        )
315
316    def test_textures_increasing_ortho(self, seed):
317        """Test M-index for textures of increasing strength."""
318        M_vals = np.empty(4)
319        factors = (2, 3, 4, 5)
320        for i, factor in enumerate(factors):
321            orientations = Rotation.from_rotvec(
322                np.stack(
323                    [
324                        [
325                            0,
326                            x * np.pi / factor - np.pi / factor / 2,
327                            x * np.pi / factor - np.pi / factor / 2,
328                        ]
329                        for x in rn.default_rng(seed=seed).random(1000)
330                    ]
331                )
332            ).as_matrix()
333            M_vals[i] = _diagnostics.misorientation_index(
334                orientations, _geo.LatticeSystem.orthorhombic
335            )
336        assert np.all(
337            np.diff(M_vals) >= 0
338        ), f"M-index values {M_vals} are not increasing"
339
340    def test_texture_girdle_ortho(self, seed):
341        """Test M-index for girdled texture."""
342        rng = rn.default_rng(seed=seed)
343        a = np.zeros(1000)
344        b = np.zeros(1000)
345        c = rng.normal(0, 1.0, size=1000)
346        d = rng.normal(0, 1.0, size=1000)
347        orientations = Rotation.from_quat(np.column_stack([a, b, c, d])).as_matrix()
348        assert np.isclose(
349            _diagnostics.misorientation_index(
350                orientations, _geo.LatticeSystem.orthorhombic
351            ),
352            0.67,
353            atol=1e-2,
354        )