gammapy / gammapy

A Python package for gamma-ray astronomy
https://gammapy.org
BSD 3-Clause "New" or "Revised" License
223 stars 194 forks source link

Access time of events through DataStoreObservation.events #922

Closed registerrier closed 7 years ago

registerrier commented 7 years ago

I am trying to understand why image_pipe is slow. While profiling the code, I saw that most of the time is spent in access to events and irfs.

If one tests the access time to the event file, one gets:

In [42]: %timeit data_store.obs(23523).events 1 loops, best of 3: 4.06 s per loop

In [43]: %timeit EventList.read("/Users/terrier/hess/fits/ash_stereo_thsq64/run023400-023599/run023523/events_023523.fits.gz") 10 loops, best of 3: 58.6 ms per loop

Why is there a factor 100 between the 2 access times? This seems to be a very serious bottleneck, no?

cdeil commented 7 years ago

@registerrier - Thanks for identifying this issue! I'll have a look on Monday.

adonath commented 7 years ago

@registerrier @cdeil

I can't reproduce this on my machine. I used the following code:

from gammapy.data import DataStore, EventList
data_store = DataStore.from_dir('$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2')

%timeit data_store.obs(23523).events
%timeit EventList.read("$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/run023400-023599/run023523/hess_events_023523.fits.gz")

The output is:

10 loops, best of 3: 36.5 ms per loop
10 loops, best of 3: 36 ms per loop
registerrier commented 7 years ago

This is because the hdu_table in GAMMAP_EXTRA is tiny. You should use a larger one e.g. HESS I data.

adonath commented 7 years ago

@registerrier Thanks, now I can reproduce it.

adonath commented 7 years ago

I used %prun to find the bottleneck:

        576668 function calls (574742 primitive calls) in 0.369 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.112    0.112    0.292    0.292 hdu_index_table.py:208(row_idx)
   124426    0.097    0.000    0.180    0.000 table.py:1194(__getitem__)
   124431    0.063    0.000    0.072    0.000 table.py:98(__getitem__)
   256366    0.021    0.000    0.021    0.000 {built-in method builtins.isinstance}
        4    0.006    0.001    0.006    0.001 data.py:111(get_readable_fileobj)
     1215    0.005    0.000    0.019    0.000 configuration.py:380(__call__)
     1484    0.004    0.000    0.007    0.000 configuration.py:481(get_config)
     2215    0.003    0.000    0.003    0.000 {method 'match' of '_sre.SRE_Pattern' objects}
       51    0.002    0.000    0.002    0.000 {method 'decompress' of 'zlib.Decompress' objects}
     2699    0.002    0.000    0.002    0.000 configobj.py:554(__getitem__)
      875    0.001    0.000    0.023    0.000 card.py:271(value)
     1215    0.001    0.000    0.001    0.000 validate.py:638(_parse_with_caching)

It seems like the majority of time is spend in the python loop in HDUIndexTable.row_index(). Maybe this should be reimplemented using some kind of dictionary lookup, that is only created once on init. Or using numpy indexing...

cdeil commented 7 years ago

@adonath - Thanks for all your work to get the Gammapy issues resolved!!!

Here's some feedback. In the end, please put what you think best, no need to discuss.

For catalog we have a similar issue. Looking at the implementation, it's quite complicated. I think at the core it's a dict cache for catalog.

And for the HDU index table, both what we have now and what you put in #945 is different, not a dict. But there are things that are cached, so in the end both for catalog and HDU index, the table should be treated as read-only. (I don't understand if there is a working cache invalidation check or not for catalog.)

Wouldn't it be simpler to have a dict cache in both cases? Keys can be tuples in Python. That would probably be fastest (and maybe also simplest in terms of code?

@adonath - If you don't want to put a dict, let me know, and I'll briefly review the current implementation you put in #945 .

cdeil commented 7 years ago

Just thought of one more possible implementation: Table itself has the option to put an index column. If there are multiple keys, one can make a hash function for both and put that as index column. I think we had this at some point, but removed it due to a bug in Table indexing, which should be no longer present with Astropy 1.3.

So @adonath - Put whatever you like for now to resolve the speed issue. And if you think it's not the best long-term solution, please open a reminder issue what you think should be put eventually, once someone has time.

adonath commented 7 years ago

Resolved in #945.