diff --git a/src/diffpy/morph/morph_io.py b/src/diffpy/morph/morph_io.py index cf194f5..680d25c 100644 --- a/src/diffpy/morph/morph_io.py +++ b/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(): diff --git a/src/diffpy/morph/morphapp.py b/src/diffpy/morph/morphapp.py index 5b03a6a..680b641 100755 --- a/src/diffpy/morph/morphapp.py +++ b/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,7 +746,7 @@ 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 @@ -738,9 +754,12 @@ def single_morph( # 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, diff --git a/src/diffpy/morph/refine.py b/src/diffpy/morph/refine.py index 8f1390d..29a7909 100644 --- a/src/diffpy/morph/refine.py +++ b/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