1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382
|
"""
Glue's fitting classes are designed to be easily subclassed for performing
custom model fitting in Glue.
See the guide on :ref:`writing custom fit plugins <fit_plugins>` for
help with using custom fitting utilities in Glue.
"""
from __future__ import absolute_import, division, print_function
import numpy as np
from glue.core.simpleforms import IntOption, Option
__all__ = ['BaseFitter1D',
'PolynomialFitter',
'AstropyFitter1D',
'SimpleAstropyGaussianFitter',
'BasicGaussianFitter']
class BaseFitter1D(object):
"""
Base class for 1D fitters.
This abstract class must be overwritten.
"""
label = "Fitter"
"""A short label for the fit, used by the GUI"""
param_names = []
"""list of parameter names that support restrictions"""
def __init__(self, **params):
self._constraints = {}
for k, v in params.items():
if k in self.param_names:
self.set_constraint(k, value=v)
else:
setattr(self, k, v)
def plot(self, fit_result, axes, x):
"""
Plot the result of a fit.
:param fit_result: The output from fit
:param axes: The Matplotlib axes to add the fit to
:param x: The values of X at which to visualize the model
:returns: A list of matplotlib artists. **This is important:**
plots will not be properly cleared if this isn't provided
"""
y = self.predict(fit_result, x)
result = axes.plot(x, y, '#4daf4a',
lw=3, alpha=0.8,
scalex=False, scaley=False)
return result
def _sigma_to_weights(self, dy):
if dy is not None:
return 1. / np.asarray(dy) ** 2
@property
def options(self):
"""
A dictionary of the current setting of each model hyperparameter.
Hyperparameters are defined in subclasses by creating class-level
:mod:`Option <glue.core.simpleforms>` attributes. This attribute
dict maps ``{hyperparameter_name: current_value}``
"""
result = []
for typ in type(self).mro():
result.extend(k for k, v in typ.__dict__.items()
if isinstance(v, Option))
return dict((o, getattr(self, o)) for o in result)
def summarize(self, fit_result, x, y, dy=None):
"""
Return a textual summary of the fit.
:param fit_result: The return value from :meth:`fit`
:param x: The x values passed to :meth:`fit`
:returns: A description of the fit result
:rtype: str
"""
return str(fit_result)
@property
def constraints(self):
"""
A dict of the constraints on each parameter in :attr:`param_names`.
Each value is itself a dict with 3 items:
:key value: The default value
:key fixed: True / False, indicating whether the parameter is fixed
:key bounds: [min, max] or None, indicating lower/upper limits
"""
result = {}
for p in self.param_names:
result[p] = dict(value=None, fixed=False, limits=None)
result[p].update(self._constraints.get(p, {}))
return result
def set_constraint(self, parameter_name, value=None,
fixed=None, limits=None):
"""
Update a constraint.
:param parameter_name: name of the parameter to update
:type parameter_name: str
:param value: Set the default value (optional)
:param limits: Set the limits to[min, max] (optional)
:param fixed: Set whether the parameter is fixed (optional)
"""
c = self._constraints.setdefault(parameter_name, {})
if value is not None:
c['value'] = value
if fixed is not None:
c['fixed'] = fixed
if limits is not None:
c['limits'] = limits
def build_and_fit(self, x, y, dy=None):
"""
Method which builds the arguments to fit, and calls that method
"""
x = np.asarray(x).ravel()
y = np.asarray(y).ravel()
if dy is not None:
dy = np.asarray(dy).ravel()
return self.fit(x, y, dy=dy,
constraints=self.constraints,
**self.options)
def fit(self, x, y, dy, constraints, **options):
"""
Fit the model to data.
*This must be overriden by a subclass.*
:param x: The x values of the data
:type x: :class:`numpy.ndarray`
:param y: The y values of the data
:type y: :class:`numpy.ndarray`
:param dy: 1 sigma uncertainties on each datum (optional)
:type dy: :class:`numpy.ndarray`
:param constraints: The current value of :attr:`constraints`
:param options: kwargs for model hyperparameters.
:returns: An object representing the fit result.
"""
raise NotImplementedError()
def predict(self, fit_result, x):
"""
Evaulate the model at a set of locations.
**This must be overridden in a subclass.**
:param fit_result: The result from the fit method
:param x: Locations to evaluate model at
:type x: :class:`numpy.ndarray`
:returns: model(x)
:rtype: :class:`numpy.ndarray`
"""
raise NotImplementedError()
class AstropyFitter1D(BaseFitter1D):
"""
A base class for wrapping :mod:`astropy.modeling`.
Subclasses must override :attr:`model_cls` :attr:`fitting_cls`
to point to the desired Astropy :mod:`model <astropy.modeling>`
and :mod:`fitter <astropy.modeling.fitting>` classes.
In addition, they should override :attr:`label` with a better label,
and :meth:`parameter_guesses` to generate initial guesses
"""
model_cls = None
"""class describing the model"""
fitting_cls = None
"""class to fit the model"""
label = "Base Astropy Fitter"
"""UI Label"""
@property
def param_names(self):
return self.model_cls.param_names
def predict(self, fit_result, x):
model, _ = fit_result
return model(x)
def summarize(self, fit_result, x, y, dy=None):
model, fitter = fit_result
result = [_report_fitter(fitter), ""]
pnames = list(sorted(model.param_names))
maxlen = max(map(len, pnames))
result.extend("%s = %e" % (p.ljust(maxlen), getattr(model, p).value)
for p in pnames)
return "\n".join(result)
def fit(self, x, y, dy, constraints):
m, f = self._get_model_fitter(x, y, dy, constraints)
dy = self._sigma_to_weights(dy)
return f(m, x, y, weights=dy), f
def _get_model_fitter(self, x, y, dy, constraints):
if self.model_cls is None or self.fitting_cls is None:
raise NotImplementedError("Model or fitting class is unspecified.")
params = dict((k, v['value']) for k, v in constraints.items())
# update unset parameters with guesses from data
for k, v in self.parameter_guesses(x, y, dy).items():
if params[k] is not None or constraints[k]['fixed']:
continue
params[k] = v
m = self.model_cls(**params)
f = self.fitting_cls()
for param_name, constraint in constraints.items():
param = getattr(m, param_name)
if constraint['fixed']:
param.fixed = True
if constraint['limits']:
param.min, param.max = constraint['limits']
return m, f
def parameter_guesses(self, x, y, dy):
"""
Provide initial guesses for each model parameter.
**The base implementation does nothing, and should be overridden**
:param x: X - values of the data
:type x: :class:`numpy.ndarray`
:param y: Y - values of the data
:type y: :class:`numpy.ndarray`
:param dy: ncertainties on Y(assumed to be 1 sigma)
:type dy: :class:`numpy.ndarray`
:returns: A dict maping ``{parameter_name: value guess}`` for each
parameter
"""
return {}
def _gaussian_parameter_estimates(x, y, dy):
amplitude = np.percentile(y, 95)
y = np.maximum(y / y.sum(), 0)
mean = (x * y).sum()
stddev = np.sqrt((y * (x - mean) ** 2).sum())
return dict(mean=mean, stddev=stddev, amplitude=amplitude)
class BasicGaussianFitter(BaseFitter1D):
"""
Fallback Gaussian fitter, for astropy < 0.3.
If :mod:`astropy.modeling` is installed, this class is replaced by
:class:`SimpleAstropyGaussianFitter`
"""
label = "Gaussian"
def _errorfunc(self, params, x, y, dy):
yp = self.eval(x, *params)
result = (yp - y)
if dy is not None:
result /= dy
return result
@staticmethod
def eval(x, amplitude, mean, stddev):
return np.exp(-(x - mean) ** 2 / (2 * stddev ** 2)) * amplitude
@staticmethod
def fit_deriv(x, amplitude, mean, stddev):
"""
Gaussian1D model function derivatives.
"""
d_amplitude = np.exp(-0.5 / stddev ** 2 * (x - mean) ** 2)
d_mean = amplitude * d_amplitude * (x - mean) / stddev ** 2
d_stddev = amplitude * d_amplitude * (x - mean) ** 2 / stddev ** 3
return [d_amplitude, d_mean, d_stddev]
def fit(self, x, y, dy, constraints):
from scipy import optimize
init_values = _gaussian_parameter_estimates(x, y, dy)
init_values = [init_values[p] for p in ['amplitude', 'mean', 'stddev']]
farg = (x, y, dy)
dfunc = None
fitparams, status, dinfo, mess, ierr = optimize.leastsq(
self._errorfunc, init_values, args=farg, Dfun=dfunc,
full_output=True)
return fitparams
def predict(self, fit_result, x):
return self.eval(x, *fit_result)
def summarize(self, fit_result, x, y, dy=None):
return ("amplitude = %e\n"
"mean = %e\n"
"stddev = %e" % tuple(fit_result))
GaussianFitter = BasicGaussianFitter
try:
from astropy.modeling import models, fitting
class SimpleAstropyGaussianFitter(AstropyFitter1D):
"""
Guassian fitter using astropy.modeling.
"""
model_cls = models.Gaussian1D
try:
fitting_cls = fitting.LevMarLSQFitter
except AttributeError: # astropy v0.3
fitting_cls = fitting.NonLinearLSQFitter
label = "Gaussian"
parameter_guesses = staticmethod(_gaussian_parameter_estimates)
GaussianFitter = SimpleAstropyGaussianFitter
except ImportError:
pass
class PolynomialFitter(BaseFitter1D):
"""
A polynomial model.
The degree of the polynomial is specified by :attr:`degree`
"""
label = "Polynomial"
degree = IntOption(min=0, max=5, default=3, label="Polynomial Degree")
def fit(self, x, y, dy, constraints, degree=2):
"""
Fit a ``degree``-th order polynomial to the data.
"""
w = self._sigma_to_weights(dy)
return np.polyfit(x, y, degree, w=w)
def predict(self, fit_result, x):
return np.polyval(fit_result, x)
def summarize(self, fit_result, x, y, dy=None):
return "Coefficients:\n" + "\n".join("%e" % coeff
for coeff in fit_result.tolist())
def _report_fitter(fitter):
if "nfev" in fitter.fit_info:
return "Converged in %i iterations" % fitter.fit_info['nfev']
return 'Converged'
__FITTERS__ = [PolynomialFitter, GaussianFitter]
|