ENH: support np.add, subtract, linspace for logarithmic units by mhvk · Pull Request #19853 · astropy/astropy

mhvk · 2026-06-02T11:12:34Z

As the title states, but note that for support to be possible, it was necessary to relax assumptions made previously about when a Magnitude is equivalent to a quantity with units of mag: before, this was the case if the physical unit was dimensionless, but now this is extended to whenever the physical unit is equivalent to dimensionless. I think in practice this is fine. EDIT: note that this also means it is a bit of an API change, now labelled as such!

This partially addresses #12548. As such, I guess one could call it a bugfix, but I think the impact goes beyond that, and there doesn't seem a pressing need.

By checking this box, the PR author has requested that maintainers do NOT use the "Squash and Merge" button. Maintainers should respect this when possible; however, the final decision is at the discretion of the maintainer that merges the PR.

github-actions · 2026-06-02T11:12:46Z

Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

Do the proposed changes actually accomplish desired goals?
Do the proposed changes follow the Astropy coding guidelines?
Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for rebase and squash.
Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the bot.
Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) before merge.

pllim

I am wary when it comes to doing math in the magnitude systems... Just some clarification questions below. Thanks!

pllim · 2026-06-02T14:59:09Z

+        lq = np.arange(1, 11.0) << u.mag("m/cm")
+        lq2 = u.Magnitude(2.0)
+        got = lq / lq2
+        exp = lq.to(u.mag) / lq2.to(u.mag)


This hurts my head lol. What does it mean physically when someone divides two magnitudes?

It doesn't mean much, but we've always had quantities where this is allowed: (10*u.mag) / (2*u.mag) == 5 and to match this before we already allowed (10*u.mag()) / (2*u.mag()) == 5. This PR extends this to also work for (5*u.mag("cm/m")) / (2*u.mag) == 5*u.one.

Now I totally agree neither of these seem all that logical, but if I wrote instead 10*u.mag == 5 * (2 * u.mag), it feels more reasonable. But it all boils down to whether one feels that mag is enough of a unit that in division one can cancel...

Note that I would do hope this will not be used, but the multiplication version is needed to get linspace to work...

pllim · 2026-06-02T15:05:28Z

+        lq2 = u.Magnitude(5.0, 2 * u.mag)
+        q1 = u.Quantity(lq1)
+        q2 = u.Quantity(lq2)
+        exp = func(lq1, lq2)


This concerns me a bit. How do we know both magnitudes are of the same zeropoints? Subtracting mags should be same as dividing linear flux, but I don't see that sort of checks are done in this test.

This only works for dimensionless (but possibly scaled) physical units. Otherwise, you cannot even make the Magnitude into a Quantity:

In [2]: u.Quantity(u.Magnitude(5.0, u.mag)) Out[2]: <Quantity 5. mag> In [3]: u.Quantity(u.Magnitude(5.0*u.m, u.mag)) UnitTypeError: Quantity instances require normal units, not <class 'astropy.units.function.logarithmic.MagUnit'> instances.

Aside: your question brings up a bug: the following should fail too:

In [4]: u.Quantity(u.Magnitude(5.0*u.m, 2*u.mag)) Out[4]: <Quantity -0.87371251 2 mag(m)>

(does not fail on current master either, sigh)

Note that I rather deal with the above crazy example separately... I raised #19867 to remind us...

pllim · 2026-06-02T15:05:51Z

        assert res.unit == u.mag()
+
+    def test_linspace(self):
+        res = np.linspace(0 * u.dex("AA"), 2 * u.dex("cm"), 11)


What does linear space mean in log units?

This gives

In [5]: np.linspace(0 * u.dex("AA"), 2 * u.dex("cm"), 11) Out[5]: <Dex [ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.] dex(Angstrom)>

which is sort-of nice to be able to do. In physical units it is logarithmic:

In [6]: np.linspace(0 * u.dex("AA"), 2 * u.dex("cm"), 11).physical Out[6]: <Quantity [1.e+00, 1.e+01, 1.e+02, 1.e+03, 1.e+04, 1.e+05, 1.e+06, 1.e+07, 1.e+08, 1.e+09, 1.e+10] Angstrom>

mhvk · 2026-06-02T15:53:53Z

@pllim - the whole thing was really weird to think through. To get the linspace example to work, we need the following to work:

start = 0*u.dex("AA")
stop = 2*u.dex("cm")
step = (stop - start) / 10
step
# <Dex 0.2 dex(cm(1/10) / Angstrom(1/10))>
offsets = np.arange(11) * step
offsets
# <Dex [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] dex>
result = offsets + start  # Actually done in-place in linspace
result[-1] = stop
result
# <Dex [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] dex(Angstrom)>

So, you see that step does exactly as you thought it should: subtraction of logarithmic quantities implies division of their physical units, and division by ten means raising the unit to the one-tenth power.

The thing tricky to see is offsets, but it becomes easier if you write it out,

np.arange(11) * (0.2 * u.dex("cm(1/10) / Angstrom(1/10)"))

Now multiplication with a number means multiplying the value and raising the physical unit to the corresponding power, so this should become a list like

[0 * dex("(cm/AA)(0/10)"), 0.2 * dex("(cm/AA)(1/10)"), 0.4 * dex("(cm/AA)(2/10)"), ..., 2 * dex("(cm/AA)(10/10)")]

As written, this is not possible to represent as a single Dex, since it is not all the same unit. What this PR does is note that the physical unit is equivalent to dimensionless, so if we first convert to that, we can do it:

(0.2 * u.dex("cm(1/10)/AA(1/10)")).to(u.dex())
# <Dex 1. dex>

After which the result is less surprising...

mhvk · 2026-06-02T15:54:24Z

p.s. Possibly I should add the above as an example to the documentation...

neutrinoceros

I'm also not very comfortable with log units arithmetics. I'll hang in there !
I got a couple questions as a warm up.

neutrinoceros · 2026-06-05T06:55:24Z

+        # as pure Quantity output?
+        if method == "__call__" and (op := ADD_SUBTRACT_UFUNCS.get(function)):
+            # Get output, if any.
+            out = kwargs.get("out", [None])[0]


I don't get why only the first element from out is used. Can you explain ?

__array_ufunc__ always passes out as a tuple with nout elements. In our case, we know there is only one output. I'll add a comment.

neutrinoceros · 2026-06-05T06:56:40Z

+            except TypeError:
+                # Operation on units failed: can happen if out got us here, but
+                # the inputs were not LogQuantity (so had no function units).
+                assert out is self


I don't get this assertion. Why is this expected ?

I put the assert here because the logic is not obvious: as the comment states, the unit operation can only lead to a TypeError if neither input unit was a LogQuantity (or, more specifically, neither unit supported addition or subtraction, which is only supported by LogUnit). But in that case, why are we in LogQuantity.__array_ufunc__? If not an input, the only option then is that out is a LogQuantity.

Does this help? Do you have a suggestion for how to make the comment clearer?

p.s. Note that the unit subtraction/addition can perhaps fail, but that would give a UnitsValueError, not a TypeError or subclass.

if not an input, the only option then is that out is a LogQuantity.

this does help a ton. ANd just to double check: do I understand correctly that __array_func__ is tried on every input as well as out with unchanged arguments (henc, out is self) ?

I think adding this explanation as an inline comment could make a life-saving difference for maintainability. I don't have a better suggestion unfortunately.

I added at least an example...

Sorry, I pushed before reading your replies. In the particular example of np.add(q1, q2, out=lq), we go directly to LogQuantity.__array_ufunc__ since LogQuantity is a subclass of Quantity (this follows the python logic for operators).

neutrinoceros · 2026-06-05T07:00:11Z

+            if result.unit != new_unit.function_unit:
+                # Ensure correct unit, as can be needed when the first
+                # argument is not a LogQuantity.
+                result <<= new_unit.function_unit


This looks symptomatic of a deeper problem, or maybe I'm just misunderstanding something, but I would not expect the result could not already have the correct unit at this point. Why is this allowed ?

The comment gave the clue, though not clear without an example.... I now have

# Ensure correct unit, as can be needed when the first # argument is not a LogQuantity, as in 1*u.dex - 10*u.mag().

For this, the call above, where the physical unit has been removed, will do 1*u.dex - 10 * u.mag, which gives -3 dex, but the unit we want is u.mag(), so we need to convert back to mag, 7.5 mag.

But I now realize I could just calculate converters, instead of letting Quantity.__array_ufunc__ take care, which is faster and probably clearer. So, let me do that...

sounds good. Please let me know when you're ready for another review !

This works by supporting np.add and np.subtract via the python operators and not falling back to regular Quantity unless truly needed.

With this, also allow output of addition or subtraction of regular quantities into logarithmic quantities, as well as the reverse (if the units allow it, of course). In the process, use the new implementation for the python operators (as is normally done for ndarray subclasses, but was not implemented yet here).

mhvk added this to the v8.1.0 milestone Jun 2, 2026

mhvk requested review from neutrinoceros and nstarman June 2, 2026 11:12

mhvk added units Enhancement labels Jun 2, 2026

mhvk force-pushed the units-logarithmic-add-subtract-ufuncs branch from eb4c16a to 35891bf Compare June 2, 2026 11:15

mhvk added the API change PRs and issues that change an existing API, possibly requiring a deprecation period label Jun 2, 2026

mhvk force-pushed the units-logarithmic-add-subtract-ufuncs branch from 35891bf to 7f81e0f Compare June 2, 2026 11:49

pllim reviewed Jun 2, 2026

View reviewed changes

neutrinoceros reviewed Jun 5, 2026

View reviewed changes

ENH: Enable np.add, np.subtract, and np.linspace.

311615a

This works by supporting np.add and np.subtract via the python operators and not falling back to regular Quantity unless truly needed.

mhvk force-pushed the units-logarithmic-add-subtract-ufuncs branch from 7f81e0f to b119b9b Compare June 5, 2026 09:51

mhvk added 2 commits June 5, 2026 12:02

DOC: changelog entry

6603b02

mhvk force-pushed the units-logarithmic-add-subtract-ufuncs branch from b119b9b to 6603b02 Compare June 5, 2026 10:02

Uh oh!

Conversation

mhvk commented Jun 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions Bot commented Jun 2, 2026

Uh oh!

pllim left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk commented Jun 2, 2026

Uh oh!

mhvk commented Jun 2, 2026

Uh oh!

neutrinoceros left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

neutrinoceros Jun 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

mhvk commented Jun 2, 2026 •

edited

Loading

neutrinoceros Jun 5, 2026 •

edited

Loading