add fast_cve to labpdfprocapp.py and some other fixes

pyproject.toml

@@ -42,6 +42,9 @@ template = "{tag}"
 dev_template = "{tag}"
 dirty_template = "{tag}"
 
+[project.scripts]
+labpdfproc = "diffpy.labpdfproc.labpdfprocapp:main"
+
 [tool.setuptools.packages.find]
 where = ["src"]  # list of folders that contain the packages (["."] by default)
 include = ["*"]  # package names should match these glob patterns (["*"] by default)

src/diffpy/labpdfproc/functions.py

@@ -1,12 +1,22 @@
 import math
+import os
 
 import numpy as np
+import pandas as pd
+from scipy.interpolate import interp1d
 
 from diffpy.utils.scattering_objects.diffraction_objects import Diffraction_object
 
 RADIUS_MM = 1
 N_POINTS_ON_DIAMETER = 300
-TTH_GRID = np.arange(1, 141, 1)
+TTH_GRID = np.arange(1, 180.1, 0.1)
+
+# pre-computed datasets for fast calculation
+MUD_LIST = [0.5, 1, 2, 3, 4, 5, 6]
+CWD = os.path.dirname(os.path.abspath(__file__))
+MULS = np.loadtxt(CWD + "/data/inverse_cve.xy")
+COEFFICIENT_LIST = np.array(pd.read_csv(CWD + "/data/coefficient_list.csv", header=None))
+INTERPOLATION_FUNCTIONS = [interp1d(MUD_LIST, coefficients, kind="quadratic") for coefficients in COEFFICIENT_LIST]
 
 
 class Gridded_circle:
@@ -27,16 +37,6 @@ def _get_grid_points(self):
         self.grid = {(x, y) for x in xs for y in ys if x**2 + y**2 <= self.radius**2}
         self.total_points_in_grid = len(self.grid)
 
-    # def get_coordinate_index(self, coordinate):  # I think we probably dont need this function?
-    #     count = 0
-    #     for i, target in enumerate(self.grid):
-    #         if coordinate == target:
-    #             return i
-    #         else:
-    #             count += 1
-    #     if count >= len(self.grid):
-    #         raise IndexError(f"WARNING: no coordinate {coordinate} found in coordinates list")
-
     def set_distances_at_angle(self, angle):
         """
         given an angle, set the distances from the grid points to the entry and exit coordinates
@@ -172,22 +172,9 @@ def get_path_length(self, grid_point, angle):
         return total_distance, primary_distance, secondary_distance
 
 
-def compute_cve(diffraction_data, mud, wavelength):
+def _cve_brute_force(diffraction_data, mud, wavelength):
     """
-    compute the cve for given diffraction data, mud and wavelength
-
-    Parameters
-    ----------
-    diffraction_data Diffraction_object
-      the diffraction pattern
-    mud float
-      the mu*D of the diffraction object, where D is the diameter of the circle
-    wavelength float
-      the wavelength of the diffraction object
-
-    Returns
-    -------
-    the diffraction object with cve curves
+    compute cve using brute-force method
 
     it is computed as follows:
     We first resample data and absorption correction to a more reasonable grid,
@@ -207,7 +194,57 @@ def compute_cve(diffraction_data, mud, wavelength):
     distances = np.array(distances) / abs_correction.total_points_in_grid
     muls = np.array(muls) / abs_correction.total_points_in_grid
     cve = 1 / muls
+    return cve
+
+
+def _cve_interp_polynomial(diffraction_data, mud, wavelength):
+    """
+    compute cve using polynomial interpolation method, raise an error if mu*D is out of the range (0.5 to 6)
+    """
+
+    if mud > 6 or mud < 0.5:
+        raise ValueError(
+            "mu*D is out of the acceptable range (0.5 to 6) for fast calculation. "
+            "Please rerun with a value within this range or use -b to enable brute-force calculation. "
+        )
+    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
+
 
+def _cve_method(diffraction_data, mud, wavelength, brute_force=False):
+    """
+    selects the appropriate CVE calculation method
+    """
+    if brute_force:
+        return _cve_brute_force(diffraction_data, mud, wavelength)
+    else:
+        return _cve_interp_polynomial(diffraction_data, mud, wavelength)
+
+
+def compute_cve(diffraction_data, mud, wavelength, brute_force=False):
+    """
+    compute the cve for given diffraction data, mud, and wavelength, using the selected method
+
+    Parameters
+    ----------
+    diffraction_data Diffraction_object
+      the diffraction pattern
+    mud float
+      the mu*D of the diffraction object, where D is the diameter of the circle
+    wavelength float
+      the wavelength of the diffraction object
+
+    Returns
+    -------
+    the diffraction object with cve curves
+
+    """
+
+    cve = _cve_method(diffraction_data, mud, wavelength, brute_force)
     orig_grid = diffraction_data.on_tth[0]
     newcve = np.interp(orig_grid, TTH_GRID, cve)
     abdo = Diffraction_object(wavelength=wavelength)
@@ -220,7 +257,6 @@ def compute_cve(diffraction_data, mud, wavelength):
         wavelength=diffraction_data.wavelength,
         scat_quantity="cve",
     )
-
     return abdo
 
 

src/diffpy/labpdfproc/labpdfprocapp.py

@@ -64,14 +64,20 @@ def get_args(override_cli_inputs=None):
         action="store_true",
         help="The absorption correction will be output to a file if this "
         "flag is set. Default is that it is not output.",
-        default="tth",
     )
     p.add_argument(
         "-f",
         "--force-overwrite",
         action="store_true",
         help="Outputs will not overwrite existing file unless --force is specified.",
     )
+    p.add_argument(
+        "-b",
+        "--brute-force",
+        action="store_true",
+        help="The absorption correction will be computed using brute-force calculation "
+        "if this flag is set. Default is using fast calculation. ",
+    )
     p.add_argument(
         "-u",
         "--user-metadata",
@@ -134,7 +140,7 @@ def main():
             metadata=load_metadata(args, filepath),
         )
 
-        absorption_correction = compute_cve(input_pattern, args.mud, args.wavelength)
+        absorption_correction = compute_cve(input_pattern, args.mud, args.wavelength, args.brute_force)
         corrected_data = apply_corr(input_pattern, absorption_correction)
         corrected_data.name = f"Absorption corrected input_data: {input_pattern.name}"
         corrected_data.dump(f"{outfile}", xtype="tth")

src/diffpy/labpdfproc/tests/test_tools.py

@@ -312,6 +312,7 @@ def test_load_metadata(mocker, user_filesystem):
             "wavelength": 0.71,
             "output_directory": str(Path.cwd().resolve()),
             "xtype": "tth",
+            "brute_force": False,
             "key": "value",
             "username": "cli_username",
             "email": "[email protected]",

src/diffpy/labpdfproc/tools.py

@@ -17,15 +17,15 @@ def set_output_directory(args):
     args argparse.Namespace
         the arguments from the parser
 
-    Returns
-    -------
-    pathlib.PosixPath that contains the full path of the output directory
-
     it is determined as follows:
     If user provides an output directory, use it.
     Otherwise, we set it to the current directory if nothing is provided.
     We then create the directory if it does not exist.
 
+    Returns
+    -------
+    pathlib.PosixPath that contains the full path of the output directory
+
     """
     output_dir = Path(args.output_directory).resolve() if args.output_directory else Path.cwd().resolve()
     output_dir.mkdir(parents=True, exist_ok=True)
@@ -110,13 +110,13 @@ def set_wavelength(args):
     args argparse.Namespace
         the arguments from the parser
 
+    we raise an ValueError if the input wavelength is non-positive
+    or if the input anode_type is not one of the known sources
+
     Returns
     -------
     args argparse.Namespace
 
-    we raise an ValueError if the input wavelength is non-positive
-    or if the input anode_type is not one of the known sources
-
     """
     if args.wavelength is not None and args.wavelength <= 0:
         raise ValueError(