np.arange and np.linspace do not work properly with log units

[olebole]

Description

When trying to create an array that is equally spaced in log units, both attempts fail:

linspace converts to a simple (non-logarithmic) quantity with dimensionless dex:

>>> import numpy as np
>>> import astropy.units as u
>>> a = np.linspace(0*u.dex(), 1*u.dex(), 5)
>>> a
<Quantity [0.  , 0.25, 0.5 , 0.75, 1.  ] dex>
>>> a.physical
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 1017, in __getattr__
    raise AttributeError(
AttributeError: 'Quantity' object has no 'physical' member

When applying a unit to dex, already linspace fails:

>>> a = np.linspace(3*u.dex(u.Angstrom), 4*u.dex(u.Angstrom), 5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<__array_function__ internals>", line 5, in linspace
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 1698, in __array_function__
    return super().__array_function__(function, types, args, kwargs)
  File "/usr/lib/python3/dist-packages/numpy/core/function_base.py", line 120, in linspace
    start = asanyarray(start) * 1.0
  File "/usr/lib/python3/dist-packages/astropy/units/function/core.py", line 582, in __mul__
    raise UnitTypeError("Cannot multiply function quantities which "
astropy.units.core.UnitTypeError: Cannot multiply function quantities which are not dimensionless with anything.

arange does something completely weird without units:

>>> np.arange(0*u.dex(), 1*u.dex(), 0.25*u.dex())
array([1.        , 1.77827941, 2.55655882, 3.33483823])

and with units, arange raises an exception as well:

>>> np.arange(3*u.dex(u.Angstrom), 4*u.dex(u.Angstrom), 0.25*u.dex())
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 893, in to_value
    scale = self.unit._to(unit)
AttributeError: 'DexUnit' object has no attribute '_to'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 1263, in __float__
    return float(self.to_value(dimensionless_unscaled))
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 896, in to_value
    value = self._to_value(unit, equivalencies)
  File "/usr/lib/python3/dist-packages/astropy/units/quantity.py", line 803, in _to_value
    return self.unit.to(unit, self.view(np.ndarray),
  File "/usr/lib/python3/dist-packages/astropy/units/function/core.py", line 256, in to
    return self.physical_unit.to(other, self.to_physical(value),
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 1135, in to
    return self._get_converter(Unit(other),
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 1066, in _get_converter
    raise exc
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 1051, in _get_converter
    return self._apply_equivalencies(
  File "/usr/lib/python3/dist-packages/astropy/units/core.py", line 1027, in _apply_equivalencies
    raise UnitConversionError(
astropy.units.core.UnitConversionError: 'Angstrom' (length) 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/units/quantity.py", line 1265, in __float__
    raise TypeError('only dimensionless scalar quantities can be '
TypeError: only dimensionless scalar quantities can be converted to Python scalars

Expected behavior

There should be a simple way to create arrays that are equally spaced in log units, or at least NotImplemented should be raised.

System Details

Linux-5.14.0-2-amd64-x86_64-with-glibc2.32
Python 3.9.9 (main, Nov 16 2021, 10:24:31)
[GCC 11.2.0]
Numpy 1.19.5
pyerfa 2.0.0.1
astropy 5.0
Scipy 1.7.1

[pllim]

Is it related to this known issue: https://docs.astropy.org/en/latest/known_issues.html#numpy-array-creation-functions-cannot-be-used-to-initialize-quantity ?

[olebole]

Partially, The docs state that linspace should work.

[pllim]

Either the doc needs updating or this is some new compatibility problem... Hopefully @mhvk et al. can advise.

[olebole]

Maybe wait until I could re-check this once we have a numpy 1.21 in Debian testing.

[pllim]

Sounds good. I added a Needs-clarification label to indicate we await your re-checking. 😸

[mhvk]

No need to wait, this does not depend on numpy, but is currently not supported. For linspace, it could be made to work, as I think the intent is quite obvious, though that then needs a custom override in quantity_helpers/function_helpers.py as right we just let numpy do the work -- the reason that that fails for log units is that internally something like start + (stop-start)*i/n is done, and multiplication with log units is tricky: it should raise the unit to some power. The reason that the first example, with u.dex(), does give a result, is that the dex unit and a functional unit are not quite the same. I think what you want here is u.dex(u.one) -- and arguably it is a bug that that is not the same as u.dex()...

For arange, overriding is harder as the inputs are normally not arrays and thus not checked for overrides (I think it is possible in principle to support np.arange(q1, q2, n, like=u.Quantity), but that has not been implemented yet).

More relevantly perhaps, what does work is

In [12]: np.logspace(0.*u.dex(u.one), 5.*u.dex(u.one), 5)
Out[12]: 
<Quantity [1.00000000e+00, 1.77827941e+01, 3.16227766e+02, 5.62341325e+03,
           1.00000000e+05]>

[mhvk]

Actually, my diagnosis of the np.linspace problem is wrong; this should work, it is just that numpy decides to multiply start and stop with 1.0 - I can see if special-casing that multiplication would make things work...

[olebole]

@pllim Just for completeness: The behavior didn't change with numpy 1.21.4.
@mhvk somehow the multiplication is also weird for unitless dex:

>>> a = 1 * u.dex()
>>> a
<Dex 1. dex>
>>> 1 * a
<Quantity 1. dex>
>>> 2 * a
<Quantity 2. dex>
>>> a + a
<Dex 2. dex>

i.e. multiplication of u.Dex with a number returns a u.Quantity.

[mhvk]

@olebole - the above are examples of what makes it hard to get the logarithmic units right. If they have a unit inside, multiplication and division make little sense, yet when dimensionless, things can work (like dB/m or mag/day). The only way I could see to solve this is to disallow multiplication/division for logarithmic quantities, but to downcast to Quantity if they were dimensionless. I'm not sure it is totally ideal, and obviously suggestions for different behaviour are very welcome!

[olebole]

Why don't make multiplication and division sense? I would formally expect, that the factor goes into the power of the physical unit:

>>> u.Dex(10 * u.Angstrom)
<Dex 1. dex(Angstrom)>
>>> u.Dex(10 * u.Angstrom) + u.Dex(10 * u.Angstrom)
<Dex 2. dex(Angstrom2)>  # This works
>>> 2 * u.Dex(10 * u.Angstrom)
# Should also be:
<Dex 2. dex(Angstrom2)>

This should have a generic solution, as the "downcasting" to Quantities is IMO ugly here:

>>> u.Dex(10*u.one)
<Dex 1. dex>
>>> u.Dex(10*u.one).physical
<Quantity 10.>
>>> u.Dex(10*u.one) + u.Dex(10*u.one)
<Dex 2. dex>
>>> (u.Dex(10*u.one) + u.Dex(10*u.one)).physical
<Quantity 100.>
>>> 2 * u.Dex(10*u.one)
<Quantity 2. dex>
>>> (2 * u.Dex(10*u.one)).physical
AttributeError: 'Quantity' object has no 'physical' member

This all works already for addition and subtraction; IMO it should just be extended to multiplication/division.
One could f.e. just add multiplication/division operators to LogUnit (maybe with a check that other are plain numbers or so):

    def __mul__(self, other):
        return self._copy(self.physical_unit**other)
        
    def __rmul__(self, other):
        return self._copy(self.physical_unit**other)

    def __div__(self, other):
        return self._copy(self.physical_unit**(1/other))

and proagate this through LogQuantity. However LogQuantity comes with the note (from you, 2014)

astropy/astropy/units/function/logarithmic.py

Lines 315 to 316 in ce74be8

	# Could add __mul__ and __div__ and try interpreting other as a power,
	# but this seems just too error-prone.

Do you have a case where this would lead to error? I would see the only major problem that multiplication is not associatve with units:

(2 ⋅ 3) dex(Å) ≠ 2 ⋅ (3 dex(Å))

but this seems obvious to me.

[olebole]

Would it be useful if I would just open a PR to play with multiplication of log units?

[mhvk]

@olebole - I think in my note I may have partially worried exactly about associativity, but perhaps you are right that it is fairly obvious (or at least something the documentation can explain). The other possible issue is with other being an array, but I think you are correct that really that should not preclude trying to support plain numbers - it should actually not be too difficult, since if one were to raise the unit to something like an array, that would already give an exception.

Anyway, if you are up for a PR, that would definitely be welcome!

[mhvk]

p.s. Just should add I'm really happy to see this used!

[mhvk]

Some more progress towards this in #12579