fix: default to brute-force computation when mud outside of range

Author: yucongalicechen

news/fastcalto7.rst

@@ -1,10 +1,11 @@
 **Added:**
 
-* Fast calculation support up to muD = 7
+* Fast calculation supports values up to muD = 7
 
 **Changed:**
 
-* Clarified error message for fast calculation to explicitly states the invalid muD value
+* Default to brute-force computation when muD < 0.5 or > 7.
+* Print a warning message instead of error, explicitly stating the input muD value
 
 **Deprecated:**
 

src/diffpy/labpdfproc/functions.py

@@ -1,4 +1,5 @@
 import math
+import warnings
 from pathlib import Path
 
 import numpy as np
@@ -16,7 +17,7 @@
 CVE_METHODS = ["brute_force", "polynomial_interpolation"]
 
 # Pre-computed datasets for polynomial interpolation (fast calculation)
-MUD_LIST = [0.5, 1, 2, 3, 4, 5, 6, 7]
+MUD_LIST = np.array([0.5, 1, 2, 3, 4, 5, 6, 7])
 CWD = Path(__file__).parent.resolve()
 MULS = np.loadtxt(CWD / "data" / "inverse_cve.xy")
 COEFFICIENT_LIST = np.array(
@@ -74,7 +75,6 @@ def _get_entry_exit_coordinates(self, coordinate, angle):
         ----------
         coordinate : tuple of floats
             The coordinates of the grid point.
-
         angle : float
             The angle in degrees.
 
@@ -90,9 +90,7 @@ def _get_entry_exit_coordinates(self, coordinate, angle):
         angle = math.radians(angle)
         xgrid = coordinate[0]
         ygrid = coordinate[1]
-
         entry_point = (-math.sqrt(self.radius**2 - ygrid**2), ygrid)
-
         if not math.isclose(angle, math.pi / 2, abs_tol=epsilon):
             b = ygrid - xgrid * math.tan(angle)
             a = math.tan(angle)
@@ -107,7 +105,6 @@ def _get_entry_exit_coordinates(self, coordinate, angle):
                 exit_point = (xexit_root1, yexit_root1)
         else:
             exit_point = (xgrid, math.sqrt(self.radius**2 - xgrid**2))
-
         return entry_point, exit_point
 
     def _get_path_length(self, grid_point, angle):
@@ -119,7 +116,6 @@ def _get_path_length(self, grid_point, angle):
         ----------
         grid_point : double of floats
             The coordinate inside the circle.
-
         angle : float
             The angle of the output beam in degrees.
 
@@ -129,7 +125,6 @@ def _get_path_length(self, grid_point, angle):
             The tuple containing three floats,
             which are the total distance, entry distance and exit distance.
         """
-
         # move angle a tad above zero if it is zero
         # to avoid it having the wrong sign due to some rounding error
         angle_delta = 0.000001
@@ -181,7 +176,9 @@ def _cve_brute_force(input_pattern, mud):
     Assume mu=mud/2, given that the same mu*D yields the same cve and D/2=1.
     """
     mu_sample_invmm = mud / 2
-    abs_correction = Gridded_circle(mu=mu_sample_invmm)
+    abs_correction = Gridded_circle(
+        n_points_on_diameter=N_POINTS_ON_DIAMETER, mu=mu_sample_invmm
+    )
     distances, muls = [], []
     for angle in TTH_GRID:
         abs_correction.set_distances_at_angle(angle)
@@ -191,7 +188,6 @@ def _cve_brute_force(input_pattern, mud):
     distances = np.array(distances) / abs_correction.total_points_in_grid
     muls = np.array(muls) / abs_correction.total_points_in_grid
     cve = 1 / muls
-
     cve_do = DiffractionObject(
         xarray=TTH_GRID,
         yarray=cve,
@@ -206,31 +202,21 @@ def _cve_brute_force(input_pattern, mud):
 
 def _cve_polynomial_interpolation(input_pattern, mud):
     """Compute cve using polynomial interpolation method,
-    raise an error if the mu*D value is out of the range (0.5 to 7).
+    default to brute-force computation if mu*D is
+    out of the range (0.5 to 7).
     """
     if mud > 7 or mud < 0.5:
-        raise ValueError(
+        warnings.warn(
             f"Input mu*D = {mud} is out of the acceptable range "
-            f"({min(MUD_LIST)} to {max(MUD_LIST)}) "
+            f"({np.min(MUD_LIST)} to {np.max(MUD_LIST)}) "
             f"for polynomial interpolation. "
-            f"Please rerun with a value within this range "
-            f"or specifying another method from {*CVE_METHODS, }."
+            f"Proceeding with brute-force computation. "
         )
-    coef1, coef2, coef3, coef4, coef5, coef6, coef7 = [
-        interpolation_function(mud)
-        for interpolation_function in INTERPOLATION_FUNCTIONS
-    ]
-    muls = np.array(
-        coef1 * MULS**6
-        + coef2 * MULS**5
-        + coef3 * MULS**4
-        + coef4 * MULS**3
-        + coef5 * MULS**2
-        + coef6 * MULS
-        + coef7
-    )
-    cve = 1 / muls
+        return _cve_brute_force(input_pattern, mud)
 
+    coeffs = np.array([f(mud) for f in INTERPOLATION_FUNCTIONS])
+    muls = np.polyval(coeffs, MULS)
+    cve = 1 / muls
     cve_do = DiffractionObject(
         xarray=TTH_GRID,
         yarray=cve,

tests/test_functions.py

@@ -105,49 +105,55 @@ def test_set_muls_at_angle(input_mu, expected_muls):
 
 
 @pytest.mark.parametrize(
-    "input_xtype, expected",
-    [
-        (
-            "tth",
+    "input_diffraction_data, input_cve_params",
+    [  # Test that cve diffraction object contains the expected info
+        # Note that all cve values are interpolated to 0.5
+        # cve do should contain the same input xarray, xtype,
+        # wavelength, and metadata
+        (  # C1: User did not specify method, default to fast calculation
             {
                 "xarray": np.array([90, 90.1, 90.2]),
-                "yarray": np.array([0.5, 0.5, 0.5]),
-                "xtype": "tth",
+                "yarray": np.array([2, 2, 2]),
             },
+            {"mud": 1, "xtype": "tth"},
         ),
-        (
-            "q",
+        (  # C2: User specified brute-force computation method
             {
-                "xarray": np.array([5.76998, 5.77501, 5.78004]),
-                "yarray": np.array([0.5, 0.5, 0.5]),
-                "xtype": "q",
+                "xarray": np.array([5.1, 5.2, 5.3]),
+                "yarray": np.array([2, 2, 2]),
             },
+            {"mud": 1, "method": "brute_force", "xtype": "q"},
+        ),
+        (  # C3: User specified mu*D outside the fast calculation range,
+            # default to brute-force computation
+            {
+                "xarray": np.array([5.1, 5.2, 5.3]),
+                "yarray": np.array([2, 2, 2]),
+            },
+            {"mud": 20, "xtype": "q"},
         ),
     ],
 )
-def test_compute_cve(input_xtype, expected, mocker):
-    xarray, yarray = np.array([90, 90.1, 90.2]), np.array([2, 2, 2])
+def test_compute_cve(mocker, input_diffraction_data, input_cve_params):
+    expected_xarray = input_diffraction_data["xarray"]
     expected_cve = np.array([0.5, 0.5, 0.5])
+    expected_xtype = input_cve_params["xtype"]
+    mocker.patch("diffpy.labpdfproc.functions.N_POINTS_ON_DIAMETER", 4)
     mocker.patch("numpy.interp", return_value=expected_cve)
     input_pattern = DiffractionObject(
-        xarray=xarray,
-        yarray=yarray,
-        xtype="tth",
+        xarray=input_diffraction_data["xarray"],
+        yarray=input_diffraction_data["yarray"],
+        xtype=input_cve_params["xtype"],
         wavelength=1.54,
         scat_quantity="x-ray",
         name="test",
         metadata={"thing1": 1, "thing2": "thing2"},
     )
-    actual_cve_do = compute_cve(
-        input_pattern,
-        mud=1,
-        method="polynomial_interpolation",
-        xtype=input_xtype,
-    )
+    actual_cve_do = compute_cve(input_pattern, **input_cve_params)
     expected_cve_do = DiffractionObject(
-        xarray=expected["xarray"],
-        yarray=expected["yarray"],
-        xtype=expected["xtype"],
+        xarray=expected_xarray,
+        yarray=expected_cve,
+        xtype=expected_xtype,
         wavelength=1.54,
         scat_quantity="cve",
         name="absorption correction, cve, for test",
@@ -156,32 +162,9 @@ def test_compute_cve(input_xtype, expected, mocker):
     assert actual_cve_do == expected_cve_do
 
 
-@pytest.mark.parametrize(
-    "inputs, msg",
-    [
-        (
-            {"mud": 10, "method": "polynomial_interpolation"},
-            f"mu*D = 10 is out of the acceptable range (0.5 to 7) "
-            f"for polynomial interpolation. "
-            f"Please rerun with a value within this range "
-            f"or specifying another method from {*CVE_METHODS, }.",
-        ),
-        (
-            {"mud": 1, "method": "invalid_method"},
-            f"Unknown method: invalid_method. "
-            f"Allowed methods are {*CVE_METHODS, }.",
-        ),
-        (
-            {"mud": 7, "method": "invalid_method"},
-            f"Unknown method: invalid_method. "
-            f"Allowed methods are {*CVE_METHODS, }.",
-        ),
-    ],
-)
-def test_compute_cve_bad(mocker, inputs, msg):
+def test_compute_cve_bad(mocker):
     xarray, yarray = np.array([90, 90.1, 90.2]), np.array([2, 2, 2])
     expected_cve = np.array([0.5, 0.5, 0.5])
-    mocker.patch("diffpy.labpdfproc.functions.TTH_GRID", xarray)
     mocker.patch("numpy.interp", return_value=expected_cve)
     input_pattern = DiffractionObject(
         xarray=xarray,
@@ -192,14 +175,21 @@ def test_compute_cve_bad(mocker, inputs, msg):
         name="test",
         metadata={"thing1": 1, "thing2": "thing2"},
     )
-    with pytest.raises(ValueError, match=re.escape(msg)):
-        compute_cve(input_pattern, mud=inputs["mud"], method=inputs["method"])
+    # Test that the function raises a ValueError
+    # when an invalid method is provided
+    with pytest.raises(
+        ValueError,
+        match=re.escape(
+            f"Unknown method: invalid_method. "
+            f"Allowed methods are {*CVE_METHODS, }."
+        ),
+    ):
+        compute_cve(input_pattern, mud=1, method="invalid_method")
 
 
 def test_apply_corr(mocker):
     xarray, yarray = np.array([90, 90.1, 90.2]), np.array([2, 2, 2])
     expected_cve = np.array([0.5, 0.5, 0.5])
-    mocker.patch("diffpy.labpdfproc.functions.TTH_GRID", xarray)
     mocker.patch("numpy.interp", return_value=expected_cve)
     input_pattern = DiffractionObject(
         xarray=xarray,