Fit the CDF with scipy.optimize.minimize

.github/workflows/ci.yml

@@ -51,7 +51,6 @@ jobs:
         test_file:
           - test_basic
           - test_king_pdf
-          - test_interpolated_king_pdf
           - test_template_smeared_king_pdf
           - test_utils
           - test_fitting

README.md

@@ -8,36 +8,35 @@ A Python library for working with King/Moffat distributions for modeling point s
 
 ## Overview
 
-The **King distribution** (also known as the **Moffat distribution**, see [wikipedia](https://en.wikipedia.org/wiki/Moffat_distribution).) is a two-parameter probability distribution commonly used to describe the point spread function (PSF) of astronomical observations. The distribution has the form:
+The **King distribution** (also known as the **Moffat distribution**, see [wikipedia](https://en.wikipedia.org/wiki/Moffat_distribution).) is a two-parameter probability distribution commonly used to describe the point spread function (PSF) of astronomical observations. This package uses the exact spherical form:
 
 ```
-f(θ | α, β) ∝ [1 + (θ/α)² / (2β)]^(-β)
+f(theta | alpha, beta) ∝ [1 + (1 - cos(theta)) / (alpha^2 * beta)]^(-beta)
 ```
 
 where:
-- **θ** is the angular separation from the source
-- **α** (alpha) is the scale parameter controlling the width of the distribution
-- **β** (beta) is the shape parameter controlling the tail behavior
+- **theta** is the angular separation from the source
+- **alpha** is the scale parameter controlling the width of the distribution
+- **beta** is the shape parameter controlling the tail behavior (must be > 1)
 
-The distribution's flexibility in modeling both the core and tail regions makes it particularly well-suited for modeling extended tails to point spread functions. These are common in many instruments in astronomy. The King distribution offers significant advantages over simpler models (e.g., Gaussian or Rayleigh) by providing independent control over the core width (α) and tail weight (β), allowing more accurate modeling of real detector responses. 
+The substitution `t = 1 - cos(theta)`, `dt = sin(theta) dtheta` yields an exact closed-form CDF, making both normalization and cumulative probabilities analytically tractable at all angular scales. For small angles, `1 - cos(theta) ~ theta^2/2`, recovering the familiar flat-sky form. In the limit beta -> inf the distribution converges to a Gaussian.
 
-Note that in the limit of β → ∞, this distribution reproduces the standard Rayleigh distribution. 
+The distribution's flexibility in modeling both the core and tail regions makes it particularly well-suited for modeling extended tails to point spread functions. These are common in many instruments in astronomy. The King distribution offers significant advantages over simpler models (e.g., Gaussian or Rayleigh) by providing independent control over the core width (alpha) and tail weight (beta), allowing more accurate modeling of real detector responses.
 
-This package implements the King/Moffat distribution on a sphere. Normalizations are numerically calculated at initialization and linearly interpolated at runtime. The current implementation includes support for
+This package implements the King/Moffat distribution on a sphere. The current implementation includes support for
 
-- Direct KingPDF evaluations with normalizations calculated at runtime
-- KingPDF evaluations with normalizations interpolated from tables calculated at intialization
+- KingPDF evaluations with exact analytic normalizations and CDFs
 - Numerically calculated "signal subtraction" PDFs
-- Preliminary healpix map convolutions and template support using spherical convolutions
+- HEALPix map convolutions and template support using spherical harmonics
 
 ## Features
 
-- **Efficient normalization**: Accurate integration over spherical caps with arbitrary cutoffs
+- **Exact normalization**: Closed-form analytic normalization and CDF via spherical geometry, valid at all angular scales
 - **Vectorized operations**: NumPy broadcasting for simultaneous evaluation at multiple points
-- **Interpolation caching**: Pre-computed normalizations for 10-200× speedup
+- **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
-- **Template smearing**: Incorporate diffuse backgrounds and Galactic plane effects
+- **Template smearing**: Incorporate diffuse backgrounds and Galactic plane effects via spherical harmonics
 - **Flexible binning**: Support for both equal-probability and explicit bin edges
 
 ## Installation
@@ -88,22 +87,6 @@ containment_68 = brentq(
 print(f"68% containment: {np.degrees(containment_68):.2f} degrees")
 ```
 
-### Interpolated PDF for Speed
-
-For applications requiring many PDF evaluations, use `InterpolatedKingPDF` which pre-computes normalizations on a grid:
-
-```python
-from kingmaker.pdf import InterpolatedKingPDF
-
-# Create interpolated version (10-200x faster)
-king_interp = InterpolatedKingPDF(angular_cutoff=np.pi)
-
-# Use identical interface
-pdf_values = king_interp.pdf(angles, alpha, beta)
-```
-
-The interpolated version provides **10-200× speedup** with **< 0.1% error** for most parameter values.
-
 ### Fitting PSF Parameters to Monte Carlo
 
 Fit King distribution parameters to simulated signal events as a function of an arbitrary set of observables:
@@ -165,24 +148,23 @@ sindec_bins, pdf_marginalized = king.marginalize(source_dec, alpha, beta)
 
 ### Template Smearing
 
-Apply Galactic or other template smearing to the PSF:
+Apply a HEALPix template (e.g. Galactic diffuse emission) to the PSF via spherical harmonic convolution:
 
 ```python
 from kingmaker.pdf import TemplateSmearedKingPDF
+import healpy as hp
 
-# Load a diffuse template (e.g., Galactic plane emission)
-template = np.load('galactic_template.npy')  # HEALPix map
+# Load a diffuse template as a HEALPix map
+skymap = hp.read_map('galactic_template.fits')
 
-# Create smeared PSF
-king_smeared = TemplateSmearedKingPDF(
-    angular_cutoff=np.pi,
-    template=template,
-    nside=128
-)
+# Create smeared PSF (pre-computes spherical harmonic expansion at init)
+king_smeared = TemplateSmearedKingPDF(skymap, angular_cutoff=np.pi)
 
-# Evaluate PDF including template effects
-ra, dec = np.radians(120), np.radians(30)  # Source position
-pdf_smeared = king_smeared.pdf(angles, ra, dec, alpha, beta)
+# Set evaluation coordinates, then convolve
+eval_ras  = np.radians([0.0, 45.0, 90.0])
+eval_decs = np.radians([0.0, 30.0, -15.0])
+king_smeared.set_coordinates(eval_decs, eval_ras)
+pdf_vals = king_smeared.convolve_at_grid_point(alpha, beta)
 ```
 
 

docs/api.rst

@@ -9,18 +9,8 @@ PDF Classes
    :nosignatures:
 
    kingmaker.pdf.KingPDF
-   kingmaker.pdf.InterpolatedKingPDF
    kingmaker.pdf.TemplateSmearedKingPDF
 
-.. autoclass:: kingmaker.pdf.KingPDF
-   :members:
-
-.. autoclass:: kingmaker.pdf.InterpolatedKingPDF
-   :members:
-
-.. autoclass:: kingmaker.pdf.TemplateSmearedKingPDF
-   :members:
-
 Distribution
 ------------
 

docs/conf.py

@@ -41,4 +41,3 @@
 exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
 
 html_theme = "furo"
-html_static_path = ["_static"]

docs/examples.rst

@@ -4,8 +4,8 @@ Examples
 The ``examples/`` directory contains Jupyter notebooks demonstrating typical workflows.
 
 `basic_demo.ipynb <https://github.com/mjlarson/kingmaker/blob/main/examples/basic_demo.ipynb>`_
-    King PDF basics, parameter effects, normalization, sampling, and speed comparisons
-    between :class:`~kingmaker.pdf.KingPDF` and :class:`~kingmaker.pdf.InterpolatedKingPDF`.
+    King PDF basics, parameter effects, normalization, sampling, and speed benchmarks
+    for :class:`~kingmaker.pdf.KingPDF`.
 
 `fitting_demo.ipynb <https://github.com/mjlarson/kingmaker/blob/main/examples/fitting_demo.ipynb>`_
     Fitting King PSF parameters to Monte Carlo simulations using

docs/index.rst

@@ -9,11 +9,15 @@ PSFs with independent control over the core width (α) and tail weight (β):
 
 .. math::
 
-   f(\theta \mid \alpha, \beta) \propto \left[1 + \frac{(\theta/\alpha)^2}{2\beta}\right]^{-\beta}
+   f(\theta \mid \alpha, \beta) \propto \left[1 + \frac{1 - \cos\theta}{\alpha^2 \beta}\right]^{-\beta}
 
-This package implements the King distribution on a sphere, including
-normalization over spherical caps, parameter fitting from simulation,
-and spherical-harmonic convolution with HEALPix template maps.
+The substitution :math:`t = 1 - \cos\theta` yields an exact closed-form CDF
+via a standard power-law integral, valid at all angular scales. For small
+angles, :math:`1 - \cos\theta \approx \theta^2/2`, recovering the flat-sky form.
+
+This package implements the King distribution on a sphere, including exact
+analytic normalization, parameter fitting from simulation, and
+spherical-harmonic convolution with HEALPix template maps.
 
 .. toctree::
    :maxdepth: 2

kingmaker/__init__.py

@@ -1,149 +1,6 @@
 # Import main classes for convenience
-from .pdf import KingPDF, InterpolatedKingPDF, TemplateSmearedKingPDF
+from .pdf import KingPDF, TemplateSmearedKingPDF
 from .fitting import KingPSFFitter
+from .wrapper import KingSpatialLikelihood
 
-'''class KingSpatialPDF:
-    def __init__(self):
-        pass
-
-    def pdf(self, x, alpha, beta):
-        return self._norm(alpha, beta) * self._unnormalized_pdf(x, alpha, beta)
-
-    def cdf(self, x, alpha, beta):
-        return self._norm(alpha, beta) * self._unnormalized_cdf(x, alpha, beta)
-
-    def _to_1d(self, *args):
-        return (np.atleast_1d(_) for _ in args)
-
-    def _check_args(self, **kwargs):
-        # Ensure all args are the same shape
-        if len(kwargs.keys()) > 1:
-            shape = list(kwargs.values())[0].shape
-            for key, val in kwargs.items():
-                if len(val.shape) > 1:
-                    message = (
-                        f"The variable {key} has more than one dimension (shape = {val.shape})."
-                        " This is not supported at the moment. Try calling again with 1d values."
-                    )
-                    raise NotImplementedError(message)
-                if val.shape != shape:
-                    message = (
-                        f"The variables {list(kwargs.keys())} have different shapes "
-                        f"{list(_.shape for _ in kwargs.values())}."
-                    )
-                    raise ValueError(message)
-        return
-
-    def _unnormalized_pdf(self, x, alpha, beta):
-        """
-        Evaluate the unnormalized radial King function (without solid angle Jacobian).
-
-        This function returns the radial shape:
-            f(x) = [1 + (x/alpha)^2 / (2*beta)]^(-beta)
-
-        Parameters
-        ----------
-        x : float or ndarray
-            Angular separation from the source, in radians with shape M.
-        alpha : float or ndarray
-            King distribution alpha parameter (scale) with shape N.
-        beta : float or ndarray
-            King distribution beta parameter (tail weight) with shape N.
-
-        Returns
-        -------
-        ndarray
-            Unnormalized King function values with units of probability/sterradian.
-            The output shape will be either N (if `x` is a scalar) or (N, M) (if `x`
-            is an ndarray).
-        """
-        print("_unnormalized_pdf")
-
-        # Check if x is a float or ndarray
-        isscalar = np.isscalar(x)
-
-        # Ensure everything is 1d
-        x, alpha, beta = self._to_1d(x, alpha, beta)
-        self._check_args(alpha=alpha, beta=beta)
-
-        # Broadcast
-        alpha = alpha[:, None]
-        beta = beta[:, None]
-        x = x[None, :]
-
-        # Calculate the King function value and return
-        return (1 + (x / alpha) ** 2 / (2 * beta)) ** -beta
-
-    def _unnormalized_cdf(self, x, alpha, beta):
-        """
-        Evaluate the CDF of the radial King function (without solid angle Jacobian).
-
-        The integral includes the sin(theta) for spherical coordinates and
-        normalizes over solid angle (i.e., integral(PDF * sin(theta) dtheta dphi = 1)).
-
-        Parameters
-        ----------
-        x : scalar or ndarray
-            The locations at which to calculate the CDF. If this is a scalar, use log-binning
-            with 1000 bins from 1e-5 to this value, assuming radians. If this is an ndarray,
-            the shape can be M, indepedent of the shape of alpha and beta.
-        alpha : float or ndarray
-            King distribution alpha parameter (scale) with shape N.
-        beta : float or ndarray
-            King distribution beta parameter (tail weight) with shape N.
-
-        Returns
-        -------
-        float or ndarray
-            King function CDF values with units of probability. If x is a single value,
-            only return the last value. Otherwise, the shape will be (N, M).
-        """
-        print("_unnormalized_cdf")
-        # Check if x is a float or ndarray
-        isscalar = np.isscalar(x)
-        if isscalar:
-            x = np.logspace(np.log10(np.pi) - 5, np.log10(x), 1000)
-
-        # Broadcast
-        x, alpha, beta = self._to_1d(x, alpha, beta)
-        alpha = alpha[:, None]
-        beta = beta[:, None]
-        x = x[None, :]
-
-        # Use the unnormalized PDF (which doesn't include sin(x))
-        unnormalized = np.array([self._unnormalized_pdf(x, a, b) for a, b in zip(alpha, beta)])
-        integrand = unnormalized * np.sin(x)
-
-        print(x.shape, integrand.shape)
-        cdf = cumulative_trapezoid(integrand, x, axis=-1, initial=0)
-
-        # If we only want a scalar value out, grab only what we need
-        if isscalar:
-            return cdf[..., -1]
-
-    def _norm(self, alpha, beta):
-        """
-        Compute the normalization constant for the King PDF over the sphere.
-
-        The integral includes the sin(theta) for spherical coordinates and
-        normalizes over solid angle (i.e., integral(PDF * sin(theta) dtheta dphi = 1)).
-
-        Parameters
-        ----------
-        alpha : float or ndarray
-            King distribution alpha parameter (scale) with shape N.
-        beta : float or ndarray
-            King distribution beta parameter (tail weight) with shape N.
-
-        Returns
-        -------
-        ndarray
-            Normalization constants such that PDF integrates to 1 over the sphere
-            with shape N.
-        """
-        print("norm")
-        cdf = self._unnormalized_cdf(np.pi, alpha, beta)
-        return cdf
-'''
-
-__all__ = ["KingPDF", "InterpolatedKingPDF", "TemplateSmearedKingPDF", "KingPSFFitter"]
+__all__ = ["KingPDF", "TemplateSmearedKingPDF", "KingPSFFitter", "KingSpatialLikelihood"]