Merge pull request #1203 from adonath/new_models_#980_and_#981

Simplified implementing new models #980

astropy/modeling/core.py

@@ -47,8 +47,10 @@
 from . import constraints
 from .utils import InputParameterError
 
+
 __all__ = ['Model', 'ParametricModel', 'PCompositeModel', 'SCompositeModel',
-           'LabeledInput', '_convert_input', '_convert_output']
+           'LabeledInput', '_convert_input', '_convert_output',
+           'Parametric1DModel']
 
 
 def _convert_input(x, pdim):
@@ -705,7 +707,7 @@ def __init__(self, transforms, inmap=None, outmap=None, n_inputs=None, n_outputs
     def inverse(self):
         try:
             transforms = [tr.inverse() for tr in self._transforms[::-1]]
-        except NotIMplementedError:
+        except NotImplementedError:
             raise
         if self._inmap is not None:
             inmap = self._inmap[::-1]
@@ -824,3 +826,48 @@ def __call__(self, *data):
             for tr in self._transforms:
                 result += tr(*data)
             return result
+
+
+class Parametric1DModel(ParametricModel):
+    """
+    Base class for one dimensional parametric models
+
+    This class provides an easier interface to defining new models.
+    Examples can be found in functional_models.py
+
+    Parameters
+    ----------
+    parameter_dict : dictionary
+        Dictionary of model parameters with initialisation values
+        {'parameter_name': 'parameter_value'}
+
+    """
+    deriv = None
+    linear = False
+
+    def __init__(self, param_dict, **cons):
+        # Get parameter dimension
+        param_dim = np.size(param_dict[self.param_names[0]])
+
+        # Initialize model parameters. This is preliminary as long there is
+        # no new parameter class. It may be more reasonable and clear to init 
+        # the parameters in the model constructor itself, with constraints etc.
+        for param_name in self.param_names:
+            setattr(self, "_" + param_name, parameters.Parameter(name=param_name,
+                            val=param_dict[param_name], mclass=self, param_dim=param_dim))
+
+        super(Parametric1DModel, self).__init__(self.param_names, n_inputs=1,
+                                                n_outputs=1, param_dim=param_dim, **cons)
+
+    def __call__(self, x):
+        """
+        Transforms data using this model.
+
+        Parameters
+        ----------
+        x : array like or a number
+            input
+        """
+        x, fmt = _convert_input(x, self.param_dim)
+        result = self.eval(x, *self.param_sets)
+        return _convert_output(result, fmt)

astropy/modeling/fitting.py

@@ -81,10 +81,6 @@ def __init__(self, model):
             self._fitpars = self.model.constraints.fitpars[:]
         else:
             self._fitpars = self.model._parameters[:]
-        if self._model.deriv is None:
-            self.dfunc = None
-        else:
-            self.dfunc = self._wrap_deriv
         self._weights = None
 
     @property
@@ -172,22 +168,22 @@ def _wrap_deriv(self, p, x, y, z=None):
         if fixed_and_tied:
             pars = self.model.constraints.modelpars
             if z is None:
-                fullderiv = self.model.deriv(pars, x, y)
+                fullderiv = np.array(self.model.deriv(x, *pars))
             else:
-                fullderiv = self.model.deriv(pars, x, y, z)
+                fullderiv = np.array(self.model.deriv(x, y, *pars))
             ind = range(len(self.model.param_names))
             for name in fixed_and_tied:
                 index = self.model.param_names.index(name)
                 ind.remove(index)
-            res = np.empty((fullderiv.shape[0], fullderiv.shape[1] - len(ind)))
-            res = fullderiv[:, ind]
+            res = np.empty((fullderiv.shape[1] - len(ind), fullderiv.shape[0]))
+            res = fullderiv[ind, :]
             return res
         else:
             pars = p[:]
             if z is None:
-                return self.model.deriv(pars, x, y)
+                return self.model.deriv(x, *pars)
             else:
-                return self.model.deriv(pars, x, y, z)
+                return self.model.deriv(x, y, *pars)
 
     def _validate_constraints(self):
         fname = self.__class__.__name__
@@ -483,7 +479,8 @@ def covar(self):
             log.info("Could not construct a covariance matrix")
             return None
 
-    def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER, epsilon=EPS):
+    def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER,
+                                epsilon=EPS, estimate_jacobian=False):
         """
         Fit data to this model.
 
@@ -505,7 +502,10 @@ def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER, epsilon=EPS):
             epsfcn is less than the machine precision, it is
             assumed that the relative errors in the functions are
             of the order of the machine precision.
-
+        estimate_jacobian : bool
+            If False (default) and if the model has a deriv method,
+            it will be used. Otherwise the Jacobian will be estimated.
+            If True, the Jacobian will be estimated in any case.
         """
         from scipy import optimize
         x = np.asarray(x, dtype=np.float)
@@ -527,9 +527,14 @@ def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER, epsilon=EPS):
             meas = np.asarray(z, dtype=np.float)
             farg = (x, y, meas)
 
+        if self._model.deriv is None or estimate_jacobian:
+            self.dfunc = None
+        else:
+            self.dfunc = self._wrap_deriv
+
         self.fitpars, status, dinfo, mess, ierr = optimize.leastsq(
             self.errorfunc, self.fitpars, args=farg, Dfun=self.dfunc,
-            maxfev=maxiter, epsfcn=epsilon, full_output=True)
+            col_deriv=1, maxfev=maxiter, epsfcn=epsilon, full_output=True)
         self.fit_info.update(dinfo)
         self.fit_info['status'] = status
         self.fit_info['message'] = mess
@@ -737,7 +742,7 @@ def errorfunc(self, fps, *args):
                     plen = sl.stop - sl.start
                     mpars.extend(mfpars[:plen])
                     del mfpars[:plen]
-            modelfit = model.eval(margs[:-1], mpars)
+            modelfit = model.eval(margs[:-1], *mpars)
             fitted.extend(modelfit - margs[-1])
         return np.ravel(fitted)