Merge pull request #2164 from kif/2158_test

Address numerical precision issue on macos/arm64

src/pyFAI/test/test_geometry.py

@@ -4,7 +4,7 @@
 #    Project: Azimuthal integration
 #             https://github.com/silx-kit/pyFAI
 #
-#    Copyright (C) 2015-2023 European Synchrotron Radiation Facility, Grenoble, France
+#    Copyright (C) 2015-2024 European Synchrotron Radiation Facility, Grenoble, France
 #
 #    Principal author:       Jérôme Kieffer ([email protected])
 #
@@ -34,7 +34,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "07/05/2024"
+__date__ = "16/05/2024"
 
 import unittest
 import random
@@ -43,7 +43,6 @@
 import itertools
 import logging
 import os.path
-import platform
 
 from . import utilstest
 logger = logging.getLogger(__name__)
@@ -54,6 +53,7 @@
 from ..detectors import detector_factory
 from ..third_party import transformations
 from .utilstest import UtilsTest
+from ..utils.mathutil import allclose_mod
 import fabio
 
 
@@ -627,72 +627,89 @@ def test_array_from_unit_chi_center(self):
         self.assertTrue(179 < r4.max() < 180, "Orientation 4 upperrange matches")
 
     def test_array_from_unit_tth_corner(self):
-        r1 = self.ai1.array_from_unit(unit="2th_deg", typ="corner")
-        r2 = self.ai2.array_from_unit(unit="2th_deg", typ="corner")
-        r3 = self.ai3.array_from_unit(unit="2th_deg", typ="corner")
-        r4 = self.ai4.array_from_unit(unit="2th_deg", typ="corner")
+        r1 = self.ai1.array_from_unit(unit="2th_rad", typ="corner")
+        r2 = self.ai2.array_from_unit(unit="2th_rad", typ="corner")
+        r3 = self.ai3.array_from_unit(unit="2th_rad", typ="corner")
+        r4 = self.ai4.array_from_unit(unit="2th_rad", typ="corner")
+
+        tth1 = r1[..., 0]
+        tth2 = r2[..., 0]
+        tth3 = r3[..., 0]
+        tth4 = r4[..., 0]
+
+        chi1 = r1[..., 1]
+        chi2 = r2[..., 1]
+        chi3 = r3[..., 1]
+        chi4 = r4[..., 1]
+
+
+        sin_chi1 = numpy.sin(chi1)
+        sin_chi2 = numpy.sin(chi2)
+        sin_chi3 = numpy.sin(chi3)
+        sin_chi4 = numpy.sin(chi4)
+
+        cos_chi1 = numpy.cos(chi1)
+        cos_chi2 = numpy.cos(chi2)
+        cos_chi3 = numpy.cos(chi3)
+        cos_chi4 = numpy.cos(chi4)
+
+        # Here we use complex numbers
+        z1 = tth1*cos_chi1 + tth1*sin_chi1 *1j
+        z2 = tth2*cos_chi2 + tth2*sin_chi2 *1j
+        z3 = tth3*cos_chi3 + tth3*sin_chi3 *1j
+        z4 = tth4*cos_chi4 + tth4*sin_chi4 *1j
+
+        # the mean is not sensitive to 2pi discontinuity in azimuthal direction
+        z1 = z1.mean(axis=-1)
+        z2 = z2.mean(axis=-1)
+        z3 = z3.mean(axis=-1)
+        z4 = z4.mean(axis=-1)
+
+        self.assertFalse(numpy.allclose(z1, z2), "orientation 1,2 differ")
+        self.assertFalse(numpy.allclose(z1, z3), "orientation 1,3 differ")
+        self.assertFalse(numpy.allclose(z1, z3), "orientation 1,3 differ")
+        self.assertFalse(numpy.allclose(z1, z4), "orientation 1,4 differ")
+        self.assertFalse(numpy.allclose(z2, z3), "orientation 2,3 differ")
+        self.assertFalse(numpy.allclose(z2, z4), "orientation 2,4 differ")
+        self.assertFalse(numpy.allclose(z3, z4), "orientation 3,4 differ")
+
+        #Check that the tranformation is OK. This is with complex number thus dense & complicated !
+        self.assertTrue(numpy.allclose(z1, -numpy.fliplr(z2.conj())), "orientation 1,2 flipped")
+        self.assertTrue(numpy.allclose(z1, -z3[-1::-1, -1::-1]), "orientation 1,3 inversed")
+        self.assertTrue(numpy.allclose(z1, numpy.flipud(z4.conj())), "orientation 1,4 flipped")
+        self.assertTrue(numpy.allclose(z2, numpy.flipud(z3.conj())), "orientation 2,3 flipped")
+        self.assertTrue(numpy.allclose(z2, -z4[-1::-1, -1::-1]), "orientation 2,4 inversion")
+        self.assertTrue(numpy.allclose(z3, -numpy.fliplr(z4.conj())), "orientation 3,4 flipped")
 
-        tth1 = r1[..., 0].mean(axis=-1)
-        chi1 = r1[..., 1].mean(axis=-1) / numpy.pi
-        tth2 = r2[..., 0].mean(axis=-1)
-        chi2 = r2[..., 1].mean(axis=-1) / numpy.pi
-        tth3 = r3[..., 0].mean(axis=-1)
-        chi3 = r3[..., 1].mean(axis=-1) / numpy.pi
-        tth4 = r4[..., 0].mean(axis=-1)
-        chi4 = r4[..., 1].mean(axis=-1) / numpy.pi
-
-        self.assertFalse(numpy.allclose(tth1, tth2), "orientation 1,2 differ tth")
-        self.assertFalse(numpy.allclose(chi1, chi2), "orientation 1,2 differ chi")
-        self.assertFalse(numpy.allclose(tth1, tth3), "orientation 1,3 differ tth")
-        self.assertFalse(numpy.allclose(chi1, chi3), "orientation 1,3 differ chi")
-        self.assertFalse(numpy.allclose(tth1, tth4), "orientation 1,4 differ tth")
-        self.assertFalse(numpy.allclose(chi1, chi4), "orientation 1,4 differ chi")
-
-        self.assertTrue(numpy.allclose(tth1, numpy.fliplr(tth2)), "orientation 1,2 flipped match tth")
-        self.assertTrue(numpy.allclose(tth1, numpy.flipud(tth4)), "orientation 1,4 flipped match tth")
-        self.assertTrue(numpy.allclose(tth2, numpy.flipud(tth3)), "orientation 2,3 flipped match tth")
-        self.assertTrue(numpy.allclose(chi2, -numpy.flipud(chi3)), "orientation 2,3 flipped match chi")
-        self.assertTrue(numpy.allclose(tth1, tth3[-1::-1, -1::-1]), "orientation 1,3 inversion match tth")
-        self.assertTrue(numpy.allclose(tth2, tth4[-1::-1, -1::-1]), "orientation 2,4 inversion match tth")
-
-        # Something fishy on mac-arm64 where this test fails ... correct on all other platforms !
-        if platform.system() == "Darwin" and platform.machine()=="arm64":
-            def angular_distance(a, b, modulo):
-                return numpy.minimum((a-b)%modulo,(b-a)%modulo)
-            self.assertLess(angular_distance(chi1 + 1, -numpy.fliplr(chi2), 1).mean(), 0.0001, "orientation 1,2 flipped match chi")
-            self.assertLess(angular_distance(chi1, -numpy.flipud(chi4), 1).mean(), 0.0001, "orientation 1,4 flipped match chi")
-            self.assertLess(angular_distance(chi1+1, chi3[-1::-1, -1::-1], 1).mean(), 0.0001, "orientation 1,3 inversion match chi")
-            self.assertLess(angular_distance(chi2 + 1, chi4[-1::-1, -1::-1], 1).mean(), 0.0001, "orientation 2,4 inversion match chi")
-        else:
-            self.assertTrue(numpy.allclose(chi1 + 1, -numpy.fliplr(chi2), atol=0.0001), "orientation 1,2 flipped match chi")
-            self.assertTrue(numpy.allclose(chi1, -numpy.flipud(chi4)), "orientation 1,4 flipped match chi")
-            self.assertTrue(numpy.allclose(chi1 + 1, chi3[-1::-1, -1::-1], atol=0.0001), "orientation 1,3 inversion match chi")
-            self.assertTrue(numpy.allclose(chi2 + 1, chi4[-1::-1, -1::-1]), "orientation 2,4 inversion match chi")
 
     def test_chi(self):
+        epsilon = 6e-3
         orient = {}
         for i in range(1, 5):
-            ai = geometry.Geometry.sload({"detector":"pilatus100k", "detector_config":{"orientation":i},
-                                           "poni1":0.01, "poni2":0.01,
-                                           "wavelength":1e-10})
+            ai = geometry.Geometry.sload({"detector":"Imxpad S10", "detector_config":{"orientation":i},
+                                           "poni1":0.005, "poni2":0.005, "wavelength":1e-10})
             chi_c = ai.center_array(unit="chi_rad") / numpy.pi
-            chi_m = ai.corner_array(unit="chi_rad")[..., 0].mean(axis=-1) / numpy.pi
+            corners = ai.corner_array(unit="r_mm")
+            corners_rad = corners[..., 0]
+            corners_ang = corners[..., 1]
+            z = corners_rad*numpy.cos(corners_ang) + corners_rad*numpy.sin(corners_ang)*1j
+            chi_m = numpy.angle(z.mean(axis=-1)) / numpy.pi
 
             orient[i] = {"ai": ai, "chi_c": chi_c, "chi_m": chi_m}
 
         for o, orien in orient.items():
-            self.assertLess(numpy.median(abs(orien["chi_m"] - orien["chi_c"])), 1e-7, f"Orientation {o} matches")
+            self.assertTrue(allclose_mod(orien["chi_m"], orien["chi_c"], 2), f"Orientation {o} matches")
             ai = orien["ai"]
-            self.assertLess(numpy.median(ai.delta_array(unit="chi_rad")) / numpy.pi, 1e-3, f"Orientation {o} delta chi small #0")
-            self.assertLess(numpy.median(ai.deltaChi()) / numpy.pi, 1e-3, f"Orientation {o} delta chi small #1")
+            self.assertLess(numpy.median(ai.delta_array(unit="chi_rad")) / numpy.pi, epsilon, f"Orientation {o} delta chi small #0")
+            self.assertLess(numpy.median(ai.deltaChi()) / numpy.pi, epsilon, f"Orientation {o} delta chi small #1")
             ai.reset()
-            self.assertLess(numpy.median(ai.delta_array(unit="chi_rad")) / numpy.pi, 1e-3, f"Orientation {o} delta chi small #2")
+            self.assertLess(numpy.median(ai.delta_array(unit="chi_rad")) / numpy.pi, epsilon, f"Orientation {o} delta chi small #2")
             ai.reset()
-            self.assertLess(numpy.median(ai.deltaChi()) / numpy.pi, 1e-3, f"Orientation {o} delta chi small #3")
+            self.assertLess(numpy.median(ai.deltaChi()) / numpy.pi, epsilon, f"Orientation {o} delta chi small #3")
             ai.reset()
             chiArray = ai.chiArray() / numpy.pi
             chi_center = orien["chi_c"]
-            self.assertTrue(numpy.allclose(chiArray, chi_center, atol=1e-5), f"Orientation {o} chiArray == center_array(chi)")
+            self.assertTrue(allclose_mod(chiArray, chi_center), f"Orientation {o} chiArray == center_array(chi)")
 
 
 class TestOrientation2(unittest.TestCase):
@@ -701,7 +718,7 @@ class TestOrientation2(unittest.TestCase):
     @classmethod
     def setUpClass(cls) -> None:
         super(TestOrientation2, cls).setUpClass()
-        p = detector_factory("Pilatus100k")
+        p = detector_factory("pilatus100k")
         c = p.get_pixel_corners()
         d1 = c[..., 1].max()
         d2 = c[..., 2].max()

src/pyFAI/test/test_utils_mathutil.py

@@ -4,7 +4,7 @@
 #    Project: Azimuthal integration
 #             https://github.com/silx-kit/pyFAI
 #
-#    Copyright (C) 2015-2021 European Synchrotron Radiation Facility, Grenoble, France
+#    Copyright (C) 2015-2024 European Synchrotron Radiation Facility, Grenoble, France
 #
 #    Principal author:       Jérôme Kieffer ([email protected])
 #
@@ -32,7 +32,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "20/02/2024"
+__date__ = "16/05/2024"
 
 import unittest
 import numpy
@@ -157,6 +157,11 @@ def test_is_far_from_group_cython(self):
         res = is_far_from_group_cython(pt, pts, dst2)
         self.assertEqual(ref, res, "cython implementation matches *is_far_from_group*")
 
+    def test_allclose_mod(self):
+        from ..utils.mathutil import allclose_mod
+        self.assertTrue(allclose_mod(numpy.arctan2(+1e-10, -1), numpy.arctan2(-1e-10, -1)),"angles matches modulo 2pi")
+
+
 def suite():
     loader = unittest.defaultTestLoader.loadTestsFromTestCase
     testsuite = unittest.TestSuite()

src/pyFAI/utils/mathutil.py

@@ -4,7 +4,7 @@
 #    Project: Fast Azimuthal integration
 #             https://github.com/silx-kit/pyFAI
 #
-#    Copyright (C) 2017-2018 European Synchrotron Radiation Facility, Grenoble, France
+#    Copyright (C) 2017-2024 European Synchrotron Radiation Facility, Grenoble, France
 #
 #    Principal author:       Jérôme Kieffer ([email protected])
 #
@@ -34,7 +34,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "20/02/2024"
+__date__ = "16/05/2024"
 __status__ = "production"
 
 import logging
@@ -394,6 +394,7 @@ def unBinning(*args, **kwargs):
     return unbinning(*args, **kwargs)
 
 
+
 def shift_fft(input_img, shift_val, method="fft"):
     """Do shift using FFTs
 
@@ -928,3 +929,13 @@ def interp_filter(ary, out=None):
     out[mask_invalid] = numpy.interp(x[mask_invalid], x[mask_valid], ary[mask_valid],
                                      left=first, right=last)
     return out
+
+
+def allclose_mod(a, b, modulo=2*numpy.pi, **kwargs):
+    """Returns True if the two arrays a & b are equal within the given
+    tolerance modulo `modulo`; False otherwise.
+
+    Thanks to "Serguei Sokol" <[email protected]>
+    """
+    di = numpy.minimum((a-b)%modulo, (b-a)%modulo)
+    return numpy.allclose(modulo*0.5, (di+modulo*0.5), **kwargs)