Add uncertainties for main parameters. (#299)

* Add uncertainties for main parameters.

* Works for all parameters now

* Fix other tests

* loosen tolerance to 1e-4 from 1e-5

* Lower tolerance more

* Remove comment

* [pre-commit.ci] auto fixes from pre-commit hooks

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: Sparks29032

docs/source/morphpy.rst

@@ -92,6 +92,11 @@ original_grid: bool
 verbose: bool
     Print additional header details to saved files. These include details about the morph
     inputs and outputs.
+uncertainty: bool
+    Estimate uncertainties for each refined morphing parameter. This is done by estimating
+    the Hessian of the fit. Caution should be taken as this is not the true uncertainty
+    of the fit, and the user should make their own judgement about what measure of uncertainty
+    to use for their particular use case.
 xmin: float
     Minimum x-value (abscissa) to use for function comparisons.
 xmax: float

docs/source/quickstart.rst

@@ -202,7 +202,13 @@ Basic diffpy.morph Workflow
 
        The optimal fit after applying the scale, smear, and stretch morphs.
 
-    9. Now, try it on your own! If you have personally collected or
+    9. We are also able to estimate the uncertainties of each of the fitted parameters.
+       This is done by using the ``uncertainty`` parameter.
+       Below we have replaced the ``-a`` from the previous step with a ``-u`` to obtain uncertainty estimates ::
+
+          diffpy.morph --scale=0.8 --smear-pdf=-0.08 --stretch=0.005 --xmin=1.5 --xmax=30 -u darkSub_rh20_C_01.gr darkSub_rh20_C_44.gr
+
+    10. Now, try it on your own! If you have personally collected or
        otherwise readily available PDF data, try this process to see if
        you can morph your PDFs to one another. Many of the parameters
        provided in this tutorial are unique to it, so be cautious about

news/uncertainty.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* New option "uncertainty" to allow user to estimate uncertainty of the refined morphing parameters.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/diffpy/morph/morph_io.py

@@ -124,9 +124,61 @@ def build_morph_inputs_container(
     return morph_inputs
 
 
+def get_terminal_morph_output(mr_copy, uncertainties):
+    morphs_out = "# Optimized morphing parameters:\n"
+    # Handle special inputs (numerical)
+    if "squeeze" in mr_copy:
+        sq_dict = mr_copy.pop("squeeze")
+        rw_pos = list(mr_copy.keys()).index("Rw")
+        morph_results_list = list(mr_copy.items())
+        for idx, _ in enumerate(sq_dict):
+            morph_results_list.insert(
+                rw_pos + idx, (f"squeeze a{idx}", sq_dict[f"a{idx}"])
+            )
+        mr_copy = dict(morph_results_list)
+
+    # Handle special inputs (functional remove)
+    func_dicts = {
+        "funcxy": [None, None],
+        "funcx": [None, None],
+        "funcy": [None, None],
+    }
+    for func in func_dicts.keys():
+        if f"{func}_function" in mr_copy:
+            func_dicts[func][0] = mr_copy.pop(f"{func}_function")
+        if func in mr_copy:
+            func_dicts[func][1] = mr_copy.pop(func)
+            rw_pos = list(mr_copy.keys()).index("Rw")
+            morph_results_list = list(mr_copy.items())
+            for idx, key in enumerate(func_dicts[func][1]):
+                morph_results_list.insert(
+                    rw_pos + idx, (f"{func} {key}", func_dicts[func][1][key])
+                )
+            mr_copy = dict(morph_results_list)
+
+    # Get uncertainties
+    if uncertainties is None:
+        morphs_out += "\n".join(
+            f"# {key} = {mr_copy[key]:.6f}" for key in mr_copy.keys()
+        )
+    else:
+        morphs_out += "\n".join(
+            f"# {key} = {mr_copy[key]:.6f}"
+            + (
+                f" +/- {uncertainties[key]:.6f}"
+                if key in uncertainties
+                else ""
+            )
+            for key in mr_copy
+        )
+
+    return morphs_out, func_dicts
+
+
 def single_morph_output(
     morph_inputs,
     morph_results,
+    uncertainties=None,
     save_file=None,
     morph_file=None,
     xy_out=None,
@@ -142,6 +194,8 @@ def single_morph_output(
         Parameters given by the user.
     morph_results: dict
         Resulting data after morphing.
+    uncertainties: dict
+        Uncertainties of all morphed parameters.
     save_file
         Name of file to print to. If None (default) print to terminal.
     morph_file
@@ -165,41 +219,7 @@ def single_morph_output(
     )
 
     mr_copy = morph_results.copy()
-    morphs_out = "# Optimized morphing parameters:\n"
-    # Handle special inputs (numerical)
-    if "squeeze" in mr_copy:
-        sq_dict = mr_copy.pop("squeeze")
-        rw_pos = list(mr_copy.keys()).index("Rw")
-        morph_results_list = list(mr_copy.items())
-        for idx, _ in enumerate(sq_dict):
-            morph_results_list.insert(
-                rw_pos + idx, (f"squeeze a{idx}", sq_dict[f"a{idx}"])
-            )
-        mr_copy = dict(morph_results_list)
-
-    # Handle special inputs (functional remove)
-    func_dicts = {
-        "funcxy": [None, None],
-        "funcx": [None, None],
-        "funcy": [None, None],
-    }
-    for func in func_dicts.keys():
-        if f"{func}_function" in mr_copy:
-            func_dicts[func][0] = mr_copy.pop(f"{func}_function")
-        if func in mr_copy:
-            func_dicts[func][1] = mr_copy.pop(func)
-            rw_pos = list(mr_copy.keys()).index("Rw")
-            morph_results_list = list(mr_copy.items())
-            for idx, key in enumerate(func_dicts[func][1]):
-                morph_results_list.insert(
-                    rw_pos + idx, (f"{func} {key}", func_dicts[func][1][key])
-                )
-            mr_copy = dict(morph_results_list)
-
-    # Normal inputs
-    morphs_out += "\n".join(
-        f"# {key} = {mr_copy[key]:.6f}" for key in mr_copy.keys()
-    )
+    morphs_out, func_dicts = get_terminal_morph_output(mr_copy, uncertainties)
 
     # Handle special inputs (functional add)
     for func in func_dicts.keys():
@@ -338,6 +358,7 @@ def multiple_morph_output(
     morph_inputs,
     morph_results,
     target_files,
+    uncertainties_dict=None,
     field=None,
     field_list=None,
     save_directory=None,
@@ -358,6 +379,8 @@ def multiple_morph_output(
         Resulting data after morphing.
     target_files: list
         Files that acted as targets to morphs.
+    uncertainties_dict: dict
+        Dictionary of uncertainties for each morph.
     save_directory
         Name of directory to save morphs in.
     field
@@ -396,11 +419,11 @@ def multiple_morph_output(
                 output = f"\n# Target: {target}\n"
             else:
                 output = f"\n# Morph: {target}\n"
-            output += "# Optimized morphing parameters:\n"
-            output += "\n".join(
-                f"# {param} = {morph_results[target][param]:.6f}"
-                for param in morph_results[target]
-            )
+
+            mr_copy = morph_results[target].copy()
+            uncertainties = uncertainties_dict[target]
+            output_body, _ = get_terminal_morph_output(mr_copy, uncertainties)
+            output += output_body
             verbose_outputs += f"{output}\n"
 
     # Get items we want to put in table

src/diffpy/morph/morphapp.py

@@ -15,6 +15,7 @@
 
 from __future__ import print_function
 
+import copy
 import sys
 from pathlib import Path
 
@@ -122,6 +123,21 @@ def morph_error(self, msg, error):
         action="store_true",
         help="Print additional header details to saved files.",
     )
+    parser.add_option(
+        "-u",
+        "--uncertainty",
+        "--estimate-uncertainty",
+        dest="estimate_uncertainty",
+        action="store_true",
+        help=(
+            "Estimate uncertainties for each refined morphing parameter. "
+            "This is done by estimating the Hessian of the fit. "
+            "Caution should be taken as this is not the true uncertainty "
+            "of the fit, and the user should make their own judgement "
+            "about what measure of uncertainty to use for their particular "
+            "use case."
+        ),
+    )
     parser.add_option(
         "--xmin",
         type="float",
@@ -734,6 +750,7 @@ def single_morph(
         refiner.residual = refiner._pearson
     if opts.addpearson:
         refiner.residual = refiner._add_pearson
+    unc = None
     if opts.refine and refpars:
         try:
             # This works better when we adjust scale and smear first.
@@ -743,17 +760,39 @@ def single_morph(
                     rptemp.append("scale")
                 refiner.refine(*rptemp)
             # Adjust all parameters
-            refiner.refine(*refpars)
+            unc = refiner.refine(*refpars, estimate_uncertainty=True)
+            # If one parameter is causing trouble with uncertainty estimate
+            # compute all uncertainties individually
+            if unc is None:
+                unc = {}
+                for param in refpars:
+                    refiner_single_param = type(refiner)(
+                        refiner.chain,
+                        refiner.x_morph,
+                        refiner.y_morph,
+                        refiner.x_target,
+                        refiner.y_target,
+                        tolerance=refiner.tolerance,
+                    )
+                    refiner_single_param.chain.config = copy.deepcopy(config)
+                    unc_param = refiner_single_param.refine(
+                        *[param], estimate_uncertainty=True
+                    )
+                    if unc_param is not None:
+                        unc.update(unc_param)
         except ValueError as e:
             parser.morph_error(str(e), ValueError)
     # Smear is not being refined, but baselineslope needs to refined to apply
     # smear
     # Note that baselineslope is only added to the refine list if smear is
     # applied
     elif "baselineslope" in refpars:
+        # Note, you cannot estimate uncertainty here as the baselineslope
+        # does not change the fit
         try:
             refiner.refine(
-                "baselineslope", baselineslope=config["baselineslope"]
+                "baselineslope",
+                baselineslope=config["baselineslope"],
             )
         except ValueError as e:
             parser.morph_error(str(e), ValueError)
@@ -839,6 +878,7 @@ def single_morph(
         io.single_morph_output(
             morph_inputs,
             morph_results,
+            uncertainties=None if opts.estimate_uncertainty is None else unc,
             save_file=opts.slocation,
             morph_file=pargs[0],
             xy_out=xy_save,
@@ -880,10 +920,12 @@ def single_morph(
     # Return different things depending on whether it is python interfaced
     if python_wrap:
         morph_info = morph_results
+        if opts.estimate_uncertainty is not None and unc is not None:
+            morph_info.update({"Uncertainties": unc})
         morph_table = numpy.array(xy_save).T
         return morph_info, morph_table
     else:
-        return morph_results
+        return morph_results, unc
 
 
 def multiple_targets(parser, opts, pargs, stdout_flag=True, python_wrap=False):
@@ -993,6 +1035,7 @@ def multiple_targets(parser, opts, pargs, stdout_flag=True, python_wrap=False):
 
     # Morph morph_file against all other files in target_directory
     morph_results = {}
+    uncs = {}
     for target_file in target_list:
         if target_file.is_file:
             # Set the save file destination to be a file within the SLOC
@@ -1002,13 +1045,11 @@ def multiple_targets(parser, opts, pargs, stdout_flag=True, python_wrap=False):
                 opts.slocation = Path(save_morphs_here).joinpath(save_as)
             # Perform a morph of morph_file against target_file
             pargs = [morph_file, target_file]
-            morph_results.update(
-                {
-                    target_file.name: single_morph(
-                        parser, opts, pargs, stdout_flag=False
-                    ),
-                }
+            morph_result, unc = single_morph(
+                parser, opts, pargs, stdout_flag=False
             )
+            morph_results.update({target_file.name: morph_result})
+            uncs.update({target_file.name: unc})
 
     target_file_names = []
     for key in morph_results.keys():
@@ -1030,6 +1071,7 @@ def multiple_targets(parser, opts, pargs, stdout_flag=True, python_wrap=False):
             morph_inputs,
             morph_results,
             target_file_names,
+            uncertainties_dict=uncs,
             save_directory=save_directory,
             morph_file=morph_file,
             target_directory=target_directory,
@@ -1187,6 +1229,7 @@ def multiple_morphs(parser, opts, pargs, stdout_flag=True, python_wrap=False):
 
     # Morph morph_file against all other files in target_directory
     morph_results = {}
+    uncs = {}
     for morph_file in morph_list:
         if morph_file.is_file:
             # Set the save file destination to be a file within the SLOC
@@ -1196,13 +1239,11 @@ def multiple_morphs(parser, opts, pargs, stdout_flag=True, python_wrap=False):
                 opts.slocation = Path(save_morphs_here).joinpath(save_as)
             # Perform a morph of morph_file against target_file
             pargs = [morph_file, target_file]
-            morph_results.update(
-                {
-                    morph_file.name: single_morph(
-                        parser, opts, pargs, stdout_flag=False
-                    ),
-                }
+            morph_result, unc = single_morph(
+                parser, opts, pargs, stdout_flag=False
             )
+            morph_results.update({morph_file.name: morph_result})
+            uncs.update({morph_file.name: unc})
 
     morph_file_names = []
     for key in morph_results.keys():
@@ -1224,6 +1265,7 @@ def multiple_morphs(parser, opts, pargs, stdout_flag=True, python_wrap=False):
             morph_inputs,
             morph_results,
             morph_file_names,
+            uncertainties_dict=uncs,
             save_directory=save_directory,
             morph_file=target_file,
             target_directory=morph_directory,