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 FOV background image modeling #503

Merged
merged 26 commits into from Apr 12, 2016

Conversation

Projects
None yet
4 participants
@JouvinLea
Contributor

JouvinLea commented Apr 6, 2016

This pull request contains:

  • improvements to the off data background maker class
  • a new image analysis class to compute stacked counts and background images for multiple observations
self.data_store = data_store
self.center = center
self.energy_band = energy_band
self.counts_image = SkyMap.empty(nxpix=1000, nypix=1000, binsz=0.01, xref=self.center.l.deg,

This comment has been minimized.

@JouvinLea

JouvinLea Apr 6, 2016

Contributor

Just start to implement the class but I have to leave now so I will continue tomorrow...
not sure If we have to define two members counts_image and bkg_image or if this is better for the method make_counts and make_bkg to return a skymap?
Then for the total_counts_map and total_bkg_map we will define two methods that will iterate on the obs_ids of the datastore,. Should I declare them as members of the class or is it better to have two methods make_total_count and make_total_bkg that will return the total sky map?

maps['counts'] = counts_image_total
maps['bkg'] = bkg_image_total
maps['exclusion'] = exclusion_mask

This comment has been minimized.

@JouvinLea

JouvinLea Apr 7, 2016

Contributor

@cdeil
I added some methods in the class so the code in test_curve is nicer...
Tell me what you think about the structure before I continue? Is this what you wanted?

@@ -499,6 +500,7 @@ def write(self, filename=None, header=None, **kwargs):
hdulist = fits.HDUList()
for name in self.get('_map_names', sorted(self)):
hdu = self[name].to_image_hdu()
hdu.name = name

This comment has been minimized.

@adonath

adonath Apr 7, 2016

Member

Actually this shouldn't be necessary, because the name is already passed to the ImageHDU object in to_image_hdu(). Can you explain why this is needed?

This comment has been minimized.

@JouvinLea

JouvinLea Apr 8, 2016

Contributor

I had a mistake without it because it was empty...

This comment has been minimized.

@adonath

adonath Apr 11, 2016

Member

This is really not important, but if you have time to find out, I'd be interested in why this error happens...because it might be a bug in .to_image_hdu() or somewhere else.

self.energy_band = energy_band
self.offset_band = offset_band
if not counts_image:
self.counts_image = SkyMap.empty(nxpix=1000, nypix=1000, binsz=0.01, xref=self.center.l.deg,

This comment has been minimized.

@adonath

adonath Apr 7, 2016

Member

I would suggest to only to initialize an empty self.maps = SkyMapCollection() and add the (i.e. counts, background, exposure etc. step by step to this container in the corresponding methods). The binning information (nxpix, nypix, binsz, etc.) shouldn't be hard-coded but rather passed to self.make_counts or be stored as a dictionary. To avoid repeating the binning information all the time you can use SkyMap.empty_like().

"""
self.make_total_counts()
self.make_total_bkg(exclusion_mask)
maps = SkyMapCollection()

This comment has been minimized.

@adonath

adonath Apr 7, 2016

Member

I think this belongs in self.__init__()

@cdeil cdeil added the feature label Apr 7, 2016

@cdeil cdeil added this to the 0.4 milestone Apr 7, 2016

def make_counts(self):
def __init__(self, center, energy_band, offset_band, data_store, counts_image=None, bkg_image=None):
self.data_store = data_store

This comment has been minimized.

@adonath

adonath Apr 8, 2016

Member

Conceptionally I find it a bit weird, that the class takes a DataStore as a member, I would rather expect an observation object or a list of grouped observations , i.e. maybe an ObservationTable? But I don't know enough about how the data structures are meant to be used in this context. @cdeil @joleroi ?

This comment has been minimized.

@cdeil

cdeil Apr 8, 2016

Member

This is an important question, how we handle observations and observation lists in algorithms.
I don't know yet what good patterns are.
LSST is passing around Butler objects (the DataStore equivalent) in high-level scripts, and then individual pieces in functions.

For now my suggestion would be that we allow to handle this differently how people like (whatever works), and then at the Gammapy coding spring discuss the patterns that have emerged and which one (or few) we want to adopt throughout.

This comment has been minimized.

@adonath

adonath Apr 8, 2016

Member

Sounds reasonable to me for now. Should we somehow keep track of those kind of questions? Maybe separated from normal issues e.g. in the wiki?

This comment has been minimized.

@cdeil

cdeil Apr 8, 2016

Member

@adonath - Great idea.
Can you please start a wiki page on such big-picture or design questions to be discussed at the f2f meeting (or at least in a larger group, outside a GH issue).

This comment has been minimized.

@joleroi

joleroi Apr 8, 2016

Contributor

FYI: The SpectrumExtraction class takes a DataStore as input because when I wrote it there was no observation object. Probably passing around an observation object would be the better solution (at least short term)

from reproject import reproject_interp
out = reproject_interp(
input_data=self.to_image_hdu(),
output_projection=refheader,
*args,
**kwargs,
**kwargs
)
map = SkyMap(name=self.name, data=out[0], wcs=self.wcs, unit=self.unit, meta=self.meta)

This comment has been minimized.

@adonath

adonath Apr 8, 2016

Member

Currently the wrong WCS is passed to SkyMap, it should be created from the new reference header. I would slightly prefer if reproject works "in place" i.e. changes the current skymap instead of returning a changed and implicit copy of it. We can easily add a .copy() method to the sky map class (which might be a good idea anyway...).

This comment has been minimized.

@JouvinLea

JouvinLea Apr 8, 2016

Contributor

you mean it shouldn't be wcs=self.wcs but wcs=reafheader , refheader being in the parameters of the method?

I don't really understand what you meant in the last part of your comments?

self.config = config
@classmethod
def from_yaml(cls, filename):

This comment has been minimized.

@adonath

adonath Apr 8, 2016

Member

I think in the end we'd like to have config I/O for image analysis tasks, but as this currently doesn't do anything it's probably OK to remove it...

self.total_bkg_image = SkyMap.empty(nxpix=1000, nypix=1000, binsz=0.01, xref=self.center.l.deg,
yref=self.center.b.deg, proj='TAN')
self.header = self.counts_image.to_image_hdu().header
self.solid_angle = Angle(0.01, "deg") ** 2

This comment has been minimized.

@adonath

adonath Apr 8, 2016

Member

You can get the solid angle from SkyMap.solid_angle where it's needed.

bkg = disk_correlate(bkg_image.data, 10)
s = significance(counts, bkg)
s_image = fits.ImageHDU(data=s, header=counts_image.to_image_hdu().header)
return s_image

This comment has been minimized.

@adonath

adonath Apr 8, 2016

Member

I'm not sure if the method should return anything at all, or just add additional maps to the self.maps attribute (which is a SkyMapCollection). The advantage of the later is, that all the available maps are stored in a single container and can be written to file at any time e.g. for debugging

@joleroi

This comment has been minimized.

Contributor

joleroi commented Apr 8, 2016

@JouvinLea Can you please add a description for this PR (and all others in the future 😉 ). It would make it easier for me to get an impression what's happening.

Example:
#479 (comment)

def __init__(self, nxpix=None, nypix=None, binsz=None, xref=None, yref=None, proj=None, coordsys=None,
energy_band=None, offset_band=None,
data_store=None, counts_image=None, bkg_image=None):
self.maps = SkyMapCollection()

This comment has been minimized.

@JouvinLea

JouvinLea Apr 11, 2016

Contributor

@adonath
I add the skymapcollection!
Do you think this is better now?

This comment has been minimized.

@adonath

adonath Apr 11, 2016

Member

Yes, I think using the SkyMapCollection is much more flexible and convenient. Looks much better to me!

from reproject import reproject_interp
out = reproject_interp(
input_data=self.to_image_hdu(),
output_projection=refheader,
*args,
**kwargs,
**kwargs

This comment has been minimized.

@cdeil

cdeil Apr 11, 2016

Member

@JouvinLea - This line was fixed in master already. Can you please rebase this PR to resolve the merge conflict?

This comment has been minimized.

@JouvinLea

JouvinLea Apr 11, 2016

Contributor

@cdeil
I have a syntax error when I add the coma

Traceback (most recent call last):
  File "test_curve.py", line 10, in <module>
    from gammapy.background import EnergyOffsetBackgroundModel
  File "/Users/jouvin/gammapy/gammapy/background/__init__.py", line 6, in <module>
    from .fov import *
  File "/Users/jouvin/gammapy/gammapy/background/fov.py", line 11, in <module>
    from ..image.utils import coordinates
  File "/Users/jouvin/gammapy/gammapy/image/__init__.py", line 5, in <module>
    from .profile import *
  File "/Users/jouvin/gammapy/gammapy/image/profile.py", line 7, in <module>
    from .maps import SkyMap
  File "/Users/jouvin/gammapy/gammapy/image/maps.py", line 350
    **kwargs,
            ^
SyntaxError: invalid syntax

This comment has been minimized.

@cdeil

cdeil Apr 11, 2016

Member

@JouvinLea This is fixed in master. If you fetch the latest master version and rebase this pull request, you should get a merge conflict for this line. If you resolve it (picking the solution from master), the error should be gone.

Let me know if this works.

self.energy_band = energy_band
self.offset_band = offset_band
self.image = SkyMap.empty(nxpix=nxpix, nypix=nypix, binsz=binsz, xref=xref,

This comment has been minimized.

@cdeil

cdeil Apr 11, 2016

Member

I think instead of re-exposing all (or part of) the SkyMap.empty parameters in the ImageAnalysis initialiser, it would be better to just take a skymap_ref as input, and let the caller construct it on the line before creating the ImageAnalysis object.

We might split this out into a SkyMapGeometry class that just has nxpix, binsz, ... as parameters, but for now I think using skymap_ref objects (that are by default zero-filled) for this purpose is fine.

def make_psf(self):
log.info('Making PSF ...')
def make_exposure(self):
log.info('Making Exposure ...')
def fit_source(self):

This comment has been minimized.

@cdeil

cdeil Apr 11, 2016

Member

@JouvinLea - I'd say remove the fit_source method here.

Modeling / fitting is already partly implemented elsewhere, and even if we do add a convenience "maker" class for it, it would be more flexible to have a separate one, no?
So the job of this ImageAnalysis class is to prepare the images that go in the likelihood fit (if that's what the user wants, it could also just go e.g. to TS map computation and source finding).

Lea Jouvin
Number of the observation
"""
self.maps["counts"] = SkyMap.empty_like(self.empty_image)

This comment has been minimized.

@adonath

adonath Apr 11, 2016

Member

Just wanted to mention that the SkyMapCollection class allows attribute access as well. i.e. self.maps.counts instead of self.maps['counts'], but I don't have a preference what to use...

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

I would suggest to use a local variable counts here or counts_map, and to only store it as a data member at the end of the method:

self.maps['counts'] = counts

Also, the log.info statement should be moved to the top of the method, and I'd suggest changing to log.debug.

self.maps["bkg"].data = scale * self.maps["bkg"].data
self.maps['exclusion'] = exclusion_mask
def make_total_bkg(self, exclusion_mask):

This comment has been minimized.

@adonath

adonath Apr 11, 2016

Member

As exclusion_mask is needed for several methods, maybe pass it already to __init__ and add it to self.maps?

Right now this function stacks the background for all observations in the data store. Is this the desired behaviour? Or should the function maybe take a list of observations IDs as well (defaulting to all observations) and allow which observations to use for the background modelling?

This comment has been minimized.

@adonath

adonath Apr 11, 2016

Member

Does make_background_normalized_offcounts actually have to be visible to the user?
Because right now make_total_bkg only does the a loop over observation IDs. But if a list of obs IDs is passed, it can still contain only one value.

This comment has been minimized.

@JouvinLea

JouvinLea Apr 12, 2016

Contributor

@adonath
Yes this is the desired behaviour! Here the bkg model is already build in the table["Acceptance"] for each observation. And we just stack the counts bkg model for each observation as for the counts map.
But what will be better is to give a data_store and an obs_table with only the observation for which we want to compute the bkg and counts map.

This comment has been minimized.

@adonath

adonath Apr 12, 2016

Member

OK, that makes sense.

counts_sum = np.sum(self.maps["counts"].data * exclusion_mask.data)
bkg_sum = np.sum(self.maps["bkg"].data * exclusion_mask.data)
scale = counts_sum / bkg_sum
self.maps["bkg"].data = scale * self.maps["bkg"].data

This comment has been minimized.

@adonath

adonath Apr 11, 2016

Member

I'm not sure if the background for a single observation should be stored in self.maps, because if you call make_total_bkg(), self.maps['bkg'] will be continuously overwritten and in the end just contain the background for the last observation in the list, which is kind of arbitrary. But I'm not yet sure what a good solution is, I have to think about that a bit more.

Lea Jouvin added some commits Apr 12, 2016

Lea Jouvin
exclusion maps + exposure
Please enter the commit message for your changes. Lines starting
Lea Jouvin
Lea Jouvin
ok
Lea Jouvin
@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 12, 2016

@cdeil
I know that some test are failing but I'm not sure this is really linked to what I implemented!
Do you think we can merge and I can continue to implement the exposure map!

maps['counts'] = counts_image_total
maps['bkg'] = bkg_image_total
maps['exclusion'] = exclusion_mask
exclusion_mask = ExclusionMask.read('$GAMMAPY_EXTRA/datasets/exclusion_masks/tevcat_exclusion.fits')

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

Why read the same exclusion mask a second time here?
Should this and the following line be removed?

@@ -183,9 +189,14 @@ def filename(self, modeltype, group_id):
Type of the background modelisation
group_id : int
number of the background model group
smooth : str

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

'str' -> 'bool'

@@ -234,20 +246,21 @@ def make_bkg_index_table(self, data_store, modeltype, out_dir_background_model=N
directory where are located the backgrounds models for each group
filename_obs_group_table : str
name of the file containing the `~astropy.table.Table` with the group infos
smooth : str

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

'str' -> 'bool'

@@ -292,6 +305,8 @@ def make_total_index_table(self, data_store, modeltype, out_dir_background_model
directory where are located the backgrounds models for each group
filename_obs_group_table : str
name of the file containing the `~astropy.table.Table` with the group infos
smooth : str

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

'str' -> 'bool'

Number of the observation
"""
self.maps["bkg"] = SkyMap.empty_like(self.empty_image)

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

Again, I would suggest to always do the computation first using a local variable, and to only store it as a data member on the last line.
I find this a little simpler, and also, if there is some error, it leaves the object in an unmodified state, which is better.

center = obs.pointing_radec.galactic
bkg_hdu = fill_acceptance_image(self.header, center, table["offset"], table["Acceptance"], self.offset_band[1])
livetime = obs.observation_live_time_duration
self.maps["bkg"].data = Quantity(bkg_hdu.data, table["Acceptance"].unit) * self.maps[

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

These three lines are a bit hard to read.
I'd suggest doing the computation using a local variable and only storing it in self.maps["bkg"] at the end.

Parameters
----------
obs_id : int
Number of the observation

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

" Number of the observation" -> "Observation ID"

Number of the observation
"""
self.make_counts(obs_id)

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

Move the scale factor computation in a separate method?
So that the caller can decide whether it should be applied or not?

If you don't think it's worth exposing as a public extra method, you could make it private (but then experts can still compute it without modifying the object if they like), i.e. _background_norm_factor.

def make_total_bkg(self):
"""Stack the total bkg from the observation in the 'DataStore'.
"""
self.maps["total_bkg"] = SkyMap.empty_like(self.empty_image)

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

Same comment as everywhere: compute first, store as data member at the end.

Disk radius in pixels.
"""
self.maps["significance"] = SkyMap.empty_like(self.empty_image)

This comment has been minimized.

@cdeil

cdeil Apr 12, 2016

Member

Same comment as above: compute first, store as data member at the end.
Move log.info statement on the first line and change to log.debug.

Also: return significance image, so that callers can get it directly on the line where it's computed if they like.
If the prefer, they can also access it from the maps data member later.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 12, 2016

@JouvinLea - I've done another round of code review ... all easy to address I think.

Is there a test that exercises ImageAnalysis.make_all?
Can you add one and add one assert, e.g. on the significance value in the pixel containing the Crab nebula?
Or do you want to do this after the HESS collab meeting and merge ASAP?

@cdeil cdeil changed the title from Some changes to the OFFBkgMaker class to Add FOV background image modeling Apr 12, 2016

Lea Jouvin
@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 12, 2016

There is no test... Yeah I would prefer to do it after the HESS collaboration! You have to remind me to do it...
I would prefer to merge ASAP and continue to work on the exposure!
@cdeil Normally I took all you comments into account, do you think we can merge?

Lea Jouvin
@cdeil

This comment has been minimized.

Member

cdeil commented Apr 12, 2016

OK, merging this now.

If there are build fails on travis-ci I can have a look tonight and fix them.

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 12, 2016

Thank you ... can't wait to use this for my analyses!

@cdeil cdeil merged commit 3b46ba1 into gammapy:master Apr 12, 2016

0 of 2 checks passed

continuous-integration/appveyor/pr Waiting for AppVeyor build to complete
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 12, 2016

hahaha thanks to be so enthusiastic for this class!!:-)

@cdeil cdeil referenced this pull request Apr 12, 2016

Closed

FOV background method #5

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