File: fitters.py

package info (click to toggle)
glueviz 0.9.1%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 17,180 kB
  • ctags: 6,728
  • sloc: python: 37,111; makefile: 134; sh: 60
file content (382 lines) | stat: -rw-r--r-- 11,728 bytes parent folder | download
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]