.. _signal-subtraction:

Signal-subtraction likelihoods and the marginalized King PDF
=============================================================

Point-source analyses in IceCube and similar experiments commonly use a
actual detector data to model their background PDFs. When this is done, the
background PDF is effectively a mixture of the true background and the signal
from the source(s) under study. This can bias the likelihood ratio test statistic
and the resulting source significance. To mitigate this, it is possible reformulate
the likelihood, explicitly recognizing the signal contribution to the background PDF.
The result is called the *signal-subtraction* (or *background-marginalized*) likelihood.

In IceCube, the typical background spatial PDF is generated by histogramming and
normalizing the detector data as a function of reconstructed declination. To
calculate the corresponding signal contribution, it is necessary to evaluate
the King PDF integrated over right ascension — the probability of observing a
reconstructed declination :math:`\delta_\mathrm{reco}` given a true source
declination :math:`\delta_\mathrm{true}`:

.. math::

   \mathcal{M}(\delta_\mathrm{reco} \mid \delta_\mathrm{true}, \alpha, \beta)
   = \int_0^{2\pi} f\!\left(\psi(\Delta\mathrm{RA}, \delta_\mathrm{reco},
     \delta_\mathrm{true}) \;\Big|\; \alpha, \beta\right) d(\Delta\mathrm{RA})

where :math:`\psi` is the angular distance between the reconstructed event
position and the source, and :math:`f` is the King PDF. Because the King
distribution depends only on the total angular offset :math:`\psi`, the
integrand is symmetric under :math:`\Delta\mathrm{RA} \to 2\pi -
\Delta\mathrm{RA}`, so the integral reduces to

.. math::

   \mathcal{M}(\delta_\mathrm{reco} \mid \delta_\mathrm{true}, \alpha, \beta)
   = 2 \int_0^{\pi} f\!\left(\psi(\Delta\mathrm{RA}, \delta_\mathrm{reco},
     \delta_\mathrm{true}) \;\Big|\; \alpha, \beta\right) d(\Delta\mathrm{RA}).

In practice, the angular cutoff is often set to values much smaller than :math:`\pi`,
so most :math:`(\delta_\mathrm{reco},\, \delta_\mathrm{true})` pairs yield
:math:`\mathcal{M} = 0`. This makes the resulting
:math:`(N_\mathrm{events} \times N_\mathrm{sources})` matrix highly sparse.

The :class:`~kingmaker.pdf.MarginalizedKingPDF` class evaluates this integral by
pre-building a four-dimensional grid over
:math:`(\delta_\mathrm{source},\, \log_{10}\alpha,\, \beta,\,
\delta_\mathrm{evt} - \delta_\mathrm{source})` at construction time, then
uses linear interpolation across parameters at runtime to return a
:class:`scipy.sparse.csr_array` of shape
:math:`(N_\mathrm{events},\, N_\mathrm{sources})`.

.. _signal-subtraction-instantiation:

Instantiation
-------------

Construct a :class:`~kingmaker.pdf.MarginalizedKingPDF` with the source
declinations and any desired grid options. The only required argument is
``source_declination``. See :ref:`signal-subtraction-grid-options` for the
full list of optional parameters.

.. code-block:: python

   import numpy as np
   from kingmaker.pdf import MarginalizedKingPDF

   source_decs = np.radians(np.linspace(-80, 80, 21))

   mkpdf = MarginalizedKingPDF(
       source_declination=source_decs,
       angular_cutoff=np.radians(10.0),
       # Optional: customize the interpolation grid.
       # points_alpha=np.logspace(np.log10(np.radians(0.05)), np.log10(np.pi), 30),
       # points_beta=np.logspace(0.01, 1, 20),
       # n_signed_delta_dec=200,
       # n_ra_bins=100,
   )

The grid is built once during :meth:`~kingmaker.pdf.MarginalizedKingPDF.__init__`
using a parallelized numba kernel. The resulting array is stored as
``mkpdf._grid`` with shape
``(n_sources, n_alpha, n_beta, n_signed_delta_dec)``.

.. _signal-subtraction-interfaces:

Choosing between ``pdf()`` and ``evaluate()``
---------------------------------------------

:class:`~kingmaker.pdf.MarginalizedKingPDF` exposes two methods for evaluating
the RA-marginalized King PDF:

:meth:`~kingmaker.pdf.MarginalizedKingPDF.pdf` ``(x, alpha, beta, source_dec)``
    Evaluates the marginalized PDF for a **single source** at an array of
    reconstructed event declinations ``x``. Returns a **dense NumPy array** of
    shape ``(n_events,)``. Use this for plotting the marginalized profile,
    quick sanity checks, or any context where you only need one source at a
    time.

:meth:`~kingmaker.pdf.MarginalizedKingPDF.evaluate` ``(source_decs, event_decs, alpha, beta)``
    Evaluates the marginalized PDF for **all (event, source) pairs** in a
    single call and returns a :class:`scipy.sparse.csr_array` of shape
    ``(n_events, n_sources)``. Events more than ``angular_cutoff`` from a
    source contribute zero and are omitted from the sparse structure, so the
    result is compact even for large event and source counts. This is the
    method to use in analyses. On repeat calls with the same event and source
    positions but different ``alpha``/``beta`` values, pass the result of a
    previous call as ``mask=`` to skip the angular-masking loop and go
    straight to the interpolation step.

In short: use ``pdf()`` for single-source plots and demos; use ``evaluate()``
for analysis-level evaluation over many sources and events.

.. _signal-subtraction-plotting:

Full-sky King PDF and its RA-marginalized profile
-------------------------------------------------

The following example compares the 2D King PDF on the sky to the 1D
marginalized profile at a fixed :math:`(\alpha, \beta)`.

.. code-block:: python

   import numpy as np
   import matplotlib.pyplot as plt
   from kingmaker.pdf import MarginalizedKingPDF
   from kingmaker.utils import angular_distance

   alpha, beta = np.radians(2.0), 2.5
   source_dec = np.radians(30.0)
   cutoff = np.radians(10.0)

   mkpdf = MarginalizedKingPDF(
       source_declination=np.array([source_dec]),
       angular_cutoff=cutoff,
   )

   # --- 2D full-sky slice: PDF as a function of (dRA, dec_reco) ---
   dra = np.linspace(-cutoff, cutoff, 300)
   dec_reco_2d = np.linspace(source_dec - cutoff, source_dec + cutoff, 300)
   DRA, DEC = np.meshgrid(dra, dec_reco_2d)
   psi = angular_distance(0.0, source_dec, DRA, DEC)
   pdf_2d = mkpdf.king.pdf_from_norm(psi, alpha, beta, mkpdf.king.norm(alpha, beta))

   # --- 1D marginalized profile ---
   dec_reco_1d = np.linspace(source_dec - cutoff, source_dec + cutoff, 200)
   alpha_arr = np.full(len(dec_reco_1d), alpha)
   beta_arr = np.full(len(dec_reco_1d), beta)
   pdf_marg = mkpdf.pdf(dec_reco_1d, alpha_arr, beta_arr, source_dec)

   fig, axes = plt.subplots(1, 2, figsize=(10, 4))

   ax = axes[0]
   im = ax.pcolormesh(
       np.degrees(dra),
       np.degrees(dec_reco_2d),
       pdf_2d,
       cmap="magma",
   )
   ax.set_xlabel(r"$\Delta\mathrm{RA}$ (deg)")
   ax.set_ylabel(r"$\delta_\mathrm{reco}$ (deg)")
   ax.set_title(r"King PDF $f(\psi\,|\,\alpha,\beta)$")
   plt.colorbar(im, ax=ax, label="PDF (sr⁻¹)")

   ax = axes[1]
   ax.plot(np.degrees(dec_reco_1d), pdf_marg)
   ax.set_xlabel(r"$\delta_\mathrm{reco}$ (deg)")
   ax.set_ylabel(r"$\mathcal{M}(\delta_\mathrm{reco}|\delta_\mathrm{true},\alpha,\beta)$ (rad⁻¹)")
   ax.set_title(r"RA-marginalized profile")

   fig.suptitle(
       rf"$\alpha={np.degrees(alpha):.1f}°$, $\beta={beta}$,"
       rf" $\delta_\mathrm{{true}}={np.degrees(source_dec):.0f}°$"
   )
   plt.tight_layout()
   plt.show()

.. _signal-subtraction-evaluation:

Evaluating the marginalized PDF over arrays of events and sources
-----------------------------------------------------------------

:meth:`~kingmaker.pdf.MarginalizedKingPDF.pdf` accepts per-event arrays of
``dec_reco``, ``alpha``, and ``beta`` together with an array of source
declinations, and returns a :class:`scipy.sparse.csr_array` of shape
``(n_events, n_sources)``.

.. code-block:: python

   import numpy as np
   from kingmaker.pdf import MarginalizedKingPDF

   rng = np.random.default_rng(42)
   n_events = 5000
   n_sources = 200

   source_decs = np.radians(np.linspace(-80, 80, n_sources))

   mkpdf = MarginalizedKingPDF(
       source_declination=source_decs,
       angular_cutoff=np.radians(10.0),
   )

   # Simulated reconstructed event parameters.
   dec_reco = np.arcsin(rng.uniform(-1, 1, n_events))
   alpha_evt = np.radians(rng.uniform(0.3, 3.0, n_events))
   beta_evt = rng.uniform(1.5, 5.0, n_events)

   pdf_sparse = mkpdf.evaluate(source_decs, dec_reco, alpha_evt, beta_evt)

   print(f"Output shape  : {pdf_sparse.shape}")
   print(f"Non-zero entries: {pdf_sparse.nnz}")
   print(f"Sparsity      : {1 - pdf_sparse.nnz / (n_events * n_sources):.3%}")

   # Extract the column for one source as a dense array:
   source_index = 100
   pdf_for_source = pdf_sparse[:, source_index].toarray().ravel()
   print(f"Max marginalized PDF for source {source_index}: {pdf_for_source.max():.4f} rad⁻¹")

   # Convert the full matrix to dense only if memory permits:
   # pdf_dense = pdf_sparse.toarray()  # shape (n_events, n_sources)

.. _signal-subtraction-grid-options:

Grid layout and configuration options
--------------------------------------

The marginalization cache is a 4D array indexed by:

#. **source_declination** — the sorted source declination(s) passed to the
   constructor.
#. **log₁₀ α** — log-spaced from ``points_alpha`` (default: 30 values,
   0.05° to π).
#. **β** — log-spaced from ``points_beta`` (default: 20 values,
   ``np.logspace(0.01, 1, 20)`` ≈ 1.023 to 10).
#. **signed_delta_dec = δ_reco − δ_source** — uniformly spaced over
   [−``angular_cutoff``, +``angular_cutoff``] (default: 200 points).

The RA integral at each grid point is computed via trapezoidal quadrature
over ``n_ra_bins + 1`` nodes on [0, π] (default: 100 intervals). All of
these can be customized at construction:

.. code-block:: python

   mkpdf = MarginalizedKingPDF(
       source_declination=np.radians([0.0, 30.0, 60.0]),
       angular_cutoff=np.radians(5.0),
       points_alpha=np.logspace(np.log10(np.radians(0.1)), np.log10(np.pi), 40),
       points_beta=np.logspace(0.01, 1, 30),
       n_signed_delta_dec=300,
       n_ra_bins=200,
   )

Finer grids improve interpolation accuracy at the cost of memory and
build time. For most use cases with realistic PSF parameters (α ≳ 0.1°,
β ∈ [1.5, 8]), the defaults are sufficient.