models subpackage in astropy

[Nadia Dencheva]

Hello,

I made a PR for a models/fitting package in Astropy:

https://github.com/astropy/astropy/pull/493

The documentation is here:

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

As Perry mentioned in a different email one of the goals is to be able to add
new models and fitting routines easily. There are not many models implemented
right now but there' s at least one of each type (and a description of how to add
new ones). I also expect to be adding new models in the future.

Please take a look if you are interested, I'd appreciate any comments.
Thanks,

Nadia

[Tom Aldcroft]

Hi Nadia,

Have you considered making the model parameter be a class object with attributes like min, max, units, value, frozen.  See [1] for an example of this in the in Sherpa fitting / modeling package.

The idea that a parameter is more than just a number is very powerful in practice.  In particular this lets you do:

- Set min and max allowed range for parameter in fitting.

- Flexibly make a parameter be varied or not varied in the fit process.  This is incredibly useful for real analysis.

In terms of coding within this paradigm, it's not much more difficult, since you can define a model property parvals that returns something like [par.value for par in model.pars].

Cheers,

Tom

[1] http://cxc.harvard.edu/sherpa/ahelp/set_source.html

[Tom Aldcroft]

Hi Nadia

You might also consider making the infrastructure to separate the optimization method (e.g. Nelder-Mead, Levenberg-Marquardt, Monte-Carlo) from the fit statistic (e.g. Chi squared, Cash (log-likelihood), C-stat, model variance, data variance, user-defined, etc).  Having this separation allows working on a wider variety of problems.  As a concrete example, low count data should not be fit using gaussian statistics like chi^2.

- Tom

On Mon, Nov 19, 2012 at 6:35 PM, Nadia Dencheva <nadia.d...@gmail.com> wrote:

[Nadia Dencheva]

Hi Tom,

There are two things here: units and constraints.

It was a deliberate decision not to attach units to parameters and let the software which uses them deal with that. It was felt this would make the software more flexible.

There's support for constraints now although I have no examples of that in the documentation (will add some soon). There are two ways to set constraints:

- pass them to the fitter, for example:

model = SomeModel(par1, par2)

fitter = SomeFitter(model, fixed=['par1'] )

- or change them interactively

fitter.constraints.fixed=['par1']

This is  a little different from sherpa but the idea was to keep constraints close to the fitting routines because htey are naturally associated with fitting. Also different fitters support different type of constraints.

>>>fitting.Constraints.fitters

{'LinearLSQFitter': ['fixed'],

 'NonLinearLSQFitter': ['fixed', 'tied'],

 'SLSQPFitter': ['bounds', 'eqcons', 'ineqcons', 'fixed', 'tied']}

There' s more to be done about constraints and they definitely need more testing but this is the general idea.

Nadia

[Nadia Dencheva]

On Tue, Nov 20, 2012 at 10:27 PM, Tom Aldcroft <aldc...@head.cfa.harvard.edu> wrote:

Hi Nadia

You might also consider making the infrastructure to separate the optimization method (e.g. Nelder-Mead, Levenberg-Marquardt, Monte-Carlo) from the fit statistic (e.g. Chi squared, Cash (log-likelihood), C-stat, model variance, data variance, user-defined, etc).  Having this separation allows working on a wider variety of problems.  As a concrete example, low count data should not be fit using gaussian statistics like chi^2.

- Tom

I agree. This would be a very good addition.

Nadia

 

[Tom Aldcroft]

On Tue, Nov 20, 2012 at 11:17 PM, Nadia Dencheva <nadia.d...@gmail.com> wrote:

Hi Tom,

There are two things here: units and constraints.

It was a deliberate decision not to attach units to parameters and let the software which uses them deal with that. It was felt this would make the software more flexible.

I was mostly thinking about units as an immutable part of the model definition that is mostly ignored by the fitting code itself.  If you envision users defining and fitting astrophysical models (e.g. Raymond-Smith plasma emission) then it becomes important to know the units of each parameter and of the model output itself.  That *can* just be handled in documentation, but with the powerful units infrastructure in astropy it seems a shame to not have the units be available as code objects that can then be manipulated by users downstream.

The other complication is the distinction between integrated and non-integrated models, which is common with spectral fitting.  Often the spectral data are integrated (e.g. summed photons within each wavelength bin) while the physical model is not (photons / angstrom or whatever).  Certainly one can put the burden on the user to deal with this, but having some awareness in the fitting package is good.

 

There's support for constraints now although I have no examples of that in the documentation (will add some soon). There are two ways to set constraints:

- pass them to the fitter, for example:

model = SomeModel(par1, par2)

fitter = SomeFitter(model, fixed=['par1'] )

How will this work in the situation of a composite model where there may be namespace collisions (e.g. a model with 4 gaussians)?

[Tom Aldcroft]

Hi Nadia,

Can you support a simple model language for creating composite models?
I think this is actually relatively easy to implement and allows very
intuitive and flexible code like:

>>> x = np.arange(1,10,.1)
>>> p1 = models.Poly1DModel(1)
>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
>>> g2 = models.Gauss1DModel(10., xsigma=1, xcen=1.2)
>>> pcomptr = p1 * g1 + g2

BTW, I don't understand the current PCompositeModel. According to the
documentation:

This is equivalent to applying the two models in parallel:
>>> y = x + (g1(x) - x) + (p1(x) - x)

That looks like "y = g1(x) + p1(x) - x", which doesn't make any sense.
I also don't understand the use of the word "parallel". One can do
simple model composition either additively or multiplicatively.
(There can also be functional composition, but that is big can of
worms that probably should not be opened).

Cheers,
Tom

On Mon, Nov 19, 2012 at 6:35 PM, Nadia Dencheva
<nadia.d...@gmail.com> wrote:

[Nadia Dencheva]

On Wed, Nov 21, 2012 at 7:12 AM, Tom Aldcroft <aldc...@head.cfa.harvard.edu> wrote:

On Tue, Nov 20, 2012 at 11:17 PM, Nadia Dencheva <nadia.d...@gmail.com> wrote:

Hi Tom,

There are two things here: units and constraints.

It was a deliberate decision not to attach units to parameters and let the software which uses them deal with that. It was felt this would make the software more flexible.

I was mostly thinking about units as an immutable part of the model definition that is mostly ignored by the fitting code itself.  If you envision users defining and fitting astrophysical models (e.g. Raymond-Smith plasma emission) then it becomes important to know the units of each parameter and of the model output itself.  That *can* just be handled in documentation, but with the powerful units infrastructure in astropy it seems a shame to not have the units be available as code objects that can then be manipulated by users downstream.

The other complication is the distinction between integrated and non-integrated models, which is common with spectral fitting.  Often the spectral data are integrated (e.g. summed photons within each wavelength bin) while the physical model is not (photons / angstrom or whatever).  Certainly one can put the burden on the user to deal with this, but having some awareness in the fitting package is good.

OK, I see what you mean. I have no strong feelings about units one way or the other. Parameters are objects and it would be relatively simple to add a 'units' attribute.

Does anyone else see any pros or cons about this?

 

There's support for constraints now although I have no examples of that in the documentation (will add some soon). There are two ways to set constraints:

- pass them to the fitter, for example:

model = SomeModel(par1, par2)

fitter = SomeFitter(model, fixed=['par1'] )

How will this work in the situation of a composite model where there may be namespace collisions (e.g. a model with 4 gaussians)?

Well, right now there's a JointFitter class which allows fitting several models simultaneously while keeping a common parameter, for example fit 4 gaussians but keep the fwhm the same for all of them. This is different from your example but I'm thinking this may be extended to work with that use case or be the basis for a more general way of setting constraints. As I said in the previous email, there's more to be done about constraints.

[Nadia Dencheva]

On Wed, Nov 21, 2012 at 7:29 AM, Tom Aldcroft <aldc...@head.cfa.harvard.edu> wrote:

Hi Nadia,

Can you support a simple model language for creating composite models?
 I think this is actually relatively easy to implement and allows very
intuitive and flexible code like:

>>> x = np.arange(1,10,.1)
>>> p1 = models.Poly1DModel(1)
>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
>>> g2 = models.Gauss1DModel(10., xsigma=1, xcen=1.2)
>>> pcomptr = p1 * g1 + g2

Currently there's no support for this.

 

BTW, I don't understand the current PCompositeModel.  According to the
documentation:

This is equivalent to applying the two models in parallel:
>>> y = x + (g1(x) - x) + (p1(x) - x)

That looks like "y = g1(x) + p1(x) - x", which doesn't make any sense.
 I also don't understand the use of the word "parallel".  One can do
simple model composition either additively or multiplicatively.
(There can also be functional composition, but that is big can of
worms that probably should not be opened).

 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. 

Cheers,

Nadia

[Peter Erwin]

On Nov 21, 2012, at 3:54 PM, Nadia Dencheva wrote:

>
>
>
> On Wed, Nov 21, 2012 at 7:12 AM, Tom Aldcroft <aldc...@head.cfa.harvard.edu> wrote:
>
>
>
> On Tue, Nov 20, 2012 at 11:17 PM, Nadia Dencheva <nadia.d...@gmail.com> wrote:
> Hi Tom,
>
> There are two things here: units and constraints.
> It was a deliberate decision not to attach units to parameters and let the software which uses them deal with that. It was felt this would make the software more flexible.
>
> I was mostly thinking about units as an immutable part of the model definition that is mostly ignored by the fitting code itself. If you envision users defining and fitting astrophysical models (e.g. Raymond-Smith plasma emission) then it becomes important to know the units of each parameter and of the model output itself. That *can* just be handled in documentation, but with the powerful units infrastructure in astropy it seems a shame to not have the units be available as code objects that can then be manipulated by users downstream.
>
> The other complication is the distinction between integrated and non-integrated models, which is common with spectral fitting. Often the spectral data are integrated (e.g. summed photons within each wavelength bin) while the physical model is not (photons / angstrom or whatever). Certainly one can put the burden on the user to deal with this, but having some awareness in the fitting package is good.
>
> OK, I see what you mean. I have no strong feelings about units one way or the other. Parameters are objects and it would be relatively simple to add a 'units' attribute.
> Does anyone else see any pros or cons about this?

I can imagine some possible issues where being overly specific about units make it more difficult --
or at least more awkward -- to re-use models for situations the original writer didn't think of.

For example, if I'm fitting a 2D function to a 2D image, I might do this to an image that has units
of counts/pixel and size units of pixels, or counts/sec/pixel and arc seconds, or electrons/sec/pixel
and kpc -- or I might do it to an image from a simulation which has units of arbitrary mass
per pixel and ...

And there are all sorts of situations where one might fit, e.g., one or more Gaussians to a 1-D dataset,
regardless of what the units are.

If the units are a kind of "decoration" -- having no effect on the fitting -- this probably isn't a problem,
just a minor annoyance for the users if they're fitting something other than the "obvious" data. If the
units have *some* influence in the fitting, then it may be a problem.


-- Peter

=============================================================
Peter Erwin Max-Planck-Insitute for Extraterrestrial
er...@mpe.mpg.de Physics, Giessenbachstrasse
tel. +49 (0)176 2481 7713 85748 Garching, Germany
fax +49 (0)89 30000 3495 http://www.mpe.mpg.de/~erwin

[Tom Aldcroft]

There is usually a natural division between mathematical models (polynomials, gaussians, etc) that are basically scale-free have no units and (astro-) physical models (plasma thermal emission, accretion disk emission) that have units and typically cannot be used in a different context.  By that I mean you cannot supply degrees fahrenheit in a thermal model where temperature in keV is expected, and the model output is always in a specific set of units.

Having said that, I was thinking about units being basically a decoration in the context of the core model fitting.

[Perry Greenfield]

I

On Nov 21, 2012, at 7:29 AM, Tom Aldcroft wrote:

> Hi Nadia,
>
> Can you support a simple model language for creating composite models?
> I think this is actually relatively easy to implement and allows very
> intuitive and flexible code like:
>
>>>> x = np.arange(1,10,.1)
>>>> p1 = models.Poly1DModel(1)
>>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
>>>> g2 = models.Gauss1DModel(10., xsigma=1, xcen=1.2)
>>>> pcomptr = p1 * g1 + g2
>
> BTW, I don't understand the current PCompositeModel. According to the
> documentation:
>
> This is equivalent to applying the two models in parallel:
>>>> y = x + (g1(x) - x) + (p1(x) - x)
>
> That looks like "y = g1(x) + p1(x) - x", which doesn't make any sense.
> I also don't understand the use of the word "parallel". One can do
> simple model composition either additively or multiplicatively.
> (There can also be functional composition, but that is big can of
> worms that probably should not be opened).
>
> Cheers,
> Tom
>

I think it is functional composition that is intended. It may be a can of worms, but it's one we have (want some worms?).

Perry

[Perry Greenfield]

One thing to keep in mind is what changes are needed to keep these options open without causing interface changes in the future, vs things we we would like it to be able to do, but it doesn't yet.

In other words, we don't need it to have all capabilities in an early release, but we sure would like to have the interface not break in adding these capabilities in the future. That's what I'm most interested in.

Perry

[Tom Aldcroft]

On Wed, Nov 21, 2012 at 10:02 AM, Nadia Dencheva <nadia.d...@gmail.com> wrote:

On Wed, Nov 21, 2012 at 7:29 AM, Tom Aldcroft <aldc...@head.cfa.harvard.edu> wrote:

Hi Nadia,

Can you support a simple model language for creating composite models?
 I think this is actually relatively easy to implement and allows very
intuitive and flexible code like:

>>> x = np.arange(1,10,.1)
>>> p1 = models.Poly1DModel(1)
>>> g1 = models.Gauss1DModel(10., xsigma=2.1, xcen=4.2)
>>> g2 = models.Gauss1DModel(10., xsigma=1, xcen=1.2)
>>> pcomptr = p1 * g1 + g2

Currently there's no support for this.

This would be a limitation in terms of modeling spectra, which is frequently done with combinations of a continuum model (powerlaw or whatever) with multiplicative absorption components and additive emission components.  If that's not the goal with this subpackage, then fine, but otherwise it would be good to plan for an architecture that will be able to support this.  One example of what I mean goes back to what I said about adding min, max, val, frozen, and link attributes to the parameter object.  This provides a natural way to specify all these, e.g.:

>>> g1.fwhm.frozen = True

>>> g2.xcen = g1.xcen  # link two parameters

I like this object-model way of controlling parameters much more than the alternative of:

>>> fitter = SomeFitter(model, fixed=['par1'])

In my work I do a lot of interactive model fitting and there is a lot of freeze/thawing of parameters and playing with constraints, etc.  A particular model is a fully contained entity that knows and remembers the state of all the parameters.  Without this, each time you want to make an incremental change to the parameter state, you need to essentially respecify the whole list.  E.g.

>>> fitter = SomeFitter(model, fixed=['nH', 'redshift', 'xcen'], minmax={'fwhm': [0.2, 0.5], 'gamma': [1.5, 2.5]})

>>> results = fitter.fit()

>>> # OK, allow more variation on gaussian fwhm and let 'xcen' be free to vary.

>>> # Now I have to repeat many of the prior constraints and define a new fitter.

>>> # The model is now partly specified by "model" and partly by the keyword args provided to SomeFitter().

>>> fitter = SomeFitter(model, fixed=['nH', 'redshift'], minmax={'fwhm': [0.2, 1.0], 'gamma': [1.5, 2.5]})

Compare:

>>> model = p1 * abs1 + g2

>>> fitter = SomeFitter(model)

>>> abs1.nH.fixed = True

>>> abs1.redshift.fixed = True

>>> g1.xcen.fixed = True

>>> g1.fwhm.min = 0.2

>>> g1.fwhm.max = 0.5

>>>  # etc

>>> results = fitter.fit()

>>> # OK, allow more variation on gaussian fwhm and let 'xcen' be free to vary.

>>> g1.xcen.fixed = False

>>> g1.fwhm.max = 1.0

>>> results = fitter.fit()

Without this attribute access then it's not as straightforward to specify a particular parameter of a particular model component (assuming that additive and multipilictive model composition is allowed). Of course different fitters don't necessarily have to support all the constraints provided by the parameter object.  When calling the fitter fit() function it could check the parameters and issue warnings or exceptions if it cannot support the constraints in the model parameters.

Now that I've written this longish screed, I think it boils down to whether the Model owns the parameter constraints or whether the Fitter owns the constraints.  I think that it's most natural for the model to own the parameter constraints since it already owns the parameters, and that the Fitter simply implements the model constraints in the fit process.

Cheers,

Tom

[Tom Aldcroft]

So (to me) that means you have models f(x) and g(x), and the Model class with then let you create a new model f(g).  If there will be support for this, then support for f(x) + g(x) and f(x) * g(x) and f(x) + g(x) * h(x) should be a piece of cake (or pumpkin pie).

 

Perry

[Perry Greenfield]

On Nov 21, 2012, at 12:53 PM, Tom Aldcroft wrote:

> I think it is functional composition that is intended. It may be a can of worms, but it's one we have (want some worms?).
>
> So (to me) that means you have models f(x) and g(x), and the Model class with then let you create a new model f(g). If there will be support for this, then support for f(x) + g(x) and f(x) * g(x) and f(x) + g(x) * h(x) should be a piece of cake (or pumpkin pie).
>

Right to the first part anyway. I'm not sure I'm that good of a cook though :-)

I could believe that fitting such composite model simultaneously would be very problematic. In these cases it is very likely that we will do the fits on the components, and then create the composite.

Perry

[Perry Greenfield]

On Nov 21, 2012, at 12:47 PM, Tom Aldcroft wrote:

>
> Without this attribute access then it's not as straightforward to specify a particular parameter of a particular model component (assuming that additive and multipilictive model composition is allowed). Of course different fitters don't necessarily have to support all the constraints provided by the parameter object. When calling the fitter fit() function it could check the parameters and issue warnings or exceptions if it cannot support the constraints in the model parameters.
>
> Now that I've written this longish screed, I think it boils down to whether the Model owns the parameter constraints or whether the Fitter owns the constraints. I think that it's most natural for the model to own the parameter constraints since it already owns the parameters, and that the Fitter simply implements the model constraints in the fit process.
>
> Cheers,
> Tom
>

That's a good question. In principle I think you are right, but my worry is that given all the varieties of possible constraints that will make the parameter class a monster (or the model classes). As long as we don't preclude eventually moving that capability into that class, I'm guessing the current approach is easier in the short run. I'm happy to have people think about and propose some parameter or model-side solutions to constraints.

Perry

[Tom Aldcroft]

The parameter class would be remain simple, and all of the hard work in implementing constraints would still be in the Fitter class just as it is now.  At a simplistic level, instead of having "fixed" and other constraint specifications as kwargs to the Fitter object initialization, take them dynamically from the model when doing the fit:

class SomeFitter(BaseFitter):

    def __init__(self, model):

        self.model = model

    def fit(self):

        pars = self.model.pars

        fixed = [par.fixed for par in pars]

        mins = [par.min for par in pars]

        maxes = [par.max for par in pars]

        links = [par.link for par in pars]  # par.link is a par object

        vals = [par.val for par in pars]

        # NOW really do the fit using vals, fixed, mins, maxes, links, etc.

        # At this time check that the requested constraints are implemented

        # by this Fitter class.

So from the code perspective it isn't a huge change, but from the user API perspective the change is significant as shown in the previous examples.  That's why I'm pressing a bit on this now rather than later.

BTW, if the fitter is prepared to handle certain complicated constraints, it is possible to specify them with this parameter object formalism.  Imagine you have a model with parameters x, y, and r2:

>>> model.r2 = model.x * model.x + model.y * model.y  # creates a link to an expression

>>> model.r2.max = 4.0

Just food for thought.  :-)  Happy Thanksgiving, the Government (and SAO) is calling it quits at 2pm today.

- Tom

[Peter Erwin]

On Nov 21, 2012, at 6:47 PM, Tom Aldcroft wrote:

> This would be a limitation in terms of modeling spectra, which is frequently done with combinations of a continuum model (powerlaw or whatever) with multiplicative absorption components and additive emission components. If that's not the goal with this subpackage, then fine, but otherwise it would be good to plan for an architecture that will be able to support this. One example of what I mean goes back to what I said about adding min, max, val, frozen, and link attributes to the parameter object. This provides a natural way to specify all these, e.g.:
>
> >>> g1.fwhm.frozen = True
> >>> g2.xcen = g1.xcen # link two parameters
>
> I like this object-model way of controlling parameters much more than the alternative of:
>
> >>> fitter = SomeFitter(model, fixed=['par1'])
>
> In my work I do a lot of interactive model fitting and there is a lot of freeze/thawing of parameters and playing with constraints, etc. A particular model is a fully contained entity that knows and remembers the state of all the parameters. Without this, each time you want to make an incremental change to the parameter state, you need to essentially respecify the whole list. E.g.
>
> >>> fitter = SomeFitter(model, fixed=['nH', 'redshift', 'xcen'], minmax={'fwhm': [0.2, 0.5], 'gamma': [1.5, 2.5]})
> >>> results = fitter.fit()

To be honest, I found the this first approach easier to follow than your second (preferred) approach,
and it seemed to involve a lot less typing ;-) (given that one can recall the initial specification line with
a few up-arrow keypresses...)

> Now that I've written this longish screed, I think it boils down to whether the Model owns the parameter constraints or whether the Fitter owns the constraints. I think that it's most natural for the model to own the parameter constraints since it already owns the parameters, and that the Fitter simply implements the model constraints in the fit process.
>

I can imagine more complicated constraints which might more easily be implemented in
the "fitter" -- e.g., wanting to ensure that parameter g1.xcen is always < g2.xcen -- which would
be difficult to specify at the level of individual component models (how is model g1 supposed
to know about the possible parameters in other models?).

But I certainly don't have any strong opinions about this; it could be that your model-based
constraints end up a better solution than "fitter-based" constraints.

cheers,

[Wolfgang Kerzendorf]

Hey guys,

First of all, thanks for the model implementation Nadia. Now straight to business: 

I agree with some of what Tom says: I have used (and created) quite a few model/fitter classes and bounds/fixed(frozen) have been immensely helpful (an implementation that I really like it pyminuit). I think they should be part of the model. This is not only useful for fitting but for input as well (e.g. I have a web interface to one of the models - I want to test if the user's inputs are out of range). 

I think that the composite models, are one step further and we should focus first on the very simple things described above (same goes for units). 

@peter.erwin: I think that even if the model has some boundaries - there's no reason (for special cases) not to implement some bound in the fitter. 

my 2 turkeys,

    Wolfgang

[Nadia Dencheva]

Hi Tom,

May be I missed this but how do you see the initial specification of constraints at the parameter level? (You agree that the heavy lifting is still going to be within the Fitter class.) In the current scheme, if I want to link parameters I would say something like:

def linkfunc(model):

    model.x1 = model.x2+model.x3

    

fitter=Fitter(model, tied={'x1': linkfunc})

or I can say:

fitter = =Fitter(model)

fitter.constraints.tied = linkfunc

So I can dynamically change the constraints but if I want to specify them initially this goes in the Fitter's initialization.

If I move them to the Parameter class, then I imagine I woul dhave to do something like:

model=SomeModel(par1, par2)

model.tied=linkfunc

fitter=Fitter(model)

So in a way this seems (at least on the surface) a bit more limiting in the sense that there;s no way to initialize the model or the fitter with the constraints.

Now, the current scheme has a limitation of not being able to specify constraints from composite models but I think this will require a change to the code and not the interface.

I'm happy to do it either way, just wondering which is more intuitive.

Another question is what to do with more general constaints. If  parameter constraints belong to the Parameter's class, where do we keep constraints on the entire solution? These are constraints in which the entire solution obeys a certain equation or an inequality relationship?

Currently everything is in one place : fitter.constraints.

One option is to make the Constraints class shared by Fitters and Models (like the Parameters class).

Nadia

[Tom Aldcroft]

On Wed, Nov 21, 2012 at 2:23 PM, Nadia Dencheva <nadia.d...@gmail.com> wrote:

Hi Tom,

May be I missed this but how do you see the initial specification of constraints at the parameter level? (You agree that the heavy lifting is still going to be within the Fitter class.) In the current scheme, if I want to link parameters I would say something like:

def linkfunc(model):

    model.x1 = model.x2+model.x3

    

fitter=Fitter(model, tied={'x1': linkfunc})

or I can say:

fitter = =Fitter(model)

fitter.constraints.tied = linkfunc

So I can dynamically change the constraints but if I want to specify them initially this goes in the Fitter's initialization.

If I move them to the Parameter class, then I imagine I woul dhave to do something like:

model=SomeModel(par1, par2)

model.tied=linkfunc

I would rather see something more explicit like:

def linkfunc(model):

    return model.x2 + model.x3

model.x1.tied = linkfunc  

# OR even simpler, but maybe not to everyone's taste

model.x1 = linkfunc

In this way when you look at this line of code you know immediately that the x1 parameter is being tied to something.  As written, your original example could be modifying any of the model parameters.  Maybe that example was not the final intended interface.

 

fitter=Fitter(model)

So in a way this seems (at least on the surface) a bit more limiting in the sense that there;s no way to initialize the model or the fitter with the constraints.

I think any of the parameter constraints we've been discussing (fixed, min, max, tied) could be put into the Model initializer in just the same way you envision for Fitter.  It's basically the same information, just located somewhat differently.  Maybe I'm not seeing some problem there that you have in mind.

 

Now, the current scheme has a limitation of not being able to specify constraints from composite models but I think this will require a change to the code and not the interface.

I'm happy to do it either way, just wondering which is more intuitive.

We've already seen some opinions on both sides.  :-)

 

Another question is what to do with more general constaints. If  parameter constraints belong to the Parameter's class, where do we keep constraints on the entire solution? These are constraints in which the entire solution obeys a certain equation or an inequality relationship?

Currently everything is in one place : fitter.constraints

Interesting question, but I guess I'm not completely clear on what you mean by "the entire solution".  Could you give an example?

- Tom

[Nadia Dencheva]

Hello,

Following the discussion in this thread, I have moved constraints from fitters to models and updated the PR :

https://github.com/astropy/astropy/pull/493

There were other suggestions for improvements but they don't require changes in the code structure and I plan to do them later (roughly in the order below):

- Add .units attribute to parameters and models.

This is useful for astrophysical models. Purely mathematical models and their parameters will have .units=None.

One thing to consider is whether this change will have any adverse effect on the future WCS package. The WCS class planned at the last meeting has a 'transform' (an instance of models.Model) and a 'coordinate_system' (an instance of a class from the coordinates package). If this plan holds, then both WCS.coordinate_system and WCS.transform will have a 'units' attribute. The WCS object should have its own units attribute which presumably points to WCS.coordinate_system.units.

I think this will work, hopefully won't be too confusing.

If there's no objection I'll add the units attribute to models and parameters.

- Add the ability to fit composite models (when possible)

- Support for model arithmetic expressions, like model1 + model2

- Separate optimization methods form fit statistic

- Add parameter confidence intervals

Thanks for all comments,

Nadia