Speed up loading of spectra from atmospheric grids.

[jluastro]

Description

Request: Speed up grid_to_spec() (or _get_spectrum()).

We have a code, SPISEA, that previously depended on PYSYNPHOT's ability to quickly read in spectra from the atmospheric model grids. Using psynphot 2.0 (on top of numpy 1.26), we could use pysynphot.Icat to load a spectrum in 0.09 seconds. (see code output below on python 3.10)

pysynphot k93models Icat + sample: 0.088233 s
pysynphot version: 2.0.0
numpy version: 1.26.4

If we upgrade python, numpy, and spisea to use stsynphot and synphot instead, we must switch to using stsynphot.grid_to_spec() instead. When we do this, it now takes 0.22 seconds to read a spectrum.... twice as slow. (see code output below on python 3.12)

stsynphot k93models grid_to_spec + flux evaluation: 0.217748 s
stsynphot version: 0.1.dev512+gf25255b3e.d20260428
numpy version: 2.4.3

I've tracked this down to two factors:

1. improved interpolation (8 spectra are loaded and then interpolated to give the final at specified values of Teff, metallicity, and log-g);
2. the support for astropy units in the new spectral objects.

On 1 -- I don't think we want to change this. So 8 spectra will still need to be loaded.

On 2 -- The addition of units is nice. The real slow down is in parsing the FITS spectral files treating their column header units as strings. There are a bunch of calls astropy.units.format.fits.parse, astropy.utils.misc.did_you_mean and difflib.get_close_matches -- ends up being 88 calls all together. I believe this is where most of the additional time is coming from.

Note that the real issue might be down in synphot (or even astropy). But I'm wondering if there is something that could be done at the stsynphot level to speed things up given that we know 100% in advance all the units of the atmospheric grid files.

I've attached a short test script below.

Additional context

Test script for timing.

def test_pysynphot_vs_stsynphot_timing(with_pysyn=None, with_stsyn=None):
    """
    Time Kurucz ``k93models`` spectrum extraction via stsynphot or pysynphot.

    Prints elapsed seconds from ``time.perf_counter()``; use ``pytest -s`` to show output.
    """
    temperature = 20000
    metallicity = 0.0
    gravity = 4.0

    if with_pysyn is None:
        with_pysyn = importlib.util.find_spec("pysynphot") is not None
    if with_stsyn is None:
        with_stsyn = importlib.util.find_spec("stsynphot") is not None
    
    if with_stsyn:
        if not os.environ.get("PYSYN_CDBS"):
            pytest.skip("PYSYN_CDBS not set; grid_to_spec needs CDBS tree")

        import stsynphot
        from stsynphot.catalog import grid_to_spec

        t0 = time.perf_counter()
        try:
            sp = grid_to_spec("k93models", temperature, metallicity, gravity)
            w = sp.waveset
            _ = sp(w)
        except Exception as exc:
            pytest.skip(f"stsynphot grid_to_spec failed: {exc}")
        elapsed = time.perf_counter() - t0
        print(f"stsynphot k93models grid_to_spec + flux evaluation: {elapsed:.6f} s")
        print(f"    stsynphot version: {stsynphot.__version__}")

        assert elapsed >= 0
        assert np.isfinite(elapsed)
        assert w.size > 0

    if with_pysyn:
        import pysynphot

        t0 = time.perf_counter()
        try:
            sp = pysynphot.Icat("k93models", temperature, metallicity, gravity)
            wave = sp.GetWaveSet()
            _ = sp.sample(wave)
        except Exception as exc:
            pytest.skip(f"pysynphot Icat failed: {exc}")
        elapsed = time.perf_counter() - t0
        print(f"pysynphot k93models Icat + sample: {elapsed:.6f} s")
        print(f"    pysynphot version: {pysynphot.__version__}")

        assert elapsed >= 0
        assert np.isfinite(elapsed)
        assert np.asarray(wave).size > 0

    print(f"   numpy version: {np.__version__}")    

    return

[pllim]

A proper profiling (https://docs.python.org/3/library/profile.html) might tell you exactly where the bottlenecks are. If https://jiffyclub.github.io/snakeviz/ still works, it builds on Python profiling and provides more user-friendly visualization.

You might also want to try with Python 3.13 or Python 3.14. Python 3.13 was advertised as performance release, though not sure how that helps in this regard, but doesn't hurt.

Thanks!

[pllim]

There are a bunch of calls astropy.units.format.fits.parse, astropy.utils.misc.did_you_mean and difflib.get_close_matches

This hinted at data files TUNITn not conforming to FITS standards.

synphot/stsynphot might also have more data validation than PySynphot, which also might cause extra run time. But then again, it is hard to tell without looking at profiling results.