The purpose of this PR is to address the basic issue raised in #2049, as well as a broader generalization of that issue, in a way that (I hope) is coherent.

Problem background

The basic issue raised in #2049 (and earlier, in different terms, by #1680) is that earlier prototypes of the modeling package were designed to support model parameters that can have array-like values. The example often raised is in what is now the AffineTransformation2D model, which takes a 2x2 matrix as one of its parameters. But it has also been requested that in principle parameter values in general should be allowed to be arrays or scalars. This is complicated somewhat by the fact that all parameter values must be compatible with each other for broadcasting; for example it has been requested that it be possible to do something like:

model = MyModel(param1=np.array([[1,2], [3.,4]]), param2=4.5)

which is fine, since here param2 is a scalar and is likely to broadcast with param1 when they are used together in any Numpy ufunc, for example. However if all parameters can be arrays, then in principle something like this is possible:

model = MyModel(param1=np.array([[1,2], [3.,4]]), param2=np.array([1, 2, 3])

except that the way most models are implemented, that won't work. On some level, if mixing parameter shapes, it's up to the user's discretion to know whether or not what they're doing is making sense. But we still want to make sure sensible error messages are raised if things go wrong.

That's the first of two major issues addressed by this PR. However, there is another factor that complicates matters:

Model sets (formerly, parameter sets)

The above issue was partially addressed by #1680. But that fix punted on the issue of supporting multiple "parameter sets" (i.e. models with param_dim=N where N > 1). The purpose of this feature is to allow a single Model instance to actually represent a whole collection of models of the same type, but with potentially different values for each fittable parameter. The use of this is primarily for fitting--for example it would enable creating a whole "set" of Gaussians, each one of which is fit to a different plane in a data cube in parallel. In many cases this would be more efficient than, for example, looping over each plan in Python in fitting separate Gaussian2D instances to them one by one.

Previously these were referred to as "parameter sets", but I've started referring to them as "model sets" since I think that makes clearer that we're really talking about a collection of related models. I don't think either terminology is more right or wrong than the other, and am open to other suggestions.

Prior to this PR there were two ways to instantiate a model with multiple model sets. The first, when instantiating a model, was to specify a keyword argument param_dim=N where N is an integer specifying the number of models in the set. Each parameter value was then required to be a list, the length of which should match the number of models:

model = MyModel(param1=[1, 2], param2=[3, 4], param_dim=2)

Equivalently, if the value given to a parameter was a Numpy array, the length of the first axis needed to match the value of param_dim:

model = MyModel(param1=np.array([[1, 2, 3], [4, 5, 6]]), param2=[3, 4], param_dim=2)

The above example would create a set of two models, each where the value of param1 is an array of shape (3,).

Unfortunately, my fix to #1680 did not effectively address the above case, where there are multiple models, while still allowing array-like values for individual parameters in each model.

The other way of making a multiple model set, which was employed by most of the models in astropy.modeling.function_models, and which I (and others) considered a misfeature was that if the first parameter of a model (typically an "amplitude" scaling coefficient) was an array, then the length of the first axis of that array automatically determined the size of the model set. That is, in the case of Gaussian1D something like:

g = Gaussian1D(amplitude=[1, 2], mean=[0.1, 0.2], stddev=[0.1, 0.2])

automatically created a set of two Gaussians. It was impossible to allow a single Gaussian whose parameters happened to be array values.

Handling inputs to model sets

Further complicating matters is how to handle inputs to Model.__call__ for instances representing a multiple model set. For example, given

g = Gaussian1D(amplitude=[1, 2], mean=[0.1, 0.2], stddev=[0.1, 0.2])

there are several use-cases for evaluation. The first, and simplest, is evaluation on a scalar value like:

g(1.2)

This case is straightforward--it simply evaluates each model in the set on x=1.2. But we also want to pass in arrays to models, and there are two different things this could represent:

1. An array whose length is the same as the number of models N, meaning map each element in this array 1-1 with a model in the set.
2. Simply evaluate all models in the set on the given array.

That is, why might want to read:

g([1.2, 2.3])

as "simultaneously evaluate the model with amplitude=1, mean=0.1, stddev=0.1 at x=1.2 and the model with amplitude=2, mean=0.2, stddev=0.2 at x=2.3".

OR, we might want to read it as "simultaneously evaluate each of x=[1.2, 2.3] across all models in the set, elementwise".

In the current implementation, the latter interpretation is the behavior for 1-D inputs, but the former interpretation is the behavior for ndim >= 3 inputs. For 2-D inputs the number of columns must match the number of models. The current rules are a reasonable compromise for simplicity, and are internally consistent. But from a user perspective they may seem confusing ("how I do need to shape my arguments for this particular model?"), and do not simultaneously allow both of the above use cases for inputs of any size/shape (scalar, or array).

Proposed solution

For a quick overview of the new proposed implementation you can read it right in this PR in this docstring: https://github.com/astropy/astropy/pull/2634/files#diff-a619f3fef9ad65068ecd1790cebc242aR16
(note, this is probably not the final resting place of that documentation; I am still in the process of updating all the documentation, and as such this PR is not quite ready to be merged)

As this is a multi-pronged problem, the solution is also multi-pronged.

Single models

The first aspect of the solution is that for single models (i.e. not model/parameter sets) everything works as before. For the most basic of one model, with scalar parameters one instantiates it like:

>>> m = MyModel(1, 2, 3)

The model can be evaluated on either a scalar or array input. The output will have the same shape as the input:

>>> m(1)  # these are just sample values for an arbitrary model
2
>>> m([[1, 2], [3, 4]])
array([[2, 4],
       [6, 8], [10, 12]])

The same applies when one or more of the parameters have array-like values. The only general requirement is that the parameter values and the inputs are all mutually compatible according to the standard Numpy array broadcasting rules. The output has the shape that results from broadcasting together all the parameters and inputs. There is at least one caveat to this approach that I will discuss below.

Model sets

Instantiating a Model that represents a set of multiple models with their own parameter values always requires stating this explicitly. O'er in #2149, as well as in other discussions it was suggested that the existing keyword argument for specifying sets, param_dim was confusing. In retrospect it makes sense if you think of a model set as an array of parameter values with an axis specifying which model a given parameter belongs to. In that sense, param_dim is giving the size of that model axis or "dimension".

This implementation preserves that concept, but does go ahead and change the name of the corresponding argument to n_models. For example:

m = MyModel([1, 2], [3, 4], [5, 6], n_models=2)

This has the same semantics as in the previous implementation. The only difference is the name, and the fact that it is required, for any model, to instantiate a model set. It also allows the parameter values to arrays (though the values for a single parameter must have the same shape across all models in the set). For example:

m = MyModel(param1=[[[1, 2], [3, 4]], [[5, 6], [7, 8]]],
            param2=[[10, 11], [12, 13]], param3=[14, 15], n_models=2)

As before, the first dimension represents what I'm calling the "model_set_axis". This is the axis of the array that determines which model a parameter values belongs to. Its length, for all the parameters, is 2. In this case the values of param1 are 2x2 arrays, the values of param2 are 1-D arrays, and the values of param3 are scalars.

Another alternative that was discussed and that looked promising was to use a boolean argument like multiple_models. If True, it would just require that all parameter values have the same length in their first dimension, and that the number of models would be determined from that common length. Unfortunately that doesn't work in the case of models that have default values for all their parameters. The main example of this in the current collection of models is the polynomial models, which have a default of zero for all their coefficients. It should be possible to write:

p = SomePolynomialType(degree=3, n_models=2)

without having to manually specify something like [0, 0, 0] for all the default parameter values. We might be able to get away with a different approach, for example, by allowing specification of the shapes of parameter values without actually specifying the values. But for now, this approach introduces the least disruption to the existing assumptions.

Evaluation

Handling of inputs to model sets is where there is some significant departure from the current implementation. The main differences are that for n_models=N with N > 1

1. By default, all inputs must be an array the first dimension/axis of which is treated as the model_set_axis, mapping input values to models in the set. For example for n_models=2 the input must be in the form[..., ...] where the ... may be either scalars or arrays so long as they broadcast sanely. One exception is made, for convenience, in the case of a single scalar input, like m(0). In this case all models in the set are evaluated at the same input value.
2. All outputs from the model are arrays, the first dimension/axis of which still corresponds to the model_set_axis, and is of length N. Therefore the inputs and outputs always have the same length in their first dimension (the subsequent dimensions may be different depending on broadcasting rules).

What, then, of wanting to evaluate a single array value across all models in the set? In order to handle this unambiguously, when evaluating the model it is necessary to supply an additional argument model_set_axis=False. For example:

>>> m = MyModel(param1=[1, 2], param2=[3, 4], param3=[5, 6], n_models=2)
>>> x = np.arange(10)
>>> m(x)
Traceback ...
InputShapeError...
>>> m(x, model_set_axis=False)
array([[10, 20, 30, 40, 50, 60, 70, 80, 90, 100],
       [100, 200, 300, 400, 500, 600, 700, 800, 900, 100]])

Trying to evaluate the 2-model set on a 10 element array would fail, because its length does not match the number of models. But supplying model_set_axis=False means "whatever the shape of this array, none of its dimensions map to a model in the model set". Which in turn implies broadcasting the value across all models. The output still has a length equal to n_models for its first dimension. But the subsequent dimensions are the result of broadcasting the input with the parameter values of each model in the set.

The following tests illustrate several more examples of what happens with parameters and inputs of different shapes: https://github.com/astropy/astropy/pull/2634/files#diff-5a6c3ec736446d944478e4ecb9cd19b4R485 But again the key point here to emphasize is that outputs have length of n_models (in their first dimension, when n_models>1).

Arbitrary model_set_axis

The previous example showed how model_set_axis can be disabled for an input, implying that it does not have an axis mapping inputs to models in a set. Alternatively, one may specify an integer for model_set_axis. By default model_set_axis=0--that is, the first dimension is equal to the number of models. But a different axis may be specified for this purpose. A model may also have a different default model_set_axis for its parameters. I won't go into much detail on this for now, as I don't think this is likely to be used very commonly, as it requires careful consideration on the part of the user. But it may be useful for some performance critical cases in order to avoid having to swap the axes of an array around, resulting in a non-contiguous input array.

If that sounds unclear, don't worry, you're not alone. Some examples of this would help, if I have time to write them...

Caveats

1. I haven't done much performance testing, but it may be much slower, especially in the case of fitting. Fortunately, once this work is merged with the work presented in Model refactoring stage 1: "evaluate" normalization #2580, that slowdown will be all but eliminated. If the regression turns out to be serious and noticeable I can easily backport the key pieces of Model refactoring stage 1: "evaluate" normalization #2580 to 0.4.x.
2. There are still some details that are glossed over regarding models with multiple inputs and multiple outputs. In particular where the number of outputs is greater than the number of inputs--I think this currently results in a crash. I don't believe there are any cases of that in the wild yet, but it's good to be aware of. [This issue was at least partially addressed by Modeling/handle more out than in #3250]
3. The Numpy broadcasting rules don't necessarily apply to how all models are evaluated. Fortunately, they do apply to most. The only model currently built into Astropy that had a problem with this is AffineTransformation2D. That model takes two inputs: an "x" coordinate and a "y" coordinate. Or it can take an array of x coordinates and an array of y coordinates. In then uses np.vstack on these inputs to create a 2 x N array that can be dotted with the transformation matrix. For now I punted on the issue by providing a flag model implementations can set that implies "I don't follow the standard rules; let me handle input validation on my own."

I'd welcome better suggestions here. But maybe it will turn out that this is good enough if there are few enough exceptions.

To do on this PR

This PR is not quite ready to be merged. If nothing else the user documentation for models requires further updating. If there are no serious objections to this general approach, we would like to include this as a fix for #2049 in Astropy 0.4. Minor objections/quibbles can either be addressed before merging if small, or in Astropy 0.4.x. This PR doesn't purport to be an end-all be-all solution for all possible cases, but it gets a good ways closer, so long as nobody opposes the general approach.

Future enhancements

One thing I still want to do with this is that make clearer the distinction between a model and a "model set", and give a little more power to the concept. Although I like using the same class for both cases, it should be clearer when a given Model represents a set of models. I also want to implement an indexing feature, allowing something like:

>>> g = Gaussian1D(amplitude=[1, 2], mean=[0, 1], stddev=[0.1, 0.5])
>>> g[0]
<Gaussian1D(amplitude=1, mean=0, stddev=0.1)>
>>> g[1]
<Gaussian1D(amplitude=2, mean=1, stddev=0.5)>

That is, indexing a model set returns a single model instance representing one of the models in that set, which can be evaluated individually, etc. Updates to the parameters in that model would also update the original model set. This enhancement is fairly easy to implement, and I wanted to have it originally as part of this PR. But it was left out for now due to time constraints, and due to the fact that some of the work I'm doing in parallel with this will make it easier to implement.