PSF component analysis with HealPix pixelization

[zimmerst]

Hi, this could be seen as clone of #203, which I just closed.

While I confirmed that gta.setup() works fine for arbitrary WCS projections with the 'component'-type summed likelihood, it seems that the current fermipy version does not support this when moving to healpix-based binning.

this is my configuration:

data:
  evfile : ft1.lst
  scfile : ft2.lst

binning:
  roiwidth   : 20.0
  binsz      : 0.1
  binsperdec : 8
  projtype   : HPX

selection :
  logemin: 2.
  logemax : 5.9030899869919438
  zmax    : 90
  evclass : 128
  evtype  : 3
  tmin    : 239557414
  tmax    : 428903014
  filter  : null
  ra      : 0
  dec     : 0

components:
  - model: {isodiff: $FERMI_DIFFUSE_DIR/P305/isotropic_p8r3_source_v2_psf0_zmax80_4years.txt}
    selection: {evtype: 4, zmax: 80}
    data: {ltcube: ltcube_00.fits}
  - model: {isodiff: $FERMI_DIFFUSE_DIR/P305/isotropic_p8r3_source_v2_psf1_zmax80_4years.txt}
    selection: {evtype: 8, zmax: 85}
    data: {ltcube: ltcube_01.fits}
  - model: {isodiff: $FERMI_DIFFUSE_DIR/P305/isotropic_p8r3_source_v2_psf2_zmax80_4years.txt}
    selection: {evtype: 16, zmax: 95}
    data: {ltcube: ltcube_02.fits}
  - model: {isodiff: $FERMI_DIFFUSE_DIR/P305/isotropic_p8r3_source_v2_psf3_zmax80_4years.txt}
    selection: {evtype: 32, zmax: 100}
    data: {ltcube: ltcube_03.fits}

gtlike:
  edisp : True
  irfs : 'P8R3_SOURCE_V2'
  edisp_disable : ['isodiff','galdiff']

model:
  src_roiwidth : 25.0
  galdiff  : '$FERMI_DIFFUSE_DIR/gll_iem_v06.fits'
  #isodiff  : 'iso_P8R2_SOURCE_V6_v06.txt'
  catalogs : ['FL8Y.fit']

tsmap:
  nthread : 10
  multithread : true

I can initialize the GTAnalysis object as usual:

2018-04-26 09:32:40 INFO    GTAnalysis.__init__(): 
--------------------------------------------------------------------------------
fermipy version 0.16.0+179.g1a8e 
ScienceTools version ScienceTools-11-07-00
WARNING: hdu= was not specified but multiple tables are present, reading in first available table (hdu=1) [astropy.io.fits.connect]
WARNING: UnitsWarning: 'photon/cm**2/MeV/s' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
WARNING: UnitsWarning: 'photon/cm**2/s' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
WARNING: UnitsWarning: 'erg/cm**2/s' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]

and then call gta.setup()

-- content clipped --
2018-04-26 09:39:50 INFO    GTAnalysis.setup(): Running setup.
2018-04-26 09:39:50 INFO    GTBinnedAnalysis.setup(): Running setup for component 00
2018-04-26 09:39:50 INFO    GTBinnedAnalysis._select_data(): Skipping data selection.
2018-04-26 09:39:50 INFO    GTBinnedAnalysis.setup(): Using external LT cube.
2018-04-26 09:39:51 INFO    GTBinnedAnalysis._create_expcube(): Skipping gtexpcube.
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
2018-04-26 09:39:51 INFO    GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps.
2018-04-26 09:39:51 INFO    GTBinnedAnalysis.setup(): Finished setup for component 00
2018-04-26 09:39:51 INFO    GTBinnedAnalysis.setup(): Running setup for component 01
2018-04-26 09:39:51 INFO    GTBinnedAnalysis._select_data(): Skipping data selection.
2018-04-26 09:39:51 INFO    GTBinnedAnalysis.setup(): Using external LT cube.
2018-04-26 09:39:52 INFO    GTBinnedAnalysis._create_expcube(): Skipping gtexpcube.
2018-04-26 09:39:52 INFO    GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps.
2018-04-26 09:39:52 INFO    GTBinnedAnalysis.setup(): Finished setup for component 01
2018-04-26 09:39:52 INFO    GTBinnedAnalysis.setup(): Running setup for component 02
2018-04-26 09:39:52 INFO    GTBinnedAnalysis._select_data(): Skipping data selection.
2018-04-26 09:39:52 INFO    GTBinnedAnalysis.setup(): Using external LT cube.
2018-04-26 09:39:53 INFO    GTBinnedAnalysis._create_expcube(): Skipping gtexpcube.
2018-04-26 09:39:53 INFO    GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps.
2018-04-26 09:39:53 INFO    GTBinnedAnalysis.setup(): Finished setup for component 02
2018-04-26 09:39:53 INFO    GTBinnedAnalysis.setup(): Running setup for component 03
2018-04-26 09:39:53 INFO    GTBinnedAnalysis._select_data(): Skipping data selection.
2018-04-26 09:39:53 INFO    GTBinnedAnalysis.setup(): Using external LT cube.
2018-04-26 09:39:54 INFO    GTBinnedAnalysis._create_expcube(): Skipping gtexpcube.
2018-04-26 09:39:54 INFO    GTBinnedAnalysis.run_gtapp(): Running gtsrcmaps.
2018-04-26 09:47:51 INFO    GTBinnedAnalysis.setup(): Finished setup for component 03
2018-04-26 09:47:51 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00.
2018-04-26 09:48:16 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 01.
2018-04-26 09:48:40 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 02.
2018-04-26 09:49:04 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 03.

but then things crash violently:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-5-1e27b2b06e27> in <module>()
----> 1 gta.setup()

/u/gl/zimmer/.local/lib/python2.7/site-packages/fermipy-0.16.0+179.g1a8e-py2.7.egg/fermipy/gtanalysis.pyc in setup(self, init_sources, overwrite, **kwargs)
   1035 
   1036         # Create likelihood
-> 1037         self._create_likelihood()
   1038 
   1039         # Determine tmin, tmax

/u/gl/zimmer/.local/lib/python2.7/site-packages/fermipy-0.16.0+179.g1a8e-py2.7.egg/fermipy/gtanalysis.pyc in _create_likelihood(self, srcmdl)
   1065         self.like.model = self.like.components[0].model
   1066         self._fitcache = None
-> 1067         self._init_roi_model()
   1068 
   1069     def _init_roi_model(self):

/u/gl/zimmer/.local/lib/python2.7/site-packages/fermipy-0.16.0+179.g1a8e-py2.7.egg/fermipy/gtanalysis.pyc in _init_roi_model(self)
   1079         proj_type = 0
   1080         for c, crd in zip(self.components, rm['components']):
-> 1081             cm = c.counts_map()
   1082             wm = c.weight_map()
   1083             cmaps += [cm]

/u/gl/zimmer/.local/lib/python2.7/site-packages/fermipy-0.16.0+179.g1a8e-py2.7.egg/fermipy/gtanalysis.pyc in counts_map(self)
   4753             z = cmap.data()
   4754             z = np.array(z).reshape(self.enumbins, self.npix, self.npix)
-> 4755             return WcsNDMap(copy.deepcopy(self.geom), z)
   4756         elif p_method == 1:  # HPX
   4757             z = cmap.data()

/nfs/farm/g/glast/u/zimmer/software/anaconda/lib/python2.7/site-packages/gammapy/maps/wcsnd.pyc in __init__(self, geom, data, dtype, meta)
     46         elif data.shape != shape[::-1]:
     47             raise ValueError('Wrong shape for input data array. Expected {} '
---> 48                              'but got {}'.format(shape, data.shape))
     49 
     50         super(WcsNDMap, self).__init__(geom, data, meta)

ValueError: Wrong shape for input data array. Expected (95582, 95582, 31) but got (31, 200, 200)

without having looked at the code in detail, the 2nd last item in the call-back suggests that the code crashed outside the HPX frame (since it never executed the elif p_method == 1: #HPX instruction)

[zimmerst]

for completeness, with this config, I get the following binning dictionary:

{'binsperdec': 8.0,
 'binsz': 0.1,
 'coordsys': 'CEL',
 'enumbins': None,
 'hpx_ebin': True,
 'hpx_order': 10,
 'hpx_ordering_scheme': 'RING',
 'npix': None,
 'proj': 'AIT',
 'projtype': 'HPX',
 'roiwidth': 20.0}

[cdeil]

Just a side note: this is very confusing, suggesting an axis ordering issue:

ValueError: Wrong shape for input data array. Expected (95582, 95582, 31) but got (31, 200, 200)

But actually it was just the error message that printed the shape with two different axis orders. I improved this in https://github.com/gammapy/gammapy/pull/1346/files#diff-e3f81d40e68a44a1d6cf5207f0e98aeaR39

It still leaves the issue in Fermipy reported by @zimmerst here. Should HPX work for that computation, or is it just not supported?

[cdeil]

@adonath and I are about to release Gammapy v0.8; it would be good to know if this issue has been resolved.

@zimmerst or @dimauromattia - Could you please have a look? Maybe even add a regression test in Fermipy that shows if there is or isn't an issue with HPX maps in gammapy.maps?

[zimmerst]

@cdeil Thanks for the update. I'll need a moment to change back to my HealPix based analysis. I had abandoned this path a few months ago so I first need to reset myself to see if this has been resolved in GammaPy.

[zimmerst]

I tried the development versions of gammapy & fermipy respectively:

print "fermipy: %s  --  gammapy: %s"%(fermipy.__version__, gammapy.__version__)
fermipy: 0.17.3+17.g6cbb  --  gammapy: 0.8.dev7133

the code still crashes when calling gta.setup() but differently than before, for the first time revealing what seems to be the underlying cause:

2018-08-31 08:25:37 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00.
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-3-1e27b2b06e27> in <module>()
----> 1 gta.setup()

/root/.local/lib/python2.7/site-packages/fermipy/gtanalysis.pyc in setup(self, init_sources, overwrite, **kwargs)
   1042 
   1043         # Create likelihood
-> 1044         self._create_likelihood()
   1045 
   1046         # Determine tmin, tmax

/root/.local/lib/python2.7/site-packages/fermipy/gtanalysis.pyc in _create_likelihood(self, srcmdl)
   1067         self._like = SummedLikelihood()
   1068         for c in self.components:
-> 1069             c._create_binned_analysis(srcmdl)
   1070             self._like.addComponent(c.like)
   1071 

/root/.local/lib/python2.7/site-packages/fermipy/gtanalysis.pyc in _create_binned_analysis(self, xmlfile, **kwargs)
   5289         self.logger.debug(kw)
   5290 
-> 5291         self._obs = ba.BinnedObs(**utils.unicode_to_str(kw))
   5292 
   5293         # Create BinnedAnalysis

/home/ScienceTools/python/BinnedAnalysis.py in __init__(self, srcMaps, expCube, binnedExpMap, irfs, phased_expmap)
     59             self.irfs = irfs
     60 
---> 61         pyLike.AppHelpers_checkExposureMap(srcMaps, binnedExpMap)
     62         # EAC switch to using AppHelpers...
     63         self.countsMap = pyLike.AppHelpers_readCountsMap(srcMaps)

/home/ScienceTools/python/pyLikelihood.py in AppHelpers_checkExposureMap(*args)
   3548 def AppHelpers_checkExposureMap(*args):
   3549   """AppHelpers_checkExposureMap(string cmapfile, string emapfile)"""
-> 3550   return _pyLikelihood.AppHelpers_checkExposureMap(*args)
   3551 
   3552 def AppHelpers_checkExposureMap_wcs(*args):

RuntimeError: HistND::setSlice:
length of input data != histogram slice size.

[zimmerst]

Addendum: in the log file from the preparation stage, I found this line:

"2018-08-30 06:40:14 DEBUG GTBinnedAnalysis._create_expcube(): Skipping local gtexpcube for HEALPix"

Can someone explain to me why there's no local gtexpcube? Shouldn't this one still be generated, regardless of the projection?

[cdeil]

That looks like a Fermi ST problem, I have no idea there.

If you can find an issue where you think the problem is in gammapy.maps, please let me know, and if you can, please try to extract a snipped of Gammapy code that lets me reproduce the issue.

[zimmerst]

@cdeil I agree... I'm trying to get more information on the issue and will keep you posted.