Skip to content

Improve FITS header performance - #8502

Merged
astrofrog merged 20 commits into
astropy:masterfrom
saimn:fits-perf-v2
Apr 19, 2019
Merged

astrofrog merged 20 commits into
astropy:masterfrom
saimn:fits-perf-v2

Conversation

@saimn

@saimn saimn commented Mar 16, 2019

Copy link
Copy Markdown
Contributor

This PR is to optimize the initial parsing of FITS headers, as described in #5593 (which could probably be closed then). As discussed in #5593, a long-standing problem with io.fits is that to access a given HDU all the previous HDU headers must be parsed. And before the lazy loading feature, all HDU headers were parsed. So this is a crucial issue for FITS files with many HDUs, and often a barrier to use them.

The main issue is that the headers must be parsed to find the structural keywords (BITPIX, NAXIS*) that are needed to know the position of the following extension. After having tried to optimize the header parsing itself (#8428), which is always good, I came to the same conclusion as @embray (#5593 (comment)): we need a fast header parser. Which this PR implements.

The _BasicHeader class is added with a "fast" parser written in Cython. It is not as fast as I would like, only twice faster as the pure Python. I think that it is possible to do better but will keep this for the future, it is already much faster that Header. This new class is used only during the construction of the HDU classes, where we need only a subset of standard keywords, so COMMENT, HIERARCH, and special keywords are not parsed. The cards are parsed only if needed, and cached. Typically only a few cards will be parsed, during the HDU class detection step. Then the header string is stored, and used to construct the usual Header object if the HDU is accessed.

One remains that remains is to know if we need to parse to full header (currently implemented), or just the first block, which would bring another interesting speedup for headers with many cards (see #5593 (comment) and following comments). We could maybe imagine having a keyword to trigger the full header parsing ?

To give some numbers, and testing on different files:

astropy/io/fits/tests/data/stddata.fits, with 2 tables and just a few
keywords, to check if there is a penalty. The PR is still faster.

❯ python -m perf timeit --fast -s "from astropy.io import fits" "fits.getheader('astropy/io/fits/tests/data/stddata.fits', ext=1)"

Before: Mean +- std dev: 1.28 ms +- 0.03 ms
After: Mean +- std dev: 1.17 ms +- 0.02 ms

❯ python -m perf timeit --fast -s "from astropy.io import fits" "fits.getheader('astropy/io/fits/tests/data/stddata.fits', ext=2)"

Before: Mean +- std dev: 2.88 ms +- 0.10 ms
After: Mean +- std dev: 2.30 ms +- 0.08 ms

Then on a quite extreme case, a file with 335 extensions, where each header
has ~1700 cards (with HIERARCH):

❯ python -m perf timeit --fast -s "from astropy.io import fits" "fits.getheader('scipost_white.fits', ext=10)"

Before: Mean +- std dev: 173 ms +- 6 ms
After: Mean +- std dev: 26.8 ms +- 1.2 ms

❯ python -m perf timeit --fast -s "from astropy.io import fits" "fits.getheader('scipost_white.fits', ext=60)"

Before: Mean +- std dev: 1.01 sec +- 0.01 sec
After: Mean +- std dev: 101 ms +- 3 ms

❯ python -m perf timeit --fast -s "from astropy.io import fits" "fits.getheader('scipost_white.fits', ext=300)"

Before: Mean +- std dev: 5.26 sec +- 0.16 sec
After: Mean +- std dev: 450 ms +- 10 ms

fits_perf_scipost

And finally with a Euclid test file, with 650 extensions but a small number of keywords (15 to 20). It seems that they have more critical cases, and use fitsio instead of Astropy because of that.

fits_perf_euclid

ping @embray - if you have some time to have a look, it would be great to have your comments!

@saimn saimn added this to the v3.2 milestone Mar 16, 2019
@saimn saimn assigned drdavella and MSeifert04 and unassigned drdavella and MSeifert04 Mar 16, 2019
@saimn
saimn requested review from MSeifert04 and drdavella March 16, 2019 00:09
Comment thread astropy/io/fits/card.py
# probably a parsing error, falling back to the internal _value
# which should be None. This may happen while calling _fix_value.
oldvalue = self._value

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

About these changes, see b135229

"""

klass = _hdu_class_from_header(cls, header)
return super().__new__(klass)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

See a9db13e, this was done twice.

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.

This is still necessary if you want to instantiate an HDU instance directly (e.g. with type.__call__) and have it return the correct HDU class based on the header. I don'e think that the existence of _readfrom_internal obviates its need because that's only used in alternative constructors for reading an HDU from data _BaseHDU.fromstring and _BaseHDU.readfrom. That's different from calling _BaseHDU.__new__ directly.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Yeah I realized this later with #8660, as some convenience methods relies on this. In #8666 I added a new class method to use the HDU class detection, i.e. _BaseHDU._from_data(data, header) instead of relying on _BaseHDU.__call__. I think it's a good compromise, it is no more possible to have an HDU class change is type based on the header which is a good thing imo, and we avoid the duplicated detection in the general case.

@codecov

codecov Bot commented Mar 16, 2019 •

Copy link
Copy Markdown

Codecov Report

Merging #8502 into master will increase coverage by <.01%.
The diff coverage is 92.7%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master    #8502      +/-   ##
==========================================
+ Coverage   86.82%   86.83%   +<.01%     
==========================================
  Files         386      386              
  Lines       58197    58267      +70     
  Branches     1060     1060              
==========================================
+ Hits        50531    50595      +64     
- Misses       7051     7057       +6     
  Partials      615      615
Impacted Files Coverage Δ
astropy/io/fits/hdu/image.py 96.04% <ø> (+0.17%) ⬆️
astropy/io/fits/hdu/table.py 93.01% <ø> (+0.13%) ⬆️
astropy/io/fits/card.py 86.47% <81.81%> (-0.22%) ⬇️
astropy/io/fits/hdu/base.py 90.31% <92.3%> (-0.04%) ⬇️
astropy/io/fits/header.py 96% <94.44%> (-0.19%) ⬇️
astropy/samp/hub.py 77.38% <0%> (-0.14%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 33db190...667cbfb. Read the comment docs.

@astrofrog astrofrog left 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.

The performance improvements look great! I'm not qualified to review this properly, but just a couple of comments about the Cython code.

# the full Header parsing
raise Exception

block_str = block.decode('ascii')

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.

Aren't headers supposed to be ASCII?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Yes, they are supposed to. But sometimes they contain unicode chars, which will fail here, and fallback to Header, which then use decode_ascii.
One other option, that I tried, would be to keep working on bytes here. And store the header string as bytes, which would save some memory and avoid decoding headers that are not used later. But there are some complications, like the need to check for padding bytes, and the gain was null... but I will probably try again this later.

Comment thread astropy/io/fits/_utils.pyx Outdated
idx = end_idx

keyword = card_image[:8].strip()
sep_idx = card_image.find(VALUE_INDICATOR)

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.

Ins't the = separator always supposed to be in a specific character position? find could be slow... Ah I see what comes below. What about, to speed things up in the most common case, making the start of the 'if' statement if card_image[8] == '=' and only using find if that fails?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Good idea, I will try this.

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.

find supports start and stop, maybe we can use these to limit the "searched region": card_image.find(VALUE_INDICATOR, 0, 9) (not sure if the indices are correct ... 😓)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Good to know. I already implemented something with @astrofrog 's suggestion (and forgot to push it ...), where I slice the first 8 characters before the .find. Not sure if it makes a big difference with using the start and stop parameters.

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.

Not sure if it makes a big difference with using the start and stop parameters.

It's probably faster because it avoids the copy. Probably not big but worth a try

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

It seems to make a small difference, always good to take 👍

Comment thread astropy/io/fits/card.py
elif self._valuestring is not None or self._image:
self._value = self._parse_value()
value = self._value
value = self._value = self._parse_value()

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.

I would prefer the old style here.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Why, multiple assignment is cool :p

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 agree with @saimn!

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.

😱

@MSeifert04 MSeifert04 Mar 18, 2019 •

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.

Maybe it's just me, but it feels like an unintended misspelling of value = self._value == self._parse_value(). Also it reminds me of C (not in a good way 😄).

Don't get me wrong, I won't start a fight about it - but this looks so weird to me because it's so rare that I see this in Python and since assignment is normally a statement not an expression it looks like it really shouldn't work.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I don't know, I frequently use this, typically to get a shorter variable name for a local block and avoid the lookup time on self. It's quite practical, and imo easier to understand that it references the same variable.

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.

Okay, so it's probably just a matter of getting used to it. I think we probably have to get used to it when Python 3.8 with PEP 572 -- Assignment Expressions comes.

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.

wow, I should read more PEP, this looks quite interesting (though I agree that feels very unusual)

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.

Ah yes, PEP 572, the one that led to Guido stepping down from BDFL! I thought it was strange at first but I constantly find myself needing it.

Comment thread astropy/io/fits/hdu/base.py Resolved
Comment thread astropy/io/fits/_utils.pyx Resolved
Comment thread astropy/io/fits/_utils.pyx Resolved
Comment thread astropy/io/fits/_utils.pyx Outdated
"""The main method allowing to parse quickly a FITS header."""

cards = OrderedDict()
cdef read_blocks = []

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.

I actually though cdef list - it's amazing that it worked 😮

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Oups! It wasn't done on purpose but yeah it seems that cython understand this. The generated code is slightly different:

+14:     cdef read_blocks = []

  __pyx_t_1 = PyList_New(0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 14, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_read_blocks = __pyx_t_1;
  __pyx_t_1 = 0;
...
+29:         read_blocks.append(block_str)

    __pyx_t_9 = __Pyx_PyObject_Append(__pyx_v_read_blocks, __pyx_v_block_str); if (unlikely(__pyx_t_9 == ((int)-1))) __PYX_ERR(0, 29, __pyx_L1_error)

and with cdef list:

+14:     cdef list read_blocks = []

  __pyx_t_1 = PyList_New(0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 14, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_read_blocks = ((PyObject*)__pyx_t_1);
  __pyx_t_1 = 0;
...
+29:         read_blocks.append(block_str)

    __pyx_t_9 = __Pyx_PyList_Append(__pyx_v_read_blocks, __pyx_v_block_str); if (unlikely(__pyx_t_9 == ((int)-1))) __PYX_ERR(0, 29, __pyx_L1_error)

@MSeifert04

Copy link
Copy Markdown
Contributor

Looks mostly good to me.

What bothers me is the performance, I would've expected it to be significantly faster, maybe it's all the string operations (slicing and joining) that slow it down. Did you do any profilings?

Since the tests are passing it's probably okay. One thing we really should do is add some docs on how this internally works.

Oh and I think some of the refactoring stuff would be good to do in separate PRs. Makes the changes a lot easier to review (for example changing the imports and the indentation of arguments/.format calls).

Comment thread astropy/io/fits/card.py
self._value = None
self._comment = None

self._valuestring = None

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 was this moved?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Because it must be defined before self.value is set. Using the new header class triggered this bug (b135229) in the getter/setter loop for .value, and the same for .comment.

Comment thread astropy/io/fits/hdu/base.py Resolved
if not block or len(block) < BLOCK_SIZE:
# header looks incorrect, raising exception to fall back to
# the full Header parsing
raise Exception

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.

Seems odd to raise a plain Exception. That way we probably swallow any exceptions that we don't expect when catching it.

I don't have a good idea though what makes sense here. Maybe ValueError or some custom exception like FastHeaderParsingFailedException? 😄

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

In any case there is a fallback to the Header parsing. I'm fine with using a custom exception but this will not make a difference, we need to catch any exception in the fast reader and then call Header as it have a better reporting of errors.

@saimn

saimn commented Mar 20, 2019

Copy link
Copy Markdown
Contributor Author

@MSeifert04 :

What bothers me is the performance, I would've expected it to be significantly faster, maybe it's all the string operations (slicing and joining) that slow it down. Did you do any profilings?

Yeah I would have expected a better performance as well, on the mockup I did it is much faster, and I don't know why. I tried several variants, keeping bytes, and nothing really helps. I think it is not critical for now, as even on small files it is faster, so I would prefer to move forward with this version and try to optimize later. Any idea or help to optimize will be greatly appreciated ;).

Since the tests are passing it's probably okay. One thing we really should do is add some docs on how this internally works.

I don't think we currently have such docs? The additional complexity is rather limited compared to the mean complexity of io.fits 🙂, the fast header is built for a short time, all in _BaseHDU._readfrom_internal, and its implementation is relatively straightforward.

Oh and I think some of the refactoring stuff would be good to do in separate PRs. Makes the changes a lot easier to review (for example changing the imports and the indentation of arguments/.format calls).

I reverted the two import changes.

if 0 < sep_idx < 8:
keyword = card_image[:sep_idx]
cards[keyword.upper()] = card_image
elif card_image == END_CARD:

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.

Maybe the END check should come before the find call? That is, first check if card_image[8:10] is the value indicator, then check if the card is the end card, then do the thing with find?

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?

@mhvk

mhvk commented Mar 23, 2019

Copy link
Copy Markdown
Contributor

@saimn - overall, this surely is the right way to go, but some questions. First, how much does the Cython buy you? It makes working on it harder for others, so it should be significant, yet since this is mostly string operations, it is not obvious to me that pre-compiling is going to help much.

Second, have you checked whether basing yourself on abc.Mapping is costly? I have no experience with this particular class, but did find that ABCs in general are not cheap: for every creation of an instance, it is checked whether the relevant methods are present (it cannot do this at the class creation level). Perhaps worth a try to just remove it?

Though really I guess some detailed profiling is needed - which perhaps is better done as follow up.

@MSeifert04

Copy link
Copy Markdown
Contributor

overall, this surely is the right way to go, but some questions. First, how much does the Cython buy you? It makes working on it harder for others, so it should be significant, yet since this is mostly string operations, it is not obvious to me that pre-compiling is going to help much.

Interesting point. I was assuming that most of the speedup is gained because of Cython. But you're right the potentially slow parts are string operations that will be (at least likely) equally fast in pure Python.

Though really I guess some detailed profiling is needed - which perhaps is better done as follow up.

👍 - but I am not sure if we should do it as follow-up, since this is adding a lot of complexity and "only" gives a factor 2 improvement it should be essential that we're on the right track.

@mhvk

mhvk commented Mar 24, 2019

Copy link
Copy Markdown
Contributor

It's true that most typical use cases with many extensions probably involve small headers, not the 1700 card ones that your PR really shines on... Also, your tests assume a user in fact knows which extension they want, where perhaps more realistically (some of) the other extensions also get looked at. Is there currently a cache for extensions, so that, e.g., if I look for extension 80 after having looked at 70, one starts at 70 and not back at 1?

@astrofrog

Copy link
Copy Markdown
Member

Just a quick note that it would probably make sense to add several benchmarks to https://github.com/astropy/astropy-benchmarks for this - maybe one simple case with one HDU, and then a couple of more complex examples?

@saimn

saimn commented Mar 25, 2019

Copy link
Copy Markdown
Contributor Author

@astrofrog - Yes, was working on it ;), see astropy/astropy-benchmarks#80

saimn added 17 commits April 18, 2019 23:14
Illustration of the issue, setting a value to a Card that was not yet
parsed change its _modified status though it should not.
```
In [1]: %astropy
Numpy 1.16.1
Astropy 3.2.dev23834

In [2]: card = fits.Card.fromstring('FOO    = 10')

In [3]: card.value = 10

In [4]: card._modified
Out[4]: True

In [5]: card = fits.Card.fromstring('FOO    = 10')

In [6]: card.value
Out[6]: 10

In [7]: card.value = 10

In [8]: card._modified
Out[8]: False
```
The HDU class detection is already done in _readfrom_internal, I cannot
see why doing it twice is necessary. And no test fail because of this :)
@saimn

saimn commented Apr 18, 2019

Copy link
Copy Markdown
Contributor Author

Rebased, please have a look @MSeifert04, @astrofrog, @mhvk :)

@astrofrog astrofrog left 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.

Looks great!

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

Projects

None yet

Development

Successfully merging this pull request may close these issues.

8 participants