[AstroPy] Composite models or input arrays

[Erik Bray]

On 08/25/2015 07:43 AM, Eric Emsellem wrote:
> Hi
>
> I just realised something odd (but nice!) about astropy.modeling.
>
> If I have a model, e.g., with 3 parameters (p1, p2, p3), I can do two
> things to get some sort of composite model :
>
> Option 1
>   >>> M1 = mymodel(p1=1,p2=1,p3=0)
>   >>> M2 = mymodel(p1=2,p2=2,p3=0)
>
>   >>> Mtotal1 = M1 + M2
>
> or
>
> Option 2
>   >>> Mtotal2 = mymodel(p1=[1,2],p2=[1,2], p3=[0,0])

Or Option 3:

MySummedModel = MyModel + MyModel

model_instance = MySummedModel(1, 1, 0, 2, 2, 0)

> In the first case, I'll have parameters set up one block after the other
> (parameters = [1, 1, 0, 2, 2, 0]) while in the second example, I'll have
> them as I did write them: parameters = [1, 2, 1, 2, 0, 0].
> In the first case, if I evaluate the model it gives the sum as written
> in Mtotal1(x) = (M1+M2)(x), while in the second case, I get [M1(x),
> M2(x)], so the broadcasting is different indeed.

That's not really the same as what I would call "broadcasting".  In the latter 
case you're treating the parameters as vectors.  This "works", but how this is 
defined, if at all, is rather dependent on the model.

> Both ways are useful. On one hand, Option 1 is the natural way to work
> out a composite model. But Option 2 is using an internal broadcasting I
> guess, hence it is just a matter of providing the right input arrays of
> parameters, and you can always do the sum after the fact by doing e.g.,
> Mtotal(x).sum() for example.

Indeed you could do something like that.  But if you want to treat a single 
model instance as a *collection* of models (of which you would then sum the 
results) there's a subtly different thing you should do.  Read on....

> However, Mtotal2 does not seem to act on arrays (again broadcasting
> issues?), so I cannot do Mtotal2(x) where x is an array, while it seems
> I can do that with Mtotal1 and it will return the right output.
>
> Is that correct, or am I on the wrong track here?

What you want is to use a model set.  This is a still poorly understood feature 
and I'm open to suggestions for how to make this more "obvious" (part of the 
difficulty in that has been supporting a multitude of different use cases 
simultaneously).  You can read about that here:

http://docs.astropy.org/en/v1.0.4/modeling/models.html

In short, a *model set* is what you want if you want to represent a collection 
of independent instances of the same model with different parameters.

When you do:

     >>> Mtotal2 = mymodel(p1=[1,2],p2=[1,2], p3=[0,0])

What you're doing is making a model where the parameters have 2D vector values. 
  Which is fine, and will work in many cases.  But the components of these 
parameters are not necessarily linearly independent, depending on how the 
particular model is defined.

On the other hand, if you specify the n_models=2 argument like:

     >>> Mtotal2 = mymodel(p1=[1, 2], p2=[1, 2], p3=[0, 0], n_models=2)

This is like an array of models of whatever type "mymodel" is, with parameters 
(1, 1, 0) and (2, 2, 0).  In other words it's more like having a list of models 
like [mymodel(1, 1, 0), mymodel(2, 2, 0)].  When these are evaluated on an input 
the broadcasting is such that they are evaluated independently (but still 
generally more efficiently than if you made a list of models and evaluated them 
one by one--especially when you're scaling up to hundreds or thousands).

One important distinction to make though, is that if you have a model set, when 
you pass in an array-like argument the default assumption is that the first axis 
of that array corresponds to the number of models in the set, as though you are 
passing in a separate input to each model.  That is, the size of the first array 
axis is assumed to be the same as the number of models in the set.  So for 
example if I have:

 >>> m = Gaussian1D([1, 2], [0, 0], [0.1, 0.1], n_models=2)

and try to evaluate it on an array of shape (3,):

 >>> m([1, 2, 3])
ERROR: ValueError: Input argument 'x' does not have the correct dimensions in 
model_set_axis=0 for a model set with n_models=2. [astropy.modeling.core]

However, as explained in the docs, the way to say "no, I really just want to 
evaluate each model in the set on the same 3 element array" is:

 >>> m([.1, .2, .3], model_set_axis=False)
array([[ 0.60653066,  0.13533528,  0.011109  ],
        [ 1.21306132,  0.27067057,  0.02221799]])

There you see that you get a result of shape (2, 3).  The rows (2) correspond to 
each model in the set, while the columns correspond to the evaluation of that 
model on each of the inputs 0.1, 0.2, and 0.3.

If it helps clarify the difference, try also using a shape (2,) array.  This 
will work with or without model_set_axis=False, but with different results in 
each case.  Without you get:

 >>> m([.1, .2])
array([ 0.60653066,  0.27067057])

This is equivalent to passing 0.1 to Gaussian1D(1, 0, 0.1) and 0.2 to 
Gaussian1D(2, 0, 0.1).

On the other hand:

 >>> m([.1, .2], model_set_axis=False)
array([[ 0.60653066,  0.13533528],
        [ 1.21306132,  0.27067057]])

This passes the input [0.1, 0.2] to Gaussian1D(1, 0, 0.1), and the same to 
Gaussian1D(2, 0, 0.1).

For your summation example this latter usecase is probably what you want. 
Although I do hope to implement some tricks to speed up compound models where 
only a single operator is used (i.e. a reduction operation like sum or product), 
for now you'll get the best performance out of a model set, and then performing 
a reduction along the desired axis:

 >>> mymodel(p1=p1s, p2=p2s, p3=p3s, n_models=n_sources)
 >>> the_summed_sources = mymodel(input, 
model_set_axis=False).sum(axis=mymodel.model_set_axis)

(where above mymodel.model_set_axis is just 0 by default)

Erik