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

CWT Algorithm

This notebook tutorial shows how to work with CWT algorithm for detecting gamma-ray sources.

You can find the docs here and source code on GitHub here for better understanding how the algorithm is constructed.

Setup

On this section we just import some packages that can be used (or maybe not) in this tutorial. You can also see the versions of the packages in the outputs below and notice that this notebook was written on Python 2.7. Don’t worry about that because the code is also Python 3 compatible.

In [1]:

# Render our plots inline
%matplotlib inline
import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = (15, 5)

In [2]:

import sys
import numpy as np
import scipy as sp

print("Python version: " + sys.version)
print("Numpy version: " + np.__version__)
print("Scipy version: " + sp.__version__)

Python version: 3.7.0 | packaged by conda-forge | (default, Nov 12 2018, 11:54:05)
[Clang 9.0.0 (clang-900.0.37)]
Numpy version: 1.15.4
Scipy version: 1.1.0

CWT Algorithm. PlayGround

First of all we import the data which should be analysied.

In [3]:

import os
from astropy.io import fits
from astropy.coordinates import Angle, SkyCoord
from gammapy.maps import Map

filename = "$GAMMAPY_DATA/fermi_survey/all.fits.gz"

counts = Map.read(filename=filename, hdu="COUNTS")
background = Map.read(filename=filename, hdu="BACKGROUND")

width = Angle([20, 10], "deg")
position = counts.geom.center_skydir
counts = counts.cutout(position=position, width=width)
background = background.cutout(position=position, width=width)

data = dict(counts=counts, background=background)

In [4]:

fig = plt.figure(figsize=(15, 3))

ax = fig.add_subplot(121, projection=data["counts"].geom.wcs)
data["counts"].plot(vmax=10, ax=ax, fig=fig)

ax = fig.add_subplot(122, projection=data["background"].geom.wcs)
data["background"].plot(vmax=10, ax=ax, fig=fig)

Out[4]:

(<Figure size 1080x216 with 2 Axes>,
 <matplotlib.axes._subplots.WCSAxesSubplot at 0x1c14c88eb8>,
 None)
../_images/notebooks_cwt_9_1.png

Let’s explore how CWT works. At first define parameters of the algorithm. An imperative parameter is kernels (detect.CWTKernels object). So we should create it.

In [5]:

# Input parameters for CWTKernels
N_SCALE = 2  # Number of scales considered.
MIN_SCALE = 6.0  # First scale used.
STEP_SCALE = 1.3  # Base scaling factor.

In [6]:

from gammapy.detect import CWTKernels

cwt_kernels = CWTKernels(
    n_scale=N_SCALE, min_scale=MIN_SCALE, step_scale=STEP_SCALE
)
print(cwt_kernels.info_table)

                     Name                              Source
---------------------------------------------- ----------------------
                              Number of scales                      2
                                 Minimal scale                    6.0
                                    Step scale                    1.3
                                        Scales              [6.  7.8]
                          Kernels approx width                     83
                            Kernels approx sum     0.9999153345163221
                            Kernels approx max  0.0015479047097235869
              Kernels base width for 6.0 scale                     63
                Kernels base sum for 6.0 scale 3.0166320971388447e-05
                Kernels base max for 6.0 scale  8.599470609575481e-05
Kernels base width for 7.800000000000001 scale                     83
  Kernels base sum for 7.800000000000001 scale  1.456141213300368e-05
  Kernels base max for 7.800000000000001 scale 3.0109136968507695e-05

Other parameters are optional, in this demonstration define them all.

In [7]:

MAX_ITER = 10  # The maximum number of iterations of the CWT algorithm.
TOL = 1e-5  # Tolerance for stopping criterion.
SIGNIFICANCE_THRESHOLD = 2.0  # Measure of statistical significance.
SIGNIFICANCE_ISLAND_THRESHOLD = (
    None
)  # Measure is used for cleaning of isolated pixel islands.
REMOVE_ISOLATED = True  # If True, isolated pixels will be removed.
KEEP_HISTORY = True  # If you want to save images of all the iterations

Let’s start to analyse input data. Import Logging module to see how the algorithm works during data analysis.

In [8]:

from gammapy.detect import CWT
import logging

logger = logging.getLogger()
logger.setLevel(logging.INFO)
cwt = CWT(
    kernels=cwt_kernels,
    tol=TOL,
    significance_threshold=SIGNIFICANCE_THRESHOLD,
    significance_island_threshold=SIGNIFICANCE_ISLAND_THRESHOLD,
    remove_isolated=REMOVE_ISOLATED,
    keep_history=KEEP_HISTORY,
)

In order to the algorithm was able to analyze source images, you need to convert them to a special format, i.e. create an CWTData object. Do this.

In [9]:

from gammapy.detect import CWTKernels, CWTData

cwt_data = CWTData(
    counts=data["counts"], background=data["background"], n_scale=N_SCALE
)

In [10]:

# Start the algorithm
cwt.analyze(cwt_data)

Results of analysis

Look at the results of CWT algorithm. Print all the images.

In [11]:

PLOT_VALUE_MAX = 5
FIG_SIZE = (15, 35)

fig = plt.figure(figsize=FIG_SIZE)
images = cwt_data.images()
for index, (name, image) in enumerate(images.items()):
    ax = fig.add_subplot(len(images), 2, index + 1, projection=image.geom.wcs)
    image.plot(vmax=PLOT_VALUE_MAX, fig=fig, ax=ax)
    plt.title(name)  # Maybe add a Name in SkyImage.plots?
../_images/notebooks_cwt_21_0.png

As you can see in the implementation of CWT above, it has the parameter keep_history. If you set to it True-value, it means that CWT would save all the images from iterations. Algorithm keeps images of only last CWT start. Let’s do this in the demonstration.

In [12]:

history = cwt.history
print(
    "Number of iterations: {0}".format(len(history) - 1)
)  # -1 because CWT save start images too

Number of iterations: 10

Let’s have a look, what’s happening with images after the first iteration.

In [13]:

N_ITER = 1
assert 0 < N_ITER < len(history)
data_iter = history[N_ITER]

fig = plt.figure(figsize=FIG_SIZE)
images_iter = data_iter.images()
for index, (name, image) in enumerate(images_iter.items()):
    ax = fig.add_subplot(
        len(images_iter), 2, index + 1, projection=image.geom.wcs
    )
    image.plot(vmax=PLOT_VALUE_MAX, fig=fig, ax=ax)
    plt.title(name)  # Maybe add a Name in SkyImage.plots?
../_images/notebooks_cwt_25_0.png

You can get the information about the one particular image in that way:

In [14]:

print(data_iter.image_info(name="approx_bkg"))

 Metrics          Source
---------- --------------------
      Name           approx_bkg
     Shape             2D image
  Variance   0.0800090238572501
      Mean   0.4230576474671889
 Max value    1.216320738970982
 Min value 0.032279497204102725
Sum values    8461.152949343777

You can also get the information about cubes. Or information about all the data.

In [15]:

print(data_iter.cube_info(name="support", per_scale=True))

Scale power  Metrics          Source
----------- ---------- -------------------
          1       Name             support
         --      Shape            2D image
         --   Variance 0.10030838999999998
         --       Mean              0.1131
         --  Max value                True
         --  Min value               False
         -- Sum values                2262
          2       Name             support
         --      Shape            2D image
         --   Variance 0.12255100000000002
         --       Mean               0.143
         --  Max value                True
         --  Min value               False
         -- Sum values                2860

In [16]:

print(data_iter.cube_info(name="support", per_scale=False))

 Metrics          Source
---------- -------------------
      Name             support
     Shape             3D cube
  Variance 0.11165319750000001
      Mean             0.12805
 Max value                True
 Min value               False
Sum values                5122

Also you can see the difference betwen the iterations in that way:

In [17]:

history = cwt.history  # get list of 'CWTData' objects
difference = (
    history[1] - history[0]
)  # get new `CWTData` obj, let's work with them

In [18]:

print(difference.cube_info("support"))

 Metrics          Source
---------- -------------------
      Name             support
     Shape             3D cube
  Variance 0.11165319750000001
      Mean             0.12805
 Max value                True
 Min value               False
Sum values                5122

In [19]:

difference.info_table.show_in_notebook()

Out[19]:
Table length=13
idxNameShapeVarianceMeanMax valueMin valueSum values
0counts2D image0.66721419750.4444524.00.08889.0
1background2D image0.16522025351138420.453299561273412373.91485793084037150.093821701373177619065.991225468248
2model2D image2.4686984587806907e-060.0004556764061886410.0154327993376764880.09.11352812377282
3approx2D image0.004022666598226336-0.008046110167776410.2713385823751633-0.2517653841639436-160.92220335552818
4approx_bkg2D image0.08000902385725010.42305764746718891.2163207389709820.0322794972041027258461.152949343777
5transform_2d2D image2.4686984587806907e-060.0004556764061886410.0154327993376764880.09.11352812377282
6model_plus_approx2D image0.004115467027140748-0.0075904337615877690.28513324259070433-0.2517653841639436-151.8086752317554
7residual2D image0.666300092318290.452040433761587823.943554761680172-0.283584814746093769040.808675231756
8maximal2D image1.2706495970464085e-060.00034396466393524450.0111684652858454630.06.87929327870489
9support_2d2D image0.139443750000000030.16751.00.03350.0
10transform_3d3D cube1.5198680587581294e-065.318961505149442e-060.011168465285845463-0.008590491273854370.2127584602059777
11error3D cube4.4263978237355114e-080.000377883070469780.00119109721270137926.826532887607761e-0515.1153228187912
12support3D cube0.111653197500000010.128051.00.05122.0
In [20]:

fig = plt.figure(figsize=FIG_SIZE)
images_diff = difference.images()
for index, (name, image) in enumerate(images_diff.items()):
    ax = fig.add_subplot(
        len(images_diff), 2, index + 1, projection=image.geom.wcs
    )
    image.plot(vmax=PLOT_VALUE_MAX, fig=fig, ax=ax)
    plt.title(name)  # Maybe add a Name in SkyImage.plots?
../_images/notebooks_cwt_35_0.png

You can save the results if you want

In [21]:

# cwt_data.write('test-cwt.fits', True)