Skip to content

A new tiled image compression submodule mostly replacing cfitsio - #14252

Merged
saimn merged 54 commits into
astropy:mainfrom
aperiosoftware:fits-byte-compress-decompress
Jan 20, 2023
Merged

saimn merged 54 commits into
astropy:mainfrom
aperiosoftware:fits-byte-compress-decompress

Conversation

@Cadair

@Cadair Cadair commented Jan 4, 2023 •

Copy link
Copy Markdown
Member

This pull request contains phase one of the work proposed in the cycle 3 funding call, and has been a collaboration between @astrofrog and I.

Background

The FITS file format provides a way of compressing multi-dimensional images to reduce their size on disk. This functionality works by splitting the image up into "tiles" compressing them using one of a selection of compression algorithms and then saving the compressed bytes into a cell in a FITS binary table.

In Astropy this functionality is exposed by CompImageHDU. When loading the data associated with a compressed image, the whole HDU is passed through to a function implemented in C called decompress_hdu (the same happens on compress when writing). This HDU object is then used to pass the whole binary table (and it's heap) to the cfitsio library which returns a complete decompressed array of the image. (Again the same happens for compress but a binary table is returned).

This functionality is the only reason astropy uses cfitsio, yes all 168,643 lines of it. Implementing this functionality this way also means that it is impossible to load and decompress only part of the image, the whole image is always fully loaded.

Description of this PR

This pull request is a big step towards implementing all the FITS tiled image compression in Python. All the functionality of decompress_hdu and compress_hdu currently provided in compressionmodule.c is replaced by equivalent Python code. The compression and decompression of the individual tiles are implemented in the style of numcodecs Codec classes (more on this below). The implementation of the compression and decompression algorithms for RICE_1, PLIO_1 and HCompress are still used from the cfitsio code, but only a small number of files from the cfitsio library are now needed.

While in this PR we do not expose any functionality for (de)compressing individual tiles, the primitives for this work now exist, as it is possible to use these Codec classes to (de)compress the bytes of an individual tile, and the C wrappers for the cfitsio functions only handle (de)compressing a single tile at a time.

Testing

To validate the implementation in this PR we needed to be able to test against FITS files not (de)compressed with Astropy. We approached this in two ways.

  1. Some canonical test files downloaded from the FITS compression pages have been added to the repo. These are small, but provide a useful "known good" file for each compression algorithm.
  2. A new heavily parametrised test_fitsio.py file has been added. This file uses the fitsio package to write or read files read or written by Astropy and then the data is compared. Running this test file (obviously) needs the fitsio package to be installed, which requires a system install of cfitsio to compile the sdist against. This test file is automatically skipped if fitsio can not be imported, but a tox factor has been added to force it to be run and to require fitsio to ensure it is always run on the CI.
    2a) For extra credit and to hunt segfaults 🔫, all these tests have been run with valgrind and no memory leaks in our code (just fitsio) were found.

New Functionality or Bug Fixes in this PR

This PR aimed to change nothing as far as users of the library are concerned. Unless someone was importing directly from astropy.io.fits.compression (which I really hope they weren't) then this change should go unnoticed. However, in the process of rewriting the compression functionality we found some things which were missing from Astropy or the implementations in cfitsio these things are:

Next Steps

This PR could have removed almost all of cextern/cfitsio but we decided that we shouldn't for the initial review (the diff is obnoxious). We can choose to either do this as the last thing before merging this PR, or in a follow up PR. We did have to already remove some code from imcompress.c here in order to get things to work.

Related Issues

This pull request addresses the first checklist item in #3895 and opens the door to some (or all) of the others.

Checklist for package maintainer(s)

This checklist is meant to remind the package maintainer(s) who will review this pull request of some common things to look for. This list is not exhaustive.

  • 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 "When to rebase and squash commits".
  • 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?
  • Is a milestone set? Milestone must be set but astropy-bot check might be missing; do not let the green checkmark fool you.
  • 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.

@Cadair
Cadair requested a review from saimn as a code owner January 4, 2023 18:49
@github-actions github-actions Bot added external PRs and issues related to external packages vendored with Astropy (astropy.extern) io.fits testing labels Jan 4, 2023

@mhvk mhvk left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Really nice. Just had a quick look as I was interested, and noticed that at least some use of .tobytes() is not needed, i.e., one can avoid the associated copies. I suspect that in cases where np.ndarray are not directly useful, a memoryview may well suffice.

buf
The decompressed buffer.
"""
cbytes = np.frombuffer(buf, dtype=np.uint8).tobytes()

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You don't need .tobytes() here (which makes a copy); gzip_decompress can handle the numpy buffer (or indeed any buffer, not sure you need even np.frombuffer, but at least that does not make a copy).

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's very weird... Frankly sounds like some kind of bug. It may be worth trying memoryview(buf) or memoryview(np.frombuffer, ...) which would also not make a copy.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The issue happens with compression - here's a MWE in case you are interested in following this up on the Numpy side:

In [1]: import numpy as np

In [2]: x = np.arange(9).reshape((3, 3)).astype(np.uint8)

In [3]: from gzip import compress

In [4]: compress(x)
Out[4]: b'\x1f\x8b\x08\x00Dk\xc6c\x02\xffc`dbfaec\xe7\x00\x00\x02C\xe1\xbc\x03\x00\x00\x00'

In [5]: compress(x.tobytes())
Out[5]: b'\x1f\x8b\x08\x00Lk\xc6c\x02\xffc`dbfaec\xe7\x00\x00\x02C\xe1\xbc\t\x00\x00\x00'

In [6]: from gzip import decompress

In [7]: decompress(compress(x))
---------------------------------------------------------------------------
BadGzipFile                               Traceback (most recent call last)
Cell In[7], line 1
----> 1 decompress(compress(x))

File /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/gzip.py:614, in decompress(data)
    612     raise BadGzipFile("CRC check failed")
    613 if length != (len(decompressed) & 0xffffffff):
--> 614     raise BadGzipFile("Incorrect length of data produced")
    615 decompressed_members.append(decompressed)
    616 data = do.unused_data[8:].lstrip(b"\x00")

BadGzipFile: Incorrect length of data produced

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is very odd. It may be a gzip that is responsible since also memoryview(x) fails. And zlib works fine.

I note that in the documentation of gzip, it states "Changed in version 3.11: Speed is improved by compressing all data at once instead of in a streamed fashion. Calls with mtime set to 0 are delegated to zlib.compress() for better speed." (and indeed, with mtime=0 it works)

I wonder if the "all data at once" is a clue here... But regardless, your example shows that it has something to do with strides and/or shape, and I worry that perhaps that might cause problems if the data are not contiguous. Indeed, even on python 3.10,

In [7]: compress(x[:, 0])
---------------------------------------------------------------------------
BufferError                               Traceback (most recent call last)
Cell In [7], line 1
----> 1 compress(x[:, 0])

File /usr/lib/python3.10/gzip.py:549, in compress(data, compresslevel, mtime)
    547 buf = io.BytesIO()
    548 with GzipFile(fileobj=buf, mode='wb', compresslevel=compresslevel, mtime=mtime) as f:
--> 549     f.write(data)
    550 return buf.getvalue()

File /usr/lib/python3.10/gzip.py:289, in GzipFile.write(self, data)
    286     length = data.nbytes
    288 if length > 0:
--> 289     self.fileobj.write(self.compress.compress(data))
    290     self.size += length
    291     self.crc = zlib.crc32(data, self.crc)

BufferError: memoryview: underlying buffer is not C-contiguous

So, perhaps the copy is worth it!?

For now, perhaps just add a comment...

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added comments in 4690dd6 (it turns out we need it for both compression and decompression)

Comment thread astropy/io/fits/_tiled_compression/codecs.py

@saimn saimn left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A first batch of comments, but I haven't really looked at the tiled_compression module yet.
Looks awesome, great work !

except ImportError:
COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = False

COMPRESSION_ENABLED = True

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why keeping this var ?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because HDUList does some horrible things using this:

compressed.COMPRESSION_ENABLED = False

Basically if you open the FITS file with disable_image_compression it will temporarily modify the global variable COMPRESSION_ENABLED. We should definitely refactor this, but for now we were trying to not refactor too much.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've now added a comment above where COMPRESSION_ENABLED is defined to explain this.

Comment thread .github/workflows/ci_workflows.yml Outdated
Comment thread .github/workflows/ci_workflows.yml Outdated
Comment thread .github/workflows/ci_workflows.yml Outdated
Comment thread astropy/io/fits/_tiled_compression/codecs.py
Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py Outdated
Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py Outdated
Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py Outdated
Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py
Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py Outdated
@astrofrog
astrofrog force-pushed the fits-byte-compress-decompress branch from 9d19273 to f6e0ce9 Compare January 16, 2023 11:21
@astrofrog

Copy link
Copy Markdown
Member

@saimn - I think I've addressed your comments so far and have tried to tidy things up a bit in places.

Once you think this is ready (which I know is probably not yet!) let me know and I can remove all the unused bundled CFITSIO code.

@astrofrog

Copy link
Copy Markdown
Member

For some reason there is an issue with gzip decompression on Python 3.11, so will need to investigate that. There are also some docs warnings that will need to be fixed.

@astrofrog

Copy link
Copy Markdown
Member

@mhvk - it seems that in Python 3.11 I need the .tobytes() in the encode method for gzip compression, otherwise things fail. I haven't had a chance to track down exactly why that is happening but just thought I'd let you know in case you have any thoughts.

@astrofrog astrofrog added Build all wheels Run all the wheel builds rather than just a selection Extra CI Run cron CI as part of PR labels Jan 16, 2023
@astrofrog

Copy link
Copy Markdown
Member

There are real failures here on s390x with PLIO, likely due to that platform being big endian, so will need to investigate that. Hopefully RTD and pre-commit will now pass.

@pllim pllim added this to the v6.0 milestone Jan 17, 2023
@astrofrog

Copy link
Copy Markdown
Member

@pllim isn't the next release 5.3? (milestone wise)

@astrofrog

Copy link
Copy Markdown
Member

@nstarman - ruff keeps trying to 'fix' this:

Screenshot 2023-01-17 at 21 53 51

What is the correct way to tell it to not touch these lines? (also I guess this is a ruff bug?)

@astrofrog
astrofrog force-pushed the fits-byte-compress-decompress branch from 693cc63 to 3f5ff25 Compare January 17, 2023 22:02
@mhvk

mhvk commented Jan 17, 2023

Copy link
Copy Markdown
Contributor

Funny, that would be correct if that line was the else clause...

@saimn saimn left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few more comments, nothing critical :)

Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py Outdated
Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py Outdated
Comment thread astropy/io/fits/_tiled_compression/tiled_compression.py
@nstarman

Copy link
Copy Markdown
Member

I'll address this more fully when I'm at my computer. For now, to find what code to ignore in a noqa, comment out -fix in pre-commit config and rerun. It will give the error code.

@nstarman

nstarman commented Jan 18, 2023 •

Copy link
Copy Markdown
Member

Back home and had a chance to look. So it looks like it's SIM401 | DictGetWithDefault | Use var = dict.get(key, "default") instead of an if block | 🛠
that ruff is trying to fix.

Our options are 3:

  1. ruff is (mostly) correct and the code should be reformatted. The else block should be removed and then we're good here, but should submit an Issue to ruff about the mistake.
  2. We just want to ignore this case, so # noqa: SIM401 should make ruff ignore that line and not auto-fix. (Probably still want to alert ruff)
  3. we turn off SIM401

I'm not so familiar with this bit of the code, but from cursory examination it looks like _header has a .get method, so option 1 is correct?

]
)

cfg["include_dirs"].append(SRC_DIR)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line should be indented I think, i.e. should be done only when not using the system cfitsio.

@saimn saimn Jan 19, 2023 •

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah actually not since you need the additional headers defined in astropy. But then some header files (fitsio2.h typically) can be found from two different locations when using system cfitsio. Not sure if that can cause problems, e.g. if system cfitsio's version is different from astropy's one ?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(from the sidelines so possibly irrelevant) - I think if astropy provides additional headers, the solution would be to have those in a different file. It may be good to check with @olebole on the system file use.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can also only comment from aside, having checked the few lines above only. My understanding of these lines is that astropy is actually not providing the include files in SRC_DIR, but uses them. In this case, I agree with @saimn that this should only happen when the astropy-provided cfitsio package is used and not a system one.

This is also the case if they are provided to other packages: if astropy was built against a system cfitsio, other packages should refer to the system API, and therefore the system include file if they are going to be linked to the same library as astropy (or otherwise uses it).

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

... but it seems that this PR proposes to use a different shared library than cfitsio (and therefor we can drop the "use-system-lib" here). In this case, I think it would be better to use another name (or location) for the include file (if it is provided to others); otherwise there would be an ambiguity and the case "I want to link to both cfitsio and astropy's module" would be very difficult.
Having a final location that allows to write

#include "astropy/fitsio2.h"

would be already sufficient.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@olebole - just to clarify:

but it seems that this PR proposes to use a different shared library than cfitsio

this isn't quite right - we are still using cftisio, just a much smaller fraction of it.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@saimn - and just going back to your original comment for this discussion, it makes sense to now keep our include dir for both the case when using our bundled CFITSIO and the system install because the header files now no longer duplicate anything in the CFITSIO library.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry I wanted to make another comment yesterday but did not have time to do so.
So actually I'm wondering if allowing to use the system cfitsio still makes sense. Mostly because we are using custom headers to export some functions that cfitsio does not export. So there is no guarantee that it will work with other versions (past or future) of cfitsio, in a sense these functions are internal and their API might change ?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In my opinion, it doesn't make sense to keep using system cfitsio here. There is no need to do this overhead for just a few lines of code.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't fully understand what is considered public/private in a C library so happy to defer to you both. It might be easier to do remove in a follow-up PR as I need to update the docs etc.

@astrofrog

astrofrog commented Jan 19, 2023 •

Copy link
Copy Markdown
Member

@saimn - I think I've addressed your latest comment about the include files!

@pllim

pllim commented Jan 20, 2023

Copy link
Copy Markdown
Member

So, you deleted cextern/trim_cfitsio.sh, so how is one supposed to grab a new CFITSIO release to be bundled since you said we are still using parts of it?

@astrofrog

Copy link
Copy Markdown
Member

Good point I will add it back and get it to copy the specific files we want instead.

@astrofrog

Copy link
Copy Markdown
Member

@pllim - I have added back the script and run it with the latest version of cfitsio to check (no updates to any of the C files as expected)

@astrofrog
astrofrog requested a review from saimn January 20, 2023 20:12
@olebole

olebole commented Jan 20, 2023

Copy link
Copy Markdown
Member

I have one late comment; not sure if it is a useful one 😀
There is still the fitsio package, which is a Python wrapper around cfitsio. Long time ago, there was the idea to base astropy.io.fits on that, so that one could also have an integrated "fast but raw" access if needed. I liked that approach, but it seems that nobody started implementing this.

This PR will probably put this idea definitely into the past.

@astrofrog

Copy link
Copy Markdown
Member

@olebole - not a bad idea, but one of the biggest blockers to being to do that is that fitsio is GPL so this would require either convincing its authors to change their license or changing astropy to GPL. The other big issue is that fitsio requires users to install cfitsio themselves, which is not necessarily trivial on all platforms. To some extent, users who care about the best possible performance right now could simply use the fitsio package as-is and not even bother with astropy.io.fits.

@astrofrog

astrofrog commented Jan 20, 2023 •

Copy link
Copy Markdown
Member

@saimn - since this uses a lot of CI resources currently because of the extra CI and all wheels, if this is ready aside from removing the option to build against the system CFITSIO, can we merge and I can then do a follow-up PR to remove the system build option (updating docs etc).

@saimn

saimn commented Jan 20, 2023

Copy link
Copy Markdown
Contributor

@astrofrog - Yep, sounds good to me. Thanks (you and @Cadair !) for the great work :)

RTD is failing but seems to be due to tutorials, maybe transient ?

embedding documentation hyperlinks...
WARNING: The following HTTPError has occurred fetching https://numpy.org/devdocs//: 404 (Not Found)

@Cadair

Cadair commented Jan 24, 2023

Copy link
Copy Markdown
Member Author

Whoop 🎉 Thanks for the reviews and comments everyone ❤️

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Build all wheels Run all the wheel builds rather than just a selection external PRs and issues related to external packages vendored with Astropy (astropy.extern) Extra CI Run cron CI as part of PR io.fits Ready-for-final-review testing

Projects

None yet

7 participants