TST: a failing doctest on MacOS + numpy 2

[neutrinoceros]

There were a number of intentional changes for MacOS arm64 in numpy 2.0, which was released 3 days ago. It seems that as a side effect of the performance optimizations conducted a table printed in docs/fitting.rst changed; # doctest: +FLOAT_CMP helpfully allows to ignore discrepancies in lower digits, but it doesn't protect against what's happening here: the number of digits printed changed. See example logs.

I tried to make this doctest more portable using the +ELLIPSIS directive and ended up with the following patch

diff --git a/docs/fitting.rst b/docs/fitting.rst
index b4a3786..4db9e35 100644
--- a/docs/fitting.rst
+++ b/docs/fitting.rst
@@ -147,16 +147,13 @@ For example, based on the spectrum defined above we can first select a region:
 
 Then estimate the line  parameters it it for a Gaussian line profile::
 
-   >>> print(estimate_line_parameters(sub_spectrum, models.Gaussian1D()))  # doctest: +FLOAT_CMP
+   >>> print(estimate_line_parameters(sub_spectrum, models.Gaussian1D()))  # doctest: +FLOAT_CMP, +ELLIPSIS
       Model: Gaussian1D
       Inputs: ('x',)
       Outputs: ('y',)
       Model set size: 1
       Parameters:
-              amplitude            mean             stddev
-                  Jy                um                um
-          ------------------ ---------------- -------------------
-          1.1845669151078486 4.57517271067525 0.19373372929165977
+      ...
 
 
 If an `~astropy.modeling.Model` is used that does not have the predefined

I could confirm it fixes the test on macos arm64, but it's inconvenient that I just ended up hiding the whole table; unfortunately my attempts to only hide parts of it did not succeed. Would this patch be accepted as is ?

[dhomeier]

I suspect that is an unfortunate coincidence in the original numbers, that mean happens to be something like 4.57517271067524999 which is then rounded to the 14th rather than the 16th significant digit. If you replace the number in the docstring with 4.5751727106752499, will that in turn break the older architectures? In that case might try to insert np.set_printoptions(floatmode='fixed') at the top...

[neutrinoceros]

If you replace the number in the docstring with 4.5751727106752499, will that in turn break the older architectures?

if we do just that, yes, it will break, because the number of - on the previous line still depends on the number of digits printed and # doctest: +FLOAT_CMP doesn't account for that.

In that case might try to insert np.set_printoptions(floatmode='fixed') at the top...

I tried it, but it doesn't seem to have the desired effect (I'm guessing the embedded table might not be composed of numpy arrays at all). In case we circle back to this idea later however, note that we should use the context manager np.printoptions instead to avoid side effects in other doctests.

[dhomeier]

In that case might try to insert np.set_printoptions(floatmode='fixed') at the top...

I tried it, but it doesn't seem to have the desired effect (I'm guessing the embedded table might not be composed of numpy arrays at all).

Unfortunate, but I also noticed it seemed to have no effect on numpy scalars like np.float64(4.5751727106752499) – don't know what exactly the table column is using for __str__ or __repr__.

In case we circle back to this idea later however, note that we should use the context manager np.printoptions instead to avoid side effects in other doctests.

I thought of that as well, but isn't every section running in its own environment anyway – at least each code-block that starts a new import numpy as np etc.? But it is reusing spectrum and other things from the previous blocks, so it's not that fine-grained. Still there are some other cases of column values printing out at different lengths above.

Last way out I can think of would be to change the spectrum parameters slightly so they don't round off to fewer digits, but that could end up in wild guessing. Still, hiding the actual table values kind of beats the whole point of the documentation in my view; I'd rather +SKIP the entire print (will it test any values with +ELLIPSIS anyway?)

[neutrinoceros]

I thought of that as well, but isn't every section running in its own environment anyway – at least each code-block that starts a new import numpy as np etc.?

Oh but I'm not guessing here: I actually tried your suggestion and found that it broke 3 other tests even with numpy 1.x I should maybe note that importing an already imported package doesn't clean up state, so it's not clear to me how doctests are separated into sessions (if at all).

[dhomeier]

Oh, and I thought you found no effect with printoptions at all – at least it is working in principle, then ;-) You are probably right, as even running sub-package tests on their own vs. the entire suite can result in very weird behaviour (and I have given up to trust things like importlib.reload sometime back in Python 2.7)...

[neutrinoceros]

So I see two ways this can go down:

1. slap a cheap workaround for it (like the one I proposed originally, or skip doctests on macOS...)
2. go the long way and figure out why the printed table isn't made of numpy objects so we can discuss if and how we can make it so (hence, enabling np.printoptions)

I wouldn't even try to go with 2) unless you greenlight it, as I'm really unsure how long this would take and if it would even work out in the end.

[pllim]

Looking here

https://github.com/astropy/specutils/blob/d16d91a077417a7aacdee06d2b259fa887908b14/specutils/fitting/fitmodels.py#L61

I think it is basically an existing astropy.modeling behavior. How did we patch the doctest for astropy.modeling? The same trick would probably work here.

👎 on blind ellipsis because it makes the example very opaque.

[dhomeier]

It is possible modeling docs simply did not have any output that happens to be rounded, at least could not find any changes along with the numpy 2 migration there. Not sure if astropy/astropy#15065 is related?

[neutrinoceros]

It is possible modeling docs simply did not have any output that happens to be rounded, at least could not find any changes along with the numpy 2 migration there

I confirm that we didn't need to do anything for modeling and I concur to your conclusion.

Not sure if https://github.com/astropy/astropy/pull/15065 is related?

It is not.

[pllim]

Looking at the log again

    -   Model: Gaussian1D
    -   Inputs: ('x',)
    -   Outputs: ('y',)
    -   Model set size: 1
    -   Parameters:
    -           amplitude            mean             stddev
    -               Jy                um                um
    -       ------------------ ---------------- -------------------
    -       1.1845669151078486 4.57517271067525 0.19373372929165977
    +Model: Gaussian1D
    +Inputs: ('x',)
    +Outputs: ('y',)
    +Model set size: 1
    +Parameters:
    +        amplitude             mean               stddev      
    +            Jy                 um                  um        
    +    ------------------ ------------------ -------------------
    +    1.1845669151078486 4.5751727106752496 0.19373372929165975

Would ellipsis work partially like this?

1.1845669151078486 4.5751727106752... 0.19373372929165...

[neutrinoceros]

I'm not saying it doesn't, but I couldn't figure it out after a 30min attempt 🤷🏻‍♂️ The real issue isn't with floating numbers though (those are handled with +FLOAT_CMP), it's with the rest of the table having a varying line length.

[pllim]

The varying line length is indeed surprising... and what does that have to do with numpy rounding? 🤯

[pllim]

estimate_line_parameters(sub_spectrum, models.Gaussian1D()) basically returns a model object. So might be possible to change the narrative to not print out the whole model, but rather specific things that the example is trying to illustrate in a more controlled manner. Maybe @rosteen can advise.

[neutrinoceros]

The varying line length is indeed surprising...

Well it's caused by a difference in rounding: we get a different number of decimals depending on the platform and in turn, the Table object's layout also varies.

what does that have to do with numpy rounding?

This is just my current conjecture but what I think is happening is that numpy 2.0 made some performance improvements specifically for macOS arm64 (using the Accelerate Apple lib), which sometimes result in slight changes in numerical results (typically 1 ULP). Ironically, the numbers contained in this Table appear to not be numpy arrays/scalars (I assume they're just plain Python floats), so we cannot easily control the precision with which they are printed.

[neutrinoceros]

So might be possible to change the narrative to not print out the whole model,

Note that the last sentence before this doctest contains at least one typo (it it), but I don't grok what the intended sentence would have been, but this might be an opportunity to address this too.

[dhomeier]

Ironically, the numbers contained in this Table appear to not be numpy arrays/scalars (I assume they're just plain Python floats), so we cannot easily control the precision with which they are printed.

As I reproduce how this is created in https://github.com/astropy/astropy/blob/5519605a3b6ec71e9ddb4c6827b8ac1658647f58/astropy/modeling/core.py#L2919-L2929 it should consist of Columns which are np.array-valued, but as Table.pformat is creating its representation per-row, this still ends up using the numpy scalar elements, which are not controllable by np.printoptions. And I don't see any options for controlling this within the table formatter either.

A somewhat hackish way to resolve this in astropy.modeling would be to pad param_table with an extra row of all e.g. 1.111111111111 (the default max. 12 digits after the decimal) and then work with param_table.pformat()[:-1] instead of str(), which should always stretch the output to the full 14 digits. But that would still run into problems for parameters represented in exponential notation...

[pllim]

I still get pinged by Actions about this once a week or so.

[rosteen]

#1159 isn't particularly elegant but avoids the failure by printing the three parameters individually, rather than printing the table with all three at once. Any objections?