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,7 +37,7 @@ 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?
+    # def get_coordinate_index(self, coordinate):
     #     count = 0
     #     for i, target in enumerate(self.grid):
     #         if coordinate == target:
@@ -172,9 +182,9 @@ def get_path_length(self, grid_point, angle):
         return total_distance, primary_distance, secondary_distance
 
 
-def compute_cve(diffraction_data, mud, wavelength):
+def compute_cve(diffraction_data, mud, wavelength, brute_force=False):
     """
-    compute the cve for given diffraction data, mud and wavelength
+    compute the cve for given diffraction data, mud and wavelength, and a boolean to determine the way to compute
 
     Parameters
     ----------
@@ -185,28 +195,41 @@ def compute_cve(diffraction_data, mud, wavelength):
     wavelength float
       the wavelength of the diffraction object
 
-    Returns
-    -------
-    the diffraction object with cve curves
-
-    it is computed as follows:
+    the brute-force method is computed as follows:
     We first resample data and absorption correction to a more reasonable grid,
     then calculate corresponding cve for the given mud in the resample grid
     (since the same mu*D yields the same cve, we can assume that D/2=1, so mu=mud/2),
     and finally interpolate cve to the original grid in diffraction_data.
+
+    Returns
+    -------
+    the diffraction object with cve curves
+
     """
 
-    mu_sample_invmm = mud / 2
-    abs_correction = Gridded_circle(mu=mu_sample_invmm)
-    distances, muls = [], []
-    for angle in TTH_GRID:
-        abs_correction.set_distances_at_angle(angle)
-        abs_correction.set_muls_at_angle(angle)
-        distances.append(sum(abs_correction.distances))
-        muls.append(sum(abs_correction.muls))
-    distances = np.array(distances) / abs_correction.total_points_in_grid
-    muls = np.array(muls) / abs_correction.total_points_in_grid
-    cve = 1 / muls
+    if brute_force:
+        mu_sample_invmm = mud / 2
+        abs_correction = Gridded_circle(mu=mu_sample_invmm)
+        distances, muls = [], []
+        for angle in TTH_GRID:
+            abs_correction.set_distances_at_angle(angle)
+            abs_correction.set_muls_at_angle(angle)
+            distances.append(sum(abs_correction.distances))
+            muls.append(sum(abs_correction.muls))
+        distances = np.array(distances) / abs_correction.total_points_in_grid
+        muls = np.array(muls) / abs_correction.total_points_in_grid
+        cve = 1 / muls
+    else:
+        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
 
     orig_grid = diffraction_data.on_tth[0]
     newcve = np.interp(orig_grid, TTH_GRID, cve)

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(