Pipeline

The orchestration layer: coerce inputs, select stamps, pick the convolution direction and basis scale, run the kernel solve and subtraction, and optionally decorrelate and score. See the usage guide for the per-knob tuning notes.

subtract

subtract(science, reference, **kwargs) -> DiffResult

Convenience wrapper: configure a :class:Subtractor and run it once.

Configuration keywords (n_max, beta, n_knots, stamp_radius, threshold_sigma, max_stamps, saturation, bright_mask, lambda_grid, clip_sigma, clip_iterations, min_stamps, cv_folds, spatial_scale, decorrelate, score, block, rescale_variance) are split from the per-call inputs (science_var, reference_var, science_mask, reference_mask, gain, read_noise, catalog, direction).

Source code in python/delta/pipeline.py

def subtract(science, reference, **kwargs) -> DiffResult:
    """Convenience wrapper: configure a :class:`Subtractor` and run it once.

    Configuration keywords (``n_max``, ``beta``, ``n_knots``, ``stamp_radius``,
    ``threshold_sigma``, ``max_stamps``, ``saturation``, ``bright_mask``,
    ``lambda_grid``, ``clip_sigma``, ``clip_iterations``, ``min_stamps``, ``cv_folds``,
    ``spatial_scale``, ``decorrelate``, ``score``, ``block``, ``rescale_variance``)
    are split from the per-call inputs
    (``science_var``, ``reference_var``, ``science_mask``, ``reference_mask``,
    ``gain``, ``read_noise``, ``catalog``, ``direction``).
    """
    config = {k: v for k, v in kwargs.items() if k in _CONFIG_KEYS}
    call = {k: v for k, v in kwargs.items() if k not in _CONFIG_KEYS}
    return Subtractor(**config).subtract(science, reference, **call)

apply

apply(
    solution: KernelSolution, science, reference, **kwargs
) -> DiffResult

Convenience wrapper: apply a frozen solution without re-fitting.

Configuration keywords (a subset of those accepted by :func:subtract: saturation, bright_mask, stamp_radius, decorrelate, score, block — the fit-only knobs are ignored) are split from the per-call inputs (science_var, reference_var, science_mask, reference_mask, gain, read_noise). The convolution direction, kernel basis, knots and coefficients all come from solution; the inputs must match solution.shape. See :meth:Subtractor.apply.

Source code in python/delta/pipeline.py

def apply(solution: KernelSolution, science, reference, **kwargs) -> DiffResult:
    """Convenience wrapper: apply a frozen ``solution`` without re-fitting.

    Configuration keywords (a subset of those accepted by :func:`subtract`:
    ``saturation``, ``bright_mask``, ``stamp_radius``, ``decorrelate``, ``score``,
    ``block`` — the fit-only knobs are ignored) are split from the per-call inputs
    (``science_var``, ``reference_var``, ``science_mask``, ``reference_mask``,
    ``gain``, ``read_noise``). The convolution direction, kernel basis, knots and
    coefficients all come from ``solution``; the inputs must match
    ``solution.shape``. See :meth:`Subtractor.apply`.
    """
    config = {k: v for k, v in kwargs.items() if k in _CONFIG_KEYS}
    call = {k: v for k, v in kwargs.items() if k not in _CONFIG_KEYS}
    return Subtractor(**config).apply(solution, science, reference, **call)

Subtractor

Subtractor(
    *,
    n_max: int = 6,
    beta: float | None = None,
    n_knots: int = 5,
    radius: int = 0,
    stamp_radius: int = 15,
    threshold_sigma: float = 5.0,
    max_stamps: int = 200,
    saturation: float = 0.0,
    match_radius: int = 2,
    fwhm_radius: int = 0,
    bright_mask: float = 0.0,
    lambda_grid: NDArray[float64] | None = None,
    clip_sigma: float = 4.0,
    clip_iterations: int = 5,
    min_stamps: int = 5,
    cv_folds: int = 5,
    spatial_scale: float | None = None,
    decorrelate: bool = False,
    score: bool = False,
    source_catalog: bool = True,
    block: int = 256,
    rescale_variance: bool = False,
)

Reusable difference-imaging configuration (SPEC §7 object-oriented path).

Source code in python/delta/pipeline.py

def __init__(
    self,
    *,
    n_max: int = 6,  # 6 is a robust default; 8 sharpens bright-star cores on
    beta: float | None = None,  # well-sampled frames (overfits sparse ones).
    n_knots: int = 5,
    radius: int = 0,
    stamp_radius: int = 15,
    threshold_sigma: float = 5.0,
    max_stamps: int = 200,
    saturation: float = 0.0,
    match_radius: int = 2,
    fwhm_radius: int = 0,
    bright_mask: float = 0.0,
    lambda_grid: NDArray[np.float64] | None = None,
    clip_sigma: float = 4.0,
    clip_iterations: int = 5,
    min_stamps: int = 5,
    cv_folds: int = 5,
    spatial_scale: float | None = None,
    decorrelate: bool = False,
    score: bool = False,
    source_catalog: bool = True,
    block: int = 256,
    rescale_variance: bool = False,
) -> None:
    self.n_max = n_max
    self.beta = beta
    self.n_knots = n_knots
    # Kernel footprint half-width; 0 auto-sizes from beta/n_max. Larger values
    # capture more of the matching-kernel wings (comparable to HOTPANTS' -r,
    # its convolution-kernel half-width).
    self.radius = radius
    self.stamp_radius = stamp_radius
    self.threshold_sigma = threshold_sigma
    self.max_stamps = max_stamps
    # Science/reference stamp cross-match tolerance in pixels; widen this if
    # inputs are astrometrically misregistered, which otherwise fails stamp
    # selection with "no stamps selected" and no obvious knob to reach for.
    self.match_radius = match_radius
    # FWHM moment window; 0 auto-sizes to max(5, stamp_radius/2).
    self.fwhm_radius = fwhm_radius
    # `saturation` is the instrument detector level: genuinely saturated
    # (non-linear) pixels, rejected from stamp selection AND masked in the
    # output. `bright_mask` is a data-driven, output-only level: it masks the
    # cores of bright (but linear) stars whose few-percent PSF-match residual
    # would otherwise dominate the difference, without removing those high-SNR
    # stars from the kernel fit. Both grow by the kernel radius in the output.
    self.saturation = saturation
    self.bright_mask = bright_mask
    self.lambda_grid = (
        np.logspace(-6.0, 6.0, 25, dtype=np.float64)
        if lambda_grid is None
        else np.asarray(lambda_grid, dtype=np.float64)
    )
    self.clip_sigma = clip_sigma
    self.clip_iterations = clip_iterations
    self.min_stamps = min_stamps
    # cv_folds >= 2 selects the spatial smoothing by leave-stamp-out k-fold CV
    # (more robust than GCV, which under-smooths correlated stamp pixels).
    self.cv_folds = cv_folds
    # Physical prior: the kernel varies on scales >> a stamp, so knots are
    # placed no finer than `spatial_scale` pixels apart (a coarse spatial
    # basis can't over-fit small-scale structure). None keeps `n_knots`.
    self.spatial_scale = spatial_scale
    self.decorrelate = decorrelate
    self.score = score
    # Source catalog from the score image (SPEC §3.7), built by default
    # whenever `score` is also on; set False to skip (e.g. for callers
    # who only want the score image itself, or care about latency).
    self.source_catalog = source_catalog
    self.block = block
    # Forces the diff-image reduced chi2 to 1 by rescaling `variance` with a
    # single robust scalar (SPEC §3.4) -- a post-hoc correction for noise
    # models (gain/read-noise, or propagated science/reference variance)
    # that under- or over-estimate the true pixel noise.
    self.rescale_variance = rescale_variance

config_cards

config_cards() -> dict[str, object]

Reproducibility cards for knobs not already in KernelSolution.header_cards().

Source code in python/delta/pipeline.py

def config_cards(self) -> dict[str, object]:
    """Reproducibility cards for knobs not already in KernelSolution.header_cards()."""
    return {
        "DLTSRAD": self.stamp_radius,
        "DLTTHSIG": self.threshold_sigma,
        "DLTMAXST": self.max_stamps,
        "DLTSAT": self.saturation,
        "DLTMRAD": self.match_radius,
        "DLTFRAD": self.fwhm_radius,
        "DLTBMASK": self.bright_mask,
        "DLTCLPSG": self.clip_sigma,
        "DLTCLPIT": self.clip_iterations,
        "DLTMINST": self.min_stamps,
        "DLTCVF": self.cv_folds,
        "DLTSSCL": self.spatial_scale if self.spatial_scale is not None else 0.0,
        "DLTDECOR": self.decorrelate,
        "DLTSCORE": self.score,
        "DLTSCAT": self.source_catalog,
        "DLTBLK": self.block,
    }

fit

fit(
    science,
    reference,
    *,
    science_var=None,
    reference_var=None,
    science_mask=None,
    reference_mask=None,
    gain=None,
    read_noise=0.0,
    catalog=None,
    direction=None,
) -> KernelSolution

Select stamps and fit the kernel/background, returning the solution.

Source code in python/delta/pipeline.py

def fit(
    self,
    science,
    reference,
    *,
    science_var=None,
    reference_var=None,
    science_mask=None,
    reference_mask=None,
    gain=None,
    read_noise=0.0,
    catalog=None,
    direction=None,
) -> KernelSolution:
    """Select stamps and fit the kernel/background, returning the solution."""
    solution, _ = self._fit(
        science,
        reference,
        science_var=science_var,
        reference_var=reference_var,
        science_mask=science_mask,
        reference_mask=reference_mask,
        gain=gain,
        read_noise=read_noise,
        catalog=catalog,
        direction=direction,
    )
    return solution

subtract

subtract(
    science,
    reference,
    *,
    science_var=None,
    reference_var=None,
    science_mask=None,
    reference_mask=None,
    gain=None,
    read_noise=0.0,
    catalog=None,
    direction=None,
) -> DiffResult

Fit and produce the difference (and optional decorrelation/score).

Source code in python/delta/pipeline.py

def subtract(
    self,
    science,
    reference,
    *,
    science_var=None,
    reference_var=None,
    science_mask=None,
    reference_mask=None,
    gain=None,
    read_noise=0.0,
    catalog=None,
    direction=None,
) -> DiffResult:
    """Fit and produce the difference (and optional decorrelation/score)."""
    start = time.perf_counter()
    solution, layers = self._fit(
        science,
        reference,
        science_var=science_var,
        reference_var=reference_var,
        science_mask=science_mask,
        reference_mask=reference_mask,
        gain=gain,
        read_noise=read_noise,
        catalog=catalog,
        direction=direction,
    )
    result = self._subtract(solution, layers)
    result.config = self.config_cards()
    result.elapsed = time.perf_counter() - start
    return result

apply

apply(
    solution: KernelSolution,
    science,
    reference,
    *,
    science_var=None,
    reference_var=None,
    science_mask=None,
    reference_mask=None,
    gain=None,
    read_noise=0.0,
) -> DiffResult

Apply a frozen solution to a new pair without re-fitting.

Runs only the full-frame subtraction (+ optional decorrelation/score) using the supplied :class:KernelSolution — the convolution direction, kernel basis, knots and coefficients all come from solution. The inputs must match solution.shape. Use this to reuse one fit across a stack of frames, or to re-apply a saved solution for QA/reproducibility.

Source code in python/delta/pipeline.py

def apply(
    self,
    solution: KernelSolution,
    science,
    reference,
    *,
    science_var=None,
    reference_var=None,
    science_mask=None,
    reference_mask=None,
    gain=None,
    read_noise=0.0,
) -> DiffResult:
    """Apply a frozen ``solution`` to a new pair without re-fitting.

    Runs only the full-frame subtraction (+ optional decorrelation/score)
    using the supplied :class:`KernelSolution` — the convolution direction,
    kernel basis, knots and coefficients all come from ``solution``. The
    inputs must match ``solution.shape``. Use this to reuse one fit across a
    stack of frames, or to re-apply a saved solution for QA/reproducibility.
    """
    start = time.perf_counter()
    layers = self._coerce_apply_layers(
        solution,
        science,
        reference,
        science_var=science_var,
        reference_var=reference_var,
        science_mask=science_mask,
        reference_mask=reference_mask,
        gain=gain,
        read_noise=read_noise,
    )
    result = self._subtract(solution, layers)
    result.config = self.config_cards()
    result.elapsed = time.perf_counter() - start
    return result

DiffResult dataclass

DiffResult(
    difference: NDArray[float32],
    variance: NDArray[float32] | None,
    mask: NDArray[uint8] | None,
    score: NDArray[float32] | None,
    solution: KernelSolution,
    psf: NDArray[float32] | None = None,
    source_catalog: NDArray | None = None,
    config: dict[str, object] = dict(),
    elapsed: float = float("nan"),
    variance_scale: float | None = None,
)

Products of a subtraction run.

header_cards

header_cards() -> dict[str, object]

Full provenance card set: fit + config + environment + runtime.

Combines :meth:KernelSolution.header_cards, the producing Subtractor's :meth:~Subtractor.config_cards, and :func:delta._provenance.environment_cards, enough on its own to reproduce this run (modulo the input pixel data itself).

Source code in python/delta/pipeline.py

def header_cards(self) -> dict[str, object]:
    """Full provenance card set: fit + config + environment + runtime.

    Combines :meth:`KernelSolution.header_cards`, the producing
    ``Subtractor``'s :meth:`~Subtractor.config_cards`, and
    :func:`delta._provenance.environment_cards`, enough on its own to
    reproduce this run (modulo the input pixel data itself).
    """
    return {
        **self.solution.header_cards(),
        **self.config,
        **environment_cards(),
        "DLTELAP": self.elapsed,
    }

write

write(
    path: str,
    overwrite: bool = False,
    *,
    compress: bool = True,
    extra_cards: dict[str, object] | None = None,
) -> None

Write products to a multi-extension FITS with provenance (needs astropy).

With compress (the default) each image layer is stored as a tile- compressed CompImageHDU (RICE), which shrinks noise-dominated difference frames by roughly an order of magnitude. A FITS primary HDU cannot itself be compressed, so the provenance header lives in a dataless primary and the difference moves to a DIFFERENCE extension. Pass compress=False for the plain layout (difference in the primary).

extra_cards adds/overrides header cards on top of :meth:header_cards — e.g. input file paths or the invoking command line, which this method has no way to know on its own.

Source code in python/delta/pipeline.py

def write(
    self,
    path: str,
    overwrite: bool = False,
    *,
    compress: bool = True,
    extra_cards: dict[str, object] | None = None,
) -> None:
    """Write products to a multi-extension FITS with provenance (needs astropy).

    With ``compress`` (the default) each image layer is stored as a tile-
    compressed ``CompImageHDU`` (RICE), which shrinks noise-dominated
    difference frames by roughly an order of magnitude. A FITS primary HDU
    cannot itself be compressed, so the provenance header lives in a
    dataless primary and the difference moves to a ``DIFFERENCE`` extension.
    Pass ``compress=False`` for the plain layout (difference in the primary).

    ``extra_cards`` adds/overrides header cards on top of
    :meth:`header_cards` — e.g. input file paths or the invoking command
    line, which this method has no way to know on its own.
    """
    from astropy.io import fits  # optional dependency

    cards = self.header_cards()
    if self.variance_scale is not None:
        cards["DLTVSCL"] = self.variance_scale
    if extra_cards:
        cards.update(extra_cards)
    layers = [("VARIANCE", self.variance), ("MASK", self.mask), ("SCORE", self.score)]
    hdus: list[fits.hdu.base._BaseHDU]
    if compress:
        primary = fits.PrimaryHDU()
        for key, value in cards.items():
            primary.header[key] = value
        hdus = [primary, fits.CompImageHDU(self.difference, name="DIFFERENCE")]
        hdus += [fits.CompImageHDU(a, name=n) for n, a in layers if a is not None]
    else:
        primary = fits.PrimaryHDU(self.difference)
        for key, value in cards.items():
            primary.header[key] = value
        hdus = [primary]
        hdus += [fits.ImageHDU(a, name=n) for n, a in layers if a is not None]
    if self.source_catalog is not None:
        hdus.append(fits.BinTableHDU(self.source_catalog, name="CATALOG"))
    fits.HDUList(hdus).writeto(path, overwrite=overwrite)

build_catalog

build_catalog(
    *,
    threshold_sigma: float = 5.0,
    threshold_sigma_dipole: float = 3.0,
    fwhm_tolerance_lo: float = 0.3,
    fwhm_tolerance_hi: float = 3.0,
    aperture_radius: int = 0,
    exclude_bad_pixels: bool = True,
    fit_psf_shape: bool = False,
    shape_chi2_tolerance: float = 3.0,
    bright_x: NDArray | None = None,
    bright_y: NDArray | None = None,
    bright_star_radius: float = 0.0,
    max_shape_fit: int = 5000,
    polarity: str = "positive",
    as_table: bool = False,
)

Rebuild the source catalog from self.score with custom thresholds (SPEC §3.7), instead of the defaults used to auto-populate self.source_catalog.

The expected PSF FWHM is measured from self.psf (the stamp built for the match filter); requires the originating subtract() / Subtractor.subtract() call to have used score=True.

fit_psf_shape=True adds a per-candidate PSF chi^2 shape fit (using self.psf and self.variance, so it stays consistent with what the matched filter itself scored with) -- see catalog.build_catalog for the false-positive rationale. bright_x/bright_y add a bright-star proximity flag; when omitted (the default) they're sourced from self.solution.stamp_x/ stamp_y (the bright stars used for the kernel fit) if available, since that's a free input -- pass empty arrays to opt out.

Source code in python/delta/pipeline.py

def build_catalog(
    self,
    *,
    threshold_sigma: float = 5.0,
    threshold_sigma_dipole: float = 3.0,
    fwhm_tolerance_lo: float = 0.3,
    fwhm_tolerance_hi: float = 3.0,
    aperture_radius: int = 0,
    exclude_bad_pixels: bool = True,
    fit_psf_shape: bool = False,
    shape_chi2_tolerance: float = 3.0,
    bright_x: NDArray | None = None,
    bright_y: NDArray | None = None,
    bright_star_radius: float = 0.0,
    max_shape_fit: int = 5000,
    polarity: str = "positive",
    as_table: bool = False,
):
    """Rebuild the source catalog from `self.score` with custom thresholds
    (SPEC §3.7), instead of the defaults used to auto-populate
    `self.source_catalog`.

    The expected PSF FWHM is measured from `self.psf` (the stamp built for
    the match filter); requires the originating ``subtract()`` /
    ``Subtractor.subtract()`` call to have used ``score=True``.

    ``fit_psf_shape=True`` adds a per-candidate PSF chi^2 shape fit (using
    ``self.psf`` and ``self.variance``, so it stays consistent with what
    the matched filter itself scored with) -- see
    ``catalog.build_catalog`` for the false-positive rationale.
    ``bright_x``/``bright_y`` add a bright-star proximity flag; when
    omitted (the default) they're sourced from ``self.solution.stamp_x``/
    ``stamp_y`` (the bright stars used for the kernel fit) if available,
    since that's a free input -- pass empty arrays to opt out.
    """
    if self.score is None or self.psf is None:
        raise ValueError(
            "build_catalog() requires a score image and PSF stamp; pass score=True "
            "to subtract()"
        )
    if fit_psf_shape and self.variance is None:
        raise ValueError("fit_psf_shape=True requires a variance map (self.variance)")
    from .catalog import build_catalog as _build_catalog

    if bright_x is None and bright_y is None and self.solution is not None:
        bright_x = self.solution.stamp_x
        bright_y = self.solution.stamp_y

    return _build_catalog(
        self.score,
        self.difference,
        mask=self.mask,
        psf=self.psf,
        threshold_sigma=threshold_sigma,
        threshold_sigma_dipole=threshold_sigma_dipole,
        fwhm_tolerance_lo=fwhm_tolerance_lo,
        fwhm_tolerance_hi=fwhm_tolerance_hi,
        aperture_radius=aperture_radius,
        exclude_bad_pixels=exclude_bad_pixels,
        fit_psf_shape=fit_psf_shape,
        shape_chi2_tolerance=shape_chi2_tolerance,
        variance=self.variance if fit_psf_shape else None,
        bright_x=bright_x,
        bright_y=bright_y,
        bright_star_radius=bright_star_radius,
        max_shape_fit=max_shape_fit,
        polarity=polarity,
        as_table=as_table,
    )