{ "cells": [ { "cell_type": "markdown", "id": "3edcfe0b-24d7-4321-b355-a6dc730c155d", "metadata": { "tags": [] }, "source": [ "# Image Deconvolution example for 511 keV with spacecraft attitude binning \n", "\n", "updated on 2025-03-19\n", "\n", "This notebook focuses on the image deconvolution with the spacecraft attitude (scatt) binning method. Using an older 511keV thin disk 3-month simulation data created for DC2, an example of the image analysis will be presented. If you have not run through 511keV-ScAtt-DataReduction.ipynb, please do that first as this analysis requires data that is binned slightly differently than the standard Galactic binning method shown in the DataIO tutorial." ] }, { "cell_type": "code", "execution_count": 1, "id": "6ac10677-610a-47f9-8915-f99850a2ae45", "metadata": {}, "outputs": [], "source": [ "import logging\n", "import sys\n", "logger = logging.getLogger('cosipy')\n", "logger.setLevel(logging.INFO)\n", "logger.addHandler(logging.StreamHandler(sys.stdout))" ] }, { "cell_type": "code", "execution_count": 2, "id": "e751bbd5", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
22:06:23 WARNING The naima package is not available. Models that depend on it will not be functions.py:48\n", " available \n", "\n" ], "text/plain": [ "\u001b[38;5;46m22:06:23\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The naima package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=454518;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=631997;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py#48\u001b\\\u001b[2m48\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:69\n", " will not be available. \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=188450;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=829818;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py#69\u001b\\\u001b[2m69\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mwill not be available. \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/numba-0.58.0-py3.10-macosx-14-arm64.egg/numba/core/decorators.py:262: NumbaDeprecationWarning: \u001b[1mnumba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.\u001b[0m\n", " warnings.warn(msg, NumbaDeprecationWarning)\n" ] }, { "data": { "text/html": [ "
WARNING The ebltable package is not available. Models that depend on it will not be absorption.py:33\n", " available \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The ebltable package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=628702;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/astromodels/functions/functions_1D/absorption.py\u001b\\\u001b[2mabsorption.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=541363;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/astromodels/functions/functions_1D/absorption.py#33\u001b\\\u001b[2m33\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/numba-0.58.0-py3.10-macosx-14-arm64.egg/numba/core/decorators.py:262: NumbaDeprecationWarning: \u001b[1mnumba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.\u001b[0m\n", " warnings.warn(msg, NumbaDeprecationWarning)\n" ] }, { "data": { "text/html": [ "
INFO Starting 3ML! __init__.py:39\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;49mINFO \u001b[0m \u001b[1;38;5;251m Starting 3ML! \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=645108;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=274244;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#39\u001b\\\u001b[2m39\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING WARNINGs here are NOT errors __init__.py:40\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m WARNINGs here are \u001b[0m\u001b[1;31mNOT\u001b[0m\u001b[1;38;5;251m errors \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=988291;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=967235;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#40\u001b\\\u001b[2m40\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING but are inform you about optional packages that can be installed __init__.py:41\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m but are inform you about optional packages that can be installed \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=916730;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=737787;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#41\u001b\\\u001b[2m41\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING to disable these messages, turn off start_warning in your config file __init__.py:44\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m \u001b[0m\u001b[1;31m to disable these messages, turn off start_warning in your config file\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=963952;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=759304;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#44\u001b\\\u001b[2m44\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
22:06:23 WARNING ROOT minimizer not available minimization.py:1345\n", "\n" ], "text/plain": [ "\u001b[38;5;46m22:06:23\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m ROOT minimizer not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=689820;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=787427;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/minimizer/minimization.py#1345\u001b\\\u001b[2m1345\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Multinest minimizer not available minimization.py:1357\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Multinest minimizer not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=875996;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=715399;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/minimizer/minimization.py#1357\u001b\\\u001b[2m1357\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING PyGMO is not available minimization.py:1369\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m PyGMO is not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=996386;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=901169;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/minimizer/minimization.py#1369\u001b\\\u001b[2m1369\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING The cthreeML package is not installed. You will not be able to use plugins which __init__.py:94\n", " require the C/C++ interface (currently HAWC) \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The cthreeML package is not installed. You will not be able to use plugins which \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=944832;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=955109;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#94\u001b\\\u001b[2m94\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mrequire the C/C++ interface \u001b[0m\u001b[1;38;5;251m(\u001b[0m\u001b[1;38;5;251mcurrently HAWC\u001b[0m\u001b[1;38;5;251m)\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Could not import plugin HAWCLike.py. Do you have the relative instrument __init__.py:144\n", " software installed and configured? \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin HAWCLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=803280;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=300317;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Could not import plugin FermiLATLike.py. Do you have the relative instrument __init__.py:144\n", " software installed and configured? \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin FermiLATLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=752487;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=576218;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING No fermitools installed lat_transient_builder.py:44\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m No fermitools installed \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=905500;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/utils/data_builders/fermi/lat_transient_builder.py\u001b\\\u001b[2mlat_transient_builder.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=845496;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/utils/data_builders/fermi/lat_transient_builder.py#44\u001b\\\u001b[2m44\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n", " performances in 3ML \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable OMP_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=620676;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=51629;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n", " performances in 3ML \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable MKL_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=44506;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=924987;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n", " performances in 3ML \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=897815;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=576834;file:///Users/ckierans/Software/COSItools/COSItools/python-env/lib/python3.10/site-packages/threeML-2.4.3.dev1-py3.10.egg/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from histpy import Histogram, HealpixAxis, Axis, Axes\n", "from mhealpy import HealpixMap\n", "from astropy.coordinates import SkyCoord, cartesian_to_spherical, Galactic\n", "\n", "from cosipy.response import FullDetectorResponse\n", "from cosipy.spacecraftfile import SpacecraftFile\n", "from cosipy.ts_map.TSMap import TSMap\n", "from cosipy.data_io import UnBinnedData, BinnedData\n", "from cosipy.image_deconvolution import SpacecraftAttitudeExposureTable, CoordsysConversionMatrix, DataIF_COSI_DC2, ImageDeconvolution\n", "\n", "# cosipy uses astropy units\n", "import astropy.units as u\n", "from astropy.units import Quantity\n", "from astropy.coordinates import SkyCoord\n", "from astropy.time import Time\n", "from astropy.table import Table\n", "from astropy.io import fits\n", "from scoords import Attitude, SpacecraftFrame\n", "\n", "#3ML is needed for spectral modeling\n", "from threeML import *\n", "from astromodels import Band\n", "\n", "#Other standard libraries\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib.gridspec import GridSpec \n", "\n", "import healpy as hp\n", "from tqdm.autonotebook import tqdm\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "id": "00f20cda-81f8-4685-b9c4-f9423e5ebcf7", "metadata": { "tags": [] }, "source": [ "## 0. Files needed for this notebook\n", "\n", "You should have used the 511keV-ScAtt-DataReduction.ipynb notebook to download and bin the data sets needed for this notebook. Specicially, you now should have:\n", "- 511keV_scatt_binning_DC2_bkg.hdf5\n", "- 511keV_scatt_binning_DC2_event.hdf5\n", "- ccm.hdf5\n", "\n", "From docs/tutorials/image_deconvolution/511keV/ScAttBinning, you will also need:\n", "- inputs_511keV.yaml\n", "- imagedeconvolution_parfile_scatt_511keV.yml\n" ] }, { "cell_type": "markdown", "id": "6c259412", "metadata": {}, "source": [ "## 1. Read the response matrix" ] }, { "cell_type": "markdown", "id": "573a7c60", "metadata": {}, "source": [ " please modify \"path_data\" corresponding to your environment." ] }, { "cell_type": "code", "execution_count": null, "id": "6c1dcdb3-0ed2-4ec6-8d47-7a629a7bd6ec", "metadata": {}, "outputs": [], "source": [ "path_data = \"/path/to/data/\"" ] }, { "cell_type": "code", "execution_count": 4, "id": "98a778c2-73cf-467b-96b6-affc42f17102", "metadata": {}, "outputs": [], "source": [ "response_path = path_data + \"SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5\"\n", "\n", "response = FullDetectorResponse.open(response_path)" ] }, { "cell_type": "code", "execution_count": 5, "id": "eab660b4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FILENAME: '/Users/ckierans/Software/COSItools/COSItools/cosipy/docs/tutorials/image_deconvolution/511keV/ScAttBinning/SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5'\n", "AXES:\n", " NuLambda:\n", " DESCRIPTION: 'Location of the simulated source in the spacecraft coordinates'\n", " TYPE: 'healpix'\n", " NPIX: 3072\n", " NSIDE: 16\n", " SCHEME: 'RING'\n", " Ei:\n", " DESCRIPTION: 'Initial simulated energy'\n", " TYPE: 'log'\n", " UNIT: 'keV'\n", " NBINS: 1\n", " EDGES: [509.0 keV, 513.0 keV]\n", " Em:\n", " DESCRIPTION: 'Measured energy'\n", " TYPE: 'log'\n", " UNIT: 'keV'\n", " NBINS: 1\n", " EDGES: [509.0 keV, 513.0 keV]\n", " Phi:\n", " DESCRIPTION: 'Compton angle'\n", " TYPE: 'linear'\n", " UNIT: 'deg'\n", " NBINS: 60\n", " EDGES: [0.0 deg, 3.0 deg, 6.0 deg, 9.0 deg, 12.0 deg, 15.0 deg, 18.0 deg, 21.0 deg, 24.0 deg, 27.0 deg, 30.0 deg, 33.0 deg, 36.0 deg, 39.0 deg, 42.0 deg, 45.0 deg, 48.0 deg, 51.0 deg, 54.0 deg, 57.0 deg, 60.0 deg, 63.0 deg, 66.0 deg, 69.0 deg, 72.0 deg, 75.0 deg, 78.0 deg, 81.0 deg, 84.0 deg, 87.0 deg, 90.0 deg, 93.0 deg, 96.0 deg, 99.0 deg, 102.0 deg, 105.0 deg, 108.0 deg, 111.0 deg, 114.0 deg, 117.0 deg, 120.0 deg, 123.0 deg, 126.0 deg, 129.0 deg, 132.0 deg, 135.0 deg, 138.0 deg, 141.0 deg, 144.0 deg, 147.0 deg, 150.0 deg, 153.0 deg, 156.0 deg, 159.0 deg, 162.0 deg, 165.0 deg, 168.0 deg, 171.0 deg, 174.0 deg, 177.0 deg, 180.0 deg]\n", " PsiChi:\n", " DESCRIPTION: 'Location in the Compton Data Space'\n", " TYPE: 'healpix'\n", " NPIX: 3072\n", " NSIDE: 16\n", " SCHEME: 'RING'\n" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "response" ] }, { "cell_type": "markdown", "id": "26d6eb3a", "metadata": {}, "source": [ "## 2. Read binned 511keV binned files (source and background)" ] }, { "cell_type": "code", "execution_count": 6, "id": "04e15347-6b38-42de-a7c5-cd99b2ae66ac", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 117 ms, sys: 671 ms, total: 788 ms\n", "Wall time: 809 ms\n" ] } ], "source": [ "%%time\n", "\n", "# background \n", "data_bkg = BinnedData(\"inputs_511keV.yaml\")\n", "data_bkg.load_binned_data_from_hdf5(\"511keV_dc2_scatt_bkg.hdf5\")\n", "\n", "## signal + background\n", "data_511keV = BinnedData(\"inputs_511keV.yaml\")\n", "data_511keV.load_binned_data_from_hdf5(\"511keV_dc2_scatt_event.hdf5\")" ] }, { "cell_type": "markdown", "id": "a409aa7b-9bd8-443b-be46-ee5a053f8349", "metadata": { "tags": [] }, "source": [ "## 3. Load the coordsys conversion matrix" ] }, { "cell_type": "code", "execution_count": 7, "id": "daaf836a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.21 s, sys: 64 ms, total: 1.28 s\n", "Wall time: 1.28 s\n" ] } ], "source": [ "%%time\n", "\n", "ccm = CoordsysConversionMatrix.open(\"ccm_dc2.hdf5\")" ] }, { "cell_type": "markdown", "id": "31ec05ad-90b7-4fad-9ad0-98cfd6483d41", "metadata": {}, "source": [ "## 4. Imaging deconvolution" ] }, { "cell_type": "markdown", "id": "6e88ca7f", "metadata": {}, "source": [ "### Brief overview of the image deconvolution\n", "\n", "Basically, we have to maximize the following likelihood function\n", "\n", "$$\n", "\\log L = \\sum_i X_i \\log \\epsilon_i - \\sum_i \\epsilon_i\n", "$$\n", "\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 (basically gamma-ray flux at $j$-th pixel)\n", "\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 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", "$$\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 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", "$$\n", "$$\n", "\\alpha^{k} < \\mathrm{max}(- \\lambda_{j}^{k} / \\delta \\lambda_{j}^{k})\n", "$$\n", "\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_{\\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", "$$\n", "$$\n", "w_j = \\left(\\sum_{i} R_{ij}\\right)^\\beta\n", "$$\n", "\n", "$\\left[ ... \\right]_{\\mathrm{gauss}}$ means that the differential image is smoothed by a gaussian filter." ] }, { "cell_type": "markdown", "id": "e0a2582e", "metadata": {}, "source": [ "## 4-1. Prepare DataInterface containing all necessary datasets" ] }, { "cell_type": "code", "execution_count": 8, "id": "de8055f7-4aab-4a17-8751-42493f9e88d6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Loading the response matrix onto your computer memory...\n", "Finished\n", "... checking the axis ScAtt of the event and background files...\n", " --> pass (edges)\n", "... checking the axis Em of the event and background files...\n", " --> pass (edges)\n", "... checking the axis Phi of the event and background files...\n", " --> pass (edges)\n", "... checking the axis PsiChi of the event and background files...\n", " --> pass (edges)\n", "...checking the axis Em of the event and response files...\n", " --> pass (edges)\n", "...checking the axis Phi of the event and response files...\n", " --> pass (edges)\n", "...checking the axis PsiChi of the event and response files...\n", " --> pass (edges)\n", "The axes in the event and background files are redefined. Now they are consistent with those of the response file.\n", "Calculating an exposure map...\n", "Finished...\n", "CPU times: user 6.67 s, sys: 1.49 s, total: 8.16 s\n", "Wall time: 8.18 s\n" ] } ], "source": [ "%%time\n", "\n", "data_interface = DataIF_COSI_DC2.load(name = \"511keV\",\n", " event_binned_data = data_511keV.binned_data,\n", " dict_bkg_binned_data = {\"albedo\": data_bkg.binned_data},\n", " rsp = response,\n", " coordsys_conv_matrix=ccm,\n", " is_miniDC2_format=False)" ] }, { "cell_type": "markdown", "id": "b1a0269e", "metadata": {}, "source": [ "### 4-2. Initialize the instance of the image deconvolution class\n", "\n", "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." ] }, { "cell_type": "markdown", "id": "79eb910c", "metadata": {}, "source": [ " please modify this parameter_filepath corresponding to your environment." ] }, { "cell_type": "code", "execution_count": 9, "id": "5fa73486", "metadata": {}, "outputs": [], "source": [ "parameter_filepath = \"imagedeconvolution_parfile_scatt_511keV.yml\"" ] }, { "cell_type": "code", "execution_count": 10, "id": "a4b47308-3e13-400d-bebc-b5d1e093201d", "metadata": {}, "outputs": [], "source": [ "image_deconvolution = ImageDeconvolution()\n", "\n", "# set data_interface to image_deconvolution\n", "image_deconvolution.set_dataset([data_interface])\n", "\n", "# set a parameter file for the image deconvolution\n", "image_deconvolution.read_parameterfile(parameter_filepath)" ] }, { "cell_type": "markdown", "id": "a2345d9d", "metadata": {}, "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, e.g., the expected counts with the initial model map, gaussian smoothing filter etc.\n", "\n", "I describe parameters in the parameter file.\n", "\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 | 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", "#### 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 length of \"energy_edges\" - 1 |\n", "\n", "#### deconvolution\n", "\n", "| Name | Unit | Description | Notes |\n", "| :---: | :---: | :---: | :---: |\n", "|algorithm | str | the name of the image deconvolution algorithm| As for now, only 'RL' is supported |\n", "|||||\n", "|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 |" ] }, { "cell_type": "code", "execution_count": 11, "id": "879053e3-ac7b-4a0a-ad58-24e3fb137065", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#### Initialization Starts ####\n", "<< Instantiating the model class AllSkyImage >>\n", "---- parameters ----\n", "coordinate: galactic\n", "energy_edges:\n", " unit: keV\n", " value:\n", " - 509.0\n", " - 513.0\n", "nside: 16\n", "scheme: ring\n", "unit: cm-2 s-1 sr-1\n", "\n", "<< Setting initial values of the created model object >>\n", "---- parameters ----\n", "algorithm: flat\n", "parameter:\n", " unit: cm-2 s-1 sr-1\n", " value:\n", " - 1e-4\n", "\n", "<< Registering the deconvolution algorithm >>\n", "Gaussian filter with FWHM of 2.0 deg will be applied to delta images ...\n", "---- parameters ----\n", "algorithm: RL\n", "parameter:\n", " acceleration: true\n", " alpha_max: 10.0\n", " background_normalization_optimization: true\n", " background_normalization_range:\n", " albedo:\n", " - 0.01\n", " - 10.0\n", " iteration_max: 10\n", " response_weighting: true\n", " response_weighting_index: 0.5\n", " save_results: false\n", " save_results_directory: ./results\n", " smoothing: true\n", " smoothing_FWHM:\n", " unit: deg\n", " value: 2.0\n", "\n", "#### Initialization Finished ####\n" ] } ], "source": [ "image_deconvolution.initialize()" ] }, { "cell_type": "markdown", "id": "5fa48a9c", "metadata": {}, "source": [ "### (You can change the parameters as follows)\n", "\n", "Note that when you modify the parameters, do not forget to run \"initialize\" again!" ] }, { "cell_type": "code", "execution_count": 12, "id": "1a658d2a-4dee-4d05-83ae-d7ac06317c73", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#### Initialization Starts ####\n", "<< Instantiating the model class AllSkyImage >>\n", "---- parameters ----\n", "coordinate: galactic\n", "energy_edges:\n", " unit: keV\n", " value:\n", " - 509.0\n", " - 513.0\n", "nside: 16\n", "scheme: ring\n", "unit: cm-2 s-1 sr-1\n", "\n", "<< Setting initial values of the created model object >>\n", "---- parameters ----\n", "algorithm: flat\n", "parameter:\n", " unit: cm-2 s-1 sr-1\n", " value:\n", " - 1e-4\n", "\n", "<< Registering the deconvolution algorithm >>\n", "Gaussian filter with FWHM of 3.0 deg will be applied to delta images ...\n", "---- parameters ----\n", "algorithm: RL\n", "parameter:\n", " acceleration: true\n", " alpha_max: 10\n", " background_normalization_optimization: true\n", " background_normalization_range:\n", " albedo:\n", " - 0.01\n", " - 10.0\n", " iteration_max: 30\n", " response_weighting: true\n", " response_weighting_index: 0.5\n", " save_results: false\n", " save_results_directory: ./results\n", " smoothing: true\n", " smoothing_FWHM:\n", " unit: deg\n", " value: 3.0\n", "\n", "#### Initialization Finished ####\n" ] } ], "source": [ "image_deconvolution.override_parameter(\"deconvolution:parameter:iteration_max = 30\")\n", "image_deconvolution.override_parameter(\"deconvolution:parameter:background_normalization_optimization = True\")\n", "image_deconvolution.override_parameter(\"deconvolution:parameter:alpha_max = 10\")\n", "image_deconvolution.override_parameter(\"deconvolution:parameter:smoothing_FWHM:value = 3.0\")\n", "\n", "image_deconvolution.initialize()" ] }, { "cell_type": "markdown", "id": "f764066e", "metadata": {}, "source": [ "## 4-3. Start the image deconvolution\n", "\n", "**On a MacBook Pro with M2 Max and 96 GB memory, it takes about 2.5 min for 30 iterations.**" ] }, { "cell_type": "code", "execution_count": 13, "id": "a57fbf71-2fcc-48c4-9ac7-4c545dca67c9", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#### Image Deconvolution Starts ####\n", "<< Initialization >>\n", "The response weighting filter was calculated.\n", "The expected count histograms were calculated with the initial model map.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "1f422577bfd94596957ee4fe6cbd7a83", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/30 [00:00, ?it/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "## Iteration 1/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 2.4783944324338596\n", " background_normalization: {'albedo': 1.0003541905620614}\n", " loglikelihood: [-1197056.4545480786]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 2/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.575421707005236\n", " background_normalization: {'albedo': 0.9922607778263346}\n", " loglikelihood: [-1148892.2494289856]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 3/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 10\n", " background_normalization: {'albedo': 0.9989836398988254}\n", " loglikelihood: [-1366763.493999292]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 4/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9807935141594855}\n", " loglikelihood: [-1130201.6042120028]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 5/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9951524418084393}\n", " loglikelihood: [-1124841.8752257773]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 6/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.996710972292622}\n", " loglikelihood: [-1123882.6965031484]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 7/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9968855622562821}\n", " loglikelihood: [-1123152.1905212747]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 8/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969083706965138}\n", " loglikelihood: [-1122553.294395423]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 9/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969134554559144}\n", " loglikelihood: [-1122056.1156774329]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 10/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969163516897256}\n", " loglikelihood: [-1121639.152483511]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 11/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969189620579678}\n", " loglikelihood: [-1121286.2048804823]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 12/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969215191961452}\n", " loglikelihood: [-1120984.8976566715]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 13/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.996924030389004}\n", " loglikelihood: [-1120725.6564115596]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 14/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969264790341887}\n", " loglikelihood: [-1120500.9815390154]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 15/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969288556146668}\n", " loglikelihood: [-1120304.963459809]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 16/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969311484060269}\n", " loglikelihood: [-1120132.8961275548]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 17/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969333516388736}\n", " loglikelihood: [-1119980.969809744]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 18/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969354690483336}\n", " loglikelihood: [-1119846.1195404662]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 19/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969374982210912}\n", " loglikelihood: [-1119725.820978066]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 20/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969394510025879}\n", " loglikelihood: [-1119618.0199778753]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 21/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969413216933766}\n", " loglikelihood: [-1119520.9964885826]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 22/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Post-processing >>\n", "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969431196567854}\n", " loglikelihood: [-1119433.3289185578]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 23/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969448441589266}\n", " loglikelihood: [-1119353.8044356154]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 24/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969465063444388}\n", " loglikelihood: [-1119281.4278713942]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 25/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969481006699753}\n", " loglikelihood: [-1119215.3384681274]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 26/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969496370088307}\n", " loglikelihood: [-1119154.8018173785]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 27/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969511182269674}\n", " loglikelihood: [-1119099.190130895]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 28/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969525475454646}\n", " loglikelihood: [-1119047.9797920918]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 29/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969539225343819}\n", " loglikelihood: [-1119000.7175890533]\n", "<< Checking Stopping Criteria >>\n", "--> Continue\n", "## Iteration 30/30 ##\n", "<< Pre-processing >>\n", "<< E-step >>\n", "<< M-step >>\n", "<< Post-processing >>\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "WARNING RuntimeWarning: invalid value encountered in divide\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "<< Registering Result >>\n", " alpha: 1.0\n", " background_normalization: {'albedo': 0.9969552449401047}\n", " loglikelihood: [-1118956.9770636673]\n", "<< Checking Stopping Criteria >>\n", "--> Stop\n", "<< Finalization >>\n", "#### Image Deconvolution Finished ####\n", "CPU times: user 18min 4s, sys: 3min 54s, total: 21min 58s\n", "Wall time: 2min 41s\n" ] } ], "source": [ "%%time\n", "\n", "image_deconvolution.run_deconvolution()" ] }, { "cell_type": "code", "execution_count": 14, "id": "cc64ea8d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[{'iteration': 1,\n", " 'model':