https://github.com/lmfit/lmfit-py
Tip revision: 555c29e63c18262a0b3ab9648c9d2e0d56f7a023 authored by Matthew Newville on 08 July 2023, 15:51:20 UTC
spellchecked whatsnew
spellchecked whatsnew
Tip revision: 555c29e
model.py
"""Implementation of the Model interface."""
from copy import deepcopy
from functools import wraps
import inspect
import json
import operator
import warnings
from asteval import valid_symbol_name
import numpy as np
from scipy.special import erf
from scipy.stats import t
import lmfit
from . import Minimizer, Parameter, Parameters, lineshapes
from .confidence import conf_interval
from .jsonutils import HAS_DILL, decode4js, encode4js
from .minimizer import MinimizerResult
from .printfuncs import ci_report, fit_report, fitreport_html_table
tiny = 1.e-15
# Use pandas.isnull for aligning missing data if pandas is available.
# otherwise use numpy.isnan
try:
from pandas import Series, isnull
except ImportError:
isnull = np.isnan
Series = type(NotImplemented)
def _align(var, mask, data):
"""Align missing data, if pandas is available."""
if isinstance(data, Series) and isinstance(var, Series):
return var.reindex_like(data).dropna()
elif mask is not None:
return var[mask]
return var
try:
import matplotlib # noqa: F401
_HAS_MATPLOTLIB = True
except Exception:
_HAS_MATPLOTLIB = False
def _ensureMatplotlib(function):
if _HAS_MATPLOTLIB:
@wraps(function)
def wrapper(*args, **kws):
return function(*args, **kws)
return wrapper
def no_op(*args, **kwargs):
print('matplotlib module is required for plotting the results')
return no_op
def get_reducer(option):
"""Factory function to build a parser for complex numbers.
Parameters
----------
option : {'real', 'imag', 'abs', 'angle'}
Implements the NumPy function with the same name.
Returns
-------
callable
See docstring for `reducer` below.
"""
if option not in ['real', 'imag', 'abs', 'angle']:
raise ValueError(f"Invalid option ('{option}') for function 'propagate_err'.")
def reducer(array):
"""Convert a complex array to a real array.
Several conversion methods are available and it does nothing to a
purely real array.
Parameters
----------
array : array-like
Input array. If complex, will be converted to real array via
one of the following NumPy functions: :numpydoc:`real`,
:numpydoc:`imag`, :numpydoc:`abs`, or :numpydoc:`angle`.
Returns
-------
numpy.array
Returned array will be purely real.
"""
if any(np.iscomplex(array)):
parsed_array = getattr(np, option)(array)
else:
parsed_array = array
return parsed_array
return reducer
def propagate_err(z, dz, option):
"""Perform error propagation on a vector of complex uncertainties.
Required to get values for magnitude (abs) and phase (angle)
uncertainty.
Parameters
----------
z : array-like
Array of complex or real numbers.
dz : array-like
Array of uncertainties corresponding to `z`. Must satisfy
``numpy.shape(dz) == numpy.shape(z)``.
option : {'real', 'imag', 'abs', 'angle'}
How to convert the array `z` to an array with real numbers.
Returns
-------
numpy.array
Returned array will be purely real.
Notes
-----
Uncertainties are ``1/weights``. If the weights provided are real,
they are assumed to apply equally to the real and imaginary parts. If
the weights are complex, the real part of the weights are applied to
the real part of the residual and the imaginary part is treated
correspondingly.
In the case where ``option='angle'`` and ``numpy.abs(z) == 0`` for any
value of `z` the phase angle uncertainty becomes the entire circle and
so a value of `math:pi` is returned.
In the case where ``option='abs'`` and ``numpy.abs(z) == 0`` for any
value of `z` the magnitude uncertainty is approximated by
``numpy.abs(dz)`` for that value.
"""
if option not in ['real', 'imag', 'abs', 'angle']:
raise ValueError(f"Invalid option ('{option}') for function 'propagate_err'.")
if z.shape != dz.shape:
raise ValueError(f"shape of z: {z.shape} != shape of dz: {dz.shape}")
# Check the main vector for complex. Do nothing if real.
if any(np.iscomplex(z)):
# if uncertainties are real, apply them equally to
# real and imaginary parts
if all(np.isreal(dz)):
dz = dz+1j*dz
if option == 'real':
err = np.real(dz)
elif option == 'imag':
err = np.imag(dz)
elif option in ['abs', 'angle']:
rz = np.real(z)
iz = np.imag(z)
rdz = np.real(dz)
idz = np.imag(dz)
# Don't spit out warnings for divide by zero. Will fix these later.
with np.errstate(divide='ignore', invalid='ignore'):
if option == 'abs':
# Standard error propagation for abs = sqrt(re**2 + im**2)
err = np.true_divide(np.sqrt((iz*idz)**2+(rz*rdz)**2),
np.abs(z))
# For abs = 0, error is +/- abs(rdz + j idz)
err[err == np.inf] = np.abs(dz)[err == np.inf]
if option == 'angle':
# Standard error propagation for angle = arctan(im/re)
err = np.true_divide(np.sqrt((rz*idz)**2+(iz*rdz)**2),
np.abs(z)**2)
# For abs = 0, error is +/- pi (i.e. the whole circle)
err[err == np.inf] = np.pi
else:
err = dz
return err
def coerce_arraylike(x):
"""
coerce lists, tuples, and pandas Series, hdf5 Groups, etc to an
ndarray float64 or complex128, but leave other data structures
and objects unchanged
"""
if isinstance(x, (list, tuple, Series)) or hasattr(x, '__array__'):
if np.isrealobj(x):
return np.asfarray(x)
if np.iscomplexobj(x):
return np.asfarray(x, dtype=np.complex128)
return x
class Model:
"""Create a model from a user-supplied model function."""
_forbidden_args = ('data', 'weights', 'params')
_invalid_ivar = "Invalid independent variable name ('%s') for function %s"
_invalid_par = "Invalid parameter name ('%s') for function %s"
_invalid_hint = "unknown parameter hint '%s' for param '%s'"
_hint_names = ('value', 'vary', 'min', 'max', 'expr')
valid_forms = ()
def __init__(self, func, independent_vars=None, param_names=None,
nan_policy='raise', prefix='', name=None, **kws):
"""
The model function will normally take an independent variable
(generally, the first argument) and a series of arguments that are
meant to be parameters for the model. It will return an array of
data to model some data as for a curve-fitting problem.
Parameters
----------
func : callable
Function to be wrapped.
independent_vars : :obj:`list` of :obj:`str`, optional
Arguments to `func` that are independent variables (default is
None).
param_names : :obj:`list` of :obj:`str`, optional
Names of arguments to `func` that are to be made into
parameters (default is None).
nan_policy : {'raise', 'propagate', 'omit'}, optional
How to handle NaN and missing values in data. See Notes below.
prefix : str, optional
Prefix used for the model.
name : str, optional
Name for the model. When None (default) the name is the same
as the model function (`func`).
**kws : dict, optional
Additional keyword arguments to pass to model function.
Notes
-----
1. Parameter names are inferred from the function arguments, and a
residual function is automatically constructed.
2. The model function must return an array that will be the same
size as the data being modeled.
3. `nan_policy` sets what to do when a NaN or missing value is
seen in the data. Should be one of:
- `'raise'` : raise a `ValueError` (default)
- `'propagate'` : do nothing
- `'omit'` : drop missing data
Examples
--------
The model function will normally take an independent variable
(generally, the first argument) and a series of arguments that are
meant to be parameters for the model. Thus, a simple peak using a
Gaussian defined as:
>>> import numpy as np
>>> def gaussian(x, amp, cen, wid):
... return amp * np.exp(-(x-cen)**2 / wid)
can be turned into a Model with:
>>> gmodel = Model(gaussian)
this will automatically discover the names of the independent
variables and parameters:
>>> print(gmodel.param_names, gmodel.independent_vars)
['amp', 'cen', 'wid'], ['x']
"""
self.func = func
if not isinstance(prefix, str):
prefix = ''
if len(prefix) > 0 and not valid_symbol_name(prefix):
raise ValueError(f"'{prefix}' is not a valid Model prefix")
self._prefix = prefix
self._param_root_names = param_names # will not include prefixes
self.independent_vars = independent_vars
self._func_allargs = []
self._func_haskeywords = False
self.nan_policy = nan_policy
self.opts = kws
# the following has been changed from OrderedSet for the time being
self.param_hints = {}
self._param_names = []
self._parse_params()
if self.independent_vars is None:
self.independent_vars = []
if name is None and hasattr(self.func, '__name__'):
name = self.func.__name__
self._name = name
def _reprstring(self, long=False):
out = self._name
opts = []
if len(self._prefix) > 0:
opts.append(f"prefix='{self._prefix}'")
if long:
for k, v in self.opts.items():
opts.append(f"{k}='{v}'")
if len(opts) > 0:
out = f"{out}, {', '.join(opts)}"
return f"Model({out})"
def _get_state(self):
"""Save a Model for serialization.
Note: like the standard-ish '__getstate__' method but not really
useful with Pickle.
"""
funcdef = None
if HAS_DILL:
funcdef = self.func
if self.func.__name__ == '_eval':
funcdef = self.expr
state = (self.func.__name__, funcdef, self._name, self._prefix,
self.independent_vars, self._param_root_names,
self.param_hints, self.nan_policy, self.opts)
return (state, None, None)
def _set_state(self, state, funcdefs=None):
"""Restore Model from serialization.
Note: like the standard-ish '__setstate__' method but not really
useful with Pickle.
Parameters
----------
state
Serialized state from `_get_state`.
funcdefs : dict, optional
Dictionary of function definitions to use to construct Model.
"""
return _buildmodel(state, funcdefs=funcdefs)
def dumps(self, **kws):
"""Dump serialization of Model as a JSON string.
Parameters
----------
**kws : optional
Keyword arguments that are passed to `json.dumps`.
Returns
-------
str
JSON string representation of Model.
See Also
--------
loads, json.dumps
"""
return json.dumps(encode4js(self._get_state()), **kws)
def dump(self, fp, **kws):
"""Dump serialization of Model to a file.
Parameters
----------
fp : file-like object
An open and `.write()`-supporting file-like object.
**kws : optional
Keyword arguments that are passed to `json.dumps`.
Returns
-------
int
Return value from `fp.write()`: the number of characters
written.
See Also
--------
dumps, load, json.dump
"""
return fp.write(self.dumps(**kws))
def loads(self, s, funcdefs=None, **kws):
"""Load Model from a JSON string.
Parameters
----------
s : str
Input JSON string containing serialized Model.
funcdefs : dict, optional
Dictionary of function definitions to use to construct Model.
**kws : optional
Keyword arguments that are passed to `json.loads`.
Returns
-------
Model
Model created from JSON string.
See Also
--------
dump, dumps, load, json.loads
"""
tmp = decode4js(json.loads(s, **kws))
return self._set_state(tmp, funcdefs=funcdefs)
def load(self, fp, funcdefs=None, **kws):
"""Load JSON representation of Model from a file-like object.
Parameters
----------
fp : file-like object
An open and `.read()`-supporting file-like object.
funcdefs : dict, optional
Dictionary of function definitions to use to construct Model.
**kws : optional
Keyword arguments that are passed to `loads`.
Returns
-------
Model
Model created from `fp`.
See Also
--------
dump, loads, json.load
"""
return self.loads(fp.read(), funcdefs=funcdefs, **kws)
@property
def name(self):
"""Return Model name."""
return self._reprstring(long=False)
@name.setter
def name(self, value):
self._name = value
@property
def prefix(self):
"""Return Model prefix."""
return self._prefix
@prefix.setter
def prefix(self, value):
"""Change Model prefix."""
self._prefix = value
self._set_paramhints_prefix()
self._param_names = []
self._parse_params()
def _set_paramhints_prefix(self):
"""Reset parameter hints for prefix: intended to be overwritten."""
@property
def param_names(self):
"""Return the parameter names of the Model."""
return self._param_names
def __repr__(self):
"""Return representation of Model."""
return f"<lmfit.Model: {self.name}>"
def copy(self, **kwargs):
"""DOES NOT WORK."""
raise NotImplementedError("Model.copy does not work. Make a new Model")
def _parse_params(self):
"""Build parameters from function arguments."""
if self.func is None:
return
kw_args = {}
keywords_ = None
# need to fetch the following from the function signature:
# pos_args: list of positional argument names
# kw_args: dict of keyword arguments with default values
# keywords_: name of **kws argument or None
# 1. limited support for asteval functions as the model functions:
if hasattr(self.func, 'argnames') and hasattr(self.func, 'kwargs'):
pos_args = self.func.argnames[:]
for name, defval in self.func.kwargs:
kw_args[name] = defval
# 2. modern, best-practice approach: use inspect.signature
else:
pos_args = []
sig = inspect.signature(self.func)
for fnam, fpar in sig.parameters.items():
if fpar.kind == fpar.VAR_KEYWORD:
keywords_ = fnam
elif fpar.kind == fpar.POSITIONAL_OR_KEYWORD:
if fpar.default == fpar.empty:
pos_args.append(fnam)
else:
kw_args[fnam] = fpar.default
elif fpar.kind == fpar.VAR_POSITIONAL:
raise ValueError(f"varargs '*{fnam}' is not supported")
# inspection done
self._func_haskeywords = keywords_ is not None
self._func_allargs = pos_args + list(kw_args.keys())
allargs = self._func_allargs
if len(allargs) == 0 and keywords_ is not None:
return
# default independent_var = 1st argument
if self.independent_vars is None:
self.independent_vars = [pos_args[0]]
# default param names: all positional args
# except independent variables
self.def_vals = {}
might_be_param = []
if self._param_root_names is None:
self._param_root_names = pos_args[:]
for key, val in kw_args.items():
if (not isinstance(val, bool) and
isinstance(val, (float, int))):
self._param_root_names.append(key)
self.def_vals[key] = val
elif val is None:
might_be_param.append(key)
for p in self.independent_vars:
if p in self._param_root_names:
self._param_root_names.remove(p)
new_opts = {}
for opt, val in self.opts.items():
if (opt in self._param_root_names or opt in might_be_param and
isinstance(val, Parameter)):
self.set_param_hint(opt, value=val.value,
min=val.min, max=val.max, expr=val.expr)
elif opt in self._func_allargs:
new_opts[opt] = val
self.opts = new_opts
if self._prefix is None:
self._prefix = ''
names = [f"{self._prefix}{pname}" for pname in self._param_root_names]
# check variables names for validity
# The implicit magic in fit() requires us to disallow some
fname = self.func.__name__
for arg in self.independent_vars:
if arg not in allargs or arg in self._forbidden_args:
raise ValueError(self._invalid_ivar % (arg, fname))
for arg in names:
if (self._strip_prefix(arg) not in allargs or
arg in self._forbidden_args):
raise ValueError(self._invalid_par % (arg, fname))
# the following as been changed from OrderedSet for the time being.
self._param_names = names[:]
def set_param_hint(self, name, **kwargs):
"""Set *hints* to use when creating parameters with `make_params()`.
The given hint can include optional bounds and constraints
``(value, vary, min, max, expr)``, which will be used by
`Model.make_params()` when building default parameters.
While this can be used to set initial values, `Model.make_params` or
the function `create_params` should be preferred for creating
parameters with initial values.
The intended use here is to control how a Model should create
parameters, such as setting bounds that are required by the mathematics
of the model (for example, that a peak width cannot be negative), or to
define common constrained parameters.
Parameters
----------
name : str
Parameter name, can include the models `prefix` or not.
**kwargs : optional
Arbitrary keyword arguments, needs to be a Parameter attribute.
Can be any of the following:
- value : float, optional
Numerical Parameter value.
- vary : bool, optional
Whether the Parameter is varied during a fit (default is
True).
- min : float, optional
Lower bound for value (default is ``-numpy.inf``, no lower
bound).
- max : float, optional
Upper bound for value (default is ``numpy.inf``, no upper
bound).
- expr : str, optional
Mathematical expression used to constrain the value during
the fit.
Example
--------
>>> model = GaussianModel()
>>> model.set_param_hint('sigma', min=0)
"""
npref = len(self._prefix)
if npref > 0 and name.startswith(self._prefix):
name = name[npref:]
if name not in self.param_hints:
self.param_hints[name] = {}
for key, val in kwargs.items():
if key in self._hint_names:
self.param_hints[name][key] = val
else:
warnings.warn(self._invalid_hint % (key, name))
def print_param_hints(self, colwidth=8):
"""Print a nicely aligned text-table of parameter hints.
Parameters
----------
colwidth : int, optional
Width of each column, except for first and last columns.
"""
name_len = max(len(s) for s in self.param_hints)
print('{:{name_len}} {:>{n}} {:>{n}} {:>{n}} {:>{n}} {:{n}}'
.format('Name', 'Value', 'Min', 'Max', 'Vary', 'Expr',
name_len=name_len, n=colwidth))
line = ('{name:<{name_len}} {value:{n}g} {min:{n}g} {max:{n}g} '
'{vary!s:>{n}} {expr}')
for name, values in sorted(self.param_hints.items()):
pvalues = dict(name=name, value=np.nan, min=-np.inf, max=np.inf,
vary=True, expr='')
pvalues.update(**values)
print(line.format(name_len=name_len, n=colwidth, **pvalues))
def make_params(self, verbose=False, **kwargs):
"""Create a Parameters object for a Model.
Parameters
----------
verbose : bool, optional
Whether to print out messages (default is False).
**kwargs : optional
Parameter names and initial values or dictionaries of
values and attributes.
Returns
---------
params : Parameters
Parameters object for the Model.
Notes
-----
1. Parameter values can be numbers (floats or ints) to set the parameter
value, or dictionaries with any of the following keywords:
``value``, ``vary``, ``min``, ``max``, ``expr``, ``brute_step``,
``is_init_value`` to set those parameter attributes.
2. This method will also apply any default values or parameter hints
that may have been defined for the model.
Example
--------
>>> gmodel = GaussianModel(prefix='peak_') + LinearModel(prefix='bkg_')
>>> gmodel.make_params(peak_center=3200, bkg_offset=0, bkg_slope=0,
... peak_amplitdue=dict(value=100, min=2),
... peak_sigma=dict(value=25, min=0, max=1000))
"""
params = Parameters()
def setpar(par, val):
# val is expected to be float-like or a dict: must have 'value' or 'expr' key
if isinstance(val, dict):
dval = val
else:
dval = {'value': float(val)}
if len(dval) < 1 or not ('value' in dval or 'expr' in dval):
raise TypeError(f'Invalid parameter value for {par}: {val}')
par.set(**dval)
# make sure that all named parameters are in params
for name in self.param_names:
if name in params:
par = params[name]
else:
par = Parameter(name=name)
par._delay_asteval = True
basename = name[len(self._prefix):]
# apply defaults from model function definition
if basename in self.def_vals:
par.value = self.def_vals[basename]
if par.value in (None, -np.inf, np.inf, np.nan):
for key, val in self.def_vals.items():
if key in name.lower():
par.value = val
# apply defaults from parameter hints
if basename in self.param_hints:
hint = self.param_hints[basename]
for item in self._hint_names:
if item in hint:
setattr(par, item, hint[item])
# apply values passed in through kw args
if basename in kwargs:
setpar(par, kwargs[basename])
if name in kwargs:
setpar(par, kwargs[name])
params.add(par)
if verbose:
print(f' - Adding parameter "{name}"')
# next build parameters defined in param_hints
# note that composites may define their own additional
# convenience parameters here
for basename, hint in self.param_hints.items():
name = f"{self._prefix}{basename}"
if name in params:
par = params[name]
else:
par = Parameter(name=name)
params.add(par)
if verbose:
print(f' - Adding parameter for hint "{name}"')
par._delay_asteval = True
for item in self._hint_names:
if item in hint:
setattr(par, item, hint[item])
if basename in kwargs:
setpar(par, kwargs[basename])
# Add the new parameter to self._param_names
if name not in self._param_names:
self._param_names.append(name)
for p in params.values():
p._delay_asteval = False
return params
def guess(self, data, x, **kws):
"""Guess starting values for the parameters of a Model.
This is not implemented for all models, but is available for many
of the built-in models.
Parameters
----------
data : array_like
Array of data (i.e., y-values) to use to guess parameter values.
x : array_like
Array of values for the independent variable (i.e., x-values).
**kws : optional
Additional keyword arguments, passed to model function.
Returns
-------
Parameters
Initial, guessed values for the parameters of a Model.
Raises
------
NotImplementedError
If the `guess` method is not implemented for a Model.
Notes
-----
Should be implemented for each model subclass to run
`self.make_params()`, update starting values and return a
Parameters object.
.. versionchanged:: 1.0.3
Argument ``x`` is now explicitly required to estimate starting values.
"""
cname = self.__class__.__name__
msg = f'guess() not implemented for {cname}'
raise NotImplementedError(msg)
def _residual(self, params, data, weights, **kwargs):
"""Return the residual.
Default residual: ``(data-model)*weights``.
If the model returns complex values, the residual is computed by
treating the real and imaginary parts separately. In this case, if
the weights provided are real, they are assumed to apply equally
to the real and imaginary parts. If the weights are complex, the
real part of the weights are applied to the real part of the
residual and the imaginary part is treated correspondingly.
Since the underlying `scipy.optimize` routines expect
``numpy.float`` arrays, the only complex type supported is
``complex``.
The "ravels" throughout are necessary to support `pandas.Series`.
"""
model = self.eval(params, **kwargs)
if self.nan_policy == 'raise' and not np.all(np.isfinite(model)):
msg = ('The model function generated NaN values and the fit '
'aborted! Please check your model function and/or set '
'boundaries on parameters where applicable. In cases like '
'this, using "nan_policy=\'omit\'" will probably not work.')
raise ValueError(msg)
diff = model - data
if diff.dtype == complex:
# data/model are complex
diff = diff.ravel().view(float)
if weights is not None:
if weights.dtype == complex:
# weights are complex
weights = weights.ravel().view(float)
else:
# real weights but complex data
weights = (weights + 1j * weights).ravel().view(float)
if weights is not None:
diff *= weights
return diff
def _strip_prefix(self, name):
npref = len(self._prefix)
if npref > 0 and name.startswith(self._prefix):
name = name[npref:]
return name
def make_funcargs(self, params=None, kwargs=None, strip=True):
"""Convert parameter values and keywords to function arguments."""
if params is None:
params = {}
if kwargs is None:
kwargs = {}
out = {}
out.update(self.opts)
# 0: if a keyword argument is going to overwrite a parameter,
# save that value so it can be restored before returning
saved_values = {}
for name, val in kwargs.items():
if name in params:
saved_values[name] = params[name].value
params[name].value = val
if len(saved_values) > 0:
params.update_constraints()
# 1. fill in in all parameter values
for name, par in params.items():
if strip:
name = self._strip_prefix(name)
if name in self._func_allargs or self._func_haskeywords:
out[name] = par.value
# 2. for each function argument, use 'prefix+varname' in params,
# avoiding possible name collisions with unprefixed params
if len(self._prefix) > 0:
for fullname in self._param_names:
if fullname in params:
name = self._strip_prefix(fullname)
if name in self._func_allargs or self._func_haskeywords:
out[name] = params[fullname].value
# 3. kwargs might directly update function arguments
for name, val in kwargs.items():
if strip:
name = self._strip_prefix(name)
if name in self._func_allargs or self._func_haskeywords:
out[name] = val
# 4. finally, reset any values that have overwritten parameter values
for name, val in saved_values.items():
params[name].value = val
return out
def post_fit(self, fitresult):
"""function that is called just after fit, can be overloaded by
subclasses to add non-fitting 'calculated parameters'
"""
pass
def _make_all_args(self, params=None, **kwargs):
"""Generate **all** function args for all functions."""
args = {}
for key, val in self.make_funcargs(params, kwargs).items():
args[f"{self._prefix}{key}"] = val
return args
def eval(self, params=None, **kwargs):
"""Evaluate the model with supplied parameters and keyword arguments.
Parameters
-----------
params : Parameters, optional
Parameters to use in Model.
**kwargs : optional
Additional keyword arguments to pass to model function.
Returns
-------
numpy.ndarray, float, int or complex
Value of model given the parameters and other arguments.
Notes
-----
1. if `params` is None, the values for all parameters are expected
to be provided as keyword arguments.
2. If `params` is given, and a keyword argument for a parameter value
is also given, the keyword argument will be used in place of the value
in the value in `params`.
3. all non-parameter arguments for the model function, **including
all the independent variables** will need to be passed in using
keyword arguments.
4. The return types are generally `numpy.ndarray`, but may depends on
the model function and input independent variables. That is, return
values may be Python `float`, `int`, or `complex` values.
"""
return coerce_arraylike(self.func(**self.make_funcargs(params, kwargs)))
@property
def components(self):
"""Return components for composite model."""
return [self]
def eval_components(self, params=None, **kwargs):
"""Evaluate the model with the supplied parameters.
Parameters
-----------
params : Parameters, optional
Parameters to use in Model.
**kwargs : optional
Additional keyword arguments to pass to model function.
Returns
-------
dict
Keys are prefixes for component model, values are value of
each component.
"""
key = self._prefix
if len(key) < 1:
key = self._name
return {key: self.eval(params=params, **kwargs)}
def fit(self, data, params=None, weights=None, method='leastsq',
iter_cb=None, scale_covar=True, verbose=False, fit_kws=None,
nan_policy=None, calc_covar=True, max_nfev=None,
coerce_farray=True, **kwargs):
"""Fit the model to the data using the supplied Parameters.
Parameters
----------
data : array_like
Array of data to be fit.
params : Parameters, optional
Parameters to use in fit (default is None).
weights : array_like, optional
Weights to use for the calculation of the fit residual [i.e.,
`weights*(data-fit)`]. Default is None; must have the same size as
`data`.
method : str, optional
Name of fitting method to use (default is `'leastsq'`).
iter_cb : callable, optional
Callback function to call at each iteration (default is None).
scale_covar : bool, optional
Whether to automatically scale the covariance matrix when
calculating uncertainties (default is True).
verbose : bool, optional
Whether to print a message when a new parameter is added
because of a hint (default is True).
fit_kws : dict, optional
Options to pass to the minimizer being used.
nan_policy : {'raise', 'propagate', 'omit'}, optional
What to do when encountering NaNs when fitting Model.
calc_covar : bool, optional
Whether to calculate the covariance matrix (default is True)
for solvers other than `'leastsq'` and `'least_squares'`.
Requires the ``numdifftools`` package to be installed.
max_nfev : int or None, optional
Maximum number of function evaluations (default is None). The
default value depends on the fitting method.
coerce_farray : bool, optional
Whether to coerce data and independent data to be ndarrays
with dtype of float64 (or complex128). If set to False, data
and independent data are not coerced at all, but the output of
the model function will be. (default is True)
**kwargs : optional
Arguments to pass to the model function, possibly overriding
parameters.
Returns
-------
ModelResult
Notes
-----
1. if `params` is None, the values for all parameters are expected
to be provided as keyword arguments. Mixing `params` and
keyword arguments is deprecated (see `Model.eval`).
2. all non-parameter arguments for the model function, **including
all the independent variables** will need to be passed in using
keyword arguments.
3. Parameters are copied on input, so that the original Parameter objects
are unchanged, and the updated values are in the returned `ModelResult`.
Examples
--------
Take ``t`` to be the independent variable and data to be the curve
we will fit. Use keyword arguments to set initial guesses:
>>> result = my_model.fit(data, tau=5, N=3, t=t)
Or, for more control, pass a Parameters object.
>>> result = my_model.fit(data, params, t=t)
"""
if params is None:
params = self.make_params(verbose=verbose)
else:
params = deepcopy(params)
# If any kwargs match parameter names, override params.
param_kwargs = set(kwargs.keys()) & set(self.param_names)
for name in param_kwargs:
p = kwargs[name]
if isinstance(p, Parameter):
p.name = name # allows N=Parameter(value=5) with implicit name
params[name] = deepcopy(p)
else:
params[name].set(value=p)
del kwargs[name]
# All remaining kwargs should correspond to independent variables.
for name in kwargs:
if name not in self.independent_vars:
warnings.warn(f"The keyword argument {name} does not " +
"match any arguments of the model function. " +
"It will be ignored.", UserWarning)
# If any parameter is not initialized raise a more helpful error.
missing_param = any(p not in params.keys() for p in self.param_names)
blank_param = any((p.value is None and p.expr is None)
for p in params.values())
if missing_param or blank_param:
msg = ('Assign each parameter an initial value by passing '
'Parameters or keyword arguments to fit.\n')
missing = [p for p in self.param_names if p not in params.keys()]
blank = [name for name, p in params.items()
if p.value is None and p.expr is None]
msg += f'Missing parameters: {str(missing)}\n'
msg += f'Non initialized parameters: {str(blank)}'
raise ValueError(msg)
# Handle null/missing values.
if nan_policy is not None:
self.nan_policy = nan_policy
mask = None
if self.nan_policy == 'omit':
mask = ~isnull(data)
if mask is not None:
data = data[mask]
if weights is not None:
weights = _align(weights, mask, data)
# If independent_vars and data are alignable (pandas), align them,
# and apply the mask from above if there is one.
for var in self.independent_vars:
if not np.isscalar(kwargs[var]):
kwargs[var] = _align(kwargs[var], mask, data)
if coerce_farray:
# coerce data and independent variable(s) that are 'array-like' (list,
# tuples, pandas Series) to float64/complex128.
data = coerce_arraylike(data)
for var in self.independent_vars:
kwargs[var] = coerce_arraylike(kwargs[var])
if fit_kws is None:
fit_kws = {}
output = ModelResult(self, params, method=method, iter_cb=iter_cb,
scale_covar=scale_covar, fcn_kws=kwargs,
nan_policy=self.nan_policy, calc_covar=calc_covar,
max_nfev=max_nfev, **fit_kws)
output.fit(data=data, weights=weights)
output.components = self.components
return output
def __add__(self, other):
"""+"""
return CompositeModel(self, other, operator.add)
def __sub__(self, other):
"""-"""
return CompositeModel(self, other, operator.sub)
def __mul__(self, other):
"""*"""
return CompositeModel(self, other, operator.mul)
def __truediv__(self, other):
"""/"""
return CompositeModel(self, other, operator.truediv)
class CompositeModel(Model):
"""Combine two models (`left` and `right`) with binary operator (`op`).
Normally, one does not have to explicitly create a `CompositeModel`,
but can use normal Python operators ``+``, ``-``, ``*``, and ``/`` to
combine components as in::
>>> mod = Model(fcn1) + Model(fcn2) * Model(fcn3)
"""
_known_ops = {operator.add: '+', operator.sub: '-',
operator.mul: '*', operator.truediv: '/'}
def __init__(self, left, right, op, **kws):
"""
Parameters
----------
left : Model
Left-hand model.
right : Model
Right-hand model.
op : callable binary operator
Operator to combine `left` and `right` models.
**kws : optional
Additional keywords are passed to `Model` when creating this
new model.
Notes
-----
The two models can use different independent variables.
"""
if not isinstance(left, Model):
raise ValueError(f'CompositeModel: argument {left} is not a Model')
if not isinstance(right, Model):
raise ValueError(f'CompositeModel: argument {right} is not a Model')
if not callable(op):
raise ValueError(f'CompositeModel: operator {op} is not callable')
self.left = left
self.right = right
self.op = op
name_collisions = set(left.param_names) & set(right.param_names)
if len(name_collisions) > 0:
msg = ''
for collision in name_collisions:
msg += (f"\nTwo models have parameters named '{collision}'; "
"use distinct names.")
raise NameError(msg)
# the unique ``independent_vars`` of the left and right model are
# combined to ``independent_vars`` of the ``CompositeModel``
if 'independent_vars' not in kws:
ivars = self.left.independent_vars + self.right.independent_vars
kws['independent_vars'] = list(np.unique(ivars))
if 'nan_policy' not in kws:
kws['nan_policy'] = self.left.nan_policy
def _tmp(self, *args, **kws):
pass
Model.__init__(self, _tmp, **kws)
for side in (left, right):
prefix = side.prefix
for basename, hint in side.param_hints.items():
self.param_hints[f"{prefix}{basename}"] = hint
def _parse_params(self):
self._func_haskeywords = (self.left._func_haskeywords or
self.right._func_haskeywords)
self._func_allargs = (self.left._func_allargs +
self.right._func_allargs)
self.def_vals = deepcopy(self.right.def_vals)
self.def_vals.update(self.left.def_vals)
self.opts = deepcopy(self.right.opts)
self.opts.update(self.left.opts)
def _reprstring(self, long=False):
return (f"({self.left._reprstring(long=long)} "
f"{self._known_ops.get(self.op, self.op)} "
f"{self.right._reprstring(long=long)})")
def eval(self, params=None, **kwargs):
"""Evaluate model function for composite model."""
return self.op(self.left.eval(params=params, **kwargs),
self.right.eval(params=params, **kwargs))
def eval_components(self, **kwargs):
"""Return dictionary of name, results for each component."""
out = dict(self.left.eval_components(**kwargs))
out.update(self.right.eval_components(**kwargs))
return out
def post_fit(self, fitresult):
"""function that is called just after fit, can be overloaded by
subclasses to add non-fitting 'calculated parameters'
"""
self.left.post_fit(fitresult)
self.right.post_fit(fitresult)
@property
def param_names(self):
"""Return parameter names for composite model."""
return self.left.param_names + self.right.param_names
@property
def components(self):
"""Return components for composite model."""
return self.left.components + self.right.components
def _get_state(self):
return (self.left._get_state(),
self.right._get_state(), self.op.__name__)
def _set_state(self, state, funcdefs=None):
return _buildmodel(state, funcdefs=funcdefs)
def _make_all_args(self, params=None, **kwargs):
"""Generate **all** function arguments for all functions."""
out = self.right._make_all_args(params=params, **kwargs)
out.update(self.left._make_all_args(params=params, **kwargs))
return out
def save_model(model, fname):
"""Save a Model to a file.
Parameters
----------
model : Model
Model to be saved.
fname : str
Name of file for saved Model.
"""
with open(fname, 'w') as fout:
model.dump(fout)
def load_model(fname, funcdefs=None):
"""Load a saved Model from a file.
Parameters
----------
fname : str
Name of file containing saved Model.
funcdefs : dict, optional
Dictionary of custom function names and definitions.
Returns
-------
Model
Model object loaded from file.
"""
m = Model(lambda x: x)
with open(fname) as fh:
model = m.load(fh, funcdefs=funcdefs)
return model
def _buildmodel(state, funcdefs=None):
"""Build Model from saved state.
Intended for internal use only.
"""
if len(state) != 3:
raise ValueError("Cannot restore Model")
known_funcs = {}
for fname in lineshapes.functions:
fcn = getattr(lineshapes, fname, None)
if callable(fcn):
known_funcs[fname] = fcn
if funcdefs is not None:
known_funcs.update(funcdefs)
left, right, op = state
if op is None and right is None:
(fname, fcndef, name, prefix, ivars, pnames,
phints, nan_policy, opts) = left
if not callable(fcndef) and fname in known_funcs:
fcndef = known_funcs[fname]
if fcndef is None:
raise ValueError("Cannot restore Model: model function not found")
if fname == '_eval' and isinstance(fcndef, str):
from .models import ExpressionModel
model = ExpressionModel(fcndef, name=name,
independent_vars=ivars,
param_names=pnames,
nan_policy=nan_policy, **opts)
else:
model = Model(fcndef, name=name, prefix=prefix,
independent_vars=ivars, param_names=pnames,
nan_policy=nan_policy, **opts)
for name, hint in phints.items():
model.set_param_hint(name, **hint)
return model
else:
lmodel = _buildmodel(left, funcdefs=funcdefs)
rmodel = _buildmodel(right, funcdefs=funcdefs)
return CompositeModel(lmodel, rmodel, getattr(operator, op))
def save_modelresult(modelresult, fname):
"""Save a ModelResult to a file.
Parameters
----------
modelresult : ModelResult
ModelResult to be saved.
fname : str
Name of file for saved ModelResult.
"""
with open(fname, 'w') as fout:
modelresult.dump(fout)
def load_modelresult(fname, funcdefs=None):
"""Load a saved ModelResult from a file.
Parameters
----------
fname : str
Name of file containing saved ModelResult.
funcdefs : dict, optional
Dictionary of custom function names and definitions.
Returns
-------
ModelResult
ModelResult object loaded from file.
"""
params = Parameters()
modres = ModelResult(Model(lambda x: x, None), params)
with open(fname) as fh:
mresult = modres.load(fh, funcdefs=funcdefs)
return mresult
class ModelResult(Minimizer):
"""Result from the Model fit.
This has many attributes and methods for viewing and working with the
results of a fit using Model. It inherits from Minimizer, so that it
can be used to modify and re-run the fit for the Model.
"""
def __init__(self, model, params, data=None, weights=None,
method='leastsq', fcn_args=None, fcn_kws=None,
iter_cb=None, scale_covar=True, nan_policy='raise',
calc_covar=True, max_nfev=None, **fit_kws):
"""
Parameters
----------
model : Model
Model to use.
params : Parameters
Parameters with initial values for model.
data : array_like, optional
Data to be modeled.
weights : array_like, optional
Weights to multiply ``(data-model)`` for fit residual.
method : str, optional
Name of minimization method to use (default is `'leastsq'`).
fcn_args : sequence, optional
Positional arguments to send to model function.
fcn_dict : dict, optional
Keyword arguments to send to model function.
iter_cb : callable, optional
Function to call on each iteration of fit.
scale_covar : bool, optional
Whether to scale covariance matrix for uncertainty evaluation.
nan_policy : {'raise', 'propagate', 'omit'}, optional
What to do when encountering NaNs when fitting Model.
calc_covar : bool, optional
Whether to calculate the covariance matrix (default is True)
for solvers other than `'leastsq'` and `'least_squares'`.
Requires the ``numdifftools`` package to be installed.
max_nfev : int or None, optional
Maximum number of function evaluations (default is None). The
default value depends on the fitting method.
**fit_kws : optional
Keyword arguments to send to minimization routine.
"""
self.model = model
self.data = data
self.weights = weights
self.method = method
self.ci_out = None
self.user_options = None
self.init_params = deepcopy(params)
Minimizer.__init__(self, model._residual, params,
fcn_args=fcn_args, fcn_kws=fcn_kws,
iter_cb=iter_cb, nan_policy=nan_policy,
scale_covar=scale_covar, calc_covar=calc_covar,
max_nfev=max_nfev, **fit_kws)
def fit(self, data=None, params=None, weights=None, method=None,
nan_policy=None, **kwargs):
"""Re-perform fit for a Model, given data and params.
Parameters
----------
data : array_like, optional
Data to be modeled.
params : Parameters, optional
Parameters with initial values for model.
weights : array_like, optional
Weights to multiply ``(data-model)`` for fit residual.
method : str, optional
Name of minimization method to use (default is `'leastsq'`).
nan_policy : {'raise', 'propagate', 'omit'}, optional
What to do when encountering NaNs when fitting Model.
**kwargs : optional
Keyword arguments to send to minimization routine.
"""
if data is not None:
self.data = data
if params is not None:
self.init_params = params
if weights is not None:
self.weights = weights
if method is not None:
self.method = method
if nan_policy is not None:
self.nan_policy = nan_policy
self.ci_out = None
self.userargs = (self.data, self.weights)
self.userkws.update(kwargs)
self.init_fit = self.model.eval(params=self.params, **self.userkws)
_ret = self.minimize(method=self.method)
self.model.post_fit(_ret)
_ret.params.create_uvars(covar=_ret.covar)
for attr in dir(_ret):
if not attr.startswith('_'):
try:
setattr(self, attr, getattr(_ret, attr))
except AttributeError:
pass
if self.data is not None and len(self.data) > 1:
dat = coerce_arraylike(self.data)
sstot = ((dat - dat.mean())**2).sum()
if isinstance(self.residual, np.ndarray) and len(self.residual) > 1:
self.rsquared = 1.0 - (self.residual**2).sum()/max(tiny, sstot)
self.init_values = self.model._make_all_args(self.init_params)
self.best_values = self.model._make_all_args(_ret.params)
self.best_fit = self.model.eval(params=_ret.params, **self.userkws)
def eval(self, params=None, **kwargs):
"""Evaluate model function.
Parameters
----------
params : Parameters, optional
Parameters to use.
**kwargs : optional
Options to send to Model.eval().
Returns
-------
numpy.ndarray, float, int, or complex
Array or value for the evaluated model.
"""
userkws = self.userkws.copy()
userkws.update(kwargs)
if params is None:
params = self.params
return self.model.eval(params=params, **userkws)
def eval_components(self, params=None, **kwargs):
"""Evaluate each component of a composite model function.
Parameters
----------
params : Parameters, optional
Parameters, defaults to ModelResult.params.
**kwargs : optional
Keyword arguments to pass to model function.
Returns
-------
dict
Keys are prefixes of component models, and values are the
estimated model value for each component of the model.
"""
userkws = self.userkws.copy()
userkws.update(kwargs)
if params is None:
params = self.params
return self.model.eval_components(params=params, **userkws)
def eval_uncertainty(self, params=None, sigma=1, **kwargs):
"""Evaluate the uncertainty of the *model function*.
This can be used to give confidence bands for the model from the
uncertainties in the best-fit parameters.
Parameters
----------
params : Parameters, optional
Parameters, defaults to ModelResult.params.
sigma : float, optional
Confidence level, i.e. how many sigma (default is 1).
**kwargs : optional
Values of options, independent variables, etcetera.
Returns
-------
numpy.ndarray
Uncertainty at each value of the model.
Notes
-----
1. This is based on the excellent and clear example from
https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#confidence-and-prediction-intervals,
which references the original work of:
J. Wolberg, Data Analysis Using the Method of Least Squares, 2006, Springer
2. The value of sigma is number of `sigma` values, and is converted
to a probability. Values of 1, 2, or 3 give probabilities of
0.6827, 0.9545, and 0.9973, respectively. If the sigma value is
< 1, it is interpreted as the probability itself. That is,
``sigma=1`` and ``sigma=0.6827`` will give the same results,
within precision errors.
3. Also sets attributes of `dely` for the uncertainty of the model
(which will be the same as the array returned by this method) and
`dely_comps`, a dictionary of `dely` for each component.
Examples
--------
>>> out = model.fit(data, params, x=x)
>>> dely = out.eval_uncertainty(x=x)
>>> plt.plot(x, data)
>>> plt.plot(x, out.best_fit)
>>> plt.fill_between(x, out.best_fit-dely,
... out.best_fit+dely, color='#888888')
"""
userkws = self.userkws.copy()
userkws.update(kwargs)
if params is None:
params = self.params
nvarys = self.nvarys
# ensure fjac and df2 are correct size if independent var updated by kwargs
feval = self.model.eval(params, **userkws)
ndata = len(feval.view('float64')) # allows feval to be complex
covar = self.covar
if any(p.stderr is None for p in params.values()):
return np.zeros(ndata)
# '0' would be an invalid prefix, here signifying 'Full'
fjac = {'0': np.zeros((nvarys, ndata), dtype='float64')}
df2 = {'0': np.zeros(ndata, dtype='float64')}
for comp in self.components:
label = comp.prefix if len(comp.prefix) > 1 else comp._name
fjac[label] = np.zeros((nvarys, ndata), dtype='float64')
df2[label] = np.zeros(ndata, dtype='float64')
# find derivative by hand!
pars = params.copy()
for i in range(nvarys):
pname = self.var_names[i]
val0 = pars[pname].value
dval = pars[pname].stderr/3.0
pars[pname].value = val0 + dval
res1 = {'0': self.model.eval(pars, **userkws)}
res1.update(self.model.eval_components(params=pars, **userkws))
pars[pname].value = val0 - dval
res2 = {'0': self.model.eval(pars, **userkws)}
res2.update(self.model.eval_components(params=pars, **userkws))
pars[pname].value = val0
for key in fjac:
fjac[key][i] = (res1[key].view('float64')
- res2[key].view('float64')) / (2*dval)
for i in range(nvarys):
for j in range(nvarys):
for key in fjac:
df2[key] += fjac[key][i] * fjac[key][j] * covar[i, j]
if sigma < 1.0:
prob = sigma
else:
prob = erf(sigma/np.sqrt(2))
scale = t.ppf((prob+1)/2.0, self.ndata-nvarys)
# for complex data, convert back to real/imag pairs
if feval.dtype in ('complex64', 'complex128'):
for key in fjac:
df2[key] = df2[key][0::2] + 1j * df2[key][1::2]
self.dely = scale * np.sqrt(df2.pop('0'))
self.dely_comps = {}
for key in df2:
self.dely_comps[key] = scale * np.sqrt(df2[key])
return self.dely
def conf_interval(self, **kwargs):
"""Calculate the confidence intervals for the variable parameters.
Confidence intervals are calculated using the
:func:`confidence.conf_interval` function and keyword arguments
(`**kwargs`) are passed to that function. The result is stored in
the :attr:`ci_out` attribute so that it can be accessed without
recalculating them.
"""
self.ci_out = conf_interval(self, self, **kwargs)
return self.ci_out
def ci_report(self, with_offset=True, ndigits=5, **kwargs):
"""Return a formatted text report of the confidence intervals.
Parameters
----------
with_offset : bool, optional
Whether to subtract best value from all other values (default
is True).
ndigits : int, optional
Number of significant digits to show (default is 5).
**kwargs : optional
Keyword arguments that are passed to the `conf_interval`
function.
Returns
-------
str
Text of formatted report on confidence intervals.
"""
return ci_report(self.conf_interval(**kwargs),
with_offset=with_offset, ndigits=ndigits)
def fit_report(self, modelpars=None, show_correl=True,
min_correl=0.1, sort_pars=False, correl_mode='list'):
"""Return a printable fit report.
The report contains fit statistics and best-fit values with
uncertainties and correlations.
Parameters
----------
modelpars : Parameters, optional
Known Model Parameters.
show_correl : bool, optional
Whether to show list of sorted correlations (default is True).
min_correl : float, optional
Smallest correlation in absolute value to show (default is 0.1).
sort_pars : callable, optional
Whether to show parameter names sorted in alphanumerical order
(default is False). If False, then the parameters will be
listed in the order as they were added to the Parameters
dictionary. If callable, then this (one argument) function is
used to extract a comparison key from each list element.
correl_mode : {'list', table'} str, optional
Mode for how to show correlations. Can be either 'list' (default)
to show a sorted (if ``sort_pars`` is True) list of correlation
values, or 'table' to show a complete, formatted table of
correlations.
Returns
-------
str
Multi-line text of fit report.
"""
report = fit_report(self, modelpars=modelpars, show_correl=show_correl,
min_correl=min_correl, sort_pars=sort_pars,
correl_mode=correl_mode)
modname = self.model._reprstring(long=True)
return f'[[Model]]\n {modname}\n{report}'
def _repr_html_(self, show_correl=True, min_correl=0.1):
"""Return a HTML representation of parameters data."""
report = fitreport_html_table(self, show_correl=show_correl,
min_correl=min_correl)
modname = self.model._reprstring(long=True)
return f"<h2>Fit Result</h2> <p>Model: {modname}</p> {report}"
def summary(self):
"""Return a dictionary with statistics and attributes of a ModelResult.
Returns
-------
dict
Dictionary of statistics and many attributes from a ModelResult.
Notes
------
1. values for data arrays are not included.
2. The result summary dictionary will include the following entries:
``model``, ``method``, ``ndata``, ``nvarys``, ``nfree``, ``chisqr``,
``redchi``, ``aic``, ``bic``, ``rsquared``, ``nfev``, ``max_nfev``,
``aborted``, ``errorbars``, ``success``, ``message``,
``lmdif_message``, ``ier``, ``nan_policy``, ``scale_covar``,
``calc_covar``, ``ci_out``, ``col_deriv``, ``flatchain``,
``call_kws``, ``var_names``, ``user_options``, ``kws``,
``init_values``, ``best_values``, and ``params``.
where 'params' is a list of parameter "states": tuples with entries of
``(name, value, vary, expr, min, max, brute_step, stderr, correl,
init_value, user_data)``.
3. The result will include only plain Python objects, and so should be
easily serializable with JSON or similar tools.
"""
summary = {'model': self.model._reprstring(long=True),
'method': self.method}
for attr in ('ndata', 'nvarys', 'nfree', 'chisqr', 'redchi', 'aic',
'bic', 'rsquared', 'nfev', 'max_nfev', 'aborted',
'errorbars', 'success', 'message', 'lmdif_message', 'ier',
'nan_policy', 'scale_covar', 'calc_covar', 'ci_out',
'col_deriv', 'flatchain', 'call_kws', 'var_names',
'user_options', 'kws', 'init_values', 'best_values'):
val = getattr(self, attr, None)
if isinstance(val, np.float64):
val = float(val)
elif isinstance(val, (np.int32, np.int64)):
val = int(val)
elif isinstance(val, np.bool_):
val = bool(val)
elif isinstance(val, bytes):
val = str(val, encoding='UTF-8')
summary[attr] = val
summary['params'] = [par.__getstate__() for par in self.params.values()]
return summary
def dumps(self, **kws):
"""Represent ModelResult as a JSON string.
Parameters
----------
**kws : optional
Keyword arguments that are passed to `json.dumps`.
Returns
-------
str
JSON string representation of ModelResult.
See Also
--------
loads, json.dumps
"""
out = {'__class__': 'lmfit.ModelResult', '__version__': '1',
'model': encode4js(self.model._get_state())}
pasteval = self.params._asteval
out['params'] = [p.__getstate__() for p in self.params.values()]
out['unique_symbols'] = {key: encode4js(pasteval.symtable[key])
for key in pasteval.user_defined_symbols()}
for attr in ('aborted', 'aic', 'best_values', 'bic', 'chisqr',
'ci_out', 'col_deriv', 'covar', 'errorbars', 'flatchain',
'ier', 'init_values', 'lmdif_message', 'message',
'method', 'nan_policy', 'ndata', 'nfev', 'nfree',
'nvarys', 'redchi', 'residual', 'rsquared', 'scale_covar',
'calc_covar', 'success', 'userargs', 'userkws', 'values',
'var_names', 'weights', 'user_options'):
try:
val = getattr(self, attr)
except AttributeError:
continue
if isinstance(val, np.bool_):
val = bool(val)
out[attr] = encode4js(val)
val = out.get('message', '')
if isinstance(val, bytes):
out['message'] = str(val, encoding='ASCII')
return json.dumps(out, **kws)
def dump(self, fp, **kws):
"""Dump serialization of ModelResult to a file.
Parameters
----------
fp : file-like object
An open and `.write()`-supporting file-like object.
**kws : optional
Keyword arguments that are passed to `json.dumps`.
Returns
-------
int
Return value from `fp.write()`: the number of characters
written.
See Also
--------
dumps, load, json.dump
"""
return fp.write(self.dumps(**kws))
def loads(self, s, funcdefs=None, **kws):
"""Load ModelResult from a JSON string.
Parameters
----------
s : str
String representation of ModelResult, as from `dumps`.
funcdefs : dict, optional
Dictionary of custom function names and definitions.
**kws : optional
Keyword arguments that are passed to `json.dumps`.
Returns
-------
ModelResult
ModelResult instance from JSON string representation.
See Also
--------
load, dumps, json.dumps
"""
modres = json.loads(s, **kws)
if 'modelresult' not in modres['__class__'].lower():
raise AttributeError('ModelResult.loads() needs saved ModelResult')
modres = decode4js(modres)
if 'model' not in modres or 'params' not in modres:
raise AttributeError('ModelResult.loads() needs valid ModelResult')
# model
self.model = _buildmodel(decode4js(modres['model']), funcdefs=funcdefs)
# params
for target in ('params', 'init_params'):
state = {'unique_symbols': modres['unique_symbols'], 'params': []}
for parstate in modres['params']:
_par = Parameter(name='')
_par.__setstate__(parstate)
state['params'].append(_par)
_params = Parameters()
_params.__setstate__(state)
setattr(self, target, _params)
for attr in ('aborted', 'aic', 'best_fit', 'best_values', 'bic',
'chisqr', 'ci_out', 'col_deriv', 'covar', 'data',
'errorbars', 'fjac', 'flatchain', 'ier', 'init_fit',
'init_values', 'kws', 'lmdif_message', 'message',
'method', 'nan_policy', 'ndata', 'nfev', 'nfree',
'nvarys', 'redchi', 'residual', 'rsquared', 'scale_covar',
'calc_covar', 'success', 'userargs', 'userkws',
'var_names', 'weights', 'user_options'):
setattr(self, attr, decode4js(modres.get(attr, None)))
self.best_fit = self.model.eval(self.params, **self.userkws)
if len(self.userargs) == 2:
self.data = self.userargs[0]
self.weights = self.userargs[1]
for parname, val in self.init_values.items():
par = self.init_params.get(parname, None)
if par is not None:
par.correl = par.stderr = None
par.value = par.init_value = self.init_values[parname]
self.init_fit = self.model.eval(self.init_params, **self.userkws)
self.result = MinimizerResult()
self.result.params = self.params
if self.errorbars and self.covar is not None:
self.uvars = self.result.params.create_uvars(covar=self.covar)
self.init_vals = list(self.init_values.items())
return self
def load(self, fp, funcdefs=None, **kws):
"""Load JSON representation of ModelResult from a file-like object.
Parameters
----------
fp : file-like object
An open and `.read()`-supporting file-like object.
funcdefs : dict, optional
Dictionary of function definitions to use to construct Model.
**kws : optional
Keyword arguments that are passed to `loads`.
Returns
-------
ModelResult
ModelResult created from `fp`.
See Also
--------
dump, loads, json.load
"""
return self.loads(fp.read(), funcdefs=funcdefs, **kws)
@_ensureMatplotlib
def plot_fit(self, ax=None, datafmt='o', fitfmt='-', initfmt='--',
xlabel=None, ylabel=None, yerr=None, numpoints=None,
data_kws=None, fit_kws=None, init_kws=None, ax_kws=None,
show_init=False, parse_complex='abs', title=None):
"""Plot the fit results using matplotlib, if available.
The plot will include the data points, the initial fit curve
(optional, with ``show_init=True``), and the best-fit curve. If
the fit model included weights or if `yerr` is specified,
errorbars will also be plotted.
Parameters
----------
ax : matplotlib.axes.Axes, optional
The axes to plot on. The default in None, which means use the
current pyplot axis or create one if there is none.
datafmt : str, optional
Matplotlib format string for data points.
fitfmt : str, optional
Matplotlib format string for fitted curve.
initfmt : str, optional
Matplotlib format string for initial conditions for the fit.
xlabel : str, optional
Matplotlib format string for labeling the x-axis.
ylabel : str, optional
Matplotlib format string for labeling the y-axis.
yerr : numpy.ndarray, optional
Array of uncertainties for data array.
numpoints : int, optional
If provided, the final and initial fit curves are evaluated
not only at data points, but refined to contain `numpoints`
points in total.
data_kws : dict, optional
Keyword arguments passed to the plot function for data points.
fit_kws : dict, optional
Keyword arguments passed to the plot function for fitted curve.
init_kws : dict, optional
Keyword arguments passed to the plot function for the initial
conditions of the fit.
ax_kws : dict, optional
Keyword arguments for a new axis, if a new one is created.
show_init : bool, optional
Whether to show the initial conditions for the fit (default is
False).
parse_complex : {'abs', 'real', 'imag', 'angle'}, optional
How to reduce complex data for plotting. Options are one of:
`'abs'` (default), `'real'`, `'imag'`, or `'angle'`, which
correspond to the NumPy functions with the same name.
title : str, optional
Matplotlib format string for figure title.
Returns
-------
matplotlib.axes.Axes
See Also
--------
ModelResult.plot_residuals : Plot the fit residuals using matplotlib.
ModelResult.plot : Plot the fit results and residuals using matplotlib.
Notes
-----
For details about plot format strings and keyword arguments see
documentation of `matplotlib.axes.Axes.plot`.
If `yerr` is specified or if the fit model included weights, then
`matplotlib.axes.Axes.errorbar` is used to plot the data. If
`yerr` is not specified and the fit includes weights, `yerr` set
to ``1/self.weights``.
If model returns complex data, `yerr` is treated the same way that
weights are in this case.
If `ax` is None then `matplotlib.pyplot.gca(**ax_kws)` is called.
"""
from matplotlib import pyplot as plt
if data_kws is None:
data_kws = {}
if fit_kws is None:
fit_kws = {}
if init_kws is None:
init_kws = {}
if ax_kws is None:
ax_kws = {}
# The function reduce_complex will convert complex vectors into real vectors
reduce_complex = get_reducer(parse_complex)
if len(self.model.independent_vars) == 1:
independent_var = self.model.independent_vars[0]
else:
print('Fit can only be plotted if the model function has one '
'independent variable.')
return False
if not isinstance(ax, plt.Axes):
ax = plt.axes(**ax_kws)
x_array = self.userkws[independent_var]
# make a dense array for x-axis if data is not dense
if numpoints is not None and len(self.data) < numpoints:
x_array_dense = np.linspace(min(x_array), max(x_array), numpoints)
else:
x_array_dense = x_array
if show_init:
y_eval_init = self.model.eval(self.init_params,
**{independent_var: x_array_dense})
if isinstance(self.model, (lmfit.models.ConstantModel,
lmfit.models.ComplexConstantModel)):
y_eval_init *= np.ones(x_array_dense.size)
ax.plot(
x_array_dense, reduce_complex(y_eval_init), initfmt,
label='initial fit', **init_kws)
if yerr is None and self.weights is not None:
yerr = 1.0/self.weights
if yerr is not None:
ax.errorbar(x_array, reduce_complex(self.data),
yerr=propagate_err(self.data, yerr, parse_complex),
fmt=datafmt, label='data', **data_kws)
else:
ax.plot(x_array, reduce_complex(self.data),
datafmt, label='data', **data_kws)
y_eval = self.model.eval(self.params, **{independent_var: x_array_dense})
if isinstance(self.model, (lmfit.models.ConstantModel,
lmfit.models.ComplexConstantModel)):
y_eval *= np.ones(x_array_dense.size)
ax.plot(x_array_dense, reduce_complex(y_eval), fitfmt, label='best fit',
**fit_kws)
if title:
ax.set_title(title)
elif ax.get_title() == '':
ax.set_title(self.model.name)
if xlabel is None:
ax.set_xlabel(independent_var)
else:
ax.set_xlabel(xlabel)
if ylabel is None:
ax.set_ylabel('y')
else:
ax.set_ylabel(ylabel)
ax.legend()
return ax
@_ensureMatplotlib
def plot_residuals(self, ax=None, datafmt='o', yerr=None, data_kws=None,
fit_kws=None, ax_kws=None, parse_complex='abs',
title=None):
"""Plot the fit residuals using matplotlib, if available.
If `yerr` is supplied or if the model included weights, errorbars
will also be plotted.
Parameters
----------
ax : matplotlib.axes.Axes, optional
The axes to plot on. The default in None, which means use the
current pyplot axis or create one if there is none.
datafmt : str, optional
Matplotlib format string for data points.
yerr : numpy.ndarray, optional
Array of uncertainties for data array.
data_kws : dict, optional
Keyword arguments passed to the plot function for data points.
fit_kws : dict, optional
Keyword arguments passed to the plot function for fitted curve.
ax_kws : dict, optional
Keyword arguments for a new axis, if a new one is created.
parse_complex : {'abs', 'real', 'imag', 'angle'}, optional
How to reduce complex data for plotting. Options are one of:
`'abs'` (default), `'real'`, `'imag'`, or `'angle'`, which
correspond to the NumPy functions with the same name.
title : str, optional
Matplotlib format string for figure title.
Returns
-------
matplotlib.axes.Axes
See Also
--------
ModelResult.plot_fit : Plot the fit results using matplotlib.
ModelResult.plot : Plot the fit results and residuals using matplotlib.
Notes
-----
For details about plot format strings and keyword arguments see
documentation of `matplotlib.axes.Axes.plot`.
If `yerr` is specified or if the fit model included weights, then
`matplotlib.axes.Axes.errorbar` is used to plot the data. If
`yerr` is not specified and the fit includes weights, `yerr` set
to ``1/self.weights``.
If model returns complex data, `yerr` is treated the same way that
weights are in this case.
If `ax` is None then `matplotlib.pyplot.gca(**ax_kws)` is called.
"""
from matplotlib import pyplot as plt
if data_kws is None:
data_kws = {}
if fit_kws is None:
fit_kws = {}
if ax_kws is None:
ax_kws = {}
# The function reduce_complex will convert complex vectors into real vectors
reduce_complex = get_reducer(parse_complex)
if len(self.model.independent_vars) == 1:
independent_var = self.model.independent_vars[0]
else:
print('Fit can only be plotted if the model function has one '
'independent variable.')
return False
if not isinstance(ax, plt.Axes):
ax = plt.axes(**ax_kws)
x_array = self.userkws[independent_var]
ax.axhline(0, **fit_kws, color='k')
y_eval = self.model.eval(self.params, **{independent_var: x_array})
if isinstance(self.model, (lmfit.models.ConstantModel,
lmfit.models.ComplexConstantModel)):
y_eval *= np.ones(x_array.size)
if yerr is None and self.weights is not None:
yerr = 1.0/self.weights
residuals = reduce_complex(self.eval()) - reduce_complex(self.data)
if yerr is not None:
ax.errorbar(x_array, residuals,
yerr=propagate_err(self.data, yerr, parse_complex),
fmt=datafmt, **data_kws)
else:
ax.plot(x_array, residuals, datafmt, **data_kws)
if title:
ax.set_title(title)
elif ax.get_title() == '':
ax.set_title(self.model.name)
ax.set_ylabel('residuals')
return ax
@_ensureMatplotlib
def plot(self, datafmt='o', fitfmt='-', initfmt='--', xlabel=None,
ylabel=None, yerr=None, numpoints=None, fig=None, data_kws=None,
fit_kws=None, init_kws=None, ax_res_kws=None, ax_fit_kws=None,
fig_kws=None, show_init=False, parse_complex='abs', title=None):
"""Plot the fit results and residuals using matplotlib.
The method will produce a matplotlib figure (if package available)
with both results of the fit and the residuals plotted. If the fit
model included weights, errorbars will also be plotted. To show
the initial conditions for the fit, pass the argument
``show_init=True``.
Parameters
----------
datafmt : str, optional
Matplotlib format string for data points.
fitfmt : str, optional
Matplotlib format string for fitted curve.
initfmt : str, optional
Matplotlib format string for initial conditions for the fit.
xlabel : str, optional
Matplotlib format string for labeling the x-axis.
ylabel : str, optional
Matplotlib format string for labeling the y-axis.
yerr : numpy.ndarray, optional
Array of uncertainties for data array.
numpoints : int, optional
If provided, the final and initial fit curves are evaluated
not only at data points, but refined to contain `numpoints`
points in total.
fig : matplotlib.figure.Figure, optional
The figure to plot on. The default is None, which means use
the current pyplot figure or create one if there is none.
data_kws : dict, optional
Keyword arguments passed to the plot function for data points.
fit_kws : dict, optional
Keyword arguments passed to the plot function for fitted curve.
init_kws : dict, optional
Keyword arguments passed to the plot function for the initial
conditions of the fit.
ax_res_kws : dict, optional
Keyword arguments for the axes for the residuals plot.
ax_fit_kws : dict, optional
Keyword arguments for the axes for the fit plot.
fig_kws : dict, optional
Keyword arguments for a new figure, if a new one is created.
show_init : bool, optional
Whether to show the initial conditions for the fit (default is
False).
parse_complex : {'abs', 'real', 'imag', 'angle'}, optional
How to reduce complex data for plotting. Options are one of:
`'abs'` (default), `'real'`, `'imag'`, or `'angle'`, which
correspond to the NumPy functions with the same name.
title : str, optional
Matplotlib format string for figure title.
Returns
-------
matplotlib.figure.Figure
See Also
--------
ModelResult.plot_fit : Plot the fit results using matplotlib.
ModelResult.plot_residuals : Plot the fit residuals using matplotlib.
Notes
-----
The method combines `ModelResult.plot_fit` and
`ModelResult.plot_residuals`.
If `yerr` is specified or if the fit model included weights, then
`matplotlib.axes.Axes.errorbar` is used to plot the data. If
`yerr` is not specified and the fit includes weights, `yerr` set
to ``1/self.weights``.
If model returns complex data, `yerr` is treated the same way that
weights are in this case.
If `fig` is None then `matplotlib.pyplot.figure(**fig_kws)` is
called, otherwise `fig_kws` is ignored.
"""
from matplotlib import pyplot as plt
if data_kws is None:
data_kws = {}
if fit_kws is None:
fit_kws = {}
if init_kws is None:
init_kws = {}
if ax_res_kws is None:
ax_res_kws = {}
if ax_fit_kws is None:
ax_fit_kws = {}
# make a square figure with side equal to the default figure's x-size
figxsize = plt.rcParams['figure.figsize'][0]
fig_kws_ = dict(figsize=(figxsize, figxsize))
if fig_kws is not None:
fig_kws_.update(fig_kws)
if len(self.model.independent_vars) != 1:
print('Fit can only be plotted if the model function has one '
'independent variable.')
return False
if not isinstance(fig, plt.Figure):
fig = plt.figure(**fig_kws_)
gs = plt.GridSpec(nrows=2, ncols=1, height_ratios=[1, 4])
ax_res = fig.add_subplot(gs[0], **ax_res_kws)
ax_fit = fig.add_subplot(gs[1], sharex=ax_res, **ax_fit_kws)
self.plot_fit(ax=ax_fit, datafmt=datafmt, fitfmt=fitfmt, yerr=yerr,
initfmt=initfmt, xlabel=xlabel, ylabel=ylabel,
numpoints=numpoints, data_kws=data_kws,
fit_kws=fit_kws, init_kws=init_kws, ax_kws=ax_fit_kws,
show_init=show_init, parse_complex=parse_complex,
title=title)
self.plot_residuals(ax=ax_res, datafmt=datafmt, yerr=yerr,
data_kws=data_kws, fit_kws=fit_kws,
ax_kws=ax_res_kws, parse_complex=parse_complex,
title=title)
plt.setp(ax_res.get_xticklabels(), visible=False)
ax_fit.set_title('')
return fig
