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

Add acceptance curve smoothing method #510

Merged
merged 5 commits into from Apr 17, 2016

Conversation

Projects
None yet
3 participants
@JouvinLea
Contributor

JouvinLea commented Apr 13, 2016

I added the smoothing for the bkg model

Lea Jouvin
@@ -8,6 +8,7 @@
from astropy.modeling.models import Gaussian1D
from astropy.table import Table
from astropy.units import Quantity
from scipy.ndimage.filters import convolve

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

Scipy is an optional dependency ... this import must be moved inside the function or method where it's used.

Also: use from scipy.ndimage import convolve instead of from scipy.ndimage.filters import convolve.

"""Smooth the bkg rate with a gaussian 1D kernel.
"""
for binE, E in enumerate(self.counts.energy[:-1]):

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

I think you're not using E, right?
Then you could just do this, no?

for idx_energy in range(len(self.counts.energy) - 1):
for binE, E in enumerate(self.counts.energy[:-1]):
counts = self.counts.data[binE, :]
Nev = np.sum(counts).value
if (Nev > 0):

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

Is this a speed optimisation to skip if Nev == 0?
Or does the code error if you execute a bin with Nev == 0?
Maybe you can leave a one-line comment like this (if this is accurate, I'm not sure):

# For zero counts, the background rate is zero and smoothing would not change it.
# For speed we're skipping the smoothing in that case
if str(ngroup) in self.models3D.keys():
self.models3D[str(ngroup)].smooth()
else:
log.info("No run in the band {}".format(ngroup))

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

"band" -> "group"
(we can discuss re-naming in the future, but for now let's be consistent and use "group" everywhere.

if str(ngroup) in self.models2D.keys():
self.models2D[str(ngroup)].smooth()
else:
log.info("No run in the band {}".format(ngroup))

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

Same here: "band" -> "group"

counts = self.counts.data[binE, :]
Nev = np.sum(counts).value
if (Nev > 0):
Np = len(counts)

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

I think the code would become easier to understand if you factor the part starting with Np = len(counts) down to the convolve call into a separate function _poisson_gauss_smooth (or something).
At the moment this method is quite complex, because you loop energy bands and then have an algorithm for each energy band and then put the result back in the 2D array. It would be better to factor the 1D algorithm out into a utility function.

So something like this (pseudo-code);

for idx in range(len(self.counts.energy)):
    counts = self.counts.data[idx]
    bkg = self.bkg.data[idx]
    bkg_smooth = _poisson_gauss_smooth(counts, bkg)
    self.bkg.data[idx] = bkg_smooth
modeltype : {'3D', '2D'}
Type of the background modelisation
Returns

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

This doesn't return anything (and shouldn't, right?).
So remove this from the docstring?

Type of the background modelisation
ngroup : int
Groups ID
Returns

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

Should this return something?
If yes, add it.
If no, remove Returns from docstring.

self.smooth_model(modeltype, ngroup)
def smooth_model(self, modeltype, ngroup):
"""Smooth the nkg model for one group

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

Typo: nkg -> bkg

@@ -198,17 +198,17 @@ def filename(self, modeltype, group_id, smooth=False):
else:
return 'background_{}_group_{:03d}_table.fits.gz'.format(modeltype, group_id)
def save_model(self, modeltype, ngroup):
def save_model(self, modeltype, ngroup, smooth=False):

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

I'm not sure it's a good idea to change the filename depending on smooth / no-smooth.

This class will get feature-richer, there will be different ways to smooth the model.
Probably filename management should be left up to the user in those cases?
(but leave as-is if you like for now)

def smooth(self):
"""Smooth the bkg rate with a gaussian 1D kernel.

This comment has been minimized.

@cdeil

cdeil Apr 13, 2016

Member

Maybe add some more info to the docstring. I don't mean a full description of the algorithm, but something like:

Calling this method modifies the ``bg_rate`` data member, replacing it with a smoothed version.

and

This method uses an adaptive Poisson method to compute the smoothing Kernel width
from the available counts (see code and inline comments for details).

@cdeil cdeil added the feature label Apr 13, 2016

@cdeil cdeil added this to the 0.4 milestone Apr 13, 2016

@cdeil cdeil self-assigned this Apr 13, 2016

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 13, 2016

@JouvinLea - I've done a review. Nothing major, I hope the suggestion to factor out the 1D smoothing algorithm in a separate function was clear).
Once those comments are addressed, this is ready to merge.

@cdeil cdeil changed the title from add smooth to Add acceptance curve smoothing method Apr 13, 2016

# in the counts histogram
Npix_sigma = np.minimum(Npix_sigma, Np / 6)
# kernel gaussian define between -3 and 3 sigma
x = np.linspace(-3, 3, 6 * Npix_sigma)

This comment has been minimized.

@adonath

adonath Apr 14, 2016

Member

I'd suggest to use the Gaussian1DKernel from astropy.convolution, this saves some code.

This comment has been minimized.

@adonath

adonath Apr 14, 2016

Member

When you care about speed, by far the best option would be to use scipy.ndimage.filters.gaussian_filter.

EDIT: Just checked, in this case with (small) 1D arrays, it doesn't make a big difference in performance, but still saves some code...

EDIT#2: There's also a scipy.ndimage.filters.gaussian_filter1d, which is slightly faster. Maybe use that?

if (Nev > 0):
Np = len(counts)
# Number of pixels per sigma of the kernel gaussian to have more than 150 events/sigma
Npix_sigma = (150 / Nev) * Np

This comment has been minimized.

@adonath

adonath Apr 14, 2016

Member

Minor: This and the following code could be reduced by using e.g. np.where() or np.clip()

@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 17, 2016

@cdeil
I think I took all you comments in to account, what do you think now?

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 17, 2016

Merging now.

Thank you!

@cdeil cdeil merged commit 30816c5 into gammapy:master Apr 17, 2016

0 of 2 checks passed

continuous-integration/appveyor/pr Waiting for AppVeyor build to complete
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment