This is a fix for #2049 .
@embray please review this

Any objections to merging this PR next week?
If anyone is planning to review it but needs more time, let me know.

@embray @astrofrog @eteq @mhvk @wkerzendorf @cdeil @mdboom

In principle this seems good to me (I can't comment on the fine details of the code since I'm not really familiar with the internals at the moment). Just to make sure I understand, param_dim is the number of sets of parameters, correct? If so, I think it should be param_size or something else because param_dim confused me - it made me think of the number of dimensions of the parameters (e.g. 2-d) not the number of parameter sets.

Also, this pull request does not currently change the documentation - should this behavior not be described in the docs?

One last thing - this needs a changelog entry.

Just to add to my previous comment, I wonder whether we want to allow param_shape instead of param_dim? That is, maybe I want to compute a grid of models varying two parameters, so it would then be convenient to do:

model = MyModel(param1=<2d array>, param2=<2d array>)

and then have model(x, y) return a 2-d array if x and y are scalars, and a 3-d array if x and y are 3-d arrays. However, this could start to get pretty complex in terms of which dimension is which, so I'm not sure whether this is sensible - just prefer to raise it now rather than later :)

Now a few questions from a user's point of view:

In [7]: g = models.Gaussian1D(amplitude=[1.,2.], mean=1., stddev=1.)
...
InputParameterError: Parameter mean value inconsistent with n_models attribute 2

This error is not obvious to me. In this situation, I would expect something like 'amplitude should be a scalar value', and from your description above, it should not be assuming there are two models, correct?

In the following example:

In [4]: g = models.Gaussian1D(amplitude=[1.,2.], mean=1., stddev=1., param_dim=2)

maybe one could broadcast the scalar parameter values to param_dim?

Finally, if doing:

In [5]: g = models.Gaussian1D(amplitude=[1.,2.], mean=[1.,3], stddev=[1.,4], param_dim=2)

users might wonder why in the following case:

In [12]: g([[1,2],[3,4]])
Out[12]: 
array([[ 1.        ,  1.93846647],
       [ 0.13533528,  1.93846647]])

there are four values returned, not 8, so this needs to be made very clear in the documentation.

I also don't like the name param_dim. If it were me, I would just have a boolean flag keyword that indicated that multiple parameter sets are present, and to take the length from the length of the the first parameter. modelsets, multimodels,... I'm not sure what name I would choose but this leads to less room for inconsistency (which is still possible if all parameters aren't consistent in length).

@astrofrog
Correct that param_dim is the number of parameter sets. Or it can be looked at as the number of model instances. For example, the statement

poly = Polynomial1D(degree=1, param_dim=2, c0=[1, 2], c1=[3, 4])

creates two polynomials. It is true that the name is unintuitive but param_shape won't be quite correct either. Internally the code keeps track of the shape of each parameter. Perhaps n_models?

Also there was a suggestion to change param_dim from int to bool to mean that multiple models/ parameter sets are being created. I am fine with that and will implement it if we agree on this.

For this example:

model = MyModel(param1=<2d array>, param2=<2d array>)

I think we should avoid such nebulous model definitions. Is there a model which doesn't "know" what its parameters are? Such generic models are useful for testing but should not be allowed or it becomes very messy.

The error message in the example

 g = models.Gaussian1D(amplitude=[1.,2.], mean=1., stddev=1.)

is a generic message on the base class. In this case amplitude is a scalar but it's not in the general case. How about

Number of values for parameter "mean"  (1) inconsistent with number of models requested (2)

In the next example

 g = models.Gaussian1D(amplitude=[1.,2.], mean=1., stddev=1.)

the current code broadcasts the parameters to param_dim=2 by repeating the values of mean and stddev. The goal of this PR is to prevent this because it leads to ambiguous interpretations of what he requested number of models is and what is the actual parameter shape. Eventually it breaks the use case of having an array as a parameter.

There are some examples of using param_dim>1 in the pages below but I'll add more

http://docs.astropy.org/en/stable/modeling/models.html

http://docs.astropy.org/en/stable/modeling/parameters.html

In the example with the input

g = models.Gaussian1D(amplitude=[1.,2.], mean=[1.,3], stddev=[1.,4], param_dim=2)
In [12]: g([[1,2],[3,4]])

2 models are created because param_dim=2. The rule for evaluation is that param_dim should match the number of rows in the input (first dimension) should match param_dim or be able to be broadcasted to param_dim. Although it looks like the current behavior returns a transposed reult (I'll check that). The rules for evaluation are described in the module docstrings of modeling.core but I see now that this is not included in the documentation.
I'll fix this.

One question that remains is what should param_dim be renamed to and should we make this a boolean?

Thanks for the comments @astrofrog .

I'm still not very familiar with the models, but worry about increasing complexity. In particular, I felt @embray's arguments against optional were well-put: I do not think one should have such optional parameters (and like the solution for the example presented of having AffineTransformation and Rotation2D models instead).

AffineTransformation also seems a good example for why array-shaped parameters are useful, so I'm fine with the PR otherwise.

@nden - just regarding this message:

Number of values for parameter "mean"  (1) inconsistent with number of models requested (2)

for this example:

 g = models.Gaussian1D(amplitude=[1.,2.], mean=1., stddev=1.)

I disagree with the error message. param_dim was not set so the requested number of models is 1. The problem is with the amplitude, not the mean. Do you see what I mean?

For the name of param_dim, I agree n_models would be better.

@astrofrog Yes, you are right. (I guess I'm partially thinking in terms of the current API which broadcasts parameters.)

Let's decide first on whether to make n_models a bool instead of an int. If it's a boolean then the error message may be

if n_models is False:
    message = "Number of values for parameter amplitude (2) inconsistent with number of models (1).
if n_models is True and number of values for each parameter is not the same:
    message = "Inconsistent parameter values; the number of values for each parameter must be the same in a model with multiple parameter sets."

Does this sound better?

@mhvk I am fine with removing optional, renaming MatrixRotation2D to RotationFromAngle2D (or some better name; any suggestions?) and redefining the class to only accept an angle. Then the use case of a rotation matrix as an argument will be covered by AffineTransform.

The question is does anyone see a use case for optional. One thing that comes to mind is the long debate we had of whether to define Gaussian1D in terms of fwhm or stddev. This would provide a way to cover both. But you are right, it adds complexity.

I will remove optional for now and we can worry about this when and if we have a case which can't be solved in any other way.

@nden - those error messages look good to me! Presumably this will only appear for models where amplitude is not expected to be e.g. a 2-element array, correct?

👍 to removing optional.

n_models=<int> does not allow fitting to a grid of pixels (e.g., JWST ramp fitting), but param_dim=<array> does.

@pllim - right, that was one of the issues I also brought up above when I suggested param_shape, which could be a tuple giving the number of dimensions and shape of the input parameters. One could then have a 'model grid'.

@nden - plan for removal of optional here sounds good; it is a different case that would be best dealt with separately (e.g., I think the gaussian stddev or fwhm case might be covered also using, e.g., a parameter-based "interpretation" -- indeed, if taken as a quantity, one can do this with the scale factor inside the unit).

As for the model name, perhaps RotateByAngle2D or Rotate2DByAngle? (I prefer this over FromAngle since it makes it more explicit what the model does)

@astrofrog @pllim I didn't elaborate on n_models being bool. I'll compare it to the current API.

Currently param_dim is int and has the meaning of number of models represented by this model instance.

poly = Polynomial1D(degree=1, c0=[1, 2], c1=[3, 4], param_dim=2)
poly.param_sets
Out[4]: array([[ 1.,  2.],
       [ 3.,  4.]])

My two polynomial instances are p1=Polynomial1D(c0=1, c1=3) and p2=Polynomial1D(c0=3, c1=4). This allows fitting those two models to two data sets.

The suggestion is to keep the functionality but change the API a bit. The new parameter will be bool and if True will have the meaning of "yes, create multiple models".
The above example becomes

poly = Polynomial1D(degree=1, c0=[1, 2], c1=[3, 4], n_models=True)
poly.param_sets
Out[4]: array([[ 1.,  2.],
       [ 3.,  4.]])

Internally the model checks that the parameter values are consistent with having multiple parameter sets. The number of models is inferred from len(par) and if all parameters dont have the same length one of the errors discussed above is raised. But at the end the two models are the same and can be evaluated and fitted to multiple data sets.
The advantage of this is there's no repetition of information.
I hope this was a clearer explanation.

@mhvk RotateByAngle2D sounds good. Thanks.

@nden, your example of model dimension of (2,) sounds fine. But things can get complicated when model dimension goes higher. How will modeling handle this call?

poly = Polynomial1D(degree=1, c0=[1, 2], c1=[[3, 4], [5,6]], n_models=True)

Should I expect it to magically broadcast c0 to match a dimension of (2, 2) like c1, or will it raise an error? What about dimension of (x, y, z, ...) -- where does it stop?

@pllim If you allow this, how do you distinguish a case where c1 is meant to be an array of two elements and c0 is meant to be a scalar?

That's my point: Guessing can have undesired effects. p1 : c0=1, c1=[3,4] does not make sense to me either.

@plim There's no guessing in this PR. Here are some use cases:

The valid case with 2 models

model = MyModel(c0=[1, 2], c1=[[3, 4], [5,6]], n_models=True)

The invalid case which raises an error

model = MyModel(c0=1, c1=[3, 4], n_models=True)

A valid case with 1 model

model = MyModel(c0=1, c1=[3, 4], n_models=False)

Another valid case with 1 model

model = MyModel(c0=[1, 2], c1=[[3, 4], [5,6]], n_models=False)

Okay, @nden, as long as it is documented and I can look that up, that's fine. Personally, I would just like to tell modeling exactly the model grid dimension that I want, but it is just a personal preference.

@nden @pllim - I think the real decision to be made here is whether we only allow 1-d sets of models, or if we want to allow model grids like @pllim was suggesting. Take the following example:

model = MyModel(c0=[[1, 2],[8,9]], c1=[[3, 4], [5,6]], n_models=True)

To me, this is ambiguous, and could mean either a 2x2 grid of models with scalar parameters:

p11: c0=1, c1=3
p12: c0=2, c1=4
p21: c0=8, c1=5
p22: c0=9, c1=6

or two models with parameters with length 2:

p1: c0=[1,2], c1=[3,4]
p2: c0=[8,9], c1=[5,6]

and what I was suggesting before is that we could distinguish these with:

model = MyModel(c0=[[1, 2],[8,9]], c1=[[3, 4], [5,6]], param_shape=(2,2))
model = MyModel(c0=[[1, 2],[8,9]], c1=[[3, 4], [5,6]], param_shape=(2,))

If we don't allow 'grids' of models with more than one dimension, then @nden's suggestion is fine. However, if we want to allow 2- or higher-dimensional 'grids' of models, then a simple bool is no longer enough.

Just to add, I'm not saying I definitely want multi-dimensional grids of models to be supported, as it might make things more complicated in terms of model evaluation, but I just want to make sure we at least discuss the possibility.

I'd like to keep it simple and only support 1D sets (it's easy enough to flatten higher dimensions for the purposes of using models, and reshaping when done). Also, I'd rather not have implied broadcasting to make all the parameters consistent. I don't see much harm in making the user explicitly make each parameter have a matching number of values for multiple parameter sets.

@perrygreenfield - that's fine by me - just wanted to make sure it was a conscious decision. The only thing I would say is that this should be made clear in the documentation. If we go down that road, then I think the argument should not be n_models, but multiple_models=True/False or similar.

multiple_models is a little confusing too, since a composite model also consists of multiple models.

I still feel like the need for this argument could be eliminated entirely if a parameter has a required dimensionality. By default it would be scalar, but if we could define something like:

class AffineTransform2D(Model):
    matrix = Parameter(..., shape=(2, 2), ...)

then if an AffineTransform2D is initialized with a list of multiple 2x2 arrays then the intent should be clear (so long as any other parameters are also passed in as lists of the same length).

I remember at one point hearing an argument that any parameter should be able to accept an array of arbitrary size/shape as its value, but it's still not clear on how that should even be interpreted (since depending on the function it might not even make sense to accept an array for arbitrary parameter values).

I agree with Perry that flat sets of parameter values are sufficient.

Also,

Models don't have true parameter validation.

Well, not automatically no, because we can't determine for an arbitrary model what parameter values are "valid". But when creating a model its parameters can be given a setter function which can be used to validate that parameter whenever values are assigned to it.

@nden - on the names: it just struck me that a better name than RotateByAngle2D is simply Rotate2D; it would seem that rotating by one angle is really what most people would expect such a function to do!

My question about this: why param_shapes is needed? Maybe I'm misunderstanding this. In the end in the actual calculation numpy will deal with this? So as I see it, you can chuck whatever parameters in any shape into this and once you evaluate the models numpy will do its broadcasting.

@embray On the topic of adding a shape to Parameter. How would this work for a model which takes an array of arbitrary shape as a parameter? An example of this is a Lookup Table distortion model defined in FITS WCS?

There's one case where the current code cannot catch a user input error:

gauss = Gaussian1D(amplitude=[1, 2], mean=1, stddev=2)

I see two options

• Since every model knows what its parameters are supposed to be, make a check that the parameters and param_dim (or whatever we end up calling it) are consistent in the __init__ .
• Rely on the sanity of users and good docstrings and don't catch this but let it fail in the eval step.

I am tempted to go with the second option but it may be a bit dangerous because there may be a case where the eval doesn't fail.

@embray Any suggestions how to fix this?

On the separation of Parameter into a ParameterDescription and Parameter.

The advantage of this is that Parameter can have then a .value and .unit and perhaps .shape or any other other attribute specific to an instance.

The advantage of this is that Parameter can have then a .value and .unit and perhaps .shape or any other other attribute specific to an instance.

It already can. Bound and unbound parameters are both handled by the same class currently and there's nothing gained from changing that.

How would this work for a model which takes an array of arbitrary shape as a parameter? An example of this is a Lookup Table distortion model defined in FITS WCS?

Perhaps in addition to shape, one could specify the number of dimensions (or just the number of dimensions--obviously in one specifies the shape then that fixes the number of dimensions as well).

Or, if we did include an option to models like has_multiple_parameter_sets (I'm not saying that's what it should be called), then if that's set to True it doesn't matter what the shape/dimensions of the parameter values are. It's just assumed that the first axis is the param sets and the subsequent axes are the dimensions of the individual param values in each set. I think this is what Perry has already suggested several times.

Though I still think there's a place for required shapes for parameter values. That way it's easy to ensure that a 2x2 array is always a 2x2 array, etc.

@embray - can this be closed now that #2634 has been merged?

Resolved by alternative in #2634.