Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Models subpackage #493

Merged
merged 76 commits into from
@nden
Collaborator

This PR is for the initial import of models and fitting routines.
The documentation is here:

https://github.com/nden/astropy/tree/master/docs/models

I have included some of the tests that I was using during development. Various parts of this package were tested against iraf, pywcs and matlab but not all tests are included here (some of the data files are large).

@eteq
Owner

For those who might not be aware, the idea to add this work to astropy came out of a discussion at the coordination meeting, and is a necessary step for the "generalized" WCS project. But It also should be very useful as a model-fitting interface that can be shared throughout astropy.

astropy/models/__init__.py
@@ -0,0 +1,7 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+import models
+import fitting
+import parameters
+import projections
+import rotations
+
@embray Collaborator
embray added a note

These should all be relative imports such as from . import models

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((94 lines not shown))
+ """
+ Parameters
+ ----------
+ model: an instance of `fitting.models.ParametricModel`
+
+ """
+ self.model = model
+ if not self.model.linear:
+ raise ModelLinearityException('Model is not linear in parameters, '
+ 'linear fit methods should not be used.')
+ self.fit_info = {'residuals': None,
+ 'rank': None,
+ 'singular_values': None,
+ 'pars': None
+ }
+ Fitter.__init__(self, model, fixed=fixed)
@embray Collaborator
embray added a note

Although not as important in single-inheritance cases, please still try to use super() instead. For example super(Fitter, self).__init__(model, fixed=fixed). See the seventh bullet point under http://astropy.readthedocs.org/en/latest/development/codeguide.html#coding-style-conventions

@nden Collaborator
nden added a note

OK, I'll change this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((3 lines not shown))
+This module provides wrappers, called Fitters, around some Numpy and Scipy
+fitting functions. All Fitters take an instance of `ParametricModel` as input
+and define a __call__ method which fits the model to the data and changes the
+model's parameters attribute. The idea is to make this extensible and allow
+users to easily add other fitters.
+
+Linear fitting is done using Numpy's linalg.lstsq function.
+There are currently two non-linear fitters which use leastsq and slsqp functions
+in scipy.optimize.
+
+"""
+from __future__ import division
+import abc
+import numpy as np
+from numpy import linalg
+from scipy import optimize
@astrofrog Owner

scipy is not a required dependency for astropy, so it should be imported only in the functions/methods that require it.

@nden Collaborator
nden added a note

Forgot about that, will make the changes

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
@@ -0,0 +1,825 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+This module provides wrappers, called Fitters, around some Numpy and Scipy
+fitting functions. All Fitters take an instance of `ParametricModel` as input
+and define a __call__ method which fits the model to the data and changes the
+model's parameters attribute. The idea is to make this extensible and allow
+users to easily add other fitters.
+
+Linear fitting is done using Numpy's linalg.lstsq function.
+There are currently two non-linear fitters which use leastsq and slsqp functions
+in scipy.optimize.
+
+"""
+from __future__ import division
@astrofrog Owner

For future-proofness, it would be best to also import print_function and change the print statements to the Python 3 syntax print().

@nden Collaborator
nden added a note

Good point! I actually haven't done much with python 3 yet. So I installed it and ran the tests. I'll change the print statements but will have to figure out why some tests are failing.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@astrofrog astrofrog referenced this pull request in astrofrog/astropy-api
Merged

Photometry API proposal #2

@astrofrog
Owner

@nden - just for information, I do plan to take a careful look at this, but it will be easier once 0.2 is out of the door. But at first glance, this looks great! I wonder whether it would already be worth isolating the WCS-specific stuff to astropy.models.wcs for clarity?

@nden
Collaborator

@astrofrog - I am going to push changes to master soon - mostly related to moving constraints from fitters to models, you may want to wait for this.
The only strictly WCS-specific modules so far are projections.py and rotations.py but models in models.py can be useful for WCS too. Since init.py will import all modules from models and models.wcs, then I'm OK with splitting them, since I won't have to remember where a specific module lives.

@taldcroft
Collaborator

@nden - sorry, I have been tied up with other things and haven't been able to look at these changes closely. From the comments and commit titles it looks like you made some changes that I'll like. :smiley: When things are quieter next week I'll get back to this.

@cdeil
Collaborator

@nden You commented here that you plan to

- Separate optimization methods form fit statistic
- Add parameter confidence intervals

It's not clear to me what the proper home for the fit statistics, optimization methods and confidence interval methods is. E.g. in Sherpa each has it's own sub-package (see here):

  • stats (fit statistics like various Likelihoods and Chi2 functions, see here)
  • optmethods (optimisation methods, see here
  • estmethods (confidence limit estimation methods, e.g. profile likelihood)

I think the fit statistics should go in astropy.stats, actually I'm implementing the Sherpa statistics (and other stats things) at the moment:
https://github.com/astropy/astropy/wiki/What-methods-do-we-want-in-astropy.stats%3F
https://groups.google.com/d/topic/astropy-dev/Zwgafam171E/discussion

Where the optimisers and confidence interval methods (well, at the moment the wrappers around scipy.optimize, not the actual implementations) should live is not clear to me. E.g. the Rolke method I propose to include in astropy.stats is really just a special case of the profile likelihood method. I guess they could go in a new astropy.optimize package or into one of the existing astropy.utils or astropy.stats or astropy.models packages, although in the astropy.models case maybe renaming it to astropy.fitting would be better, because e.g. I only expect to find models in astropy.models. Anyways, this is "just" naming, and I don't have a strong opinion what we should do, I just wanted to mention the question here.

@keflavich
Collaborator

I'm very interesting in using the astropy modeling tools, so if and when this gets merged, I'll work on a tutorial for it based on blackbody fitting - I'd like to translate this: https://code.google.com/p/agpy/source/browse/trunk/agpy/blackbody.py to astropy, so that it uses astropy's units framework + fitting framework.

But keeping the linked example in mind, I think it would be worthwhile to make "mpfit to astropy.models" and "lmfit-py to astropy.models" guides. I haven't looked at the PR closely yet, but I'd really like to know how it differs from https://github.com/newville/lmfit-py.

@nden
Collaborator

@keflavich It seems this is similar to lmfit-py in many ways. I think the main difference is that we have a formalized Model class and the Parameters (functionally similar to lmfit-py.Parameters) is attached to it but shared between models and fitters. This makes it possible to create composite models. astropy.models supports the same constraints, although after looking at lmfit-py I see that the bounds are calculated in a different way. It probably differs in other details as well. Also, the goal in astropy.models is to create a formal separation of optimization methods and statistics methods so that they can be combined arbitrarily.
I looked briefly at blackbody.py and I think it should be straightforward to convert it to astropy.models.

@nden
Collaborator

@cdeil Sorry for not responding for such a long time - I was working on other things while waiting for the 0.2 release and this fell through the cracks.

Having separate astropy subpackages for optimization, statistics and estimation methods may make sense in the future. But I think we should try to avoid package proliferation unless necessary. If code functionality is used by several packages it may make sense to separate it. I am not too happy about the name too but I can't come up with a better one. The initial name was "fitting" but at the last coordination meeting it was decided to rename it to "models". I have no strong feeling either way.

@astrofrog
Owner

@taldcroft - can you review this in the near future to see if it addresses your concerns?

I'll also review this within the next week.

@taldcroft
Collaborator

@astrofrog - I'll plan to review this within the next week as well.

astropy/models/constraints.py
@@ -0,0 +1,181 @@
+from __future__ import division, print_function
+import numpy as np
+class Constraints(object):
+ """
+ Fitting constraints
+
+ """
+
+ def __init__(self, model, fixed={}, tied={}, bounds={},
+ eqcons=[], ineqcons=[]):
+ """
+ Parameters
+ ----------
+ fitter: object which supports fitting
@mdboom Collaborator
mdboom added a note

It's annoying, but there needs to be a space on both sides of the colon in numpydoc parameter lists, otherwise it doesn't see the separation between the name of the argument and the type.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/constraints.py
((12 lines not shown))
+ Parameters
+ ----------
+ fitter: object which supports fitting
+ fixed: dict
+ (parameter_name: True/False}
+ parameters to be held fixed during fitting
+ tied: dict
+ keys are parameter names
+ values are callable/function providing a relationship
+ between parameters. Currently the callable takes a model
+ instance as an argument.
+ In the example below xcen is tied to the value of xsigma
+
+ def tie_center(model):
+ xcen = 50*model.xsigma
+ return xcen
@mdboom Collaborator
mdboom added a note

I'm not sure this kind of formatting works here -- we'll have to check. This might have to be moved to an examples section below.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/constraints.py
((2 lines not shown))
+import numpy as np
+class Constraints(object):
+ """
+ Fitting constraints
+
+ """
+
+ def __init__(self, model, fixed={}, tied={}, bounds={},
+ eqcons=[], ineqcons=[]):
+ """
+ Parameters
+ ----------
+ fitter: object which supports fitting
+ fixed: dict
+ (parameter_name: True/False}
+ parameters to be held fixed during fitting
@mdboom Collaborator
mdboom added a note

I don't quite understand what fixed is supposed to be from this docstring.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/constraints.py
((6 lines not shown))
+
+ """
+
+ def __init__(self, model, fixed={}, tied={}, bounds={},
+ eqcons=[], ineqcons=[]):
+ """
+ Parameters
+ ----------
+ fitter: object which supports fitting
+ fixed: dict
+ (parameter_name: True/False}
+ parameters to be held fixed during fitting
+ tied: dict
+ keys are parameter names
+ values are callable/function providing a relationship
+ between parameters. Currently the callable takes a model
@mdboom Collaborator
mdboom added a note

This needs some punctuation, since Sphinx won't preserve the line endings here. As in:

Keys are parameter names; values are callables/functions providing a 
relationship between parameters.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/constraints.py
((27 lines not shown))
+ return xcen
+
+ tied ={'xcen':tie_center}
+ bounds: dict
+ keys: parameter names
+ values: list of length 2 giving the desired range for hte parameter
+ eqcons: list
+ A list of functions of length n such that
+ eqcons[j](x0,*args) == 0.0 in a successfully optimized
+ problem.
+ ineqcons: list
+ A list of functions of length n such that
+ ieqcons[j](x0,*args) >= 0.0 in a successfully optimized
+ problem.
+ """
+ self.model = model
@mdboom Collaborator
mdboom added a note

Should we (for future flexibility) store model in a private member with a public property (like all the other members in this class)?

@nden Collaborator
nden added a note

@mdboom The guide says "Classes should either use direct variable access, or python’s property mechanism ...". Am I wrong in thinking that a property makes sense mostly when something more than just returning or setting the variable needs to be done? The guide for me in this code was to always set model parameters with properties since they are the variables users will interact with. Constraints and fitters keep an instamce of the model as an attribute because they need to manipulate its parameters but users don't need to access constraints.model or fitter.model. I don't see a benefit in making 'model' a property. A private variable - may be but this seems a bit arbitrary. Does this make sense?

@mdboom Collaborator
mdboom added a note
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((31 lines not shown))
+ }
+
+class ModelLinearityException(Exception):
+ """
+ Called when a linear model is passed to a non-linear fitter and vice versa.
+ """
+ pass
+
+class Fitter(object):
+ """
+ Base class for all fitters.
+
+ The purpose of this class is to deal with constraints
+ """
+ def __init__(self, model):
+ __metaclass__ = abc.ABCMeta
@mdboom Collaborator
mdboom added a note

I think the __metaclass__ declaration needs to be at the top level of the class. i.e. outside of the __init__ function. All this does is assign it to a local variable, if I understand correctly.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((33 lines not shown))
+class ModelLinearityException(Exception):
+ """
+ Called when a linear model is passed to a non-linear fitter and vice versa.
+ """
+ pass
+
+class Fitter(object):
+ """
+ Base class for all fitters.
+
+ The purpose of this class is to deal with constraints
+ """
+ def __init__(self, model):
+ __metaclass__ = abc.ABCMeta
+
+ self.model = model
@mdboom Collaborator
mdboom added a note

As above, maybe this should be private with a property?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@mdboom mdboom commented on the diff
astropy/models/fitting.py
((97 lines not shown))
+ except NotImplementedError:
+ pass
+ self._fitpars[:] = fps
+ self.model.parameters = fps
+
+ def _set_bounds(self, pars):
+ """
+ This method is to be implemented by subcclasses of Fitter if necessary
+
+ Diferent fitting algorithms deal with bounds in a different way.
+ For example, the SLSQP algorithm accepts bounds as input while
+ the leastsq algorithm does not handle bounds at all and they are
+ dealt with in a separate method.
+
+ """
+ raise NotImplementedError
@mdboom Collaborator
mdboom added a note

raise NotImplementedError()

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((105 lines not shown))
+
+ Diferent fitting algorithms deal with bounds in a different way.
+ For example, the SLSQP algorithm accepts bounds as input while
+ the leastsq algorithm does not handle bounds at all and they are
+ dealt with in a separate method.
+
+ """
+ raise NotImplementedError
+
+ def _wrap_deriv(self, p, x, y, z=None):
+ """
+ Wraps the method calculating the Jacobian of the function to
+ account for model constraints
+
+ Currently the only fitter that uses a derivative is the
+ NonLinearLSQFitter. This wrapper may neeed to be revised when
@mdboom Collaborator
mdboom added a note

Put backticks around identifiers so they will be hyperlinked. e.g.

`NonLinearLSQFitter`
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((106 lines not shown))
+ Diferent fitting algorithms deal with bounds in a different way.
+ For example, the SLSQP algorithm accepts bounds as input while
+ the leastsq algorithm does not handle bounds at all and they are
+ dealt with in a separate method.
+
+ """
+ raise NotImplementedError
+
+ def _wrap_deriv(self, p, x, y, z=None):
+ """
+ Wraps the method calculating the Jacobian of the function to
+ account for model constraints
+
+ Currently the only fitter that uses a derivative is the
+ NonLinearLSQFitter. This wrapper may neeed to be revised when
+ when other fitters using function derivative are added or when
@mdboom Collaborator
mdboom added a note

when is repeated here.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((150 lines not shown))
+ return self.model.deriv(pars, x, y)
+ else:
+ return self.model.deriv(pars, x, y, z)
+
+ if z is None:
+ res = self.model.deriv(pars, x, y)
+ else:
+ res = self.model.deriv(pars, x, y, z)
+ return res
+
+ def _validate_constraints(self):
+ fname = self.__class__.__name__
+ try:
+ c = constraintsdef[fname]
+ except KeyError:
+ print("%s does not support fitting with constraints" % fname)
@mdboom Collaborator
mdboom added a note

Rather than using print, this should probably be a warning, so it can be logged and controlled etc. Also, it seems like the purpose of this is to explain why the KeyError is raised. Maybe it makes sense to just raise a different error with this message in it?

@nden Collaborator
nden added a note

You are right, this is to explain the KeyError, so I changed it to raise a custom Error.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/constraints.py
((44 lines not shown))
+ self._fixed = fixed
+ self._tied = tied
+ self._bounds = bounds
+ self.set_range(bounds)
+ self._eqcons = eqcons
+ self._ineqcons = ineqcons
+ if any(self._fixed.values()) or any(self._tied.values()):
+ self._fitpars = self._model_to_fit_pars()
+ else:
+ self._fitpars = self.model._parameters[:]
+
+ def __str__(self):
+ fmt = ""
+ if any(self._fixed.values()):
+ fixedstr = [par for par in self._fixed if self._fixed[par]]
+ fmt += "fixed=%s, " % str(fixedstr)
@mdboom Collaborator
mdboom added a note

Our astropy coding guidelines suggest using new-style formatting (str.format) rather than % whenever possible. (There are some weird side cases where we don't because the behavior has changed, but for new code like this we should use format, I think).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((18 lines not shown))
+
+from .util import pmapdomain
+
+__all__ = ['LinearLSQFitter', 'NonLinearLSQFitter', 'SLSQPFitter',
+ 'JointFitter']
+
+MAXITER = 100
+EPS = np.sqrt(np.finfo(float).eps)
+
+# supported constraints
+constraintsdef = {'NonLinearLSQFitter': ['fixed', 'tied', 'bounds'],
+ 'SLSQPFitter': ['bounds', 'eqcons', 'ineqcons', 'fixed', 'tied'],
+ 'LinearLSQFitter': ['fixed'],
+ }
+
+class ModelLinearityException(Exception):
@mdboom Collaborator
mdboom added a note

I've been bitten by this myself -- most of the built-in Python exceptions are named SomethingError, not SomethingException, so we should probably follow that convention.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((228 lines not shown))
+ if y is None:
+ d = self.model.deriv(x)
+ else:
+ d = self.model.deriv(x, y)
+ fixed = [name for name in self.model.constraints.fixed if
+ self.model.constraints.fixed[name]]
+ ind = range(len(self.model.parnames))
+ for name in fixed:
+ index = self.model.parnames.index(name)
+ ind.remove(index)
+ res = d[:, ind]
+ return res
+
+ def __call__(self, x, y, z=None, w=None, rcond=None):
+ """
+ Parameters
@mdboom Collaborator
mdboom added a note

There should be a short description of what this does here. i.e., why would one "call" this instance?

@mdboom Collaborator
mdboom added a note

I see the answer is documented in the models.py module docstring -- that's probably ok for the detailed explanation, but maybe all of these just say "Transforms data using the fitter." or something generic like that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((245 lines not shown))
+ x: array
+ input coordinates
+ y: array
+ input coordinates
+ z: array (optional)
+ input coordinates
+ w: array (optional)
+ weights
+ rcond: float, optional
+ Cut-off ratio for small singular values of `a`.
+ Singular values are set to zero if they are smaller than `rcond`
+ times the largest singular value of `a`.
+
+ """
+ multiple = False
+ x = np.asarray(x) + 0.0
@mdboom Collaborator
mdboom added a note

What is the + 0.0 here for? Could it be handled by explicitly setting a dtype in np.asarray?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((253 lines not shown))
+ rcond: float, optional
+ Cut-off ratio for small singular values of `a`.
+ Singular values are set to zero if they are smaller than `rcond`
+ times the largest singular value of `a`.
+
+ """
+ multiple = False
+ x = np.asarray(x) + 0.0
+ y = np.asarray(y) + 0.0
+
+ if self.model.ndim == 2 and z is None:
+ raise ValueError("Expected x, y and z for a 2 dimensional model.")
+
+ if z is None:
+ if x.shape[0] != y.shape[0]:
+ raise ValueError("x and y should have the same length")
@mdboom Collaborator
mdboom added a note

Do they need to be the same length, or is it good enough that they are broadcastable? Could it be made such?

@nden Collaborator
nden added a note

If x and y are 1D arrays I think it makes sense to require them to be the same size. The intention here is to be able to fit arrays of shapes (n,) and (n,m) where m different sets of y are fit to the same x. But I can see that the error message is not very helpful. May be this sounds better?
"Expected the measured and model data to have the same size"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/fitting.py
((388 lines not shown))
+ self.fitpars = fps
+ meas = args[0]
+ if self.weights is None:
+ return np.ravel(self.model(*args[1:]) - meas)
+ else:
+ return np.ravel(self.weights * (self.model(*args[1:]) - meas))
+
+ def _set_bounds(self, fitpars):
+ if any(c != (-1E12, 1E12) for c in
+ self.model.constraints.bounds.values()):
+ bounds = [self.model.constraints.bounds[par] for
+ par in self.model.parnames]
+ [setattr(self.model, name, par if par>b[0] else b[0]) for
+ name, par, b in zip(self.model.parnames, fitpars, bounds)]
+ [setattr(self.model, name, par if par<b[1] else b[1]) for
+ name, par, b in zip(self.model.parnames, fitpars, bounds)]
@mdboom Collaborator
mdboom added a note

This is probably a matter of personal style, but I'm generally not a fan of list comprehensions that don't get assigned to anything. I'd rather just see this as a regular for loop, particularly since both of these are iterating over the same thing.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/models.py
((47 lines not shown))
+from . import parameters
+from . import constraints
+from .util import pmapdomain, InputParametersException, comb
+
+__all__ = ['ChebyshevModel', 'Gauss1DModel', 'Gauss2DModel', 'ICheb2DModel',
+ 'ILegend2DModel', 'LegendreModel', 'Poly1DModel', 'Poly2DModel',
+ 'ScaleModel', 'ShiftModel', 'SIPModel',
+ 'PCompositeModel', 'SCompositeModel', 'LabeledInput']
+
+def _convert_input(x, pdim):
+ """
+ format the input into appropriate shape
+
+ 'N' - normal (don't do anything to the output)
+ 'T' - transposed
+ 'S' - scalar
@mdboom Collaborator
mdboom added a note

This should be reformatted in numpydoc format.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/models.py
((117 lines not shown))
+ "have the required shape" % name)
+ else:
+ setattr(self, '_'+name, par)
+ self._parameters = parameters.Parameters(self,
+ self.parnames,
+ paramdim=self.paramdim)
+ else:
+ setattr(self, '_'+name, val)
+ else:
+ par = parameters._Parameter(name, val, self, self.paramdim)
+ if not getattr(self, name).parshape == par.parshape:
+ raise InputParametersException(
+ "Input parameter '%s' does not "
+ "have the required shape" % name)
+ else:
+ setattr(self, '_'+name, par)
@mdboom Collaborator
mdboom added a note

This getpar/setpar stuff could be done perhaps more easily with a descriptor class. You basically create a class with __get__ and __set__ members, and then you can assign instances of that class to another class like you would a property. I think that's more straightforward than this plus the lambda/property magic used to define them below.

See this: http://docs.python.org/2/howto/descriptor.html

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/models.py
((210 lines not shown))
+ """
+ raise NotImplementedError
+
+ def add_model(self, newtr, mode):
+ """
+ Parameters
+ ----------
+ newtr: an instance of a subclass of Model
+ mode: string
+ 'parallel', 'serial', 'p' or 's'
+ a flag indicating whether to combine the models
+ in series or in parallel
+
+ Returns
+ -------
+ an instance of CompositeModel
@mdboom Collaborator
mdboom added a note

The Returns section needs to be in the same parameter format, i.e.:

model : CompositeModel
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/models.py
((828 lines not shown))
+ x: array, of minimum dimensions 1
+
+ Note: See the module docstring for rules for model evaluation.
+ """
+ x, fmt = _convert_input(x, self.paramdim)
+ result = self.horner(x, self.psets)
+ return _convert_output(result, fmt)
+
+class Poly2DModel(PModel):
+ """
+ 2D Polynomial model
+
+ Represents a general polynomial of degree n:
+
+ P(x,y) = c0_0 + c1_0*x + ...+ cn_0*x**n + c0_1*y + ...+ c0_n*y**n +
+ c1_1*x*y + c1_2*x*y**2 + ... + c1_(n-1)*x*y**(n-1)+ ... + c(n-1)_1*x**(n-1)*y
@mdboom Collaborator
mdboom added a note

This would be nice as a math directive so it gets rendered nicely as LaTeX in the docs.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/parameters.py
((47 lines not shown))
+
+ """
+ def __init__(self, name, val, mclass, paramdim, fixed=False, tied=False,
+ minvalue=None, maxvalue=None):
+ """
+ Parameters
+ ----------
+ name: string
+ parameter name
+ val: number or an iterable of numbers
+ mclass: object
+ an instance of a Model class
+ paramdim: int
+ parameter dimension
+ """
+ self.paramdim = paramdim
@mdboom Collaborator
mdboom added a note

private member with property?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@mdboom mdboom commented on the diff
astropy/models/projections.py
((14 lines not shown))
+.. [2] Calabretta, M.R., Greisen, E.W., 2002, A&A, 395, 1077 (Paper II)
+
+"""
+from __future__ import division, print_function
+from .models import Model
+from . import parameters
+import numpy as np
+
+projcodes = ['TAN', 'AZP', 'SZP', 'STG', 'SIN', 'ARC', 'ZPN', 'ZEA', 'AIR',
+ 'CYP', 'CEA', 'MER']
+
+__all__ = ['Pix2Sky_AZP', 'Sky2Pix_AZP', 'Pix2Sky_CAR', 'Sky2Pix_CAR',
+ 'Pix2Sky_CEA', 'Sky2Pix_CEA', 'Pix2Sky_COP', 'Sky2Pix_COP',
+ 'Pix2Sky_CYP', 'Sky2Pix_CYP', 'Pix2Sky_MER', 'Sky2Pix_MER',
+ 'Pix2Sky_SIN', 'Sky2Pix_SIN', 'Pix2Sky_STG', 'Sky2Pix_STG',
+ 'Pix2Sky_TAN', 'Sky2Pix_TAN']
@mdboom Collaborator
mdboom added a note

I understand these three-letter projection codes come from FITS WCS and are probably familiar to many, but I for one always have trouble remembering them. Maybe we should have aliases with longer names, e.g.:

Pix2Sky_Tangent = Pix2Sky_TAN
@mdboom Collaborator
mdboom added a note

Also -- question as to whether we want to use the word sky or world here. I know wcslib uses sky, but I think @astrofrog at one point explained that a better word would be world, but I leave that up for others to decide.

@nden Collaborator
nden added a note

I agree in many cases it's appropriate to have 'world' instead of 'sky'. In this case though the world coordinates are really sky coordinates. But I don't have strong feelings about this and would also leave it to other people to decide.

@wkerzendorf Collaborator

I prefer world as in an n-dimensional theory simulation where axes map to e.g. density, temperature and hydrogen abundance - sky would not make a lot of sense.

@nden Collaborator
nden added a note

In my mind 'world' in the context of WCS means a physical quantity. The classes in projections.py represent geometrical relations. It just happens that 'sphere' in astronomy is 'sky'. We could abstract the names even more and name them something like 'pix_to_sphere' or 'plain_to_sphere' but certainly not 'world'. Unless I'm mistaken...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@mdboom mdboom commented on the diff
astropy/models/tests/data/idcompspec.fits
@@ -0,0 +1,36 @@
+# Thu 16:08:18 31-Mar-2011
@mdboom Collaborator
mdboom added a note

This file has a fits extension, but doesn't look like a FITS file. Is that deliberate?

@nden Collaborator
nden added a note

It's the output of IRAF's identify task. The naming convention for their database file is to prepend 'id' to the fits file name. It doesn't really matter what the file name is for this tests since I just parse the IRAF database file to compare the fit results. I can change it if it's too confusing.

@mdboom Collaborator
mdboom added a note

Ok -- it's probably fine in that case if that's the convention. @eteq: As the "data" guy, what do you think?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/tests/testfitters.py
@@ -0,0 +1,159 @@
+"""
+Module to test fitting routines
+"""
+from __future__ import division
+import os
+from .. import models, fitting
+from . import irafutil
+import numpy as np
+from numpy import linalg
+from numpy.testing import utils
+from scipy import optimize
+from data import dpath
@mdboom Collaborator
mdboom added a note

To load data files, astropy has utils.data.get_pkg_data_filename and friends. I have nothing against with this approach here, but maybe we should use utils.data to be consistent with other packages?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
@@ -0,0 +1,102 @@
+.. _fitting:
+
+*******
+Fitting
+*******
+
+This module provides wrappers, called Fitters, around some Numpy and Scipy
+fitting functions. All Fitters take an instance of
+**`models.ParametricModel`** as input and define a **__call__** method
@mdboom Collaborator
mdboom added a note

When an identifier is in both ** and backticks, it doesn't properly create the link. I think you just want backticks for all of these.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
@@ -0,0 +1,102 @@
+.. _fitting:
+
+*******
+Fitting
+*******
+
+This module provides wrappers, called Fitters, around some Numpy and Scipy
+fitting functions. All Fitters take an instance of
+**`models.ParametricModel`** as input and define a **__call__** method
+which fits the model to the data and changes **model.parameters**
+attribute. The idea is to make this extensible and allow users to easily add
+other fitters.
+
+Linear fitting is done using Numpy's **linalg.lstsq** function.
@mdboom Collaborator
mdboom added a note

If you do

`~numpy.linalg.lstsq`

here, it will link to the function in Numpy's documentation and only display the last part of the name, which is what I think you want.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
@@ -0,0 +1,155 @@
+**************************
+Models (`astropy.models`)
+**************************
+
+Introduction
+============
+The **`models`** and **`fitting`** modules described here are designed to work as
@mdboom Collaborator
mdboom added a note

This paragraph does a good job of displaying how astropy.models is designed. There should be a paragraph before this describing what it does and why someone would want to use it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@mdboom mdboom commented on the diff
docs/models/index.rst
((52 lines not shown))
+>>> ch1.parameters
+[0.0, 0.0, 0.0, 0.0]
+>>> ch1.parameters = [1,2,3,4]
+>>> ch1.parameters
+[1.0, 2.0, 3.0, 4.0]
+>>> print ch1
+Model: ChebyshevModel
+Dim: 1
+Degree: 3
+Parameter sets: 1
+Parameters:
+ c0: [1.0]
+ c1: [2.0]
+ c2: [3.0]
+ c3: [4.0]
+>>> y = ch1(x)
@mdboom Collaborator
mdboom added a note

This comes many lines after where x is defined. Maybe move the definition of x down here?

@mdboom Collaborator
mdboom added a note

Then, it would be illustrative to print out what y is at this point.

@nden Collaborator
nden added a note

y is too big and, at least to me, not very illustrative. I'd rather add a picture like here

http://ssb.stsci.edu/doc/jwst_dev/jwstlib.models.doc/html/models.html#models-examples

using

.. figure:: images/gauss1D_eval_psets.png
:scale: 75 %

Would this work in this setup?

@mdboom Collaborator
mdboom added a note

Sure -- I agree a plot would be even better. No reason we can't have lots of plots in the docs wherever they make sense.

@nden Collaborator
nden added a note

@mdboom I can't get images to scale using the :scale: option. Does this work?
.. image:: images/p1d_5psets.png
:scale: 25 %

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((55 lines not shown))
+>>> ch1.parameters
+[1.0, 2.0, 3.0, 4.0]
+>>> print ch1
+Model: ChebyshevModel
+Dim: 1
+Degree: 3
+Parameter sets: 1
+Parameters:
+ c0: [1.0]
+ c1: [2.0]
+ c2: [3.0]
+ c3: [4.0]
+>>> y = ch1(x)
+>>> n=np.random.randn(90)
+>>> ny = y + n
+
@mdboom Collaborator
mdboom added a note

Why are random numbers being added to y here?

@nden Collaborator
nden added a note

I just add some noise to the data in this way.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((72 lines not shown))
+Fit a Chebyshev polynomial to the data
+
+>>> ch2 = models.ChebyshevModel(3)
+>>> chfit = fitting.LinearLSQFitter(ch2)
+>>> chfit(x,y)
+>>> ch2.parameters
+[0.7863153, 1.7473515, 2.8038203, 3.9106717]
+
+- Working with 2D models
+
+First create some data to be fitted with a 2D polynomial
+
+>>> x, y = np.mgrid[:10, :10]
+>>> def poly2(x, y):
+ return 1+2*x+3*x**2+4*y+5*y**2+6*x*y
+z = poly2(x, y)
@mdboom Collaborator
mdboom added a note

Should this have >>> in front?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@mdboom mdboom commented on the diff
docs/models/index.rst
((116 lines not shown))
+<Gauss1DModel(amplitude= [10.0],xcen= [4.2000000000000002],xsigma= [2.1000000000000001],paramdim=1)>
+>>> y = g1(x)
+>>> n = np.random.randn(90)
+>>> ny = y + n
+>>> gfit = fitting.NonLinearLSQFitter(g1)
+>>> gfit(x, ny)
+>>> print g1
+Model: Gauss1DModel
+Dim: 1
+Degree: N/A
+Parameter sets: 1
+Parameters:
+ amplitude: [10.141697966089579]
+ xcen: [4.2140429078454309]
+ xsigma: [2.0780002458907352]
+
@mdboom Collaborator
mdboom added a note

If I understand correctly, this is two examples, each with a few parts. If that's the case, each example should have a section header, rather than a bullet point, since it's not clear when this gets rendered that these are all independent examples.

@nden Collaborator
nden added a note

I moved stuff around and made section headers.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@mdboom
Collaborator

@nden: I've gone through and made some stylistic comments. I just commented on the first instance of most things I came across. Also, you probably want to run everything through PEP8 to find some more formatting glitches.

I'll probably have some comments later about how all the pieces fit together once I've digested it in my mind, but I think this looks really good. Thanks.

docs/models/index.rst
((85 lines not shown))
+ 3.97988800e+03 4.17389200e+03 4.37406400e+03 4.58050000e+03
+ 4.79329600e+03 5.01254800e+03 5.23835200e+03 5.47080400e+03
+ 5.71000000e+03 5.95603600e+03 6.20900800e+03 6.46901200e+03
+ 6.73614400e+03 7.01050000e+03 7.29217600e+03 7.58126800e+03
+ 7.87787200e+03 8.18208400e+03 8.49400000e+03 8.81371600e+03
+ 9.14132800e+03 9.47693200e+03 9.82062400e+03 1.01725000e+04
+ 1.05326560e+04 1.09011880e+04 1.12781920e+04 1.16637640e+04
+ 1.20580000e+04 1.24609960e+04 1.28728480e+04 1.32936520e+04
+ 1.37235040e+04 1.41625000e+04 1.46107360e+04 1.50683080e+04
+ 1.55353120e+04 1.60118440e+04]
+
+
+Add some noise
+
+>>> n=np.random.randn(90)
+>>> ny = y + n
@mdboom Collaborator
mdboom added a note

I wasn't clear in my original comment -- the thing I found puzzling in this example is that ny isn't subsequently used later on -- it's overwritten in a later example, but it's not clear what the purpose of this line is, even if it's clear what it does.

@nden Collaborator
nden added a note

Ah, duh, yes, I want to fit the noisy data!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@wkerzendorf wkerzendorf referenced this pull request in astropy/specutils
Closed

Spectrum1d 2nd draft #6

@hamogu hamogu referenced this pull request from a commit in hamogu/specutils
@hamogu hamogu Three Readers for Spectrum1D f14ec0d
@nden
Collaborator

@mdboom I think I addressed all comments you made. Thanks for reviewing this. If you get a chance please look at models.models._ParameterProperty since I haven't used descriptors before.

@wkerzendorf
Collaborator

@nden Hey Nadia, this is very close to merge, right? I think there's lots in there which we could use in Spectrum1D (specutils). Great job!!

@nden
Collaborator

@wkerzendorf Hi Wolfgang! Have you tried it for specutils? So far the comments have been mostly on style. It'd be nice if it gets some real testing before merging.

@keflavich
Collaborator
docs/models/models.rst
((96 lines not shown))
+ [ 0., 1., 2., 3., 4.],
+ [ 0., 1., 2., 3., 4.]])
+>>> print y.shape
+(10,5)
+
+- Create and evaluate a parallel composite model
+
+>>> x = np.arange(1,10,.1)
+>>> p1 = models.Poly1DModel(1)
+>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
+>>> pcomptr = models.PCompositeModel([g1, p1])
+>>> y = pcomptr(x)
+
+This is equivalent to applying the two models in parallel:
+
+>>> y = x + (g1(x) - x) + (p1(x) - x)
@taldcroft Collaborator

We discussed this before, and I still don't know exactly what PCompositeModel does. The explanation in email was:

  • PCompositeModel applies each model to the input data, takes the difference between the result and the input and adds the sum of all deltas to the input. This is a common operation with distortion corrections and the above equation is meant to represent these operations, not to be mathematically elegant.

Would it be possible to write down a mathematically (or programmatically) correct statement of what PCompositeModel does? When you write y = x + (g1(x) - x) ... my parser just shuts down because y and x may have entirely different units so the subtraction / addition operation is not meaningful. So what I want to know is what is the model output y as a function of x given the two functions g1 and p1. Assume that the data being fit with the model is called y_data.

@nden Collaborator
nden added a note

The initial context of this package was WCS support and this particular model may not be applicable to astrophysical models. The use case it covers is a detector with a distortion described with more than one model, e.g. one component is a polynomial model and another component is a filter dependent displacement represented as a lookup table. So 'x' and 'y' are coordinates in slightly different coordinate systems and in this case the above statement represents the mathematical operation. As an example, in astropy.wcs this is represented through sip_pix2foc() and p4_pix2foc() methods. The disadvantage there is that this is hardwired and other models cannot be used. The aim here is to make it more general and allow other types of mathematical (unitless) models.
When units are added to parameters I think there will be an optional None value for parameters of mathematical models.
But I see your point about astrophysical models. May be there needs to be a check for the type of models added in this way and a rejection mechanism if any of them are astrophysical models.
Or would it be better to rename this class to something which indicates its intended use?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((114 lines not shown))
+
+>>> x, y = np.mgrid[:10, :10]
+>>> off = models.ShiftModel(-3.2)
+>>> poly2 = models.Poly2DModel(2)
+>>> scomptr = models.SCompositeModel([off, poly2], inmap=[['x'], ['x', 'y']], outmap=[['x'], ['z']])
+
+The above composite transform will apply an inplace shift to x, followed by a 2D
+polynomial and will save the result in an array, labeled 'z'.
+To evaluate this model use a LabeledInput object
+
+>>> ado = models.LabeledInput([x, y], ['x', 'y'])
+>>> result = scomptr(ado)
+
+The output is also a LabeledInput object and the result is stored in label 'z'.
+
+>>> print result
@taldcroft Collaborator

You could probably get the idea across with a 5x5 model instead of 10x10, making the print statements a bit more compact and readable.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((104 lines not shown))
+>>> p1 = models.Poly1DModel(1)
+>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
+>>> pcomptr = models.PCompositeModel([g1, p1])
+>>> y = pcomptr(x)
+
+This is equivalent to applying the two models in parallel:
+
+>>> y = x + (g1(x) - x) + (p1(x) - x)
+
+In more complex cases the input and output may be mapped to transformations:
+
+>>> x, y = np.mgrid[:10, :10]
+>>> off = models.ShiftModel(-3.2)
+>>> poly2 = models.Poly2DModel(2)
+>>> scomptr = models.SCompositeModel([off, poly2], inmap=[['x'], ['x', 'y']], outmap=[['x'], ['z']])
+
@taldcroft Collaborator

So within this current framework would it be possible to eventually have something like SumCompositeModel and ProductCompositeModel models? In these cases the model outputs would be g1(x) + p1(x) and g1(x) * p1(x). This would make the model language useful for astrophysical data modeling.

@nden Collaborator
nden added a note

I have a toy implementation of SumCompositeModel. I think it would be straightforward to create these for the purpose of evaluating them.
Do you also expect to fit instances of SumComposite and ProductComposite models?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@wkerzendorf
Collaborator

@nden: Well it would be very useful to use the Chebychev and Legendre models. @hamogu has written code that imports quite a bit of that stuff from fits files, but for now it just uses this in functions. I think using these models is the right way forward and also a quick way forward.

When looking over the code, we should use the excellent np.polynomial.Legendre, np.polynomial.Polynomial, and np.polynomial.Chebychev classes instead of doing the calculations ourselves.

@nden
Collaborator

@wkerzendorf Can you elaborate a bit on why we should use the numpy polynomial classes instead of these? The advantage of implementing them here is to use this interface and framework. Whenever possible for low level computation, I agree we should use numpy I just don't consider this low level.

@wkerzendorf
Collaborator

@nden: well I think we should use this interface and framework - but when evaluating the Chebychev and Legendre polynomials use numpy's already built-in functions. They are much more tested than ours.

@jehturner
Collaborator

It's very good to see this is still active and getting positive feedback. I noticed a while ago that this module requires Python 2.7, unlike the rest of AstroPy (otherwise I'd probably have more feedback by now, as my laptop disk is full and our 2.7 build is perpetually almost stable...) -- will that be the final version requirement?

@nden
Collaborator

@jehturner Ah, thanks for reminding me. I think the only dependence on 2.7 is the OrderedDict but Astropy has a port back to 2.6. I'll switch the code to that and let you know.

@nden
Collaborator

@jehturner The dependence on OrderedDict from python 2.7 is eliminated now.

@jehturner
Collaborator
@nden
Collaborator

@wkerzendorf I agree in principle with you but am a bit hesitant to switch to numpy's polynomial functions (we are not talking about the classes, right?) right now for several reasons. The 2D Chebyshev and Legendre polynomials in numpy were released for the first time in numpy v 1.7. The implementation in this package, while going into astropy now, has been used for some time locally so it has some real testing. I see that the 2D algorithms are different (which is not surprising), but I'll have to install numpy 1.7 and test performance. As for the 1D algorithms, most people will agree that the Clenshaw algorithm is best and this is what both packages use. In fact they look awfully similar ;).
Of course there's the question about polynomial support in previous versions of numpy.

@perrygreenfield
@wkerzendorf
Collaborator

@nden and @perrygreenfield I agree that this should not keep the package from merging. On a related note, I also think that we should merge this PR sooner rather than later and "fix" some of the problems later (as some problems will only show up when we use the models package).

I'm sure that the algorithms are similar or even the same as in numpy - however, I think that in general, we should avoid maintaining additional code that someone else is already maintaining.

astropy/models/fitting.py
((515 lines not shown))
+ """
+ Parameters
+ ----------
+ model : a fittable :class: `models.ParametricModel`
+ model to fit to data
+
+ Raises
+ ------
+ ModelLinearityError
+ A linear model is passed to a nonlinear fitter
+
+ """
+ super(SLSQPFitter, self).__init__(model)
+ if self.model.linear:
+ raise ModelLinearityError('Model is linear in parameters, '
+ 'consider using linear fitting methods.')
@taldcroft Collaborator

Shouldn't this be a warning instead of raising an exception? Maybe the user has a reason to use this minimization method, even if the model is linear.

@taldcroft Collaborator

All other cases of raising an exception where a warning is sufficient should be changed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@taldcroft
Collaborator

I saw a nice talk at PyCon by Alex Martelli "Good enough" is good enough!. With that (and the one comment about exception vs. warnings), I'm good with merging this.

@nden
Collaborator

@taldcroft Thanks for the review and all the suggestions. Now I'll have to see this talk :)
I was experimenting with fitting composite models and think this can be supported with the current framework.
A number of details have to be decided on, like parameter name conflicts and others, but it is possible.

@astrofrog
Owner

@nden - I've finally got a chance to review this, so a few comments will follow :-) First, when I try importing astropy.models, I get:

In [1]: from astropy.models import models, fitting
ERROR: ImportError: cannot import name OrderedDict [IPython.core.interactiveshell]
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-22b53997cd5b> in <module>()
----> 1 from astropy.models import models, fitting

/Users/tom/Library/Python/3.2/lib/python/site-packages/astropy/models/__init__.py in <module>()
      1 # Licensed under a 3-clause BSD style license - see LICENSE.rst
----> 2 from . import models
      3 from . import fitting
      4 from . import parameters
      5 from . import projections

/Users/tom/Library/Python/3.2/lib/python/site-packages/astropy/models/models.py in <module>()
     43 import operator
     44 import abc
---> 45 from ..utils.compat import OrderedDict
     46 import numpy as np
     47 from . import parameters

ImportError: cannot import name OrderedDict

because I think the import should be

from ..utils.compat.odict import OrderedDict

I'm also seeing a more minor issue when building Astropy, but this should probably be fixed:

Can't parse docstring in build/lib.macosx-10.6-x86_64-3.2/astropy/models/models.py line 1822: TokenError: ('EOF in multi-line statement', (2, 0))
Can't parse docstring in build/lib.macosx-10.6-x86_64-3.2/astropy/models/models.py line 1823: TokenError: ('EOF in multi-line statement', (2, 0))

However, I can't figure out at first glance what is wrong with the docstring.

docs/models/index.rst
((27 lines not shown))
+`~astropy.models.SCompositeModel`, so that the output values of one model are
+used as input values to another. It is also possible to form a new model by
+combining models in parallel (each model is evaluated separately with the
+original input and the deltas are summed), `~astropy.models.PCompositeModel`.
+Since models may have multiple input values, machinery is provided that allows
+assigning outputs from one model into the appropriate input of another in a
+flexible way, `~astropy.models.models.LabeledInput`. Finally, it is permitted
+to combine any number of models using all of these mechanisms simultaneously.
+A composite model can be used to make further composite models. The goal
+is to eventually provide a rich toolset of models and fitters such that most users will
+not need to define new model classes, nor special purpose fitting
+routines (but not making that hard to do if it is necessary).
+
+While this system is initially being developed to support WCS models,
+its generality extends well beyond WCS cases, and thus the initial
+examples are not specifically geared to WCS problems.
@astrofrog Owner

I think this introduction is too long and gives too many details that users shouldn't have to worry about - in particular, we don't need to mention that this is being developed for WCS, because it looks like it can already be used with common generic models. This section should also maybe include a short paragraph on what this gives over e.g. curve_fit (obviously astropy.models is much for flexible and powerful, but we should mention this).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((80 lines not shown))
+>>> ny = y + n
+
+Fit a Chebyshev polynomial to the data
+
+>>> ch2 = models.ChebyshevModel(3)
+>>> chfit = fitting.LinearLSQFitter(ch2)
+>>> chfit(x, ny)
+>>> ch2.parameters
+[0.957, 1.931, 3.029, 4.305]
+
+.. figure:: images/cheb_fit.png
+ :scale: 25 %
+
+Fit a data set with a gaussian model.
+
+>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
@astrofrog Owner

I think that the gaussian example, being simpler, should come first.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@astrofrog astrofrog commented on the diff
docs/models/index.rst
((122 lines not shown))
+
+Fit a 2D polynomial to the data
+
+>>> p2 = models.Poly2DModel(2)
+>>> print p2
+Model: Poly2DModel
+Dim: 2
+Degree: 2
+Parameter sets: 1
+Parameters:
+ c0_0: [0.0]
+ c1_0: [0.0]
+ c2_0: [0.0]
+ c0_1: [0.0]
+ c0_2: [0.0]
+ c1_1: [0.0]
@astrofrog Owner

It might be worth showing users how to figure out what these parameters mean.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((127 lines not shown))
+Model: Poly2DModel
+Dim: 2
+Degree: 2
+Parameter sets: 1
+Parameters:
+ c0_0: [0.0]
+ c1_0: [0.0]
+ c2_0: [0.0]
+ c0_1: [0.0]
+ c0_2: [0.0]
+ c1_1: [0.0]
+>>>pfit = fitting.LinearLSQFitter(p2)
+>>>n = np.random.randn(100)
+>>>n.shape = (10, 10)
+>>>pfit(x, y, z+n)
+>>>fitter = fitting.LinearLSQFitter(self.model)
@astrofrog Owner

fitter doesn't get used - should this line be removed? Also, can you add a space after >>> on these few lines?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((129 lines not shown))
+Degree: 2
+Parameter sets: 1
+Parameters:
+ c0_0: [0.0]
+ c1_0: [0.0]
+ c2_0: [0.0]
+ c0_1: [0.0]
+ c0_2: [0.0]
+ c1_1: [0.0]
+>>>pfit = fitting.LinearLSQFitter(p2)
+>>>n = np.random.randn(100)
+>>>n.shape = (10, 10)
+>>>pfit(x, y, z+n)
+>>>fitter = fitting.LinearLSQFitter(self.model)
+>>> p2.parameters
+[0.6354845, 2.016544, 3.0035796, 4.0907439, 4.989999, 6.000127]
@astrofrog Owner

I think that for this example, it would be really cool to use real data and fit it to model the sky background (with plots to demonstrate). I can provide such an example as a PR once this is merged. SImilarly, we could fit a 1D gaussian above to real spectral data (I can also provide).

@keflavich Collaborator

If you're looking for inspiration or templates, you can look here:
http://pyspeckit.bitbucket.org/html/sphinx/examples.html, e.g. https://github.com/keflavich/pyspeckit/blob/master/doc/example_hcop.rst. Data there is free to use if you want.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
@@ -0,0 +1,146 @@
+.. _models:
+
+******
+Models
+******
+
+The base class of all models is `~astropy.models.models.Model`, however fittable
+models should subclass `~astropy.models.models.ParametricModel`. Parametric
+models can be linear or nonlinear in a regression analysis sense.
+
+All models provide a `__call__` method which performs the transformation in a
+purely mathematical way, i.e. the models are unitless. In addition, when possible the
@astrofrog Owner

I think that talking about __call__ doesn't really fit in here, because users shouldn't care about it (however, you could show examples of a model being called). The only place it really makes sense to mention hidden methods is when discussing creating custom models.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/parameters.rst
((4 lines not shown))
+Parameters
+**********
+
+Parameters are used in three different contexts within this package:
+communicating with fitters, model evaluation and getting values to/from users.
+
+Models maintain a list of parameter names, `parnames`. Single parameters are list-like
+objects, instances of `~astropy.models.parameters._Parameter`. Simple mathematical operations can be
+performed with them. The preferred way for users to interact with models is through
+individual parameters.
+
+The goal of this package is, when possible, to allow simultaneous model evaluation
+and fitting with multiple parameter sets. Because of this, all models have a `psets`
+attribute, an array of shape `(len(parnames), paramdim)`, where `paramdim` is the number of
+parameter sets. Typically the array is of type float but can become an object array in
+some cases. `Psets` is used for model evaluation.
@astrofrog Owner

I think Psets should be all lowercase

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/parameters.rst
((30 lines not shown))
+parameters and cannot be updated directly.
+
+Parameters Examples
+-------------------
+
+- Polynomial models are created by default with all coefficients set to 0.
+
+>>> p1=models.Poly1DModel(degree=4)
+>>> p1.parnames
+['c0', 'c1', 'c2', 'c3', 'c4']
+>>> p1.parameters
+[0.0, 0.0, 0.0, 0.0, 0.0]
+
+- Coefficients can be set using the **`parameters`** attribute
+
+>>> p1.parameters = [0, 1, 2, 3, 4]
@astrofrog Owner

This doesn't work for me:

In [11]: p1.parameters = [0,1,2,3,4]
ERROR: InputParameterError [astropy.models.parameters]
---------------------------------------------------------------------------
InputParameterError                       Traceback (most recent call last)
<ipython-input-11-b31655922c88> in <module>()
----> 1 p1.parameters = [0,1,2,3,4]

/Users/tom/Library/Python/3.2/lib/python/site-packages/astropy/models/models.py in parameters(self, value)
    390         elif isinstance(value, (list, np.ndarray)):
    391             _val = parameters._tofloat(value)
--> 392             if self._parameters._is_same_length(_val):
    393                 self._parameters._changed = True
    394                 self._parameters[:] = _val

/Users/tom/Library/Python/3.2/lib/python/site-packages/astropy/models/parameters.py in _is_same_length(self, newpars)
    291 
    292         """
--> 293         parsize = _tofloat(newpars).size
    294         return parsize == self.__len__()
    295 

/Users/tom/Library/Python/3.2/lib/python/site-packages/astropy/models/parameters.py in _tofloat(value)
     35     else:
     36         raise InputParameterError(
---> 37             "Don't know how to convert parameter to float")
     38     return _value
     39 

InputParameterError: Don't know how to convert parameter to float
@nden Collaborator
nden added a note

@astrofrog I can't reproduce this. What versions of numpy and python are you using?

@astrofrog Owner

Huh, I don't see this anymore - ignore me!

@nden Collaborator
nden added a note

@astrofrog No, it was real. There was a problem with python 3 - some of the functions in the operator module are moved to the numbers module in python 3.
Anyway, I think I addressed all your comments. Previously you suggested to move the WCS related models into a wcs subpackage. I don't always see a clear distinction - some of the general models can be used for WCS purposes. Instead, how about pulling all base classes and classes which set up the framework in general into a separate module, "core.py", and keeping the rest of the models in models.py. In the future, when we have more models, we can rearrange the modules again.

I had very good comments and the documentation improved significantly. But I really think the only way we can know if this is useful is to start using it with real data. I'm sure this will bring also other documentation shortcomings.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
@@ -0,0 +1,102 @@
+.. _fitting:
+
+*******
+Fitting
+*******
+
+This module provides wrappers, called Fitters, around some Numpy and Scipy
+fitting functions. All Fitters take an instance of
+`models.ParametricModel` as input and define a `__call__` method
+which fits the model to the data and changes `model.parameters
@astrofrog Owner

I think that __call__ shouldn't be mentioned here (just say that they can be used as a function or something like that)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@astrofrog
Owner

@nden - this is great, and will be a very useful module! (I can think of many applications in my work)

I left some detailed comments above, mostly about documentation. I think that we need to be careful to separate documentation for users who just want to use it, and implementation details that only matter for developers or users who want to extend the functionality by creating their own models. I think that the 'getting started' examples should, if possible, use real astronomy data, but I've volunteered to provide this once the PR is merged. Also, as I mentioned above, I don't think users need to know that the module was originally intended for WCS.

Once the issues mentioned above are fixed, I think that we should consider merging soon so that we can start using this, as it will allow us to better determine if anything in the interface needs to be improved.

@keflavich keflavich referenced this pull request in pyspeckit/pyspeckit
Open

Incorporate astropy.models #1

@nden nden commented on the diff
astropy/models/tests/test_projections.py
((5 lines not shown))
+import os.path
+import pytest
+import numpy as np
+from numpy.testing import utils
+from .. import projections
+from ...io import fits
+from ... import wcs
+from ...utils.data import get_pkg_data_filename
+
+class TestProjections(object):
+ """
+ Test composite models evaluation in series
+ """
+ def setup_class(self):
+ self.w = wcs.WCS()
+ test_file = get_pkg_data_filename(os.path.join('data', '1904-66_AZP.fits'))
@nden Collaborator
nden added a note

@mdboom Mike, can you tell why this fails in Travis (works on my desktop). Instead of finding the local file, it tries to download it. (I must be blind ...)

@mdboom Collaborator
mdboom added a note

This is a wild guess, but perhaps that data file isn't getting installed? I don't see a setup_package.py file here which would normally be required to specify which data files get installed.

@mdboom Collaborator
mdboom added a note

Having just checked this out, indeed that's what it is. Look at the wcs module for an example of how to install data. It's also documented here: https://astropy.readthedocs.org/en/v0.2/development/building_packaging.html

@nden Collaborator
nden added a note

I missed that. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@astrofrog
Owner

Regarding the current Travis failures, I think that any test requiring scipy should be skipped if scipy is not present (see astropy.cosmology for an example of this). Basically, just make sure that there are no failures if scipy is not installed.

@nden
Collaborator

@astrofrog Tests pass now with one exception - a sphinx build fails with a KeyError in the threading module.

docs/models/index.rst
((28 lines not shown))
+===============
+
+All examples assume the following modules have been imported
+
+>>> import numpy as np
+>>> from astropy.models import models, fitting
+
+Working with 1D models
+======================
+
+Fit a data set with a gaussian model.
+
+>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
+>>> g1
+<Gauss1DModel(amplitude= [10.0],xcen= [4.2000000000000002],xsigma= [2.1000000000000001],paramdim=1)>
+>>> y = g1(x)
@cdeil Collaborator
cdeil added a note

x is not defined here.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((27 lines not shown))
+Getting Started
+===============
+
+All examples assume the following modules have been imported
+
+>>> import numpy as np
+>>> from astropy.models import models, fitting
+
+Working with 1D models
+======================
+
+Fit a data set with a gaussian model.
+
+>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
+>>> g1
+<Gauss1DModel(amplitude= [10.0],xcen= [4.2000000000000002],xsigma= [2.1000000000000001],paramdim=1)>
@cdeil Collaborator
cdeil added a note

Super minor comment: typically in Python there should be no space between parameter=value and there should be a space after commas:

<Gauss1DModel(amplitude=[10.0], xcen=[4.2000000000000002], xsigma=[2.1000000000000001], paramdim=1)>

If you change the code to print it this way it should be updated everywhere in the docs.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((66 lines not shown))
+[1.0, 2.0, 3.0, 4.0]
+>>> print ch1
+Model: ChebyshevModel
+Dim: 1
+Degree: 3
+Parameter sets: 1
+Parameters:
+ c0: [1.0]
+ c1: [2.0]
+ c2: [3.0]
+ c3: [4.0]
+>>> y = ch1(x)
+
+Add some noise
+
+>>> n=np.random.randn(90)
@cdeil Collaborator
cdeil added a note

There should be spaces around operators to be PEP8 compliant:

n = np.random.randn(90)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((105 lines not shown))
+Fit a 2D polynomial to the data
+
+>>> p2 = models.Poly2DModel(2)
+>>> print p2
+Model: Poly2DModel
+Dim: 2
+Degree: 2
+Parameter sets: 1
+Parameters:
+ c0_0: [0.0]
+ c1_0: [0.0]
+ c2_0: [0.0]
+ c0_1: [0.0]
+ c0_2: [0.0]
+ c1_1: [0.0]
+>>>pfit = fitting.LinearLSQFitter(p2)
@cdeil Collaborator
cdeil added a note

There should be a space after >>> on this and the following lines:

>>> pfit = fitting.LinearLSQFitter(p2)

I'll shut up now about whitespace ... if you care you can look through the docs.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((77 lines not shown))
+>>> y = ch1(x)
+
+Add some noise
+
+>>> n=np.random.randn(90)
+>>> ny = y + n
+
+Fit a Chebyshev polynomial to the data
+
+>>> ch2 = models.ChebyshevModel(3)
+>>> chfit = fitting.LinearLSQFitter(ch2)
+>>> chfit(x, ny)
+>>> ch2.parameters
+[0.957, 1.931, 3.029, 4.305]
+
+.. figure:: images/cheb_fit.png
@cdeil Collaborator
cdeil added a note

Maybe draw the data as points, i.e. not connected by lines?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/algorithms.rst
@@ -0,0 +1,56 @@
+.. _algorithms:
+
+**********
+Algorithms
+**********
+
+Univariate Polynomial Evaluation
+################################
+
+* The evaluation of 1D polynomial uses Horner's algorithm.
@cdeil Collaborator
cdeil added a note

polynomial -> polynomials

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/algorithms.rst
@@ -0,0 +1,56 @@
+.. _algorithms:
+
+**********
+Algorithms
+**********
+
+Univariate Polynomial Evaluation
+################################
+
+* The evaluation of 1D polynomial uses Horner's algorithm.
+
+* The evaluation of 1D Chebyshev and Legendre polynomials uses Clenshaw algorithm.
@cdeil Collaborator
cdeil added a note

"uses Clenshaw" -> "uses the Clenshaw"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/algorithms.rst
((37 lines not shown))
+ 2. Set :math:`r_0` to :math:`c_{\alpha(0)}`, where c is a list of coeeficients for each multiindex in inverse lexical order.
+
+ 3. For each monomial, n, in the polynomial:
+
+ - determine :math:`k = max \{1 \leq j \leq di: \alpha(n)_j \neq \alpha(n-1)_j\}`
+
+ - Set :math:`r_k := l_k(x)* (r_0 + r_1 + \dots + r_k)`
+
+ - Set :math:`r_0 = c_{\alpha(n)}, r_1 = \dots r_{k-1} = 0.`
+
+ 4. return :math:`r_0 + \dots + r_{di}`
+
+* The evaluation of multivariate Chebyshev and Legendre polynomials uses a
+ variation of the above Horner's scheme, in which every Legendre or Chebyshev
+ function is considered a separate variable. In this case the length
+ of the :math:`\alpha` indices tuple is equal to the number of functions in x plus teh number of functions in y. In addition the Chebyshev and Legendre functions are cached for efficiency.
@cdeil Collaborator
cdeil added a note

"teh" -> "the"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
@@ -0,0 +1,101 @@
+.. _fitting:
+
+*******
+Fitting
+*******
+
+This module provides wrappers, called Fitters, around some Numpy and Scipy
+fitting functions. All Fitters can be called as funcitons. They take an instance of
@cdeil Collaborator
cdeil added a note

"funcitons" -> "functions"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
@@ -0,0 +1,101 @@
+.. _fitting:
+
+*******
+Fitting
+*******
+
+This module provides wrappers, called Fitters, around some Numpy and Scipy
+fitting functions. All Fitters can be called as funcitons. They take an instance of
+`models.ParametricModel` as input and modily `model.parameters`
@cdeil Collaborator
cdeil added a note

"modily model.parameters" -> "modify the model.parameters"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
((69 lines not shown))
+>>> y = p1(x)
+>>> p1.c1.fixed = True
+>>> pfit=fitting.LinearLSQFitter(p1)
+>>> pfit(x, y)
+>>> p1.psets
+array([[ 5.50225913, 5.50225913],
+ [ 2. , 2. ],
+ [ 3.17551299, 3.17551299]])
+
+
+- Parameters can be tied. This can be done in two ways:
+
+>>> def tiedfunc(g1):
+ ... xcen = 3*g1.xsigma[0]
+ ... return xcen
+>>> g1 = models.Gauss1D(amplitude=10., xcen=3, xsigma=.5, tied={'xcen':tiedfunc})
@cdeil Collaborator
cdeil added a note

NameError: Gauss1D -> Gauss1DModel

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
((73 lines not shown))
+>>> p1.psets
+array([[ 5.50225913, 5.50225913],
+ [ 2. , 2. ],
+ [ 3.17551299, 3.17551299]])
+
+
+- Parameters can be tied. This can be done in two ways:
+
+>>> def tiedfunc(g1):
+ ... xcen = 3*g1.xsigma[0]
+ ... return xcen
+>>> g1 = models.Gauss1D(amplitude=10., xcen=3, xsigma=.5, tied={'xcen':tiedfunc})
+
+or
+
+>>> g1 = models.Gauss1D(amplitude=10., xcen=3, xsigma=.5)
@cdeil Collaborator
cdeil added a note

NameError: Gauss1D -> Gauss1DModel

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
((80 lines not shown))
+
+>>> def tiedfunc(g1):
+ ... xcen = 3*g1.xsigma[0]
+ ... return xcen
+>>> g1 = models.Gauss1D(amplitude=10., xcen=3, xsigma=.5, tied={'xcen':tiedfunc})
+
+or
+
+>>> g1 = models.Gauss1D(amplitude=10., xcen=3, xsigma=.5)
+>>> g1.xcen.tied = tiedfunc
+>>> gfit = fitting.NonLinearLSQFitter(g1)
+
+
+- Print a list of available fitting constraints
+
+>>> fitting.Constraints.fitters
@cdeil Collaborator
cdeil added a note
AttributeError: 'module' object has no attribute 'Constraints'

Add import astropy.models.constraints?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/fitting.rst
((64 lines not shown))
+>>> p1.parameters=[1,1,2,2,3,3]
+>>> p1.psets
+array([[ 1., 1.],
+ [ 2., 2.],
+ [ 3., 3.]])
+>>> y = p1(x)
+>>> p1.c1.fixed = True
+>>> pfit=fitting.LinearLSQFitter(p1)
+>>> pfit(x, y)
+>>> p1.psets
+array([[ 5.50225913, 5.50225913],
+ [ 2. , 2. ],
+ [ 3.17551299, 3.17551299]])
+
+
+- Parameters can be tied. This can be done in two ways:
@cdeil Collaborator
cdeil added a note

I'm not familiar with that term. Can you give a short definition what a "tied parameter" is here?

@cdeil Collaborator
cdeil added a note

It's clear to me now ... still, maybe you add a short explanation like tied, i.e. forced to always have the same value ?

@eteq Owner
eteq added a note

@nden - I agree with @cdeil - can you just add a sentence explaining what "tied" means?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((5 lines not shown))
+******
+
+The base class of all models is `~astropy.models.models.Model`, however fittable
+models should subclass `~astropy.models.models.ParametricModel`. Parametric
+models can be linear or nonlinear in a regression analysis sense.
+
+To evaluate a model, it is called like a function. When possible the
+transformation is done using multiple parameter sets, `psets`.
+The number of parameter sets is stored in an attribute `paramdim`.
+
+Parametric models also store a flat list of all parameters as an instance of
+`~astropy.models.parameters.Parameters`. When fitting, this list-like object is
+modified by a subclass of `~astropy.fitting.fitting.Fitter`. When fitting nonlinear models,
+the values of the parameters are used as initial guesses by the fitting class.
+
+Models have dimensions, `ndim` attribute, which show how many coordinates the
@cdeil Collaborator
cdeil added a note

"have dimensions" -> "have a dimensions"
"show" -> "shows"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((10 lines not shown))
+
+To evaluate a model, it is called like a function. When possible the
+transformation is done using multiple parameter sets, `psets`.
+The number of parameter sets is stored in an attribute `paramdim`.
+
+Parametric models also store a flat list of all parameters as an instance of
+`~astropy.models.parameters.Parameters`. When fitting, this list-like object is
+modified by a subclass of `~astropy.fitting.fitting.Fitter`. When fitting nonlinear models,
+the values of the parameters are used as initial guesses by the fitting class.
+
+Models have dimensions, `ndim` attribute, which show how many coordinates the
+model expects as an input. All models expect coordinates as separate arguments.
+For example a 2D model expects x and y to be passed separately,
+e.g. as two arrays or two lists. When a model has multiple parameter sets and x, y are
+2D arrays, the model is evaluated with each of the parameter sets and the same x, y as
+input. The shape of the output array is (paramdim, xsh, ysh) where paramdim is the number
@cdeil Collaborator
cdeil added a note

Use more explicit, less cryptic x_shape and y_shape instead of xsh, ysh.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((90 lines not shown))
+ [ 0., 1., 2., 3., 4.],
+ [ 0., 1., 2., 3., 4.],
+ [ 0., 1., 2., 3., 4.],
+ [ 0., 1., 2., 3., 4.],
+ [ 0., 1., 2., 3., 4.],
+ [ 0., 1., 2., 3., 4.],
+ [ 0., 1., 2., 3., 4.]])
+>>> print y.shape
+(10,5)
+
+- Create and evaluate a parallel composite model
+
+>>> x = np.arange(1,10,.1)
+>>> p1 = models.Poly1DModel(1)
+>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
+>>> pcomptr = models.PCompositeModel([g1, p1])
@cdeil Collaborator
cdeil added a note

pcomptr is a bit cryptic. Maybe it means "parallel composite transformation"?
How about calling it composite_model or c1 instead?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((108 lines not shown))
+This is equivalent to applying the two models in parallel:
+
+>>> y = x + (g1(x) - x) + (p1(x) - x)
+
+In more complex cases the input and output may be mapped to transformations:
+
+>>> x, y = np.mgrid[:5, :5]
+>>> off = models.ShiftModel(-3.2)
+>>> poly2 = models.Poly2DModel(2)
+>>> scomptr = models.SCompositeModel([off, poly2], inmap=[['x'], ['x', 'y']], outmap=[['x'], ['z']])
+
+The above composite transform will apply an inplace shift to x, followed by a 2D
+polynomial and will save the result in an array, labeled 'z'.
+To evaluate this model use a LabeledInput object
+
+>>> ado = models.LabeledInput([x, y], ['x', 'y'])
@cdeil Collaborator
cdeil added a note

No idea what ado stands for. Use a more descriptive name in this example?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/parameters.rst
((13 lines not shown))
+individual parameters.
+
+The goal of this package is, when possible, to allow simultaneous model evaluation
+and fitting with multiple parameter sets. Because of this, all models have a `psets`
+attribute, an array of shape `(len(parnames), paramdim)`, where `paramdim` is the number of
+parameter sets. Typically the array is of type float but can become an object array in
+some cases. `psets` is used for model evaluation.
+
+In addition, all models maintain an attribute, `parameters`, an instance of
+`~astropy.parameters.Parameters`. This is a flat list of
+parameter values which fitters update. It serves as a communication tool between fitters
+and models.
+
+Individual parameters, `psets` and the flat list of parameter values are kept in sync.
+Single parameters are updated through properties. An update to a single parameter
+trigers an update to `psets` and `parameters`. Single parameters are updated
@cdeil Collaborator
cdeil added a note

trigers -> triggers

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/parameters.rst
((28 lines not shown))
+trigers an update to `psets` and `parameters`. Single parameters are updated
+after a change to `parameters`. `psets` are always constructed on demand from single
+parameters and cannot be updated directly.
+
+Parameters Examples
+-------------------
+
+- Polynomial models are created by default with all coefficients set to 0.
+
+>>> p1=models.Poly1DModel(degree=4)
+>>> p1.parnames
+['c0', 'c1', 'c2', 'c3', 'c4']
+>>> p1.parameters
+[0.0, 0.0, 0.0, 0.0, 0.0]
+
+- Coefficients can be set using the **`parameters`** attribute
@cdeil Collaborator
cdeil added a note

Maybe not make parameters bold here to be consistent with the rest of the docs?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/parameters.rst
((34 lines not shown))
+
+- Polynomial models are created by default with all coefficients set to 0.
+
+>>> p1=models.Poly1DModel(degree=4)
+>>> p1.parnames
+['c0', 'c1', 'c2', 'c3', 'c4']
+>>> p1.parameters
+[0.0, 0.0, 0.0, 0.0, 0.0]
+
+- Coefficients can be set using the **`parameters`** attribute
+
+>>> p1.parameters = [0, 1, 2, 3, 4]
+>>> p1.parameters
+[0.0, 1.0, 2.0, 3.0, 4.0]
+
+- It is possible to set the coefficients passing a dictionary to the `pars` keyword
@cdeil Collaborator
cdeil added a note

You are not showing the pars keyword in the following example.

I tried it and it doesn't work here:


In [111]: ch2 = models.ICheb2DModel(xdeg=2, ydeg=3, pars=coeff)
ERROR: AssertionError [astropy.models.models]
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-111-792e3098de53> in <module>()
----> 1 ch2 = models.ICheb2DModel(xdeg=2, ydeg=3, pars=coeff)

/Users/deil/code/astropy/astropy/models/models.pyc in __init__(self, xdeg, ydeg, xdomain, xwindow, ydomain, ywindow, paramdim, **pars)
   1074                                            xdomain=xdomain, ydomain=ydomain,
   1075                                            xwindow=xwindow, ywindow=ywindow,
-> 1076                                            paramdim=paramdim, **pars)
   1077 
   1078 

/Users/deil/code/astropy/astropy/models/models.pyc in __init__(self, xdeg, ydeg, xdomain, xwindow, ydomain, ywindow, paramdim, **pars)
    577                 print("Creating a model with {0} parameter sets\n".format(lenpars))
    578                 paramdim = lenpars
--> 579             self._validate_pars(**pars)
    580             self.set_coeff(pardim=paramdim, **pars)
    581         super(IModel, self).__init__(self.parnames, ndim=2, outdim=1,

/Users/deil/code/astropy/astropy/models/models.pyc in _validate_pars(self, **pars)
    617     def _validate_pars(self, **pars):
    618         numcoeff = self.get_numcoeff()
--> 619         assert(len(pars) == numcoeff)
    620 
    621     def _invlex(self):

AssertionError: 
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/parameters.rst
((48 lines not shown))
+
+- It is possible to set the coefficients passing a dictionary to the `pars` keyword
+
+>>> ch2 = models.ICheb2DModel(xdeg=2,ydeg=3, paramdim=2)
+>>> coeff = {}
+>>> for i, j in zip(ch2.parnames, range(len(ch2.parnames))):
+ coeff[i] = [j, j+10]
+>>> ch2 = models.ICheb2DModel(xdeg=2, ydeg=3, **coeff)
+>>> ch2.psets
+array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
+ [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]])
+
+
+- or directly, using keyword arguments
+
+>>> ch2 = models.ICheb2DModel(xdeg=2, ydeg=3, 'c0_0'=[0, 10], 'c0_1'=[3, 13],
@cdeil Collaborator
cdeil added a note

In [115]: ch2 = models.ICheb2DModel(xdeg=2, ydeg=3, 'c0_0'=[0, 10], 'c0_1'=[3, 13],'c0_2'=[6, 16], 'c0_3'=[9, 19],'c1_0'=[1, 11], 'c1_1'=[4, 14],'c1_2'=[7, 17], 'c1_3'=[10, 20,], 'c2_0'=[2, 12], 'c2_1'=[5, 15],'c2_2'=[8, 18], 'c2_3'=[11, 21])
  File "<ipython-input-115-c79fc0414ab1>", line 1
    ch2 = models.ICheb2DModel(xdeg=2, ydeg=3, 'c0_0'=[0, 10], 'c0_1'=[3, 13],'c0_2'=[6, 16], 'c0_3'=[9, 19],'c1_0'=[1, 11], 'c1_1'=[4, 14],'c1_2'=[7, 17], 'c1_3'=[10, 20,], 'c2_0'=[2, 12], 'c2_1'=[5, 15],'c2_2'=[8, 18], 'c2_3'=[11, 21])
SyntaxError: keyword can't be an expression
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@cdeil
Collaborator

@nden Thanks for writing this!
I'm looking forward to starting to use this and hope it gets merged soon ... sorry for commenting so late.

I tried out the tutorial and left some minor inline comments, a few of the examples don't work directly.

For the new.rst tutorial I think it would be better to (maybe additionally) give the reader a working example that she can simply copy & paste to try it out or as a template to get started on her own model. This doesn't fit too well with the narrative, maybe the example can be a separate .py file that you link to at the top?

Is it correct that integration of models over bins and convolution are not available at the moment? Will it be possible to integrate this with the current API?

docs/models/new.rst
((21 lines not shown))
+
+::
+
+ class Gauss1DModel(ParametricModel):
+ parnames = ['amplitude', 'xcen', 'xsigma']
+
+
+As a minimum the __init__ method takes all parameters and the number of parameter sets, `paramdim`:
+
+::
+
+ def __init__(self, amplitude, xcen, xsigma, paramdim=1)
+ self.linear = False
+ self.ndim = 1
+ self.outdim = 1
+ self._amplitude = parameters._Parameter(name='amplitude', val=amplitude, mclass=self, paramdim=paramdim)
@cdeil Collaborator
cdeil added a note

Rename _Parameter to Parameter?
A user-defined model is quite common and it shows up in the docs, so I would say it should be part of the public API and not internal as the underscore suggests.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/new.rst
((17 lines not shown))
+Pass the list of parameter names and the number of parameter sets to the base
+class. Note, that if the method which evaluates the model cannot work
+with multiple parameter sets, `paramdim` should not be given as an argument
+in the __init__ method. The default for `paramdim` is set in the base class to 1.
+
+::
+
+ class Gauss1DModel(ParametricModel):
+ parnames = ['amplitude', 'xcen', 'xsigma']
+
+
+As a minimum the __init__ method takes all parameters and the number of parameter sets, `paramdim`:
+
+::
+
+ def __init__(self, amplitude, xcen, xsigma, paramdim=1)
@cdeil Collaborator
cdeil added a note

Missing semicolon.
Also when copy & pasting the indentation is lost.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/new.rst
((58 lines not shown))
+ return params[0] * np.exp((-(1/(params[2]**2)) * (x-params[1])**2))
+
+
+The `deriv` method takes as input all coordinates as separate arguments.
+There is an option to compute numerical derivatives for nonlinear models
+in which case the `deriv` method should return None.
+
+Finally, the __call__ method takes input coordinates as separate arguments.
+It reformats them (if necessary) and calls the `eval` method to perform the
+model evaluation using model.psets as parameters.
+The reason there is a separate `eval` method is to allow fitters to call the `eval`
+method with different parameters which is necessary for fitting with constraints.
+
+::
+
+ def __call__(self, x):
@cdeil Collaborator
cdeil added a note

This looks like boilerplate code the user has to copy & paste for most models.
Can't this default implementation be in the ParametercModel?

@nden Collaborator
nden added a note

It is in Model but is an abstract method. It may be possible to move the implementation there as well but currently there are few models in the package and some models may need more complex call methods. We can always move it later.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/new.rst
((19 lines not shown))
+with multiple parameter sets, `paramdim` should not be given as an argument
+in the __init__ method. The default for `paramdim` is set in the base class to 1.
+
+::
+
+ class Gauss1DModel(ParametricModel):
+ parnames = ['amplitude', 'xcen', 'xsigma']
+
+
+As a minimum the __init__ method takes all parameters and the number of parameter sets, `paramdim`:
+
+::
+
+ def __init__(self, amplitude, xcen, xsigma, paramdim=1)
+ self.linear = False
+ self.ndim = 1
@cdeil Collaborator
cdeil added a note

I tried implementing my own model:

class Line(models.ParametricModel):
    parnames = ['slope', 'intercept']    
    def __init__(self, slope, intercept, paramdim=1):
        self.linear = True
        self.ndim = 1
        self.outdim = 1
        self._slope = parameters._Parameter(name='slope', val=slope, mclass=self, paramdim=paramdim)
        self._intercept = parameters._Parameter(name='intercept', val=intercept, mclass=self, paramdim=paramdim)
        models.ParametricModel.__init__(self, self.parnames, paramdim=paramdim)
    def eval(self, x, params):
        return params[0] * x + params[1]
    def __call__(self, x):
        x, format = _convert_input(x, self.paramdim)
        result = self.eval(x, self.psets)
        return _convert_output(result, format)

model = Line(slope=1, intercept=1)

I get an AttributeError: can't set attribute for self.ndim = 1 and if I comment out that line also for self.outdim = 1.
If I comment out both of these lines I get TypeError: __init__() takes at least 4 arguments (3 given) for this line: models.ParametricModel.__init__(self, self.parnames, paramdim=paramdim)

@nden Collaborator
nden added a note

@cdeil I realize now that some of the documentation did not keep up with changes in the code triggered by reviews, new.rst is one of them. There was a request to make ndim and outdim properties and I also moved them to the base class. So they are passed now to the base class constructor. Here's how your example looks now:

from astropy.models import models, parameters

class Line(models.ParametricModel):
parnames = ['slope', 'intercept']

def init(self, slope, intercept, paramdim=1):
self.linear = True
self.slope = parameters.Parameter(name='slope', val=slope, mclass=self, paramdim=paramdim)
self._intercept = parameters._Parameter(name='intercept', val=intercept, mclass=self, paramdim=paramdim)
models.ParametricModel.__init
(self, self.parnames, ndim=1, outdim=1, paramdim=paramdim)

def eval(self, x, params):
    return params[0] * x + params[1]
def __call__(self, x):
    x, format = models._convert_input(x, self.paramdim)
    result = self.eval(x, self.psets)
    return models._convert_output(result, format)

model = Line(slope=1, intercept=1)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@cdeil
Collaborator

Here's the notebook where I tried to make a user-defined model and fit statistic and the errors I ran into:
http://nbviewer.ipython.org/5335831

@nden
Collaborator

@cdeil Thanks for the detailed notes. I fixed all problems you pointed out. Still need to recreate the figures without lines but some on the mailing list actually offered to write examples with real data, so I'll wait on that. Also I'm not sure I'm convinced that _Parameter should be renamed to Parameter. It's true there will be many custom classes but only those who write them have to use this class so in that sense they are developers. What do other think?

I can' t view your notebook, I don't know why. Was there anything else there (other than the Line class)? If so could you send me the errors directly in an email, or possibly send me the notebook.

@cdeil
Collaborator

Hmm ... it's a public gist at https://gist.github.com/cdeil/5335831 and I can see the notebook fine ... maybe it was a glitch and you can try again?

I still can't get the user-defined model to work. Here's where I am now:

import numpy as np
from astropy.models import models, parameters, fitting

class Line(models.ParametricModel):
    parnames = ['slope', 'intercept']
    def __init__(self, slope, intercept, paramdim=1):
        self.linear = True
        self._slope = parameters._Parameter(name='slope', val=slope, mclass=self, paramdim=paramdim)
        self._intercept = parameters._Parameter(name='intercept', val=intercept, mclass=self, paramdim=paramdim)
        models.ParametricModel.__init__(self, self.parnames, ndim=1, outdim=1, paramdim=paramdim)
    def eval(self, x, params):
        return params[0] * x + params[1]
    def __call__(self, x):
        x, format = models._convert_input(x, self.paramdim)
        result = self.eval(x, self.psets)
        return models._convert_output(result, format)

xdata = np.array([0, 1, 2, 3, 4, 5])
ydata = np.array([1, 1, 5, 7, 8, 12])
sigma = np.array([1, 2, 1, 2, 1, 2], dtype=float)

model = Line(slope=1, intercept=1)
fitter = fitting.LinearLSQFitter(model)
fitter(xdata, ydata, weights=sigma ** (-1))

I get

AttributeError: 'Line' object has no attribute 'deriv'

and then if I add self.deriv = None in the constructor a similar error for domain and then for window.

I'll just wait for you to push the documentation update and then try it again.

@nden
Collaborator

@cdeil All fittable models should have a deriv method. The non-linear fitters have the option to estimate the derivatives if a method was not provided and deriv is None. The linear fitters expect a method which creates the Vandermond matrix. So in general, the developer of a class must provide a deriv method.

There are other issues and I am hoping to get an opinion on the best way to resolve them. The linear fitter is written in a way which assumes that some kind of a polynomial model is going to be fit. Since some polynomials are defined on specific domains and these are important when fitting, the polynomial classes have a 'domain' and a 'window' attributes. All polynomial classes have a base class PModel but I realize now that the domain and window attributes are in the individual subclasses (Poly1DModel, Cheb1DModel, ...) instead of PModel. There are two ways to fix this: either move domain and window to PModel or check for these attributes in the LinearLSQFitter and remap only if the attributes are present. I think the second is probably better. Any thoughts on that?

Also , there's a rank check in the LinearLSQFitter which assumes that model._order is present. I guess I would get rid of this.

In any case, here's a working example of line.py:

from astropy.models import models, parameters
import numpy as np

class Line(models.PModel):
    parnames = ['slope', 'intercept']

def init(self, slope, intercept, paramdim=1):
    self.linear = True 
    self._slope = parameters._Parameter(name='slope', val=slope, mclass=self, paramdim=paramdim)
    self._intercept = parameters._Parameter(name='intercept', val=intercept, mclass=self, paramdim=paramdim)
    models.ParametricModel.__init__(self, self.parnames, ndim=1, outdim=1, paramdim=paramdim)
    self.domain = [-1, 1]
    self.window = [-1, 1]
    self._order = 2

def eval(self, x, params):
    return params[0] * x + params[1]

def call(self, x):
    x, format = models._convert_input(x, self.paramdim)
    result = self.eval(x, self.psets)
    return models._convert_output(result, format)

def deriv(self, x):
    res = np.ones((len(x), 2))
    res[:, 1] = x.copy()
    return res
@wkerzendorf wkerzendorf referenced this pull request in sunpy/sunpy
Closed

Change PyFITS dependency to Astropy #425

@cdeil
Collaborator

@nden I have been playing some more with the latest version of astropy.models in this PR in this IPython notebook: gist html

With your helpful comments I was able to make the UserLineModel run (cells 12 and 13).
But the result I get isn't correct, so I am still doing something wrong.
I do get correct and matching results using the built-in astropy polynomial model (cell 11) and also another fitting package, iminuit (cells 10).

Another thing I tried is a chi2 fit with a non-linear model, namely a 1D Gauss. The result I get from astropy (cell 23) doesn't match the iminuit results (cell 22).

There's a few TODOs in the notebook of things I think should already be possible with astropy.models, e.g. I would like to be able to do a Poisson likelihood fit (see cells 26 to 30 for what I mean, which gives an example and reference results with iminuit).

@cdeil
Collaborator

Some more general feedback:

  • Add unit tests for user-defined models (linear and non-linear) and fit statistics (chi2 and e.g. cash (see cell 29 here).
  • new.rst should contain full working examples.
  • astropy.models.rotations should be described in the docs, at least mentioned for now or user's won't find it.
  • Currently some big pieces are still missing:
    • model integration over bins
    • model convolution
    • symmetric errors (covariance matrix as as inverse Hesse matrix)
    • asymmetric errors (profile likelihood) Would it make sense to have a fitting API document like there is for some other parts of astropy? This would help me a lot to better get the "big picture", i.e. see how realistic use cases might look like in the end. E.g. I have this one concrete suggestions, but I'm not sure if it makes sense:
      • The concepts of fit statistic and optimiser are combined in the Fitter class. I think it could be better to have fit statistic separate like e.g. sherpa or iminuit do. I can make a case for this, but basically I think it's more modular and will make user-defined fit statistics and fitters easier to implement.

I'll now browse through the code and leave some comments.

@cdeil
Collaborator

@nden I've made a pull request with some changes: nden#1

You should also run autopep8, there's quite a few PEP8 warnings:
https://gist.github.com/cdeil/5361924

It will change a lot of code though ... so most of the inline code comments on github will be lost.

Side remark: Can you put your code examples here into triple-backticks so that it is shows as properly-formatted and -indented code? See "fenced code blocks" here.

astropy/models/parameters.py
((95 lines not shown))
+ self._max = maxvalue
+
+ @property
+ def paramdim(self):
+ return self._paramdim
+
+ @paramdim.setter
+ def model(self, val):
+ self._paramdim = val
+
+ @property
+ def mclass(self):
+ return self._mclass
+
+ @mclass.setter
+ def model(self, val):
@cdeil Collaborator
cdeil added a note

Is it a problem that this is a duplicated signature, i.e. this class has two setters called model?

    @paramdim.setter
    def model(self, val):
        self._paramdim = val

    @mclass.setter
    def model(self, val):
        self._mclass = val
@astrofrog Owner

This is normal syntax, even though pyflakes doesn't seem to recognize it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/util.py
@@ -0,0 +1,56 @@
+"""
+This module provides utility functions for the models package
+"""
+from __future__ import division, print_function
+import numpy as np
+
+class InputParameterError(Exception):
@cdeil Collaborator
cdeil added a note

There's nothing astropy.model package specific about this InputParameterError.
Maybe see what other sub-packages use and unify?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/util.py
((5 lines not shown))
+import numpy as np
+
+class InputParameterError(Exception):
+ """
+ Called when there's a problem with input parameters.
+ """
+ def __init__(self, message):
+ self._message = message
+
+ def __str__(self):
+ return self._message
+
+__all__ = ['pmapdomain', 'comb', 'InputParameterError']
+
+
+def pmapdomain(oldx, domain, window):
@cdeil Collaborator
cdeil added a note

What does the p stand for. Maybe change to parameter_map_domain if that's what it means?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@cdeil
Collaborator

At the moment there's ChebyshevModel and ICheb2DModel, where I stands for IRAF.

Before I looked at the docstring I thought the I in ICheb2DModel means Integral, i.e. that this model integrates over bins.

If it is important to keep the reference to IRAF I would suggest renaming to IRAFChebyshev2DModel, but hopefully in a few years there won't be anyone left who knows what IRAF is :-), so how about just Chebyshev1DModel and Chebyshev2DModel, which would be consistent with e.g. Gauss1DModel and Gauss2DModel?

@cdeil cdeil commented on the diff
astropy/models/tests/test_projections.py
@@ -0,0 +1,85 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+Test sky projections defined in WCS Paper II
+"""
+from __future__ import division
+import os.path
+import numpy as np
+from numpy.testing import utils
+from .. import projections
+from ...io import fits
+from ... import wcs
+from ...utils.data import get_pkg_data_filename
+from ...tests.helper import pytest
+
+def b(s):
@cdeil Collaborator
cdeil added a note

Is this encoding really necessary?
I think it shouldn't be because you just use the public wcs package API?

@nden Collaborator
nden added a note

I think after #956 is merged this won't be necessary.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/util.py
((24 lines not shown))
+ Parameters
+ ----------
+ oldx : array
+ original coordinates
+ domain : list or tuple of length 2
+ function domain
+ window : list or tuple of length 2
+ range into which to map the domain
+ """
+ domain = np.array(domain, dtype=np.float64)
+ window = np.array(window, dtype=np.float64)
+ scl = (window[1]-window[0])/(domain[1]-domain[0])
+ off = (window[0]*domain[1] - window[1]*domain[0])/(domain[1]-domain[0])
+ return off + scl*oldx
+
+def comb(N, k):
@cdeil Collaborator
cdeil added a note

Maybe implement like this to be more readable?

from math import factorial

def comb(N, k):
    return factorial(N) / factorial(k) / factorial(N - k)

You chose this name because it is the same as scipy.misc.comb, right?
http://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.comb.html
While it's nice having the same name for the same thing, I think it would be better to give a longer, more descriptive name here,
when I'm somewhere else and I see comb in the code I have to go look at the documentation to see what it means.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/models.py
((1,302 lines not shown))
+ if not the Jacobian will be estimated.
+
+ Returns
+ -------
+ model : Gauss1DModel
+ 1D Gaussian
+ """
+ self._amplitude = parameters._Parameter('amplitude', amplitude, self, 1)
+ if xsigma is None and fwhm is None:
+ raise InputParameterError(
+ "Either fwhm or xsigma must be specified")
+ if xsigma is not None:
+ xsigmaval = xsigma
+ else:
+ try:
+ xsigmaval = 0.42466 * fwhm
@cdeil Collaborator
cdeil added a note

Is it possible to define this constant in astropy.constants?
Here's an analytical formula(to get correct results up to rounding error precision):

In [12]: import numpy as np
In [13]: sigma_to_fwhm = np.sqrt(8 * np.log(2))         
In [14]: fwhm_to_sigma = 1 / sigma_to_fwhm
In [15]: sigma_to_fwhm
Out[15]: 2.3548200450309493
In [16]: fwhm_to_sigma
Out[16]: 0.42466090014400953

I frequently use this constant in my own analysis code, so having it available in astropy.constants would be nice.
I can make a separate pull request if someone tells me where it would fit best.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
astropy/models/models.py
((1,269 lines not shown))
+ "Expected input arrays to have the same shape"
+ invcoeff = self.invlex_coeff()
+ x, fmt = _convert_input(x, self.paramdim)
+ y, fmt = _convert_input(y, self.paramdim)
+ result = self.imhorner(x, y, invcoeff)
+ return _convert_output(result, fmt)
+
+class Gauss1DModel(ParametricModel):
+ """
+
+ Implements 1D Gaussian model
+
+ """
+ parnames = ['amplitude', 'xcen', 'xsigma']
+ def __init__(self, amplitude, xcen, fwhm=None, xsigma=None,
+ fjac=None, **cons):
@cdeil Collaborator
cdeil added a note

Should the cons parameter be explained in the docstring here?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@nden
Collaborator

@cdeil Thanks, this is useful feedback. I'll go over your comments today but to start I looked at the gaussian fitting example.

I am not familiar with Minuit. I wonder if the difference comes from using different algorithms. The NonLinearLSQFitter uses the Levenberg-Marquardt algorithm of scipy.optimize.leastsq. Do you know what Minuit is using in this case? Below is a plot similar to your cell 20 but I added the results from Minuit and Models. Actually I ran models.fitter twice: first, providing a function to calculate the derivatives (the default for Gauss1DModel) and second, telling the fitter to estimate the derivatives by passing fjac='estimated' to Gauss1DModel. This plot shows several interesting things. First, in this simple case, the algorithm should not be confused and fitting with estimated and computed derivatives should give the same results (and in this case the fit with the estimated derivatives is the correct one). I must have introduced a bug when I moved constraints from fitting to models. Will look into this.
But more importantly how do we tell which result is "more correct"? It seems to me the best fit in this case is given by models with estimated derivatives. Any thoughts on that?

BTW, the NonLinearLSQFitter has a covar property which returns the covariance matrix (computed using the Hessian matrix).

gauss_fit_comparison

@keflavich
Collaborator

@nden, @cdeil - I'm a latecomer to this, but perhaps it would make more sense to benchmark against lmfit-py, which also uses scipy.optimize.leastsq (optionally) or mpfit.py http://cars9.uchicago.edu/software/python/mpfit.html, which implements the Levenberg-Marquardt algorithm directly.

In terms of assessing which is "more correct", I think it would be sensible to compare both the uncertainties and the fit parameters, so maybe plot the amplitude + error & width + error rather than the fits. In terms of absolute correctness, it would be good to compare to the parameter space generated by a MCMC method (e.g., pymc).

...I'll try to work on that today but no promises.

@cdeil
Collaborator

Here's a paper describing the methods used by MINUIT.
I never had the time to learn exactly how the minimisation step (called MIGRAD) works.

But if the minimisation and derivative estimation works (and they should in such a simple case) results will not depend of the minimisation and derivative estimation method used, it doesn't matter which software you use.
But sure, we can also compare against sherpa or lmfit and e.g. try the lmfit test cases with astropy.models.

From the plot you show I think "Minuit fit" and "Models fit with deriv" is incorrect. I can have a look tomorrow and see what goes wrong with Minuit. I also have other test-cases of nonlinear fits which I'll try with astropy.
If you have working examples besides what is in this pull request, maybe you can post them somewhere else (e.g. as a gist) so that I can check them out. But more examples / unit tests in astropy.models itself would definitely be good.

@wkerzendorf
Collaborator

@nden: I think it would be very useful to have the base of the models class accessible very soon. So would it be possible to split up this PR into different ones. I have the feeling the discussion on fitting might still take a while (and is important).

@cdeil
Collaborator

Like @wkerzendorf I also think we should merge this pull request as soon as possible so that more people will try it out and contribute / give feedback. The 0.3 release is not imminent, so we can improve the package (and adjust the API if we think it can be improved) in the coming weeks.

IMO this is ready to be merged today.

@keflavich
Collaborator
@hamogu
Collaborator
@nden
Collaborator

@cdeil I fixed the deriv method and 'domain' and 'window' attributes in the Line example above. Now it agrees with your output from Minuit.

When I wrote the polynomial classes I tried (also following what was done in numpy) to write them in a uniform way and so they all have 'domain' and 'window' attributes and the LinearLSQFitter expects these. I am starting to think this is too confusing. It may be better to keep these attributes only in the models for which it really matters (Chebyshev and Legendre sofar) and have the LinearLSQFitter check for the presence of the attributes. If no one objects I will make this change after I merge your PR.

@astrofrog
Owner

I agree about merging!

@nden - is this ready to be merged on your side?

@eteq - any objections?

@eteq
Owner

I've been waiting for this so I can do a final review before merge - I will do that right now.

@nden
Collaborator
@astrofrog
Owner

Ok, so once @eteq has reviewed this, we can merge it. @eteq - if there are any comments that are non-trivial but also non-critical, maybe you can open issues for those?

@eteq
Owner

Yep, that's what I was thinking - only small changes for this PR.

docs/models/index.rst
@@ -0,0 +1,145 @@
+**************************************
+Models And Fitting (`astropy.models`)
+**************************************
+
+.. _designed: design.rst
+
+Introduction
+============
+`~astropy.models` provides a framework for representing models and
+performing model evaluation and fitting. It supports 1D and 2D models
+and fitting with parameter constraints.
+
+It is `designed`_ to be easily extensible and flexible.
@eteq Owner
eteq added a note

The link for `designed`_ here doesn't work - I see what you meant to do, but this isn't how Sphinx documentation works (this generates a link to the rst file, and you want a link to another document that may or may not be an html file). You'll want to get rid of the .. _designed: design.rst at the top of this file once you've fixed this, as well.

Fortunately, there's a better way to do this in sphinx (see http://sphinx-doc.org/markup/inline.html#cross-referencing-arbitrary-locations if you want more detail). Use the :ref: directive. That is, change this to :ref:`designed` (or :ref:models-design if you take my comment in design.rst about changing the label name.) Note that this will replace the text with the section title, though, so what you probably want is actually something like :ref:`designed <models-design>`

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((9 lines not shown))
+`~astropy.models` provides a framework for representing models and
+performing model evaluation and fitting. It supports 1D and 2D models
+and fitting with parameter constraints.
+
+It is `designed`_ to be easily extensible and flexible.
+Models do not reference fitting algorithms explicitely
+(though exceptions are sometimes necessary) and new fitting
+algorithms may be added without changing the existing models.
+In addition models can be combined in different ways using a machinery
+that allows assigning outputs from one model into the appropriate input
+of another in a flexible way, `~astropy.models.models.LabeledInput`.
+The goal is to eventually provide a rich toolset of models and fitters
+such that most users will not need to define new model classes, nor
+special purpose fitting routines (but not making that hard to do if it is necessary).
+
+This is a work in progress but the basic framework is in place and usable
@eteq Owner
eteq added a note

I think this sentence is unnecessary - no need to short-change yourself :) Unless, of course, you want to give warning like we did in coordinates that the API might change in future versions. If so, this is probably better done as a .. warning:: directive.

@nden Collaborator
nden added a note

A lot of new features were requested/proposed (the list is in https://gist.github.com/nden/5382534).
All I meant here is that feature-wise it is incomplete but the basis for the future functionality is in place and we should make sure nothing in the current framework prevents us from adding these features. If we ever get past the documentation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/index.rst
((78 lines not shown))
+
+Add some noise
+
+>>> n = np.random.randn(90)
+>>> ny = y + n
+
+Fit a Chebyshev polynomial to the data
+
+>>> ch2 = models.Chebyshev1DModel(3)
+>>> chfit = fitting.LinearLSQFitter(ch2)
+>>> chfit(x, ny)
+>>> ch2.parameters
+[0.957, 1.931, 3.029, 4.305]
+
+.. figure:: images/cheb_fit.png
+ :scale: 25 %
@eteq Owner
eteq added a note

I like this plot, but it would be much better to implement this using the :plot: directive, so we don't have binary files in the repository (and the figure stays up-to-date with future changes to the package). If you're not familiar with :plot:, see http://matplotlib.org/sampledoc/extensions.html#inserting-matplotlib-plots . I'd also say you probably want to use :include-source: so that people reading the documentation understand how they should plot their models.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/design.rst
@@ -0,0 +1,32 @@
+.. _designed:
@eteq Owner
eteq added a note

These labels are not necessary unless you intend to use the :ref: directive to point to them. So you can feel free to remove them unless you add in some :ref:s.

Also, these labels are shared globally across all of the astropy documentation. So I recommend using a slightly more descriptive name to avoid collisions with other subpackages. Something like .. _models_design.

This comment applies to all the other .rst files with labels, as well.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((51 lines not shown))
+ [ 0.127398, 0.084932]])
+
+Evaluate the model on one data set
+
+>>> y = g1(x)
+>>> print y.shape
+(90, 2)
+
+or two data sets (any other number would be an error)
+
+>>> y = g1(np.array([x, x]).T)
+>>> print y.shape
+(90, 2)
+
+.. figure:: images/gauss1D_eval_psets.png
+ :scale: 25 %
@eteq Owner
eteq added a note

This should also be replaced with :plot: and the png file removed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/models.rst
((65 lines not shown))
+.. figure:: images/gauss1D_eval_psets.png
+ :scale: 25 %
+
+- Evaluating polynomial models with multiple parameter sets with one input data set creates multiple output data sets
+
+>>> p1 = models.Poly1DModel(1, paramdim=5)
+>>> len(p1.parameters)
+10
+>>> p1.c1 = [0, 1, 2, 3, 4]
+>>> p1.psets
+array([[ 0., 0., 0., 0., 0.],
+ [ 0., 1., 2., 3., 4.]])
+>>> y = p1(x)
+
+.. image:: images/p1d_5psets.png
+ :scale: 25 %
@eteq Owner
eteq added a note

:plot: here too.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
docs/models/new.rst
((15 lines not shown))
+if not it should subclass `models.Model`. If the model takes parameters,
+their names are stored in a list as a class attribute named `parnames`.
+Pass the list of parameter names and the number of parameter sets to the base
+class. Note, that if the method which evaluates the model cannot work
+with multiple parameter sets, `paramdim` should not be given as an argument
+in the __init__ method. The default for `paramdim` is set in the base class to 1.
+
+::