2D models created with units cannot be called with a dimensionless 0.0

[olebole]

Description

If an argument in a 2D model (that was created with units) is zero, one gets an exception.

Expected behavior

I would expect that a zero (0.0) does not need to be armed with a unit, since this is the case in other places as well, f.e.

>>> 1*u.arcsec + 0
<Quantity 1. arcsec>

Actual behavior

>>> from astropy.modeling.models import Gaussian2D
>>> import astropy.units as u
>>> gaus = Gaussian2D(x_stddev=1*u.arcsec, y_stddev=1*u.arcsec)
>>> gaus(0, 0)
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/astropy/units/quantity_helper/helpers.py", line 32, in get_converter
    scale = from_unit._to(to_unit)
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 950, in _to
    raise UnitConversionError(
astropy.units.core.UnitConversionError: 'Unit("1 / arcsec2")' is not a scaled version of 'Unit(dimensionless)'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/astropy/units/quantity_helper/helpers.py", line 147, in helper_dimensionless_to_dimensionless
    return ([get_converter(unit, dimensionless_unscaled)],
  File "/usr/lib/python3/dist-packages/astropy/units/quantity_helper/helpers.py", line 34, in get_converter
    return from_unit._apply_equivalencies(
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 886, in _apply_equivalencies
    raise UnitConversionError(
astropy.units.core.UnitConversionError: '1 / arcsec2' and '' (dimensionless) are not convertible

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/astropy/modeling/core.py", line 401, in __call__
    new_call = make_function_with_signature(
  File "/usr/lib/python3/dist-packages/astropy/modeling/core.py", line 380, in __call__
    return super(cls, self).__call__(*inputs, **kwargs)
  File "/usr/lib/python3/dist-packages/astropy/modeling/core.py", line 919, in __call__
    return generic_call(self, *new_args, **kwargs)
  File "/usr/lib/python3/dist-packages/astropy/modeling/core.py", line 4121, in generic_call
    outputs = self.evaluate(*chain(inputs, parameters))
  File "/usr/lib/python3/dist-packages/astropy/modeling/functional_models.py", line 374, in evaluate
    return amplitude * np.exp(-((a * xdiff ** 2) + (b * xdiff * ydiff) +
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 457, in __array_ufunc__
    converters, unit = converters_and_unit(function, method, *inputs)
  File "/usr/lib/python3/dist-packages/astropy/units/quantity_helper/converters.py", line 166, in converters_and_unit
    converters, result_unit = ufunc_helper(function, *units)
  File "/usr/lib/python3/dist-packages/astropy/units/quantity_helper/helpers.py", line 150, in helper_dimensionless_to_dimensionless
    raise UnitTypeError("Can only apply '{}' function to "
astropy.units.core.UnitTypeError: Can only apply 'exp' function to dimensionless quantities

While this is usually easy to avoid, it is somehow unexpected.

Steps to Reproduce

See above.

System Details

• Linux-5.10.0-2-amd64-x86_64-with-glibc2.31
• Python 3.9.1+ (default, Feb 5 2021, 13:46:56) [GCC 10.2.1 20210110]
• Numpy 1.19.4
• astropy 4.2
• Scipy 1.6.0
• Matplotlib 3.3.4

(all from Debian Testing)

[pllim]

a zero (0.0) does not need to be armed with a unit

Not necessarily. In magnitude, zero mag in one system isn't zero in another system.

[mhvk]

@olebole - the issue likely is that a unit is attached to the input because the model has units, and for a plain 0 that would be dimensionless. And 1 * u.arcsec + 0 * u.one also fails.

In principle, it may be possible to change this, though there is some advantage to being explicit (as @pllim notes, 0 does not always have obvious meaning). But I fear that I'm not well enough versed in astropy.model to know where to do that....

I should perhaps add that the main reason we allow 0 to have any unit is for comparisons: q < 0, etc.; we had quite a bit of discussion before agreeing that overall the benefits outweighed breaking "explicit is better than implicit".

[nden]

The intention with units of models is to be explicit. If a model was initialized with units it expects inputs with units.
I need to check if this is well documented though.
Perhaps we can close the issue?

[olebole]

@nden I think 0 is a special case, as usually the unit doesn't give additional information: 0" is the same as 0°, and I would rather support the "be strict in what you return, be generous in what you accept" rule. Requiring to be explicit in units when referring to the coordinate origin seems a bit too strict to me.

I'd think the question is rather how much effort it takes to solve this; I agree that it is a rather unimportant problem. But when it takes more time to document the behaviour properly than to fix it, maybe the fix is better?

[nden]

@olebole This has come up in different contexts (not specifically the case with 0 as input) - should one be able to fit a model with units to arrays or fit a model without units to quantities, and the corresponding evaluation questions. I think the general rule needs to be documented .
Are you defending the specific case that 0 should be a valid input to models with units and we should make a special case for it?

[olebole]

I think, 0 (and ±∞) is a special case, because it is not uncommon to use it without a unit. For other values, I think the unit should be mandatory, for a simple reason: For example, a Gaussian with σ=1nm is the same as a Gaussian with σ=10Å; so I would expect that they behave the same. The value of gaus(2.2) therefore would depend on the (physically irrelevant) detail with which unit gaus was created -- here explicit is definitely better than implicit. But for gaus(0.0) this does not matter.

[WilliamJamieson]

@nden I think 0 is a special case, as usually the unit doesn't give additional information: 0" is the same as 0°, and I would rather support the "be strict in what you return, be generous in what you accept" rule. Requiring to be explicit in units when referring to the coordinate origin seems a bit too strict to me.

The problem is that there are exceptions to 0 being a "special case", for example temperature and magnitude (as @pllim noted). Models (really ~astropy.units) has no way to distinguish if one is in an exception to this special case or not (indeed it is technically possible for someone to define a unit with similar properties). So one will get an UnitConversionError because there is no way for astropy to be certain that the "intended unit" of a dimensionless 0 will correspond to the 0 of the units involved in the model (0°C is not 0K so what does 0 in temperature "mean").

In most cases, the unit is dropped from a 0, as a notational connivence because we normally work with units which are all convertible to on-another with [linear maps)(https://en.wikipedia.org/wiki/Linear_map). This means 0 is always fixed (f(0*x) = 0*f(x) = 0). In reality, that 0 value and its unit cancel out because of this fact. However, this is not generally true.

As you do point out, in the majority of cases we do have this nice property for units so there are a few approaches to allowing a dimensionless 0 to be used in a model which defines units:

1. Just allow a dimensionless 0 input and return a dimensionless output (all the units will need to be dropped inside the model). In this case, modeling should raise a warning to the user that they should not be doing this.
2. Allow the user the ability to easily create a dimensionless copy of an existing model. Then the user can use a dimensionless 0.

I think the best thing to do is to leave models as they are and just add a note to the documentation. In any case this should be better documented.

[olebole]

After re-reading the stacktrace, I think the problem is not in the model, but in the units themself:

>>> import numpy as np
>>> import astropy.units as u
>>> np.exp(0)
1.0
>>> 1*u.arcsec**-2 + 0
<Quantity 1. 1 / arcsec2>
>>> np.exp(0*u.arcsec**-2)
Traceback (most recent call last):
  File "[…]/units/quantity_helper/helpers.py", line 32, in get_converter
    scale = from_unit._to(to_unit)
  File "[…]/units/core.py", line 950, in _to
    raise UnitConversionError(
astropy.units.core.UnitConversionError: 'Unit("1 / arcsec2")' is not a scaled version of 'Unit(dimensionless)'

During handling of the above exception, another exception occurred:
[…]

I think that the confusion is that some operations (like +) work with a unitless zero, while others don't. So, the unit seems to know by itself whether it is "linear" or not. This is however anyway broken:

>>> import numpy as np
>>> import astropy.units as u
>>> 1*u.Celsius + 22*u.Celsius
<Quantity 23. deg_C>
>>> 1*u.Celsius - 22*u.Celsius
<Quantity -21. deg_C>

The + shouldn't work; Celsius is not additive, and the difference should be Kelvin or just degrees.

[mhvk]

@olebole - the units system always respects units, so 0 arcsec is always an angle. The only exception is for 0 (or infinity) without a unit - it is treated as being able to have an arbitrary unit. I should perhaps note again that the main reason we implemented this is to allow comparisons, i.e., things like velocity < 0.