cdeil · Dec 9, 2013 · Dec 9, 2013 · Dec 9, 2013 · Dec 9, 2013 · Dec 9, 2013
diff --git a/docs/modeling/custom_model.py b/docs/modeling/custom_model.py
@@ -0,0 +1,28 @@
+__doctest_requires__ = {('.'): ['matplotlib', 'scipy']}
+import numpy as np
+import matplotlib.pyplot as plt
+from astropy.modeling.models import custom_model_1d
+from astropy.modeling.fitting import NonLinearLSQFitter
+
+# Define model
+@custom_model_1d
+def sum_of_gaussians(x, amplitude1=1., mean1=-1., sigma1=1.,
+                     amplitude2=1., mean2=1., sigma2=1.):
+    return (amplitude1 * np.exp(-0.5 * ((x - mean1) / sigma1) ** 2) +
+            amplitude2 * np.exp(-0.5 * ((x - mean2) / sigma2) ** 2))
+
+# Generate fake data
+np.random.seed(0)
+x = np.linspace(-5., 5., 200)
+m_ref = sum_of_gaussians(amplitude1=2., mean1=-0.5, sigma1=0.4,
+                         amplitude2=0.5, mean2=2., sigma2=1.0)
+y = m_ref(x) + np.random.normal(0., 0.1, x.shape)
+
+# Fit model to data
+m_init = sum_of_gaussians()
+fit = NonLinearLSQFitter()
+m = fit(m_init, x, y)
+
+# Plot the data and the best fit
+plt.plot(x, y, 'o', color='k')
+plt.plot(x, m(x), color='r', lw=2)
diff --git a/docs/modeling/index.rst b/docs/modeling/index.rst
@@ -177,7 +177,8 @@ Using `modeling`
    parameters
    models
    fitting
-   new
+   new_model
+   new_fitter
    algorithms
    design
 

diff --git a/docs/modeling/line_model.py b/docs/modeling/line_model.py
@@ -0,0 +1,28 @@
+import numpy as np
+from astropy.modeling import models, Parameter, format_input
+
+class LineModel(models.PolynomialModel):
+    slope = Parameter('slope')
+    intercept = Parameter('intercept')
+    linear = True
+
+    def __init__(self, slope, intercept, param_dim=1, **constraints):
+        super(LineModel, self).__init__(slope=slope, intercept=intercept,
+                                        param_dim=param_dim, **constraints)
+        self.domain = [-1, 1]
+        self.window = [-1, 1]
+        self._order = 2
+
+    @staticmethod
+    def eval(x, slope, intercept):
+        return slope * x + intercept
+
+    @staticmethod
+    def fit_deriv(x, slope, intercept):
+        d_slope = x
+        d_intercept = np.ones_like(x)
+        return [d_slope, d_intercept]
+
+    @format_input
+    def __call__(self, x):
+        return self.eval(x, *self.param_sets)
diff --git a/docs/modeling/new_fitter.py b/docs/modeling/new_fitter.py
@@ -0,0 +1,27 @@
+__doctest_requires__ = {('Chi2Fitter.*'): ['scipy']}
+from astropy.modeling.fitting import Fitter
+
+class Chi2Fitter(Fitter):
+    """Chi-square (a.k.a. least squares) fitter.
+
+    Wraps the `scipy.optimize.minimize` function.
+    """
+    supported_constraints = ['bounds', 'fixed']
+
+    def errorfunc(self, fitparams, *args):
+        """Compute chi-square fit statistic."""
+        model, x, y = args
+        self._fitter_to_model_params(model, fitparams)
+        return ((y - model(x)) ** 2).sum()
+
+    def __call__(self, model, x, y , maxiter=100):
+        """Minimize objective function."""
+        from scipy.optimize import minimize
+        self._validate_constraints(model)
+        model_copy = model.copy()
+        x0, _ = model_copy._model_to_fit_params()
+        fitparams = minimize(self.errorfunc, x0=x0,
+                            args=(model_copy, x, y),
+                            iter=maxiter)
+        self._fitter_to_model_params(model_copy, fitparams)
+        return model_copy
diff --git a/docs/modeling/new_fitter.rst b/docs/modeling/new_fitter.rst
@@ -0,0 +1,67 @@
+***************************
+Defining a new Fitter class
+***************************
+
+If you would like to use a fit statistic or minimizer different from the ones built into `~astropy.modeling`,
+you have to create a subclass of `~astropy.modeling.fitting.Fitter`. 
+
+Here is an example how to define a least squares fitter for nonlinear models
+by wrapping the `~scipy.optimize.minimize`  minimisation function:
+
+.. literalinclude:: new_fitter.py
+   :linenos:
+
+The user calls ``Chi2Fitter.__call__`` with a given ``model`` and data ``x`` and ``y``,
+which checks and reformats the inputs a bit and then calls `~scipy.optimize.minimize`,
+which repeatedly calls ``Chi2Fitter.errorfunc`` in each fit iteration. 
+
+Let's go through it line by line:
+
+* Line TODO: Create ``Chi2Fitter`` as a sub-class of `~astropy.modeling.fitting.Fitter`.
+* Line TODO: Declare the supported contraints of this fitter as a list of strings.
+  TODO: link to complete list and explanation of supported constraints types.
+
+The special method ``Chi2Fitter.errorfunc`` computes the chi-square fit statistic:
+
+* The ``fitparams`` array contains the parameter values in the current fit iteration.
+* The ``args`` parameter tuple can be used to pass other parameters needed to compute the
+  fit statistic, in this case ``args[0] = model``, ``args[1] = x`` and ``args[2] = y``.
+* Line TODO: The `~astropy.modeling.fitting.Fitter._fitter_to_model_params` helper function is used
+  to set the model parameters.
+* Line TODO: The model is evaluated and the fit statistic
+  (in this case the sum of squared residuals) computed and returned.
+
+The special method ``Chi2Fitter.__call__`` runs the minimisation:
+
+* The main inputs are the ``model`` and the data ``x`` and ``y``.
+  Extra arguments like ``maxiter`` to control the fitting can be added. 
+* Line TODO: Import the `~scipy.optimize.fmin_slsqp` minimisation function.
+  Lines TODO to TODO: Checks if the model is fittable and has valid constraints.
+* Line TODO: Make a copy of the model, because fitters are supposed to leave the state of the
+  input model unchanged and return a copy with optimised parameters.
+* Line TODO: Call the `~astropy.modeling.fitting.Fitter._model_to_fit_params` helper function
+  to extract the input parameter values as a numpy array ``x0``.
+* Line TODO: Call `~scipy.optimize.fmin_slsqp`, passing the fit statistic function ``self.errorfunc``,
+  the initial parameter values ``x0``, the ``model`` and data ``x`` and ``y`` as well a the ``iter`` parameter
+  that defines the maximum allowed number of ``self.errorfunc`` evaluations. 
+* Line TODO: Call the `~astropy.modeling.fitting.Fitter._fitter_to_model_params` helper function
+  to copy the ``fitparams`` returned by ``fmin_slsqp`` into the ``model_copy``.
+
+Usage
+-----
+
+You can use this user-defined `Chi2Fitter` just like the built-in astropy fitters:: 
+
+   import numpy as np
+   from astropy.modeling import models
+
+   g = models.Gaussian1D(amplitude=1.2, mean=0.9, stddev=0.5)
+
+   # Generate fake data
+   np.random.seed(0)
+   x = np.linspace(-5., 5., 100)
+   y = g(x) + np.random.normal(0., 0.2, x.shape)
+
+   # Fit the data using a Gaussian
+   fitter = Chi2Fitter()
+   g = fitter(g, x, y, maxiter=50)
diff --git a/docs/modeling/new.rst → docs/modeling/new_model.rst b/docs/modeling/new.rst → docs/modeling/new_model.rst
@@ -1,5 +1,5 @@
 **************************
-Defining New Model Classes
+Defining a new model class
 **************************
 
 This document describes how to add a model to the package or to create a
@@ -8,44 +8,16 @@ write an eval function which evaluates the model.  If the model is fittable, a
 function to compute the derivatives with respect to parapemeters is required
 if a linear fitting algorithm is to be used and optional if a non-linear fitter is to be used.
 
-
 Custom 1D models
 ----------------
 
 For 1D models, the `~astropy.modeling.functional_models.custom_model_1d`
 decorator is provided to make it very easy to define new models. The following
 example demonstrates how to set up a model consisting of two Gaussians:
 
-.. plot::
+.. plot:: modeling/custom_model.py
    :include-source:
 
-    import numpy as np
-    from astropy.modeling.models import custom_model_1d
-    from astropy.modeling.fitting import NonLinearLSQFitter
-
-    # Define model
-    @custom_model_1d
-    def sum_of_gaussians(x, amplitude1=1., mean1=-1., sigma1=1.,
-                            amplitude2=1., mean2=1., sigma2=1.):
-        return (amplitude1 * np.exp(-0.5 * ((x - mean1) / sigma1)**2) +
-                amplitude2 * np.exp(-0.5 * ((x - mean2) / sigma2)**2))
-
-    # Generate fake data
-    np.random.seed(0)
-    x = np.linspace(-5., 5., 200)
-    m_ref = sum_of_gaussians(amplitude1=2., mean1=-0.5, sigma1=0.4,
-                             amplitude2=0.5, mean2=2., sigma2=1.0)
-    y = m_ref(x) + np.random.normal(0., 0.1, x.shape)
-
-    # Fit model to data
-    m_init = sum_of_gaussians()
-    fit = NonLinearLSQFitter()
-    m = fit(m_init, x, y)
-
-    # Plot the data and the best fit
-    plt.plot(x, y, 'o', color='k')
-    plt.plot(x, m(x), color='r', lw=2)
-
 A step by step definition of a 1D Gaussian model
 ------------------------------------------------
 
@@ -157,85 +129,6 @@ with constraints.::
 A full example of a LineModel
 -----------------------------
 
-::
-
-    from astropy.modeling import models, Parameter, format_input
-    import numpy as np
-
-    class LineModel(models.PolynomialModel):
-        slope = Parameter('slope')
-        intercept = Parameter('intercept')
-        linear = True
-
-    def __init__(self, slope, intercept, param_dim=1, **constraints):
-        super(LineModel, self).__init__(slope=slope, intercept=intercept,
-                                        param_dim=param_dim, **constraints)
-        self.domain = [-1, 1]
-        self.window = [-1, 1]
-        self._order = 2
-
-    @staticmethod
-    def eval(x, slope, intercept):
-        return slope * x + intercept
-
-    @staticmethod
-    def fit_deriv(x, slope, intercept):
-        d_slope = x
-        d_intercept = np.ones_like(x)
-        return [d_slope, d_intercept]
-
-    @format_input
-    def __call__(self, x):
-        return self.eval(x, *self.param_sets)
-
+.. literalinclude:: line_model.py
+   :linenos:
 
-***************************
-Defining New Fitter Classes
-***************************
-
-This section describes how to add a new nonlinear fitting algorithm to this
-package or write a user-defined fitter.  In short, one needs to define an error
-function and a ``__call__`` method and define the types of constraints which
-work with this fitter (if any).
-
-The details are described below using scipy's SLSQP algorithm as an example.
-The base class for all fitters is `~astropy.modeling.fitting.Fitter`::
-
-    class SLSQPFitter(Fitter):
-        supported_constraints = ['bounds', 'eqcons', 'ineqcons', 'fixed', 'tied']
-
-        def __init__(self):
-            super(SLSQPFitter,self).__init__()
-
-All fitters take a model (their ``__call__`` method modifies the model's
-parameters) as their first argument.
-
-Next, the error function takes a list of parameters returned by an iteration of
-the fitting algorithm and input coordinates, evaluates the model with them and
-returns some type of a measure for the fit.  In the example the sum of the
-squared residuals is used as a measure of fitting.::
-
-    def errorfunc(self, fps, *args):
-        model = args[0]
-        meas = args[-1]
-        model.fitparams(fps)
-        res = self.model(*args[1:-1]) - meas
-        return np.sum(res**2)
-
-The ``__call__`` method performs the fitting. As a minimum it takes all
-coordinates as separate arguments. Additional arguments are passed as
-necessary.::
-
-    def __call__(self, model, x, y , maxiter=MAXITER, epsilon=EPS):
-        if model.linear:
-                raise ModelLinearityException(
-                    'Model is linear in parameters; '
-                    'non-linear fitting methods should not be used.')
-        model_copy = model.copy()
-        init_values, _ = model_copy._model_to_fit_params()
-        self.fitparams = optimize.fmin_slsqp(self.errorfunc, p0=init_values,
-                                             args=(y, x),
-                                             bounds=self.bounds,
-                                             eqcons=self.eqcons,
-                                             ineqcons=self.ineqcons)
-        return model_copy