Skip to content
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

Simplify PSF stacking and containment radius computation #2019

Merged
merged 12 commits into from Feb 2, 2019

Conversation

@adonath
Copy link
Member

@adonath adonath commented Jan 31, 2019

This PR simplifies the PSF stacking and containment radii computation code for EnergyDependdentTablePSF. It also fixes #2008.

Copy link
Contributor

@registerrier registerrier left a comment

Thanks @adonath ! Some minor comments.

Loading

@@ -598,31 +603,34 @@ def info(self):
"""Print basic info"""
print(str(self))

def plot_psf_vs_rad(self, energies=[1e4, 1e5, 1e6]):
def plot_psf_vs_rad(self, energies=None, ax=None, **kwargs):
Copy link
Contributor

@registerrier registerrier Feb 1, 2019

Choose a reason for hiding this comment

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

Is energy=None is really safe/meaningful? Isn't it better not to put any default?

Loading

Copy link
Member Author

@adonath adonath Feb 2, 2019

Choose a reason for hiding this comment

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

Mutable objects such as dicts, lists or "Quantity" should be avoided as default arguments, because they are created when the function definition is executed and so are shared between multiple instances. The patter we adopted in Gammapy is to put None instead and put the default in the first line of the function if energies is None: .... I think it's nice to have a default, but I don't insist to keep it.

Loading

psf_value = self.evaluate(energy=energies, rad=rad).to_value("deg-2")
containment = cumtrapz(y=rad * psf_value, x=rad.to_value("deg"), axis=1)

with np.errstate(divide="ignore", invalid="ignore"):
Copy link
Contributor

@registerrier registerrier Feb 1, 2019

Choose a reason for hiding this comment

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

Should the renormalization be the default behavior? Or should it be done here, instead of __init__ for instance?

Loading

Copy link
Member Author

@adonath adonath Feb 2, 2019

Choose a reason for hiding this comment

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

I think this is basically an open question. I guess ideally the PSF would be normalized already on the DL3 data level and only re-normalized in Gammapy, when the normalization changed because of some computation (such as broadening or clipping the PSF). Right now PSFs are not necessarily normalized and we need to do it here to get the correct containment.

Loading

from astropy.io import fits
from astropy.units import Quantity
from astropy import units as u
Copy link
Contributor

@registerrier registerrier Feb 1, 2019

Choose a reason for hiding this comment

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

Is using u.Quantity instead of importing Quantity preferable for some reason?

Loading

Copy link
Member Author

@adonath adonath Feb 2, 2019

Choose a reason for hiding this comment

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

from astropy import units as u is the recommended import pattern for astropy.units, which we already use in most parts of Gammapy. To avoid the second import I typically use u.Quantity(). It doesn't really matter, but slightly prefer, when it's consistent.

Loading

psf.containment_radius(1 * u.TeV, 0 * u.deg, 0.5),
rtol=1e-3,
psf.containment_radius(1 * u.TeV, theta=0 * u.deg, fraction=0.5),
rtol=1e-2,
Copy link
Contributor

@registerrier registerrier Feb 1, 2019

Choose a reason for hiding this comment

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

This will require a code rebase.

Loading

Copy link
Member Author

@adonath adonath Feb 2, 2019

Choose a reason for hiding this comment

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

Done

Loading

def _get_1d_psf_values(self, energy_index):
"""Get 1-dim PSF value array.
# TODO: make this public once the API for IRF stacking is defined
def _stack(self, psf):
Copy link
Contributor

@registerrier registerrier Feb 1, 2019

Choose a reason for hiding this comment

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

Why make it private?

Loading

Copy link
Member Author

@adonath adonath Feb 2, 2019

Choose a reason for hiding this comment

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

Done

Loading

psfs = [self.table_psf_at_energy(energy, interp_kwargs) for energy in energies]
rad = [psf.containment_radius(fraction) for psf in psfs]
return Quantity(rad)
# oversample for better precision
Copy link
Contributor

@registerrier registerrier Feb 1, 2019

Choose a reason for hiding this comment

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

Shouldn't the containment_radius of the TablePSF be adapted too, to follow the same approach?
This would completely remove the need for spline parametrization.

Loading

Copy link
Member Author

@adonath adonath Feb 2, 2019

Choose a reason for hiding this comment

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

I agree this is a good suggestion, but I'd prefer to do this in a separate PR.

Loading

@adonath adonath self-assigned this Feb 2, 2019
@adonath adonath added the cleanup label Feb 2, 2019
@adonath adonath added this to the 0.11 milestone Feb 2, 2019
@adonath adonath merged commit 38bc7fa into gammapy:master Feb 2, 2019
4 checks passed
Loading
@adonath adonath mentioned this pull request Feb 5, 2019
@adonath adonath changed the title Simplify PSF stacking and containment radii code Simplify PSF stacking and containment radius computation Mar 26, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

2 participants