DC2 Image Analysis, 511 keV, Image Deconvolution

updated on 2024-01-30 (the commit 26cfdeacb25335bd511a91c4f8a29bdeb36408f2)

This notebook focuses on the image deconvolution with the spacecraft attitude (scatt) binning method for DC2. Using the 511 keV thin disk 3-month simulation data created for DC2, an example of the image analysis will be presented. If you have not run through 511keV-DC2-ScAtt-DataReduction.ipynb, please see it first.

[1]:

from histpy import Histogram, HealpixAxis, Axis, Axes
from mhealpy import HealpixMap
from astropy.coordinates import SkyCoord, cartesian_to_spherical, Galactic

from cosipy.response import FullDetectorResponse
from cosipy.spacecraftfile import SpacecraftFile
from cosipy.ts_map.TSMap import TSMap
from cosipy.data_io import UnBinnedData, BinnedData
from cosipy.image_deconvolution import SpacecraftAttitudeExposureTable, CoordsysConversionMatrix, DataLoader, ImageDeconvolution

# cosipy uses astropy units
import astropy.units as u
from astropy.units import Quantity
from astropy.coordinates import SkyCoord
from astropy.time import Time
from astropy.table import Table
from astropy.io import fits
from scoords import Attitude, SpacecraftFrame

#3ML is needed for spectral modeling
from threeML import *
from astromodels import Band

#Other standard libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

import healpy as hp
from tqdm.autonotebook import tqdm

%matplotlib inline


WARNING: version mismatch between CFITSIO header (v4) and linked library (v4.01).

Welcome to JupyROOT 6.24/06

19:41:37 WARNING   The naima package is not available. Models that depend on it will not be         functions.py:48
                  available

         WARNING   The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it  functions.py:69
                  will not be available.

19:41:38 WARNING   The ebltable package is not available. Models that depend on it will not be     absorption.py:36
                  available

         WARNING   We have set the min_value of K to 1e-99 because there was a postive transform   parameter.py:704

         WARNING   We have set the min_value of K to 1e-99 because there was a postive transform   parameter.py:704

         WARNING   We have set the min_value of K to 1e-99 because there was a postive transform   parameter.py:704

         WARNING   We have set the min_value of K to 1e-99 because there was a postive transform   parameter.py:704

         WARNING   We have set the min_value of F to 1e-99 because there was a postive transform   parameter.py:704

         WARNING   We have set the min_value of K to 1e-99 because there was a postive transform   parameter.py:704

19:41:38 INFO      Starting 3ML!                                                                     __init__.py:35

         WARNING   WARNINGs here are NOT errors                                                      __init__.py:36

         WARNING   but are inform you about optional packages that can be installed                  __init__.py:37

         WARNING    to disable these messages, turn off start_warning in your config file            __init__.py:40

         WARNING   Multinest minimizer not available                                           minimization.py:1357

         WARNING   PyGMO is not available                                                      minimization.py:1369

         WARNING   The cthreeML package is not installed. You will not be able to use plugins which  __init__.py:94
                  require the C/C++ interface (currently HAWC)

         WARNING   Could not import plugin HAWCLike.py. Do you have the relative instrument         __init__.py:144
                  software installed and configured?

         WARNING   Could not import plugin FermiLATLike.py. Do you have the relative instrument     __init__.py:144
                  software installed and configured?

         WARNING   No fermitools installed                                              lat_transient_builder.py:44

         WARNING   Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal         __init__.py:387
                  performances in 3ML

         WARNING   Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal         __init__.py:387
                  performances in 3ML

         WARNING   Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal     __init__.py:387
                  performances in 3ML

0. Files needed for this notebook

From wasabi - cosi-pipeline-public/COSI-SMEX/DC2/Responses/SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5

From docs/tutorials/image_deconvolution/511keV/ScAttBinning - inputs_511keV_DC2.yaml - imagedeconvolution_parfile_scatt_511keV.yml

As outputs from the notebook 511keV-DC2-ScAtt-DataReduction.ipynb - 511keV_scatt_binning_DC2_bkg.hdf5 - 511keV_scatt_binning_DC2_event.hdf5 - ccm.hdf5

1. Read the response matrix

please modify “path_data” corresponding to your environment.

[2]:

path_data = "path/to/data/"

[3]:

response_path = path_data + "SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5"

response = FullDetectorResponse.open(response_path)

[4]:

response

[4]:

FILENAME: '/Users/yoneda/Work/Exp/COSI/cosipy-2/data_challenge/DC2/prework/data/Responses/SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5'
AXES:
  NuLambda:
    DESCRIPTION: 'Location of the simulated source in the spacecraft coordinates'
    TYPE: 'healpix'
    NPIX: 3072
    NSIDE: 16
    SCHEME: 'RING'
  Ei:
    DESCRIPTION: 'Initial simulated energy'
    TYPE: 'log'
    UNIT: 'keV'
    NBINS: 1
    EDGES: [509.0 keV, 513.0 keV]
  Em:
    DESCRIPTION: 'Measured energy'
    TYPE: 'log'
    UNIT: 'keV'
    NBINS: 1
    EDGES: [509.0 keV, 513.0 keV]
  Phi:
    DESCRIPTION: 'Compton angle'
    TYPE: 'linear'
    UNIT: 'deg'
    NBINS: 60
    EDGES: [0.0 deg, 3.0 deg, 6.0 deg, 9.0 deg, 12.0 deg, 15.0 deg, 18.0 deg, 21.0 deg, 24.0 deg, 27.0 deg, 30.0 deg, 33.0 deg, 36.0 deg, 39.0 deg, 42.0 deg, 45.0 deg, 48.0 deg, 51.0 deg, 54.0 deg, 57.0 deg, 60.0 deg, 63.0 deg, 66.0 deg, 69.0 deg, 72.0 deg, 75.0 deg, 78.0 deg, 81.0 deg, 84.0 deg, 87.0 deg, 90.0 deg, 93.0 deg, 96.0 deg, 99.0 deg, 102.0 deg, 105.0 deg, 108.0 deg, 111.0 deg, 114.0 deg, 117.0 deg, 120.0 deg, 123.0 deg, 126.0 deg, 129.0 deg, 132.0 deg, 135.0 deg, 138.0 deg, 141.0 deg, 144.0 deg, 147.0 deg, 150.0 deg, 153.0 deg, 156.0 deg, 159.0 deg, 162.0 deg, 165.0 deg, 168.0 deg, 171.0 deg, 174.0 deg, 177.0 deg, 180.0 deg]
  PsiChi:
    DESCRIPTION: 'Location in the Compton Data Space'
    TYPE: 'healpix'
    NPIX: 3072
    NSIDE: 16
    SCHEME: 'RING'

2. Read binned 511keV binned files (source and background)

[5]:

%%time

#  background
data_bkg = BinnedData("inputs_511keV_DC2.yaml")
data_bkg.load_binned_data_from_hdf5("511keV_scatt_binning_DC2_bkg.hdf5")

##  signal + background
data_511keV = BinnedData("inputs_511keV_DC2.yaml")
data_511keV.load_binned_data_from_hdf5("511keV_scatt_binning_DC2_event.hdf5")

CPU times: user 149 ms, sys: 806 ms, total: 955 ms
Wall time: 958 ms

3. Load the coordsys conversion matrix

[6]:

%%time

ccm = CoordsysConversionMatrix.open("ccm.hdf5")

CPU times: user 1.63 s, sys: 77.8 ms, total: 1.71 s
Wall time: 1.72 s

4. Imaging deconvolution

Brief overview of the image deconvolution

Basically, we have to maximize the following likelihood function

\[\log L = \sum_i X_i \log \epsilon_i - \sum_i \epsilon_i\]

\(X_i\): detected counts at \(i\)-th bin ( \(i\) : index of the Compton Data Space)

\(\epsilon_i = \sum_j R_{ij} \lambda_j + b_i\) : expected counts ( \(j\) : index of the model space)

\(\lambda_j\) : the model map (basically gamma-ray flux at \(j\)-th pixel)

\(b_i\) : the background at \(i\)-th bin

\(R_{ij}\) : the response matrix

Since we have to optimize the flux in each pixel, and the number of parameters is large, we adopt an iterative approach to find a solution of the above equation. The simplest one is the ML-EM (Maximum Likelihood Expectation Maximization) algorithm. It is also known as the Richardson-Lucy algorithm.

\[\lambda_{j}^{k+1} = \lambda_{j}^{k} + \delta \lambda_{j}^{k}\]

\[\delta \lambda_{j}^{k} = \frac{\lambda_{j}^{k}}{\sum_{i} R_{ij}} \sum_{i} \left(\frac{ X_{i} }{\epsilon_{i}} - 1 \right) R_{ij}\]

We refer to \(\delta \lambda_{j}^{k}\) as the delta map.

As for now, the two improved algorithms are implemented in COSIpy.

Accelerated ML-EM algorithm (Knoedlseder+99)

\[\lambda_{j}^{k+1} = \lambda_{j}^{k} + \alpha^{k} \delta \lambda_{j}^{k}\]

\[\alpha^{k} < \mathrm{max}(- \lambda_{j}^{k} / \delta \lambda_{j}^{k})\]

Practically, in order not to accelerate the algorithm excessively, we set the maximum value of \(\alpha\) (\(\alpha_{\mathrm{max}}\)). Then, \(\alpha\) is calculated as:

\[\alpha^{k} = \mathrm{min}(\mathrm{max}(- \lambda_{j}^{k} / \delta \lambda_{j}^{k}), \alpha_{\mathrm{max}})\]

Noise damping using gaussian smoothing (Knoedlseder+05, Siegert+20)

\[\lambda_{j}^{k+1} = \lambda_{j}^{k} + \alpha^{k} \left[ w_j \delta \lambda_{j}^{k} \right]_{\mathrm{gauss}}\]

\[w_j = \left(\sum_{i} R_{ij}\right)^\beta\]

\(\left[ ... \right]_{\mathrm{gauss}}\) means that the differential image is smoothed by a gaussian filter.

4-1. Prepare DataLoader containing all neccesary datasets

[7]:

dataloader = DataLoader.load(data_511keV.binned_data,
                             data_bkg.binned_data,
                             response,
                             ccm,
                             is_miniDC2_format = False)

[8]:

dataloader._modify_axes()


WARNING FutureWarning: Note that _modify_axes() in DataLoader was implemented for a temporary use. It will be removed in the future.


WARNING UserWarning: Make sure to perform _modify_axes() only once after the data are loaded.

... checking the axis ScAtt of the event and background files...
    --> pass (edges)
    --> pass (unit)
... checking the axis Em of the event and background files...
    --> pass (edges)
    --> pass (unit)
... checking the axis Phi of the event and background files...
    --> pass (edges)
    --> pass (unit)
... checking the axis PsiChi of the event and background files...
    --> pass (edges)
    --> pass (unit)
...checking the axis Em of the event and response files...
    --> pass (edges)
...checking the axis Phi of the event and response files...
    --> pass (edges)
...checking the axis PsiChi of the event and response files...
    --> pass (edges)
The axes in the event and background files are redefined. Now they are consistent with those of the response file.

(In the future, we plan to remove the method “_modify_axes.”)

4-2. Load the response file

The response file will be loaded on the CPU memory. It requires a few GB. In the actual COSI satellite analysis, the response could be much larger, perhaps ~1TB wiht finer bin size.

So loading it on the memory might be unrealistic in the future. The optimized (lazy) loading would be a next work.

[9]:

%%time

dataloader.load_full_detector_response_on_memory()

CPU times: user 13.5 s, sys: 1.48 s, total: 15 s
Wall time: 15 s

Here, we calculate an auxiliary matrix for the deconvolution. Basically, it is a response matrix projected into the model space (\(\sum_{i} R_{ij}\)). Currently, it is mandatory to run this command for the image deconvolution.

[10]:

dataloader.calc_image_response_projected()

... (DataLoader) calculating a projected image response ...

4-3. Initialize the instance of the image deconvolution class

First, we prepare an instance of the ImageDeconvolution class and then register the dataset and parameters for the deconvolution. After that, you can start the calculation.

please modify this parameter_filepath corresponding to your environment.

[11]:

parameter_filepath = "imagedeconvolution_parfile_scatt_511keV.yml"

[12]:

image_deconvolution = ImageDeconvolution()

# set dataloader to image_deconvolution
image_deconvolution.set_data(dataloader)

# set a parameter file for the image deconvolution
image_deconvolution.read_parameterfile(parameter_filepath)

data for image deconvolution was set ->  <cosipy.image_deconvolution.data_loader.DataLoader object at 0x2b65478e0>
parameter file for image deconvolution was set ->  imagedeconvolution_parfile_scatt_511keV.yml

Initialize image_deconvolution

In this process, a model map is defined following the input parameters, and it is initialized. Also, it prepares ancillary data for the image deconvolution, e.g., the expected counts with the initial model map, gaussian smoothing filter etc.

I describe parameters in the parameter file.

model_property

Name	Unit	Description	Notes
coordinate	str	the coordinate system of the model map	As for now, it must be ‘galactic’
nside	int	NSIDE of the model map	it must be the same as NSIDE of ‘lb’ axis of the coordinate conversion matrix
scheme	str	SCHEME of the model map	As for now, it must be ‘ring’
energy_edges	list of float [keV]	The definition of the energy bins of the model map	As for now, it must be the same as that of the response matrix

model_initialization

Name	Unit	Description	Notes
algorithm	str	the method name to initialize the model map	As for now, only ‘flat’ can be used
parameter_flat:values	list of float [cm-2 s-1 sr-1]	the list of photon fluxes for each energy band	the length of the list should be the same as the length of “energy_edges” - 1

deconvolution

Name	Unit	Description	Notes
algorithm	str	the name of the image deconvolution algorithm	As for now, only ‘RL’ is supported

parameter_RL:iteration	int	The maximum number of the iteration
parameter_RL:acceleration	bool	whether the accelerated ML-EM algorithm (Knoedlseder+99) is used
parameter_RL:alpha_max	float	the maximum value for the acceleration parameter
parameter_RL:save_results_each_iteration	bool	whether an updated model map, detal map, likelihood etc. are saved at the end of each iteration
parameter_RL:response_weighting	bool	whether a delta map is renormalized based on the exposure time on each pixel, namely \(w_j = (\sum_{i} R_{ij})^{\beta}\) (see Knoedlseder+05, Siegert+20)
parameter_RL:response_weighting_index	float	\(\beta\) in the above equation
parameter_RL:smoothing	bool	whether a Gaussian filter is used (see Knoedlseder+05, Siegert+20)
parameter_RL:smoothing_FWHM	float, degree	the FWHM of the Gaussian in the filter
parameter_RL:background_normalization_fitting	bool	whether the background normalization factor is optimized at each iteration	As for now, the single background normalization factor is used in all of the bins
parameter_RL:background_normalization_range	list of float	the range of the normalization factor	should be positive

[13]:

image_deconvolution.initialize()

#### Initialization ####
1. generating a model map
---- parameters ----
coordinate: galactic
energy_edges:
- 509.0
- 513.0
nside: 16
scheme: ring

2. initializing the model map ...
---- parameters ----
algorithm: flat
parameter_flat:
  values:
  - 1e-4

3. registering the deconvolution algorithm ...
... calculating the expected events with the initial model map ...
... calculating the response weighting filter...
... calculating the gaussian filter...

---- parameters ----
algorithm: RL
parameter_RL:
  acceleration: true
  alpha_max: 10.0
  background_normalization_fitting: false
  background_normalization_range:
  - 0.01
  - 10.0
  iteration: 10
  response_weighting: true
  response_weighting_index: 0.5
  save_results_each_iteration: false
  smoothing: true
  smoothing_FWHM: 2.0
  smoothing_max_sigma: 10.0

#### Done ####

(You can change the parameters as follows)

Note that when you modify the parameters, do not forget to run “initialize” again!

[14]:

image_deconvolution.override_parameter("deconvolution:parameter_RL:iteration = 30")
image_deconvolution.override_parameter("deconvolution:parameter_RL:background_normalization_fitting = True")
image_deconvolution.override_parameter("deconvolution:parameter_RL:alpha_max = 10")
image_deconvolution.override_parameter("deconvolution:parameter_RL:smoothing_FWHM = 3.0")

image_deconvolution.initialize()

#### Initialization ####
1. generating a model map
---- parameters ----
coordinate: galactic
energy_edges:
- 509.0
- 513.0
nside: 16
scheme: ring

2. initializing the model map ...
---- parameters ----
algorithm: flat
parameter_flat:
  values:
  - 1e-4

3. registering the deconvolution algorithm ...
... calculating the expected events with the initial model map ...
... calculating the response weighting filter...
... calculating the gaussian filter...

---- parameters ----
algorithm: RL
parameter_RL:
  acceleration: true
  alpha_max: 10
  background_normalization_fitting: true
  background_normalization_range:
  - 0.01
  - 10.0
  iteration: 30
  response_weighting: true
  response_weighting_index: 0.5
  save_results_each_iteration: false
  smoothing: true
  smoothing_FWHM: 3.0
  smoothing_max_sigma: 10.0

#### Done ####

4-5. Start the image deconvolution

With MacBook Pro with M1 Max and 64 GB memory, it takes about 6 minutes for 30 iterations.

[15]:

%%time

all_results = image_deconvolution.run_deconvolution()

#### Deconvolution Starts ####

  Iteration 1/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step


WARNING RuntimeWarning: invalid value encountered in divide

--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 2.5039460293538105
    loglikelihood: -1563364.0277526558
    background_normalization: 1.0048700233481955
  Iteration 2/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1519357.4702937151
    background_normalization: 0.9944142064277177
  Iteration 3/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 2.4670760938135765
    loglikelihood: -1499867.5506138196
    background_normalization: 0.999275691887223
  Iteration 4/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1496920.474980764
    background_normalization: 1.0004892236020582
  Iteration 5/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 4.546207203696152
    loglikelihood: -1490909.3204384344
    background_normalization: 0.9998689870447892
  Iteration 6/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1490365.0435509903
    background_normalization: 0.9995258381190871
  Iteration 7/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 2.854692362839457
    loglikelihood: -1489307.7214813665
    background_normalization: 0.9997388449308033
  Iteration 8/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1489058.487761884
    background_normalization: 0.9998124108372027
  Iteration 9/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 2.2283812062183777
    loglikelihood: -1488593.384611151
    background_normalization: 0.9997745302553334
  Iteration 10/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1488431.7662045578
    background_normalization: 0.999764145247152
  Iteration 11/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.2705691250781777
    loglikelihood: -1488249.8944230902
    background_normalization: 0.9997686118565604
  Iteration 12/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1488124.8362333952
    background_normalization: 0.9997696073518008
  Iteration 13/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1488012.3045336306
    background_normalization: 0.9997697876916849
  Iteration 14/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.2230022900243092
    loglikelihood: -1487888.5786615435
    background_normalization: 0.9997700189565418
  Iteration 15/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1487798.520730182
    background_normalization: 0.9997703163980328
  Iteration 16/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.8098974214440344
    loglikelihood: -1487652.4673017936
    background_normalization: 0.999770562358397
  Iteration 17/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1487583.1765680623
    background_normalization: 0.999771015377049
  Iteration 18/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 5.618448128695514
    loglikelihood: -1487253.9933618857
    background_normalization: 0.9997712604348642
  Iteration 19/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1487215.6271944637
    background_normalization: 0.9997727102976933
  Iteration 20/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 4.541756790045432
    loglikelihood: -1487060.3103523117
    background_normalization: 0.999772909371205
  Iteration 21/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1487034.1422262355
    background_normalization: 0.9997739086816684
  Iteration 22/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1487009.4576823683
    background_normalization: 0.9997741100024176
  Iteration 23/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1486986.1429297891
    background_normalization: 0.9997742950835806
  Iteration 24/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1486964.0949290707
    background_normalization: 0.9997744737637747
  Iteration 25/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 4.967042011343623
    loglikelihood: -1486864.6152302232
    background_normalization: 0.9997746469610125
  Iteration 26/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1486848.7990036565
    background_normalization: 0.9997754874062363
  Iteration 27/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.6330713142086324
    loglikelihood: -1486824.3117899313
    background_normalization: 0.9997756319817698
  Iteration 28/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1486810.3517182083
    background_normalization: 0.9997758601840663
  Iteration 29/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> continue
--> registering results
--> showing results
    alpha: 1.1820824055223498
    loglikelihood: -1486794.6074471278
    background_normalization: 0.9997759950994234
  Iteration 30/30
--> pre-processing
--> E-step
... skip E-step ...
--> M-step
--> post-processing
... calculating the expected events with the updated model map ...
--> checking stopping criteria
--> --> stop
--> registering results
--> showing results
    alpha: 1.0
    loglikelihood: -1486781.9555685548
    background_normalization: 0.9997761488793757
#### Done ####

CPU times: user 33min 1s, sys: 3min 35s, total: 36min 37s
Wall time: 5min 41s

[16]:

import pprint

pprint.pprint(all_results)

[{'alpha': <Quantity 2.50394603>,
  'background_normalization': 1.0048700233481955,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8a2f0>,
  'iteration': 1,
  'loglikelihood': -1563364.0277526558,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8a380>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b89f00>},
 {'alpha': 1.0,
  'background_normalization': 0.9944142064277177,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x10368e1d0>,
  'iteration': 2,
  'loglikelihood': -1519357.4702937151,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x2b657bc70>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b88850>},
 {'alpha': <Quantity 2.46707609>,
  'background_normalization': 0.999275691887223,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8bd60>,
  'iteration': 3,
  'loglikelihood': -1499867.5506138196,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8ba00>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b89990>},
 {'alpha': 1.0,
  'background_normalization': 1.0004892236020582,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8ad70>,
  'iteration': 4,
  'loglikelihood': -1496920.474980764,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x2b6568820>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b88430>},
 {'alpha': <Quantity 4.5462072>,
  'background_normalization': 0.9998689870447892,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8afb0>,
  'iteration': 5,
  'loglikelihood': -1490909.3204384344,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8afe0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b89780>},
 {'alpha': 1.0,
  'background_normalization': 0.9995258381190871,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8bf40>,
  'iteration': 6,
  'loglikelihood': -1490365.0435509903,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8b970>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b88880>},
 {'alpha': <Quantity 2.85469236>,
  'background_normalization': 0.9997388449308033,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8a410>,
  'iteration': 7,
  'loglikelihood': -1489307.7214813665,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8ac50>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8ae30>},
 {'alpha': 1.0,
  'background_normalization': 0.9998124108372027,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07790>,
  'iteration': 8,
  'loglikelihood': -1489058.487761884,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e076d0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07f10>},
 {'alpha': <Quantity 2.22838121>,
  'background_normalization': 0.9997745302553334,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b89600>,
  'iteration': 9,
  'loglikelihood': -1488593.384611151,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8ae60>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b88340>},
 {'alpha': 1.0,
  'background_normalization': 0.999764145247152,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b896f0>,
  'iteration': 10,
  'loglikelihood': -1488431.7662045578,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8ad10>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8b0d0>},
 {'alpha': <Quantity 1.27056913>,
  'background_normalization': 0.9997686118565604,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07bb0>,
  'iteration': 11,
  'loglikelihood': -1488249.8944230902,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e05a50>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07670>},
 {'alpha': 1.0,
  'background_normalization': 0.9997696073518008,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b88640>,
  'iteration': 12,
  'loglikelihood': -1488124.8362333952,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b88cd0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b882b0>},
 {'alpha': 1.0,
  'background_normalization': 0.9997697876916849,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8b0a0>,
  'iteration': 13,
  'loglikelihood': -1488012.3045336306,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8a4a0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b89840>},
 {'alpha': <Quantity 1.22300229>,
  'background_normalization': 0.9997700189565418,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8bc10>,
  'iteration': 14,
  'loglikelihood': -1487888.5786615435,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8aad0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b89540>},
 {'alpha': 1.0,
  'background_normalization': 0.9997703163980328,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e075e0>,
  'iteration': 15,
  'loglikelihood': -1487798.520730182,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e05ae0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x10367dbd0>},
 {'alpha': <Quantity 1.80989742>,
  'background_normalization': 0.999770562358397,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b88d30>,
  'iteration': 16,
  'loglikelihood': -1487652.4673017936,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8a2c0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07880>},
 {'alpha': 1.0,
  'background_normalization': 0.999771015377049,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07df0>,
  'iteration': 17,
  'loglikelihood': -1487583.1765680623,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07b20>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x2b6569b10>},
 {'alpha': <Quantity 5.61844813>,
  'background_normalization': 0.9997712604348642,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b8ace0>,
  'iteration': 18,
  'loglikelihood': -1487253.9933618857,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x2b6579690>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b029b0>},
 {'alpha': 1.0,
  'background_normalization': 0.9997727102976933,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07c40>,
  'iteration': 19,
  'loglikelihood': -1487215.6271944637,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e06950>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b024d0>},
 {'alpha': <Quantity 4.54175679>,
  'background_normalization': 0.999772909371205,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b018d0>,
  'iteration': 20,
  'loglikelihood': -1487060.3103523117,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x2b6569780>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b000d0>},
 {'alpha': 1.0,
  'background_normalization': 0.9997739086816684,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07460>,
  'iteration': 21,
  'loglikelihood': -1487034.1422262355,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x373e07700>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b016f0>},
 {'alpha': 1.0,
  'background_normalization': 0.9997741100024176,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b02bc0>,
  'iteration': 22,
  'loglikelihood': -1487009.4576823683,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b00c40>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b01a80>},
 {'alpha': 1.0,
  'background_normalization': 0.9997742950835806,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b005e0>,
  'iteration': 23,
  'loglikelihood': -1486986.1429297891,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b013f0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b03ac0>},
 {'alpha': 1.0,
  'background_normalization': 0.9997744737637747,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b00f40>,
  'iteration': 24,
  'loglikelihood': -1486964.0949290707,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b01150>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b021d0>},
 {'alpha': <Quantity 4.96704201>,
  'background_normalization': 0.9997746469610125,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b00100>,
  'iteration': 25,
  'loglikelihood': -1486864.6152302232,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b002b0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b020b0>},
 {'alpha': 1.0,
  'background_normalization': 0.9997754874062363,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b03730>,
  'iteration': 26,
  'loglikelihood': -1486848.7990036565,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b03ca0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b02ce0>},
 {'alpha': <Quantity 1.63307131>,
  'background_normalization': 0.9997756319817698,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b02fe0>,
  'iteration': 27,
  'loglikelihood': -1486824.3117899313,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b02200>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b035b0>},
 {'alpha': 1.0,
  'background_normalization': 0.9997758601840663,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b03490>,
  'iteration': 28,
  'loglikelihood': -1486810.3517182083,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b011b0>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b02ad0>},
 {'alpha': <Quantity 1.18208241>,
  'background_normalization': 0.9997759950994234,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b03c10>,
  'iteration': 29,
  'loglikelihood': -1486794.6074471278,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b01030>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b00c70>},
 {'alpha': 1.0,
  'background_normalization': 0.9997761488793757,
  'delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b02e90>,
  'iteration': 30,
  'loglikelihood': -1486781.9555685548,
  'model_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x374b01330>,
  'processed_delta_map': <cosipy.image_deconvolution.modelmap.ModelMap object at 0x2b65790c0>}]

5. Analyze the results

Examples to see/analyze the results are shown below.

Log-likelihood

Plotting the log-likelihood vs the number of iterations

[17]:

x, y = [], []

for result in all_results:
    x.append(result['iteration'])
    y.append(result['loglikelihood'])

plt.plot(x, y)
plt.grid()
plt.xlabel("iteration")
plt.ylabel("loglikelihood")
plt.show()

../../../../_images/tutorials_image_deconvolution_511keV_ScAttBinning_511keV-DC2-ScAtt-ImageDeconvolution_35_0.png

Alpha (the factor used for the acceleration)

Plotting \(\alpha\) vs the number of iterations. \(\alpha\) is a parameter to accelerate the EM algorithm (see the beginning of Section 4). If it is too large, reconstructed images may have artifacts.

[18]:

x, y = [], []

for result in all_results:
    x.append(result['iteration'])
    y.append(result['alpha'])

plt.plot(x, y)
plt.grid()
plt.xlabel("iteration")
plt.ylabel("alpha")
plt.show()

../../../../_images/tutorials_image_deconvolution_511keV_ScAttBinning_511keV-DC2-ScAtt-ImageDeconvolution_37_0.png

Background normalization

Plotting the background nomalization factor vs the number of iterations. If the background model is accurate and the image is reconstructed perfectly, this factor should be close to 1.

[19]:

x, y = [], []

for result in all_results:
    x.append(result['iteration'])
    y.append(result['background_normalization'])

plt.plot(x, y)
plt.grid()
plt.xlabel("iteration")
plt.ylabel("background_normalization")
plt.show()

../../../../_images/tutorials_image_deconvolution_511keV_ScAttBinning_511keV-DC2-ScAtt-ImageDeconvolution_39_0.png

The reconstructed images

[20]:

def plot_reconstructed_image(result, source_position = None): # source_position should be (l,b) in degrees
    iteration = result['iteration']
    image = result['model_map']

    for energy_index in range(image.axes['Ei'].nbins):
        map_healpxmap = HealpixMap(data = image[:,energy_index], unit = image.unit)

        _, ax = map_healpxmap.plot('mollview')

        _.colorbar.set_label(str(image.unit))

        if source_position is not None:
            ax.scatter(source_position[0]*u.deg, source_position[1]*u.deg, transform=ax.get_transform('world'), color = 'red')

        plt.title(label = f"iteration = {iteration}, energy_index = {energy_index} ({image.axes['Ei'].bounds[energy_index][0]}-{image.axes['Ei'].bounds[energy_index][1]})")

Plotting the reconstructed images in all of the energy bands at the 20th iteration

[21]:

iteration = 19

plot_reconstructed_image(all_results[iteration])

../../../../_images/tutorials_image_deconvolution_511keV_ScAttBinning_511keV-DC2-ScAtt-ImageDeconvolution_43_0.png

An example to plot the image in the log scale

[22]:

iteration_idx = 29

result = all_results[iteration_idx]

iteration = result['iteration']
image = result['model_map']

data = image[:,0]
data[data <= 0 * data.unit] = 1e-12 * data.unit

hp.mollview(data, min = 1e-5, norm ='log', unit = str(data.unit), title = f'511 keV image at {iteration}th iteration', cmap = 'magma')

plt.show()

../../../../_images/tutorials_image_deconvolution_511keV_ScAttBinning_511keV-DC2-ScAtt-ImageDeconvolution_45_0.png

[ ]: