Concept: pseudocode proposal for wavelength calibration infrastructure

[kecnry]

FOR DISCUSSION: NOT INTENDED TO MERGE

This shows the basic class-infrastructure and syntax that could be designed to support 1d and 2d wavelength calibration, both with manual line input and auto-identifying from a line-list. See also #71 for an older proposal that served as some inspiration here. See the diff for the class architecture itself, and I'll just copy the user-facing proposed syntax here as well:

# various supported syntaxes for creating a calibration object:
cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine(5500, 20), CalibrationLine(6200, 32)))
cal1d = WavelengthCalibration1D(spec1d_lamp, (CalibrationLine.by_name('H_alpha', 20, 'uphill'), CalibrationLine.by_name('Lyman_alpha', 32, 'max', {'range': 10})))
cal1d = WavelengthCalibration1D(spec1d_lamp, ((5500, 20), ('Lyman_alpha', 32)))
cal1d = WavelengthCalibration1D(spec1d_lamp, ((5500, 20), (6000, 30)), default_refinement_method='gaussian')
cal1d = WavelengthCalibration1D.autoidentify(spec1d_lamp, common_stellar_line_list)

# once we have that object, we can access the refined line wavelengths (if applicable)
# the resulting WCS, and apply it to a spectrum (either the same or separate from the spectrum used for fitting)
plt.plot(cal1d.refined_wavelengths)
targ_wcs = cal1d.wcs
spec1d_targ1_cal = cal1d(spec1d_targ1)
spec1d_targ2_cal = cal1d(spec1d_targ2)


# we could either start a 2D calibration the same way (except also passing the trace) or by passing
# the refined lines from the 1D calibration as the starting point
cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, ((5500, 20), (6000, 30)))
cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, cal1d.refined_lines)
spec2d_targ_cal = cal2d(spec2d_targ, trace_targ?)

Open questions:

• How should we handle parameter names for wavelength vs frequency (or use something general like "spectral_axis")?
• This proposal makes use of cached-properties... do we want to implement our own setter/getters and either not use dataclass or start moving away from them package-wide?
• Are we happy with the call to calibration returning spectrum1d instead of WCS or a separate WavelengthSolution object (if we stick with dataclasses or don't use caching, then we might want to have a separate object returned to avoid having to recalculate the solution when applying the same calibration to multiple specta)? See relevant lines in @eteq's original proposal
• Do we want any smoothing or model fitting in the cross-dispersion direction?

Other considerations:

• 2D refinement for non-flat traces will need to handle PART of shifted trace going out of bounds
• Extraction methods will now need to collapse according to 2D wavelength solution, if available. This might take significant work and require interpolation methods, etc.

[kecnry]

Do we want to support per-line or global refinement options, or both?

[camipacifici]

In my mind, the user would do a first pass automatically and then refine the lines that need refinement. I wonder if the refinement to the line centroid should happen earlier and not while the wavelength solution is calculated. To follow Tim's process (and if I understand this code correctly), this represents step 3 (fit model), while the centroid refinement happens in step 2 (identify).

[kecnry]

That's right, this proposal merges Tim's steps 2 and 3 into a single call, but does allow the user to inspect the results from step 2 and modify the inputs. My reasoning for having these as a single class is:

1. I don't expect the wavelength solution to be expensive at all, so having it computed and then modifying the inputs shouldn't be a computational "burden". And although both steps are covered in the same class, the user technically can defer step 3 by asking for cal1d.refined_lines and modifying those inputs as necessary before calling cal1d(spec_targ).
2. The 2D case (reidentify, "step 4") requires information about how to refine the lines at any given "row", and I think it's elegant to have similar inputs to the 1D and 2D cases, so I like having refinement information passed to the 1D case as well. Alternatively, we could split the 2D case into two steps: the first that just refines the positions across all the rows and returns those positions and the second that then fits a wavelength solution.

In other words, I would suggest we either merge steps 2 & 3 from Tim's summary below or also split step 4 into "reidentify lines across rows" and "fit 2D model".

[tepickering]

In my mind, the user would do a first pass automatically and then refine the lines that need refinement. I wonder if the refinement to the line centroid should happen earlier and not while the wavelength solution is calculated. To follow Tim's process (and if I understand this code correctly), this represents step 3 (fit model), while the centroid refinement happens in step 2 (identify).

this is correct. my step 2 has a step 2a: find lines and centroid them and step 2b: match lines to lines with known wavelengths. centroids shouldn't need to be tweaked. there will be iteration between steps 2b and 3 to find new matches as the wavelength solution is refined and then refining the solution further with more matched lines.

[kecnry]

@tepickering - I think the steps you state as "match lines to lines with known wavelengths. centroids shouldn't need to be tweaked. there will be iteration between steps 2b and 3 to find new matches as the wavelength solution is refined and then refining the solution further with more matched lines." would fall within this proposal's autoidentify classmethod. That is imagined to take a line list instead of a manual list of wavelength-pixel pairs, find the lines in the spectrum, match to the lines in the list, and create that list of wavelength-pixel pairs automatically. At that point they can still be refined when passing on to the 2D step.

So modifying the original example:

cal1d = WavelengthCalibration1D.autoidentify(spec1d_lamp, common_stellar_line_list)
cal1d.lines  # list of wavelength-pixel pairs determined automatically by identifying and matching lines between the observations and the line-list.  User can plot, inspect, modify, etc.

cal2d = WavelengthCalibration2D(spec2d_lamp, trace_lamp, cal1d.lines)
calibrated_spec2d = cal2d(spec2d_source)

Does that sound reasonable to you or is there something missing or in the wrong order for this workflow?

[kecnry]

@tepickering - just wanted to revist this... do you think the modification about would accomplish your use-case? Or if not, do you have suggestions for how to alter the API in order to support your steps better?

[kecnry]

Does 2D calibration need to support different trace for lamp and source? If so, how would that work?

[eteq]

I think the answer is "no" here in one 2D case, but see the broader comment I'm about to leave for a bit more nuance
in the other case. ("no" means "what you have here is fine")

[codecov]

Codecov Report

Merging #152 (3df4342) into main (05a174d) will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##             main     #152   +/-   ##
=======================================
  Coverage   75.46%   75.46%           
=======================================
  Files           9        9           
  Lines         693      693           
=======================================
  Hits          523      523           
  Misses        170      170           

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[tepickering]

let me break down the process at a higher level so that we're all on the same page and can have a clearer picture of how things need to fit together for various use cases:

1. Extraction: the first step is to extract a 1D calibration spectrum from a 2D image of a calibration source (lamp, sky, laser comb, etc) using a Trace. the output should be a Spectrum1D with spectral axis units of u.pixel.
2. Identify: find the peaks within the extracted calibration spectrum and match them to lines with known wavelengths. the basic way to do this is to manually match a handful of lines. more sophisticated approaches involve pattern matching against a provided line list or cross-correlating against archived wavelength solutions. i thought we had more utilities for loading line lists from NIST and elsewhere, but it looks like that work never got merged. pypeit has some well-vetted line lists that we could provide access to. likewise the IRAF line lists that some observatories (like the one i work at) still use as references.
3. Fit Model: use the line matches from the previous step to fit a wavelength solution model, identify more lines in the calibration spectrum, redo the fit possibly with revisions (e.g. increase order, change model type, etc), and iterate as needed until the desired goodness of fit is reached. the output is fundamentally a mapping of pixel->wavelength. this could simply be a WCS object, but we may want our own class to wrap WCS and simplify common operations we'll need (e.g. simple shifts).
4. Reidentify: use output from the previous step to identify calibration lines in other areas of the 2D calibration image or in other calibration images and re-fit the wavelength solution at different positions. basically, go from a single wavelength solution to a map of wavelength solutions as a function of image and image position. the simple case of a single source in a single slit does not require this step. cases that have multiple sources, slits, or orders in a single 2D image will require it. also useful for re-using wavelength solutions between different sets of data with similar configurations.
5. Apply and Tweak: use a wavelength solution to update the WCS for an extracted source's Spectrum1D. shifts can occur between the times when calibration and target images are obtained so we'll want to have ways to tweak wavelength solutions based on, for example, sky lines.

[eteq]

My own answers to @kecnry's open questions. Note I have not yet read @tepickering's comment to try to keep my perspective independent, so I might revise some of this immediately after writing it after I've read his 😉

How should we handle parameter names for wavelength vs frequency (or use something general like "spectral_axis")?

I think general - i.e. use "spectral_axis" or just "world_coord" or "physical_value", and also, very critically, they must be Quantity objects, not raw numbers. E.g., CalibrationLine(5500*u.angstrom, 20*u.pixel) (although the u.pixel is perhaps optional since that's the only valid unit for "pixel number"). I'd insist on that even if they're all the same physical type, because we absolutely cannot assume people only want to use only angstroms.

This proposal makes use of cached-properties... do we want to implement our own setter/getters and either not use dataclass or start moving away from them package-wide?

I don't see why we can't do both? I.e., dataclass by default, but cache only when needed? (a la https://stackoverflow.com/questions/51079503/dataclasses-and-property-decorator )

Are we happy with the call to calibration returning spectrum1d instead of WCS or a separate WavelengthSolution object (if we stick with dataclasses or don't use caching, then we might want to have a separate object returned to avoid having to recalculate the solution when applying the same calibration to multiple specta)? See relevant lines in @eteq's original proposal

I am very against returning a full Spectrum1D every time, because some applications (e.g. the forward modeling approach below) only need the WCS, not the spectrum itself. But I don't object to there being a few different ways of calling, where some yield the WCS, and others are a convenience that yields a Spectrum1D with the WCS already applied. And that convenience one should probably be the default that we point people to in the docs, now that you say it @kecnry. But there does need to be a way to get just the WCS out that is functionally identical (i.e., the "default method" would just call that, and then be a trivial addition of attaching the WCS to the spectrum).

Do we want any smoothing or model fitting in the cross-dispersion direction?

I don't completely get this question, but I think the answer is yes - so my comment at the end for how I think about this.

Other considerations:
2D refinement for non-flat traces will need to handle PART of shifted trace going out of bounds

I think it's probably ok to just fit part of the trace that's in-bounds. That is, the model ends up just unconstraned for pixels where the trace fell off the edge (and might fall into the wildly-extrapolated failure mode of polynomials, etc, but that's fine)

Extraction methods will now need to collapse according to 2D wavelength solution, if available. This might take significant work and require interpolation methods, etc.

Definitely this is an important and highly non-trivial thing. I think it's ok right now to basically say "all the extraction methods don't work with 2D solutions right now", and we can work on that later, as there's a whole pile of literature and experience on this topic we'd need to explore more. But for many science cases one doesn't care because you want to forward model the whole spectrum anyway, so the user doesn't need to extract.

[kecnry]

Thanks for the comments/thoughts @eteq

I think general - i.e. use "spectral_axis" or just "world_coord" or "physical_value", and also, very critically, they must be Quantity objects, not raw numbers

spectral_axis for a single value feels wrong to me (just grammatically feels like it should take an array or the index of the axis, not a single value along the axis), but I do like your other suggestions as alternates to consider. I'd be fine with accepting raw numbers (at least for pixel) and assume default units, but entirely agree that we must accept valid quantity objects.

But there does need to be a way to get just the WCS out that is functionally identical (i.e., the "default method" would just call that, and then be a trivial addition of attaching the WCS to the spectrum).

Ok, I think that is consistent with the proposal here, but correct me if I'm misunderstanding. cal1d.wcs returns the WCS itself. cal1d(spec_targ) calls cal1d.wcs internally and applies it to the passed spectrum. That WCS is cached so that additional calls to cal1d(another_spec) do not repeat anything computationally expensive.

Do we want any smoothing or model fitting in the cross-dispersion direction?
I don't completely get this question...

Let me try to explain a little better before breaking out a whiteboard. At each "row" (ok, stepping in the cross-dispersion axis away from the trace), the pixel positions of the lines are refined, giving a pixel coordinate at each row for each line. Are these allowed to be entirely independent of each other as just raw data points? Does each row get an independent 1D wavelength solution or does the whole image get a 2D model fitted/applied? Or are the wavelengths smoothed or have a separate 1D model fit to them before the wavelength solution is fitted?

I think it's probably ok to just fit part of the trace that's in-bounds. That is, the model ends up just unconstraned for pixels where the trace fell off the edge

Works for me, and since your answer to the following question is to not handle 2D solutions during extraction yet, we also won't have to handle this scenario... yet. But down the road we may need to think about how extraction will work on a 2D spectrum where some regions have an undefined solution.

Either way, as we scope out the work, we at least need to make sure we have the checks to raise an error if passing a 2D wavelength solution into extraction.

[eteq]

Now some separate high-level thoughts:

I found this API's naming a bit un-intuitive because I don't see the calibration as intrinsically 1D or 2D (because it's always 2D) but rather the resulting WCS is 1D or 2D. So I think we want the two cases to be like this:

Case 1: follow-the-spectrum-trace:

1. Extract a spectrum
2. Use the exact same trace to do an extraction on the arc
3. Produce the (x)>(wavelength) mapping for that extraction
4. Fit a 1D WCS

Case 2: 2D wavelength solution

1. Start with some sort of trace, which might be a spectral extraction trace, or might be something else, maybe even just a straight one-pixel box along the rows
2. Apply that to the arc
3. Find a 1D wavelength solution a la your WavelengthCalibration1D
4. use autoidentify to produce a pile of (x,y) -> (wavelength) mappings
5. fit a 2D WCS to the above

So the 2D solution isn't a completely separate thing, but is rather an extension of the 1D case. You'll also note that neither require an extracted spectrum, just a trace, and some use cases it's very very important that you not expect an extracted spectrum.

Relatedly, but separately, I think too much is happening in the initializer in this API. That is, it looks like you are suggesting that all of the steps above happen on initialization, and the WCS is re-computed when something changes. I think that's not ideal because frequently one wants to apply the same identify steps to a bunch of spectra, and I think it's better to create an object that operates on the spectra, rather than one that contains the spectrum. So instead I'd modify your API to do something like:

cal1d = WavelengthCalibration1D(CalibrationLine(5500*u.angstrom, 20), CalibrationLine(6200*u.angstrom, 32), onedextractorobject), wcs_model=Polynomial1D(degree=5))
newthingie = cal1d.attach_to_spectrum(spec2d_lamp, trace)  # whether this produces a new WavlenegthCalibration1D or some other object I don't really care
newthingie.plot_1d()  # plots in matplotlib, nothing too fancy
specviz.fancy_wavecal_plot_magicalness(newthingie)  # this we'd implement in matplotlib
newthingie.add_more_lines(...)
newthingie.move_lines(...)
wcs1d = newthingie.make_1dwcs() 
spec_with_wl = newthingie.apply_wavlenegth_solution(spec1d_extracted)

that's for my case 1. Case 2 would add this:

cal2d = newthingie.autoidentify(... some autoidentify parameters...)
cal2d.choose_2dWCS_model(2DSpline(...))
cal2d.make_2dwcs()  # fits the model

I'm less confident that cal2d api there makes sense, but hopefully you get the gist that the one builds from the other, they can't be treated independently.

EDIT: I added onedextractorobject above, which is a critical thing I realized I was missing while responding below.

[eteq]

(Oops, sorry, I think our posts collided there, @kecnry - I was writing that before seeing yours)

[eteq]

To respond to @kecnry last post:

spectral_axis for a single value feels wrong to me (just grammatically feels like it should take an array or the index of the axis, not a single value along the axis), but I do like your other suggestions as alternates to consider. I'd be fine with accepting raw numbers (at least for pixel) and assume default units, but entirely agree that must accept valid quantity objects.

👍 to most of this but I am very very very strongly against "assume default units". History has shown that way leads to madness and sadness.

cal1d.wcs returns the WCS itself. ...

Sorry, yes you're right, I missed the cal1d.wcs in the first read, soyes I think that's all fine (although with the caveat that I don't like all the work happening in the initializer as detailed above, but I think what you're saying still works fine with my alternative). I think maybe I'm against having a __call__ at all, and instead favor just .wcs and .apply_to_spectrum() or similar, but that's not a terribly strong opinion.

Let me try to explain a little better before breaking out a whiteboard. ...

Hopefully my post above cleared a lot of this up, but my thinking is that yes, we definitely do need to support this stuff, but it would be the same extraction parameters used to extract the arc in the first place. That is, this class just uses the extraction object to extract from the arc and that would do any smoothing, etc. The simplest version (which would probably be the default?) would be a straight line along the rows, but the user can provide really whatever they want and autoidentify just marches it up and down the spectrum while keeping track of the relevant pixel coordinates.

Works for me, ...

👍

[eteq]

Ohh, something I just realized: when I was saying autoidentify above, I really meant reidentify in @tepickering's parlance. Reading more closely I get now that's not what you meant @kecnry. Just substitute "autoidentify" -> "reidentify" in all of my messages above. I get now that you'd want to "autoidentify" as an alternative to the interactive mode I proposed above, and reidentify might or might not do that.

And also I agree with all of @tepickering's points/terminology, although I'm not convinced we want step 5, since in my mind "apply" means "put WCS where there wasn't one before", not "update an existing one with small changes"

[kecnry]

I'll need to think carefully through your proposed changes before replying to those, but I just wanted to clarify that I was thinking of having the refined_lines and wcs cached properties contain all of the "work", not the initializer. Those caches would then be cleared if any of the inputs were changed by the user and would have to do the work again on the next time they're called. This allows the user to access any of those intermediate steps (refined line locations, wcs, spectrum with wcs applied) without extra overhead.

[eteq]

Oh, gotcha, @kecnry . I agree with that in a general sense, but my quibble is: PEP8 (see note 2 under https://peps.python.org/pep-0008/#designing-for-inheritance) says properties should not be computationally expensive. Instead you should use methods for anything that's going to be expensive. I get this is a bit of an edge case given the caching, but my argument here is that we can work around that by instead designing an API where it's explicit that the computation happens when you apply the object to a specific datasets.

(But I'm not dead set on this - PEP8 also says "However, know when to be inconsistent – sometimes style guide recommendations just aren’t applicable. ")

[eteq]

Reading more carefully, I like how close #152 (comment) and #152 (comment) in practice even if framed differently. That gives me some confidence this is a good path to go down

[tepickering]

Ohh, something I just realized: when I was saying autoidentify above, I really meant reidentify in @tepickering's parlance. Reading more closely I get now that's not what you meant @kecnry. Just substitute "autoidentify" -> "reidentify" in all of my messages above. I get now that you'd want to "autoidentify" as an alternative to the interactive mode I proposed above, and reidentify might or might not do that.

in my parlance, identify is required when there is no prior knowledge of a wavelength solution and one needs to be bootstrapped. reidentify/autoidentify requires a input wavelength solution to do the line matching. identify is not always required. prior knowledge could come from a spectrograph specification, image headers, or archived solutions and be fed directly info reidentify/autoidentify.

And also I agree with all of @tepickering's points/terminology, although I'm not convinced we want step 5, since in my mind "apply" means "put WCS where there wasn't one before", not "update an existing one with small changes"

again, i was conflating a step 5a: apply calibration WCS to target spectra and step 5b: update/tweak target spectra WCS. "apply" could be "put WCS where one wasn't before", but i think it's more accurately update detector WCS to physical WCS, at least in the API we're discussing here. conceptually, step 5b is just using lines within the target spectrum (e.g. sky lines) and doing another iteration of steps 3 and 4 with them and then 5a to update the WCS. so maybe we don't need this separated out, programatically. 5a might just be target.wcs = new_wcs and 5b is just a re-iteration of previous steps that is useful for some cases, but not others. however, i wanted to put into words everything that's often needed to achieve science-grade wavelength calibration, at least from a ground-based perspective.

i do want to note that as a pypeit user, i have often wanted to use sky lines to perform step 5b or include it in pypeit's process, but its much harder and more complicated than it needs to be. so it's a bit of a personal peeve that this be more easily possible.

[kecnry]

@eteq - I'm revisiting this whole conversation and think we're all mostly (?) on the same page, but looks like a few items you brought up were lost in all the threads and not discussed (if there's anything else that you think was not resolved, .

I don't see the calibration as intrinsically 1D or 2D (because it's always 2D) but rather the resulting WCS is 1D or 2D

The 1D and 2D in this PR refer to the input spectrum. In the first case the spectrum is either a single row/trace of a 2D spectrum or an already extracted spectrum, and in the 2nd (2D) case, the calibration "walks" up and down a 2D spectral image. Because of the separate inputs and logic under-the-hood, it made sense to me to separate them, but the names themselves could easily be swapped.

I think that's not ideal because frequently one wants to apply the same identify steps to a bunch of spectra, and I think it's better to create an object that operates on the spectra, rather than one that contains the spectrum. So instead I'd modify your API to do something like:

cal1d = WavelengthCalibration1D(CalibrationLine(5500*u.angstrom, 20), CalibrationLine(6200*u.angstrom, 32), onedextractorobject), wcs_model=Polynomial1D(degree=5))
newthingie = cal1d.attach_to_spectrum(spec2d_lamp, trace)  # whether this produces a new WavlenegthCalibration1D or some other object I don't really care
newthingie.plot_1d()  # plots in matplotlib, nothing too fancy
specviz.fancy_wavecal_plot_magicalness(newthingie)  # this we'd implement in matplotlib
newthingie.add_more_lines(...)
newthingie.move_lines(...)
wcs1d = newthingie.make_1dwcs() 
spec_with_wl = newthingie.apply_wavlenegth_solution(spec1d_extracted)

To me this feels much more "interactive" and would feel a bit out of place with the current APIs of other steps in specreduce, and maybe more natural as the API within jdaviz's plugin which would wrap the calibration in specreduce. In my mind, the user interacts with specreduce objects by either changing inputs and re-running or possibly by revising the input "parameters" on an already instantiated object and then able to re-run the actual computation/calibration (in other words, no specific methods for interactively plotting/modifying individual lines, etc). But I'm also open to the idea of trying to start heading in this general direction if that is the consensus... in which case we might want to also think about what that means for trace/background/extract.

[PatrickOgle]

I don't have much to add to this discussion, which seems quite nuanced. As long as we capture the following workflow,
I think we will be in good shape:

1. Inputs: 1D or 2D calibration spectrum, line list
2. Algorithm: a) line-finding + b) Manual or auto line identification + c) Iterative fit of wavelength vs. pixel to some tolerance
3. Output: 1D/2D WCS or spectral_axis
   Apart from that, I will be happy to test-drive a PR for usability and accuracy, once there is a draft.

[rosteen]

Is there an expectation here that we will limit this to linear solutions to the dispersion, or do we need to account for non-linear spectral WCS solutions as well? I've been doing a bit of reading up on spectral WCS representations as part of this, and the non-linear case turns out to be more complicated than I expected. As far as I can tell you can't just encode, e.g., a second-degree polynomial instead of a linear solution. Do we want to use GWCS for this, which may be more flexible or have better support in this regard?

[tepickering]

yes, we will need to use GWCS for this. we will need the ability to do distortion corrections of arbitrary order as well as support discontiguous WCS (e.g. IFU slices, echelle orders, multi-slit). the gemini pipeline code, DRAGONS, uses GWCS internally and may be a good source of examples to glean from.