Implementation of the KDE/photospline PSF (9yr/11yr NT samples) as sp…

…atial PDF for flarestack (#333) * implementation of KDE generated MC PDF in flarestack * add conda environment file to install photospline * display conda and python version in separate step * use photospline v2.1.0 * dont activate base env * try activating environment manually * specify shell to correctly activate conda env * implementation now supports gamma-dependence of KDE PSF * fix mypy / black issues * cleanup benchmark analysis files * add photospline installation instructions * cleanup benchmark scripts * try pinning urllib 2.0.6 --------- Co-authored-by: JannisNe <necker@physik.hu-berlin.de>
icecube · Nov 30, 2023 · 4f56077 · 4f56077
1 parent 956b8ca
commit 4f56077
Show file tree

Hide file tree

Showing 16 changed files with 1,733 additions and 935 deletions.
diff --git a/.github/workflows/continous_integration.yml b/.github/workflows/continous_integration.yml
@@ -33,27 +33,44 @@ jobs:
       # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
       - uses: actions/checkout@v3
 
-      # Set up the python versions
-      - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+      # Set up the python versions and install conda dependencies
+      - name: Install conda dependencies
+        uses: conda-incubator/setup-miniconda@v2
         with:
           python-version: ${{ matrix.python-version }}
+          environment-file: conda_env.yml
+          activate-environment: flarestack
+          auto-activate-base: false
+
+      # Display the python version and conda info
+      - name: Display Python version and Conda info
+        shell: bash -el {0}
+        run: |
+            python --version
+            conda --version
+            conda config --show-sources
+            conda config --show
+            conda list
+            conda info -a
 
      # Setting up dependencies
       - name: Install dependencies
+        shell: bash -el {0}
         run: |
           python -m pip install --upgrade pip
           python -m pip install poetry
           python -m poetry install 
 
       # Runs a set of commands using the runners shell
       - name: Test the code
+        shell: bash -el {0}
         run: |
           poetry run coverage run --concurrency=multiprocessing -m unittest discover tests/
           poetry run coverage combine
 
       - name: Run Coveralls
         if: ${{ success() }}
+        shell: bash -el {0}
         run: |
           poetry run coveralls --service=github
         env:
@@ -63,10 +80,12 @@ jobs:
         run: echo "Tag is ${{ github.ref }}, Deploy is ${{ startsWith(github.ref, 'refs/tags/') && matrix.python-version == 3.9}}"
 
       - name: Build a binary wheel and a source tarball
+        shell: bash -el {0}
         run: python -m poetry build
 
       - name: Publish distribution 📦 to PyPI
         if: ${{ startsWith(github.ref, 'refs/tags/') && success() && matrix.python-version == 3.9}}
         env:
           POETRY_PYPI_TOKEN_PYPI: ${{ secrets.PYPI_API_TOKEN }}
+        shell: bash -el {0}
         run: poetry publish
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,4 @@
 flarestack.egg-info
 *.pyc
-docs/build
+docs/build
+scratch/
diff --git a/README.md b/README.md
@@ -76,6 +76,18 @@ yourself, and still enable you to import flarestack.
 If you do want to contribute to _flarestack_, you can check out some guidelines [here](https://github.com/icecube/flarestack/blob/master/CONTRIBUTING.md).
 
 
+### Using KDE PDFs
+If you want to use the KDE spatial PDFs you need to install `photospline` using `conda`.
+`flarestack` ships with a `conda_env.yml` file that pins the correct `photospline` version. 
+Find your `flarestack` directory (either in your `site-packages` directory if you followed Option A or in your `git` clone if you followed Option B) and execute:
+
+```shell
+conda env create -f conda_env.yml
+```
+
+Alternatively, You can try and install the latest `photospline` version with as described [here](https://anaconda.org/conda-forge/photospline).
+
+
 ### Right, anyway, I've now downloaded *flarestack*. Can I use it right away?
 
 You can get started with *flarestack* immediately using public IceCube datasets provided as part of the code. You can simply run scripts such as those under `/flarestack/analyses/`, and do your science!

diff --git a/conda_env.yml b/conda_env.yml
@@ -0,0 +1,6 @@
+name: flarestack
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - photospline=2.1.0
diff --git a/docs/source/setup.md b/docs/source/setup.md
@@ -60,6 +60,18 @@ yourself, and still enable you to import flarestack.
 If you do want to contribute to _flarestack_, you can check out some guidelines [here](https://github.com/icecube/flarestack/blob/master/CONTRIBUTING.md).
 
 
+### Using KDE PDFs
+If you want to use the KDE spatial PDFs you need to install `photospline` using `conda`.
+`flarestack` ships with a `conda_env.yml` file that pins the correct `photospline` version. 
+Find your `flarestack` directory (either in your `site-packages` directory if you followed Option A or in your `git` clone if you followed Option B) and execute:
+
+```shell
+conda env create -f conda_env.yml
+```
+
+Alternatively, You can try and install the latest `photospline` version with as described [here](https://anaconda.org/conda-forge/photospline).
+
+
 ### Right, anyway, I've now downloaded *flarestack*. Can I use it right away?
 
 You can get started with *flarestack* immediately using public IceCube datasets provided as part of the code. You can simply run scripts such as those under `/flarestack/analyses/`, and do your science!

diff --git a/flarestack/analyses/benchmarks/benchmark_point_source_sensitivity_10yr.py b/flarestack/analyses/benchmarks/benchmark_point_source_sensitivity_10yr.py
@@ -24,28 +24,32 @@
     "time_pdf_name": "steady",
 }
 
-llh_time = {
-    "time_pdf_name": "steady",
+injection_spatial = {
+    "spatial_pdf_name": "circular_gaussian",
 }
 
 inj_kwargs = {
     "injection_energy_pdf": injection_energy,
     "injection_sig_time_pdf": injection_time,
+    "injection_spatial_pdf": injection_spatial,
 }
 
+
+llh_time = injection_time
+llh_spatial = injection_spatial
 llh_energy = injection_energy
 
 llh_kwargs = {
     "llh_name": "standard",
     "llh_energy_pdf": llh_energy,
+    "llh_spatial_pdf": llh_spatial,
     "llh_sig_time_pdf": llh_time,
     "llh_bkg_time_pdf": {"time_pdf_name": "steady"},
     "negative_ns_bool": True,
 }
 
 name = "analyses/benchmarks/ps_sens_10yr"
 
-# sindecs = np.linspace(0.90, -0.90, 3)
 sindecs = np.linspace(0.90, -0.90, 19)
 # sindecs = np.linspace(0.5, -0.5, 3)
 #
@@ -60,7 +64,7 @@
 
     subname = name + "/sindec=" + "{0:.2f}".format(sindec) + "/"
 
-    scale = flux_to_k(reference_sensitivity(sindec)) * 3
+    scale = (flux_to_k(reference_sensitivity(sindec)) * 3)[0]
 
     mh_dict = {
         "name": subname,
@@ -70,14 +74,15 @@
         "inj_dict": inj_kwargs,
         "llh_dict": llh_kwargs,
         "scale": scale,
-        "n_trials": 5000,
+        "n_trials": 1000,
+        "trials_per_task": 50,
         "n_steps": 15,
     }
 
     job_id = analyse(
         mh_dict,
         cluster=cluster,
-        n_cpu=1 if cluster else 32,
+        n_cpu=1 if cluster else 16,
         h_cpu="23:59:59",
         ram_per_core="8.0G",
     )
@@ -136,7 +141,7 @@
 ax1.set_xlim(xmin=-1.0, xmax=1.0)
 # ax1.set_ylim(ymin=1.e-13, ymax=1.e-10)
 ax1.grid(True, which="both")
-ax1.semilogy(nonposy="clip")
+ax1.semilogy(nonpositive="clip")
 ax1.set_ylabel(r"Flux Strength [ GeV$^{-1}$ cm$^{-2}$ s$^{-1}$ ]", fontsize=12)
 
 ax2 = plt.subplot2grid((4, 1), (3, 0), colspan=3, rowspan=1, sharex=ax1)

diff --git a/flarestack/analyses/benchmarks/benchmark_point_source_sensitivity_11yr_NT_KDE+gamma.py b/flarestack/analyses/benchmarks/benchmark_point_source_sensitivity_11yr_NT_KDE+gamma.py
@@ -0,0 +1,195 @@
+import numpy as np
+from flarestack import ResultsHandler, MinimisationHandler
+from flarestack.data.icecube import nt_v005_p01
+from flarestack.shared import plot_output_dir, flux_to_k
+from flarestack.utils.prepare_catalogue import ps_catalogue_name
+from flarestack.icecube_utils.reference_sensitivity import (
+    reference_sensitivity,
+    reference_discovery_potential,
+)
+import matplotlib.pyplot as plt
+from flarestack import analyse, wait_cluster
+import logging
+import sys, pickle
+
+logging.basicConfig(level=logging.INFO)
+
+# Initialise Injectors/LLHs
+
+injection_energy = {
+    "energy_pdf_name": "power_law",
+    "gamma": 2.0,
+}
+
+injection_time = {
+    "time_pdf_name": "steady",
+}
+
+injection_spatial = {
+    "spatial_pdf_name": "northern_tracks_kde",
+    "spatial_pdf_data": "/lustre/fs22/group/icecube/data_mirror/northern_tracks/version-005-p01/KDE_PDFs_v007/IC86_pass2/sig_E_psi_photospline_v006_4D.fits",
+}
+
+inj_kwargs = {
+    "injection_energy_pdf": injection_energy,
+    "injection_sig_time_pdf": injection_time,
+    "injection_spatial_pdf": injection_spatial,
+}
+
+
+llh_time = injection_time
+llh_spatial = injection_spatial
+llh_energy = injection_energy
+
+llh_kwargs = {
+    "llh_name": "standard_kde_enabled",
+    "llh_energy_pdf": llh_energy,
+    "llh_spatial_pdf": llh_spatial,
+    "llh_sig_time_pdf": llh_time,
+    "llh_bkg_time_pdf": {"time_pdf_name": "steady"},
+    "negative_ns_bool": True,
+}
+
+name = "analyses/benchmarks/ps_sens_10yr_NT_KDE+gamma_v2"
+
+# sindecs = np.linspace(0.90, -0.90, 3)
+sindecs = np.linspace(0.95, -0.05, 11)
+# sindecs = np.linspace(0.5, -0.5, 3)
+#
+analyses = []
+
+cluster = True
+
+job_ids = []
+
+for sindec in sindecs:
+    cat_path = ps_catalogue_name(sindec)
+
+    subname = name + "/sindec=" + "{0:.2f}".format(sindec) + "/"
+
+    scale = (flux_to_k(reference_sensitivity(sindec)) * 2)[0]
+
+    mh_dict = {
+        "name": subname,
+        "mh_name": "fixed_weights",
+        "dataset": nt_v005_p01,
+        "catalogue": cat_path,
+        "inj_dict": inj_kwargs,
+        "llh_dict": llh_kwargs,
+        "scale": scale,
+        "n_trials": 2000,
+        "n_steps": 15,
+    }
+
+    if len(sys.argv) > 1 and sys.argv[1] == "--run-jobs":
+        job_id = analyse(
+            mh_dict,
+            cluster=cluster,
+            n_cpu=1 if cluster else 16,
+            h_cpu="23:59:59",
+            ram_per_core="8.0G",
+            remove_old_logs=False,
+            trials_per_task=50,
+        )
+        job_ids.append(job_id)
+
+    analyses.append(mh_dict)
+
+if len(sys.argv) > 1 and sys.argv[1] == "--run-jobs":
+    wait_cluster(job_ids)
+
+
+sens = []
+sens_err = []
+disc_pots = []
+
+for rh_dict in analyses:
+    rh = ResultsHandler(rh_dict)
+    sens.append(rh.sensitivity)
+    sens_err.append(rh.sensitivity_err)
+    disc_pots.append(rh.disc_potential)
+
+with open(plot_output_dir(name) + "/PS10yrKDE+gamma.pkl", "wb") as pf:
+    pickle.dump(
+        {
+            "sin_dec": sindecs,
+            "sensitivity": list(sens),
+            "sens_error": [list(x) for x in list(sens_err)],
+            "disc_potential": list(disc_pots),
+        },
+        pf,
+    )
+
+sens_err = np.array(sens_err).T
+
+# sens = reference_sensitivity(sindecs, sample="7yr")
+# disc_pots = reference_discovery_potential(sindecs, sample="7yr")
+# sens_err = 0.1*sens
+
+plot_range = np.linspace(-0.99, 0.99, 1000)
+
+plt.figure()
+ax1 = plt.subplot2grid((4, 1), (0, 0), colspan=3, rowspan=3)
+ax1.plot(
+    sindecs,
+    reference_sensitivity(sindecs, sample="10yr"),
+    color="blue",
+    label=r"10-year Point Source analysis",
+)
+# ax1.plot(sindecs, reference_sensitivity(sindecs, sample="7yr"), color="green",
+#          label=r"7-year Point Source analysis")
+
+# ax1.plot(sindecs, sens, color='orange', label="Flarestack")
+
+ax1.errorbar(
+    sindecs, sens, yerr=sens_err, color="orange", label="Sensitivity", marker="o"
+)
+
+ax1.plot(
+    sindecs,
+    reference_discovery_potential(sindecs, sample="10yr"),
+    color="blue",
+    linestyle="--",
+)
+
+ax1.plot(
+    sindecs, disc_pots, color="orange", linestyle="--", label="Discovery Potential"
+)
+
+ax1.set_xlim(xmin=-1.0, xmax=1.0)
+# ax1.set_ylim(ymin=1.e-13, ymax=1.e-10)
+ax1.grid(True, which="both")
+ax1.semilogy(nonpositive="clip")
+ax1.set_ylabel(r"Flux Strength [ GeV$^{-1}$ cm$^{-2}$ s$^{-1}$ ]", fontsize=12)
+
+ax2 = plt.subplot2grid((4, 1), (3, 0), colspan=3, rowspan=1, sharex=ax1)
+
+sens_ratios = np.array(sens) / reference_sensitivity(sindecs, sample="10yr")
+sens_ratio_errs = sens_err / reference_sensitivity(sindecs, sample="10yr")
+
+disc_ratios = np.array(disc_pots) / reference_discovery_potential(
+    sindecs, sample="10yr"
+)
+
+ax2.errorbar(sindecs, sens_ratios, yerr=sens_ratio_errs, color="red", marker="o")
+ax2.scatter(sindecs, disc_ratios, color="k")
+ax2.plot(sindecs, disc_ratios, color="k", linestyle="--")
+ax2.set_ylabel(r"ratio", fontsize=12)
+ax2.set_xlabel(r"sin($\delta$)", fontsize=12)
+
+plt.suptitle("Point Source Sensitivity (10 year)")
+#
+ax1.set_xlim(xmin=-1.0, xmax=1.0)
+
+xticklabels = ax1.get_xticklabels()
+plt.setp(xticklabels, visible=False)
+plt.subplots_adjust(hspace=0.001)
+
+ax1.legend(loc="upper right", fancybox=True, framealpha=1.0)
+
+savefile = plot_output_dir(name) + "/PS10yrKDE+gamma.pdf"
+
+logging.info(f"Saving to {savefile}")
+
+plt.savefig(savefile)
+plt.close()