Fixed incorrect string representation of floats in Cards #14508

kYwzor · 2023-03-09T11:08:51Z

Description

This pull request is to address a bug that caused io.fits.Card to format floats incorrectly, which could force the associated comments to be truncated.

EDIT:

github-actions · 2023-03-09T11:09:27Z

mhvk

I think we may need upper case. Some other comments that are less directly related to your PR.

mhvk · 2023-03-09T14:14:33Z

astropy/io/fits/card.py

-    if "." not in value_str and "E" not in value_str:
-        value_str += ".0"
-    elif "E" in value_str:
+    value_str = str(value)


I think we may have to do str(value).upper() to ensure we still have an E - that may well be FITS standard.

I've changed the if and .split() to check for "e" instead of "E" since that is how python represents the exponent. I believe the upper case E is indeed part of the FITS standard, however that is not an issue since I did not change the following line:

astropy/astropy/io/fits/card.py

Line 1314 in 60835aa

value_str = f"{significand}E{sign}{int(exponent):02d}"

This will force the end result to always use the upper case E.

Good point. So, then it is more relevant if we actually can remove that whole stanza (which I think is likely to be the case).

Sounds a bit risky but we can give it a try. I notice that the code specifically has int(exponent), but I don't think it's possible that Python ever returns a fractional exponent (1e0.5 isn't even legal syntax). Also, 1e-9 is converted to "1e-09" so it seems Python automatically pads to at least 2 digits. The only way this would fail is if indeed the representation is different in different platforms.

This part looks very sketchy to me:

astropy/astropy/io/fits/card.py

Lines 1315 to 1316 in 60835aa

elif "." not in value_str:

value_str += ".0"

The only way this if statement would be called is if a integer reached this point, but I don't think that could ever happen because we only call this function in these places:

astropy/astropy/io/fits/card.py

Lines 1287 to 1292 in 60835aa

elif isinstance(value, (float, np.floating)):

return f"{_format_float(value):>20}"

elif isinstance(value, (complex, np.complexfloating)):

val_str = f"({_format_float(value.real)}, {_format_float(value.imag)})"

return f"{val_str:>20}"

I've checked and using both Python built-in complex numbers and numpy's complex numbers, .real and .imag are float even if we use 1j+1. In other words, I don't see how we could get a number that does not have a decimal point while also not using an exponent.

Agreed, that can be removed too! It was needed only with the format string:

In [2]: f"{5.0:16G}" Out[2]: ' 5'

mhvk · 2023-03-09T14:16:38Z

astropy/io/fits/card.py

-    elif "E" in value_str:
+    value_str = str(value)
+
+    if "e" in value_str:
        # On some Windows builds of Python (and possibly other platforms?) the


I doubt this is still true.

Do you mean the comments regarding the Windows builds of Python?

Yes, I think python now no longer relies on anything OS-dependent for typesetting float values (especially since we no longer use the format string).

astropy/io/fits/card.py

kYwzor · 2023-03-09T14:57:57Z

Does anyone know why Read the Docs does not like me using "io.fits.Card" nor "Card" on the changelog?

docs/changes/io.fits/14508.bugfix.rst

Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com>

kYwzor · 2023-03-09T17:47:07Z

I tried to be a bit more aggressive with the changes and only truncate if we are using a scientific notation. I assumed that the string representation of floats in Python couldn't go over 20 characters in instances where scientific notation isn't used. However, this is apparently not true (for example, -0.0021428662759603867 and 0.0014104245585713215 both get represented exactly as they are written), so the old code did truncate these. So I will revert that part to the old way of truncating for now, at least while we don't opt for rounding instead. I do wonder though: would it be possible to have a float that is 22 characters long without scientific notation and just a single decimal place? Because if so, with the way the strings are being truncated, we could hypothetically truncate the decimal point (which would break the standard).

mhvk · 2023-03-09T18:56:48Z

Probably best to leave the truncation, rounding for another time. Sounds like standard notation can give 22 characters, negating your second example:

In [3]: len(str(-0.0014104245585713215))
Out[3]: 22

We could only make this 21 by removing the leading 0, so not a real solution.

kYwzor · 2023-03-09T19:38:36Z

Sorry, maybe my comment wasn't clear. What I meant is that from my testing I figured out that floats can surpass 20 chars even when not using scientific notation (which I initially didn't think possible). My question now is: is there a float for which Python gives a string representation that is 1) standard notation (not scientific), 2) more than 20 characters long and 3) has just a single decimal place? Because if such a float exists, the truncation code would break, since it would delete the decimal places, making it an invalid float in FITS.

mhvk · 2023-03-09T21:19:03Z

@kYwzor - I think the truncation itself is safe, just (slightly) wrong. Probably best to keep that for another PR, as you suggested earlier...

mhvk

This looks all good to me - let's handle the truncation separately.

Approving now, but will wait with merging to give @saimn a chance to have a look.

saimn

Changes look good, thanks for the PR @kYwzor (and for the review @mhvk :))

keflavich · 2023-03-13T16:48:34Z

This broke some tests downstream in spectral-cube because we're using wcs_out.wcs.compare(result.wcs.wcs). I'm still trying to ID the issue specifically, but I may request a change.

keflavich · 2023-03-13T17:04:02Z

OK I think I see the problem - the error noted in #14507 is possibly not a bug. It is not possible to round-trip WCS->header->WCS with this PR in place.

However, I can't produce a MWE; I don't know exactly what wcs.compare is doing, and in the case where it's returning False, the WCSes are by-eye identical. I'm still digging.

keflavich · 2023-03-13T17:09:31Z

MWE:

from astropy.wcs import WCS
ww = WCS(naxis=1)
ww.wcs.cdelt[0] = 1288.2149687900003
ww2 = WCS(ww.to_header())
ww.wcs.compare(ww.wcs) # True
ww.wcs.compare(ww2.wcs) # False

keflavich · 2023-03-13T17:18:47Z

Could we please add the above example as a regression test and revert this until we come up with a solution?

keflavich · 2023-03-13T17:21:02Z

To expand a bit: it looks like we have a problem with round-tripping from header->wcs->header identified in #14507, but the solution here breaks the round-trip from wcs->header->wcs.

pllim · 2023-03-13T17:33:13Z

PR to revert at #14524

dhomeier · 2023-03-13T18:22:29Z

(ww.wcs.cdelt[0] - ww2.wcs.cdelt[0]) / ww.wcs.cdelt[0]
1.7650289815899323e-16

is below float64 resolution of 1e-15, so it may not be reasonable to expect this can roundtrip. In fact that example also fails in 5.2, even when setting ww.wcs.cdelt[0] = 1288.214968790003 (i.e. more than resolution above 1288.21496879).
But before and after this PR

fits.Card('CDELT1', ww.wcs.cdelt[0], 'Coordinate increment at reference point')
('CDELT1', 1288.2149687900003, 'Coordinate increment at reference point')
ww.to_header().cards[2]
('CDELT1', 1288.21496879, 'Coordinate increment at reference point')

as the latter is created by WCS.to_header(), with a default precision of {:22.14G} set by WCSLIB (WCSHDO_P14) – changing it to the desired precision here would require

 ww.to_header(relax=wcs.WCSHDO_P17).cards[2]
('CDELT1', 1288.2149687900003, 'Coordinate increment at reference point')

does not seem like it falls within the scope of this PR or the discussion to be had on https://github.com/astropy/astropy/pull/14508/files#r1133599288.

pllim · 2023-03-13T18:26:53Z

So does it mean we need to revert or not?

keflavich · 2023-03-13T18:27:20Z

Ah, I'm afraid you're right, my MWE does not succeed on 5.2 either, it fails on both, so it is not a good MWE. The test that fails in spectral-cube has a little more going on in it; I thought I had narrowed it down to the right spot. I used git bisect to pin the spectral-cube failure down to 3d8a2c6, but I can't tell yet why spectral-cube's test succeeds on 5.2.1 while my MWE still fails. I'll try to produce a better one.

keflavich · 2023-03-13T18:51:36Z

@pllim we still do need to revert, I think - just because we haven't found where this breaks the spectral-cube test does not mean that this is a safe change.

I'm finding it extremely difficult to produce a MWE that passes on 5.2.1 and fails on main, but spectral-cube's test definitely does.

pllim · 2023-03-13T18:53:56Z

Are you sure it is this PR?

keflavich · 2023-03-13T19:00:34Z

Sure as I can be. I ran git bisect and it points to the commit I referenced above. I verified that spectral-cube still passes when built against last night's nightly build.

pllim · 2023-03-13T19:06:50Z

I'd rather see that better MWE first. Or you can also propose a revert PR from your fork.

mhvk · 2023-03-13T19:14:04Z

Hmm, sounds like a wcs.is_close (or a loosening in wcs.compare) is called for... It does seem that @kYwzor's comment is right, that the fact that it did work is more "bug-for-bug compatibility" than anything else.

keflavich · 2023-03-13T19:20:27Z

I've given up on the MWE. I have reproduced the path as precisely as possible and cannot get the MWE to pass with 5.2.1 and fail with main. I'm just changing the comparison to wcs.compare(..., tolerance=1e-12) and moving on.

mhvk · 2023-03-13T19:29:21Z

Ah, so it is like .is_close() already - certainly, sounds like a sensible solution to just allow for some floating point errors.

dhomeier · 2023-03-13T19:32:03Z

Yes; it is correct that that commit has changed the output of

# 5.2.1
fits.card._format_value(1288.2149687900003)
'       1288.21496879'
# 5.3dev
fits.card._format_value(1288.2149687900003)
'  1288.2149687900003'

which happened to just match the precision limitation of wcs.to_header in that specific case (but not even for 1288.21496879003 or 1288.214968790003). So wcs.compare did fail in those examples and still does; seems not a helpful solution to re-introduce the precision loss in fits.card here, but it is a puzzle where card._format_value would have been used in that spectral-cube test so that it could just make up for the WCS roundtrip bug.

dhomeier · 2023-03-13T19:38:30Z

I'm just changing the comparison to wcs.compare(..., tolerance=1e-12) and moving on.

Alternatively you might still try the header_out.update(wcs_out.to_header(relax=wcs.WCSHDO_P17)) option from above; still unclear why the default precision has to be lower (WCSLIB design decision?).
That version does pass with astropy 5.2.1 anyway (which I really was not sure of at this point ;-).

dhomeier · 2023-03-13T21:33:46Z

@keflavich it seems the spectral-cube bug is already introduced when reading the input file, which actually has

CDELT3 = 1.28821496879E+00 /

but is parsed (with any astropy version) as

fits.open('adv.fits')[0].header['CDELT3']
1.28821496879
SpectralCube.read('adv.fits')._header['CDELT3']
1.28821496879
SpectralCube.read('adv.fits').header['CDELT3']
1.2882149687900002

Likely there is some floating point error introduced in BaseSpectralCube.header, perhaps in
https://github.com/radio-astro-tools/spectral-cube/blob/8747c377e190c0f5d193f9dd9916e1fe4b4af1c1/spectral_cube/spectral_cube.py#L2516-L2521
(although in the present case there should be no unit conversion) that is creating the odd values here.

Fixed incorrect string representation of floats in Cards

3d8a2c6

kYwzor requested a review from saimn as a code owner March 9, 2023 11:08

github-actions bot added the io.fits label Mar 9, 2023

kYwzor added 6 commits March 9, 2023 11:10

Fixed incorrect variable name

7566bc2

Added changelog entry

a65f074

Fixed wrong check for exponent symbol

89d4691

Fix test_floating_point_string_representation_card assertion

03730f3

Fixed incorrect changelog

50c51e3

Fixed incorrect changelog reference

60835aa

mhvk reviewed Mar 9, 2023

View reviewed changes

pllim added the Bug label Mar 9, 2023

pllim added this to the v5.3 milestone Mar 9, 2023

pllim reviewed Mar 9, 2023

View reviewed changes

docs/changes/io.fits/14508.bugfix.rst Outdated Show resolved Hide resolved

kYwzor and others added 3 commits March 9, 2023 15:01

Fix Card reference in changelog

b3cc12c

Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com>

Simplify handling of scientific notation in cards

631b607

Fix incorrect assignment

13db404

Revert truncation

7054d58

mhvk approved these changes Mar 9, 2023

View reviewed changes

saimn approved these changes Mar 13, 2023

View reviewed changes

saimn merged commit 386b52e into astropy:main Mar 13, 2023
20 of 21 checks passed

keflavich mentioned this pull request Mar 13, 2023

Expand test coverage radio-astro-tools/spectral-cube#854

Merged

pllim mentioned this pull request Mar 13, 2023

Revert "Fixed incorrect string representation of floats in Cards" #14524

Closed

	elif isinstance(value, (float, np.floating)):
	return f"{_format_float(value):>20}"

	elif isinstance(value, (complex, np.complexfloating)):
	val_str = f"({_format_float(value.real)}, {_format_float(value.imag)})"
	return f"{val_str:>20}"

Fixed incorrect string representation of floats in Cards #14508

Fixed incorrect string representation of floats in Cards #14508

Conversation

kYwzor commented Mar 9, 2023 • edited by pllim

Description

github-actions bot commented Mar 9, 2023

mhvk left a comment

Choose a reason for hiding this comment

mhvk Mar 9, 2023

Choose a reason for hiding this comment

kYwzor Mar 9, 2023 • edited

Choose a reason for hiding this comment

mhvk Mar 9, 2023

Choose a reason for hiding this comment

kYwzor Mar 9, 2023 • edited

Choose a reason for hiding this comment

kYwzor Mar 9, 2023

Choose a reason for hiding this comment

mhvk Mar 9, 2023

Choose a reason for hiding this comment

mhvk Mar 9, 2023

Choose a reason for hiding this comment

kYwzor Mar 9, 2023

Choose a reason for hiding this comment

mhvk Mar 9, 2023

Choose a reason for hiding this comment

kYwzor commented Mar 9, 2023

kYwzor commented Mar 9, 2023

mhvk commented Mar 9, 2023

kYwzor commented Mar 9, 2023

mhvk commented Mar 9, 2023

mhvk left a comment

Choose a reason for hiding this comment

saimn left a comment

Choose a reason for hiding this comment

keflavich commented Mar 13, 2023

keflavich commented Mar 13, 2023

keflavich commented Mar 13, 2023

keflavich commented Mar 13, 2023

keflavich commented Mar 13, 2023

pllim commented Mar 13, 2023

dhomeier commented Mar 13, 2023

pllim commented Mar 13, 2023

keflavich commented Mar 13, 2023

keflavich commented Mar 13, 2023

pllim commented Mar 13, 2023

keflavich commented Mar 13, 2023

pllim commented Mar 13, 2023

mhvk commented Mar 13, 2023

keflavich commented Mar 13, 2023

mhvk commented Mar 13, 2023

dhomeier commented Mar 13, 2023

dhomeier commented Mar 13, 2023 • edited

dhomeier commented Mar 13, 2023 • edited

kYwzor commented Mar 9, 2023 •

edited by pllim

kYwzor Mar 9, 2023 •

edited

kYwzor Mar 9, 2023 •

edited

dhomeier commented Mar 13, 2023 •

edited

dhomeier commented Mar 13, 2023 •

edited