Models should support Quantity

@nden

I just noticed the following behavior:

In [1]: from astropy.modeling import models
In [2]: m = models.Linear1DModel(slope=1,intercept=0)
In [3]: m([1,2.] * u.m)
Out[3]: array([ 1.,  2.])

even though you give a Quantity in, a regular array comes out.

Also,

In [28]: m = models.Linear1DModel(slope=1,intercept=2 * u.m)
InputParameterError: Don't know how to convert parameter of <class 'astropy.units.quantity.Quantity'> to float

So even if I wanted to use quantities as the variable in the model, I can't do it because the parameters are forced to be a float.

So it would be much more valuable if parameters could be a quantity where q.isscalar == True, and the varaible in the model can be any sort of quantity.

cc @nden @adrn @astrofrog

@cdeil

I have a "toy" implementation of this issue but would like to get others' opinion on how this is supposed to work.

More specifically, I think Parameter should not subclass Quantity but rather models should be able to accept as input parameters instances of Quantity. Models which can accept quantities need to define a default_unit and default_equivalencies for their parameters. Passing quantities as parameters when initializing a model should be optional, i.e. models should be able to work with numbers. If quantities are used internally the model should work in units of default_unit.

When a model is initialized with quantities, it will return a quantity when evaluated, when it is initialized with numbers, it will have the current behaviour. For this reason it is necessary for models to define output_units.

How does this sound?
@cdeil @plim @adrn @wkerzendorf

@nden

@nden , could it be easier on your life if modeling converts everything to Quantity at the beginning, using default_unit (which can be unitless) if only numbers are given? That way, you do not have to keep doing things like

if isinstance(x, u.Quantity):
    y = x.value
else:
    y = x

I find that I have to do that in pysynphot for functions that are supposed to support both normal arrays and Quantity arrays. Although I do not know of the performance impact if your calculations are done with Quantity instead of normal arrays. Also, Quantity does not yet support all math functions, .e.g, log.

@plim

@plim I know quantities support most math operations but within modeling/fitting it's much simpler to work with numbers. So right at the beginning I convert everything to numbers (after unit conversion) and at the very end attach output_units to the final answer (if the model was initialised with quantities). In addition, since it is not possible to add a number to a quantity, looking up the value will be necessary anyway.

Also we should consider what this means for composite models.

The one good thing working in Quantity is that it keeps track of the units for you. In your implementation, what happens when user gives output_units that does not make sense given the input Quantity?

output_units is not a user argument. It is defined in the model and is consistent with default_unit. If the user prefers different output units, they will have to do result.to(units.MyUinit)

Ah, okay. 👍

I know quantities support most math operations but within modeling/fitting it's much simpler to work with numbers.

If so, I would consider it undesired behaviour in Quantity! My sense would be to approach it from the opposite end: convert both parameters and inputs to Quantity automatically; with that, outputs would automatically be Quantities too.

Also, Quantity does not yet support all math functions, .e.g, log.

It does, as long as the units are dimensionless; are there models already where one wants log of a dimensional quantity?

We will never know what people need modeling for. But we cannot tell them not to use it because log does not work, or something like that.

I'd think it is desirable to impose astrophysically meaningful restrictions, as they catch unclear thinking and/or bugs. More practically, why not implement something that is physically meaningful and see whether problems arise? (For the particular case of log, the problem will go away with the introduction of functional quantities anyway -- which I hope to have a PR for soon after 0.3 is released).

@nden

Using Quantity internally for models wouldn't work in the current implementation. This is because the float values of multiple parameters are often combined in a single array for use with the fitting algorithms. There's no sensible way to do this with Quantity. Instead, parameters should only be returned as Quantities to the user when they do something like model.theta (where theta is some parameter).

The only issue I have with this currently is that right now the Parameter class, while neither an array nor a number, does act as a proxy that encapsulates info about the parameter, but can still be used in expressions like an array or number:

>>> g = models.Gaussian1D(1, 0, 0.5)
>>> g.amplitude
Parameter('amplitude', value=1.0)
>>> type(g.amplitude)
<class 'astropy.modeling.parameters.Parameter'>
>>> g.amplitude * np.array([1, 2, 3, 4])
array([ 1.,  2.,  3.,  4.])

This works even if there are multiple parameter sets, in which case each parameter is returned as an array:

g = models.Gaussian1D([1, 2], [0, 2], [0.5, 0.1])
>>> g.mean
Parameter('mean', value=array([ 0.,  2.]))
>>> g.mean * 3
array([ 0.,  6.])

But you can't do something like g.mean.mean() because it does not proxy all of the ndarray methods or attributes. You can however do:

>>> g.mean.value.mean()
1.0

Parameter.value gives direct access to the underlying data object. Maybe that would be sufficient for Quantities. If nothing else Parameter could maybe grow a unit attribute, or just use default_unit as @nden suggested.

If default_unit is unspecified then a parameter would just use the input quantity's unit as the "default" default_unit. Obviously any attempt to reassign the value of that parameter would have to have compatible units. Makes sense to me...

#1796 shows what I am proposing.

To write an astrophysical model one has to define default_unit (and possibly default_equivalencies) for each parameter, together with output_units. It imposes restrictions if quantities are passed but does not prevent one from using simply numbers.

A couple of silly examples to look at are the AZP projection:

https://github.com/nden/astropy/blob/720a4b17027107ed080e4d5f1ba848d8e5a4531a/astropy/modeling/projections.py#L73

or the silly implementation of Gaussian1D.
The examples are not particularly meaningful but test the functionality.

@eteq

I would like to raise again the question about using quantities within models and understand how others see this working.

To summarize what was proposed earlier in the PR:

add default_unit attribute to the Parameter class
add default_input_unit and default_output_unit attributes to the Model class
when initializing a model with quantities, convert them to the default units of the parameters
if input is quantities, convert it to default_input_unit and return the result as a quantity with default_output_unit
internally models work with numbers

This scheme makes sense for a model like black body but I am not sure how it will work correctly for any of the current models in modeling. If we take as an example Gaussian1D, there is no one set of "correct" units with which to define the model since gaussians may be used for many types of data.

Now let's suppose we take a different approach and instead of defining default units when the model class is written, we allow units to be defined when the model is initialized. Then models will have to make sure their parameter units are consistent and their input and output units are consistent with their parameter units. Is this always possible though?

I am leaning towards having the above scheme implemented but using it only for well defined physical models, like black body. Is there a better solution? Can someone give examples of other models that can take advantage of the above scheme?

@eteq @astrofrog @mhvk @embray @perrygreenfield @cdeil @pllim

@nden

@nden et al., for those without well-defined units, can we just use u.dimensionless_unscaled as defaults? Users can always change it to something meaningful that suit their science via keywords to constructor, or something like that.

For the current models in astropy I think the only cases where quantities could be used are angle parameters like e.g. theta in Gaussian2D or the parameters of RotateCelestial2Native.

I'm not sure the API and / or implementation would simplify if quantities were used for those. E.g. we probably soon want to have covariance matrices to represent fit errors, and if some parameters have units and others don't or different parameters have different units, what type of object is the covariance matrix then?

I think we should be a bit careful where to use quantities in modeling. If we seriously consider turning all parameters into quantities we would probably need a prototype to see if that works out well for the API and implementation.

@nden

@nden - thanks for the great summary. I think your example of the Gaussian is a very good one, in that, yes, many models typically will have units, but, no, these will not be known until initialisation.

For those types of models, I think the only sensible solution is to write them such that no conversion of quantity at all is done, i.e., if one initializes the Gaussian with a sigma with some unit, the user has to make sure the independent variable has the same unit. Since quantities now work just like numbers, this means the code for Gaussian would remain trivially simple.

Overall, two solutions would seem possible.

Along the lines you suggest, for parameters where the physical meaning is obvious beforehand, you have a default_unit=<something> while for ones where it only becomes clear at initialisation, one may have default_unit=None. For the implementation, though, my suggestion would not be to do the equivalent of real_par = u.Quantity(input, default_unit).value, but rather stick with the quantity in the code (i.e., omit the .value), at least down to the deepest level where one can. I am now so used to unit conversions catching mistakes, that I'd really love for it to be used consistently everywhere!
Be bold, and go Quantity all the way, i.e., models that need something with a given physical unit, should accept only quantities. With that, the whole concept of default_unit can be abandoned.

Now, @cdeil is correct that some aspects are less easily captured using quantities (polynomial coefficients, covariance matrices). Mostly, I think this can be captured by not insisting that there are together in a single ndarray (for which I think there is, in principle, little need).

@embray

@embray - Is the actual topic of this issue of your roadmap of modeling work you are planning? That is, a lot was discussed here, but the key feature is to let models accept Quantity parameters and inputs/outputs.

I don't think this will get resolved before 1.0, so removing milestone.

🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧🚧
^Big line of barriers to make it easier to find the resumed discussion on this already quite long issue.

I'm working on this now and hope to have an early prototype soon. Below I'm going to try to address some of the questions/comments from the earlier discussion on this issue.

One point I want to make early on is that the first prototype probably won't address every possible issue that one could raise with this, and that's fine. I think the most important thing will be to get something working that can evaluate models with quantities. I believe it will be important for Models to continue to work either with or without Quantities. To be clear, if Quantities are used at all (that as, any one of the inputs or parameters have units) then it will be required that all inputs and parameters are Quantities (even if their unit is 'dimensionless')--this is necessary to have any hope for consistency. Allowing both is important for at least three related reasons:

Backwards compatibility must be maintained.
Not all use cases call for Quantities.
For corner cases where Quantities don't immediately make sense to use, we can still fall back on purely numeric results.

@mhvk

Okay, I'm going to respond to a few comments from the original discussion (before things got side-tracked) which I have now re-re-re-reviewed. First up:

@mhvk in #1464 (comment)

I know quantities support most math operations but within modeling/fitting it's much simpler to work with numbers.

If so, I would consider it undesired behaviour in Quantity! My sense would be to approach it from the opposite end: convert both parameters and inputs to Quantity automatically; with that, outputs would automatically be Quantities too.

I'm still trying to figure out exactly what to do here. I want to actually avoid using Quantity arithmetic for all evaluation of the models, as this could get unbearably slow for fitting. On the other hand I also want to be able to determine the correct output units for some evaluation (where the required units haven't already been fixed), so I think it will also be necessary to pass the values as Quantities through the evaluate function to get the correct units. Perhaps this could be the default for simple evaluation, while for fitting we can drop the units. Note that in either case quantities will all be converted to the same set of base units for internal storage in the model's parameters array, so the arithmetic will work out correctly when calling the evaluate without units. On the other hand, maybe I could go ahead and enable different storage mechanisms for different parameters. For example, just keep each parameter internally as a Quantity to start with, but convert to the flattened array storage for fitting purposes. That isn't my initial plan and is not something I've tried yet, but we could resort to that if the benefits are clear.

In short, there needs to be some balance found between keeping simple one-time evaluation efficient, and also making evaluation in loops efficient. That will take some playing around.

In #1464 (comment)

I think your example of the Gaussian is a very good one, in that, yes, many models typically will have units, but, no, these will not be known until initialisation.

For those types of models, I think the only sensible solution is to write them such that no conversion of quantity at all is done, i.e., if one initializes the Gaussian with a sigma with some unit, the user has to make sure the independent variable has the same unit. Since quantities now work just like numbers, this means the code for Gaussian would remain trivially simple.

Some details aside, that's still the idea. Little to nothing should have to be changed in the evaluate methods in order to support units. While there are some things we'll be able to do in terms of validation, at the end of the day it will be up to users to make sure their inputs are consistent. Maybe in the future we could add some additional constraints to parameters regarding units. For example we could have some way of specifying (on Gaussian1D, for example) that we know a priori that the independent variable, mean, and stddev must all have the same (or compatible, up to equivalencies) units. This won't be in the initial prototype though, unless it turns out to be more important than I think it will be.

@nden

Responding to myself in #1464 (comment)

The only issue I have with this currently is that right now the Parameter class, while neither an array nor a number, does act as a proxy that encapsulates info about the parameter, but can still be used in expressions like an array or number. ... But you can't do something like g.mean.mean() because it does not proxy all of the ndarray methods or attributes. Parameter.value gives direct access to the underlying data object. Maybe that would be sufficient for Quantities. If nothing else Parameter could maybe grow a unit attribute, or just use default_unit as @nden suggested.

It seems like this will work for all Quantity arithmetic now, so long as the Parameter has a .unit attribute.

And relatedly, to @eteq in #1464 (comment)

So I'm willing to accept @mhvk's idea of QuantityParameter or similar, but only if there's a clear reason why it won't work to just use Quantity all the time, by just auto-converting to a dimensionless Quantity if you do something like Gaussian1D(amplitude=1, stddev=5, mean=6563). It seems like that addresses the problem of how to deal with something like polynomials - as long as the input is dimensionless (or a raw array), the output will be dimensionless. So then it seems like we get the best of both worlds.

Converting everything to dimensionless Quantities, in all cases, is a bit too much trouble, at least for now. However, when all parameters are explicitly given as Quantities, it should be possible to use those parameters in any arithmetic with other Quantities (since the parameters will have .unit attributes). It might not work transparently equivalently to a Quantity in all cases as it won't have all the ndarray methods. But if we really want to we could maybe have it proxy any ndarray/Quantity methods where applicable.

@cdeil

To @cdeil in #1464 (comment)

I'm not sure the API and / or implementation would simplify if quantities were used for those. E.g. we probably soon want to have covariance matrices to represent fit errors, and if some parameters have units and others don't or different parameters have different units, what type of object is the covariance matrix then?

Yes, if you need a covariance matrix that might be one of those cases where one should not use Quantities. However, this might not be entirely meaningless. As I responded to @mhvk above, all units will be converted to a common base (if possible--if not I'm not sure what the best action is) for fitting. So for something like a covariance matrix we could have an auxiliary list somewhere specifying what units were used for each of the variables.

Somewhat beyond the scope of models+Quantities, but this leads me to think it might be useful to have some kind of array type in which each element can have its own unit. This would not be terribly efficient or useful for vectorized operations, but it might be convenient just as a user interface. This would not be entirely unlike an interface I've already been working on for working with Model parameters more conveniently, though I've set that work aside for now.

@astrofrog

And to @astrofrog in #1464 (comment)

The model requires a parameter in certain units - for example a position angle in degrees. Then when it defines the Parameter it can specify unit=u.deg, and when someone sets the parameter, whatever is passed, it is converted to a float value in degrees. If the model and the user access the parameter value differently, we can make it return the value for the model, and the quantity for the user.

This speaks to a point I've been struggling with a little bit. Right now there are some models that require inputs in degrees (the namely the rotation models, and other models that have a rotation parameter like Gaussian2D. There's nothing special about degrees in this example; it's just one that's already applicable to many models. Right now, the models that require values in degrees will take just any array of floats and assume they have been correctly specified in degrees. For the sake of backwards compatibility, I want to keep this functionality for now, at least on models that already use this. However, in the long term, I want to change this so that if a model specifies unit=u.degree for some Parameter, then the user must supply a Quantity for that parameter, with units compatible with degrees. It might also be important to specify on some parameters what base unit that parameter should be converted to for evaluation. For degrees we might want to always store in radians, for example. It may be necessary in some cases to make this explicit; not sure yet though.

The model does not require a parameter in certain units - in that case, the Parameter class can just split up the value and the unit, and again when the model requests the value, the value is given, and when the user requests it, the quantity is returned.

Not always so simple. If the model is supplied parameters in different, but equivalent units then we still need to convert them to common units before storing them away so that evaluation can proceed correctly. But the original units could still be saved too, and returned to that unit when the user requests the parameter.

@embray

@embray - great to see your thoughts! Two answers and one comment:

It seems like this will work for all Quantity arithmetic now, so long as the Parameter has a .unit attribute.

Yay to ducktyping!

Somewhat beyond the scope of models+Quantities, but this leads me to think it might be useful to have some kind of array type in which each element can have its own unit.

Maybe just an ndarray with `dtype='O'? I know this is possible because annoyingly, currently,

In [2]: np.ma.array([1., 2.]) * u.m
Out[2]: 
masked_array(data = [<Quantity 1.0 m> <Quantity 2.0 m>],
             mask = False,
       fill_value = ?)

Finally, a comment on speed: with __numpy_ufunc__, speed differences will become substantially less since each calculation will just be done once (instead of twice, once with the wrong units, and once with the right ones). But like you, I hope that in the end so much of the code can be agnostic about whether or not input is in quantities or in numbers that it will not be too hard to support both explicitly, in which case pure arrays can still be used for cases where speed is of the essence.

@mhvk

@mhvk Yeah, an object array might be fine for that, actually. I was still thinking it might be nice to have a special array subclass that, if nothing else, presents an array of mixed-unit Quantities more nicely (cut out all <Quantity ...> noise, so that you just have something like [1.0 m, 2.0 m, 0.3 Jy] or whatever.

@astrofrog

@astrofrog , can we close this now?

Navigation Menu

Repository navigation

Models should support Quantity #1464

48 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Models should support Quantity #1464

Description

Activity

nden commented on Nov 15, 2013

pllim commented on Nov 15, 2013

nden commented on Nov 15, 2013

pllim commented on Nov 15, 2013

nden commented on Nov 15, 2013

pllim commented on Nov 15, 2013

mhvk commented on Nov 15, 2013

pllim commented on Nov 15, 2013

mhvk commented on Nov 15, 2013

embray commented on Nov 15, 2013

nden commented on Nov 15, 2013

nden commented on Feb 25, 2014

pllim commented on Feb 25, 2014

cdeil commented on Feb 25, 2014

mhvk commented on Feb 26, 2014

48 remaining items

eteq commented on Aug 8, 2014

astrofrog commented on Jan 14, 2015

embray commented on May 14, 2015

embray commented on May 14, 2015

embray commented on May 14, 2015

embray commented on May 14, 2015

embray commented on May 14, 2015

mhvk commented on May 15, 2015

embray commented on May 15, 2015

pllim commented on Jun 15, 2017

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions