feat: increase fast cal to muD=7

news/fastcalto7.rst

@@ -0,0 +1,24 @@
+**Added:**
+
+* Fast calculation supports values up to muD = 7
+
+**Changed:**
+
+* 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:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/diffpy/labpdfproc/data/coefficient_list.csv

@@ -1,5 +1,7 @@
--6.543824057313183,2.517676049711881e-10,-430.98705781061807,-1390.1242894055792,-2233.6043531698424,-2624.226082870084,-2590.22123249732
-5.020666247920394,-4.676545705328227e-10,875.3743680589598,2744.521375373735,4366.230491983095,5113.620084734874,5051.51673764446
-1.8615172482459748,3.25649449911046e-10,-663.8228448830648,-2027.2101338882571,-3194.0427904869357,-3728.3876362721985,-3684.6430045258744
--1.9688714316525813,0.9999999998992586,224.64320300800793,666.0546805817288,1038.3830530062712,1207.4338980088214,1193.1790131535502
-0.9819818977427489,1.1680695590205036e-11,-28.533774174029173,-82.17596308875525,-126.6770208129477,-146.65755463691826,-144.8506009433783
+-20619.128648244743,2.4364997877410417e-06,26505.585137008606,-337279.3213902739,-1056051.2792994028,-1815539.1160187968,-2372171.980894541,-2642993.4538858426
+56203.20048346139,-6.774323967911372e-06,-51225.01460831775,1005023.6254387972,3051609.2101441943,5201294.212075991,6773264.754466731,7538833.212217793
+-63802.106117440504,7.845612740225726e-06,33089.880661840034,-1242527.9972772636,-3668761.3028854067,-6202335.40819644,-8050708.344579743,-8951452.327576341
+38603.18842240787,-4.844608638276231e-06,-4028.752491140328,816333.2796219026,2349291.871863084,3940828.5623758254,5099123.959731602,5663734.051678175
+-13126.533425725584,1.6822282153058894e-06,-4395.4012467221355,-300757.7094990732,-845218.597424347,-1407243.8630987857,-1815244.4369415767,-2014105.8565458413
+2378.155272695758,0.9999996885549859,1912.0755691304305,58944.34388638723,162014.93660541167,267802.2088119374,344395.5526651557,381710.2334140182
+-178.70587585316136,2.4018054840276848e-08,-234.76082555771987,-4803.1698085743365,-12928.521798999534,-21220.235622246797,-27207.092578054173,-30121.285267280207

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]
+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,28 +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 6).
+    default to brute-force computation if mu*D is
+    out of the range (0.5 to 7).
     """
-    if mud > 6 or mud < 0.5:
-        raise ValueError(
-            f"mu*D is out of the acceptable range (0.5 to 6) "
+    if mud > 7 or mud < 0.5:
+        warnings.warn(
+            f"Input mu*D = {mud} is out of the acceptable range "
+            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. "
         )
-    coeff_a, coeff_b, coeff_c, coeff_d, coeff_e = [
-        interpolation_function(mud)
-        for interpolation_function in INTERPOLATION_FUNCTIONS
-    ]
-    muls = np.array(
-        coeff_a * MULS**4
-        + coeff_b * MULS**3
-        + coeff_c * MULS**2
-        + coeff_d * MULS
-        + coeff_e
-    )
-    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": 7, "method": "polynomial_interpolation"},
-            f"mu*D is out of the acceptable range (0.5 to 6) "
-            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,