equivalent_width calculation: performance issues, and more

[cartilage-ftw]

I'm trying to calculate the equivalent width of a few absorption lines in a high res stellar spectrum my mentor provided me (size: 2.9 mb). I have a pretty simple piece of code for calculating equivalent of a Spectrum1D loaded from a fits file

import astropy.units as u
import matplotlib.pyplot as plt

from specutils import Spectrum1D
from specutils import SpectralRegion
from specutils.analysis import equivalent_width

spec = Spectrum1D.read('../data/arcturus.fits')

#plt.plot(spec.spectral_axis, spec.flux, 'k-')
#plt.axvline(x=5462.960, ls=':')
#plt.xlabel('Wavelength ($\AA$)')
#plt.ylabel('Normalized Flux')

def calc_width_milliaa(spectrum, start, end, continuum=1):
    """
    Parameters
    -----------
    spectrum: the Spectrum1D object
    continuum: the value of the (normalized) continuum above the absorption dip
    start: wavelength range start (in Angstroms)
    end: wavelength range end (in Angstroms)

    Returns
    ----------
    Equivalent width in milli Angstroms.
    """
    return 1000*equivalent_width(spec, continuum, regions=SpectralRegion(start*u.AA, end*u.AA))

calc_width_milliaa(spec, 5462.727, 5463.117, continuum=0.993)

But EW calculation of one line in the region described in the last line of the code above takes ~180 seconds (benchmarked). For the spectrum demonstrated in the specutils documentation, a similar calculation is.. almost instantaneous. There are some ~150 or so odd lines I'd like to calculate for a certain element. Is there something that I'm doing wrong or is it supposed to be this slow for higher resolution spectra?

Also, is there a way to actually see what it's doing, as you can in IRAF by interactively fitting a Gaussian/Voigt/Lorentzian as one wishes? (In either specutils or any other well known astropy package, or for that purpose, any modern astronomy software).

My package versions

astropy                       4.2.1
pytest-astropy                0.8.0
pytest-astropy-header         0.1.2
specutils                     1.2

Python version: Python 3.8.5

Let me know if any additional information is needed.

[keflavich]

Could you supply the .fits file so we can reproduce the issue?

[keflavich]

Regarding plotting and fitting line profiles, check out https://specutils.readthedocs.io/en/stable/fitting.html and https://docs.astropy.org/en/latest/api/astropy.modeling.functional_models.Voigt1D.html#astropy.modeling.functional_models.Voigt1D

[cartilage-ftw]

Regarding plotting and fitting line profiles, check out https://specutils.readthedocs.io/en/stable/fitting.html and https://docs.astropy.org/en/latest/api/astropy.modeling.functional_models.Voigt1D.html#astropy.modeling.functional_models.Voigt1D

Thank you so much for sharing the links. I will check these to see if I can get the same level of interactivity as IRAF offers. Sadly, even in 2021 some of the leaders in the field use the 35 years old tool IRAF because of the convenience it provides. I hope some day we will have a replacement for that.

I'm attaching a .zip containing the .fits file below.

arcturus.zip

[rosteen]

If you're looking for an interactive tool, I'd recommend checking out Jdaviz, which includes a configuration (Specviz) for working with spectra that is based on specutils.

Back to the issue at hand - I was able to confirm that I'm seeing bad performance as well, and was able to pinpoint it to the extract_region call that happens as part of the equivalent width calculation. I'm continuing to look into it...

Edit: every spectrum.wcs.world_to_pixel() call is taking 40 seconds which is...clearly not ideal.

[rosteen]

@eteq @keflavich Some more detail on the exact cause of the slowdown:

This is an IRAF file, which Spectrum1D.read() initializes from by creating a spectral_axis array and initializing the Spectrum1D object with that. When initializing with a spectral axis Spectrum1D creates a GWCS SpectralTabular1D lookup table for the WCS. Using world_to_pixel on this lookup table appears to scale with N, so 1 million data points takes 10x longer (~53 seconds in my test with a dummy spectrum) than for 100,000 data points (~5 seconds with a dummy spectrum).

I see a few possibilities here:

1. We could have Spectrum1D do something different to create a WCS when a Spectrum1D is initialized without one.
2. We could do something more manual/custom in extract_region in these cases that is faster than wcs.world_to_pixel
3. Fix whatever GWCS is doing here
4. Come up with a sorting algorithm that scales better than N and win some sort of prize and/or large sum of money (if I'm remembering algorithms correctly...)

Thoughts? I'm still kind of shocked that this takes so long. I mean (np.abs(spectral_axis - 5080*u.AA)).argmin() on a million data point spectral axis takes 10 ms.

[keflavich]

There are two solutions here: (1) if this is a linear WCS, we shouldn't be using a lookuptable. It's probably a linear WCS? (2) The scaling is not really a problem, I think (N is good scaling). The absolute time is worrying. 5 seconds for 100k is very, very slow. So I think point (3) above is the issue.

[eteq]

After some additional sleuthing, I think I see why this is ending up as a lookup table - the key point is here: https://github.com/astropy/specutils/blob/71b7f465a3523702f81dccf938a7b796bb723194/specutils/io/default_loaders/wcs_fits.py#L437 where basically the multispec format gets translated into a wavelength array. More specifically, https://github.com/astropy/specutils/blob/71b7f465a3523702f81dccf938a7b796bb723194/specutils/io/default_loaders/wcs_fits.py#L506 shows where the code creates a model and then uses that model to create the spectral_axis for lookup.

So I think "all" that's needed is to update https://github.com/astropy/specutils/blob/71b7f465a3523702f81dccf938a7b796bb723194/specutils/io/default_loaders/wcs_fits.py#L506 and a few other places in the same file to return a gwcs that uses the model itself instead of producing a lookup table. As per @keflavich's 1, this particular spectrum is not linear, but it's a pretty simple polynomial so should not require a lookup table.

That said, I think it's also surprisingly slow that a 100k point lookup table is that slow! So I think we should fix this immediate issue by my suggestion above but we should also try to fix the GWCS interpolation at some point (but that's probably a bigger job...).

[eteq]

Also one other piece of data for troubleshooting: pixel_to_world is much more reasonable - ~3 ms for the ~1mil element spectral axis. Still slower than I thought (np.searchsorted takes only a few microseconds), but still in the realm of "not a bottleneck", unlike 5 seconds

[eteq]

Aha! After a deep dive I think I have a solution for why the table is slow even thought I still don't fully understand the cause. I can reproduce the slowness minimally like so:

from astropy import modeling
pts = np.arange(1e6)*u.angstrom
lut = np.arange(1e6)
modeling.tabular.Tabular1D(pts, lut)(5000*u.angstrom)

that's what the code around https://github.com/astropy/specutils/blob/8c3a375/specutils/utils/wcs_utils.py#L207 is producing. Interestingly,

modeling.tabular.Tabular1D(pts.value, lut)(5000*u.angstrom)

is 83 ms, which is still a bit slow but much more reasonable. So there most be something strange happening in the quantity conversion. So I'm guessing the quickest thing there is to do a pre-conversion in the backwards transform...

[nmearl]

Here's some profiling. It certainly seems to be due to the unit machinery, and more specifically the world_to_pixel transformations (as can perhaps be seen more clearly in the call graph at the bottom).

+----------------------------------------------------+------------+-----------+---------------+
|                        Name                        | Call Count | Time (ms) | Own Time (ms) |
+----------------------------------------------------+------------+-----------+---------------+
| _get_physical_type_id                              |   12196496 |     37384 |         18335 |
| is_equivalent                                      |    9145058 |     81194 |         12649 |
| <function Quantity.__array_ufunc__ at 0x11c01fe50> |    3048388 |     10489 |         10489 |
| __array_ufunc__                                    |    3048388 |     41701 |          8505 |
| <built-in method builtins.getattr>                 |   42738968 |     13848 |          7978 |
| _normalize_equivalencies                           |    9145061 |     14329 |          7139 |
| __array_ufunc__                                    |    3048320 |     50491 |          6756 |
| __array_finalize__                                 |    6096650 |     18566 |          6547 |
| _new_view                                          |    3048473 |    116738 |          6164 |
| _is_equivalent                                     |    6096803 |     44461 |          6029 |
| __array_finalize__                                 |    6096650 |      9991 |          5691 |
| converters_and_unit                                |    3048388 |     14409 |          4937 |
| decompose                                          |    9149301 |      8377 |          4323 |
| <genexpr>                                          |    9144990 |     68938 |          4295 |
| <listcomp>                                         |   12196496 |      6142 |          4254 |
| <method 'view' of 'numpy.ndarray' objects>         |    9146497 |     12976 |          3889 |
| _set_unit                                          |    3048363 |     86755 |          3534 |
| <built-in method builtins.any>                     |    6102024 |     73168 |          3429 |
| to_value                                           |    6097169 |      4931 |          3251 |
| _normalize_equivalencies                           |    9145068 |      3239 |          3239 |
| decompose                                          |    9150012 |      4058 |          3038 |
| <built-in method numpy.array>                      |    3050360 |      2586 |          2586 |
| <built-in method builtins.max>                     |       1595 |     87592 |          2570 |
| <built-in method builtins.min>                     |      49734 |     86822 |          2553 |
| __array_finalize__                                 |    6097245 |      2544 |          2543 |
| equivalencies                                      |    9145068 |      2461 |          2461 |
| <built-in method builtins.isinstance>              |   18534946 |      2372 |          2370 |
| quantity_iter                                      |    3048358 |    118797 |          2060 |
| <method 'view' of 'numpy.generic' objects>         |    3048313 |      2041 |          2041 |
| name                                               |   12212052 |      1890 |          1890 |
+----------------------------------------------------+------------+-----------+---------------+