PERF: use scaler class instead of lambda for unit scale conversions - #20183
PERF: use scaler class instead of lambda for unit scale conversions#20183mhvk wants to merge 33 commits into
Conversation
|
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.
|
|
p.s. Forgot to mention: there is another part still to come, that helps with larger arrays (>10**5 elements), by allowing numpy to do the scaling as part of the iterator. Since that avoids a memory copy, that has quite a big impact for those larger arrays (and is somewhat helped by this PR). |
|
I'll move this to draft while I try to figure out why builds fail on older python and windows/mac/aarch64... |
743cb88 to
6b35323
Compare
This is slower than lambda, so not great.
Also ensures that we do not promote complex64 to complex128.
This means we can also get rid of the fake imaginary components.
Agreed. I guess targeting 8.1, like this? |
wait the phrasing is either misleading or revealing a misunderstanding; this isn't a choice we can make at the module level: if a single module in astropy isn't compliant, we cannot ship |
| #include <stddef.h> // for offsetof() | ||
|
|
||
| // Once we're at 3.12, move SCALER_TP_FLAGS to its use, and remove this whole block. | ||
| #if Py_LIMITED_API + 0 >= 0x030C0000 |
There was a problem hiding this comment.
I saw that done in CPython, but to be honest have no idea why. Maybe simply if someone by mistake simply does #define Py_LIMITED_API? Which I guess would only be relevant for header files, not C code...
| sc1 = Scaler(10.0) | ||
| sc1_2 = Scaler(10.0) | ||
| assert sc1_2 == sc1 | ||
| assert not (sc1_2 != sc1) # noqa: SIM202 |
There was a problem hiding this comment.
possibly naive question: isn't __ne__ automatically derived for any object subclass ?
There was a problem hiding this comment.
I'm not sure that is the case in C -- I get passed Py_EQ or Py_NE and deal with them more or less separately.
A bit of both: I started not following the limited API, but then remembered vaguely that it was good to do so, so followed it. This is partially (but not only) why there is such a large hack for python 3.11. In the next PR, I fear the need to work around a numpy issue that was only solved in 2.6 (i.e., -dev) means I have to calculate offsets in structures... But I do like the idea of the limited API very much, in particular that one can have just one package for all supported versions of python and numpy -- saves a lot of compute, bandwidth and storage if nothing else! |
6fa80c8 to
be73429
Compare
|
Since I'm only middling at C, I asked Claude Opus 5 to It found the bump to numpy v2.0 :) Click to expandBug review — astropy PR #20183PR: PERF: use scaler class instead of lambda for unit scale conversions Method: the PR's C sources were checked out and built as a standalone extension Blockers1. Use-after-free in
|
| input | old (scale * value) |
new |
|---|---|---|
np.arange(3, dtype='f2') |
float16 |
float64 |
np.float16(1) |
float16 |
float16 |
Note the array and scalar paths now also disagree with each other: f2 scalars still return f16
because they take the O_factor Python-float path.
4. Complex arrays give different answers depending on memory layout
astropy/units/_scaler/src/scaler.c:74-82
The contiguous path reinterprets complex data as 2N reals and uses the real multiply loop; the
non-contiguous fallback promotes the factor to complex. For inf/nan these disagree:
c = np.array([complex(1, np.inf), complex(np.inf, 2)])
numpy c * 10.0 : [nan+infj inf+nanj]
scaler contiguous : [10.+infj inf+20.j]
scaler non-contiguous : [nan+infj inf+nanj]The real-loop answer is arguably the better one, but a single function returning two different
results for the same values based on strides is not defensible. Pick one and make both paths agree.
5. Scaler(10.0, **{}) raises SystemError
astropy/units/_scaler/src/scaler.c:254-259
if (nargs != 1 || kwds != NULL) {
char *kwlist[] = {"", NULL};
PyArg_ParseTupleAndKeywords(args, kwds, "d:Scaler", kwlist, &factor);
return NULL;
}A non-NULL but empty kwargs dict is treated as an error. The parser is then called with valid
input, succeeds, and the function returns NULL with no exception set:
SystemError: <class 'astropy.units._scaler.Scaler'> returned NULL without setting an exception
Fix: guard on kwds != NULL && PyDict_Size(kwds) > 0.
6. Multiplication order flipped: value * scale instead of scale * value
Old _get_converter returned lambda val: scale * _condition_arg(val). Every new path does
PyNumber_Multiply(obj, factor). For any duck-array reaching unit.to(other, value) this now
invokes __mul__ where it used to invoke __rmul__:
scaler(obj) -> left-mul # old: 10.0 * obj -> right-mul
Commutative for ndarray / dask / masked arrays, but it is an unannounced API-surface change and
deserves a changelog line at minimum.
7. make_converter changes equivalency conversion results
astropy/units/core.py:484
- return func(_condition_arg(v) / scale1) * scale2
+ return scaler2(func(scaler1(v))) # scaler1 multiplies by 1/scale1v / scale1 and v * (1/scale1) are not the same in floating point. Over
np.linspace(1, 2, 10001):
scale1 |
differing results |
|---|---|
3.0 |
3342 / 10001 |
7.0 |
3573 / 10001 |
299792458.0 (c) |
1361 / 10001 |
1.602176634e-19 |
1159 / 10001 |
86400 * 365.25 |
6252 / 10001 |
All 1-ulp, but astropy has ulp-level tests in coordinates/spectral. Separately, 1/scale1
overflows to inf for subnormal scale1 (1/1e-320 == inf) where v / scale1 was fine.
8. Stray fence breaks the docs build
docs/changes/units/20183.perf.rst:6
The file ends with a literal ``` line. Sphinx emits
"Inline literal start-string without end-string", and astropy's docs build runs with -W.
Design / lower severity
9. The limited-API guard is inverted
astropy/units/_scaler/src/scaler_limited_api_workarounds.h:6-14
#if defined(Py_LIMITED_API) && Py_LIMITED_API + 0 >= 0x030D0000
/* full flags, incl. Py_TPFLAGS_BASETYPE */
#else
#define PyType_GetModuleByDef(type, unused) PyType_GetModule(type)
/* flags without BASETYPE */
#endifsetup_package.py never defines Py_LIMITED_API and never passes py_limited_api=True, so
astropy's actual build takes the #else branch — replacing the perfectly available
PyType_GetModuleByDef with PyType_GetModule and dropping Py_TPFLAGS_BASETYPE. Confirmed on a
non-limited build:
class Sub(Scaler): pass
TypeError: type 'astropy.units._scaler.Scaler' is not an acceptable base type
Wants #if !defined(Py_LIMITED_API) || Py_LIMITED_API >= 0x030D0000.
Also: if the entire limited-API workaround file is dead code in-tree, that is worth stating
explicitly — the PR description's "w/o limited API" benchmark numbers imply it was meant to be on.
Dormant behind the same flag: Scaler_new ignores its type argument when factor == 1.0 and
returns the base-class singleton, so Sub(1.0) would return a plain Scaler if BASETYPE is
ever enabled.
10. Free-threading
No Py_mod_gil slot, and O_factor / A_factor / A_factor_f are lazily initialised with no
synchronisation (scaler.c:100-110, scaler.c:154-162) — a benign-ish leak race, but there is a
weekly 3.13t job running with PYTHON_GIL=0 that will at least start warning.
11. PyArray_ISALIGNED is not checked
astropy/units/_scaler/src/scaler.c:149 hands PyArray_DATA straight to the ufunc inner loop
after checking only ISONESEGMENT and ISNOTSWAPPED. The ufunc machinery normally buffers
unaligned arrays before calling the loop; this bypasses that. Produced correct results on arm64 in
testing, but it is a real hazard on strict-alignment targets.
12. Undocumented (probably good) behaviour change
Scaler(10.)([True, False]) now returns [10., 0.] where _condition_arg raised ValueError.
Bool arrays already worked, so this removes an old inconsistency — it just needs a changelog line.
13. Nits
Py_tp_docslot discardsconst— compiler warning on every build.PyModule_GetStateresult unchecked atscaler.c:426.Scaler_clear's comment says "called from finalize and dealloc", but there is no
Py_tp_finalizeslot.#endif // _SCALER_PY311_WORKAROUNDS_Hdoes not match the guard name
_SCALER_LIMITED_API_WORKAROUNDS_H.
Summary
| # | Finding | Severity |
|---|---|---|
| 1 | Use-after-free in Scaler_dealloc (segfault) |
Blocker |
| 2 | NPY_TARGET_VERSION bump drops numpy 1.25/1.26 | Blocker |
| 3 | float16 arrays upcast to float64 |
High |
| 4 | Complex inf/nan differs by memory layout | High |
| 5 | Scaler(x, **{}) → SystemError |
Medium |
| 6 | __mul__ instead of __rmul__ |
Medium |
| 7 | Equivalency conversions change by 1 ulp | Medium |
| 8 | Stray ``` fails the -W docs build |
Medium |
| 9 | Limited-API guard inverted | Low |
| 10 | No Py_mod_gil; unsynchronised lazy caches |
Low |
| 11 | PyArray_ISALIGNED unchecked |
Low |
| 12 | Bool lists silently accepted now | Low |
| 13 | Nits | Trivial |
|
What I find useful with reports from Claude like this is to ask it to actually generate some regression tests for these issues - if you provide those it can reduce effort on the part of the PR author (obviously not all things can be written as a test, but it is less work to just read a concise well written regression test and decide how to fix it) |
|
@nstarman - thanks, that LLM report is both interesting and useful. Fixed (EDIT: moved 2 to not followed)
Not followed:
Not completely sure:
@astrofrog - as you have a bit more experience with C, would you be able to review? (note that it is the follow-up that is truly worth it, doubling speed for large arrays...). |
be73429 to
de4a7ab
Compare
I generally agree, but didn't want to come back to this PR like here's 500 lines fixing everything and a full regression test suite. I could if requested. |
|
@nstarman - I think this was just right, as I definitely wanted to think about the problems. Thanks! |
With the inspectable C-based scaler class, which has optimizations for scalars and contiguous arrays, conversions become quite a bit faster ``` import astropy.units as u import numpy as np c = u.km.get_converter(u.m) %timeit c(10.0) 45.7 ns -> 14.1 ns b = np.arange(10.0) %timeit c(b) 392 ns -> 84 ns ``` (With the static type one gets an ~10 ns extra speed-up, but that does not seem worth it in the scheme of things.) Note that while this propagates to `Quantity` operations, the effect is only at the 10% level since there are other quite large overheads: ``` q1 = 1 * u.m q2 = np.arange(10.0) << u.km %timeit q1 + q2 4.19 μs -> 3.66 μs q3 = np.arange(30.0)*u.deg %timeit np.sin(q3) 4.12 μs -> 3.55 μs ``` Also: time needed for construction of the scaler does not change (for that, would need to do caching): ``` %timeit u.km.get_converter(u.m) 704 ns -> 704 ns ```
In particular, let astropy decide on Py_LIMITED_API and NPY_TARGET_VERSION.
de4a7ab to
f11e8eb
Compare
|
I rebased now the caching is merged, and edited the top comment to give combined benefit as well: ~30% reduction in execution time for conversion of small quantity arrays. I'll rebase again once python 3.12 is the minimum version (see #20223), as I then can remove a lot of the work-arounds. |
5 participants
This PR introduces a C-based scaler class, which has two advantages over the lambda that is currently used in unit conversions: it is inspectable, and it is faster, with in particular nice optimizations for scalars and contiguous arrays:
(With the static type one gets an ~10 ns extra speed-up, but that does not seem worth it in the scheme of things.)
Note that while this propagates to
Quantityoperations, the effect is only at the 10% level since there are other quite large overheads:Also: time needed for construction of the scaler does not change (for that, would need to do caching):
Finally, @neutrinoceros - I developed the class separately from astropy, and kept its source and the tests separate here, to help the eventual move of C based code outside of astropy.