#!/usr/bin/env python
##############################################################################
#
# diffpy.srfit by DANSE Diffraction group
# Simon J. L. Billinge
# (c) 2008 The Trustees of Columbia University
# in the City of New York. All rights reserved.
#
# File coded by: Chris Farrow, Pavol Juhas
#
# See AUTHORS.txt for a list of people who contributed.
# See LICENSE_DANSE.txt for license information.
#
##############################################################################
"""The FitResults and ContributionResults classes for storing results of
a fit.
The FitResults class is used to display the current state of a
FitRecipe. It stores the state, and uses it to calculate useful
statistics, which can be displayed on screen or saved to file.
"""
from __future__ import print_function
__all__ = ["FitResults", "ContributionResults", "initializeRecipe"]
import re
from collections import OrderedDict
import numpy
from diffpy.srfit.util import _DASHEDLINE
from diffpy.srfit.util import sortKeyForNumericString as numstr
from diffpy.srfit.util.inpututils import inputToString
from diffpy.utils._deprecator import build_deprecation_message, deprecated
fitresults_base = "diffpy.srfit.fitbase.FitResults"
removal_version = "4.0.0"
formatResults_dep_msg = build_deprecation_message(
fitresults_base,
"formatResults",
"get_results_string",
removal_version,
)
printResults_dep_msg = build_deprecation_message(
fitresults_base,
"printResults",
"print_results",
removal_version,
)
saveResults_dep_msg = build_deprecation_message(
fitresults_base,
"saveResults",
"save_results",
removal_version,
)
resultsDictionary_dep_msg = build_deprecation_message(
"diffpy.srfit.fitbase",
"resultsDictionary",
"get_results_dictionary",
removal_version,
new_base="diffpy.srfit.fitbase.FitResults",
)
initializeRecipe_dep_msg = build_deprecation_message(
"diffpy.srfit.fitbase",
"initializeRecipe",
"initialize_recipe_with_results",
removal_version,
new_base="diffpy.srfit.fitbase.FitRecipe",
)
[docs]
class FitResults(object):
"""Class for processing, presenting and storing results of a fit.
Attributes
----------
recipe : FitRecipe
The recipe from which the results were generated.
cov : numpy.ndarray or None
The covariance matrix of the refined variables. None if unavailable.
conresults : collections.OrderedDict[str, ContributionResults]
The ordered mapping of FitContribution name → ContributionResults.
derivstep : float
The fractional step size used for numerical derivatives (default 1e-8).
varnames : list[str]
The names of refined variables in the recipe.
varvals : numpy.ndarray
The optimized values of the refined variables.
varunc : numpy.ndarray or None
The estimated standard uncertainties of the variables. None if invalid.
showfixed : bool
The flag indicating whether to show the fixed variables in the
formatted output (default True).
fixednames : list[str]
The names of variables held fixed during refinement.
fixedvals : numpy.ndarray
The values of the fixed variables.
showcon : bool
The flag indicating whether to show the constrained parameters
in the formatted output (default False).
connames : list[str]
The names of constrained parameters.
convals : numpy.ndarray
The values of constrained parameters.
conunc : numpy.ndarray or None
The uncertainties of constrained parameters. None if unavailable.
residual : float
The scalar residual value of the recipe.
penalty : float
The penalty contribution to the residual from restraints.
chi2 : float
The chi-squared value of the fit.
cumchi2 : numpy.ndarray
The cumulative chi-squared as a function of data index.
rchi2 : float
The reduced chi-squared of the fit.
rw : float
The weighted R-factor of the fit.
cumrw : numpy.ndarray
The cumulative weighted R-factor as a function of data index.
messages : list[str]
The informational or warning messages associated with the results.
precision : int
The number of digits used when formatting numeric output (default 8).
_dcon : numpy.ndarray
The jacobian of constraint equations with respect to variables.
Used internally for uncertainty propagation.
Each of these attributes, except the recipe, are created or updated when
the update method is called.
"""
def __init__(self, recipe, update=True, showfixed=True, showcon=False):
"""Initialize the attributes.
Parameters
----------
recipe : FitRecipe
The recipe containing the results.
update : bool
The flag indicating whether to do an immediate update
(default True).
showfixed : bool
The flag indicating whether to show fixed variables in the
output (default True).
showcon : bool
The flag indicating whether to show constraint values in
the output (default False).
"""
self.recipe = recipe
self.conresults = OrderedDict()
self.derivstep = 1e-8
self.varnames = []
self.varvals = []
self.varunc = []
self.fixednames = []
self.fixedvals = []
self.connames = []
self.convals = []
self.conunc = []
self.cov = None
self.residual = 0
self.penalty = 0
self.chi2 = 0
self.rchi2 = 0
self.rw = 0
self.precision = 8
self._dcon = []
self.messages = []
self.showfixed = bool(showfixed)
self.showcon = bool(showcon)
if update:
self.update()
return
[docs]
def update(self):
"""Update the results according to the current state of the
recipe."""
# Note that the order of these operations are chosen to reduce
# computation time.
recipe = self.recipe
if not recipe._contributions:
return
# Make sure everything is ready for calculation
recipe._prepare()
# Store the variable names and values
self.varnames = recipe.get_names()
self.varvals = recipe.get_values()
fixedpars = recipe._tagmanager.union(recipe._fixedtag)
fixedpars = [p for p in fixedpars if not p.constrained]
self.fixednames = [p.name for p in fixedpars]
self.fixedvals = [p.value for p in fixedpars]
# Store the constraint information
self.connames = [con.par.name for con in recipe._oconstraints]
self.convals = [con.par.getValue() for con in recipe._oconstraints]
if self.varnames:
# Calculate the covariance
self._calculate_covariance()
# Get the variable uncertainties
self.varunc = [
self.cov[i, i] ** 0.5 for i in range(len(self.varnames))
]
# Get the constraint uncertainties
self._calculate_constraint_uncertainties()
# Store the fitting arrays and metrics for each FitContribution.
self.conresults = OrderedDict()
for con, weight in zip(
recipe._contributions.values(), recipe._weights
):
self.conresults[con.name] = ContributionResults(con, weight, self)
# Calculate the metrics
res = recipe.residual()
self.residual = numpy.dot(res, res)
self._calculate_metrics()
# Calculate the restraints penalty
w = self.chi2 / len(res)
self.penalty = sum([r.penalty(w) for r in recipe._restraintlist])
return
def _calculate_covariance(self):
"""Calculate the covariance matrix.
This is called by update. This code borrowed from PARK. It finds
the pseudo-inverse of the Jacobian using the singular value
decomposition.
"""
try:
J = self._calculate_jacobian()
u, s, vh = numpy.linalg.svd(J, 0)
self.cov = numpy.dot(vh.T.conj() / s**2, vh)
except numpy.linalg.LinAlgError:
self.messages.append("Cannot compute covariance matrix.")
lvarvals = len(self.varvals)
self.cov = numpy.zeros((lvarvals, lvarvals), dtype=float)
return
def _calculate_jacobian(self):
"""Calculate the Jacobian for the fitting.
Adapted from PARK. Returns the derivative wrt the fit variables
at point p.
This also calculates the derivatives of the constrained
parameters while we're at it.
Numeric derivatives are calculated based on step, where step is
the portion of variable value. E.g. step = dv/v.
"""
recipe = self.recipe
step = self.derivstep
# Make sure the input vector is an array
pvals = numpy.asarray(self.varvals)
# Compute the numeric derivative using the center point formula.
delta = step * pvals
# Center point formula:
# df/dv = lim_{h->0} ( f(v+h)-f(v-h) ) / ( 2h )
#
r = []
# The list of constraint derivatives with respect to variables
# The forward difference would be faster, but perhaps not as accurate.
conr = []
for k, v in enumerate(pvals):
h = delta[k]
pvals[k] = v + h
rk = self.recipe.residual(pvals)
# The constraints derivatives
cond = []
for con in recipe._oconstraints:
con.update()
cond.append(con.par.getValue())
pvals[k] = v - h
rk -= self.recipe.residual(pvals)
# FIXME - constraints are used for vectors as well!
for i, con in enumerate(recipe._oconstraints):
con.update()
val = con.par.getValue()
if numpy.isscalar(val):
cond[i] -= con.par.getValue()
cond[i] /= 2 * h
else:
cond[i] = 0.0
conr.append(cond)
pvals[k] = v
r.append(rk / (2 * h))
# Reset the constrained parameters to their original values
for con in recipe._oconstraints:
con.update()
self._dcon = numpy.vstack(conr).T
# return the jacobian
jac = numpy.vstack(r).T
return jac
def _calculate_metrics(self):
"""Calculate chi2, cumchi2, rchi2, rw and cumrw for the
recipe."""
cumchi2 = numpy.array([], dtype=float)
# total weighed denominator for the ratio in the Rw formula
yw2tot = 0.0
numpoints = 0
for con in self.conresults.values():
cc2w = con.weight * con.cumchi2
c2last = cumchi2[-1:].sum()
cumchi2 = numpy.concatenate([cumchi2, c2last + cc2w])
yw2tot += con.weight * (con.chi2 / con.rw**2)
numpoints += len(con.x)
chi2 = cumchi2[-1:].sum()
cumrw = numpy.sqrt(cumchi2 / yw2tot)
rw = cumrw[-1:].sum()
numpoints += len(self.recipe._restraintlist)
rchi2 = chi2 / (numpoints - len(self.varnames))
self.chi2 = chi2
self.rchi2 = rchi2
self.rw = rw
self.cumchi2 = cumchi2
self.cumrw = cumrw
return
def _calculate_constraint_uncertainties(self):
"""Calculate the uncertainty on the constrained parameters."""
vu = self.varunc
# sig^2(c) = sum_i sum_j sig(v_i) sig(v_j) (dc/dv_i)(dc/dv_j)
# sig^2(c) = sum_i sum_j [sig(v_i)(dc/dv_i)][sig(v_j)(dc/dv_j)]
# sig^2(c) = sum_i sum_j u_i u_j
self.conunc = []
for dci in self._dcon:
# Create sig(v_i) (dc/dv_i) array.
u = dci * vu
# The outer product is all possible pairings of u_i and u_j
# uu_ij = u_i u_j
uu = numpy.outer(u, u)
# Sum these pairings to get sig^2(c)
sig2c = sum(uu.flatten())
self.conunc.append(sig2c**0.5)
return
[docs]
def get_results_string(self, header="", footer="", update=False):
"""Format the results and return them in a string.
This function is called by ``print_results`` and
``save_results``. Overloading the formatting here will change
all three functions.
Parameters
----------
header : str
The header to add to the output (default "")
footer : str
The footer to add to the output (default "")
update : bool
The flag indicating whether to call ``update()`` (default
False).
Returns
-------
str
The string containing the formatted results.
"""
if update:
self.update()
lines = []
corrmin = 0.25
p = self.precision
pe = "%-" + "%i.%ie" % (p + 6, p)
pet = "%" + ".%ie" % (p,)
# Check to see if the uncertainty values are reliable.
certain = True
for con in self.conresults.values():
if (con.dy == 1).all():
certain = False
break
# User-defined header
if header:
lines.append(header)
if not certain:
err_msg = (
"Some quantities invalid due to missing profile uncertainty"
)
if err_msg not in self.messages:
self.messages.append(err_msg)
lines.extend(self.messages)
# Overall results
err_msg = "Overall"
if not certain:
err_msg += " (Chi2 and Reduced Chi2 invalid)"
lines.append(err_msg)
lines.append(_DASHEDLINE)
formatstr = "%-14s %.8f"
lines.append(formatstr % ("Residual", self.residual))
lines.append(
formatstr % ("Contributions", self.residual - self.penalty)
)
lines.append(formatstr % ("Restraints", self.penalty))
lines.append(formatstr % ("Chi2", self.chi2))
lines.append(formatstr % ("Reduced Chi2", self.rchi2))
lines.append(formatstr % ("Rw", self.rw))
# Per-FitContribution results
if len(self.conresults) > 1:
keys = list(self.conresults.keys())
keys.sort(key=numstr)
lines.append("")
err_msg = "Contributions"
if not certain:
err_msg += " (Chi2 and Reduced Chi2 invalid)"
lines.append(err_msg)
lines.append(_DASHEDLINE)
formatstr = "%-10s %-42.8f"
for name in keys:
res = self.conresults[name]
lines.append("")
namestr = name + " (%f)" % res.weight
lines.append(namestr)
lines.append("-" * len(namestr))
lines.append(formatstr % ("Residual", res.residual))
lines.append(formatstr % ("Chi2", res.chi2))
lines.append(formatstr % ("Rw", res.rw))
# The variables
if self.varnames:
lines.append("")
err_msg = "Variables"
if not certain:
err_msg2 = "Uncertainties invalid"
err_msg += " (%s)" % err_msg2
lines.append(err_msg)
lines.append(_DASHEDLINE)
varnames = self.varnames
varvals = self.varvals
varunc = self.varunc
varlines = []
w = max(map(len, varnames))
w = str(w + 1)
# Format the lines
formatstr = "%-" + w + "s " + pe + " +/- " + pet
for name, val, unc in zip(varnames, varvals, varunc):
varlines.append(formatstr % (name, val, unc))
varlines.sort()
lines.extend(varlines)
# Fixed variables
if self.showfixed and self.fixednames:
varlines = []
lines.append("")
lines.append("Fixed Variables")
lines.append(_DASHEDLINE)
w = max(map(len, self.fixednames))
w = str(w + 1)
formatstr = "%-" + w + "s " + pet
for name, val in zip(self.fixednames, self.fixedvals):
varlines.append(formatstr % (name, val))
varlines.sort()
lines.extend(varlines)
# The constraints
if self.connames and self.showcon:
lines.append("")
err_msg = "Constrained Parameters"
if not certain:
err_msg += " (Uncertainties invalid)"
lines.append(err_msg)
lines.append(_DASHEDLINE)
w = 0
keys = []
vals = {}
for con in self.conresults.values():
for i, loc in enumerate(con.conlocs):
names = [obj.name for obj in loc]
name = ".".join(names)
w = max(w, len(name))
val = con.convals[i]
unc = con.conunc[i]
keys.append(name)
vals[name] = (val, unc)
keys.sort(key=numstr)
w = str(w + 1)
formatstr = "%-" + w + "s %- 15f +/- %-15f"
for name in keys:
val, unc = vals[name]
lines.append(formatstr % (name, val, unc))
# Variable correlations
lines.append("")
corint = int(corrmin * 100)
err_msg = "Variable Correlations greater than %i%%" % corint
if not certain:
err_msg += " (Correlations invalid)"
lines.append(err_msg)
lines.append(_DASHEDLINE)
tup = []
cornames = []
n = len(self.varnames)
for i in range(n):
for j in range(i + 1, n):
name = "corr(%s, %s)" % (varnames[i], varnames[j])
val = self.cov[i, j] / (self.cov[i, i] * self.cov[j, j]) ** 0.5
if abs(val) > corrmin:
cornames.append(name)
tup.append((val, name))
tup.sort(key=lambda vn: abs(vn[0]))
tup.reverse()
if cornames:
w = max(map(len, cornames))
w = str(w + 1)
formatstr = "%-" + w + "s %.4f"
for val, name in tup:
lines.append(formatstr % (name, val))
else:
lines.append("No correlations greater than %i%%" % corint)
# User-defined footer
if footer:
lines.append(footer)
out = "\n".join(lines) + "\n"
return out
[docs]
def print_results(self, header="", footer="", update=False):
"""Format and print the results.
Parameters
----------
header : str
The header to add to the output (default "")
footer : str
The footer to add to the output (default "")
update : bool
The flag indicating whether to call ``update()`` (default
False).
"""
print(self.get_results_string(header, footer, update).rstrip())
return
[docs]
@deprecated(printResults_dep_msg)
def printResults(self, header="", footer="", update=False):
"""This function has been deprecated and will be removed in version
4.0.0.
Please use diffpy.srfit.fitbase.FitResults.print_results
instead.
"""
self.print_results(header, footer, update)
return
def __str__(self):
"""Return the formatted results string."""
return self.get_results_string()
[docs]
def save_results(self, filename, header="", footer="", update=False):
"""Format and save the results.
Parameters
----------
filename : str
The name of the save file.
header : str
The header to add to the output (default "")
footer : str
The footer to add to the output (default "")
update : bool
The flag indicating whether to call ``update()`` (default
False).
"""
# Save the time and user
from getpass import getuser
from time import ctime
myheader = "Results written: " + ctime() + "\n"
myheader += "produced by " + getuser() + "\n"
header = myheader + header
res = self.get_results_string(header, footer, update)
f = open(filename, "w")
f.write(res)
f.close()
return
[docs]
@deprecated(saveResults_dep_msg)
def saveResults(self, filename, header="", footer="", update=False):
"""This function has been deprecated and will be removed in version
4.0.0.
Please use diffpy.srfit.fitbase.FitResults.save_results
instead.
"""
self.save_results(filename, header, footer, update)
return
[docs]
def get_results_dictionary(self):
"""Get a dictionary of results, with variable names and values,
and overall metrics.
Returns
-------
results_dict : dict
The dictionary containing the variable names and values,
and overall metrics, from the FitResults.
"""
parameter_names = self.varnames
parameter_values = self.varvals
results_dict = dict(zip(parameter_names, parameter_values))
results_dict.update(
{
"Residual": self.residual,
"Contributions": self.residual - self.penalty,
"Restraints": self.penalty,
"Chi2": self.chi2,
"Reduced Chi2": self.rchi2,
"Rw": self.rw,
}
)
return results_dict
# End class FitResults
[docs]
class ContributionResults(object):
"""Class for processing, storing FitContribution results.
This does not store the FitContribution.
Attributes
----------
y : numpy.ndarray or None
The FitContribution's profile over the calculation range
(default None).
dy : numpy.ndarray or None
The uncertainty in the FitContribution's profile over the
calculation range (default None).
x : numpy.ndarray or None
The numpy array of the calculated independent variable for the
FitContribution (default None).
ycalc : numpy.ndarray or None
The numpy array of the calculated signal for the FitContribution
(default None).
residual : float
The scalar residual of the FitContribution.
chi2 : float
The chi2 of the FitContribution.
cumchi2 : numpy.ndarray
The cumulative chi2 of the FitContribution.
rw : float
The Rw of the FitContribution.
cumrw : numpy.ndarray
The cumulative Rw of the FitContribution.
weight : float
The weight of the FitContribution in the recipe.
conlocs : list
The location of the constrained parameters in the
FitContribution (see the
RecipeContainer._locate_managed_object method).
convals : list
The values of the constrained parameters.
conunc : list
The uncertainties in the constraint values.
"""
def __init__(self, con, weight, fitres):
"""Initialize the attributes.
Parameters
----------
con : FitContribution
The FitContribution to summarize.
weight : float
The weight of the FitContribution in the recipe.
fitres : FitResults
The FitResults instance containing this
ContributionResults.
"""
self.x = None
self.y = None
self.dy = None
self.ycalc = None
self.residual = 0
self.chi2 = 0
self.rw = 0
self.weight = 0
self.conlocs = []
self.convals = []
self.conunc = []
self._init(con, weight, fitres)
return
def _init(self, con, weight, fitres):
"""Initialize the attributes, for real."""
# Note that the order of these operations is chosen to reduce
# computation time.
if con.profile is None:
return
recipe = fitres.recipe
# Store the weight
self.weight = weight
# First the residual
res = con.residual()
self.residual = numpy.dot(res, res)
# The arrays
self.x = numpy.array(con.profile.x)
self.y = numpy.array(con.profile.y)
self.dy = numpy.array(con.profile.dy)
self.ycalc = numpy.array(con.profile.ycalc)
# The other metrics
self._calculate_metrics()
# Find the parameters
for i, constraint in enumerate(recipe._oconstraints):
par = constraint.par
loc = con._locate_managed_object(par)
if loc:
self.conlocs.append(loc)
self.convals.append(fitres.convals[i])
self.conunc.append(fitres.conunc[i])
return
# FIXME: factor rw, chi2, cumrw, cumchi2 to separate functions.
def _calculate_metrics(self):
"""Calculate chi2 and Rw of the recipe."""
# We take absolute values in case the signal is complex
num = numpy.abs(self.y - self.ycalc)
y = numpy.abs(self.y)
chiv = num / self.dy
self.cumchi2 = numpy.cumsum(chiv**2)
# avoid index error for empty array
self.chi2 = self.cumchi2[-1:].sum()
yw = y / self.dy
yw2tot = numpy.dot(yw, yw)
if yw2tot == 0.0:
yw2tot = 1.0
self.cumrw = numpy.sqrt(self.cumchi2 / yw2tot)
# avoid index error for empty array
self.rw = self.cumrw[-1:].sum()
return
# End class ContributionResults
@deprecated(resultsDictionary_dep_msg)
def resultsDictionary(results):
"""This function has been deprecated and will be removed in version
4.0.0.
Please use
diffpy.srfit.fitbase.FitResults.get_results_dictionary instead.
Get dictionary of results from file.
This reads the results from file and stores them in a dictionary to be
returned to the caller. The dictionary may contain non-result entries.
Parameters
----------
results : str or file-like
The open file-like object, name of a file that contains
results from FitResults, or a string containing fit results.
Returns
-------
dict
The mapping of result names to their string values.
"""
resstr = inputToString(results)
rx = {
"f": r"[+-]? *(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?",
"n": r"[a-zA-Z_]\w*",
}
pat = r"(%(n)s)\s+(%(f)s)" % rx
matches = re.findall(pat, resstr)
# Prefer the first match
matches.reverse()
mpairs = dict(matches)
return mpairs
[docs]
@deprecated(initializeRecipe_dep_msg)
def initializeRecipe(recipe, results):
"""This function has been deprecated and will be removed in version
4.0.0.
Please use
diffpy.srfit.fitbase.FitRecipe.initialize_recipe_with_results
instead.
Initialize the variables of a recipe from a results file.
This reads the results from file and initializes any variables (fixed or
free) in the recipe to the results values. Note that the recipe has to be
configured, with variables. This does not reconstruct a FitRecipe.
Parameters
----------
recipe : FitRecipe
The configured recipe with variables.
results : str or file-like
The open file-like object, name of a file that contains
results from FitResults, or a string containing fit results.
Raises
------
AttributeError
If no results can be found in ``results``.
"""
mpairs = resultsDictionary(results)
if not mpairs:
raise AttributeError("Cannot find results")
# Get variable names
names = recipe._parameters.keys()
for vname in names:
value = mpairs.get(vname)
if value is not None:
var = recipe.get(vname)
var.value = float(value)
return