fastqso -- fast QSO fitter

[moustakas]

[WIP]

@paulmartini @stephjuneau @arjundey @yzfkieran817 @abrodze and many others

At long last, I have an initial implementation of a "fast" QSO fitter called fastqso, building on work that I started two years ago (e.g., the templates have had AGN models in them for some time) and the "wish list" in #112.

Quoting from the new continuum.qso_continuum_fastspec docstring:

    The QSO continuum model has three components:

    - Power law: F_λ ∝ (λ/1450 Å)^(α-2) with IGM attenuation applied, where α
      is the spectral index in F_ν ∝ ν^{-α}.  The fitted amplitude A_PL is the
      observed flux at 1450 Å rest-frame in units of 10^{-17} erg/s/cm²/Å.
    - UV Fe emission: the Vestergaard & Wilkes (2001) template convolved to
      velocity dispersion σ_Fe, resampled to the template wavelength grid.
    - IR AGN torus: the Nenkova+08 template; its amplitude A_torus is derived
      from energy balance (absorbed UV/optical luminosity re-emitted in the IR).

    Fitting proceeds in four steps:

    1. σ_Fe is determined by a chi2 grid scan restricted to the Fe-template
       wavelength window (1075–3090 Å rest), with two nuisance power-law columns
       marginalizing over the continuum.  The minimum is refined with a 1D
       parabola fit in log-σ space.
    2. The aperture correction is estimated from the ratio of broadband
       photometry to a dust-free spectral model at α=1.5.
    3. α and τ_V are optimized jointly via VARPRO: an outer grid on α (11 points
       over [−0.5, 3.0]) with an inner Brent scalar minimization over τ_V ∈ [0, 5].
       For each (α, τ_V) evaluation, dust attenuation and energy-balance torus
       emission are applied via :meth:`ContinuumTools.attenuate`, and (A_PL, A_Fe)
       are solved by NNLS on the joint spectro+photometric system.
    4. A final NNLS solve at (α_best, τ_V,best) yields the model arrays and the
       torus amplitude from energy balance.

    All parameter uncertainties are estimated by Monte Carlo: the full four-step
    pipeline is repeated for ``nmonte`` noise realizations drawn from the spectral
    and photometric inverse variances.

Here's an example run (6.79 sec for 50 Monte Carlo fitting iterations and 0.84 sec with no Monte Carlo):

% time fastqso -o fqso.fits /dvs_ro/cfs/cdirs/desi/spectro/redux/loa/healpix/main/dark/323/32363/redrock-main-dark-32363.fits --targetids=39628262946641205
INFO:fastspecfit.py:114:parse: fastspec -o fqso.fits /dvs_ro/cfs/cdirs/desi/spectro/redux/loa/healpix/main/dark/323/32363/redrock-main-dark-32363.fits --targetids=39628262946641205
INFO:fastspecfit.py:325:fastspec: fsftime 0.82 sec for sc_data_init at 2026-05-24T12:59:22.780
INFO:fastspecfit.py:338:fastspec: Cached stellar templates /dvs_ro/cfs/cdirs/desi/public/external/templates/fastspecfit/2.0.0/ftemplates-chabrier-2.0.0.fits
INFO:fastspecfit.py:339:fastspec: Cached emission-line table /global/common/software/desi/users/ioannis/fastspecfit/py/fastspecfit/data/emlines.ecsv
INFO:fastspecfit.py:340:fastspec: Cached photometric filters and parameters /global/common/software/desi/users/ioannis/fastspecfit/py/fastspecfit/data/legacysurvey-dr9.yaml
INFO:fastspecfit.py:341:fastspec: Cached cosmology table /global/common/software/desi/users/ioannis/fastspecfit/py/fastspecfit/data/desi_fiducial_cosmology.dat
INFO:fastspecfit.py:342:fastspec: Cached Inoue+2014 IGM attenuation parameters.
INFO:io.py:726:gather_metadata: specprod=loa, coadd_type=healpix, survey=main, program=dark, healpix=32363
INFO:io.py:1057:read: Reading 1 spectrum from /dvs_ro/cfs/cdirs/desi/spectro/redux/loa/healpix/main/dark/323/32363/coadd-main-dark-32363.fits
INFO:io.py:1165:read: fsftime 0.34 sec for read_spectra [nobj=1, nfiles=1] at 2026-05-24T12:59:26.365
INFO:fastspecfit.py:192:fastspec_one: Continuum- and emission-line fitting object 0 [targetid 39628262946641205, seed 2032329983, z=1.799320].
INFO:continuum.py:2117:qso_continuum_fastspec: Median spectral S/N_b=5.38 S/N_r=7.62 S/N_z=7.32
INFO:continuum.py:2562:continuum_specfit: Smooth continuum correction: b=4.319% r=-0.867% z=4.833%
INFO:continuum.py:2786:continuum_specfit: pl_slope=0.55, pl_amplitude=6.85+/-0.0698 (10^-17 erg/s/cm2/A at 1450 A), torus_amplitude=1.28e+04+/-1.79e+03
INFO:continuum.py:2791:continuum_specfit: fe_vdisp=1927+/-475 km/s, fe_amplitude=0.0301+/-0.00414, tauv=0.017+/-0.004
INFO:continuum.py:2796:continuum_specfit: Median aperture correction 1.247 [1.049-1.904].
INFO:emlines.py:1681:emline_specfit: Skipping broad-line fitting (no broad Balmer lines in the spectral range).
INFO:emlines.py:1735:emline_specfit: delta(v) UV=-476.7+/-1027.9 km/s Balmer broad=0.0 km/s narrow=-256.1 km/s
INFO:emlines.py:1745:emline_specfit: sigma UV=2455+/-1247 km/s Balmer broad=0 km/s narrow=86 km/s
INFO:fastspecfit.py:236:fastspec_one: fsftime 6.36 sec for fastspec_one [targetid=39628262946641205] at 2026-05-24T12:59:32.771
INFO:fastspecfit.py:477:fastspec: fsftime 6.79 sec for fit_all [nobj=1,6.79s/obj/core] at 2026-05-24T12:59:33.194
INFO:io.py:1662:write: Writing 1 object to fqso.fits
INFO:io.py:1733:write_fastspecfit: fsftime 1.83 sec for write_fastspecfit [file=fqso.fits] at 2026-05-24T12:59:35.024

And here's the current version of the QA:

For comparison, here's the same object fitted with fastspec and the nominal set of SPS models:

This definitely isn't ready for primetime---and there's a lot I'm going to need some help with---but initial reactions are welcome.

[coveralls]

Coverage Report for CI Build 26371427444

Coverage decreased (-3.6%) to 62.118%

Details

• Coverage decreased (-3.6%) from the base build.
• Patch coverage: 435 uncovered changes across 4 files (133 of 568 lines covered, 23.42%).
• 2 coverage regressions across 2 files.

Uncovered Changes

File	Changed	Covered	%
py/fastspecfit/continuum.py	355	56	15.77%
py/fastspecfit/qa.py	145	32	22.07%
py/fastspecfit/io.py	50	28	56.0%
py/fastspecfit/fastspecfit.py	3	2	66.67%

Coverage Regressions

2 previously-covered lines in 2 files lost coverage.

File	Lines Losing Coverage	Coverage
py/fastspecfit/io.py	1	70.14%
py/fastspecfit/qa.py	1	62.11%

---

Coverage Stats

Relevant Lines:	7365
Covered Lines:	4575
Line Coverage:	62.12%
Coverage Strength:	0.62 hits per line

---

💛 - Coveralls

[moustakas]

One important unresolved issue is that the IR photometry may not have enough signal-to-noise to break the degeneracy between the UV/optical power law parameterization and the dust attenuation. For example, in this QSO the WISE photometry is disregarded entirely--

[Swayamtrupta]

1. Possible inclusion of a Balmer continuum component (as an optional parameter, similar to PyQSOFit). PyQSOFit uses the approach from Dietrich et al. (2002). Here is the code snippet from the pyqsofit.py:

def Balmer_conti(self, xval, pp):
        """Fit the Balmer continuum from the model of Dietrich+02"""
        # xval = input wavelength, in units of A
        # pp=[norm, Te, tau_BE] -- in units of [--, K, --]
        xval = xval * u.AA
        lambda_BE = 3646.  # A
        bb_lam = BlackBody(pp[1] * u.K, scale=1.0 * u.erg / (u.cm ** 2 * u.AA * u.s * u.sr))
        bbflux = bb_lam(xval).value * 3.14  # in units of ergs/cm2/s/A
        tau = pp[2] * (xval.value / lambda_BE) ** 3
        result = pp[0] * bbflux * (1 - np.exp(-tau))
        ind = np.where(xval.value > lambda_BE, True, False)
        if ind.any() == True:
            result[ind] = 0
        return result

2. Another aspect (which isn't clear to me, sorry) is how the choice of number of components (esp. for broad lines) is made in fastqso. Would a BIC-like criterion help to break the degeneracies?
3. A parameter I like to look at is frac_host (F_host/F_AGN) or frac_AGN at some specific wavelength (or wavelength ranges). For example, PyQSOFit reports this at 5100 Angstroms when the host fitting is on.