Repository navigation
Vectorize string output for Angle.to_string - #20130
veyron-kairo wants to merge 1 commit 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.
|
69ce3be to
cfb70b6
Compare
|
@taldcroft @mhvk — this is the coordinates follow-up to the It's the same idea as before: the sexagesimal math in Roughly ~10x faster for whenever you have a chance to look please review this once . Thanks! |
This comment was marked as resolved.
This comment was marked as resolved.
This comment was marked as resolved.
This comment was marked as resolved.
|
Thanks @pllim — fair on both counts. On the number of tests: to_string wasn't under-tested before. These are new because the PR adds a second, vectorized path alongside the existing per-element one, and the goal is byte-identical output — so most of them just pin "fast path == old path" across the options (separators, padding, fields, precision, sign, latex/unicode, plus NaN/inf and the fallbacks). But you're right that it's more than needed: I ended up checking that same equivalence twice — once through to_string and again directly against the private helper — and over-parametrized the direct one. I've consolidated it down to one equivalence test through the public API plus a small direct/fallback test. On the oldest-deps failures: those are all coming from that over-parametrized test, and it's a mistake in the test rather than an incompatibility. The fast path deliberately falls back to the per-element code on numpy < 2.1 (np.strings.zfill is buggy there), so the helper returns None — and the test was comparing against that None instead of skipping. Trimming it and adding a numpy < 2.1 skip clears all 96. I'll push that shortly. |
|
Thanks @pllim, you're right that there was redundancy. I was checking the same fast-vs-per-element equivalence both through |
|
devdeps failures are unrelated, and benchmark decrease might also be unrelated |
There was a problem hiding this comment.
@veyron-kairo - this looks great to me. It would be best to get @mhvk to have a look since I'm not a coordinates maintainer. Although I looked through all the code and it looked about right, I didn't do a deep line-by-line. The testing looks solid to me though.
I pushed one commit to re-organize the code a bit to put all the sexagesmal processing (the fast-path or the per-element path setup) within the same initial if-block. In the normal fast-path case then this avoid running unnecessary per-element code. This change is pretty straightforward but if you have concerns we can discuss and back it out.
|
Sorry for the slow reply, I was away for a few days. Thanks for the review, and for pushing that commit. I pulled it down and looked it over, and I have no concerns at all. Folding the fast path and the per-element setup into the one sexagesimal if-block reads a lot better than the split version I had, and it also means the fast path stops building the I re-ran the coordinates tests locally on top of your commit and they all pass. The only red job is the allowed-failure devdeps one that @pllim already flagged as unrelated. @mhvk, no rush, but this is ready whenever you have time for a look. |
mhvk
left a comment
mhvk
left a comment
There was a problem hiding this comment.
I was reviewing this in some detail and lhaving mostly comments about readability (replace np.where usage, use a + b instead of np.strings.add(a, b)when I started wondering whether there really isn't more overlap with how it is done inastropy.time. There, I really liked how your code used the existing format strings. I think it would be much nicer to do that here too (effectively using septo create a fixed format string), and perhaps factor out the common code to something inastropy/utils`.
But looking at Time, I was reminded also that there we use an erfa routine to calculate parts. In contrast, here we do things ourselves, which is silly and risky--indeed, as you'll see in one of the in-line comments, I now noticed we round seconds incorrectly to 60!!
A quick look at pyerfa shows there is in fact no need to do things ourselves; one can just convert the angle to radian and use their routines (which are ufuncs wrapped around C code, so pretty fast):
- a2af(ndp, angle): Decompose radians Alternativelyinto degrees, arcminutes, arcseconds, fraction.
- a2tf(ndp, angle): Decompose radians into hours, minutes, seconds, fraction.
Alternatively, one could also divide either the degrees or hours by 24 and use
- d2tf(ndp, days): Decompose days to hours, minutes, seconds, fraction.
With those, I think a refactoring out common code with Time makes even more sense. The only downside of those routines is that they use int32 for the degrees/hours, so the maximum is smaller, 2**31.
But that is still unlikely enough that falling back to the slow route is fine.
Which brings me to a second thought: is there really a reason to keep the slow route at all (once NUMPY 2.0 is dropped...)? Speed is not so important for small arrays, and for the large numbers it is also not obvious: there are really only 15 significant digits in a float, and if you have minutes and seconds, you already used about 4 with precision=0, so really giving more than 9 digits for degrees or hours is pretty pointless. I'd personally be quite happy with deciding that the maximum number of digits for hour/deg is 8 and replacing anything that doesn't fit with an exponential representation.
| ) | ||
|
|
||
| # Carry rounding upwards through the fields, as in the per-element version. | ||
| rounding_thresh = 60.0 - (10.0 ** -(8 if precision is None else precision)) |
There was a problem hiding this comment.
Yikes, this is not your problem, but is definitely wrong:
In [8]: Angle("1d2m59.49999s").to_string(precision=0)
Out[8]: np.str_('1d03m00s')
There is a divide by 2 missing...
There was a problem hiding this comment.
Confirmed, thanks for spotting it. It is not rare either: over 40,000 random angles it hits 322 at precision=0 and 30 at precision=1. Adding the factor of 2 fixes it, and I checked it against correct rounding over a 180,000 value sweep.
I have left it out of this branch for the moment, since it is the scoping question in my comment above and I did not want to pre-empt your answer. The fix, a regression test and a changelog entry are ready to push wherever you would like them.
There was a problem hiding this comment.
Update on this one: it went into #20227, which is now merged, so the fix and its regression test are on main rather than in this branch.
| else: | ||
| d = np.where(m >= 30.0, d + 1.0, d) | ||
|
|
||
| # Degrees have a variable number of digits, so size the buffer from the |
There was a problem hiding this comment.
I think this should be done before the degrees are even calculated, on angle itself.
There was a problem hiding this comment.
Done. Taking the mask from angle also let the second np.where on the minutes further down go away.
One subtlety worth noting: the sign has to be taken before the non-finite values are zeroed, otherwise -inf loses it. The existing tests caught that.
| rounding_thresh = 60.0 - (10.0 ** -(8 if precision is None else precision)) | ||
| if fields == 3: | ||
| carry = s >= rounding_thresh | ||
| s = np.where(carry, 0.0, s) |
There was a problem hiding this comment.
I think this reads better as,
m[carry] += 1
s[carry] = 0
Since it is in-place, this will be faster for large arrays (say, >1e4 or 1e5).
Same below, also for the string operations.
There was a problem hiding this comment.
Done, and applied to the string concatenation below as well.
For what it is worth, the commit as a whole is worth about 5% on a million elements (424 ms to 402 ms) and nothing measurable at 10k, so it is mostly a readability win at the sizes people actually use.
|
Thanks @mhvk, this is a really useful review. On the seconds rounding: confirmed, and it's not specific to the fast path.
On erfa: I hadn't looked at Before I start moving things around, two questions on scope:
On dropping the slow path once numpy 2.0 goes: I like the idea, and the 8 digit cap with exponential beyond it sounds reasonable to me. I'd rather do that as a separate step after the erfa work, so the behaviour change isn't mixed in with the refactor. I'll do the in-place updates and move the finiteness check onto |
|
I did a bit of measuring while waiting, in case any of it is useful. The rounding problem turns out to be more common than I expected. Over 40,000 random angles at Your point about erfa being the safer foundation also holds up in a way I hadn't expected. For 67.94486111111111 degrees, our own decomposition gives 41.49999999999295 seconds where The one thing I wanted to raise before starting is that moving to erfa isn't output preserving, over and above fixing the bug. erfa breaks ties away from zero while Python's
So something entered as One small practical note on the erfa routines: the fraction field is None of this changes my mind that erfa is the right way to go. It mainly affects how the work is best split up and what the changelog needs to say. |
|
Pushed the in-place rounding, the I've left the seconds rounding alone for now, since that's one of the questions above and I didn't want to pre-empt your answer on whether it belongs here or in its own PR. It's ready to go either way. |
|
@veyron-kairo - on separate PRs, I think I agree that a small PR just fixing the rounding error is better, as it allows us also to put at test in place for it. It will also make review of this PR easier. On the other hand, I don't want to make extra work, so I'm OK to go with what you decide. |
|
p.s. On |
|
Quick update: I have the erfa version working locally and will push it in a bit, once I have tidied up the commits. First, a correction to my numbers above. The 22% I quoted was measured with So the real change is smaller than I made it sound. Comparing the branch against main over 1.2 million formatted values (both units, precisions None/0/1/3/6, all three field counts), 0.911% differ. The worst case is Something else fell out of that comparison which I do not think has come up yet: The angle is 285d29m41.76s, so it should be 285d. The current code carries the seconds into the minutes at 30 first, which takes the minutes to exactly 30, and then rounds the minutes into the degrees at 30 as well, so it rounds twice. It affects 295 of 40,000 random angles. erfa gets it right because it is handed the whole angle and a single resolution. Would you like that in the same small bugfix PR as the seconds rounding? It is the same kind of mistake and it seems odd to fix one and leave the other. Last thing, on scope. To use a format string here the way For what it is worth, on a million angles |
|
@veyron-kairo - wow, that's very nice! Your plan sounds good: let's do the fix first (both seconds and minutes), and then the larger one to vectorize. I see no problem at all with touching |
|
Thanks, that all makes sense. I'll put the fix up first with both the seconds and the minutes, then rebase this one on top of it. On where the shared code goes, I ended up giving it a small private module of its own rather than |
|
My sense is to put it in |
|
Actually, let me retract that: why don't we first just extract more of the code in the right function in Also, yet another thing, I'm not sure why in the slow path we do not use python's round instead of that whole |
|
Good call on Both fixes are in, PR coming shortly. And noted on keeping the extraction inside |
20b3084 to
1e579a1
Compare
|
Pushed the erfa version. I kept the extraction inside I also rebased on main, so the seconds and minutes rounding fix from #20227 now sits underneath this. That covers the On a million angles |
mhvk
left a comment
mhvk
left a comment
There was a problem hiding this comment.
apologies for not responding earlier, but with term rapidly approaching, I have only limited time for review. As you'll see inline, my main request now is to keep things simpler.
| ) | ||
| if result is not None: | ||
| if alwayssign: | ||
| result = np.where( |
There was a problem hiding this comment.
I would still prefer to use direct operators over np.strings.add...
There was a problem hiding this comment.
Done, the latex delimiters are "$" + result + "$" now. The alwayssign one has gone entirely, since the sign is added where the fields come from, as you suggest below.
| ) | ||
| if format == "latex": | ||
| result = np.strings.add("$", np.strings.add(result, "$")) | ||
| is_nan = np.isnan(values) |
There was a problem hiding this comment.
Why not just:
results[np.isnan(values)] = "nan"
There was a problem hiding this comment.
Done. It fits in place because the non-finite values were written as inf first, so the dtype is always wide enough for nan.
| result = np.where(is_nan, "nan", result) | ||
| # NumPy string ops can collapse a scalar to a 0-d string; make | ||
| # sure we return a bare Python string in that case, as before. | ||
| result = np.asarray(result) |
There was a problem hiding this comment.
It is not. The fast path only runs when there are more than _VECTORIZE_MIN_SIZE elements, so the result is always an array. The asarray and the ndim check on the return are both gone.
| if angle.size and np.abs(angle).max() >= _ERFA_MAX_FIELD: | ||
| return None | ||
|
|
||
| _, parts = erfa.d2tf(_COARSE_RESOLUTION.get(fields, ndp), angle / 24.0) |
There was a problem hiding this comment.
I think we might as well keep the sign here, and use it below, i.e., do sign, parts = ...
I also think we should pass alwayssign down to here since then we can just deal with adding the sign once, using
if not alwayssign:
sign[sign == "+"] = ""
There was a problem hiding this comment.
Done, it is sign, parts = erfa.d2tf(...) now and alwayssign is passed down, so the sign goes on in one place and the extra pass in core.py is gone.
The one thing that still needs the sign bit is negative zero: ERFA reports it as +, and to_string has always written it as -0, so sign[negative] = "-" stays.
| # its own column and stripped, rather than going into the template. | ||
| degrees = parts["h"] | ||
| width = len(str(int(degrees.max()))) if degrees.size else 1 | ||
| out = _format_template(f"{{deg:0{width}d}}", {"deg": degrees}, degrees.shape) |
There was a problem hiding this comment.
I'm very confused why there are two calls to _format_template rather than just one that uses minutes and seconds as well (as appropriate for fields).
There was a problem hiding this comment.
The leading field is the only one with a variable number of digits. _format_template writes fixed width columns, so the degrees come out as 005 and the zeros have to come off afterwards, and stripping them off the joined string is not safe: with an empty separator and zero degrees, the 000 runs straight into the 05 of the minutes and the strip eats that zero too. So the leading field gets its own call, and the minutes and seconds, which keep their zeros, go in the second one.
There is a way to do it in one call: build the whole template, then take the leading zeros off with np.strings.slice and a per element start worked out from the digit count of the degrees. I tried it and the output is identical, it costs about 6% (143 ms now against 151 ms for a million elements), but np.strings.slice is NumPy 2.3 while the rest of the fast path only needs 2.1, so on 2.0 to 2.2 it would mean no speedup at all. That is the only reason I did not go with it. Happy to switch if you would rather have the single call.
| # large for its 32-bit fields, and its integer fields have no room for the | ||
| # non-finite values, all of which are left to the slower code below. | ||
| if 0 <= ndp <= 9 and abs(angle) < _ERFA_MAX_FIELD: | ||
| _, parts = erfa.d2tf(_COARSE_RESOLUTION.get(fields, ndp), angle / 24.0) |
There was a problem hiding this comment.
I think performance of conversion to string is not worth duplicating code unless it necessary, so I'd prefer to just calling the string version. With good use of operators instead of np.strings functions, it might in fact work almost as well anyway.
There was a problem hiding this comment.
I tried it, and the fixed cost of going through the array version for one value is high: the helper goes from 3.1 to 31 us, Angle.to_string on a scalar from 7.6 to 36 us, and an eight element array (just under the vectorization threshold) from 34 to 255 us. That is slower than main as well, where the helper is 3.8 us. It is not the np.strings calls, it is the half dozen array operations and the two buffers, so operators do not win it back.
What is actually duplicated is the guard and the erfa.d2tf call, four lines or so. What follows differs anyway, since one builds a Python str and the other writes into a code point buffer.
Going the other way and taking ERFA back out of the scalar path is not an option, since then the two disagree. Over 30,000 random angles they differ in the last digit for about 0.4% of them at the default precision, and on ties at precision=0 and 1, and it is ERFA that is right there: -247.51375 deg is exactly -247d30m49.5s, which ERFA rounds to 50 and the float decomposition gives as 49.
So I would rather keep it, but it is your call, and if you still prefer the one implementation I will make the change.
|
sorry for the delay @mhvk , got busy in internals and exams , didn't got enough time to resolve this ! |
|
Sorry for the slow reply on this one. I have pushed the changes for the straightforward points, and left notes on the two about the structure of the helper, where I would rather have your call than guess. |
b5f6bc1 to
0d7f7a4
Compare
|
@veyron-kairo - just to let you know my teaching started this week and that means I'll be slow to review; may have some time after next week. |
|
Thanks for letting me know @mhvk , please take you time 😊 |
|
@mhvk - just a ping, might you have time to review this? |
|
Thanks for the ping. I fear that my tardiness in replying is an answer in itself -- my teaching leaves little time this fall, and I don't think I'll get to this soon, sadly. |
|
@veyron-kairo - can you add an AI disclosure to the description to comply with current policy? It helps the review to know what model(s) you used to help out with this PR. |
|
@taldcroft , did updated the pr description with Ai disclosure |
|
This is not independent, but I asked Claude Fable 5.1 for a review and there are some issues that should be addressed or at least discussed. Thanks for this, the speedup for large arrays is real and the approach of building a template once and filling it with array ops is sound. I checked out the branch, ran the coordinates tests (734 passed), and compared the output against 1. Rounding of the last digit changed (blocking)The new scalar path (and therefore the array path) goes through Comparing against
The second is an exact value (187d28.4999999999997m) that This changes I think the vectorized path needs to reproduce the old decomposition with numpy ( 2. Masked Angles lose their mask on the fast path
a = Angle(Masked(np.linspace(0, 100, 12), mask=[True, False] * 6), 'deg')
a.to_string()
# main: MaskedNDArray([———, '9d05m27.27272727s', ———, ...])
# PR: array(['0d00m00s', '9d05m27.27272727s', '18d10m54.54545455s', ...])
a[:3].to_string()
# PR: MaskedNDArray([———, '9d05m27.27272727s', ———]) (loop path, still correct)So the same array formats differently depending on its length. 3. Overflow guard is applied before roundingThe Angle(2**31 - 0.5, 'deg').to_string(fields=1)
# RuntimeWarning: invalid value encountered in d2tf (an error under astropy's filterwarnings)
# returns '2147483647' on arm64 (saturating cast); on x86 the cast of 2^31 is undefined
# main returns '2147483648'Same on the array path. Guard with a margin, or check the returned 4. A
|
|
@taldcroft , thanks for pointing out the issues I'll fix it ASAP |
0d7f7a4 to
292f7c0
Compare
|
@taldcroft all nine were real, all fixed now. First a correction to what I told @mhvk earlier: -247.51375 is not exactly -247d30m49.5s, it is -247d30m49.49999999995s, so main's 49 was right and ERFA's 50 was not. I had compared against the decimal literal instead of the float. With Fraction on the binary value, d2tf gets the last digit wrong on 267 of 50,000 angles at the default precision, the decomposition on none. So ERFA is out. The array path is back on the modf decomposition, printf does the seconds, and the carry comes from the 60 it can produce. Output matches main byte for byte over 960 keyword combinations and 12,005 values. A million angles take 0.39 s against 4.4 s on main, SkyCoord.to_string 0.83 s against 23 s. The rest are in as well: masked angles keep their mask (SkyCoord stays on the per-element path, the mask has to sit inside the joined string), the int64 boundary is tested now, a sep with braces falls back to the loop, NaN is written in the helper, the size check uses >=, and the new rounding test pins your three values plus two more and fails on the old branch. The cross-package import I left alone, @mhvk wanted the extraction kept inside astropy.time for now. |
|
@veyron-kairo - with apologies, I'm just copying another Claude review. I honestly don't have time to give this the attention it needs, nor am I a coordinates maintainer. But hopefully this helps in getting the PR to a good place where @mhvk can give a real review. The Claude review is basically good now with only some small cleanups that are not strictly required.
@veyron-kairo - thanks, this is a big improvement. I reviewed Earlier review items
What I ranI compared this branch against its merge base, using a separate script from the tests in this PR:
All valid numeric output was identical, and warnings matched in a separate run on NumPy 2.2.3. The only differences were items 1 and 2 below, plus an empty I also confirmed your correction on ERFA. Checked with Speed, compared to
The array path is first faster at 8 elements (48 µs against 67 µs). Small items
Before merging
Disclosure: this review was prepared with assistance from Claude Code (Opus 5.5). All results above were run against the branch and its merge base with |
Angle.to_string formatted every element through np.vectorize, so large arrays paid the full Python overhead per value. Build the strings with array operations instead, reusing the fixed-width string builder from the Time formatting, which is moved into module-level helpers in astropy.time.formats for this. The fields are decomposed and rounded exactly as in the per-element formatter, so the output is identical. Small arrays, NumPy < 2.1, degrees of 1e18 and up, and separators the template cannot take still go through the loop. SkyCoord.to_string formats both angles for the whole array and joins them, except for masked coordinates. Co-authored-by: Tom Aldcroft <[email protected]> Signed-off-by: Shridhar Panigrahi <[email protected]>
2ae5f0b to
803e5d5
Compare
|
@taldcroft thanks a lot for going through it again, really appreciate it.
I've squashed everything into one commit and kept you as co-author for your part. Will update the PR description next. @mhvk whenever you get time, the two things left are dropping ERFA and whether the scalar path should call the array one. |
Fixes #19992.
Angle.to_stringformatted every element in a Python loop (np.vectorize), so large arrays were slow. This builds the sexagesimal strings with array operations instead, the same way #19942 did forTime.TimeString._value_fastis moved into_format_templateinastropy.time.formats(kept private there, per @mhvk), and coordinates uses it.1e18and up, and asepthat is not a string or contains a{name:spec}field still use the loop. Masked angles keep their mask.SkyCoord.to_stringformats both angles for the whole array and joins them, except for masked coordinates.For 1,000,000 elements on NumPy 2.4.6:
Angle.to_string(sep="dms")SkyCoord.to_string("hmsdms")Scalars are unchanged. Tests compare the new path against the per-element one, including NaN, inf, masked input and every fallback. Changelog and What's New entry included.
AI Disclosure
I used Claude (Opus) to help write the tests in this PR. I went through all of it myself, ran the tests and the benchmarks locally, and I'm responsible for every change here.