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 mean PSF computation #549

Merged
merged 28 commits into from Jun 13, 2016

Conversation

Projects
None yet
2 participants
@JouvinLea
Contributor

JouvinLea commented Jun 4, 2016

I implemented this in gammapy one month ago but I wanted to test it on some picture before to present my idea!!This is computation of the mean of the psf for a set of runs in order to be able to do some morphology fit!
I think there is a lot to change to make it easier to use it but we will discuss this next week.

See you on Monday,
Lea

@cdeil cdeil added the feature label Jun 6, 2016

@cdeil cdeil added this to the 0.5 milestone Jun 6, 2016

@cdeil cdeil self-assigned this Jun 6, 2016

@cdeil cdeil modified the milestones: PyGamma1, 0.5 Jun 6, 2016

@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Jun 7, 2016

@cdeil
I changed a little the computation of the mean psf I think this is a lot more nicer now and it uses the EnergyDependantTablePSF and PSFTable!
See you tomorow

@@ -211,6 +224,44 @@ def excess_image(self):
total_excess.data = self.maps["counts"].data - self.maps["bkg"].data
self.maps["excess"] = total_excess
def make_mean_psf(self, region_center, spectral_index=2.3):

This comment has been minimized.

@cdeil

cdeil Jun 8, 2016

Member

Changes:

  • make_mean_psf -> make_psf
  • region_center -> pos
  • spectral_index can be removed
  • Update docstring
  • Remove if-else. Either don't do any check, or raise an error. Normal code (what's now in the else block) shouldn't be indented.

@cdeil cdeil changed the title from Image pipe meanpsf to Add mean PSF computation Jun 8, 2016

@@ -325,3 +384,35 @@ def excess_image(self):
total_excess = SkyMap.empty_like(self.empty_image)
total_excess.data = self.maps["counts"].data - self.maps["bkg"].data
self.maps["excess"] = total_excess
def make_mean_psf_tab(self, region_center, spectral_index=2.3):

This comment has been minimized.

@cdeil

cdeil Jun 8, 2016

Member

This is the example we wrote this morning how I would have implemented the mean PSF computation:
https://gist.github.com/cdeil/8f61fc0cb0604c9115942da283da9552

As we said, maybe this methods should be moved to a DataStoreObservation method.
And the mean PSF computation across observations should be separate from the integration of the PSF for a given energy band (that's the last step, a spectral model is only needed there).

Can you add an example or test to this PR, so that I can try out this code.

@JouvinLea JouvinLea force-pushed the JouvinLea:ImagePipe_meanpsf branch from a4d035d to 3952261 Jun 9, 2016

@requires_data('gammapy-extra')
def test_make_psftable():

This comment has been minimized.

@JouvinLea

JouvinLea Jun 9, 2016

Contributor

@cdeil
I don't have great ideas for testing this mean psf

This comment has been minimized.

@cdeil

cdeil Jun 10, 2016

Member

Please put this test for the mean psf:

psf1 = PSF from run 1
psf2 = PSF from run 2
psf_mean = mean PSF for both runs

# Add assertions that all three PSFs are similar, e.g.
# - compute the R68 and R95 in the 1 - 10 TeV energy band and see what you get.
# If the values are all similar and sensible, just add asserts on exactly those values (as a regression test).

# Maybe also assert something about the exposure attribute of the mean PSF ... shouldn't it be exactly or approximately the sum of the exposures of the two PSFs?

That's enough of a test for the mean PSF for now.

This comment has been minimized.

@JouvinLea

JouvinLea Jun 10, 2016

Contributor

@cdeil
hi,
I took two observations of the Crab in gammapy extra: 23523 and 23526 for testinf the mean psf computation.
The problem is that for these two observations the psfvalue you obtain for the EnergydependantTablePsf (psf1 and psf2) are really different so this is weird knowing that the Crab offset in the FOV for these two observations is almost equal (0.49 and 0.51 degree…). I think that there is either a problem on the interpolation of the parameters of the tripple gauss for these runs in your exporter or in the definition of the tripple gauss in gammapy. The difference between the two psf is superior to 10 order of magnitude so there is obviously a bug...

When I am doing the same test for the two same observations using the PA data and their king profile for the psf, it works fine. The mean psf value integrated between 1 and 10 TeV of the two runs is well between the values of each run and the R68 of the psf of each observation and of the psfmean is almost equal.

So I think what I coded is ok but for the moment I can’t do any test with the data of HAP in gammapy-extra …

This comment has been minimized.

@cdeil

cdeil Jun 10, 2016

Member

Can you please still put this test?
If there's an issue with our PSF, we'll fix it, we shouldn't ignore it.
You can add pytest.mark.xfail for now if it fails for an unrelated reason, so that this pull request can go in.

There's a few other inline comments from me that aren't addressed though, so this isn't ready to be merged yet.

This comment has been minimized.

@JouvinLea

JouvinLea Jun 11, 2016

Contributor

@cdeil
ok but I mean shouldn't we try to undersand why the HAP gaussian psf are so different for these two runs?

This comment has been minimized.

@cdeil

cdeil Jun 11, 2016

Member

ok but I mean shouldn't we try to undersand why the HAP gaussian psf are so different for these two runs?

Yes, we should!
(but I could do that and that doesn't have to hold up this pull request)

Do you have an example code snippet I can run to start investigating (that e.g. prints some numbers or makes a plot, e.g. in the examples folder in this branch, or in a separate https://gist.github.com/) or should I write my own?

This comment has been minimized.

@JouvinLea

JouvinLea Jun 11, 2016

Contributor

ok!!
I don't really have an example. But if you take the code in test_observationlist.py and try to access for examples psf1.psf_value or psf2.psf_value, you will see that there are completely different!!

sig1=Angle(np.array([ 0.02858012,  0.1211779 ,  0.0515614 ]), 'deg')
norm1=np.array([ 0.1735659 ,  0.41401832,  0.41241578])
 sig2=Angle(np.array([ 0.02559773,  0.10916024,  0.04580572]), 'deg')
 norm2=np.array([ 0.17220046,  0.40929013,  0.41850941])
filename = 'fov_bg_maps_test.fits'
mosaic_images.make_images(make_background_image=True, for_integral_flux=True, radius=10., make_psf = True,
region_center=center)
import IPython; IPython.embed()

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

Remove import IPython; IPython.embed()

Parameters
----------
source_position : `~astropy.coordinates.SkyCoord`
Coordinates of the interest region for which we want to calculate the psf

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

Suggest: Source position as description.

@@ -584,9 +585,94 @@ def make_exposure_image(self, fov, energy_range):
"""
raise NotImplementedError
def make_psf(self, source_position, energy=None, theta=None):
"""Evaluating the energydependent psf for a given energy array or theta array

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

Suggest: Make energy-dependent PSF for a given source position. as description.

----------
source_position : `~astropy.coordinates.SkyCoord`
Coordinates of the interest region for which we want to calculate the psf
energy: `~astropy.units.Quantity`

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

Missing space:

energy : `~astropy.units.Quantity`
arf = self.aeff.evaluate(offset=offset, energy=energy).to("cm2")
exposure = arf
exposure *= self.observation_live_time_duration
energy_dependant_psftab = EnergyDependentTablePSF(energy=energy, offset=theta, exposure=exposure,

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

Suggest other names:

  • psf_tab -> psf_value
  • energy_dependant_psftab -> psf
psf_value=psf_tab_tot.T)
return psf_total
def make_psftable(self, source_position, energy, theta=None, spectral_index=2.3):

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

This is OK, no?

obs.psf.to_table_psf(energy_band)

Then you don't need this method.
Remove?

@requires_data('gammapy-extra')
def test_makepsf():
"""Test creating a datastore as subset of another datastore"""

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

Remove docstring ... old copy & paste remnant?

class ObservationList(list):
"""List of `~gammapy.data.DataStoreObservation`
Could be extended to hold a more generic class of observations
"""
def make_energydependant_psf(self, source_position, energy=None, theta=None):

This comment has been minimized.

@cdeil

cdeil Jun 9, 2016

Member

Rename make_energydependant_psf -> make_psf

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 9, 2016

This is the test I would have written:

First test:

psf_mean = ObservationList.make_psf(two runs)
psf1 = Observation.make_psf(run number 1)
psf2 = Observation.make_psf(run number 2)
# Assert value at some offset / energy and some containment radius
# Should all be similar
def make_images(self, make_background_image=False, bkg_norm=True, spectral_index=2.3, for_integral_flux=False,
radius=10):
def make_images(self, make_background_image=False, bkg_norm=True, spectral_index=2.3,
radius=10, make_psf=False, region_center=None):

This comment has been minimized.

@cdeil

cdeil Jun 10, 2016

Member

I think the make_psf and region_center arguments should now be removed from this make_images method (and the docstring as well), because they aren't used any more?

Returns
-------

This comment has been minimized.

@cdeil

cdeil Jun 10, 2016

Member

Document what this make_1d_exposure method returns.

Is it really an exposure? Exposure computation doesn't need an assumed spectrum, right?
So maybe it's exposure (= area x time) times flux, i.e. expected counts?
Maybe pick a better name?

@cdeil cdeil modified the milestones: PyGamma1, 0.5 Jun 11, 2016

@JouvinLea JouvinLea force-pushed the JouvinLea:ImagePipe_meanpsf branch from 47e4c4c to 1584fe6 Jun 11, 2016

from ...data import DataStore, ObservationList
@pytest.mark.xfail

This comment has been minimized.

@cdeil

cdeil Jun 11, 2016

Member

I'm getting this error:

gammapy/data/tests/test_observationlist.py:11: in <module>
    @pytest.mark.xfail
E   NameError: name 'pytest' is not defined

You have to put this import at the top of the file:

from astropy.tests.helper import pytest
obslist = ObservationList(list)
psf_tot = obslist.make_psf(source_position=center, energy=energy)
psf_tot_int = psf_tot.table_psf_in_energy_band(energy_band, spectral_index=2.3)
assert_allclose(psf_tot_int.containment_radius(0.68), psf1_int.containment_radius(0.68), rtol=1e-3)

This comment has been minimized.

@cdeil

cdeil Jun 11, 2016

Member

After fixing the import, I get this test fail:

______________________________________________________________ test_make_psftable _______________________________________________________________

    @requires_data('gammapy-extra')
    def test_make_psftable():
        center = SkyCoord(83.63, 22.01, unit='deg')
        store = gammapy_extra.filename("datasets/hess-crab4-hd-hap-prod2")
        data_store = DataStore.from_dir(store)
        obs1 = data_store.obs(23523)
        obs2 = data_store.obs(23526)
        energy = EnergyBounds.equal_log_spacing(1, 10, 100, "TeV")
        energy_band = Energy([energy[0].value, energy[-1].value], energy.unit)
        psf1 = obs1.make_psf(source_position=center, energy=energy, theta=None)
        psf2 = obs2.make_psf(source_position=center, energy=energy, theta=None)
        psf1_int = psf1.table_psf_in_energy_band(energy_band, spectral_index=2.3)
        psf2_int = psf2.table_psf_in_energy_band(energy_band, spectral_index=2.3)
        list = [obs1, obs2]
        obslist = ObservationList(list)
        psf_tot = obslist.make_psf(source_position=center, energy=energy)
        psf_tot_int = psf_tot.table_psf_in_energy_band(energy_band, spectral_index=2.3)
>       assert_allclose(psf_tot_int.containment_radius(0.68), psf1_int.containment_radius(0.68), rtol=1e-3)
E       AssertionError: 
E       Not equal to tolerance rtol=0.001, atol=0
E       
E       (mismatch 100.0%)
E        x: Angle(<Angle 0.12588112016412742 deg>)
E        y: Angle(<Angle 0.11645282474146176 deg>)

gammapy/data/tests/test_observationlist.py:29: AssertionError

My recommendation would be that you change the test to this (with the imports at the top of the file):

    from astropy.coordinates import Angle
    from astropy.tests.helper import assert_quantity_allclose
    # Check that the mean PSF is consistent with the individual PSFs
    # (in this case the R68 of the mean PSF is in between the R68 of the individual PSFs)
    assert_quantity_allclose(psf1_int.containment_radius(0.68), Angle(0.11645282474146176, 'deg'))
    assert_quantity_allclose(psf2_int.containment_radius(0.68), Angle(0.14232666985400985, 'deg'))
    assert_quantity_allclose(psf_tot_int.containment_radius(0.68), Angle(0.12588112016412742, 'deg'))

I.e. put assert_quantity_allclose statements against fixed values as regression tests, to establish the current results (so that we note when something changes or the method is improved in the future).

OK?

This comment has been minimized.

@cdeil

cdeil Jun 11, 2016

Member

The R68 asserts I posted pass on my machine.
Maybe remove the xfail from this test?
(and if the tests further down fail, uncomment those with a coment # TODO: check why these numbers don't make sense?

Then at least the mean PSF computation is executed and tested.

Lea Jouvin
@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Jun 12, 2016

@cdeil
The passed Travis, so can we merge?
Cheers Lea

Lea Jouvin
if not theta:
theta = self.psf.to_table_psf(theta=offset).offset
psf_value = self.psf.to_table_psf(theta=offset).evaluate(energy)
arf = self.aeff.evaluate(offset=offset, energy=energy).to("cm2")

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Remove to('cm^2') here?
Seems like a weird place to enforce units, no?

If you do want to enforce units, put an extra line

exposure = exposure.to('cm^2 s')

below ... the units of observation_live_time_duration could change as well (e.g. hour or seconds).
But I don't think we need to enforce units here at all, do we?

Returns
-------
energy_dependant_psftab : `~gammapy.irf.EnergyDependentTablePSF`

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Rename energy_dependant_psftab -> psf.
(that name just duplicates the info that is given by the type EnergyDependentTablePSF on the same line)

Returns
-------
energy_dependant_psftab : `~gammapy.irf.EnergyDependentTablePSF`

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Rename energy_dependant_psftab -> psf.
(that name just duplicates the info that is given by the type EnergyDependentTablePSF on the same line)

psf2 = obs2.make_psf(source_position=center, energy=energy, theta=None)
psf1_int = psf1.table_psf_in_energy_band(energy_band, spectral_index=2.3)
psf2_int = psf2.table_psf_in_energy_band(energy_band, spectral_index=2.3)
list = [obs1, obs2]

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Replace this:

list = [obs1, obs2]
obslist = ObservationList(list)

with this:

obslist = ObservationList([obs1, obs2])

(list is a Python built-in type, and you shouldn't overwrite it by assigning [obs1, obs2] to it ... that can cause very weird issues further down and it's also something PyCharm or other static code analysers will warn you about not to do.)

Returns
-------
offset_tab: `numpy.array`

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

I think it would be nicer if you return a Table with columns offset and npred instead of two arrays separately.

eref = EnergyBounds(self.energy_band).log_centers
spectrum = (energy_bin / eref) ** (-spectral_index)
offset_tab = Angle(np.linspace(self.offset_band[0].value, self.offset_band[1].value, 10), self.offset_band.unit)
arf = self.aeff.evaluate(offset=offset_tab, energy=energy_bin).to("cm2").T

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Remove .to("cm2") here.
Again, if you do want to force units, do it on a separate line at the end of the function for the quantity you return.
But I don't think we should be enforcing units in Gammapy in computations, just on I/O.

This comment has been minimized.

@JouvinLea

JouvinLea Jun 13, 2016

Contributor

@cdeil
I agree with you!
Do you think that I should enforce the unity in the method make_imagesof the class MosaicImages when I am saving the exposure map as a skymap collection? But maybe it's better to leave it as m^2? It's just that you don't have unit in a skymap so when I read the skymap collection I wanted it to be in cm^2 in order to have flux in cm-2?

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

I haven't used SkyMapCollection with units yet, so I don't know yet how we should handle this.
I would start by not enforcing units, and only when it becomes clear that enforcing units is helpful add it in some places.

This comment has been minimized.

@JouvinLea

JouvinLea Jun 13, 2016

Contributor

@cdeil
The problem here for example is that you don't know if the exposure map is in m2 ou cm2 and this is problematic when you want to compute the flux?

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Why? If we use quantities everywhere and propagate units, the flux will come out with some flux units.
It could be m^-2 aor 'cm^-2' if you don't pin the output units.

But we should stop discussion this here ... that's a question we have to decide for 100s or functions in Gammapy and should discuss in a separate issue and then handle in a uniform way.

I just would prefer with starting more or less with a codebase in Gammapy where we always either:

  1. don't pin units, output units depend on input units, or
  2. pin units at the very end of the function, on the quantity that's returned.

If you do one of those two coding styles here it's fine with me. OK?

This comment has been minimized.

@JouvinLea

JouvinLea Jun 13, 2016

Contributor

@cdeil
Ok let's stop this discussion in this PR. But just my point is that for the moment you can't attach unit to skymap because it's failing when writting to fits with the Value error. Therefore when the user will read again the exposure map, he will not know the unit... so the fact that the exposure is in m2 or cm2 is absolutely not a pb, the problem is that the skymap doesn't have the information when you write it....

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

@JouvinLea - I've split out the skymap units problem into a separate Github issue:
#575
You're right that this isn't working correctly yet.

@cdeil

This comment has been minimized.

Member

cdeil commented Jun 13, 2016

@JouvinLea - Thanks for the update. Tests now passed on travis-ci. I've left some more inline comments.

My main suggestion here would be to add one more test that call the make_psf for a single observation, to make sure that the handling of user-supplied energy and offset as well as the default values work (it's not clear to me from reading the code that this is the case).
I.e. ideally compute the psf 4 times with:

  • energy=None, theta=None
  • energy=something, theta=None
  • energy=None, theta=something
  • energy=something, theta=something
    and assert on the shape and one pixel value of the returned psf data members (i.e. the energy, theta, psf_value and exposure there).

This is a good example how you can parametrize a test, so that you don't have to write a lot of boilerplate code:

@pytest.mark.parametrize("pars,results",[

@JouvinLea - If you don't want to add this, let me know, I can also do it.

Lea Jouvin added some commits Jun 13, 2016

Lea Jouvin
Lea Jouvin
@@ -301,6 +319,8 @@ def make_images(self, make_background_image=False, bkg_norm=True, spectral_index
if make_background_image:
self.maps["bkg"] = total_bkg
self.maps["exposure"] = total_exposure
self.maps["exposure"].unit = obs_image.maps["exposure"].unit

This comment has been minimized.

@JouvinLea

JouvinLea Jun 13, 2016

Contributor

@cdeil:
If you add the unit to the Skymap exposure, you get this error when you save it as a fits with the write method:

ValueError: Illegal value: Unit("cm2 s TeV")

I don't know why?

data_store = data_manager['hess-crab4-hd-hap-prod2']
t = data_store.data_summary([23523, 23592])
assert t[0]['events'] == 620975
assert t[1]['edisp_2d'] == 28931
t = data_store.data_summary([23523, 23592], summed=True)
assert t[0]['psf_3gauss'] == 6042
@pytest.mark.parametrize("pars,result", [(dict(energy=None, theta=None),

This comment has been minimized.

@JouvinLea

JouvinLea Jun 13, 2016

Contributor

@cdeil:
I add this test for a single obs, is this what you wanted?

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Yes, this is exactly what I had in mind.
Thank you!

Can you try to make the content in parametrize a little more readable (i.e. not so far on the right side)?
Maybe just adding a newline after the opening [ would be helpful, i.e. something like this:

@pytest.mark.parametrize("pars,result", [
    "put stuff here"
])
@@ -301,6 +319,8 @@ def make_images(self, make_background_image=False, bkg_norm=True, spectral_index
if make_background_image:
self.maps["bkg"] = total_bkg
self.maps["exposure"] = total_exposure
#self.maps["exposure"].unit = obs_image.maps["exposure"].unit

This comment has been minimized.

@JouvinLea

JouvinLea Jun 13, 2016

Contributor

@cdeil, @adonath
There is a error here if I give a unit to the exposure map when you write the SkyMap:

ValueError: Illegal value: Unit("cm2 s TeV").
Returns
-------
table: `astropy.table.Table`

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Missing whitespace before the colon, should be:

table : `astropy.table.Table`

Please always check the HTML version of new docstrings you've written locally using:

python setup.py build_sphinx
open docs/_build/html/index.html
Lea Jouvin
data_store = DataStore.from_dir(store)
obs1 = data_store.obs(23523)
psf = obs1.make_psf(source_position=center, energy=pars["energy"], theta=pars["theta"])
assert_allclose(len(psf.offset), result["theta_shape"])

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Can you please change all these len function calls to .shape attribute calls. I.e. instead of len(psf.offset) call psf.offset.shape and check that.

This is a better check, because len of a multi-dimensional array will return the length of the first (or last, not sure) dimension, and checking the shape explicitly tests more.

This comment has been minimized.

@JouvinLea

JouvinLea Jun 13, 2016

Contributor

@cdeil :
ok this is what I did for the psf_value that is a multidimensional array but offset is 1D as energy or exposure so why is it important?

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Just to make sure they are 1d arrays.
(I didn't have much time this morning when reading your code, but somehow it wasn't obvious to me that it works in all cases and returns correct PSFs, with energy and offset 1D and psf_value 2D.)
Can't hurt to assert on shape here, no?

Lea Jouvin
if for_integral_flux:
norm = np.sum(spectrum * energy_band)
npred_tab /= norm
table = Table()

This comment has been minimized.

@cdeil

cdeil Jun 13, 2016

Member

Can you change to this?

from astropy.table import QTable # put this at the top of the file
table = QTable()
table['theta'] = offset_tab
table['npred'] = npred_tab

QTable is a table that can store quantities as columns.
And there's no need to create Column objects first, they will be created on table['colname'] = ....

Lea Jouvin
@cdeil

This comment has been minimized.

Member

cdeil commented Jun 13, 2016

Test fails on travis-ci are unrelated.

This looks good to me ... merging this now.

Thanks!

@cdeil cdeil merged commit 3df2a01 into gammapy:master Jun 13, 2016

1 of 2 checks passed

continuous-integration/travis-ci/pr The Travis CI build failed
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
@cdeil

This comment has been minimized.

Member

cdeil commented Jun 13, 2016

I've added a changelog entry and done some minor code cleanup in 88d01dd .

One thing I changed was to rename "source_position" to "position" ... sometimes one wants to compute the PSF for any kind of position, there doesn't have to be a source present, and just "position" is shorter.

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