A comment unrelated to the structural changes in this PR: can we add a description string to Parameter?

E.g. in #1048 it's not clear to me what scale is ... at the moment I have to look at the eval code to see that it is what I usually call amplitude or norm (see e.g. powlaw1d parameter names in Sherpa).
There's been so much discussion about names ... I don't want to do this for PowerLaw here ... it's just an example that shows that setting a Parameter description (default empty string) can be useful.

I have to admit that I don't fully understand the advantages / disadvantages of the change proposed here.

Just one observation: intuitively I would have thought a Parameter stores it's value and a Model has (i.e. stores references to) parameters. Why is it the other way around here? The model stores the parameter value and the parameter has a reference to the model.

Here's a high-level description of what RooFit, the most advanced (and most popular / used I think) modelling / fitting framework used by physicists does: PDF
Parameters store their value and models have pointers to parameters (I didn't check the code if Parameters actually have links back to models, but I don't think so).

Maybe before merging this we can increase the number of unit tests to make sure this works well for all common use cases, e.g. for linked parameters and constraints and composed models?

@piti118 has a lot of experience from writing iminuit and probfit.
Piti, if you have time to look at the emerging astropy.modeling and offer your perspective on the design, that would be great.

Regarding the addition of a description attribute to Parameters, that's a suggestion I had buried somewhere in the wall of text that explains this PR. This will certainly make it easier to implement such a feature.

In this design the model still stores references to parameters. But the value referenced by the model is just a floating point value or an array or whatever the parameter's value is. The Parameter object just wraps the reference to this value and gets any other metadata about itself from the model instance that it's bound to. This makes it easier, meanwhile, for the model to access the lighter-weight unadorned parameter value for use in calculations internally. In other words, the Parameter wrapper is just a convenience for users and doesn't really add anything internally.

This could be changed of course--the model could store a Parameter object internally and have that carry around all the constraints for that parameter as well. But I guess I didn't like that pieces of information were being copied around between different objects.

Either way, another point in this PR is that there needs to be a Parameter descriptor that is distinct from the Parameter instance int that it works like a property--it defines what parameters are available on a model, and enables things like assigning a parameter value as a float, for example, without having to manually create a Parameter instance.

I still need to look more at this since I don't fully understand the Parameter implementation yet. I like a lot of things about the new Parameters. Is it necessary though to have a metaclass with the sole purpose of generating the list of parameter names? This looks like an overkill. The current approach requires a list of parameter names and constructs parameters from them. The proposal here does the opposite. Another thing I am concerned with is how this will work when fitting composite models. I am not saying it won't be possible with this proposal, I don't understand everything yet.
@cdeil I think we should merge all tests that you wrote to give a better idea of what is expected to work. I will try to check in my changes to constraints and polynomials as well.
@iguananaut I am curious why you think adding properties to all models with a descriptor is a bad idea. Models may have many parameters. To have to write the same kind of property for every parameter would be tedious.

I think if you want a pretty good demonstration of the advantages of this approach take a look at the roations module:

https://github.com/iguananaut/astropy/blob/97addf39bcf69b4740967131f762144ac6c54155/astropy/modeling/rotations.py#L29

vs.

https://github.com/iguananaut/astropy/blob/master/astropy/modeling/rotations.py#L29

It requires a lot less work to implement a model, and with fewer opportunities for mistakes.

This also addresses #980 on some level, albeit indirectly. It would be considerably easier now to write a sort of model_factory() function that generates a Model subclass from a function. It could generate the parameters either be introspection on that function's arguments, or through additional arguments to the model_factory() specifying details about the parameters. I don't think it's a particularly necessary thing to have so long as it's easy to create subclasses, but it might not hurt either.

We avoid this kind of declaring parameter in iminuit and probfit. Though this makes things very clear think it's just too repetitive.

iminuit and probfit take advantage of introspection for funccode.co_varname (see iminuit.describe) so that user do not need to declare parameter like that.

@nden The approach used here of using Parameter as a descriptor is more powerful than simply providing a list of parameter names. For example, previously if you wanted to provide a validation function for a parameter you had to supply it to a separate class-level check_params dict. Now it's just passed in as the setter keyword argument to the Parameter descriptor. Another use case is the one @cdeil brought up of adding a description argument to each parameter. If parameters are just described by a list of parameter names, then to add this you might also need a class-level dict of parameter descriptions or something like that. In this view using a metaclass is not at all excessive--it's just the right way to do this sort of thing.

Compare, for example, this (with what the description feature might look like added in to each example):

class MatrixRotation2D(Model):
    def __init__(self, rotmat=None, angle=None):
        if rotmat is None and angle is None:
            raise InputParameterError("Expected at least one argument - "
                                           "a rotation matrix or an angle")
        if rotmat is not None:
            self._validate_rotmat(rotmat)
            self._rotmat = Parameter('rotmat', np.asarray(rotmat) + 0.,
                                      self, 1, description='arbitrary rotation matrix')
            super(MatrixRotation2D, self).__init__(param_names=['rotmat'], n_inputs=1,
                                                                                n_outputs=1, param_dim=1)
        else:
            self._validate_angle(angle)
            self._angle = Parameter('angle', np.deg2rad(angle), self, 1, description='rotation angle')
            super(MatrixRotation2D, self).__init__(param_names=[], n_inputs=1,
                                                                            n_outputs=1,param_dim=1)
            self.param_names = ['angle']
            self._rotmat = Parameter('rotmat', self._compute_matrix(angle),
                                      self, 1, description='arbitrary rotation matrix')
        self._n_inputs = self._rotmat[0].shape[0]
        self._n_outputs = self.n_inputs
        self._parcheck = {'rotmat': self._validate_rotmat,
                                     'angle': self._validate_angle}

    @property
    def angle(self):
        return Parameter('angle', np.rad2deg(self._angle), self, param_dim=1, description='rotation angle')

    @angle.setter
    def angle(self, val):
        self._angle = Parameter('angle', np.deg2rad(val), self, param_dim=1,  description='rotation angle')

to this:

class MatrixRotation2D(Model):
    angle = Parameter('angle', getter=np.deg2rad, setter=_validate_angle,
                        description='rotation angle')
    matrix = Parameter('matrix', setter=_validate_matrix,
                         description='arbitrary rotation matrix')

    def __init__(self, matrix=None, angle=None):
        if matrix is None and angle is None:
            raise InputParameterError("Expected at least one argument - "
                                           "a rotation matrix or an angle")
        if matrix is not None:
            super(MatrixRotation2D, self).__init__(n_inputs=1, n_outputs=1,
                                                   param_dim=1)
            self.matrix = np.asarray(matrix) + 0.0
        else:
            self.param_names = ['angle']
            n_inputs = self._matrix[0].shape[0]
            super(MatrixRotation2D, self).__init__(n_inputs=n_inputs,
                                                   n_outputs=n_inputs,
                                                   param_dim=1)
            self.angle = angle
            self.matrix = self._compute_matrix(angle)

@piti118 It's not terribly repetitive if subclassing existing models. It also provides a capability for a great deal more expressivity such as attaching descriptions, custom getter/setters, and default constraints to each parameter. That said, I've already written here that using introspection is also perfectly valid as a means of quickly generating new models.

In fact, fitters could do just like iminuit appears to do and, when passed in a function, create a Model instance out of it via introspection. But for common functions like those that astropy.modeling already provides a little additional descriptive power is nice too.

Now I have time to write a long reply. Here is what I was considering whether to go RooFit style or PyMinuit style.

1. Performance: RooFit style of setting parameters --> call compute() performance would be extremely bad because of one need to do a bunch of name look ups when trying to compute a function. Method call in python is extremely slow. This is crucial for fitting since it will be called again and again. But, I'd say that this kind of performance hit could be avoid even with RooFit style if one think about this a little bit. This name look up minimization should be thought of thoroughly when designing the interface.

Subclassing would create another performance problem too. I think it does a bit more expensive name look up.

2. Learning curve and Repetitiveness:

IMO PyMinuit style has much less learning code. User can get started right away. There is nothing user need to know about the internal of PyMinuit to get started. def xxx: return yyy works magically out of the box. However this can be alleviate if tutorial is written nicely though.

Repetitiveness should be avoid. I absolutely abhor RooRealVar etc. I hate having to tell it again what the parameter names are (and which order).

3. Manipulating Combining and Sharing parameters:

This could be done with both styles though. RooFit is obvious on how to implements this. Probfit relies on FakeFunc object to replicate func_code.

However, this should be obvious from user point of view which one they share too. I think this could be made obvious in both. probfit use same-name policy. But if one rely on shared parameter to be the same object, it should be obvious to user that sometime event hough parameters have the same name they aren't really shared in the fit.

4. Built-in constraint: This is obvious for RooFit style but this could be done in iminuit style too. I haven't implemented this yet but the idea is that one can basically add any field to an object, including function objects, on the fly:

x = lambda p: (p-2)**2
x.constraints = {'p':(0,1)}

5. Yet... if we have competing packages we might want to use different approach so we can learn which one actually work best.:smile:

One suggestion though Parameter constructor could take another parameter as an object so that it automatiaclly accumulate the parameters

class Gauss:
    def __init__(self):
        mu = Parameter('mu',self)
        ...

#in Parameter
class Parameter:
    def __init__(self, name, acc=None):
        acc.param_name.append(name)
        ...etc...        

Under this proposal users can get/set parameter values using normal attribute access. For example:

>>> g = Gaussian1DModel(amplitude=1.0, mean=0.5, stddev=0.3)
>>> g.mean
0.5
>>> g.amplitude = 4.0

etc. This is slower because it's going through the descriptor lookup process, but the advantage is that it can also transform the value, such as how the Rotation models covert between degrees and radians. But then in the __call__() implementation that does the actual computation it's not necessary to look up those parameters over and over again through the descriptor since they're all stored directly in the model as semi-private attributes:

def __call__(self, x):
    amplitude = self._amplitude
    mean = self._mean
    stddev = self._stddev
    return ...

@nden's implementation already does something like this by shoving each parameter into an array and then using it in Numpy expressions, so it doesn't really matter much either way how the parameters are stored in the object-oriented interface.

That's not what I meant though. The problem is that chi^2 minimization will need to call this __call__ many times and it will do name this self._name look up everytime.

def __call__(self, x):
    amplitude = self._amplitude
    mean = self._mean
    stddev = self._stddev
    return ...

I see what you're saying. This can still be gotten around though by having a means on the model to return a list of parameters (the fitter only needs to look this up once):

def __call__(self, x, **parameters):

and then the fitter can call this repeatedly, passing in the parameter set each time and avoiding the repeated attribute lookups.

This has little to do with my implementation v.s. the existing one in the modeling package, though the point you bring up here is a a good one to keep in mind whichever way we go. The main point of my changes here is to just make it easier to define the parameters on a function when using a class is desired. Taking that heavier-weight class and turning into a function that's fast to evaluate is trivial--the function itself could even be written in Cython or something.

Since a lot of the discussion now is about speed ... it's probably easiest to just give it a try (see #1095) to make sure the design does eventually support fast modeling / fitting.

In the current implementation parameters is a list and this is what the fitter uses in the call method. If the parameters list is used to evaluate the model it will be fast and this is how it initially worked. However, the current implementation uses param_sets to do the evaluation and (may be ironically) this was done in part to support faster evaluation of multiple models and partly to support constraints in fitters which don't naturally support them. An option in the current implementation would be to pick the cases when only one model is used without constraints and use the parameters list to do the fitting without going through the descriptor. The main block to using only the parameters list is having fitters which don't support constraints. With the future work to separate optimization algorithms from fit statistic this may be remedied.

Just FYI. Dictionary argument unpack is quite expensive though; you probably want to do tuple unpack instead.

def __call__(self, x, *arg)

iminuit bottleneck right now is actually the tuple unpacking.

Doesn't necessarily even need to do that if the number of parameters is fixed. For models with a variable number of parameters hypothetically one can auto-generate a function with a fixed number of arguments once all the parameters are fixed (though it would have to be updated every time the parameters are changed). But that sort of strikes me as overoptimization, at least for this package, at least for now.

Just giving my 2 cents though. I'd recommend you thought about the whole process how the fitter will call the function; how will it manipulate the parameters. I do not know what kind of fit astro physicist typically do. If it's just binned chi^2, it probably doesn't really matter but when you have to do things like unbinned LH, these name lookup will make the fitting package unusable. The proposal is a subset of what RooFit does in RooRealVar etc. RooFit has limit and unit built in to parameter too IIRC.

Here is a good concern: how do you combine models.
If I need to fit the sum of two gaussian with shared width but not the mean to a dataset.

How would this model do it: sharing Parameter object or sharing the name?

@piti118 I appreciate your feedback on this, but most of the points you're bringing up are outside the scope of this PR, the purpose of which is to offer some, what I think are improvements, to the existing API we have implemented in Astropy. This is a high-level object oriented API meant to be user-facing.

The details of, for example, how to efficiently fit parameters to a function is something that still needs to be explored but those are lower-level details than this PR is meant to address. Combining models is mostly outside this PR's scope as well. Whether this PR is accepted either in part or in whole, the Astropy modeling package has to address the issues you're bringing up. I would propose moving this discussion to astropy-dev mailinglist where I'm sure more of your input would be greatly appreciated. But for now I'd like to limit discussion on this PR to comparing it with the current implementation in master and the relative weaknesses or strengths of both versions of the API.

I have looked now at Erik's proposal in some more detail and ran some tests.

First, this is not an API only change. The fundamental idea on which this package was based has changed. And this is that Parameters is a list of float values which is shared between models and fitters. This ensures that the iteration of the minimizer is as fast as possible. Individual parameters can be updated after the minimizer is done. So in this respect @piti118 's comments are quite relevant to this discussion. In my opinion the Parameters class implementation should stay as it is now since the implementation proposed here would degrade performance and does not improve the API.

The part which changes the API is the Parameter class. I think this implementation is an improvement over the previous one. What I like most about it is the way it returns the value of the parameter. Also it may be a bit easier to define parameters for new classes this way but that is to be seen. I would like to see a full (final?) implementation of the Parameter class, one which supports constraints as well.

I played a little bit with a "mixed" implementation - one which uses the Parameter class from this proposal and the Parameters class from the earlier proposed implementation. This may be a good way to go if the new Parameter class supports constraints and if most think it provides an improved API.

I didn't do any performance testing on the Parameters list, but that makes sense that it's too slow. I still find it problematic, however, that parameter values are copied between the Model object and the Parameters object. The problem being solved there is that it's too fragile to have to keep these two objects in sync with each other, when really the purpose of the Parameters object is just to provide a view of the parameter values that's easier for the fitting routines to use.

One other possible implementation I mentioned, I think, in commit in which I changed the Parameters class, is that it all the parameter values should just be stored in an array in the first place on the Model object. That is, by default Model._parameters would just be an Numpy array of the parameters in the appropriate order. Returning the parameter values via attribute access would just get them from slices out of that array.

I'll see if I can come up with an implementation to demonstrate this idea.

Either way, as @nden wrote the changes to Parameters and to Parameter (singular) that I made here are basically separable and perhaps should have been made as separate PRs.

@iguananaut: So I put up a new API PR up at astropy/astropy-api#9 and feel that I'm proposing exactly the same thing as you. Do you mind having a quick look?
I also saw that @cdeil asked for a description parameter and I support this +1. I worked with the models for WCS purposes and they are functional but a little clunky.

I rebased my original PR on top of the (more or less) latest master. I didn't change or updated anything so several tests should be broken. Now I can begin fixing the things that need fixing to work with the several new tests (clean up handling of constraints, mainly) and then I need to address the performance concerns, especially, though not limited to the Model.parameters attribute.

Ack, ignore the last push of this branch. I accidentally overwrote the version of the branch that I had rebased on master due to git not being configured as I intended on this system. Hoping I can fix this but I don't know...

@iguananaut May the git-god be with you and help you in these times of trouble.

I've made a PR against this one that fixes the sphinx warnings ... it was just one docstring.

One of the examples on this docs page doesn't work any more:

In [39]: %paste
>>> y = p1(x)
>>> p1.c1.fixed = True
>>> pfit = fitting.LinearLSQFitter(p1)
>>> pfit(x, y)
>>> p1.param_sets
## -- End pasted text --
WARNING: The fit may be poorly conditioned
 [astropy.modeling.fitting]
ERROR: ValueError: could not broadcast input array from shape (4) into shape (5) [astropy.modeling.fitting]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-39-d33e58f16a0c> in <module>()
      2 p1.c1.fixed = True
      3 pfit = fitting.LinearLSQFitter(p1)
----> 4 pfit(x, y)
      5 p1.param_sets

/Users/deil/code/astropy2/astropy/modeling/fitting.py in __call__(self, x, y, z, weights, rcond)
    484             warnings.warn("The fit may be poorly conditioned\n",
    485                           AstropyUserWarning)
--> 486         self.fitparams = lacoef.flatten()
    487 
    488 

/Users/deil/code/astropy2/astropy/modeling/fitting.py in fitparams(self, fps)
    148         """
    149 
--> 150         self._fitparams[:] = fps
    151 
    152         if any(self.fixed) or any(self.tied):

ValueError: could not broadcast input array from shape (4) into shape (5)

Do we need to modify the example or is this a real issue?

@cdeil The failure of that example is documented in #1624. I already have a fix for it locally. I'm also going through the docs and updating them to better reflect what's in the PR. I'll merge in your doc fixes as well since it looks like you got some things I missed. Then once Travis is done testing I'll merge this.

sounds good, @embray ! Maybe you can open an issue about the parameters performance concerns? If consensus is acheived by next week, a separate PR can remove it after this is merged, or if not it will just get pushed to 0.4.

I made a few more fixes to the documentation. Gonna wait for the Travis build to finish and then I'll merge.

The auto-generated API documentation still needs work (true for much of Astropy, unfortunately). A particular problem I noticed, however, that in the documentation for each model, the docstring used for each Parameter is just a repeat of the entire Parameter class's docstring. Obviously individual Parameter descriptors should have their own documentation. @wkerzendorf had previously suggested adding a description argument for Parameter. So after this is merged I will go ahead and do that--the description can override individual parameters' __doc__ so as to fix this problem. I don't have time to go filling in descriptions of every parameter of every model (but if anyone would like to do that, that would be awesome); in the meantime the description will default to blank for each parameter.

Another thing that needs to be fixed is the boilerplate text describing constraints that's appended to each model's docstring. It should be added to the list of class arguments and not just appended to the end of each class's docstring. I'll see if I can figure something out for that, but if not I might just disable it for now, as it can be confusing in some cases.

Sounds good!

@embray - Travis is 8 builds behind, so it may take a few hours to get done. Should whoever notices first that the build passes merge, or do you need to do any final checks?

It's fine as long as the build passes. I'm still concerned about whether or not the docs build will pass.

Looks like it passed, so I'm pulling the trigger. Thanks for all the contributions!

🎆