Skip to content

Commit

Permalink
minor changes on the example notebooks of the image deconvolution
Browse files Browse the repository at this point in the history
  • Loading branch information
@hiyoneda
hiyoneda committed Jan 30, 2024
1 parent 26cfdea commit ffc19ef
Show file tree
Hide file tree
Showing 10 changed files with 376 additions and 315 deletions.
103 changes: 59 additions & 44 deletions ...rials/image_deconvolution/511keV/GalacticCDS/511keV-DC2-Galactic-ImageDeconvolution.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,12 @@
"source": [
"# DC2 Image Analysis, 511 keV, Image Deconvolution using CDS in the Galactic coordinate system\n",
"\n",
"updated on 2024-01-23 (the commit fe724f0f8d317ae347596708e6c50db6899c7e89)\n",
"updated on 2024-01-30 (the commit 26cfdeacb25335bd511a91c4f8a29bdeb36408f2)\n",
"\n",
"This notebook focuses on the image deconvolution with the Compton data space (CDS) in the Galactic coordinate system.\n",
"Using the 511keV thin disk 3month simulation data created for DC2, an example of the image analysis will be presented."
"An example of the image analysis will be presented using the 511keV thin disk 3-month simulation data created for DC2.\n",
"\n",
"In DC2, we have two options on the coordinate system to describe the Compton scattering direction ($\\chi\\psi$) in the image deconvolution. Please also check the notes written in 511keV-DC2-ScAtt-DataReduction.ipynb."
]
},
{
Expand Down Expand Up @@ -67,9 +69,12 @@
"\n",
"From wasabi\n",
"- cosi-pipeline-public/COSI-SMEX/DC2/Responses/PointSourceReponse/psr_gal_511_DC2.h5.gz (please gunzip it)\n",
" - a pre-computed 511 keV line response file rotated into the Galactic coordinate system\n",
" - a pre-computed 511 keV line response file converted into the Galactic coordinate system\n",
"- cosi-pipeline-public/COSI-SMEX/DC2/Data/Sources/511_thin_disk_3months.fits.gz\n",
"- cosi-pipeline-public/COSI-SMEX/DC2/Data/Backgrounds/albedo_photons_3months_unbinned_data.fits.gz\n",
" - In this notebook, only the albedo gamma-ray background is considered for a tutorial.\n",
" - If you want to consider all of the background components, please replace it with cosi-pipeline-public/COSI-SMEX/DC2/Data/Backgrounds/total_bg_3months_unbinned_data.fits.gz\n",
" - Note that total_bg_3months_unbinned_data.fits.gz is 14.15 GB.\n",
"\n",
"From docs/tutorials/image_deconvolution/511keV/GalacticCDS\n",
"- inputs_511keV_DC2.yaml\n",
Expand Down Expand Up @@ -203,7 +208,7 @@
"id": "4eb8577f-d394-49b9-a13f-a527d4512f77",
"metadata": {},
"source": [
"Convert the data into sparse matrices / add the signal with the background"
"Convert the data into sparse matrices & add the signal to the background"
]
},
{
Expand Down Expand Up @@ -263,7 +268,9 @@
"id": "0e7bb933-0ec0-47af-a18c-ac241abfea82",
"metadata": {},
"source": [
"(In DC2, the number of time bins is 1 because all of events are projected into a signle galactic CDS. In the future, we may also introduce time/scatt binning in the galactic method.)"
"In DC2, the number of time bins should be 1 when you perform the image deconvolution using the galactic CDS.\n",
"It is because the pre-computed response files in the galactic coordinate have no time axis, and all of the events are assumed to be projected into a single galactic CDS.\n",
"In the future, we plan to introduce more flexible binning."
]
},
{
Expand Down Expand Up @@ -373,10 +380,10 @@
"# 3. Prepare a 'fake' coordsys conversion matrix\n",
"\n",
"The coordsys conversion matrix was initially introduced to convert a model map in the Galactic coordinates into the local coordinates.\n",
"In the case on this notebook, the CDS is in the Galactic coordinates, thus ideally we do not have to convert the coordinates of the model map.\n",
"However, as for now, the code for the image deconvolution was mainly developed for the CDS in the local coordinates, so here a ‘fake’ coordinate conversion matrix will be created. \n",
"Then, the same code can be used for both methods. \n",
"We will consider to remove this fake coordinate conversion matrix in the future."
"In the case of this notebook, the CDS is in the Galactic coordinates; thus, ideally, we do not have to convert the coordinates of the model map.\n",
"However, as for now, the code for the image deconvolution was mainly developed for the CDS in the local coordinates and requires the conversion matrix, so here, we generate a ‘fake’ coordinate conversion matrix, which is an unit matrix. \n",
"Then, the same code can be applied for both methods. \n",
"We will consider removing this fake coordinate conversion matrix in the future."
]
},
{
Expand Down Expand Up @@ -453,7 +460,9 @@
{
"cell_type": "markdown",
"id": "6e88ca7f",
"metadata": {},
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"## Brief overview of the image deconvolution\n",
"\n",
Expand All @@ -463,44 +472,50 @@
"\\log L = \\sum_i X_i \\log \\epsilon_i - \\sum_i \\epsilon_i\n",
"$$\n",
"\n",
"$X_i$: detected counts ( $i$ : index of the Compton Data Space)\n",
"$X_i$: detected counts at $i$-th bin ( $i$ : index of the Compton Data Space)\n",
"\n",
"$\\epsilon_i = \\sum_j R_{ij} \\lambda_j + b_i$ : expected counts ( $j$ : index of the model space)\n",
"\n",
"$\\lambda_j$ : the model map\n",
"$\\lambda_j$ : the model map (basically gamma-ray flux at $j$-th pixel)\n",
"\n",
"$b_i$ : the background at the index $i$\n",
"$b_i$ : the background at $i$-th bin\n",
"\n",
"$R_{ij}$ : the response matrix\n",
"\n",
"Since we have to optimize the flux in each pixel, and the number of parameters are large, we adopt an iterative approach to find a solution of the above equation. The simplest one is ML-EM (maximum likelihood expectation maximazation) algorithm.\n",
"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.\n",
"\n",
"$$\n",
"\\lambda_{j}^{k+1} = \\lambda_{j}^{k} + \\delta \\lambda_{j}^{k} \\\\\n",
"\\lambda_{j}^{k+1} = \\lambda_{j}^{k} + \\delta \\lambda_{j}^{k}\n",
"$$\n",
"$$\n",
"\\delta \\lambda_{j}^{k} = \\frac{\\lambda_{j}^{k}}{\\sum_{i} R_{ij}} \\sum_{i} \\left(\\frac{ X_{i} }{\\epsilon_{i}} - 1 \\right) R_{ij} \n",
"$$\n",
"\n",
"We refer to $\\delta \\lambda_{j}^{k}$ as the delta map.\n",
"\n",
"As for now, the two improved algorithms are implemented.\n",
"As for now, the two improved algorithms are implemented in COSIpy.\n",
"\n",
"- Accelerated ML-EM algorithm (Knoedlseder+99)\n",
"\n",
"$$\n",
"\\lambda_{j}^{k+1} = \\lambda_{j}^{k} + \\alpha^{k} \\delta \\lambda_{j}^{k} \\\\\n",
"\\lambda_{j}^{k+1} = \\lambda_{j}^{k} + \\alpha^{k} \\delta \\lambda_{j}^{k}\n",
"$$\n",
"$$\n",
"\\alpha^{k} < \\mathrm{max}(- \\lambda_{j}^{k} / \\delta \\lambda_{j}^{k})\n",
"$$\n",
"\n",
"Paractically, in order not to accelerate the algorithm excessively, we set the maximu value of $\\alpha$. Thus, $\\alpha$ can be determined as:\n",
"Practically, in order not to accelerate the algorithm excessively, we set the maximum value of $\\alpha$ ($\\alpha_{\\mathrm{max}}$). Then, $\\alpha$ is calculated as:\n",
"\n",
"$$\n",
"\\alpha^{k} = \\mathrm{min}(\\mathrm{max}(- \\lambda_{j}^{k} / \\delta \\lambda_{j}^{k}), \\alpha_{max})\n",
"\\alpha^{k} = \\mathrm{min}(\\mathrm{max}(- \\lambda_{j}^{k} / \\delta \\lambda_{j}^{k}), \\alpha_{\\mathrm{max}})\n",
"$$\n",
"\n",
"- Noise damping using gaussian smoothing (Knoedlseder+05, Siegert+20)\n",
"\n",
"$$\n",
"\\lambda_{j}^{k+1} = \\lambda_{j}^{k} + \\alpha^{k} \\left[ w_j \\delta \\lambda_{j}^{k} \\right]_{\\mathrm{gauss}} \\\\\n",
"\\lambda_{j}^{k+1} = \\lambda_{j}^{k} + \\alpha^{k} \\left[ w_j \\delta \\lambda_{j}^{k} \\right]_{\\mathrm{gauss}}\n",
"$$\n",
"$$\n",
"w_j = \\left(\\sum_{i} R_{ij}\\right)^\\beta\n",
"$$\n",
"\n",
Expand All @@ -527,7 +542,7 @@
"dataloader.event_dense = event\n",
"dataloader.bkg_dense = bkg\n",
"\n",
"# the loaded response matrix should be assigned to both full_detector_response/image_response_dense.\n",
"# the loaded response matrix should be assigned to both full_detector_response/image_response_dense in the Galactic CDS method.\n",
"dataloader.full_detector_response = image_response\n",
"dataloader.image_response_dense = image_response \n",
"\n",
Expand Down Expand Up @@ -588,15 +603,15 @@
"id": "241505ad",
"metadata": {},
"source": [
"(the following function is tentetive. It will be removed in the future.)"
"(In the future, we plan to remove the method \"_modify_axes.\")"
]
},
{
"cell_type": "markdown",
"id": "5bc6a570",
"metadata": {},
"source": [
"Calculate an auxiliary matrix for the deconvolution (mandatory)"
"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."
]
},
{
Expand Down Expand Up @@ -628,7 +643,7 @@
"source": [
"## 4-3. Initialize the instance of the image deconvolution class\n",
"\n",
"First we prepare an instance of ImageDeconvolution class, and then, resister the dataset, parameters for the deconvolution. After that, you can start the calculation."
"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."
]
},
{
Expand Down Expand Up @@ -681,42 +696,42 @@
"source": [
"### Initialize image_deconvolution\n",
"\n",
"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\n",
"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.\n",
"\n",
"I describe the parameters in the parameter file.\n",
"I describe parameters in the parameter file.\n",
"\n",
"#### The description of parameters: model_property\n",
"#### model_property\n",
"\n",
"| Name | Unit | Description | Notes |\n",
"| :---: | :---: | :---: | :---: |\n",
"| coordinate | str | the coordinate system of the model map | As for now, it must be 'galactic' |\n",
"| nside | int | NSIDE of the model map | As for now, it must be the same as that of the response matrix |\n",
"| nside | int | NSIDE of the model map | it must be the same as NSIDE of 'lb' axis of the coordinate conversion matrix|\n",
"| scheme | str | SCHEME of the model map | As for now, it must be 'ring' |\n",
"| 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 |\n",
"\n",
"#### The description of parameters: model_initialization\n",
"#### model_initialization\n",
"\n",
"| Name | Unit | Description | Notes |\n",
"| :---: | :---: | :---: | :---: |\n",
"| algorithm | str | the method name to initialize the model map | As for now, only 'flat' can be used |\n",
"| 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 number of energy_edges - 1\" |\n",
"| 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 |\n",
"\n",
"#### The description of parameters: deconvolution\n",
"#### deconvolution\n",
"\n",
"| Name | Unit | Description | Notes |\n",
"| :---: | :---: | :---: | :---: |\n",
"|algorithm | str | the name of the image deconvolution algorithm| As for now, only 'RL_memsave' is supported |\n",
"|algorithm | str | the name of the image deconvolution algorithm| As for now, only 'RL' is supported |\n",
"|||||\n",
"|parameter_RL_memsave:iteration | int | The maximum number of the iteration | |\n",
"|parameter_RL_memsave:acceleration | bool | whether the accelerated ML-EM algorithm (Knoedlseder+99) is used | |\n",
"|parameter_RL_memsave:alpha_max | float | the maximum value for the acceleration | |\n",
"|parameter_RL_memsave:save_results_each_iteration | bool | whether a updated model map, detal map, likelihood etc. are save at the end of each iteration | |\n",
"|parameter_RL_memsave:response_weighting | bool | whether a factor $w_j = (\\sum_{i} R_{ij})^{\\beta}$ for weighting the delta image is introduced (see Knoedlseder+05, Siegert+20) | |\n",
"|parameter_RL_memsave:response_weighting_index | float | $\\beta$ in the above equation | |\n",
"|parameter_RL_memsave:smoothing | bool | whether a Gaussian filter is used (see Knoedlseder+05, Siegert+20) | |\n",
"|parameter_RL_memsave:smoothing_FWHM | float, degree | the FWHM of the Gaussian in the filter | |\n",
"|parameter_RL_memsave:background_normalization_fitting | bool | whether the background normalization is optimized at each iteration | As for now, the same single background normalization factor is used in all of the time bins |\n",
"|parameter_RL_memsave:background_normalization_range | list of float | the range of the normalization factor | should be positive |"
"|parameter_RL:iteration | int | The maximum number of the iteration | |\n",
"|parameter_RL:acceleration | bool | whether the accelerated ML-EM algorithm (Knoedlseder+99) is used | |\n",
"|parameter_RL:alpha_max | float | the maximum value for the acceleration parameter | |\n",
"|parameter_RL:save_results_each_iteration | bool | whether an updated model map, detal map, likelihood etc. are saved at the end of each iteration | |\n",
"|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) | |\n",
"|parameter_RL:response_weighting_index | float | $\\beta$ in the above equation | |\n",
"|parameter_RL:smoothing | bool | whether a Gaussian filter is used (see Knoedlseder+05, Siegert+20) | |\n",
"|parameter_RL:smoothing_FWHM | float, degree | the FWHM of the Gaussian in the filter | |\n",
"|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 |\n",
"|parameter_RL:background_normalization_range | list of float | the range of the normalization factor | should be positive |"
]
},
{
Expand Down Expand Up @@ -2075,7 +2090,7 @@
"metadata": {},
"source": [
"# 5. Analyze the results\n",
"Below examples to see/analyze the results are shown."
"Examples to see/analyze the results are shown below."
]
},
{
Expand Down Expand Up @@ -2240,7 +2255,7 @@
"id": "4b8cdf58",
"metadata": {},
"source": [
"Plotting the reconstructed images in all of the energy bands at the 20th iteration"
"Plotting the reconstructed images in all of the energy bands at the 50th iteration"
]
},
{
Expand Down
Loading

0 comments on commit ffc19ef

Please sign in to comment.