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 example script for cube background models #338

Merged
merged 19 commits into from Sep 3, 2015

Conversation

Projects
None yet
2 participants
@mapazarr
Member

mapazarr commented Aug 25, 2015

The script uses the same model (except for absolute normalization) to produce a true cube bg model and a reco cube bg model. The models can be used to test the cube bg model production and can be compared to each other using the plot_bg_cube_model_comparison.py example script.

@mapazarr mapazarr force-pushed the mapazarr:plots_for_last_gsoc_post branch from dafe1fa to 8f49ec9 Aug 25, 2015

@mapazarr

This comment has been minimized.

Member

mapazarr commented Aug 25, 2015

Und noch einen Kaffee! :-)
Mergen?

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 26, 2015

@mapazarr – Thanks for adding this!

As with your other recent PRs I'm OK with merging this without much review, but still I'd like to do some review to avoid accumulating more and more technical debt in Gammapy, i.e. things that need to be reviewed / tested / cleaned up later.

@cdeil cdeil added the feature label Aug 26, 2015

@cdeil cdeil added this to the 0.4 milestone Aug 26, 2015

@cdeil cdeil self-assigned this Aug 26, 2015

@cdeil cdeil changed the title from Added example script to produce true and reco bg cube model. to Add example script to produce true and reco bg cube model Aug 26, 2015

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 26, 2015

Please don't use minus - in Python file names, use underscore _ instead.

outdir = OUTDIR
overwrite = OVERWRITE
# create output folder

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

These ~10 lines of code have appeared now a few times in your code I think.
Please make it a utility function.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

@@ -809,6 +809,38 @@ def plot_spectrum(self, coord, ax=None, style_kwargs=None):
plt.close(fig)
return ax
@property
def integral(self):
"""Integral of the cube (`~astropy.units.Quantity`)"""

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

What is the "Integral of the cube"?

I think it's predicted counts?

If so you could call it npred or npred_total (see
https://gammapy.readthedocs.org/en/latest/api/gammapy.data.SpectralCube.html#gammapy.data.SpectralCube.integral_flux_image and
https://gammapy.readthedocs.org/en/latest/api/gammapy.data.compute_npred_cube.html),
or at least mention this in the docstring.

Also maybe add "scalar" to make it clearer?

(scalar `~astropy.units.Quantity`)

This comment has been minimized.

@mapazarr

mapazarr Aug 26, 2015

Member

It depends on the contents of the cube. If it is counts, it retruns a scalar with the total number of counts, but if it is for instance the bg rate it will return the number of counts corrected by the livetime in units of 1/s.

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

OK, then just say that in the docstring, e.g. "the integral returned has the unit of the values in the cube times solid angle times energy". You could even mention the one or two common examples that we have.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

az_range=Angle([0, 360], 'degree'),
alt_range=Angle([45, 90], 'degree'),
date_range=(Time('2010-01-01T00:00:00',
format='isot', scale='utc'),

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Remove format='isot', scale='utc' ... UTC is the default and the format it's important.

This comment has been minimized.

@mapazarr

mapazarr Aug 26, 2015

Member

format='isot', scale='utc' is the standard we defined in file:///home/mapaz/astropy/development_code/gammapy/docs/_build/html/time/index.html#time-handling-in-gammapy
It might be the default in Astropy right now (I'm not sure). But if it changes at some point for some reason, we would have to change all entries throughout Gammapy, to be consistent with our standard. Isn't it better to just explicitly give the format and scale?

BTW, what do you mean by `it's important'?

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Let's adopt the Astropy default for Gammapy and trust that they won't change it.
(It's highly unlikely ... it would break a lot of user codes if they did ... satellites falling out of the sky ... doomsday.)

I wanted to write "the format is not important", because it's not used for computations.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done (everywhere in gammapy, again).
You'll be responsible for ordering the doomsday :-)

"""Integral of the cube images (`~astropy.units.Quantity`)
Calculate the integral of each energy bin (slice) in the
cube. Returns an array of integrals.

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Again, from the method name and docstring it's not obvious to me what this does.

It looks to me like this is computing a 1-dim npred array for each energy bin?

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done (similar to the previous one).

Starting date for random generation of observation start time.
dateend : `~astropy.time.Time`, optional
Ending date for random generation of observation start time.
az_range : `~astropy.coordinatesAngle`, optional

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Missing dot: ~astropy.coordinatesAngle -> ~astropy.coordinates.Angle

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Good catch!
Fixed.

az_range : `~astropy.coordinatesAngle`, optional
Azimuth angle range (start, end) for random generation of
observation pointing positions.
alt_range : `~astropy.coordinatesAngle`, optional

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Missing dot again.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Fixed.

datestart=None, dateend=None,
az_range=Angle([0, 360], 'degree'),
alt_range=Angle([45, 90], 'degree'),
date_range=(Time('2010-01-01T00:00:00',

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Simplify:

date_range=Time(['2010-01-01', '2015-01-01'])

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

It is already done... You commented on the old implementation.

Starting date for random generation of observation start time.
dateend : `~astropy.time.Time`, optional
Ending date for random generation of observation start time.
az_range : `~astropy.coordinatesAngle`, optional

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Again: missing docs (also two lines down)

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Dots. Fixed.

@@ -311,8 +311,7 @@ def stack_observations(fits_path, outdir, overwrite, method='default'):
bg_cube_model = make_bg_cube_model(observation_table, fits_path, method)
# save model to file
outfile = outdir +\
'/bg_cube_model_group{}'.format(group)
outfile = outdir + '/bg_cube_model_group{}'.format(group)

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Use os.path.join ... this doesn't work on Windows. 😄

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done. Maybe we need an issue 'support windows'? ;-)
And a travis-ci test :-D

This comment has been minimized.

@cdeil

cdeil Sep 3, 2015

Member

Oh ... we will support Windows!
And, to avoid the appearance of bias, we shall support Apple watch and MS DOS as well!

# create output folder
outdir = OUTDIR + '/reco'
if not os.path.isdir(outdir):

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Again the same code to make a folder ... use utility function.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

CHANGES.rst Outdated
- Fix TS map boundary handling [#332] (Axel Donath)
- Production of true/reco bg cube models should use the same model. [#335] (Manuel Paz Arribas)
- Fix TS map boundary handling [#332] (Axel Donath)
- Added example script to produce true and reco bg cube model. [#338] (Manuel Paz Arribas)

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Let's use a uniform format for GH issue titles and changelog. Active voice ("Add" instead of "Added") and no dot at the end of the line.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

overwrite = OVERWRITE
# create output folder

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Use utility function.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

def create_dummy_observation_grouping():
"""Define dummy observation grouping."""
altitude_edges = ALT_RANGE

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

I think this function can be shortened by half:

def create_dummy_observation_grouping():
    alt_ax = ObservationGroupAxis('ALT', ALT_RANGE, 'bin_edges')
    az_ax = ObservationGroupAxis('AZ', AZ_RANGE, 'bin_edges')
    obs_groups = ObservationGroups([alt_ax, az_ax])
    obs_groups.obs_groups_table['GROUP_ID'][0] = GROUP_ID

    return obs_groups

The following functions can also be simplified / improved.

Just in case you say "this is just an example and not worth the effort of code review / polishing", I'd say "especially for examples it's important to make them nice and clean" ... it's how users / co-developers learn and what determines whether they like Gammapy or now.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

obs_groups.write(outfile)
# use binning from CubeBackgroundModel.define_cube_binning
detx_range = (Angle(-0.07, 'radian').to('degree'),

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

If you're using the same binning for x and y, just have one det_range variable and use it twice in the call below.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

# 0: do not normalize
# 1: normalize w.r.t. cube integral
# 2: normalize w.r.t images integral (normalize each image on its own)
input_dir1 = '/home/mapaz/astropy/working_dir/gammapy_scripts/20150819_ready_to_merge_PR319/bg_cube_models_gammapy_a_la_michi'

This comment has been minimized.

@cdeil

cdeil Aug 26, 2015

Member

Hard-coded PATHs make it hard and unlikely that others will run the example.
So it would be much better if you make this a command line option.
For quick scripts it's OK to just use input_dir = sys.argv[1] ... no need to use an ArgumentParser, so it's just one or two extra lines to change this.

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Done.

@cdeil

This comment has been minimized.

Member

cdeil commented Aug 27, 2015

@mapazarr – There's a merge conflict (probably in the changelog) now that #337 is merged. Please rebase against master.

@mapazarr mapazarr force-pushed the mapazarr:plots_for_last_gsoc_post branch from 3fbfe33 to f1e6bee Sep 3, 2015

@mapazarr

This comment has been minimized.

Member

mapazarr commented Sep 3, 2015

Ok, now true and reco are very similar:
https://hess-confluence.desy.de/confluence/display/HESS/HESS+Open+Source+Tools+-+HOWTO+create+cube+background+models+using+Gammapy#HESSOpenSourceTools-HOWTOcreatecubebackgroundmodelsusingGammapy-Comparisontrue/reco

bg_cube_model_comparison_group27_sim_improved_true_method

I'm not sure if we still need to investigate the differences. @cdeil what is your opinion?

For me, I think I have reached my goal in this PR: have similar true/reco models. I will implement the comments from the code review tomorrow.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Sep 3, 2015

@cdeil there seems to be some incompatibility issues between numpy 1.6 or 1.7 and astropy in the travis-ci allowed failure jobs.

I killed the jobs, because they were stuck and taking forever: > 25 min.

@cdeil

This comment has been minimized.

Member

cdeil commented Sep 3, 2015

I think the best is to finish this PR this week and leave the investigation of the remaining differences true / reco to the future. IMO this is just task number 21 😄 in #333 ... I'll try to work on that myself in the coming months or get a student.

@cdeil

This comment has been minimized.

Member

cdeil commented Sep 3, 2015

@cdeil there seems to be some incompatibility issues between numpy 1.6 or 1.7 and astropy in the travis-ci allowed failure jobs.

I'm short on time this week, my suggestion would be to just ignore this.
(or file a separate issue if you like.)

@mapazarr

This comment has been minimized.

Member

mapazarr commented on gammapy/background/cube.py in e0a3b68 Sep 3, 2015

This takes care that MeV units will be used for the bg rates per default. There is though a flag to avoid this if desired (turned off per default).

This comment has been minimized.

Member

cdeil replied Sep 3, 2015

This seems like a random location to force units to a default.
Shouldn't this happen at the time of object construction or read / write?

This comment has been minimized.

Member

mapazarr replied Sep 3, 2015

It is not random: this is the function where the bg cube gets its bg rate units, and is called each time a cube is created.
If we put it in read or write, the cube might still be using other units in the meanwhile. I think read should read the cube as it is in the file, and write should only write, not interpret the results.

But now that you mention it, this is not the right place, since this function is inside the generic Cube container class, and we don't want to force MeV for all cubes, right? I'll move the units forcing to the make_bg_cube_model function.

if do_not_force_mev_units:
bin_volume = delta_energy*(delta_y*delta_x).to('sr')
else:
bin_volume = delta_energy.to('MeV')*(delta_y*delta_x).to('sr')

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

This takes care that MeV units will be used for the bg rates per default. There is though a flag to avoid this if desired (turned off per default).

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Further discussion on this topic should go here: e0a3b68

@@ -12,7 +12,7 @@
]
def make_bg_cube_model(observation_table, fits_path, method='default'):
def make_bg_cube_model(observation_table, fits_path, method='default', do_not_force_mev_units=False):

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

OK, I moved the enforcing of the units here. This seems the right place, since it is the function that creates the cubes.

Set to `True` to use the same energy units as the energy
binning for the bg rate. (Only applicable to the *default*
method; the *michi* always enforces `MeV` in the bg rate
units.)

This comment has been minimized.

@cdeil

cdeil Sep 3, 2015

Member

Can you make this an option unit='1 / (MeV sr s)' that is applied for any method?
Simpler, more consistent (across methods) and more flexible, no?

This comment has been minimized.

@mapazarr

mapazarr Sep 3, 2015

Member

Well, actually the whole idea is to be LESS flexible, since we are enforcing some particular units.

In addition, the michi method should always use these units, since we want it to be as close to @mimayer's results as possible, in order to compare, so I think is is fine how it is.

The user can always unset this forcing of the units, by setting the varaible to True, in which case, the same units as the ones defined by the binning will be used (the most natural choice). If he wants something else, he will have to do it manually.

So, I'm leaving this as it is.

@cdeil

This comment has been minimized.

Member

cdeil commented Sep 3, 2015

@mapazarr – I'm going to bed in 10 min. Do you have any question here or are you all set?

Feel free to merge if you think this is ready or let me know if you want me to have another look tomorrow morning ...

@mapazarr

This comment has been minimized.

Member

mapazarr commented Sep 3, 2015

Ok, I addressed all comments, and travis-ci is green, so I'm merging this now.

mapazarr added a commit that referenced this pull request Sep 3, 2015

Merge pull request #338 from mapazarr/plots_for_last_gsoc_post
Add example script to produce true and reco bg cube model

@mapazarr mapazarr merged commit a8acc4a into gammapy:master Sep 3, 2015

1 check passed

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

This comment has been minimized.

Member

mapazarr commented Sep 6, 2015

@cdeil I forgot one last thing. I'd like to add the plot comparing true/reco:

bg_cube_model_comparison_group27_sim_improved_true_method

to the high-level docs, at the end of:

http://gammapy.readthedocs.org/en/latest/background/make_models.html#datasets-for-testing

I think it would help illustrate the capabilities of bg cube production. The file is 260 kB; is it ok, or is it too big?

@cdeil

This comment has been minimized.

Member

cdeil commented Sep 6, 2015

@mapazarr – That is a nice Figure and it would be good to have in the docs.

But indeed, I'd prefer to keep the gammapy repo small and to have all example datasets and images such as this one in gammapy-extra. So please put this on the list of things to do for background modeling.

If you like, you can already put the image here for now:
https://github.com/gammapy/gammapy-extra/tree/master/figures
and then we'll properly include it in the Sphinx build later, when gammapy-extra is integrated better.

@mapazarr

This comment has been minimized.

Member

mapazarr commented Sep 7, 2015

Ok, adressing the figure issue in #347.

@cdeil cdeil changed the title from Add example script to produce true and reco bg cube model to Add example script to produce cube background models Sep 16, 2015

@cdeil cdeil changed the title from Add example script to produce cube background models to Add example script for cube background models Sep 16, 2015

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