New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement TS map computation in Cython #252

Merged
merged 2 commits into from Apr 13, 2015

Conversation

Projects
None yet
3 participants
@adonath
Member

adonath commented Apr 6, 2015

This implements cash and model functions for TS computation in cython. Additionally it adresses a few minor issues (storing the sign of the TS values correctly, computing sqrt(TS)).

Also removes overlap.py.

@@ -16,7 +16,7 @@
from astropy.nddata.utils import extract_array
from astropy.io import fits
from ..stats import cash
from _test_statistics import _cash, __amplitude_bounds, _cash_sum, _f_cash_root

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

This fails on Python 3:
https://travis-ci.org/gammapy/gammapy/jobs/57322817#L1785
I didn't look in detail, but maybe it's just that you need to use a relative import here?

@cython.cdivision(True)
@cython.boundscheck(False)
def __amplitude_bounds(np.ndarray[np.float_t, ndim=2] counts,

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

Why double underscore?
I know single underscore means "private" ... if this is roughly the same can you use single underscore here, too?

sum = 0
for j in range(nj):
for i in range(ni):
sum += (model[j, i] * (counts[j, i] / (x * FLUX_FACTOR * model[j, i]

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

I always find such long formulae hard to read and prefer temp variables over line breaks.
Would that be OK with you here to split the formula in two lines and introduce a temp variable?
Speed should be the same because the C compiler optimises this, no?

@@ -307,6 +309,10 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
threshold : float (None)
If the TS value corresponding to the initial flux estimate is not above
this threshold, the optimizing step is omitted to save computing time.
sqrt : bool

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

Maybe a bit longer, more descriptive parameter name like return_sqrt_map?

@@ -307,6 +309,10 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
threshold : float (None)
If the TS value corresponding to the initial flux estimate is not above
this threshold, the optimizing step is omitted to save computing time.
sqrt : bool
Whether to return a sqrt(TS) map.
clean : bool

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

Maybe a bit longer, more descriptive parameter name like clean_edges or apply_edge_cleaning?

return TSMapResult(ts=TS, amplitude=amplitudes,
# Compute edge cleaned sqrt TS map
sqrt_TS = np.where(TS > 0, np.sqrt(TS), -np.sqrt(-TS))

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

Many users will wonder what a negative sqrt_TS means.
Can you give this formula in the function docstring or even make a section explaining it in the high-level docs and then reference to it from function docstrings that use this?

result[name] = upsample_2N(result[name], 2, order=order)
assert_allclose([[99], [99]], np.where(result.ts == result.ts.max()))
assert_allclose(1705.840212274973, result.ts[99, 99], rtol=1e-3)
assert_allclose(1675.174869581638, result.ts[99, 99], rtol=1e-3)

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

Do you understand these changes? Wild guess: is it because you use 32 bit floats in Cython?
Are the new numbers more correct or the old ones?

if len(results) > 1:
for scale, result in zip(scales, results):
filename_ = filename.replace('.fits', '_{0:.3f}.fits'.format(scale))
logging.info('Writing {0}'.format(os.path.join(folder, filename_)))

This comment has been minimized.

@cdeil

cdeil Apr 6, 2015

Member

No need to put {0}, you can just write {} because we don't support Python 2.6 and in Python 2.7 and 3 one doesn't have to give the index of the argument to put in str.format.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 6, 2015

@adonath – Thanks!

I know you've done benchmarking on this and using Cython is warranted here?
Can you include whatever script or notebook you have in https://github.com/gammapy/gammapy-benchmarks (maybe in a new notebooks folder) and link to it from here?
Or maybe it's even possible to write an asv benchmark for this that would work before / after this change and would show a drop in execution time with this PR? That would be a nice graph to have for future discussions where to use Cython in Gammapy.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 7, 2015

FYI: I just now wrote up the Gammapy policy for C / Cython / Numba:
https://gammapy.readthedocs.org/en/latest/development/index.html#when-to-use-c-or-cython-or-numba-for-speed

@adonath – OK?

@cdeil cdeil changed the title from Added Cython cash and model function and minor fixes to Implement TS map computation in Cython Apr 8, 2015

@cdeil cdeil added the feature label Apr 8, 2015

@cdeil cdeil added this to the 0.2 milestone Apr 8, 2015

@adonath adonath force-pushed the adonath:ts_map_performance_improvement branch from 6835cb2 to 03c9a82 Apr 13, 2015

@adonath

This comment has been minimized.

Member

adonath commented Apr 13, 2015

Maybe it would make sense to define a naming convention for cython functions? E.g. they should be private i.e. start with _ and then end with _cython (That's what I adopted for the TS functions...)

@adonath adonath force-pushed the adonath:ts_map_performance_improvement branch from 03c9a82 to 0b0d3b2 Apr 13, 2015

cimport numpy as np
cimport cython
cdef np.float_t FLUX_FACTOR

This comment has been minimized.

@cdeil

cdeil Apr 13, 2015

Member

Put 1e-12 on this line

-----
Negative :math:`TS` values are defined as following:
.. math::

This comment has been minimized.

@cdeil

cdeil Apr 13, 2015

Member

Dedent one level

amplitudes = np.zeros(TS.shape)
niter = np.zeros(TS.shape)
TS[j, i] = [_[0] for _ in results]
amplitudes[j, i] = [_[1] for _ in results]
niter[j, i] = [_[2] for _ in results]
return TSMapResult(ts=TS, amplitude=amplitudes,
# Compute negative TS values

This comment has been minimized.

@cdeil

cdeil Apr 13, 2015

Member

compute -> handle

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 13, 2015

I think the example notebook doesn't work any more:
https://gist.github.com/cdeil/f5217525788215fd7824
Can you try and if needed adapt to these changes?

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 13, 2015

Please also try out and update the example here:
https://gammapy.readthedocs.org/en/latest/detect/index.html#command-line-tool

$ gammapy-ts-image all.fits.gz --threshold 0 --scale 0
usage: gammapy-ts-image [-h] [--psf PSF] [--morphology MORPHOLOGY]
                        [--width WIDTH] [--scales SCALES [SCALES ...]]
                        [--downsample DOWNSAMPLE] [--residual] [--model MODEL]
                        [--threshold THRESHOLD] [--overwrite]
                        input_file output_file
gammapy-ts-image: error: the following arguments are required: output_file

@cdeil cdeil referenced this pull request Apr 13, 2015

Closed

Further TS map improvements #223

6 of 10 tasks complete
@cdeil

This comment has been minimized.

Member

cdeil commented Apr 13, 2015

@adonath Something unrelated I noticed while preparing the Gammapy 0.2 release.

This external link is broken:
http://home.strw.leidenuniv.nl/~oberg/Pulsars/Verbunt_Heise_NeutronStars.pdf
https://github.com/gammapy/gammapy/search?utf8=%E2%9C%93&q=verbunt

Could you please remove it or, if something equivalent exists, change the reference?

@adonath adonath force-pushed the adonath:ts_map_performance_improvement branch 2 times, most recently from e2a9353 to 295b04c Apr 13, 2015

@adonath adonath force-pushed the adonath:ts_map_performance_improvement branch from 295b04c to c4f90fb Apr 13, 2015

@coveralls

This comment has been minimized.

coveralls commented Apr 13, 2015

Coverage Status

Coverage increased (+1.32%) to 46.65% when pulling 295b04c on adonath:ts_map_performance_improvement into d89ced4 on gammapy:master.

@coveralls

This comment has been minimized.

coveralls commented Apr 13, 2015

Coverage Status

Coverage increased (+1.32%) to 46.65% when pulling c4f90fb on adonath:ts_map_performance_improvement into d89ced4 on gammapy:master.

cdeil added a commit that referenced this pull request Apr 13, 2015

Merge pull request #252 from adonath/ts_map_performance_improvement
Implement TS map computation in Cython

@cdeil cdeil merged commit 4d063c0 into gammapy:master Apr 13, 2015

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
@cdeil

This comment has been minimized.

Member

cdeil commented Apr 13, 2015

Thanks!

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