Works for all parameters now

Author: Sparks29032

docs/source/morphpy.rst

@@ -86,6 +86,11 @@ get_diff: 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,53 +124,7 @@ def build_morph_inputs_container(
     return morph_inputs
 
 
-def single_morph_output(
-    morph_inputs,
-    morph_results,
-    uncertainties=None,
-    save_file=None,
-    morph_file=None,
-    xy_out=None,
-    verbose=False,
-    stdout_flag=False,
-):
-    """Helper function for printing details about a single morph.
-    Handles both printing to terminal and printing to a file.
-
-    Parameters
-    ----------
-    morph_inputs: dict
-        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
-        Name of the morphed function file. Required when printing to a
-        non-terminal file.
-    param xy_out: list
-        List of the form [x_morph_out, y_morph_out]. x_morph_out is a List of
-        r values and y_morph_out is a List of gr values.
-    verbose: bool
-        Print additional details about the morph when True (default False).
-    stdout_flag: bool
-        Print to terminal when True (default False).
-    """
-
-    print(uncertainties)
-
-    # Input and output parameters
-    morphs_in = "\n# Input morphing parameters:\n"
-    morphs_in += (
-        "\n".join(
-            f"# {key} = {morph_inputs[key]}" for key in morph_inputs.keys()
-        )
-        + "\n"
-    )
-
-    mr_copy = morph_results.copy()
+def get_terminal_morph_output(mr_copy, uncertainties):
     morphs_out = "# Optimized morphing parameters:\n"
     # Handle special inputs (numerical)
     if "squeeze" in mr_copy:
@@ -202,7 +156,7 @@ def single_morph_output(
                 )
             mr_copy = dict(morph_results_list)
 
-    # Normal inputs
+    # Get uncertainties
     if uncertainties is None:
         morphs_out += "\n".join(
             f"# {key} = {mr_copy[key]:.6f}" for key in mr_copy.keys()
@@ -218,6 +172,56 @@ def single_morph_output(
             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,
+    verbose=False,
+    stdout_flag=False,
+):
+    """Helper function for printing details about a single morph.
+    Handles both printing to terminal and printing to a file.
+
+    Parameters
+    ----------
+    morph_inputs: dict
+        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
+        Name of the morphed function file. Required when printing to a
+        non-terminal file.
+    param xy_out: list
+        List of the form [x_morph_out, y_morph_out]. x_morph_out is a List of
+        r values and y_morph_out is a List of gr values.
+    verbose: bool
+        Print additional details about the morph when True (default False).
+    stdout_flag: bool
+        Print to terminal when True (default False).
+    """
+
+    # Input and output parameters
+    morphs_in = "\n# Input morphing parameters:\n"
+    morphs_in += (
+        "\n".join(
+            f"# {key} = {morph_inputs[key]}" for key in morph_inputs.keys()
+        )
+        + "\n"
+    )
+
+    mr_copy = morph_results.copy()
+    morphs_out, func_dicts = get_terminal_morph_output(mr_copy, uncertainties)
+
     # Handle special inputs (functional add)
     for func in func_dicts.keys():
         if func_dicts[func][0] is not None:
@@ -356,6 +360,7 @@ def multiple_morph_output(
     morph_inputs,
     morph_results,
     target_files,
+    uncertainties_dict=None,
     field=None,
     field_list=None,
     save_directory=None,
@@ -376,6 +381,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
@@ -415,11 +422,16 @@ 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]
-            )
+            # 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
 
@@ -747,6 +748,25 @@ def single_morph(
                 refiner.refine(*rptemp)
             # Adjust all parameters
             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
@@ -886,10 +906,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):
@@ -999,6 +1021,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
@@ -1008,13 +1031,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():
@@ -1036,6 +1057,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,
@@ -1193,6 +1215,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
@@ -1202,13 +1225,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():
@@ -1230,6 +1251,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,

src/diffpy/morph/refine.py

@@ -204,12 +204,22 @@ def refine(self, *args, estimate_uncertainty=False, **kw):
         self._update_chain(vals)
 
         if estimate_uncertainty:
+            par_names = list(self.pars)
             if hess_inv_sol is None:
-                warnings.warn(
-                    "Warning: Could not estimate "
-                    "uncertainty as estimated Hessian is singular.",
-                    UserWarning,
-                )
+                if len(par_names) == 1:
+                    warnings.warn(
+                        "Warning: Could not estimate "
+                        f"uncertainty of {par_names[0].lower()} parameter "
+                        "as estimated Hessian is singular.",
+                        UserWarning,
+                    )
+                else:
+                    warnings.warn(
+                        "Warning: Certain parameters may be "
+                        "missing uncertainties as estimated Hessian "
+                        "is singular.",
+                        UserWarning,
+                    )
                 return None
             dof = len(fvec) - len(sol)
             if dof <= 0:
@@ -220,7 +230,23 @@ def refine(self, *args, estimate_uncertainty=False, **kw):
                     UserWarning,
                 )
                 return None
-            par_names = list(self.pars)
+            # Handle squeeze morph params
+            if "squeeze" in par_names and "squeeze" in config.keys():
+                squeeze_par_names = [
+                    f"squeeze a{i}"
+                    for i in range(len(config["squeeze"].keys()))
+                ]
+                squeeze_idx = par_names.index("squeeze")
+                par_names[squeeze_idx : squeeze_idx + 1] = squeeze_par_names
+            # Handle funcx, funcy, funcxy morph params
+            func_types = ["funcx", "funcy", "funcxy"]
+            for func in func_types:
+                if func in par_names and func in config.keys():
+                    func_par_names = [
+                        f"{func} {fparam}" for fparam in config[func].keys()
+                    ]
+                    func_idx = par_names.index(func)
+                    par_names[func_idx : func_idx + 1] = func_par_names
             cov_sol = hess_inv_sol * dot(fvec, fvec) / dof
             unc = sqrt(diag(cov_sol))