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

Improve npred cube functionality #151

Merged
merged 7 commits into from Aug 20, 2014

Conversation

Projects
None yet
2 participants
@ellisowen
Contributor

ellisowen commented Jul 24, 2014

Implementation which I intend to use to produce significance images (currently this is still running on my machine).

The unit tests will need much improvement soon as I just implemented some quick ones. Looking at the outputs with ds9 using the FermiGalacticCenter files seemed reasonable, but obviously not rigorous yet.

Don't know if you were aiming to merge this asap, so did a pull request just in case.

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 24, 2014

No need to merge this ASAP.
I won't have time to trying it out today, so just continue on your machine for now (and add tests that show that this is doing what it's supposed to do) and I'll review this on Monday.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jul 25, 2014

To do (Ellis):

  • Bug fix
  • Include documents of tests against gtmodel for Vela region

@ellisowen ellisowen referenced this pull request Jul 26, 2014

Merged

Add Fermi Vela dataset #156

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jul 29, 2014

Pushed for discussion & review... I think this will need PR #156 operational to work

@@ -361,15 +432,22 @@ def __repr__(self):
def compute_npred_cube(flux_cube, exposure_cube, energy_bounds):
"""TODO
"""Computes predicated counts cube in energy bins.

This comment has been minimized.

@cdeil

cdeil Jul 30, 2014

Member

predicated -> predicted

False: returns error if a file exists of the same name in the
output directory.
Returns

This comment has been minimized.

@cdeil

cdeil Jul 30, 2014

Member

No need to add a Returns section for functions / methods that don't return anything.

@cdeil

This comment has been minimized.

Member

cdeil commented Jul 30, 2014

I can review this tonight.

There's a few TODO left in the code that should be resolved before this can be merged.

And for your example scripts, please add an empty line every once in a while where a new step starts (e.g. compute A, compute B, plot C) ... I find it hard to get a quick overview of what this code does if it's so dense.

"""Runs commands to produce convolved predicted counts map in current directory
"""
from gammapy.spectral_cube import compute_npred_cube, GammaSpectralCube, convolve_npred_cube

This comment has been minimized.

@cdeil

cdeil Jul 30, 2014

Member

Sort imports.

First scientific Python packages like numpy and matplotlib.
Then astropy
Then gammapy.

ratio = np.nan_to_num(model / gtmodel)
# Plotting
vmin, vmax = 0, 1
fig, axes = plt.subplots(nrows=1, ncols=3)

This comment has been minimized.

@cdeil

cdeil Jul 30, 2014

Member

I think we should start using APLPy for sky image plotting to have sky coordinate labels.

The Sphinx build on travis-ci as well as the docs build on readthedocs are already configured such that this should work.

There might be some plotting applications that require using WCSAxes directly instead of using the higher-level APLPy.

This comment has been minimized.

@ellisowen

ellisowen Jul 31, 2014

Contributor

A number of changes would be required to get APLPy to work here & would take me a little time to get going. Can we keep using matplotlib here for now, and then switch to APLPy in a future pull request? Or should this be a priority?

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Jul 31, 2014

I think this addresses your comments.

The Travis CI failure is due to two time out errors, so I guess I don't need to worry about this?

@@ -140,32 +140,76 @@ def test_spatial_coordinate_images(self):
# TODO assert the four corner values
@pytest.mark.skipif('not HAS_SCIPY')

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

reproject is an optional dependency ... please add a skipif to any tests that use it.
Note that you can "stack" decorators, i.e. simply have two lines:

@pytest.mark.skipif('not HAS_SCIPY')
@pytest.mark.skipif('not HAS_REPROJECT')

This comment has been minimized.

@ellisowen

ellisowen Aug 6, 2014

Contributor

I changed this in the line below so it didn't spot the change in the diff here

# Check npredicted is same order of magnitude of true counts
assert_allclose(expected_sum, actual_sum, rtol=1)
# PSF convolve the npred cube
npred_cube_convolved = convolve_npred_cube(npred_cube, max_offset=3,

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

I see this test error locally:
https://gist.github.com/cdeil/ca385fe7271a6a9bfdb2

Looks like the signature is

def convolve_npred_cube(npred_cube, psf_object, max_offset, resolution=1):

i.e. you have to pass a psf object?

I don't know why this didn't show up on travis-ci ... can you rebase this PR on master to make sure that's not the issue?

from gammapy.irf import EnergyDependentTablePSF
# Reads in data
background_file = get_fermi_diffuse_background_model(filename='gll_iem_v05_rev1.fit')

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

gll_iem_v05_ref1.fit is a 500 MB file.

Can we avoid this in the Sphinx docs build?
I think travis-ci will fail always and often to download it, like it did here

Can you cut out a small part using ftcopy like I did here
and put the resulting file (which should be a few MB at most) here?

If it's still large (say more than 5 MB), you could just use a cutout of one of the older, lower-resolution diffuse model files?

This comment has been minimized.

@ellisowen

ellisowen Aug 5, 2014

Contributor

The cutout file is about half a MB, so will put it into the vela_region folder of the gammapy-extra repo in a new PR

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 4, 2014

Can you please move the tutorial in a sub-folder docs/tutorials/npred?
(we really should do this with all tutorials ... as this grows it will otherwise be hard to know which .rst or .py files belong together.

reprojected_cube : GammaSpectralCube
Cube spatially reprojected to the reference cube.
"""
from reproject.interpolation import reproject_celestial

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

The reproject.interpolation.reproject_celestial function will be renamed and I think should be considered an implementation detail.

Can you use the top-level reproject function instead?

new_cube = np.zeros((cube.shape[0], reference.shape[1],
reference.shape[2]))
energy_slices = np.arange(cube.shape[0])
for i in energy_slices:

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

reproject should work for cubes.
So do you really need to loop over energy slices manually here?
Please try to get this to work, but if you can't, leave a TODO comment here so that we don't forget to make this work in the future.

Parameters
----------
reference_cube : GammaSpectralCube

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

Always put single ticks around class names so that a html link is generated in the html docs.

Returns
-------
reprojected_cube : GammaSpectralCube

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

same here ... add single ticks around class name.

Quantitatively, we may compare these against the true counts observed by Fermi LAT in this region for the same parameters:
* True counts: 261
* Fermi Tools gtmodel predicted counts: 352

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

If you have 261 true counts and 352 predicted counts from gtmodel, obviously the gtmodel npred image is incorrect.

I guess the reason is that you don't have enough energy bins in the npred computation?
Can't you increase those by a factor of a few to get a correct npred image?

This comment has been minimized.

@ellisowen

ellisowen Aug 4, 2014

Contributor

So, re-running gtmodel with 100 energy bins in the npred computation (instead of 3) gives 285 counts, which seems much more reasonable...

I guess this improved gtmodel result should also be included in gammapy-extra (?) as it would be nice to include a better gtmodel example in the tutorial... should I replace the existing gtmodel npred cube in gammapy-extra, or provide this as an extra one?

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

Replace the existing one. But maybe 10 or 20 energy bins in the get model computation should be enough to get the same result as with 100 say within 1% or so, no?

This comment has been minimized.

@ellisowen

ellisowen Aug 4, 2014

Contributor

Yeah, I went for overkill just to be sure. I will replace the files with ones for 20 energy bins (gives npred counts within 1% of the 100 energy bins npred counts)

@cdeil cdeil referenced this pull request Aug 4, 2014

Merged

Add code to make model images from a source catalog #160

5 of 5 tasks complete
repro_bg_cube = background_model.reproject_to(exposure_cube)
# Define energy band required for output
energies = Quantity([10000, 500000], 'MeV')

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

Use 'GeV' instead of 'MeV' ... Easier to read if fewer zeros.

npred_cube = compute_npred_cube(repro_bg_cube, exposure_cube, energies)
# Convolve with Energy-dependent Fermi LAT PSF
psf_object = EnergyDependentTablePSF.read(FermiGalacticCenter.filenames()['psf'])

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

It probably doesn't matter much, but why not add a PSF file to the Vela example to be sure it's correct?

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

'Psf_object' -> 'psf'

counts_wcs = WCS(fits.open(counts_file)[0].header)
counts_cube = GammaSpectralCube(data=Quantity(counts_data, ''), wcs=counts_wcs,
energy=energies)
counts_cube = counts_cube.reproject_to(npred_cube)

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

Reprojecting counts is always problematic because you get non-integer output and Poisson stats are incorrect ... It should not be necessary here, I.e. Counts and exposure should have the same spatial binning, no?

This comment has been minimized.

@ellisowen

ellisowen Aug 4, 2014

Contributor

It's not necessary here specifically; I did need it in a few early test cases where the spatial binning was not the same, and left it in to keep the example more general... I guess people can figure out themselves whether or not this is required in their case though, so this isn't needed any more.

energy=energies)
counts_cube = counts_cube.reproject_to(npred_cube)
counts = np.nan_to_num(counts_cube.data[0])

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

There shouldn't be nan values in counts or model ... Why did you add this?

This comment has been minimized.

@ellisowen

ellisowen Aug 4, 2014

Contributor

Again, as per previous comment. 'NaN' is returned by reproject for the 'outside' areas if a larger region is reprojected to a smaller region. I will remove this as it's probably not actually helpful.

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.3, 0.025, 0.4])
fig.colorbar(im, cax=cbar_ax)
a = fig.get_axes()[0]

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

I'd still prefer you use APLPy, but if you stick with matplotlib use a for loop or list comprehension here.

This comment has been minimized.

@ellisowen

ellisowen Aug 5, 2014

Contributor

I still can't get APLPy to work here without writing a fits file to disk, which probably isn't desirable.

The error I get is:

/home/ellis/.local/lib/python2.7/site-packages/APLpy-0.9.12.dev620-py2.7.egg/aplpy/wcs_util.pyc in wcs_pix2sky(self, x, y, origin)
157 return yw, xw
158 else:
--> 159 return AstropyWCS.wcs_pix2sky(self, x, y, origin)
160 else:
161 coords = []

AttributeError: type object 'WCS' has no attribute 'wcs_pix2sky'

In the meantime, I will stick with matplotlib with a loop or list of some kind, like you said

This comment has been minimized.

@cdeil

cdeil Aug 5, 2014

Member

OK to stick with matplotlib here ... I can try to change the plotting to APLPy or WCSAxes here once I'm able to run this code.

Checks
------
For small regions, the predicted counts cube and significance images may be checked against the gtmodel output. The Vela region shown above is taken in this example in one

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

Use ~80 char line length limit also in RST files.

.. plot:: tutorials/npred_convolved.py
Quantitatively, we may compare these against the true counts observed by Fermi LAT in this region for the same parameters:

This comment has been minimized.

@cdeil

cdeil Aug 4, 2014

Member

Remove 'qualitatively' ...this is a quantitative comparison.

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 5, 2014

@ellisowen - I've merged the other related PRs ... can you update this one tomorrow?

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 5, 2014

Yes, should be updated by first thing tomorrow

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 5, 2014

This might need some more review ... I couldn't try this out yet because my internet connection is super-bad and the 500 MB diffuse model file download always failed.
So if you have a working version where the examples in the Sphinx docs build work with the small file, please push to Github and let me know.

energy_bounds : `~astropy.units.Quantity`
An array of Quantities specifying the edges of the energy bins
required for the predicted counts cube.
integral_resolution : int (optional)

This comment has been minimized.

@cdeil

cdeil Aug 18, 2014

Member

Let's try to use consistent names ... maybe energy_band and energy_bins as in the GammaSpectralCube.integral_flux_image method?

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 18, 2014

It's great that you now have consistent npred between your code and gtmodel!

I've left a few minor inline comments, but this is basically ready once you fix the issue mentioned at gammapy/gammapy-extra#6 (comment) .

This will change some of the unit tests reference values, so please make sure to run python setup.py test --remote-data and python setup.py build_sphinx to make sure this is OK.

About the model images you show above ... can you change the color scale for the model images (not the ratio image)? It's all blue now and I can't see anything.
You might also want to use projection_type ='bicubic' ... that should give smoother and better results (and be a bit slower, but that doesn't matter here).

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 18, 2014

Some of the Travis CI failures appear to relate to stuff that isn't changed in this PR, the rest should pass when PR 6 in gammapy-extra (vela corrections) is merged.

Otherwise, is this OK?

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 19, 2014

Restarted travis-ci after the other pull request has been merged and the unrelated error has been fixed ... let's see.

@@ -188,7 +188,7 @@ class FermiVelaRegion(object):
* 10 GeV < Energy < 500 GeV in 3 equal log10 bins
* CenterPosition: GLON = 263.05836967702709; GLAT = 3.9298511274632784
* Coordinate System: GALACTIC
* Image Size: 5 x 5 degrees (Search center 3 degrees)
* Image Size: 5 x 5 degrees (Search center 4 degrees)

This comment has been minimized.

@cdeil

cdeil Aug 19, 2014

Member

I think the search radius is 5 deg actually.
And the GLAT given above is incorrect.

I'd suggest to just add a link to the README in the gammapy-extra repo to avoid duplicating this info and having to check / adapt it in multiple places.

TODO
npred_cube : GammaSpectralCube
Predicted counts cube in energy bins.
psf_object : EnergyDependentTablePSF

This comment has been minimized.

@cdeil

cdeil Aug 19, 2014

Member

Remove _object ... everything in Python is an object, so putting this in the name is not useful.

Please check the rendered html version for the docstrings you write, in this case you need this to generate an HTML link:

psf : `~gammapy.irf.EnergyDependentTablePSF`
return npred_cube
def convolve_npred_cube(npred_cube, max_offset, resolution=1):
"""TODO
def convolve_npred_cube(npred_cube, psf_object, max_offset, resolution=1):

This comment has been minimized.

@cdeil

cdeil Aug 19, 2014

Member

There's nothing specific about this function why it would only work for npred cubes.
We'll also apply it to real counts cubes or to estimated background cubes (not necessarily the npred from a flux model ... could be from the data).

Actually, why didn't we make this a class method, i.e.

psf = EnergyDependentTablePSF(...)
cube = Spectralcube(...)
cube.convolve(psf)

(if there is a reason, please change at least the function name convolve_npred_cube to convolve_cube and the first parameter name from npred_cube to cube.)

This comment has been minimized.

@ellisowen

ellisowen Aug 19, 2014

Contributor

I changed the name as you suggest.

We didn't make this a class method as you planned to make big changes to the spectral_cube API in the near future (i.e. have separate classes for Spectral cubes, counts cubes, exposure cubes, etc) and it isn't currently clear how best this will fit in.

def convolve_npred_cube(npred_cube, max_offset, resolution=1):
"""TODO
def convolve_npred_cube(npred_cube, psf_object, max_offset, resolution=1):
"""Convolves a predicted counts cube in energy bins with the Fermi PSF.

This comment has been minimized.

@cdeil

cdeil Aug 19, 2014

Member

"the Fermi PSF" -> "an energy-dependent PSF"
(we'll use the same PSF class for HESS or CTA)

return hdu_list
def writeto(self, filename, clobber=False):

This comment has been minimized.

@cdeil

cdeil Aug 19, 2014

Member

I still think this method is superfluous, because

cube.to_fits().writeto('cube.fits')

is almost as short as

cube.writeto('cube.fits')

But if you absolutely want to keep it, please implement it as the one line

self.to_fits.writeto(filename, clobber)

instead of re-writing the three lines from to_fits.

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 19, 2014

OK ... finally I was able to build the Sphinx docs for this ... I've put them here for reference:
http://www.mpi-hd.mpg.de/personalhomes/deil/public/gammapy/issue_151/tutorials/npred/predicted_counts.html

  • I thought now the predicted counts agreed completely ... now I see again gtmodel_npred = 265 and gammapy_npred = 281. Are these simply outdated numbers or are results non-consistent again?
  • I don't understand the content of any of the images shown. The background model should look like what ds9 gammapy-extra/datasets/vela_region/background_vela.fits shows like, no? (see attached screenshot). Using the same color scale for the npred images and ratio image is confusing ... I'd suggest plotting the ratio image separately with it's own colorbar.
  • The gtmodel image looks smoother ... does the PSF convolution work properly?
  • The significance images don't show Vela? In any case I'd suggest removing them from this tutorial or putting it at the end, because the npred computation is a separate first step and the computation of a significance image is a separate later step.
  • Don't repeat the dataset parameters (the bullet list) here ... this is hard to maintain ... instead link to this README file.

screen shot 2014-08-19 at 14 20 03

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 19, 2014

I know we've discussed this before and there was some problem, but for sky images in the docs we really should use aplpy.FITSFigure.
It should just work if you pass it a filename or an HDU.
Could you please give this another try and paste the exact code / error message here if you can't get this to work?

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 19, 2014

The sphinx docs look quite different to mine... a couple of questions in line below:
(The other stuff I will address and push up a new version this evening).


OK ... finally I was able to build the Sphinx docs for this ... I've put them here for reference:
http://www.mpi-hd.mpg.de/personalhomes/deil/public/gammapy/issue_151/tutorials/npred/predicted_counts.html

  • I thought now the predicted counts agreed completely ... now I see again gtmodel_npred = 265 and gammapy_npred = 281. Are these simply outdated numbers or are results non-consistent again?

These values have never agreed exactly - from what I understand, the method I implement differs from the method used by gtmodel (I don't understand the full details of how this works) - so we thought before values within 10% would be fine?

  • The gtmodel image looks smoother ... does the PSF convolution work properly?

My PSF parameters may be different - I looked for a few hours yesterday to try and find out what kernel size, correction radius etc I should use, but didn't find this info yet... do you know what would be best?

  • The significance images don't show Vela? In any case I'd suggest removing them from this tutorial or putting it at the end, because the npred computation is a separate first step and the computation of a significance image is a separate later step.

This was due to the URL being wrong in gammapy.datasets.load (I couldn't check this before while the gammapy-extra PR wasn't merged). This is now fixed and shows correct images.

  • I know we've discussed this before and there was some problem, but for sky images in the docs we really should use aplpy.FITSFigure.
    It should just work if you pass it a filename or an HDU.
    Could you please give this another try and paste the exact code / error message here if you can't get this to work?

The issue here is that aplpy doesn't work on my laptop (normally it just freezes and then throws some time out error or memory error after 10 mins or so). This is why we decided to make a second PR in the days when I'm back in HD & can work on my desktop, and stick with matplotlib for now as I need this code working & merged asap.

@ellisowen ellisowen force-pushed the ellisowen:npred_cube2 branch from e89a700 to 5b676c5 Aug 19, 2014

@ellisowen ellisowen force-pushed the ellisowen:npred_cube2 branch from 5b676c5 to 5c706ff Aug 19, 2014

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 19, 2014

OK, so I made the changes you suggested apart from a couple of changes to the tutorial page (mentioned above) which I would prefer to address in a future PR when this is less time-critical, if that's OK?

I also left a couple of questions above.

The Travis CI failure appears to relate to a time-out error in build_sphinx as far as I can tell.

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 20, 2014

The reason I got weird results before was that I had old example data files in my cache ... after rm -r ~/.astropy I get sensible images.

There's a few things I'd like to improve about this tutorial (e.g. use aplpy for plotting) and the code (I plan to add a Dataset class that combine counts, exposure and psf). And we should try to get really consistent npred results (within 1%).
But it's better to do this later in separate PRs ... this is ready to be merged.

Thanks for your hard work and patience with the review of this PR!

cdeil added a commit that referenced this pull request Aug 20, 2014

@cdeil cdeil merged commit 9e4164e into gammapy:master Aug 20, 2014

1 check failed

continuous-integration/travis-ci The Travis CI build failed
Details
@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 20, 2014

Great! - I agree that future PRs will be better to make these improvements, which I will work on when I am back.

@cdeil cdeil changed the title from Npred cube2 to Improve npred cube functionality Apr 8, 2015

@cdeil cdeil added the feature label Apr 8, 2015

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

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