[AstroPy] Composite models or input arrays
Erik Bray
embray at stsci.edu
Tue Aug 25 11:17:21 EDT 2015
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
More information about the AstroPy
mailing list