From ffc19ef6dfc4ce3c9459c37ff8173b8d8b2a3013 Mon Sep 17 00:00:00 2001 From: Hiroki Yoneda Date: Tue, 30 Jan 2024 17:00:43 +0100 Subject: [PATCH] minor changes on the example notebooks of the image deconvolution --- ...1keV-DC2-Galactic-ImageDeconvolution.ipynb | 103 ++++++++++------- .../511keV-DC2-ScAtt-DataReduction.ipynb | 62 +++++++--- .../511keV-DC2-ScAtt-ImageDeconvolution.ipynb | 80 +++++++------ .../511keV-DC2-ScAtt-Upsampling.ipynb | 33 +++--- .../511keV/ScAttBinning/mk_ccm_upsampling.py | 2 +- ...Crab-DC2-Galactic-ImageDeconvolution.ipynb | 101 ++++++++++------- .../Crab-DC2-ScAtt-DataReduction.ipynb | 52 ++++++--- .../Crab-DC2-ScAtt-ImageDeconvolution.ipynb | 82 +++++++------- .../Crab-DC2-ScAtt-Upsampling.ipynb | 69 +++-------- .../GRB/miniDC2-GRB-ImageDeconvolution.ipynb | 107 +++++++++--------- 10 files changed, 376 insertions(+), 315 deletions(-) diff --git a/docs/tutorials/image_deconvolution/511keV/GalacticCDS/511keV-DC2-Galactic-ImageDeconvolution.ipynb b/docs/tutorials/image_deconvolution/511keV/GalacticCDS/511keV-DC2-Galactic-ImageDeconvolution.ipynb index 288236ac..fda3fd46 100644 --- a/docs/tutorials/image_deconvolution/511keV/GalacticCDS/511keV-DC2-Galactic-ImageDeconvolution.ipynb +++ b/docs/tutorials/image_deconvolution/511keV/GalacticCDS/511keV-DC2-Galactic-ImageDeconvolution.ipynb @@ -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." ] }, { @@ -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", @@ -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" ] }, { @@ -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." ] }, { @@ -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." ] }, { @@ -453,7 +460,9 @@ { "cell_type": "markdown", "id": "6e88ca7f", - "metadata": {}, + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, "source": [ "## Brief overview of the image deconvolution\n", "\n", @@ -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", @@ -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", @@ -588,7 +603,7 @@ "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.\")" ] }, { @@ -596,7 +611,7 @@ "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." ] }, { @@ -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." ] }, { @@ -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 |" ] }, { @@ -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." ] }, { @@ -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" ] }, { diff --git a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-DataReduction.ipynb b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-DataReduction.ipynb index 72f8e519..77d42138 100644 --- a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-DataReduction.ipynb +++ b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-DataReduction.ipynb @@ -7,11 +7,30 @@ "source": [ "# DC2 Image Analysis, 511 keV, Data Reduction\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 how to produce the binned datasets with the spacecraft attitude (scatt) binning method for DC2.\n", - "Using the 511keV thin disk 3month simulation data created for DC2, an example of the image analysis will be presented.\n", - "After running through this notebook, you can go to the next notebook, namely 511keV-DC2-ScAtt-ImageDeconvolution.ipynb.\n", + "Using the 511keV thin disk 3-month simulation data created for DC2, an example of the image analysis will be presented.\n", + "After running through this notebook, you can go to the next notebook, 511keV-DC2-ScAtt-ImageDeconvolution.ipynb.\n", + "\n", + "### Notes on the coordinate system of Compton data space in the image deconvolution ###\n", + "\n", + "We have two options on the coordinate system to describe the Compton scattering direction ($\\chi\\psi$) with, namely the Galactic coordinate or the detector coordinate.\n", + "\n", + "Using the Galactic coordinate is intuitive, and the spectral fitting adopts this coordinate. Thus, we suppose that Galactic coordinate should be adopted also for image deconvolution eventually. However, in this case, we need to convert the detector response into the Galactic coordinate for each pixel in the sky because the response matrix is described in the detector coordinate. As for now, it takes a long time to compute it. Thus, the pre-computed converted response are provided in DC2 for several main sources (511 keV, Al-26, Ti-44, continuum). The pre-computed responses assume that we analyze 3-month data without extracting some time intervals, and the pixel resolution of the model map is already fixed in them. While there is less flexibility in binning/modeling, it is relatively fast to perform the image deconvolution in DC2 since the most computationally heavy part, the coordinate conversion of the response, can be skipped.\n", + "\n", + "Using the detector coordinates for Compton data space may not be so intuitive. However, the advantage is that we do not have to convert the response matrix. Instead, we will convert the model map into the detector coordinate. Because the model map generally has a much smaller data size than the response, we can compute this coordinate conversion quickly. \n", + "\n", + "The disadvantage of this method is that we need more bins due to continuous pointing changes of the COSI satellite. Since COSI is an all-sky monitoring satellite with ∼90-minute orbits, it changes its pointing by ∼4 degrees every minute. Thus, in this case, we need to divide the data into several bins so that astronomical sources can be considered at rest in the detector coordinate for each bin within the COSI's angular resolution. The straightforward way could be to divide the data every $\\sim$15 seconds, considering that the COSI's angular resolution is an order of degrees. However, we need $5\\times10^5$ time bins for 3-month observations, which makes the event histogram very huge. To avoid this issue, the spacecraft attitude (scatt) binning method is introduced. Instead of binning data over time, we first analyze the satellite attitude and find the time intervals when the satellite has almost the same attitude within the angular resolution. Then, we assign the events in such intervals into the same CDS. In the DC2 simulation, the orbit inclination is assumed to be 0 degrees. In this case, the number of the scatt bins becomes 100-1000, which makes the computation more executable. With this method, at least in DC2, we can perform the image deconvolution using the original response matrix and have flexibilities to change binning/modeling, e.g., the pixel resolution can be changed in a relatively easy way.\n", + "\n", + "While both methods have pros and cons, our baseline is to eventually use the Galactic coordinate. But we still need to carefully investigate how they will be scaled with longer exposure, finer pixel resolution, etc. Thus, we provide the notebooks of both methods for the image deconvolution in DC2.\n", + "\n", + "For the Crab image analysis, the following notebooks are based on the scatt binning method\n", + "- ScAttBinning/511keV-DC2-ScAtt-DataReduction.ipynb\n", + "- ScAttBinning/511keV-DC2-ScAtt-ImageDeconvolution.ipynb\n", + "- ScAttBinning/511keV-DC2-ScAtt-Upsampling.ipynb\n", + "\n", + "GalacticCDS/511keV-DC2-Galactic-ImageDeconvolution.ipynb uses the galactic coordinate.\n", "\n", "If you want to know about the other analysis, e.g., the spectral analysis, you can see the notebooks in docs/tutorials/spectral_fits." ] @@ -385,13 +404,16 @@ "# 0. Prepare the data\n", "Before running the cells, please download the files needed for this notebook. You can get them from wasabi. \n", "\n", - "Actually, the data reduction is not optimized and takes hours depending on your environments. So I skip this process.\n", + "Basically, the data reduction from raw tra files may take hours depending on your environments. So we can skip this process.\n", "Please download the following data files and then run the following cells.\n", "\n", "From wasabi\n", "- cosi-pipeline-public/COSI-SMEX/DC2/Responses/SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5\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", "- cosi-pipeline-public/COSI-SMEX/DC2/Data/Orientation/20280301_3_month.ori\n", "\n", "From docs/tutorials/image_deconvolution/511keV/ScAttBinning\n", @@ -403,8 +425,7 @@ "id": "8462d0dc", "metadata": {}, "source": [ - "You can download the data and detector response from wasabi. You can skip this cell if you already have downloaded the files.\n", - "Note that the response is not public yet (2023-11-21)." + "You can download the data and detector response from wasabi. You can skip this cell if you already have downloaded the files." ] }, { @@ -1012,7 +1033,7 @@ "id": "e9306cf5", "metadata": {}, "source": [ - "SpacecraftAttitudeExposureTable can produce SpacecraftAttitudeMap that has an exposure time in each Z- and X-poiting pixels." + "SpacecraftAttitudeExposureTable can produce SpacecraftAttitudeMap that has an exposure time in each Z- and X-pointing pixels." ] }, { @@ -1117,16 +1138,15 @@ "source": [ "# 2. Calculate the coordinate conversion matrix\n", "\n", - "\n", "CoordsysConversionMatrix.spacecraft_attitude_binning_ccm can produce the coordinate conversion matrix for the spacecraft attitude binning.\n", "\n", - "In this calculation, the dwell time map is calculated for each model pixel and each scatt_binning_index.\n", + "In this calculation, we calculate the exposure time map in the detector coordinate for each model pixel and each scatt_binning_index. We refer to it as the dwell time map.\n", "\n", - "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table), and then the dwell time map is calculated once for ach model pixel and each scatt_binning_index.\n", + "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table) for each scatt_binning_index, and then the dwell time map is calculated assuming the averaged pointings for each model pixel and each scatt_binning_index.\n", "\n", - "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing in the exposure table, and then the calculated dwell time maps are summed up. \n", + "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing (basically every 1 second), and then the calculated dwell time maps are summed up for each model pixel and each scatt_binning_index. \n", "\n", - "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate but it takes a long time to calculate it." + "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate, but it takes a very long time to calculate it." ] }, { @@ -1209,7 +1229,7 @@ "source": [ "# 3. produce the binned data\n", "\n", - "Using the exposure table, we can produce the binned data." + "Using the exposure table, we can produce the binned data. Note that here we generate the binned histogram manually. We consider implementing this functionality in the DataIO class in the future." ] }, { @@ -1372,6 +1392,22 @@ "binned_event.write(\"511keV_scatt_binning_DC2_event.hdf5\")\n", "binned_bkg.write(\"511keV_scatt_binning_DC2_bkg.hdf5\")" ] + }, + { + "cell_type": "markdown", + "id": "aa3e61c2-976a-4b08-befb-b607ed510442", + "metadata": {}, + "source": [ + "**You can move on the next notebook (511keV-DC2-ScAtt-ImageDeconvolution.ipynb).**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00723a19-9304-4625-94e1-998247098e0c", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-ImageDeconvolution.ipynb b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-ImageDeconvolution.ipynb index adca7508..2d25a249 100644 --- a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-ImageDeconvolution.ipynb +++ b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-ImageDeconvolution.ipynb @@ -9,10 +9,10 @@ "source": [ "# DC2 Image Analysis, 511 keV, Image Deconvolution\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 spacecraft attitude (scatt) binning method for DC2.\n", - "Using the 511 keV thin disk 3month simulation data created for DC2, an example of the image analysis will be presented.\n", + "Using the 511 keV thin disk 3-month simulation data created for DC2, an example of the image analysis will be presented.\n", "If you have not run through 511keV-DC2-ScAtt-DataReduction.ipynb, please see it first." ] }, @@ -577,44 +577,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", @@ -696,7 +702,7 @@ "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.\")" ] }, { @@ -706,9 +712,9 @@ "source": [ "## 4-2. Load the response file\n", "\n", - "The response file will be loaded on the CPU memory. It requires a few GB. In the actuall analysis, the response will be much larger, ~TB. \n", + "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. \n", "\n", - "So loading it on the memory might be unrealistic. The optimized (lazy) loading would be a next work." + "So loading it on the memory might be unrealistic in the future. The optimized (lazy) loading would be a next work." ] }, { @@ -737,7 +743,7 @@ "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." ] }, { @@ -765,7 +771,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." ] }, { @@ -818,42 +824,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 |" ] }, { @@ -1760,7 +1766,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." ] }, { @@ -1852,7 +1858,7 @@ "source": [ "## Background normalization\n", "\n", - "Plotting the background nomalization factor vs the number of iterations. If the backgroud model is accurate and the image is reconstructed perfectly, this factor should be close to 1." + "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." ] }, { diff --git a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-Upsampling.ipynb b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-Upsampling.ipynb index 2ea9cafe..577ddc28 100644 --- a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-Upsampling.ipynb +++ b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/511keV-DC2-ScAtt-Upsampling.ipynb @@ -7,30 +7,31 @@ "source": [ "# DC2 Image Analysis, 511keV, Upsampling\n", "\n", - "updated on 2024-01-23 (the commit fe724f0f8d317ae347596708e6c50db6899c7e89)\n", + "updated on 2024-01-30 (the commit 26cfdeacb25335bd511a91c4f8a29bdeb36408f2)\n", "\n", - "This notebook explains image reconstruction when the pixel resolution of the model map is set higher than the response matrix.\n", + "This notebook explains image reconstruction using the pixel resolution of the model map finer than that of the response matrix.\n", "\n", "Note that this notebook is advanced. It is assumed that you have already performed the two notebooks (511keV-DC2-ScAtt-DataReduction.ipynb, 511keV-DC2-ScAtt-ImageDeconvolution.ipynb).\n", "\n", "## Point\n", "\n", - "In the current implementation, the pixel size of the model map can be differnt from the that of the response. The model pixel size is used in the following instances:\n", + "In the current implementation, the pixel size of the model map can be differnt from that of the response matrix. The model pixel size is used in the following instances:\n", "\n", "- coordsys_conv_matrix\n", "- image_deconvolution\n", "\n", "Thus, make sure that NSIDE in these instances must be the same. In this notebook, I present the case with NSIDE = 32 in the model map.\n", "\n", - "If the NSIDE for the exposure table is not changed, you do not have to create the binned data.\n", + "When we convert the model map in the galactic coordinate to the detector coordinate, the pixel size will be downscaled so as the converted model map has the same pixel resolution matching the detector response.\n", + "Thus, using finer resolution in the model space does not improve the angular resolution in principle, while the reconstructed image will be smoother.\n", "\n", - "There are three different NSIDE in priciple:\n", + "There are three different NSIDE defined in the analysis:\n", "\n", "- NSIDE for the pixel resolution of the model (coordsys_conv_matrix, image_deconvolution)\n", "- NSIDE for the pixel resolution of the response/data/background CDS (full_detector_response, inputs_511keV_DC2.yaml)\n", "- NSIDE for the pixel resolution of the spacecraftattitude binning (exposure_table)\n", "\n", - "Usually, these three values can be set equal." + "Normally, these three values are set equal, but in principle they can be different." ] }, { @@ -95,14 +96,10 @@ }, { "cell_type": "markdown", - "id": "2a7ca026", + "id": "7d93d4e9-d70f-41b5-93b6-fa8c556403c8", "metadata": {}, "source": [ "# 0. Prepare the data\n", - "Before running the cells, please download the files needed for this notebook. You can get them from wasabi. \n", - "\n", - "Actually, the data reduction is not optimized and takes hours depending on your environments. So I skip this process.\n", - "Please download the following data files and then run the following cells.\n", "\n", "From wasabi\n", "- cosi-pipeline-public/COSI-SMEX/DC2/Responses/SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5\n", @@ -651,13 +648,13 @@ "\n", "CoordsysConversionMatrix.spacecraft_attitude_binning_ccm can produce the coordinate conversion matrix for the spacecraft attitude binning.\n", "\n", - "In this calculation, the dwell time map is calculated for each model pixel and each scatt_binning_index.\n", + "In this calculation, we calculate the exposure time map in the detector coordinate for each model pixel and each scatt_binning_index. We refer to it as the dwell time map.\n", "\n", - "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table), and then the dwell time map is calculated once for ach model pixel and each scatt_binning_index.\n", + "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table) for each scatt_binning_index, and then the dwell time map is calculated assuming the averaged pointings for each model pixel and each scatt_binning_index.\n", "\n", - "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing in the exposure table, and then the calculated dwell time maps are summed up. \n", + "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing (basically every 1 second), and then the calculated dwell time maps are summed up for each model pixel and each scatt_binning_index. \n", "\n", - "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate but it takes a long time to calculate it.\n", + "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate, but it takes a very long time to calculate it.\n", "\n", "**With NSIDE = 32, this process may take a few hours. I also prepared a python script to create the conversion matrix. If the notebook does not work, you can use mk_ccm_upsampling.py**" ] @@ -1025,7 +1022,7 @@ "id": "89aeb1ad", "metadata": {}, "source": [ - "**Do not forget to make sure that NSIDE of the model map is modified**" + "**Do not forget to make sure that NSIDE of the model map is modified to 32**" ] }, { @@ -1651,7 +1648,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." ] }, { @@ -1761,7 +1758,7 @@ "source": [ "## Background normalization\n", "\n", - "Plotting the background nomalization factor vs the number of iterations. If the backgroud model is accurate and the image is reconstructed perfectly, this factor should be close to 1." + "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." ] }, { diff --git a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/mk_ccm_upsampling.py b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/mk_ccm_upsampling.py index 2d3c362a..be3320f5 100644 --- a/docs/tutorials/image_deconvolution/511keV/ScAttBinning/mk_ccm_upsampling.py +++ b/docs/tutorials/image_deconvolution/511keV/ScAttBinning/mk_ccm_upsampling.py @@ -1,7 +1,7 @@ from cosipy.response import FullDetectorResponse from cosipy.image_deconvolution import SpacecraftAttitudeExposureTable, CoordsysConversionMatrix -full_detector_response_filename = "/Users/yoneda/Work/Exp/COSI/cosipy-2/data_challenge/DC2/prework/data/Responses/SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5" +full_detector_response_filename = "/Users/yoneda/Work/Exp/COSI/cosipy-2/data_challenge/DC2/prework/data/Responses/SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5" # Please replace this with your file path full_detector_response = FullDetectorResponse.open(full_detector_response_filename) exposure_table = SpacecraftAttitudeExposureTable.from_fits("exposure_table_nside32.fits") diff --git a/docs/tutorials/image_deconvolution/Crab/GalacticCDS/Crab-DC2-Galactic-ImageDeconvolution.ipynb b/docs/tutorials/image_deconvolution/Crab/GalacticCDS/Crab-DC2-Galactic-ImageDeconvolution.ipynb index 519ca8a5..64f0dde4 100644 --- a/docs/tutorials/image_deconvolution/Crab/GalacticCDS/Crab-DC2-Galactic-ImageDeconvolution.ipynb +++ b/docs/tutorials/image_deconvolution/Crab/GalacticCDS/Crab-DC2-Galactic-ImageDeconvolution.ipynb @@ -9,10 +9,12 @@ "source": [ "# DC2 Image Analysis, Crab, 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 Crab 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 the Crab 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 Crab-DC2-ScAtt-DataReduction.ipynb." ] }, { @@ -67,9 +69,12 @@ "\n", "From wasabi\n", "- cosi-pipeline-public/COSI-SMEX/DC2/Responses/PointSourceReponse/psr_gal_continuum_DC2.h5.zip (please unzip it)\n", - " - a pre-computed continuum response file rotated into the Galactic coordinate system\n", + " - a pre-computed continuum response file converted into the Galactic coordinate system\n", "- cosi-pipeline-public/COSI-SMEX/DC2/Data/Sources/crab_3months_unbinned_data.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/Crab/GalacticCDS\n", "- inputs_Crab_DC2.yaml\n", @@ -204,7 +209,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" ] }, { @@ -264,7 +269,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." ] }, { @@ -374,10 +381,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." ] }, { @@ -464,44 +471,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", @@ -528,7 +541,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", @@ -589,7 +602,7 @@ "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.\")" ] }, { @@ -597,7 +610,7 @@ "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." ] }, { @@ -629,7 +642,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." ] }, { @@ -682,42 +695,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 |" ] }, { @@ -2112,7 +2125,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." ] }, { @@ -2277,7 +2290,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" ] }, { @@ -2470,7 +2483,7 @@ "source": [ "## Spectrum\n", "\n", - "Plotting the gamma-ray spectrum at 20th interation. The photon flux at each energy band shown here is calculated as the accumulation of the flux values in all pixel at each energy band." + "Plotting the gamma-ray spectrum at the 50th iteration. The photon flux at each energy band shown here is calculated as the accumulation of the flux values in all pixels at each energy band." ] }, { diff --git a/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-DataReduction.ipynb b/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-DataReduction.ipynb index 1365dde0..1e88ea8c 100644 --- a/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-DataReduction.ipynb +++ b/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-DataReduction.ipynb @@ -7,11 +7,28 @@ "source": [ "# DC2 Image Analysis, Crab, Data Reduction\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 how to produce the binned datasets with the spacecraft attitude (scatt) binning method for DC2.\n", - "Using the Crab 3month simulation data created for DC2, an example of the image analysis will be presented.\n", - "After running through this notebook, you can go to the next notebook, namely Crab-DC2-ScAtt-ImageDeconvolution.ipynb.\n", + "This notebook focuses on how to produce the binned datasets with the spacecraft attitude (scatt) binning method for DC2. An example of the image analysis will be presented using the Crab 3-month simulation data created for DC2. After running through this notebook, you can go to the next notebook, Crab-DC2-ScAtt-ImageDeconvolution.ipynb.\n", + "\n", + "### Notes on the coordinate system of Compton data space in the image deconvolution ###\n", + "\n", + "We have two options on the coordinate system to describe the Compton scattering direction ($\\chi\\psi$) with, namely the Galactic coordinate or the detector coordinate.\n", + "\n", + "Using the Galactic coordinate is intuitive, and the spectral fitting adopts this coordinate. Thus, we suppose that Galactic coordinate should be adopted also for image deconvolution eventually. However, in this case, we need to convert the detector response into the Galactic coordinate for each pixel in the sky because the response matrix is described in the detector coordinate. As for now, it takes a long time to compute it. Thus, the pre-computed converted response are provided in DC2 for several main sources (511 keV, Al-26, Ti-44, continuum). The pre-computed responses assume that we analyze 3-month data without extracting some time intervals, and the pixel resolution of the model map is already fixed in them. While there is less flexibility in binning/modeling, it is relatively fast to perform the image deconvolution in DC2 since the most computationally heavy part, the coordinate conversion of the response, can be skipped.\n", + "\n", + "Using the detector coordinates for Compton data space may not be so intuitive. However, the advantage is that we do not have to convert the response matrix. Instead, we will convert the model map into the detector coordinate. Because the model map generally has a much smaller data size than the response, we can compute this coordinate conversion quickly. \n", + "\n", + "The disadvantage of this method is that we need more bins due to continuous pointing changes of the COSI satellite. Since COSI is an all-sky monitoring satellite with ∼90-minute orbits, it changes its pointing by ∼4 degrees every minute. Thus, in this case, we need to divide the data into several bins so that astronomical sources can be considered at rest in the detector coordinate for each bin within the COSI's angular resolution. The straightforward way could be to divide the data every $\\sim$15 seconds, considering that the COSI's angular resolution is an order of degrees. However, we need $5\\times10^5$ time bins for 3-month observations, which makes the event histogram very huge. To avoid this issue, the spacecraft attitude (scatt) binning method is introduced. Instead of binning data over time, we first analyze the satellite attitude and find the time intervals when the satellite has almost the same attitude within the angular resolution. Then, we assign the events in such intervals into the same CDS. In the DC2 simulation, the orbit inclination is assumed to be 0 degrees. In this case, the number of the scatt bins becomes 100-1000, which makes the computation more executable. With this method, at least in DC2, we can perform the image deconvolution using the original response matrix and have flexibilities to change binning/modeling, e.g., the pixel resolution can be changed in a relatively easy way.\n", + "\n", + "While both methods have pros and cons, our baseline is to eventually use the Galactic coordinate. But we still need to carefully investigate how they will be scaled with longer exposure, finer pixel resolution, etc. Thus, we provide the notebooks of both methods for the image deconvolution in DC2.\n", + "\n", + "For the Crab image analysis, the following notebooks are based on the scatt binning method\n", + "- ScAttBinning/Crab-DC2-ScAtt-DataReduction.ipynb\n", + "- ScAttBinning/Crab-DC2-ScAtt-ImageDeconvolution.ipynb\n", + "- ScAttBinning/Crab-DC2-ScAtt-Upsampling.ipynb\n", + "\n", + "GalacticCDS/Crab-DC2-Galactic-ImageDeconvolution.ipynb uses the galactic coordinate.\n", "\n", "If you want to know about the other analysis, e.g., the spectral analysis, you can see the notebooks in docs/tutorials/spectral_fits." ] @@ -385,13 +402,16 @@ "# 0. Prepare the data\n", "Before running the cells, please download the files needed for this notebook. You can get them from wasabi. \n", "\n", - "Actually, the data reduction is not optimized and takes hours depending on your environments. So I skip this process.\n", + "Basically, the data reduction from raw tra files may take hours depending on your environments. So we can skip this process.\n", "Please download the following data files and then run the following cells.\n", "\n", "From wasabi\n", "- cosi-pipeline-public/COSI-SMEX/DC2/Responses/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5 (please unzip it)\n", "- cosi-pipeline-public/COSI-SMEX/DC2/Data/Sources/crab_3months_unbinned_data.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", "- cosi-pipeline-public/COSI-SMEX/DC2/Data/Orientation/20280301_3_month.ori\n", "\n", "From docs/tutorials/image_deconvolution/Crab/ScAttBinning\n", @@ -403,8 +423,7 @@ "id": "8462d0dc", "metadata": {}, "source": [ - "You can download the data and detector response from wasabi. You can skip this cell if you already have downloaded the files.\n", - "Note that the response is not public yet (2023-11-21)." + "You can download the data and detector response from wasabi. You can skip this cell if you already have downloaded the files." ] }, { @@ -1097,16 +1116,15 @@ "source": [ "# 2. Calculate the coordinate conversion matrix\n", "\n", - "\n", "CoordsysConversionMatrix.spacecraft_attitude_binning_ccm can produce the coordinate conversion matrix for the spacecraft attitude binning.\n", "\n", - "In this calculation, the dwell time map is calculated for each model pixel and each scatt_binning_index.\n", + "In this calculation, we calculate the exposure time map in the detector coordinate for each model pixel and each scatt_binning_index. We refer to it as the dwell time map.\n", "\n", - "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table), and then the dwell time map is calculated once for ach model pixel and each scatt_binning_index.\n", + "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table) for each scatt_binning_index, and then the dwell time map is calculated assuming the averaged pointings for each model pixel and each scatt_binning_index.\n", "\n", - "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing in the exposure table, and then the calculated dwell time maps are summed up. \n", + "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing (basically every 1 second), and then the calculated dwell time maps are summed up for each model pixel and each scatt_binning_index. \n", "\n", - "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate but it takes a long time to calculate it." + "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate, but it takes a very long time to calculate it." ] }, { @@ -1351,10 +1369,18 @@ "binned_bkg.write(\"Crab_scatt_binning_DC2_bkg.hdf5\")" ] }, + { + "cell_type": "markdown", + "id": "d8290e73-3a82-4308-ab2f-96da00f67bf5", + "metadata": {}, + "source": [ + "**You can move on the next notebook (Crab-DC2-ScAtt-ImageDeconvolution.ipynb).**" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "12cc3b12", + "id": "7a6d5643-989b-4efc-b30a-1676cfdd02f7", "metadata": {}, "outputs": [], "source": [] diff --git a/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-ImageDeconvolution.ipynb b/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-ImageDeconvolution.ipynb index 310ff820..675c3d35 100644 --- a/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-ImageDeconvolution.ipynb +++ b/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-ImageDeconvolution.ipynb @@ -9,10 +9,10 @@ "source": [ "# DC2 Image Analysis, Crab, Image Deconvolution\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 spacecraft attitude (scatt) binning method for DC2.\n", - "Using the Crab 3month simulation data created for DC2, an example of the image analysis will be presented.\n", + "Using the Crab 3-month simulation data created for DC2, an example of the image analysis will be presented.\n", "If you have not run through Crab-DC2-ScAtt-DataReduction.ipynb, please see it first." ] }, @@ -590,44 +590,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", @@ -709,7 +715,7 @@ "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.\")" ] }, { @@ -719,9 +725,9 @@ "source": [ "## 4-2. Load the response file\n", "\n", - "The response file will be loaded on the CPU memory. It requires a few GB. In the actuall analysis, the response will be much larger, ~TB. \n", + "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. \n", "\n", - "So loading it on the memory might be unrealistic. The optimized (lazy) loading would be a next work." + "So loading it on the memory might be unrealistic in the future. The optimized (lazy) loading would be a next work." ] }, { @@ -750,7 +756,7 @@ "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." ] }, { @@ -782,7 +788,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." ] }, { @@ -835,42 +841,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 |" ] }, { @@ -1602,7 +1608,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." ] }, { @@ -1712,7 +1718,7 @@ "source": [ "## Background normalization\n", "\n", - "Plotting the background nomalization factor vs the number of iterations. If the backgroud model is accurate and the image is reconstructed perfectly, this factor should be close to 1." + "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." ] }, { @@ -2191,7 +2197,7 @@ "source": [ "## Spectrum\n", "\n", - "Plotting the gamma-ray spectrum at 20th interation. The photon flux at each energy band shown here is calculated as the accumulation of the flux values in all pixel at each energy band." + "Plotting the gamma-ray spectrum at 11th interation. The photon flux at each energy band shown here is calculated as the accumulation of the flux values in all of the pixels at each energy band." ] }, { diff --git a/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-Upsampling.ipynb b/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-Upsampling.ipynb index bbadd065..423b6aba 100644 --- a/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-Upsampling.ipynb +++ b/docs/tutorials/image_deconvolution/Crab/ScAttBinning/Crab-DC2-ScAtt-Upsampling.ipynb @@ -7,30 +7,31 @@ "source": [ "# DC2 Image Analysis, Crab, Upsampling\n", "\n", - "updated on 2024-01-23 (the commit fe724f0f8d317ae347596708e6c50db6899c7e89)\n", + "updated on 2024-01-30 (the commit 26cfdeacb25335bd511a91c4f8a29bdeb36408f2)\n", "\n", - "This notebook explains image reconstruction when the pixel resolution of the model map is set higher than the response matrix.\n", + "This notebook explains image reconstruction using the pixel resolution of the model map finer than that of the response matrix.\n", "\n", "Note that this notebook is advanced. It is assumed that you have already performed the two notebooks (Crab-DC2-ScAtt-DataReduction.ipynb, Crab-DC2-ScAtt-ImageDeconvolution.ipynb).\n", "\n", "## Point\n", "\n", - "In the current implementation, the pixel size of the model map can be differnt from the that of the response. The model pixel size is used in the following instances:\n", + "In the current implementation, the pixel size of the model map can be differnt from that of the response matrix. The model pixel size is used in the following instances:\n", "\n", "- coordsys_conv_matrix\n", "- image_deconvolution\n", "\n", "Thus, make sure that NSIDE in these instances must be the same. In this notebook, I present the case with NSIDE = 16 in the model map.\n", "\n", - "If the NSIDE for the exposure table is not changed, you do not have to create the binned data.\n", + "When we convert the model map in the galactic coordinate to the detector coordinate, the pixel size will be downscaled so as the converted model map has the same pixel resolution matching the detector response.\n", + "Thus, using finer resolution in the model space does not improve the angular resolution in principle, while the reconstructed image will be smoother.\n", "\n", - "There are three different NSIDE in priciple:\n", + "There are three different NSIDE defined in the analysis:\n", "\n", "- NSIDE for the pixel resolution of the model (coordsys_conv_matrix, image_deconvolution)\n", - "- NSIDE for the pixel resolution of the response (full_detector_response)\n", + "- NSIDE for the pixel resolution of the response/data/background CDS (full_detector_response, inputs_Crab_DC2.yaml)\n", "- NSIDE for the pixel resolution of the spacecraftattitude binning (exposure_table)\n", "\n", - "Usually, these three values can be set equal." + "Normally, these three values are set equal, but in principle they can be different." ] }, { @@ -99,10 +100,6 @@ "metadata": {}, "source": [ "# 0. Prepare the data\n", - "Before running the cells, please download the files needed for this notebook. You can get them from wasabi. \n", - "\n", - "Actually, the data reduction is not optimized and takes hours depending on your environments. So I skip this process.\n", - "Please download the following data files and then run the following cells.\n", "\n", "From wasabi\n", "- cosi-pipeline-public/COSI-SMEX/DC2/Responses/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5 (please unzip it)\n", @@ -117,35 +114,6 @@ "- Crab_scatt_binning_DC2_event.hdf5" ] }, - { - "cell_type": "markdown", - "id": "8462d0dc", - "metadata": {}, - "source": [ - "You can download the data and detector response from wasabi. You can skip this cell if you already have downloaded the files.\n", - "Note that the response is not public yet (2023-11-21)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5173596b", - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "\n", - "header = \"AWS_ACCESS_KEY_ID=GBAL6XATQZNRV3GFH9Y4 AWS_SECRET_ACCESS_KEY=GToOczY5hGX3sketNO2fUwiq4DJoewzIgvTCHoOv aws s3api get-object\"\n", - "\n", - "os.system(header + \" --bucket cosi-pipeline-public --key COSI-SMEX/DC2/Responses/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5.zip --endpoint-url=https://s3.us-west-1.wasabisys.com SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5.zip\")\n", - "\n", - "os.system(header + \" --bucket cosi-pipeline-public --key COSI-SMEX/DC2/Data/Sources/crab_3months_unbinned_data.fits.gz --endpoint-url=https://s3.us-west-1.wasabisys.com crab_3months_unbinned_data.fits.gz\")\n", - "\n", - "os.system(header + \" --bucket cosi-pipeline-public --key COSI-SMEX/DC2/Data/Backgrounds/albedo_photons_3months_unbinned_data.fits.gz --endpoint-url=https://s3.us-west-1.wasabisys.com albedo_photons_3months_unbinned_data.fits.gz\")\n", - "\n", - "os.system(header + \" --bucket cosi-pipeline-public --key COSI-SMEX/DC2/Data/Orientation/20280301_3_month.ori --endpoint-url=https://s3.us-west-1.wasabisys.com 20280301_3_month.ori\")" - ] - }, { "cell_type": "markdown", "id": "dc91fb24", @@ -681,16 +649,15 @@ "source": [ "# 2. Calculate the coordinate conversion matrix\n", "\n", - "\n", "CoordsysConversionMatrix.spacecraft_attitude_binning_ccm can produce the coordinate conversion matrix for the spacecraft attitude binning.\n", "\n", - "In this calculation, the dwell time map is calculated for each model pixel and each scatt_binning_index.\n", + "In this calculation, we calculate the exposure time map in the detector coordinate for each model pixel and each scatt_binning_index. We refer to it as the dwell time map.\n", "\n", - "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table), and then the dwell time map is calculated once for ach model pixel and each scatt_binning_index.\n", + "If use_averaged_pointing is True, first the averaged Z- and X-pointings are calculated (the average of zpointing or xpointing in the exposure table) for each scatt_binning_index, and then the dwell time map is calculated assuming the averaged pointings for each model pixel and each scatt_binning_index.\n", "\n", - "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing in the exposure table, and then the calculated dwell time maps are summed up. \n", + "If use_averaged_pointing is False, the dwell time map is calculated for each attitude in zpointing and xpointing (basically every 1 second), and then the calculated dwell time maps are summed up for each model pixel and each scatt_binning_index. \n", "\n", - "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate but it takes a long time to calculate it." + "In the former case, the computation is fast but may lose the angular resolution. In the latter case, the conversion matrix is more accurate, but it takes a very long time to calculate it." ] }, { @@ -801,9 +768,9 @@ "id": "4ae2fcdb", "metadata": {}, "source": [ - "# 3. produce the binned data\n", + "# 3. Load the binned data\n", "\n", - "skip this section mostly" + "Since NSIDE of exposure_table on this notebook is the same as that in Crab-DC2-ScAtt-DataReduction.ipynb, you can use the files generated before." ] }, { @@ -992,7 +959,7 @@ "source": [ "## 4-4. modify the parameters\n", "\n", - "**NSIDE for the model map must be the same as nside_model**" + "**Do not forget to make sure that NSIDE of the model map is modified to 16**" ] }, { @@ -1637,7 +1604,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." ] }, { @@ -1759,7 +1726,7 @@ "source": [ "## Background normalization\n", "\n", - "Plotting the background nomalization factor vs the number of iterations. If the backgroud model is accurate and the image is reconstructed perfectly, this factor should be close to 1." + "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." ] }, { @@ -1968,7 +1935,7 @@ "source": [ "## Spectrum\n", "\n", - "Plotting the gamma-ray spectrum at 20th interation. The photon flux at each energy band shown here is calculated as the accumulation of the flux values in all pixel at each energy band." + "Plotting the gamma-ray spectrum at 20th iteration. The photon flux at each energy band shown here is calculated as the accumulation of the flux values in all pixel at each energy band." ] }, { diff --git a/docs/tutorials/image_deconvolution/GRB/miniDC2-GRB-ImageDeconvolution.ipynb b/docs/tutorials/image_deconvolution/GRB/miniDC2-GRB-ImageDeconvolution.ipynb index cd67f3c6..6ca311e4 100644 --- a/docs/tutorials/image_deconvolution/GRB/miniDC2-GRB-ImageDeconvolution.ipynb +++ b/docs/tutorials/image_deconvolution/GRB/miniDC2-GRB-ImageDeconvolution.ipynb @@ -9,9 +9,9 @@ "source": [ "# GRB image analysis (miniDC2)\n", "\n", - "updated on 2024-01-23 (commit fe724f0f8d317ae347596708e6c50db6899c7e89)\n", + "updated on 2024-01-30 (commit 26cfdeacb25335bd511a91c4f8a29bdeb36408f2)\n", "\n", - "This notebook focuses on the image analysis with the image deconvolution code with the commit b7f97aafea23ef88466e62f45adfbdde223a4be9. Using the GRB simulation data created for miniDC2, an example of the image analysis will be presented.\n", + "Using the GRB simulation data created for miniDC2, an example of the image analysis will be presented.\n", "\n", "If you want to know about the other analysis, e.g., the spectral analysis, you can see the notebooks in docs/tutorials/spectral_fits." ] @@ -21,7 +21,7 @@ "id": "2bc243c8", "metadata": {}, "source": [ - "### Note that when the headline is inside parentheses it is not necessary to run the following cell, which is prepared for readers to understand the code clearly" + "**Note that it is not necessary to run the following cell when the headline is inside parentheses. These cells are prepared for readers to understand the code more clearly**" ] }, { @@ -391,9 +391,6 @@ "\n", "Before running the cells, please download the files needed for this notebook. You can get them from wasabi. \n", "\n", - "Actually, the data reduction is not optimized and takes hours depending on your environments. So I skip this process.\n", - "Please download the following data files and then run the following cells.\n", - "\n", "From wasabi\n", "- cosi-pipeline-public/ComptonSphere/mini-DC2/FlatContinuumIsotropic.LowRes.binnedimaging.imagingresponse.area.nside8.cosipy.h5.zip (please unzip it)\n", "- cosi-pipeline-public/ComptonSphere/mini-DC2/bkg_binned_data_full.hdf5\n", @@ -628,7 +625,7 @@ "id": "26f17d4c", "metadata": {}, "source": [ - "## Modify the axis (this is just for this case!)\n", + "## Modify the axis\n", "\n", "Here the time axis in the data and background files are modified as a single time bin. This is because the current code requires the same time intervals in both files.\n", "\n", @@ -636,7 +633,7 @@ "The background files is renormalized because the background is 2-hour data while the source data is 2-s duration.\n", "\n", "\n", - "Such file modification may be confusing, so it will be modified in the future." + "Such a procedure might be confusing, but it will be improved in the future, for example, by introducing a user-friendly background generator." ] }, { @@ -715,15 +712,15 @@ "source": [ "# 2. Calculate the coordinate conversion matrix\n", "\n", - "Here we calculate the dwell time map on each sky pixel and each time bin, and then combine them as a coordinate conversion matrix (ccm).\n", - "\n", "CoordsysConversionMatrix.spacecraft_time_binning_ccm can produce the ccm for the time binning.\n", "\n", - "The ccm $C^{lb, \\nu\\lambda}_{t}$ is a three-dimentional matrix with the axes of 'lb', 'Time' and 'NuLambda'.\n", + "Here we calculate the dwell time map on each sky pixel and each time bin, and then combine them as a coordinate conversion matrix (ccm).\n", + "\n", + "The ccm $C^{lb, \\nu\\lambda}_{t}$ is a three-dimensional matrix with the axes of 'lb', 'Time' and 'NuLambda'.\n", "\n", - "$C^{lb, \\nu\\lambda}_{t}$ is the dwell time map on the local coordinate pixel $\\nu\\lambda$ for the sky pixel $lb$ (in the galactic coordinate) during the time bin $t$.\n", + "$C^{lb, \\nu\\lambda}_{t}$ is the exposure time on the pixel $\\nu\\lambda$ on the detector coordinate for the model pixel $lb$ (in the galactic coordinate) during the time bin $t$.\n", "\n", - "Effectively, $C^{lb, \\nu\\lambda}_{t}$ acts as the coordinate conversion matrix of the sky map from the galactic coordinate to the local coordinate for each time bin." + "By multiplying $C^{lb, \\nu\\lambda}_{t}$ with the model map, it can be converted into the detector coordinate for each time bin." ] }, { @@ -731,7 +728,7 @@ "id": "47b489df", "metadata": {}, "source": [ - "## Read the orientation file and extract the orientation information around the GRB interval" + "## Read the orientation file and extract the orientation information around the GRB event" ] }, { @@ -952,44 +949,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", @@ -1023,7 +1026,7 @@ "id": "b23f1fbe", "metadata": {}, "source": [ - "DataLoader is a data container for the image deconvolution. It also calculates several matrices for the analysis." + "DataLoader is a data container for the image deconvolution. It also calculates several auxiliary matrices for the analysis." ] }, { @@ -1033,9 +1036,9 @@ "source": [ "## 4-2. Load the response file\n", "\n", - "The response file will be loaded on the CPU memory. It requires a few GB. In the actuall analysis, the response will be much larger, ~TB. \n", + "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. \n", "\n", - "So loading it on the memory might be unrealistic. The optimized (lazy) loading would be a next work." + "So loading it on the memory might be unrealistic in the future. The optimized (lazy) loading would be a next work." ] }, { @@ -1078,7 +1081,7 @@ "id": "0c58e43e", "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." ] }, { @@ -1099,14 +1102,6 @@ "dataloader.calc_image_response_projected()" ] }, - { - "cell_type": "markdown", - "id": "241505ad", - "metadata": {}, - "source": [ - "(the following function is tentetive. It will be removed in the future.)" - ] - }, { "cell_type": "code", "execution_count": 23, @@ -1162,7 +1157,7 @@ "source": [ "## 4-4. 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." ] }, { @@ -1215,42 +1210,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 |" ] }, { @@ -2171,7 +2166,7 @@ "source": [ "# 5. Analyze the results\n", "\n", - "Below examples to see/analyze the results are shown." + "Examples to see/analyze the results are shown below." ] }, {