Add uncertainties for main parameters.

Author: Sparks29032

src/diffpy/morph/morph_io.py

@@ -127,6 +127,7 @@ def build_morph_inputs_container(
 def single_morph_output(
     morph_inputs,
     morph_results,
+    uncertainties=None,
     save_file=None,
     morph_file=None,
     xy_out=None,
@@ -142,6 +143,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
@@ -156,6 +159,8 @@ def single_morph_output(
         Print to terminal when True (default False).
     """
 
+    print(uncertainties)
+
     # Input and output parameters
     morphs_in = "\n# Input morphing parameters:\n"
     morphs_in += (
@@ -198,9 +203,20 @@ def single_morph_output(
             mr_copy = dict(morph_results_list)
 
     # Normal inputs
-    morphs_out += "\n".join(
-        f"# {key} = {mr_copy[key]:.6f}" for key in mr_copy.keys()
-    )
+    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
+        )
 
     # Handle special inputs (functional add)
     for func in func_dicts.keys():

src/diffpy/morph/morphapp.py

@@ -109,6 +109,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",
@@ -721,6 +736,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.
@@ -730,17 +746,20 @@ def single_morph(
                     rptemp.append("scale")
                 refiner.refine(*rptemp)
             # Adjust all parameters
-            refiner.refine(*refpars)
+            unc = refiner.refine(*refpars, estimate_uncertainty=True)
         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)
@@ -825,6 +844,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,

src/diffpy/morph/refine.py

@@ -15,7 +15,9 @@
 """refine -- Refine a morph or morph chain
 """
 
-from numpy import array, concatenate, dot, exp, ones_like
+import warnings
+
+from numpy import array, concatenate, diag, dot, exp, ones_like, sqrt
 from scipy.optimize import leastsq
 from scipy.stats import pearsonr
 
@@ -134,7 +136,7 @@ def _add_pearson(self, pvals):
         res = concatenate([res1, res2])
         return res
 
-    def refine(self, *args, **kw):
+    def refine(self, *args, estimate_uncertainty=False, **kw):
         """Refine the chain.
 
         Additional arguments are used to specify which parameters are to be
@@ -145,7 +147,9 @@ def refine(self, *args, **kw):
 
         This uses the leastsq algorithm from scipy.optimize.
 
-        This returns the final scalar residual value.
+        If estimate_uncertainty is True, return an estimated uncertainty
+        for each parameter.
+        Otherwise, return the final scalar residual value (used for testing).
         The parameters from the fit can be retrieved from the config
         dictionary of the morph or morph chain.
 
@@ -180,7 +184,7 @@ def refine(self, *args, **kw):
                 initial.append(val)
                 self.flat_to_grouped[len(initial) - 1] = (p, None)
 
-        sol, cov_sol, infodict, emesg, ier = leastsq(
+        sol, hess_inv_sol, infodict, emesg, ier = leastsq(
             self.residual,
             array(initial),
             full_output=True,
@@ -199,7 +203,30 @@ def refine(self, *args, **kw):
             vals = [vals]
         self._update_chain(vals)
 
-        return dot(fvec, fvec)
+        if estimate_uncertainty:
+            if hess_inv_sol is None:
+                warnings.warn(
+                    "Warning: Could not estimate "
+                    "uncertainty as estimated Hessian is singular.",
+                    UserWarning,
+                )
+                return None
+            dof = len(fvec) - len(sol)
+            if dof <= 0:
+                warnings.warn(
+                    "Warning: Could not estimate "
+                    "uncertainty as the number of fit parameters "
+                    "exceeds the number of data points.",
+                    UserWarning,
+                )
+                return None
+            par_names = list(self.pars)
+            cov_sol = hess_inv_sol * dot(fvec, fvec) / dof
+            unc = sqrt(diag(cov_sol))
+
+            return dict(zip(par_names, unc))
+        else:
+            return dot(fvec, fvec)
 
 
 # End class Refiner