I cannot judge this well enough to do a thorough review, but like where it is going very much. I'll think a bit more about the arbitrary unit -- probably more helpful!

@embray - can we remove the 1.1 milestone? (to try and do the feature freeze this week)

Don't think this can realistically go in 1.1, removing milestone.

This comment ignores all the detailed discussion between @mhvk and @embray and just answers in general to the "big block of text" on the top.

The concept looks workable to me. I'm not saying "good" because it's another layer of a fairly complex interface. I understand that in the most general case, defining a model is a complex undertaking that requires a complex interface. I assume that is the reason, why non of the other fitting packages that I know of has something like this. (Some domain specific packages like sherpa use some simpler internal tracking of units to distingiush counts/s/keV from counts/s/bin, but that's not really transparent for meaningful for the user.)

That said, this interface covers all the use cases I can think of for now, most importantly, it is opt-in for all models so far. Even going forward, I think it's important to be able to offer an opt-out of units option for all models, even those where a default units makes total sense (blackbody). Astropy is still unique in it's use of units and interaction of astropy.modelling with other packages requires that is all works on plain floats or numpy arrays. I understand that even with a requires unit I could always write a wrapper function to add the required units before I call the astropy model and strip them off again after, but I see that as deterring people form an otherwise nice framework.
That's not a pushback against this feature, just a feature request for an easy opt-out should there ever be models with hard required units.

@hamogu Thanks for the feedback. I forget I mention this anywhere else, but one thing I've been meaning to add to this is an option, defined at the class-level, for how strict it is about units. I think some models might absolutely require units, where others could be more flexible in what it allows. There were several possible settings for this. I could also be fine with a user-supplied argument though to completely ignore all units handling.

@embray - I'd be interested in helping make this happen. Are you still going to have time to work on this once you have adequate feedback? I'm going to review the concepts here and provide some feedback in my next comment. If you don't have time to pursue this further, or not enough time to do it all, I'm happy to help implement some of the remaining things.

@embray - with apologies for the delay, here is some detailed feedback on the high-level ideas here. I've tested out the PR but haven't looked at any code at this point. As mentioned in my previous comment, I'd be very interested in making sure this happens and can commit some time to working on it. If you no longer have time, let me know, and I can take over this branch and work from it - otherwise I can also do PRs to your branch if you want to share the work.

In some cases below, I'm just summarizing some of the changes with examples for anyone else who wants to review this (and also will help figure out if I'm misunderstanding anything)

Parameters with units

I think your approach to setting parameters to Quantities is good. Essentially, the units get set at initialization time:

In [1]: g = Gaussian1D(amplitude=1 * u.Jy)

In [2]: g
Out[2]: <Gaussian1D(amplitude=1.0 Jy, mean=0.0, stddev=1.0)>

It's then not possible to change the units afterwards:

In [3]: g.amplitude = 3 * u.joule

In [4]: g
Out[4]: <Gaussian1D(amplitude=3.0 Jy, mean=0.0, stddev=1.0)>

However, this does illustrate that we need better validation of input units when we set the parameters (the above should have raised an error). Assuming we would raise an error here, this seems fine to me.

In some cases, some units should be linked, so for example the units for mean and stddev should be the same. We'll need a way to specify that, since at the moment I can set inconsistent units:

In [5]: g = Gaussian1D(amplitude=1 * u.Jy, mean=2 * u.km, stddev=4 * u.s)

In [6]: g
Out[6]: <Gaussian1D(amplitude=1.0 Jy, mean=2.0 km, stddev=4.0 s)>

I do realize that this is just a WIP, and not expected to all work yet, but just using it as an illustration of what we need to consider. And you did mention this as a to-do:

More sophisticated unit specs for parameters--currently parameters can be given a required unit but don't support fancier unit specs like inputs and outputs do. For example, it might be useful to require one parameter's units to be equivalent to another parameter's units.

so +1 to this.

Any active equivalency can be taken into account, but what about also allowing equivalencies to be attached to parameters? For example, one could imagine wanting to have a spectral parameter which always includes the spectral equivalency - for example for a blackbody model but you wouldn't want the equivalency to apply to other parameters.

Evaluating models

As you mentioned, you can set input_units and output_units, which are set to the name of the parameter from which to use the units that any input should be in when evaluating the models. So in the case of the Gaussian:

In [7]: g.input_units
Out[7]: 'mean'

In [8]: g.output_units
Out[8]: 'amplitude'

This makes sense to me, and I think we should also not accept multiple incompatible units as input since it wouldn't make sense (and equivalencies take care of the cases where different units could make sense).

Evaluating a model using a value with the wrong units leads to an error:

In [9]: g = Gaussian1D(amplitude=1 * u.Jy, mean=2 * u.micron, stddev=0.1 * u.micron)

In [10]: g(3)
...
UnitConversionError: Units of input 'x', (dimensionless), could not be converted to required input units of micron (length)

whereas using the correct units or equivalent ones works:

In [11]: g(3 * u.nm)
Out[11]: <Quantity 2.5204894036926934e-87 Jy>

In [12]: u.set_enabled_equivalencies(u.spectral())
Out[12]: <astropy.units.core._UnitContext at 0x10aa63b38>

In [13]: g(300 * u.THz)
Out[13]: <Quantity 1.799785424849269e-22 Jy>

As mentioned above, the issue here is the enabled equivalencies apply to all parameters, whereas one might only want it to apply to one.

For more complex units, I like the ability to define more complex input/output units. I would personally prefer if we recommended that these were defined using properties instead of lambda functions, for readability, e.g.

@property
def input_units(self):
    return self.intercept.unit / self.slope.unit

instead of

input_units = lambda slope, intercept: intercept.unit / slope.unit

Fitting models with quantities

I would advocate actually not working on this as part of this PR but doing it in a separate PR. For now, it would be sufficient to return an error saying that fitting models with quantities is not yet supported. And it would then be trivial later to do another PR where if we pass a model with units to a fitter, things would indeed be converted to some common base units as you suggest.

Compound models

Creating compound models from models with units will place more constraints on what operations can be done on two models. This is a good thing. For example, it will not be possible to add two models together unless all their output_units are compatible.

I agree, this sounds like a good approach

Overriding input/output units

Possible support for per-instance unit spec overrides? Right now input_units and output_units are defined at the class-level. But it may also be useful to write code that defines these at the instance level. This would allow writing code that requires input to a particular model to always be of some unit, without having to subclass the model type.

This appears to already be possible?

In [14]: g.input_units = u.s

In [15]: g(3 * u.nm)
...
UnitConversionError: Units of input 'x', nm (length), could not be converted to required input units of s (time)

Equivalencies

Clearer support for unit equivalencies. Right now equivalencies will just work both when instantiating and evaluating a model. However, it might be useful to have equivalencies baked into a model somehow. For example, when evaluating a model it might be good to automatically apply all equivalencies that were active when the model was instantiated.

Agree - and as mentioned above, equivalencies should be specifiable at the parameter level.

Parameter.value

Open question: Should Parameter.value return a Quantity or the raw value? Currently something like model.amplitude.value or model.amplitude.default does not return a Quantity. Instead the Parameter object is a Quantity-like object in that it has a .unit attribute (if units are defined on a parameter). Sometimes it might be convenient to get a Quantity, but on the other hand the current interface feels more consistent with the existing Quantity interface.

I think .value has to return the value (and .unit the unit) otherwise it will be confusing. We can add .quantity to return it as a quantity?

@astrofrog , officially @embray does not have any time to work on this anymore, so I would say you can "take over" but give him credit to the existing work so far.

@astrofrog

However, this does illustrate that we need better validation of input units when we set the parameters (the above should have raised an error). Assuming we would raise an error here, this seems fine to me.

Yep, that should be fixed (I'm surprised it allowed that...?)

In some cases, some units should be linked, so for example the units for mean and stddev should be the same. We'll need a way to specify that, since at the moment I can set inconsistent units:

Yes, this was on my to do list. It's currently possible to link the units of outputs to parameters and vice-versa. But it should also be possible to link the units of individual parameters to each other via some kind of unit constraint specified when defining the parameters.

@astrofrog

Any active equivalency can be taken into account, but what about also allowing equivalencies to be attached to parameters?

👍

@astrofrog

I would personally prefer if we recommended that these were defined using properties instead of lambda functions

Using a a property doesn't really work in the context of the current implementation, which IIRC wants to sometimes check unit consistency in cases where a self might not exist. That said there's no reason to use property, and writing

def input_inputs(...):
    ...

is no different than using input_units = lambda ... You don't even have to put @staticmethod if you don't want to since the metaclass knows what to do with this, but maybe putting @staticmethod also makes it more clear that it is not an instance method. That should be allowed if it isn't already.

I prefer the lambda expression for short expressions because it makes clear that input_units is generally treated as a "class attribute" even though it can also be a function. A strange but true fact.

@astrofrog

I think .value has to return the value (and .unit the unit) otherwise it will be confusing. We can add .quantity to return it as a quantity?

I think that makes some sense.

Fitting could be added now, but it would be very slow. My idea for fitting would be to perform all unit conversions first, and then do the fit using the bare float values, and only restore units at the end. Otherwise you'd have to strip units before putting them through the optimizer, then somehow re-attach the units when the models evaluate is called to update the residuals, and it will be quite a slog.

I could try rebasing this if you want, but as @pllim said I don't have time to otherwise do work on it right now.

@embray

Fitting could be added now, but it would be very slow. My idea for fitting would be to perform all unit conversions first, and then do the fit using the bare float values, and only restore units at the end. Otherwise you'd have to strip units before putting them through the optimizer, then somehow re-attach the units when the models evaluate is called to update the residuals, and it will be quite a slog.

I agree, this makes sense. Maybe one could imagine having a context manager internally that temporarily makes evaluate return plain values in a specified set of units, and we can use that context manager inside the fitters?

I could try rebasing this if you want, but as @pllim said I don't have time to otherwise do work on it right now.

Shall I just take over this branch and open another pull request to continue it, since I'll have time to work on it over the next couple of weeks?

Shall I just take over this branch and open another pull request to continue it, since I'll have time to work on it over the next couple of weeks?

👍

All right, in that case I am now working on this branch, and will open a new pull request with additional changes!

@astrofrog is working on this over at #4615, so closing this one.