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.

[1]:

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

[2]:

from astropy.io import fits
from astropy.coordinates import Angle, SkyCoord
from gammapy.maps import Map
from gammapy.detect import CWTKernels, CWT, CWTData

CWT Algorithm

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

[3]:

counts = Map.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz")

background = Map.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-background.fits.gz"
)

maps = {"counts": counts, "background": background}

[4]:

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

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

ax = fig.add_subplot(122, projection=maps["background"].geom.wcs)
maps["background"].plot(vmax=8, ax=ax);
../_images/notebooks_cwt_9_0.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.

[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.

[6]:

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

[6]:
Table length=13
NameSource
str46str32
Number of scales2
Minimal scale6.0
Step scale1.3
Scales[6. 7.8]
Kernels approx width83
Kernels approx sum0.9999153345163221
Kernels approx max0.0015479047097235869
Kernels base width for 6.0 scale63
Kernels base sum for 6.0 scale3.0166320971388447e-05
Kernels base max for 6.0 scale8.599470609575481e-05
Kernels base width for 7.800000000000001 scale83
Kernels base sum for 7.800000000000001 scale1.456141213300368e-05
Kernels base max for 7.800000000000001 scale3.0109136968507695e-05

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

[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.

[8]:

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.

[9]:

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

[10]:

# Start the algorithm
cwt.analyze(cwt_data)

Results of analysis

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

[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?

plt.subplots_adjust(hspace=0.4)
../_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.

[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.

[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?

plt.subplots_adjust(hspace=0.4)
../_images/notebooks_cwt_25_0.png

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

[14]:

data_iter.image_info(name="approx_bkg")

[14]:
Table length=7
MetricsSource
str10str20
Nameapprox_bkg
Shape2D image
Variance0.06792041028304846
Mean0.3450471835001858
Max value1.3431556509385707
Min value0.025070552405371926
Sum values27603.774680014867

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

[15]:

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

[15]:
Table masked=True length=14
Scale powerMetricsSource
int64str10str19
1Namesupport
--Shape2D image
--Variance0.08244281109375
--Mean0.0906625
--Max valueTrue
--Min valueFalse
--Sum values7253
2Namesupport
--Shape2D image
--Variance0.10431465234374998
--Mean0.1183125
--Max valueTrue
--Min valueFalse
--Sum values9465
[16]:

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

[16]:
Table length=7
MetricsSource
str10str19
Namesupport
Shape3D cube
Variance0.09356986234375002
Mean0.1044875
Max valueTrue
Min valueFalse
Sum values16718

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

[17]:

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

[18]:

difference.cube_info("support")

[18]:
Table length=7
MetricsSource
str10str19
Namesupport
Shape3D cube
Variance0.09356986234375002
Mean0.1044875
Max valueTrue
Min valueFalse
Sum values16718
[19]:

difference.info_table.show_in_notebook()

[19]:
Table length=13
idxNameShapeVarianceMeanMax valueMin valueSum values
0counts2D image0.68441189750.4085539.00.032684.0
1background2D image0.088803403195964910.356857903749868252.45129036903381350.0762577354907989528548.632299989462
2model2D image8.921001389341455e-060.00056713013956863420.056259393496281560.045.37041116549074
3approx2D image0.0138816258886271530.050778161816230560.935265850323209-0.220590620943778044062.2529452984445
4approx_bkg2D image0.067920410283048460.34504718350018581.34315565093857070.02507055240537192627603.774680014867
5transform_2d2D image8.921001389341455e-060.00056713013956863420.056259393496281560.045.37041116549074
6model_plus_approx2D image0.0143920039096260640.051345291955799190.9915252438194906-0.220590620943778044107.623356463935
7residual2D image0.62161152639357820.357204708044200838.02873658648009-0.923947594132767128576.376643536063
8maximal2D image3.9089915789042415e-060.000384855301411285340.0403273033377930660.030.788424112902828
9support_2d2D image0.111234437343749980.12748751.00.010199.0
10transform_3d3D cube3.3095296022395382e-061.3167352804236912e-050.040327303337793066-0.0097946299032343162.106776448677906
11error3D cube3.499614659470687e-080.000338744908420145240.00114967824192518316.1270922334229e-0554.199185347223235
12support3D cube0.093569862343750020.10448751.00.016718.0
[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?

    plt.subplots_adjust(hspace=0.4)
../_images/notebooks_cwt_35_0.png

You can save the results if you want

[21]:

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