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 tests for spectral_cube.integral_flux_image #188

Merged
merged 6 commits into from Sep 16, 2014

Conversation

Projects
None yet
2 participants
@ellisowen
Contributor

ellisowen commented Aug 26, 2014

Adds the unit test changes suggested in PR #178 comment.

@ellisowen ellisowen force-pushed the ellisowen:issue_178_test branch from 0b442cb to 289bb4d Aug 26, 2014

@@ -473,6 +473,7 @@ def compute_npred_cube(flux_cube, exposure_cube, energy_bins,
''.format(flux_cube.data.shape[1:], exposure_cube.data.shape[1:]))
energy_centers = energy_bin_centers_log_spacing(energy_bins)
energy_centers.to('GeV')

This comment has been minimized.

@ellisowen

ellisowen Aug 26, 2014

Contributor

The new bug fix and unit tests flagged up this bug too, which I fix here.

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

I think this line is a no-op:

>>> from astropy.units import Quantity
>>> e = Quantity(10, 'TeV')
>>> e
<Quantity 10.0 TeV>
>>> e.to('GeV')
<Quantity 10000.0 GeV>
>>> e
<Quantity 10.0 TeV>

Why do you think it's a bug?
If it is, please add a test and fix it.

This comment has been minimized.

@ellisowen

ellisowen Aug 26, 2014

Contributor

At some point just the value is called - and I required this to be in GeV if I remember correctly (in some instances Quantities aren't supported - e.g ImageHDUs, so I couldn't find a way of avoiding this). If I don't force the conversion in this line, the energy will just be passed in as a value in potentially the wrong energy units.

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

But energy_centers.value doesn't change if you call energy_centers.to('GeV'), no?
You'd have to write this I think!?

energy_centers = energy_centers.to('GeV')

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

Actually the idiom we should use if we need a given unit for some quantity is

energy_centers = Quantity(energy_centers, 'GeV')

but in this case please add a comment line why 'GeV' is needed (is it the assumed internal unit in some part of the class)?

This comment has been minimized.

@ellisowen

ellisowen Aug 26, 2014

Contributor

Ah yes, energy_centers = energy_centers.to('GeV') is what I should have. And yes, 'GeV' is assumed somewhere in the class.

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

The fact that your tests passed on travis-ci tells you that this is not properly tested yet.
And I don't understand where we have the internal assumption that energy needs to be GeV ... please track this down and add a comment (or even better ... do the unit conversion at the point where it's really needed instead).

This comment has been minimized.

@ellisowen

ellisowen Aug 26, 2014

Contributor

A more careful look tells me this line isn't required at all. I was mislead to believe this was a bug by some out of date code (I had run setup.py install --user but thought I had run setup.py develop --user). Apologies!

actual = new_image.data
assert_allclose(actual, expected)
# Test CDELT1

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

Reduce this repetetive code a bit, e.g.

expected = [('CDELT1', 0.5),
                    ...
                  ]
for key, value in expected:
    assert_allclose(np.abs(image.header[key]), value)
expected = 5.2481972772213124e-02
assert_allclose(actual, expected)
# Test integral flux for energy bands with units.

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

if it fits well on one line, no need to create temp expected and actual variables.
In this case simply write:

assert_allclose(image.data, image2.data)

@cdeil cdeil changed the title from Implements additional unit tests for integral_flux_image in spectral_cub... to Add more unit tests for spectral_cube.integral_flux_image Aug 26, 2014

@@ -229,7 +266,7 @@ def test_analytical_npred_cube():
# Exposure = 1, so solid angle only factor which varies.
# Result should be 0.5 * 1 * solid_angle_array from integrating analytically
energies = Quantity([1, 2], 'GeV')
energies = Quantity([1, 2], 'MeV')

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

Why can you change the unit here but tests still pass on travis-ci ... is this not covered by unit tests?
If so, please add a test for this.

This comment has been minimized.

@ellisowen

ellisowen Aug 26, 2014

Contributor

The test fails with 'GeV' as the unit, as the energy should be 'MeV'. (It was passing before because there was no conversion being done and 'MeV' was just assumed somewhere along the line)

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

But then you need to add a conversion to MeV somewhere, no?
(I don't see any relevant code change in this PR, only the test change)

This comment has been minimized.

@cdeil

cdeil Aug 26, 2014

Member

I'll be offline now and have another look at this PR tomorrow.

This comment has been minimized.

@ellisowen

ellisowen Aug 26, 2014

Contributor

So, I've had look into the code more carefully: The 'MeV' is just specific to this test which shouldn't be energy-scale independent to test units. It correctly fails with 'GeV' where it didn't before because of the previous bug fix (before this, it assumed the 'GeV' energy was 'MeV' - and in the npred code & elsewhere I used 'MeV' as the internal assumption because this is the energy unit usually used by cubes produced by the Fermi science tools).

Running the test on the version of the code in master currently fails, but this test is skipped by Travis CI because of the reproject dependency, which is why Travis CI passes.

This comment has been minimized.

@cdeil

cdeil Aug 27, 2014

Member

travis-ci does install reproject and should run any unit tests that require reproject :
https://travis-ci.org/gammapy/gammapy/jobs/33666780#L618

This comment has been minimized.

@ellisowen

ellisowen Aug 27, 2014

Contributor

I still think Travis CI skips these tests:

gammapy/data/tests/test_spectral_cube.py ssssssssssssss

There is also a scipy dependency on them (due to interpolation I think), so if reproject is OK, does Travis CI skip scipy tests?

@ellisowen ellisowen force-pushed the ellisowen:issue_178_test branch from 289bb4d to 824b178 Aug 27, 2014

actual = self.spectral_cube.integral_flux_image(energy_band).data[0, 0]
assert_quantity(actual, expected, rtol=1e-3)
assert_allclose(actual, expected, rtol=1e-3)

This comment has been minimized.

@cdeil

cdeil Aug 27, 2014

Member

expected is a Quantity here, no?
If yes, please keep using assert_quantity instead ... this way also the units are tested.

What is it's value and what are it's units?

This comment has been minimized.

@ellisowen

ellisowen Aug 27, 2014

Contributor

expected is a Quantity, however actual is not (because the ImageHDU returned by integral_flux_image no longer returned quantities).

actual = array(1.8630308796498186e-07)
expected = <Quantity 1.8642901977727886e-07 1 / (cm2 s sr)>

I have changed to an assert_quantity to test the units too as you suggest.

@ellisowen ellisowen force-pushed the ellisowen:issue_178_test branch from 824b178 to 81e72a6 Aug 27, 2014

@ellisowen ellisowen referenced this pull request Aug 28, 2014

Merged

Add iterative kernel background estimator #186

5 of 5 tasks complete
# Compute image with fraction of pixel area at equator
glat = coordinates(image)[1]
area_fraction = np.cos(np.radians(glat))
result = area_fraction * equator_area
result = area_fraction * equator_area.to('sr')
return result

This comment has been minimized.

@ellisowen

ellisowen Aug 28, 2014

Contributor

I moved this change from PR #186 so that the changes were easier to keep track of. Also changed corresponding unit tests to call the correct units.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 28, 2014

It was getting hard to keep track of the small fixes to ensure they were consistent, so I have moved them to this PR.

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Aug 28, 2014

I think these corrections are all consistent now. The tests seem to pass when I expect them too, and the results in the npred_cube are unchanged - again, as I would expect. Does this seem OK to you?

@@ -22,7 +22,7 @@ def prepare_images():
background_model = SpectralCube.read(background_file)
exposure_cube = SpectralCube.read(exposure_file)
# Add correct units
exposure_cube.data = Quantity(exposure_cube.data, '1/(cm2 deg2 s MeV)')
exposure_cube.data = Quantity(exposure_cube.data, '1/(cm2 sr s MeV)')

This comment has been minimized.

@cdeil

cdeil Aug 31, 2014

Member

I think the Fermi exposure cubes have units cm2 s and not 1/(cm2 sr s MeV), no?
(so was incorrect before and still is)

assert_allclose(actual, expected)
# Test integral flux for energy bands with units.
energy_band_check = Quantity([1000, 10000], 'MeV')

This comment has been minimized.

@cdeil

cdeil Aug 31, 2014

Member

Write this as 1 to 10 GeV ... easier to read if there's fewer 0s and everything should still work (that's the whole point of using quantities ... results don't depend on the units of inputs)

This comment has been minimized.

@ellisowen

ellisowen Sep 16, 2014

Contributor

This was the test to show the bug was fixed (where Quantity([1000, 10000], 'MeV') was giving a different result to Quantity([1, 10], 'GeV')) so in this case, I will have to keep the 0s

def test_solid_angle_image(self):
actual = self.spectral_cube.solid_angle_image[10][30]
expected = Quantity(0.24999762018018362, 'steradian')
assert_quantity(actual, expected)
expected = Quantity(0.24999762018018362, 'deg2')

This comment has been minimized.

@cdeil

cdeil Aug 31, 2014

Member

Can you compute the expected value here so that this is clear that the expected value is OK?

expected = Quantity(self.spectral_cube.wcs.wcs.cdelt[:-1].prod(), 'deg2')

ellisowen added some commits Sep 16, 2014

@ellisowen

This comment has been minimized.

Contributor

ellisowen commented Sep 16, 2014

I think this addresses your comments

cdeil added a commit that referenced this pull request Sep 16, 2014

Merge pull request #188 from ellisowen/issue_178_test
Add more unit tests for spectral_cube.integral_flux_image

@cdeil cdeil merged commit 3e14c12 into gammapy:master Sep 16, 2014

1 check passed

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

@cdeil cdeil changed the title from Add more unit tests for spectral_cube.integral_flux_image to Add tests for spectral_cube.integral_flux_image Apr 8, 2015

@cdeil cdeil added the tests label Apr 8, 2015

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

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