Add --smear-function option and Lorentzian smearing (#309)

* Add Lorentzian smearing

* More verbose test descriptors

* Fix headers for tests

Author: Sparks29032

docs/source/morphpy.rst

@@ -139,17 +139,28 @@ squeeze: list of float
     When this parameter is given, hshift is disabled.
     When n>1, stretch is disabled.
 smear: float
-    Smear the peaks with a Gaussian of width smear.
+    Smear the peaks with a function of width SMEAR.
 
+    The smearing function is chosen by the --smear-function option
+    (default is Gaussian if that option is not enabled).
     This is done by convolving the function with a Gaussian
     with standard deviation smear. If both smear and
     smear_pdf are used, only smear_pdf will be
     applied.
 smear_pdf: float
     Convert PDF to RDF. Then, smear peaks with a Gaussian
-    of width smear_pdf. Convert back to PDF. If both smear and
+    of width smear_pdf.
+    The smearing function is chosen by the --smear-function option
+    (default is Gaussian if that option is not enabled).
+    Convert back to PDF. If both smear and
     smear_pdf are used, only smear_pdf will be
     applied.
+smear_func: str
+    Choose the function for the smear morph.
+    Only used if --smear or --smear-pdf is enabled.
+    Available options: Gaussian (default) and Lorentzian.
+    If Gaussian, the SMEAR parameter is (+/-) the standard deviation.
+    If Lorentzian, the SMEAR parameter is (+/-) the half width half maximum.
 slope: float
     Slope of the baseline used in converting from PDF to RDF.
 

news/general-smear.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* New option `--smear-function` allows user to select either gaussian or lorentzian smearing. The default is still gaussian, so no workflow changes are needed for those who are using the gaussian smearing.
+
+**Changed:**
+
+* The choice of smear function is now saved to the file header when `--verbose` is enabled.
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/diffpy/morph/morph_io.py

@@ -39,6 +39,7 @@ def build_morph_inputs_container(
     stretch,
     smear_pdf,
     smear,
+    smear_func,
     hshift,
     vshift,
     squeeze,
@@ -104,11 +105,27 @@ def build_morph_inputs_container(
         smear_in = smear
     else:
         smear_in = smear_pdf
-    morph_inputs = {
-        "scale": scale_in,
-        "stretch": stretch_in,
-        "smear": smear_in,
-    }
+    if smear_func is not None:
+        smear_func = smear_func.lower()
+    if smear_func is None or smear_func == "gaussian":
+        smear_func_in = "gaussian"
+    elif smear_func == "lorentzian":
+        smear_func_in = "lorentzian"
+    else:
+        smear_func_in = "unknown"
+    if smear_in is None:
+        morph_inputs = {
+            "scale": scale_in,
+            "stretch": stretch_in,
+            "smear": smear_in,
+        }
+    else:
+        morph_inputs = {
+            "scale": scale_in,
+            "stretch": stretch_in,
+            "smear": smear_in,
+            "smear-function": smear_func_in,
+        }
 
     if squeeze_poly_deg < 0:
         hshift_in = hshift
@@ -126,6 +143,11 @@ def build_morph_inputs_container(
 
 def get_terminal_morph_output(mr_copy, uncertainties):
     morphs_out = "# Optimized morphing parameters:\n"
+
+    # Handle special inputs (strings)
+    if "smear_func" in mr_copy:
+        mr_copy.pop("smear_func")
+
     # Handle special inputs (numerical)
     if "squeeze" in mr_copy:
         sq_dict = mr_copy.pop("squeeze")

src/diffpy/morph/morphapp.py

@@ -245,9 +245,11 @@ def morph_error(self, msg, error):
         type="float",
         metavar="SMEAR",
         help=(
-            "Smear the peaks with a Gaussian of width SMEAR. "
+            "Smear the peaks with a function of width SMEAR. "
+            "The smearing function is chosen by the --smear-function option "
+            "(default is Gaussian if that option is not enabled). "
             "This is done by convolving the function with a "
-            "Gaussian with standard deviation SMEAR. "
+            "function with width SMEAR. "
             "If both --smear and --smear-pdf are enabled, "
             "only --smear-pdf will be applied."
         ),
@@ -258,12 +260,29 @@ def morph_error(self, msg, error):
         metavar="SMEAR",
         help=(
             "Convert PDF to RDF. "
-            "Then, smear peaks with a Gaussian of width SMEAR. "
+            "Then, smear peaks with a function of width SMEAR. "
+            "The smearing function is chosen by the --smear-function option "
+            "(default is Gaussian if that option is not enabled). "
             "Convert back to PDF. "
             "If both --smear and --smear-pdf are enabled, "
             "only --smear-pdf will be applied."
         ),
     )
+    group.add_option(
+        "--smear-function",
+        "--smear-func",
+        dest="smear_func",
+        metavar="SMEARFUNCTION",
+        help=(
+            "Choose the function for the smear morph. "
+            "Only used if --smear or --smear-pdf is enabled. "
+            "Available options: Gaussian (default) and Lorentzian. "
+            "If Gaussian, the SMEAR parameter is the standard "
+            "deviation. "
+            "If Lorentzian, the SMEAR parameter is the half width "
+            "half maximum."
+        ),
+    )
     group.add_option(
         "--slope",
         type="float",
@@ -665,13 +684,17 @@ def single_morph(
         config["stretch"] = stretch_in
         refpars.append("stretch")
     # Smear
+    smear_func = "gaussian"
+    if opts.smear_func is not None:
+        smear_func = opts.smear_func.lower()
     if opts.smear_pdf is not None:
         smear_in = opts.smear_pdf
         chain.append(helpers.TransformXtalPDFtoRDF())
         chain.append(morphs.MorphSmear())
         chain.append(helpers.TransformXtalRDFtoPDF())
         refpars.append("smear")
         config["smear"] = smear_in
+        config["smear_func"] = smear_func
         # Set baselineslope if not given
         config["baselineslope"] = opts.baselineslope
         if opts.baselineslope is None:
@@ -682,6 +705,7 @@ def single_morph(
         chain.append(morphs.MorphSmear())
         refpars.append("smear")
         config["smear"] = smear_in
+        config["smear_func"] = smear_func
     # Shift
     # Only enable hshift is squeeze is not enabled
     shift_morph = None
@@ -825,6 +849,7 @@ def single_morph(
         opts.stretch,
         opts.smear_pdf,
         opts.smear,
+        opts.smear_func,
         opts.hshift,
         opts.vshift,
         opts.squeeze,
@@ -1060,6 +1085,7 @@ def multiple_targets(parser, opts, pargs, stdout_flag=True, python_wrap=False):
         opts.stretch,
         opts.smear_pdf,
         opts.smear,
+        opts.smear_func,
         opts.hshift,
         opts.vshift,
         opts.squeeze,
@@ -1254,6 +1280,7 @@ def multiple_morphs(parser, opts, pargs, stdout_flag=True, python_wrap=False):
         opts.stretch,
         opts.smear_pdf,
         opts.smear,
+        opts.smear_func,
         opts.hshift,
         opts.vshift,
         opts.squeeze,

src/diffpy/morph/morphs/morphsmear.py

@@ -37,7 +37,7 @@ class MorphSmear(Morph):
     yinlabel = LABEL_RR
     xoutlabel = LABEL_RA
     youtlabel = LABEL_RR
-    parnames = ["smear"]
+    parnames = ["smear", "smear_func"]
 
     def morph(self, x_morph, y_morph, x_target, y_target):
         """Resample arrays onto specified grid."""
@@ -51,10 +51,23 @@ def morph(self, x_morph, y_morph, x_target, y_target):
         r = self.x_morph_in
         rr = self.y_morph_in
         r0 = r[len(r) // 2]
-        gaussian = numpy.exp(-0.5 * ((r - r0) / self.smear) ** 2)
+
+        function = None
+        # FIXME: Add a condition here that checks if smear_func is a function.
+        # FIXME: If smear_func is a function, set function = self.smear_func.
+        if self.smear_func is None or self.smear_func == "gaussian":
+            function = numpy.exp(-0.5 * ((r - r0) / self.smear) ** 2)
+        elif self.smear_func.lower() == "lorentzian":
+            function = numpy.abs(self.smear) / ((r - r0) ** 2 + self.smear**2)
+        else:
+            raise ValueError(
+                "Could not process smearing function. "
+                "Smearing function must be either "
+                "'gaussian' or 'lorentzian'."
+            )
 
         # Get the full convolution
-        c = numpy.convolve(rr, gaussian, mode="full")
+        c = numpy.convolve(rr, function, mode="full")
         # Find the centroid of the RDF, we don't want this to change from the
         # convolution.
         x1 = numpy.arange(len(rr), dtype=float)
@@ -69,7 +82,7 @@ def morph(self, x_morph, y_morph, x_target, y_target):
         rrbroad = numpy.interp(x1, xc, c)
 
         # Normalize so that the integrated magnitude of the RDF doesn't change.
-        rrbroad /= sum(gaussian)
+        rrbroad /= sum(function)
 
         self.y_morph_out = rrbroad
 

tests/test_morph_func.py

@@ -93,6 +93,7 @@ def test_smear_with_morph_func():
     x_morph = x_target.copy()
     y_morph = np.exp(-0.5 * ((x_morph - r0) / sigma0) ** 2)
     cfg = morph_default_config(smear=0.1, scale=1.1, stretch=0.1)  # off init
+    cfg.update({"smear_func": None})
     morph_rv = morph(x_morph, y_morph, x_target, y_target, **cfg)
     morphed_cfg = morph_rv["morphed_config"]
     # verified they are morphable

tests/test_morphapp.py

@@ -243,6 +243,21 @@ def test_parser_errors(self, capsys, setup_parser):
             single_morph(self.parser, opts, pargs, stdout_flag=False)
         assert "a could not be converted to float." in str(excinfo.value)
 
+        # Non-sensical smearing function
+        opts, pargs = self.parser.parse_args(
+            [
+                f"{nickel_PDF}",
+                f"{nickel_PDF}",
+                "--smear",
+                "0",
+                "--smear-function",
+                "not-a-smearing-function",
+            ]
+        )
+        with pytest.raises(ValueError) as excinfo:
+            single_morph(self.parser, opts, pargs, stdout_flag=False)
+        assert "Could not process smearing function." in str(excinfo.value)
+
     def test_morphsequence(self, setup_morphsequence):
         # Parse arguments sorting by field
         opts, pargs = self.parser.parse_args(