This is a fixed-text formatted version of a Jupyter notebook.

3D simulation and fitting

This tutorial shows how to do a 3D map-based simulation and fit.

For a tutorial on how to do a 3D map analyse of existing data, see the analysis_3d tutorial.

This can be useful to do a performance / sensitivity study, or to evaluate the capabilities of Gammapy or a given analysis method. Note that is is a binned simulation as is e.g. done also in Sherpa for Chandra, not an event sampling and anbinned analysis as is done e.g. in the Fermi ST or ctools.

Imports and versions

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord, Angle
from gammapy.utils.random import get_random_state
from gammapy.irf import EffectiveAreaTable2D, EnergyDispersion2D, EnergyDependentMultiGaussPSF, Background3D
from gammapy.maps import WcsGeom, MapAxis, WcsNDMap, Map
from gammapy.spectrum.models import PowerLaw
from gammapy.image.models import SkyGaussian
from gammapy.cube.models import SkyModel, SkyModels
from gammapy.cube import MapFit, MapEvaluator, PSFKernel
from gammapy.cube import make_map_exposure_true_energy, make_map_background_irf
In [3]:
!gammapy info --no-envvar --no-dependencies --no-system

Gammapy package:

        path                   : /Users/deil/work/code/gammapy/gammapy
        version                : 0.8.dev6700
        githash                : 1c8db0e996201e378004568a57ffc159d6c6e415

Simulate

In [4]:
def get_irfs():
    """Load CTA IRFs"""
    filename = '$GAMMAPY_EXTRA/datasets/cta-1dc/caldb/data/cta//1dc/bcf/South_z20_50h/irf_file.fits'
    psf = EnergyDependentMultiGaussPSF.read(filename, hdu='POINT SPREAD FUNCTION')
    aeff = EffectiveAreaTable2D.read(filename, hdu='EFFECTIVE AREA')
    edisp = EnergyDispersion2D.read(filename, hdu='ENERGY DISPERSION')
    bkg = Background3D.read(filename, hdu='BACKGROUND')
    return dict(psf=psf, aeff=aeff, edisp=edisp, bkg=bkg)

irfs = get_irfs()
In [5]:
import os
os.environ['GAMMAPY_EXTRA']
Out[5]:
'/Users/deil/work/code/gammapy-extra'
In [6]:
# Define sky model to simulate the data
spatial_model = SkyGaussian(
    lon_0='0.2 deg',
    lat_0='0.1 deg',
    sigma='0.3 deg',
)
spectral_model = PowerLaw(
    index=3,
    amplitude='1e-11 cm-2 s-1 TeV-1',
    reference='1 TeV',
)
sky_model = SkyModel(
    spatial_model=spatial_model,
    spectral_model=spectral_model,
)
print(sky_model)
SkyModel

spatial_model = SkyGaussian

Parameters:

         name   value   error unit min max frozen
        ----- --------- ----- ---- --- --- ------
        lon_0 2.000e-01   nan  deg nan nan  False
        lat_0 1.000e-01   nan  deg nan nan  False
        sigma 3.000e-01   nan  deg nan nan  False

spectral_model = PowerLaw

Parameters:

           name     value   error       unit         min    max frozen
        --------- --------- ----- --------------- --------- --- ------
            index 3.000e+00   nan                       nan nan  False
        amplitude 1.000e-11   nan 1 / (cm2 s TeV)       nan nan  False
        reference 1.000e+00   nan             TeV 0.000e+00 nan   True

In [7]:
# Define map geometry
axis = MapAxis.from_edges(
    np.logspace(-1., 1., 10), unit='TeV', name='energy',
)
geom = WcsGeom.create(
    skydir=(0, 0), binsz=0.02, width=(5, 4),
    coordsys='GAL', axes=[axis],
)
In [8]:
# Define some observation parameters
# Here we just have a single observation,
# we are not simulating many pointings / observations
pointing = SkyCoord(1, 0.5, unit='deg', frame='galactic')
livetime = 1 * u.hour
offset_max = 2 * u.deg
offset = Angle('2 deg')
In [9]:
exposure = make_map_exposure_true_energy(
    pointing=pointing,
    livetime=livetime,
    aeff=irfs['aeff'],
    geom=geom,
)
exposure.slice_by_idx({'energy': 3}).plot(add_cbar=True);
../_images/notebooks_simulate_3d_12_0.png
In [10]:
background = make_map_background_irf(
    pointing=pointing,
    livetime=livetime,
    bkg=irfs['bkg'],
    geom=geom,
)
background.slice_by_idx({'energy': 3}).plot(add_cbar=True);
../_images/notebooks_simulate_3d_13_0.png
In [11]:
psf = irfs['psf'].to_energy_dependent_table_psf(theta=offset)
psf_kernel = PSFKernel.from_table_psf(
    psf,
    geom,
    max_radius=0.3 * u.deg,
)
psf_kernel.psf_kernel_map.sum_over_axes().plot(stretch='log');
WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
../_images/notebooks_simulate_3d_14_1.png
In [12]:
edisp = irfs['edisp'].to_energy_dispersion(offset=offset)
edisp.plot_matrix();
/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.6/site-packages/astropy/units/quantity.py:639: RuntimeWarning: invalid value encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
../_images/notebooks_simulate_3d_15_1.png
In [13]:
%%time
# The idea is that we have this class that can compute `npred`
# maps, i.e. "predicted counts per pixel" given the model and
# the observation infos: exposure, background, PSF and EDISP
evaluator = MapEvaluator(
    model=sky_model,
    exposure=exposure,
    background=background,
    psf=psf_kernel,
)
CPU times: user 9 µs, sys: 0 ns, total: 9 µs
Wall time: 12.9 µs
In [14]:
# Accessing and saving a lot of the following maps is for debugging.
# Just for a simulation one doesn't need to store all these things.
# dnde = evaluator.compute_dnde()
# flux = evaluator.compute_flux()
npred = evaluator.compute_npred()
npred_map = WcsNDMap(geom, npred)
In [15]:
npred_map.sum_over_axes().plot(add_cbar=True);
../_images/notebooks_simulate_3d_18_0.png
In [16]:
# The npred map contains negative values, this is probably a bug in the PSFKernel application
npred[npred<0] = 0
In [17]:
# This one line is the core of how to simulate data when
# using binned simulation / analysis: you Poisson fluctuate
# npred to obtain simulated observed counts.
# Compute counts as a Poisson fluctuation
rng = get_random_state(42)
counts = rng.poisson(npred)
counts_map = WcsNDMap(geom, counts)
In [18]:
counts_map.sum_over_axes().plot();
../_images/notebooks_simulate_3d_21_0.png

Fit

Now let’s analyse the simulated data. Here we just fit it again with the same model we had before, but you could do any analysis you like here, e.g. fit a different model, or do a region-based analysis, …

In [19]:
# Define sky model to fit the data
spatial_model = SkyGaussian(
    lon_0='0 deg',
    lat_0='0 deg',
    sigma='1 deg',
)
spectral_model = PowerLaw(
    index=2,
    amplitude='1e-11 cm-2 s-1 TeV-1',
    reference='1 TeV',
)
model = SkyModel(
    spatial_model=spatial_model,
    spectral_model=spectral_model,
)

model.parameters.set_parameter_errors({
    'lon_0': '0.1 deg',
    'lat_0': '0.1 deg',
    'sigma': '0.1 deg',
    'index': '0.1',
    'amplitude': '1e-12 cm-2 s-1 TeV-1',
})

# model.parameters['sigma'].min = 0
print(model.parameters['sigma'])
print(model)
Parameter(name='sigma', value=1.0, unit='deg', min=nan, max=nan, frozen=False)
SkyModel

spatial_model = SkyGaussian

Parameters:

         name   value   error unit min max frozen
        ----- --------- ----- ---- --- --- ------
        lon_0 0.000e+00   nan  deg nan nan  False
        lat_0 0.000e+00   nan  deg nan nan  False
        sigma 1.000e+00   nan  deg nan nan  False

spectral_model = PowerLaw

Parameters:

           name     value   error       unit         min    max frozen
        --------- --------- ----- --------------- --------- --- ------
            index 2.000e+00   nan                       nan nan  False
        amplitude 1.000e-11   nan 1 / (cm2 s TeV)       nan nan  False
        reference 1.000e+00   nan             TeV 0.000e+00 nan   True

In [20]:
%%time
fit = MapFit(
    model=model,
    counts=counts_map,
    exposure=exposure,
    background=background,
    psf=psf_kernel,
)

fit.fit()

FCN = 267764.03642114025 TOTAL NCALL = 215 NCALLS = 215
EDM = 0.00014033799852801926 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 par_000_lon_0 0.208187 0.00933751 0 0
2 par_001_lat_0 0.0874705 0.00932412 0 0
3 par_002_sigma 0.292286 0.00656484 0 0
4 par_003_index 3.01174 0.0319113 0 0
5 par_004_amplitude 9.62712e-12 5.10933e-13 0 0
6 par_005_reference 1 0 0 0 0.0 FIXED

CPU times: user 2min 27s, sys: 18.3 s, total: 2min 45s
Wall time: 41.8 s
In [21]:
print('True values:\n\n{}\n\n'.format(sky_model.parameters))
print('Fit result:\n\n{}\n\n'.format(model.parameters))
True values:

ParameterList
Parameter(name='lon_0', value=0.2, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='lat_0', value=0.1, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='sigma', value=0.3, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='index', value=3.0, unit='', min=nan, max=nan, frozen=False)
Parameter(name='amplitude', value=1e-11, unit='1 / (cm2 s TeV)', min=nan, max=nan, frozen=False)
Parameter(name='reference', value=1.0, unit='TeV', min=0.0, max=nan, frozen=True)

Covariance:
None


Fit result:

ParameterList
Parameter(name='lon_0', value=0.20818674984823562, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='lat_0', value=0.08747047208677965, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='sigma', value=0.29228569404349286, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='index', value=3.011740728194449, unit='', min=nan, max=nan, frozen=False)
Parameter(name='amplitude', value=9.627120472938576e-12, unit='1 / (cm2 s TeV)', min=nan, max=nan, frozen=False)
Parameter(name='reference', value=1.0, unit='TeV', min=0.0, max=nan, frozen=True)

Covariance:
[[ 8.71890016e-05 -1.59075005e-06  9.96408069e-07 -1.27215882e-06
   1.70076868e-17  0.00000000e+00]
 [-1.59075005e-06  8.69392353e-05 -7.39577817e-07 -4.07548557e-07
  -3.08985589e-17  0.00000000e+00]
 [ 9.96408069e-07 -7.39577817e-07  4.30970756e-05 -4.25486329e-06
   1.07553805e-15  0.00000000e+00]
 [-1.27215882e-06 -4.07548557e-07 -4.25486329e-06  1.01833295e-03
  -1.36070442e-14  0.00000000e+00]
 [ 1.70076868e-17 -3.08985589e-17  1.07553805e-15 -1.36070442e-14
   2.61052580e-25  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]]


In [22]:
# TODO: show e.g. how to make a residual image

iminuit

What we have done for now is to write a very thin wrapper for http://iminuit.readthedocs.io/ as a fitting backend. This is just a prototype, we will improve this interface and add other fitting backends (e.g. Sherpa or scipy.optimize or emcee or …)

As a power-user, you can access fit.iminuit and get the full power of what is developed there already. E.g. the fit.fit() call ran Minuit.migrad() and Minuit.hesse() in the background, and you have access to e.g. the covariance matrix, or can check a likelihood profile, or can run Minuit.minos() to compute asymmetric errors or …

In [23]:
# Check correlation between model parameters
# As expected in this simple case,
# spatial parameters are uncorrelated,
# but the spectral model amplitude and index are correlated as always
fit.minuit.print_matrix()
+
par_000_lon_0
par_001_lat_0
par_002_sigma
par_003_index
par_004_amplitude
par_000_lon_0 1.00 -0.02 0.02 -0.00 0.00
par_001_lat_0 -0.02 1.00 -0.01 -0.00 -0.01
par_002_sigma 0.02 -0.01 1.00 -0.02 0.32
par_003_index -0.00 -0.00 -0.02 1.00 -0.83
par_004_amplitude 0.00 -0.01 0.32 -0.83 1.00
In [24]:
# You can use likelihood profiles to check if your model is
# well constrained or not, and if the fit really converged
fit.minuit.draw_profile('par_002_sigma');
../_images/notebooks_simulate_3d_29_0.png