PAHFIT / pahfit

Model Decomposition for Near- to Mid-Infrared Spectroscopy of Astronomical Sources
https://pahfit.readthedocs.io/
18 stars 26 forks source link

Multi + glob-style instrument specification for stitched spectra fitting #229

Closed jdtsmith closed 1 year ago

jdtsmith commented 2 years ago

This supports fitting stitched spectra by allowing the user to specify an iterable of instruments instead of just a single value per input spectrum:

ins=("jwst.miri.mrs.ch1.A","jwst.miri.mrs.ch1.B")

and/or glob-style match strings:

ins="jwst.miri.mrs.ch*.[AB]"

It works by identifying all overlapping instrument segments at the given wavelengths, computing all relevant resolutions, and using the average value with lower and upper bounds set to the min and max predicted resolution. Lines which intersect only one segment have masked bounds (i.e., are fixed).

Note that, if multiple segments are provided or matched:

  1. the return value of wave_ranges() can be a list of (2-element) lists.
  2. fwhm() and resolution() return (n, 3) masked arrays (given n input wavelengths). The additional 2 values specify lower and upper bounds.

Also note that instruments() now takes an optional match argument, which can be used to test a given match:

In [399]: instruments(('jwst.nirspec.*.medium','jwst.miri.mrs.ch*.[AB]'))
Out[399]: 
['jwst.nirspec.g140.medium',
 'jwst.nirspec.g235.medium',
 'jwst.nirspec.g395.medium',
 'jwst.miri.mrs.ch1.A',
 'jwst.miri.mrs.ch1.B',
 'jwst.miri.mrs.ch2.A',
 'jwst.miri.mrs.ch2.B',
 'jwst.miri.mrs.ch3.A',
 'jwst.miri.mrs.ch3.B',
 'jwst.miri.mrs.ch4.A',
 'jwst.miri.mrs.ch4.B']
codecov-commenter commented 2 years ago

Codecov Report

Merging #229 (8e92c05) into master (661e32b) will not change coverage. The diff coverage is n/a.

@@           Coverage Diff           @@
##           master     #229   +/-   ##
=======================================
  Coverage   43.91%   43.91%           
=======================================
  Files          12       12           
  Lines         822      822           
=======================================
  Hits          361      361           
  Misses        461      461           

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

drvdputt commented 2 years ago

I have merged this into my clone of the #222 branch.

Things are a bit fiddly because some functions return different things depending on whether there's one or multiple instrument segments specified. (wave_range and fwhm being the most important ones). Currently making some extra functions to work around this. Things like a wavelength range check, and a fwhm recommendation function that gives one number + fixed or not + upper and lower bounds (the latter are already returned by the current implementation).

jdtsmith commented 2 years ago

Thanks; happy to simplify to make it easier to use. I left the default case (one segment) returning just a normal array, and multi-segments loosening the bounds. If you think it's better I can always return bounded (n,3) fwhm, etc., where the bounds will be masked for single segments.

Also, thoughts on wave_range? When there are multiple, I could "merge" them using min/max, but what if the user passes in non-continuous spectra, like jwst.miri.mrs.*.[AB]?

drvdputt commented 2 years ago

I have worked around it for now. But I feel having a consistent behavior independent of the number of segments would be more maintainable.

As for wave_range, instead of returning the range explicitly, it might be more interesting to just have some functions that tell you what the deal is for a given wavelength.

For example I created this function

def test_waves_in_any_segment(wave_micron, segments):
    """Test if each given wavelength is covered by at least one instrument

    Parameters
    ----------
    wave_micron: numpy array of wavelengths in micron (cannot be quantity)
    segments: list of segments

    Returns
    -------
    np.array of bool, True where wave_micron[i] is in the range of any of the segments."""

For non-continuous spectra, a similar check could be done to remove any wavelength points that fall in between segments. The details of that can be discussed (do we want to stop everything and raise an error, or just discard the offending wavelength points and continue).

jdtsmith commented 2 years ago

it might be more interesting to just have some functions that tell you what the deal is for a given wavelength.

I was actually thinking this same thing. I think it should also have a keyword like near which, if enabled, computes the FWHM at any of the segment ends, and uses it to estimate if you are within a fixed multiple of fwhm (2.25 is what I've been using), either direction (basically dilates the wave_ranges by 2.25 FWHMs before the test). I'm happy to adapt this for inclusion.

For non-continuous spectra, a similar check could be done to remove any wavelength points that fall in between segments. The details of that can be discussed (do we want to stop everything and raise an error, or just discard the offending wavelength points and continue).

I think don't worry about it. If segments map to a non-continuous set, so be it, and we will only tell you if you are "in" (or "near") at least one of them. You can then ask for the fwhm and be given a (possible bounded) value. Right now fwhm for multi-segments only returns values for "in". Need to evaluate this for "near" as well.

jdtsmith commented 2 years ago

Please checkout the new within_segment function to see how it works, and give it a check over. For unifying the output, what do you think about adding an optional flag as_bounded that always returns bounded output, like in the most recent commit. This way interactive use is still convenient for the common case of "plot resolution vs. wavelength" etc.

jdtsmith commented 1 year ago

I added a new check_range function that performs the "out of bounds" warning/error checks, as discussed in #230, which this addresses.

jdtsmith commented 1 year ago

Note that check_range is never called automatically, so we should do so first thing upon fit/guess/etc.

jdtsmith commented 1 year ago

Since instrument.py isn't yet used in master, I'm going to merge this and then cut the dev:api+yaml "trunk" branch.