Unify XSpectrum types

[icweaver]

Doc preview: https://juliaastro.org/Spectra.jl/previews/PR24/

Closes: #25, #37

What's changed

• Following the architecture of specutils, our Spectrum, EchelleSpectrum, and a new IFUSpectrum for spectral cube support, are all just Spectrum{S, F, M, N} <: AbstractSpectrum{S, F}s now, with their dimensionality specified in the type domain:

  const SingleSpectrum = Spectrum{S, F, 1, 1} where {S, F}
  const EchelleSpectrum = Spectrum{S, F, 2, 2} where {S, F}
  const IFUSpectrum = Spectrum{S, F, 1, 3} where {S, F}

  where S is the eltype of the spectral axis with dimensionality M, and F is the eltype of the flux axis with dimensionality N.

• More of the array interface is implemented now so things like this will work:

  julia> spec_1D = spectrum([20, 40, 120, 160, 200], [1, 3, 6, 7, 20])
  SingleSpectrum(Int64, Int64)
    spectral axis (5,): 20 .. 200
    flux axis (5,): 1 .. 20
    meta: Dict{Symbol, Any}()
  
  julia> spec_1D[2:4]
  SingleSpectrum(Int64, Int64)
    spectral axis (3,): 40 .. 160
    flux axis (3,): 3 .. 7
    meta: Dict{Symbol, Any}()
  
  julia> spec_1D[2:4] |> spectral_axis
  3-element Vector{Int64}:
    40
   120
   160

• Merged common.jl into Spectra.jl (for my own sanity)

• Renamed the Analysis page of the docs to Utilities and moved blackbody over to it

• Shelved some analysis functionality (e.g., continuum fitting and equivalent width support) in favor of using downstream packages like SpectralFitting.jl

• Generalized terms in preparation for introducing X-ray support in a future PR. e.g., wave --> spectral_axis and flux --> spectral-axis. These are now exported. Wavelength and energy now both referenced instead of just wavelength

• Filled in some more API documentation

To-do

• array ops
• better (more?) show methods
• begin generalizing for X-ray spectra
• tests
• docs

Future PR(s)

• enforce that spectral axes are compatible before doing algebra between multiple Spectrums Spectra algebra #43
• build out X-ray support Unify XSpectrum types #24 (comment)
• update plot recipes Explore moving out units and/or plotting #31
• delete or move out old code (commented out for now) Trim dead code #42

[codecov]

Codecov Report

❌ Patch coverage is 88.23529% with 14 lines in your changes missing coverage. Please review.
✅ Project coverage is 90.00%. Comparing base (88edb4d) to head (d147e6a).
⚠️ Report is 34 commits behind head on main.

Files with missing lines	Patch %	Lines
src/spectrum_ifu.jl	60.00%	8 Missing ⚠️
src/Spectra.jl	94.91%	3 Missing ⚠️
src/spectrum_echelle.jl	90.00%	2 Missing ⚠️
src/plotting.jl	66.66%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             main      #24       +/-   ##
===========================================
+ Coverage   74.28%   90.00%   +15.71%     
===========================================
  Files           8        7        -1     
  Lines         140      150       +10     
===========================================
+ Hits          104      135       +31     
+ Misses         36       15       -21     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[cgarling]

We did this in PhotometricFilters.jl, see here.

Is it useful for getindex to return a Spectrum object and the metadata? I would expect something like

Base.getindex(f::Spectrum, i::Int) = (wave(f)[i], flux(f)[i])

[icweaver]

Thanks! Yea, I liked how this keeps things consistent with EchelleSpcetrum's behavior.

Actually...

This may go nowhere, but I think I want to start investigating #25 (comment) now

[icweaver]

Ok, I think this might be going somewhere. Updated the parent comment with a summary

[cgarling]

If you want to cover specutils' case 2 (multi-dimensional flux, but same spectral axis along each) you may need to add a type parameter for the dimensionality of the wavelength array as well. I think case 2 is for IFU data (??) so maybe

mutable struct Spectrum{W<:Number, F<:Number, N, M} <: AbstractSpectrum{W, F}
    wave::AbstractArray{W, N}
    flux::AbstractArray{F, M}
    meta::Dict{Symbol,Any}
end

const SingleSpectrum = Spectrum{W, F, 1, 1} where {W, F}
const IFUSpectrum = Spectrum{W, F, 1, 2} where {W, F}
const EchelleSpectrum = Spectrum{W, F, 2, 2} where {W, F}

We can also worry about this later if you want

[icweaver]

oh I like this a lot. Added some methods to start playing with this. There's something really satisfying about being able to do, say:

julia> spec_IFU = spectrum([20, 40, 120, 160, 200], rand(10, 5))
IFUSpectrum(Int64, Float64)
  # orders: 10
  wave: (20, 200)
  meta: Dict{Symbol, Any}()

julia> spec_IFU[3, begin:4]
SingleSpectrum(Int64, Float64)
  wave: (20, 200)
  flux: (0.12659667784419948, 0.3942549580490494)
  meta: Dict{Symbol, Any}(:Order => 3)

to get the first 4 flux measurements from the 3rd order of an IFUSpectrum. I imagine this could be opened up a bit more to allow for more uniform indexing across the other spectrum types. I don't do space rainbow science much these days, so does this sound like something that would be useful to have, or is this trying to do too much?

[cgarling]

Actually I think an IFU with a uniform wavelength (dispersion) sampling would be a 3-D cube (2 spatial axes across the IFU field, one wavelength axis) so M = 3, see here for example.

I think it's a good idea, not sure what the best API is. I would probably put the wavelength axis first in the flux matrix because the most common operation is probably selecting a 1-D spectrum from the cube based on pixel.

julia> wave, flux = [20, 40, 120, 160, 200], rand(5, 10, 10);

julia> spec_ifu = spectrum(wave, flux);

julia> spec_ifu[:, 1, 1] # Selects full 1-D spectrum from spatial pixel with index (1, 1)

My intuition is that a scalar wavelength index would return a basic array (i.e, an image of the flux at a fixed wavelength).

julia> spec_IFU[3, begin:4, begin:4] isa Matrix{Float64}
true

julia> spec_IFU[3, begin:4, begin:4] == flux[3, begin:4, begin:4]
true

I think returning a spectrum type only makes sense if the index is a range for the wavelength. I'm really not sure though, I don't really work on this either

[icweaver]

Oh cool! Thanks for the reference, I've never worked with IFUs before, so this was a really interesting read. I think "spaxel" is my new favorite word 😂. I see that I was mixing concepts between echelles and ifus when I was talking about orders, sorry about that. The split to 3D really helped clarify things for me. How does this new design look? I also cleaned up the show methods a bit (still need to update tests once this settles though)

julia> using Spectra

julia> wave, flux = [20, 40, 120, 160, 200], rand(5, 10, 6);

julia> spec_ifu = spectrum(wave, flux)
IFUSpectrum(Int64, Float64)
  wave ((5,)): 20 .. 200
  flux ((5, 10, 6)): 0.6211135421955702 .. 0.4976591674645773
  meta: Dict{Symbol, Any}()

julia> spec_ifu[:, 1, 1]
SingleSpectrum(Int64, Float64)
  wave ((5,)): 20 .. 200
  flux ((5,)): 0.6211135421955702 .. 0.9046136130205523
  meta: Dict{Symbol, Any}()

julia> spec_ifu[:, begin:4, begin:3]
IFUSpectrum(Int64, Float64)
  wave ((5,)): 20 .. 200
  flux ((5, 4, 3)): 0.6211135421955702 .. 0.8211549470853685
  meta: Dict{Symbol, Any}()

julia> spec_ifu[3, begin:4, begin:3]
4×3 Matrix{Float64}:
 0.184341  0.942616   0.546656
 0.206153  0.754342   0.120533
 0.186781  0.0043856  0.0914326
 0.601994  0.352066   0.82047

julia> spec_ifu[3, begin:4, begin:3] == flux[3, begin:4, begin:3]
true

[fjebaker]

Hello @icweaver ! I've had a look at this PR, and include below some thought about how this would work with SpectralFitting.jl. These are just comments and observations towards some kind of base spectra package, but don't let me stop you with this PR.

TL;DR

For this to be useful for SpectralFitting.jl, I'd need at the very least:

# Returns the domain of the spectrum in units of energy (eV or KeV)
# This may have to be bin widths for X-ray too
Spectra.energy(::AbstractSpectrum)

It would probably also be good to have an interface for converting the units of both wavelength / energy, but also any (specific) fluxes. There's a table for how to do the flux conversions I have printed out on my desk that I can share here too (once I find the PDF of it again).

There's a few other caveats, but I can sort those within SpectralFitting.jl too if I were to depend on Spectra.jl. See below for those.

Details

I keep things brief here. X-ray astronomy (and other high energy fields) typically operate on spectra that are binned into channels by the measurement device. We usually deal with photon counts (as we don't observe many photons) that we can measure the energy (from the discrete channel) of individually.

The data is typically arranged with three vectors (after some processing):

• $E_{\text{min},i}$, the low energy of bin $i$.
• $E_{\text{max},i}$, the high energy of bin $i$. It is generally the case that $E_{\text{max},i} = E_{\text{min},i+1}$.
• $F_{i}$, the flux, usually in either photon counts per second per cm^2 of the detector effective area.

There's usually also an instrument response matrix that is carried around with the raw channel data that can be used to calculate the above three vectors, and is also needed for modelling purposes (akin to a PSF but which tabulates the probability of detection of a photon of true energy $E$ in a channel $C(E)$ ).

So the proposed changes in this PR aren't really suitable for SpectralFitting's alignment with X-ray astronomy at the moment, but would be perfect for other wavelengths.

For X-rays, we need spectra where we can:

• Have the domain be in units of energy (typically ev or keV or higher).
• Be able to represent contiguously binned data as a spectrum, so where there are bin edges instead of single wavelengths associated with each flux measurement.
• Be flexible about units being in energy flux or photon count flux.

There's also a point about uncertainties in the spectra, as X-ray have (co)variances that are not necessarily Gaussian, but that could be done by convention in the metadata.

In SpectralFitting, I currently have an AbstractLayout type that is concretely implemented as either ContiguouslyBinned() or OneToOne() (poorly named, there's an open issue about that) that is compile time queryable via the Holy Trait's pattern. Things can then check if a spectrum is one thing or the other, and approximately convert it to the type of data they want to deal with.

To make an example:

• I have an X-ray spectrum that is contiguously binned and want to use it in an SED fit. There is an automatic system that converts the X axis into the average energy in each bin and sets the standard deviation to the bin width.
• Vice versa.

Something like that would be very useful in Spectra.jl. Not everything would need to support both types of data, and to some extent that conversion can be automated provided it can be done without losing information.

I'd also encourage not using terms like wave and flux for the fields. Flux is perhaps alright, but wave would not extend to X-ray data or to multi-dimensional time-series spectra either. In SpectralFitting.jl I've used domain instead, but that's only because I couldn't think of a better term. I'm planning to use spectral_axis and flux_axis going forwards.

Virtual Observatory and SSA

I don't know if we want to also take a leaf out of the VO's book and make our spectral types directly based on the Spectral Data Model. See, for example, 3.1 of SpectrumDM. They obviously are much more complete than we would need to be to get started, but they've thought through most of the problems we would need to face in a base spectra package. The table in 4.2 in flux (as spectral intensity) and all the metadata and units that would need to include is a good thing to look at.

To reiterate what I said in the call though -- I would love a base spectra package, and I hope this doesn't look like I'm trying to poke holes in what you're doing here. Let me know if there's concrete ways you'd like me to contribute to this (very happy to open PRs also) so that we can get this moving 👍

[fjebaker]

Sorry after writing the above I noticed that wave and flux can have different dimensions, so that would sort the bin low / high bit very nicely too!

[cgarling]

Appreciate your thorough review Fergus! What you're doing with SpectralFitting.jl looks great.

Just FYI for everyone involved, I will be very busy for the next month or two at least so I can't commit to doing much for this right now. I also don't really do observational spectroscopy myself (my interest is in loading/manipulating theoretical SEDs) so I am not really a great domain expert here anyhow.

[icweaver]

For sure, thanks for all of your current input, Chris!

And Fergus, no apologies necessary at all. This is exactly the kind of expert feedback we need, thank you! I'm down to see if we can get some X-ray support going in this PR if you think there might be something there.

Funnily enough, the XSpectrum title I used for this PR as a placeholder for Single, Echelle, and IFU spectrums looks to be pretty similar to the name of the X-ray object that specutils is also working on. I hope that means we might be on to something, lol.

I'm still getting up to speed, but would a design like this be starting to go in the direction that you are thinking?

julia> using Spectra
Precompiling Spectra finished.
  1 dependency successfully precompiled in 3 seconds. 114 already precompiled.

julia> const XraySpectrum = Spectrum{S, F, 2, 1} where {S, F}
Spectrum{S, F, 2, 1} where {S, F}

julia> spec = Spectrum([2:2:10;; 4:2:12], 1:5, Dict())
Spectrum{Int64, Int64, 2, 1}([2 4; 4 6; … ; 8 10; 10 12], 1:5, Dict{Symbol, Any}())

julia> energy_low, energy_high = eachcol(spectral_axis(spec))
2-element ColumnSlices{Matrix{Int64}, Tuple{Base.OneTo{Int64}}, SubArray{Int64, 1, Matrix{Int64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 [2, 4, 6, 8, 10]
 [4, 6, 8, 10, 12]

julia> counts = flux_axis(spec)
1:5

Or maybe even subtypee XraySpectrum further so that it can support your ContiguouslyBinned and OneToOne cases. Would it make sense to store the response matrix separately in the meta field? Sort of a separate issue, but I guess we could also switch from dicts to named tuples too. Could we also use the Holy Trait design you mentioned to help validate construction (i.e., checking axis sizes make sense for the relevant spectrum type)?

At any rate, I quite like your suggested wave => spectral_axis and flux => flux_axis name change. Besides being more general, it also reduces the chance of us cluttering the namespace with common names if we wanted to export them.

I'd be happy to help support any additional PRs you might have in mind if I can. I've started working with SpectralFitting.jl, and working with the sample s54405.pha dataset has been really helpful!

[fjebaker]

I'm still getting up to speed, but would a design like this be starting to go in the direction that you are thinking?

Yes, that would certainly get me some of the way.

So thinking aloud for a few other bits:

• Uncertainties in flux can be encoded in a Measurements.jl array, which would then also do error propagation under the array interface operations. I can pull out the uncertainties when needed, and keep track in the metadata whether it's Gaussian or Poisson statistics.
• When it comes to fitting, there's a few other utilities that I make available (such as masking out certain bins in a fit), but I can use a wrapper for that as I already do.

Instrument responses and the ancillary files is another thing. Since they are big and tend not to change for a particular instrument (modulo updated calibrations), I have kept those seperate from spectra in the past. That does mean my Spectrum type actually keeps channel information on the spectral axis instead of energies, as you need the response files to work out what the energies are.

So at the moment I'm thinking I could do something like

using Spectra
const BinnedSpectrum = Spectrum{S, F, 2, 1} where {S, F}

struct SpectralData # a dataset type specific to SpectralFitting
    source::BinnedSpectrum
    background::BinnedSpectrum
    response::ResponseMatrix # specific to SF
    ancillary::AncillaryResponse # specific to SF
    # .. other SF specific fields
end

it's then easy for someone to take a SpectralData and get Spectra.jl spectra from it, and the other way round.

For the response files, I wrote most of an OGIP parser for SpectralFitting.jl. It doesn't have any dependencies beyond a fits reader of some kind (will be FITSFiles.jl in the future). Would this be better served migrated over to this package? If not, then Spectra.jl will probably be of limited use to X-ray observers, as the spectra would only ever contain channel information in the spectral axis and not physical units. That would also mean needing some kind of redistribution matrix and ancillary response types in Spectra.jl as well. I don't think that would be out of scope for this package, but I am happy to be told otherwise :)

I can elaborate on this if you like, or can start throwing PRs this way if you're happy for it to be included.

Could we also use the Holy Trait design you mentioned to help validate construction (i.e., checking axis sizes make sense for the relevant spectrum type)?

Yes absolutely. However, there would be some overlap at the moment. E.g., the ContiguouslyBinned and OneToOne traits are implicit in the Spectrum{S,F,M,N} -- if N == 1 and M == 2, I know it's ContiguouslyBinned and we could create a trait as such -- but If M == 2 and N != 1, I don't know if its an echelle spectrum or have some timing axis or even something else, in which case there's not enough information in the type to assign it an unambiguous trait.

One thing I found useful to do on this front is to have a Tag parameter in the type. So

abstract type AbstractSpectrum{S,F,Tag}

By default, the Tag can just be Nothing to suggest that this has no additional anything associated with it. But I can then tag a spectrum as Binned or Echelle or otherwise, and dispatch special methods for it where needed. This could lift some of the trait degeneracies. It also means a downstream user can create a tag if needed and overwrite specific functions for their needs too, which has occasionally been very useful in SF.jl. The downside is some compile invalidations if used excessively, and it's another parameteric type to carry around...

[icweaver]

Oh, this is getting exciting. To start preparing, I've just updated this PR with some very basic generalizations:

• wave --> spectral_axis, flux --> flux_axis
• W --> S in type parameters
• Discussion of energies included in addition to wavelengths in docs, and in unit handling
• Shelved some analysis (equivalent_width, line_flux, and continuum(!)) in favor of using downstream packages like SF.jl. Is there a better place for blackbody?

To help keep things organized, would merging this PR once the above looks good, and then having you open a separate PR so we can start exploring incorporating an updated AbstractSpectrum{S, F, Tag} scheme on top of this sound reasonable to you? I image that could also be a good PR to explore using the meta field more generally, if at all.

[icweaver]

For the response files, I wrote most of an OGIP parser for SpectralFitting.jl. It doesn't have any dependencies beyond a fits reader of some kind (will be FITSFiles.jl in the future). Would this be better served migrated over to this package?

TIL about OGIP, thanks for the link! I'm still reading up on how specutils handles data loaders, but it looks like moving this parsing functionality here would be in scope, yea? Would the idea be to also try and replace the SF.Spectrum with a Spectra.AbstractSpectrum that will be defined here?

Overall, since we're excising some of the old analysis bits here (e.g., continuum, equivalent_width, and friends) in favor of SF.jl and/or other downstream packages, this makes me a bit more optimistic that introducing some general purpose file I/O functionality like your OGIP parser won't lead to too much feature creep

[icweaver]

This could lift some of the trait degeneracies. It also means a downstream user can create a tag if needed and overwrite specific functions for their needs too, which has occasionally been very useful in SF.jl.

haha, yea, I had a feeling that just trying to dispatch based on array dimension may have been trying to be a bit too clever. Looking forward to exploring the tag solution more!

[fjebaker]

To help keep things organized, would merging this PR once the above looks good, and then having you open a separate PR so we can start exploring incorporating an updated AbstractSpectrum{S, F, Tag} scheme on top of this sound reasonable to you?

Sounds good! I'd like to play around with this in SF.jl a bit too and start putting a branch together using Spectra.jl to provide the spectral interface, to see what it look likes and what we have (inevitably) missed.

Would the idea be to also try and replace the SF.Spectrum with a Spectra.AbstractSpectrum that will be defined here?

Yes, that would be ideal. Then the spectra types are really owned by Spectra.jl and SpectralFitting.jl can solely be a set of utilities and abstractions for model fitting.

this makes me a bit more optimistic that introducing some general purpose file I/O functionality like your OGIP parser won't lead to too much feature creep

In that case I'll start preparing an OGIP branch for Spectra.jl as well, to move the code over from SF.jl :)

[icweaver]

Perfect. I'll give this PR another once-over then merge things in