Fixes for signal subtraction

README.md

@@ -35,7 +35,7 @@ This package implements the King/Moffat distribution on a sphere. The current im
 - **Vectorized operations**: NumPy broadcasting for simultaneous evaluation at multiple points
 - **JIT-compiled kernels**: All inner-loop functions compiled and cached with numba for near-native performance
 - **Multi-dimensional fitting**: Parameterize PSF as function of energy, declination, angular error, etc.
-- **Signal subtraction**: Built-in RA marginalization for likelihood-based analyses
+- **Signal subtraction**: RA-averaged King PDF for signal-subtraction likelihoods
 - **Template smearing**: Incorporate diffuse backgrounds and Galactic plane effects via spherical harmonics
 - **Flexible binning**: Support for both equal-probability and explicit bin edges
 
@@ -134,15 +134,15 @@ beta_value = beta_interp(test_point)
 
 ### Signal-Subtracted Likelihood
 
-`MarginalizedKingPDF` pre-computes the King PDF integrated over right ascension
-on a grid and exposes two methods for different use cases:
+`MarginalizedKingPDF` pre-computes the King PDF averaged over right ascension
+(sr⁻¹) on a grid and exposes two methods for different use cases:
 
-- **`pdf(x, alpha, beta, source_dec)`** — evaluates the marginalized PDF for a
+- **`pdf(x, alpha, beta, source_dec)`** — evaluates the RA-averaged PDF for a
   *single* source at a given set of event declinations. Returns a dense NumPy
   array. Intended for plotting, diagnostics, and demonstrations.
 
 - **`evaluate(source_decs, event_decs, alpha, beta)`** — evaluates the
-  marginalized PDF for *all* (event, source) pairs at once and returns a
+  RA-averaged PDF for *all* (event, source) pairs at once and returns a
   `scipy.sparse.csr_array` of shape `(n_events, n_sources)`. The sparse format
   is efficient because events far from a source (beyond `angular_cutoff`)
   contribute zero and are omitted. This is the method to use in analyses.
@@ -183,7 +183,7 @@ print(f"Nonzero fraction: {pdf_matrix.nnz / pdf_matrix.shape[0] / pdf_matrix.sha
 pdf_matrix_2 = mkpdf.evaluate(source_decs, ev_decs, alpha_evt, beta_evt, mask=pdf_matrix)
 ```
 
-`KingSpatialLikelihood` integrates the marginalized path directly. Pass
+`KingSpatialLikelihood` supports the RA-averaged path directly. Pass
 `enable_marginalization=True` and a source-declination array at construction;
 `set_events` then precomputes one sparse matrix per spectral index so that
 `evaluate_marginalized_pdf` can interpolate in O(nnz) without additional grid
@@ -214,7 +214,7 @@ wrapper.set_events(data_events, source_ras=np.array([src_ra]),
 pdf_matrix = wrapper.evaluate_pdf(data_events, gamma=2.0)
 pdf_values = pdf_matrix.toarray().ravel()
 
-# RA-marginalized PDF: sparse (n_events, n_sources) matrix.
+# RA-averaged PDF (sr⁻¹): sparse (n_events, n_sources) matrix.
 pdf_sparse = wrapper.evaluate_marginalized_pdf(data_events, gamma=2.0)
 ```
 

docs/examples.rst

@@ -9,9 +9,9 @@ Jupyter notebook in ``examples/`` for additional plots and benchmarks.
 King PDF basics
 ----------------
 
-:class:`~kingmaker.pdf.KingPDF` evaluates the PDF/CDF (see :doc:`quickstart`),
-draws samples, and marginalizes over right ascension for signal-subtraction
-likelihoods.
+:class:`~kingmaker.pdf.KingPDF` evaluates the PDF/CDF (see :doc:`quickstart`)
+and draws samples; :class:`~kingmaker.pdf.MarginalizedKingPDF` averages it over
+right ascension for signal-subtraction likelihoods.
 
 **Sampling angular offsets**
 
@@ -29,15 +29,15 @@ likelihoods.
    # (ra, dec) for a given true source position.
    psi = king.sample(10_000, alpha, beta, n_grid=10_000)
 
-**Marginalizing over right ascension**
+**Averaging over right ascension**
 
-Signal-subtraction likelihood terms need the King PDF marginalized over
-right ascension: the probability of an event reconstructing at a given
-declination, for a source at a given true declination.
+Signal-subtraction likelihood terms need the King PDF averaged over right
+ascension: the signal density (sr⁻¹) at a reconstructed declination after RA
+scrambling, for a source at a given true declination.
 :class:`~kingmaker.pdf.MarginalizedKingPDF` exposes two methods for this:
 
 - :meth:`~kingmaker.pdf.MarginalizedKingPDF.pdf` ``(x, alpha, beta, source_dec)``
-  evaluates the marginalized profile for a **single source** and returns a
+  evaluates the RA-averaged profile for a **single source** and returns a
   **dense array**. Intended for plotting and single-source checks.
 - :meth:`~kingmaker.pdf.MarginalizedKingPDF.evaluate` ``(source_decs, event_decs, alpha, beta)``
   evaluates all (event, source) pairs at once and returns a
@@ -75,7 +75,7 @@ declination, for a source at a given true declination.
 - ``source_declination`` / ``points_alpha`` / ``points_beta`` /
   ``n_signed_delta_dec`` / ``n_ra_bins`` (:class:`~kingmaker.pdf.MarginalizedKingPDF`
   constructor): control the source positions and interpolation grid for the
-  RA-marginalized PDF. See :doc:`signal_subtraction` for details.
+  RA-averaged PDF. See :doc:`signal_subtraction` for details.
 
 `basic_demo.ipynb <https://github.com/mjlarson/kingmaker/blob/main/examples/basic_demo.ipynb>`_
     Parameter effects, normalization checks, and sampling/evaluation speed
@@ -180,7 +180,7 @@ End-to-end point-source likelihood
 :class:`~kingmaker.fitting.KingPSFFitter` and
 :class:`~kingmaker.pdf.KingPDF` behind a single interface: fit (or load
 cached fit results) once, then evaluate the PDF per-event many times across
-trials. It also optionally integrates
+trials. It also optionally wraps
 :class:`~kingmaker.pdf.MarginalizedKingPDF` for signal-subtraction
 likelihoods via :meth:`~kingmaker.wrapper.KingSpatialLikelihood.evaluate_marginalized_pdf`.
 
@@ -201,14 +201,14 @@ above:
        parametrization_bins=parametrization_bins,
        spectral_indices=[1.0, 2.0, 3.0, 4.0],
        cache_parameters=False,
-       # Enable the RA-marginalized path for signal-subtraction likelihoods.
+       # Enable the RA-averaged path for signal-subtraction likelihoods.
        enable_marginalization=True,
        marginalization_source_decs=np.array([source_dec]),
        marginalization_angular_cutoff=np.radians(10.0),
    )
 
    # Per trial: cache per-event parameters once, then evaluate as needed.
-   # set_events precomputes both the standard and marginalized PDF matrices.
+   # set_events precomputes both the standard and RA-averaged PDF matrices.
    wrapper.set_events(
        data_events,
        source_ras=np.array([source_ra]),
@@ -221,7 +221,7 @@ above:
    pdf_matrix_steeper = wrapper.evaluate_pdf(data_events, gamma=2.5)  # interpolated
    pdf_values = pdf_matrix.toarray().ravel()  # dense 1-D for downstream use
 
-   # RA-marginalized PDF — sparse (n_events, n_sources) matrix.
+   # RA-averaged PDF (sr⁻¹) — sparse (n_events, n_sources) matrix.
    # Cheap O(nnz) interpolation on precomputed .data arrays; no interpn call.
    marg_matrix = wrapper.evaluate_marginalized_pdf(data_events, gamma=2.0)
    marg_steeper = wrapper.evaluate_marginalized_pdf(data_events, gamma=2.5)
@@ -240,12 +240,12 @@ above:
   :class:`~kingmaker.fitting.KingPSFFitter`.
 - ``enable_marginalization`` / ``marginalization_source_decs``: set
   ``enable_marginalization=True`` and supply a source-declination array to
-  activate the marginalized path. At construction,
+  activate the RA-averaged path. At construction,
   :class:`~kingmaker.pdf.MarginalizedKingPDF` builds its 4D interpolation
   grid (expensive, done once). ``set_events`` then precomputes one sparse
   matrix per spectral index; ``evaluate_marginalized_pdf`` interpolates
   between them in O(nnz) without further grid lookups.
-- ``marginalization_angular_cutoff``: cutoff for the marginalized PDF,
+- ``marginalization_angular_cutoff``: cutoff for the RA-averaged PDF,
   independent of the point-source ``angular_cutoff``. Defaults to the same
   value as ``angular_cutoff`` if not set.
 - **Gotcha:** ``evaluate_pdf`` and ``evaluate_marginalized_pdf`` both require

docs/signal_subtraction.rst

@@ -1,7 +1,7 @@
 .. _signal-subtraction:
 
-Signal-subtraction likelihoods and the marginalized King PDF
-=============================================================
+Signal-subtraction likelihoods and the RA-averaged 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
@@ -14,34 +14,36 @@ The result is called the *signal-subtraction* (or *background-marginalized*) lik
 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}`:
+the King PDF averaged over right ascension — the signal density (per steradian)
+at reconstructed declination :math:`\delta_\mathrm{reco}` after scrambling in
+right ascension, 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},
+   = \frac{1}{2\pi} \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
+position and the source, and :math:`f` is the King PDF. Like :math:`f`,
+:math:`\mathcal{M}` integrates to one over the sphere, so the density in
+reconstructed declination alone is :math:`2\pi \cos\delta_\mathrm{reco}\,\mathcal{M}`.
+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 average 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},
+   = \frac{1}{\pi} \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
+The :class:`~kingmaker.pdf.MarginalizedKingPDF` class evaluates this average 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
@@ -71,7 +73,7 @@ full list of optional parameters.
        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),
+       # points_beta=np.geomspace(1.01, 1000, 60),
        # n_signed_delta_dec=200,
        # n_ra_bins=100,
    )
@@ -87,17 +89,17 @@ Choosing between ``pdf()`` and ``evaluate()``
 ---------------------------------------------
 
 :class:`~kingmaker.pdf.MarginalizedKingPDF` exposes two methods for evaluating
-the RA-marginalized King PDF:
+the RA-averaged King PDF:
 
 :meth:`~kingmaker.pdf.MarginalizedKingPDF.pdf` ``(x, alpha, beta, source_dec)``
-    Evaluates the marginalized PDF for a **single source** at an array of
+    Evaluates the RA-averaged 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,
+    shape ``(n_events,)``. Use this for plotting the RA-averaged 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
+    Evaluates the RA-averaged 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
@@ -112,11 +114,11 @@ for analysis-level evaluation over many sources and events.
 
 .. _signal-subtraction-plotting:
 
-Full-sky King PDF and its RA-marginalized profile
--------------------------------------------------
+Full-sky King PDF and its RA-averaged profile
+---------------------------------------------
 
 The following example compares the 2D King PDF on the sky to the 1D
-marginalized profile at a fixed :math:`(\alpha, \beta)`.
+RA-averaged profile at a fixed :math:`(\alpha, \beta)`.
 
 .. code-block:: python
 
@@ -141,7 +143,7 @@ marginalized profile at a fixed :math:`(\alpha, \beta)`.
    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 ---
+   # --- 1D RA-averaged 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)
@@ -164,8 +166,8 @@ marginalized profile at a fixed :math:`(\alpha, \beta)`.
    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")
+   ax.set_ylabel(r"$\mathcal{M}(\delta_\mathrm{reco}|\delta_\mathrm{true},\alpha,\beta)$ (sr⁻¹)")
+   ax.set_title(r"RA-averaged profile")
 
    fig.suptitle(
        rf"$\alpha={np.degrees(alpha):.1f}°$, $\beta={beta}$,"
@@ -176,10 +178,10 @@ marginalized profile at a fixed :math:`(\alpha, \beta)`.
 
 .. _signal-subtraction-evaluation:
 
-Evaluating the marginalized PDF over arrays of events and sources
------------------------------------------------------------------
+Evaluating the RA-averaged PDF over arrays of events and sources
+----------------------------------------------------------------
 
-:meth:`~kingmaker.pdf.MarginalizedKingPDF.pdf` accepts per-event arrays of
+:meth:`~kingmaker.pdf.MarginalizedKingPDF.evaluate` 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)``.
@@ -214,7 +216,7 @@ declinations, and returns a :class:`scipy.sparse.csr_array` of shape
    # 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⁻¹")
+   print(f"Max RA-averaged PDF for source {source_index}: {pdf_for_source.max():.4f} sr⁻¹")
 
    # Convert the full matrix to dense only if memory permits:
    # pdf_dense = pdf_sparse.toarray()  # shape (n_events, n_sources)
@@ -224,32 +226,39 @@ declinations, and returns a :class:`scipy.sparse.csr_array` of shape
 Grid layout and configuration options
 --------------------------------------
 
-The marginalization cache is a 4D array indexed by:
+The interpolation grid 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:
+#. **β** — log-spaced from ``points_beta`` (default: 60 values,
+   ``np.geomspace(1.01, 1000, 60)``, matching the
+   :class:`~kingmaker.fitting.KingPSFFitter` bounds).
+#. **signed_delta_dec = δ_reco − δ_source** — asinh-spaced over
+   [−``angular_cutoff``, +``angular_cutoff``] (default: 200 points): uniform
+   within half the smallest ``points_alpha`` of zero, logarithmic beyond.
+
+The RA average at each grid point is computed via trapezoidal quadrature
+over ``n_ra_bins + 1`` asinh-spaced nodes spanning the RA half-width of
+``angular_cutoff`` around each source (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),
+       points_beta=np.geomspace(1.01, 1000, 90),
        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.
+build time. For most use cases with realistic PSF parameters (α ≳ 0.1°),
+the defaults are sufficient.
+
+Queries with α or β outside the grid are clamped to the nearest grid edge
+and emit a :class:`RuntimeWarning`; widen ``points_alpha`` or
+``points_beta`` if this appears.