Native vs EasySpec showdown¶

The whole point of the shelf is that for every operation you have a choice of backends - the native spectro-kernel implementation and (where available) a wrapper around a pro library like EasySpec. This walk-through runs the same operations through both paths on the same input, and shows you how to read the disagreement.

Pick whichever side you trust more after this notebook. Or run both in production and cross-check on suspect frames.

1. The synthetic test set¶

A 3-step setup: build three bias frames, three darks, three flats, plus a science frame, all as on-disk FITS so both paths see the same bytes.

import os, tempfile
import numpy as np
from astropy.io import fits

WORK = tempfile.mkdtemp(prefix="shelf-shootout-")
os.makedirs(os.path.join(WORK, "biases"))
os.makedirs(os.path.join(WORK, "darks"))
os.makedirs(os.path.join(WORK, "flats"))

def _save_set(prefix, level, noise, n=3, shape=(64, 64)):
    rng = np.random.default_rng(hash(prefix) & 0xFFFFFFFF)
    for i in range(n):
        data = rng.normal(level, noise, shape).astype(np.float64)
        path = os.path.join(WORK, prefix, f"{prefix.rstrip('s')}_{i}.fits")
        fits.PrimaryHDU(data).writeto(path)

_save_set("biases", level=10.0, noise=0.5)
_save_set("darks",  level=2.0,  noise=0.3)
_save_set("flats",  level=1000.0, noise=10.0)

rng = np.random.default_rng(0)
science = (
    rng.normal(80.0, 1.5, (64, 64)) + 10.0    # signal + bias
).astype(np.float64)
science_path = os.path.join(WORK, "science.fits")
fits.PrimaryHDU(science).writeto(science_path)

2. Build the master bias - both ways¶

from spectro_kernel import WorkContext, run_algorithm
from spectro_kernel.io import read_fits
from spectro_kernel.types import ImageFrame
import glob, numpy as np

def _load_frames(paths):
    return [
        ImageFrame(
            data=np.asarray(fits.getdata(p), dtype=np.float64),
            frame_type="bias",
        )
        for p in paths
    ]

bias_paths = sorted(glob.glob(os.path.join(WORK, "biases", "*.fits")))

# --- native (bias_combine reads ctx.images, median by default) ----
ctx_native = WorkContext(images=_load_frames(bias_paths))
run_algorithm("bias_combine", ctx_native, {"method": "median"})
native_master = ctx_native.extras["master_bias"]

# --- easyspec wrapper ---------------------------------------------
ctx_easy = WorkContext()
run_algorithm("bias_combine_easyspec", ctx_easy, {"bias_paths": bias_paths})
easy_master = ctx_easy.extras["master_bias"]

print(f"native median = {float(np.median(native_master.data)):.3f}")
print(f"easy   median = {float(np.median(easy_master.data)):.3f}")
print(f"|delta| pixelwise (median) = {float(np.median(np.abs(native_master.data - easy_master.data))):.4f}")

For a uniform median stack on the same input, the two masters should agree to the last digit. If they do not, look at:

Method - easyspec's mean differs from the native mean only by edge cases (NaN handling, integer truncation).
Trim - easyspec optionally trims to a region of interest; the native algorithm does not. If you crop manually before, the masters match.

3. Run the SNR comparator on the science frame¶

sci = read_fits(science_path)
ctx_sci = WorkContext(spectrum=type(sci)(
    np.arange(sci.npix, dtype=np.float64) if hasattr(sci, "npix") else
    np.arange(64*64, dtype=np.float64),
    sci.flux if hasattr(sci, "flux") else science.ravel(),
))
# Note: bias_combine produced an ImageFrame; for SNR we'd usually run on a
# 1D extracted spectrum. The point here is to illustrate compare_snr_methods.

For a proper SNR comparison, run the comparator on a real 1D extraction:

ctx_sci = WorkContext(spectrum=read_fits("/data/already_reduced_1d.fits"))
run_algorithm("compare_snr_methods", ctx_sci)
print(ctx_sci.extras["snr_methods"])
print(f"spread = {ctx_sci.metrics['snr_methods_spread']:.2f}")

Typical output: the three SNR estimates land within 5-10% of each other on a clean spectrum. A spread above 30% usually means one of:

A strong line is inside the linear-fit window.
The edge SNR window catches the noisy ends of the detector.
The DER_SNR median is being pulled by a low-flux dip.

That disagreement is itself diagnostic.

4. Three line profiles, one comparator¶

ctx_line = WorkContext(spectrum=read_fits("/data/some_star_with_halpha.fits"))
run_algorithm("normalize_polynomial", ctx_line, {"order": 3})
run_algorithm(
    "compare_line_fits",
    ctx_line,
    {"line_center_angstrom": 6562.79, "window_angstrom": 40.0},
)

comp = ctx_line.extras["line_fit_comparison"]
for name, fit in comp["results"].items():
    print(f"{name:25s} R²={fit['r_squared']:.4f}  FWHM={fit['fwhm_angstrom']:.2f} Å")
print(f"best profile: {comp['best_profile']}")

On a stellar absorption line, a Voigt usually wins by 1-3% in R² over a pure Gaussian thanks to the wings; on instrument-broadened nebular emission a Gaussian often suffices. The comparator removes the guesswork.

5. Reduction chain - native algo where it exists, EasySpec elsewhere¶

The natural pattern: pick native algorithms when they exist (you save on the EasySpec overhead - the staging directory, the subprocess-like state, the plot suppression), reach for the EasySpec wrapper for steps that do not have a native equivalent yet (extract_spectrum_easyspec, wavelength_calibrate_easyspec, flux_calibrate_easyspec).

ctx = WorkContext(images=_load_frames(bias_paths))
run_algorithm("bias_combine", ctx, {"method": "median"})              # native
# load_more_images...
run_algorithm("dark_combine_easyspec", ctx, {"dark_dir": "/data/darks"})  # easyspec (no native)
# … flat combine, applies, extract …

A ProcessingStep is recorded in ctx.history for every step regardless of backend, so the provenance is preserved. Grep ctx.history for the algorithm names you used; you'll see them interleaved.

What you have gained¶

A practical answer to the "native or easyspec?" question, recorded directly in your audit trail. The shelf is not a marketing choice - it's a way to run both implementations on the same bytes when the result matters, and pick the one that matches your reference.

Going further¶

Apply the same pattern to your own pipeline: every time you have an EasySpec wrapper and a native algorithm for the same operation, write the two-line shootout. Commit the result as a unit test (assert delta < tolerance) so future versions can't drift.
Add a custom backend (a fork of EasySpec, a different community library) by writing one more wrapper file - it joins the shelf automatically.
Pair with Claude + MCP end-to-end: ask the agent "compare the native and easyspec implementations of bias_combine on these frames and tell me if they agree" - it does.