Fitting overlapping spectral segments

[karllark]

Currently, if we have overlapping spectral segments first the spectral segments are merged and then the resulting spectrum fit. This averages/ignores the fact that often the spectral resolutions between segments can be different (different instruments, different gratings, etc.).

One idea to dealing with spectral segments in PAHFIT is to simultaneously fit the segments while fixing the common parameters. In the case of an unresolved gas emission line, this would be the line area and central wavelength (assuming the wavelength solutions are well known). The JointFitter in astropy models may provide the capability for what we need. As an example notebook showing a JointFitter result for two mocked spectral segments with 1 line. https://gist.github.com/karllark/6598fedb8129860c3c17d4ec579589b8

[karllark]

Notebook updated to use a custom Gaussian1D that as (area, mean, stddev) as parameters (area replaces amplitude). The two spectral segments have a factor of 2 different resolutions.

[karllark]

One nice thing about the JointFitter (that looks like a bad thing at first glance), is that each spectral segment requires a separate compound model. This is a good thing in that we can define the compound model for each segment to not include lines that are well outside the spectral range of the segment, hence we will not have to pay for the computation of calculating lines that are not in the spectral segment.

To account for the different spectral resolutions, we will just have setup the lines in each spectral segment to have the fwhm that corresponds with the known spectral resolution and fix the stddev/fwhm parameter for each unresolved line.

[jdtsmith]

Very cool! Thanks for looking into this. I was just sketching out a plan for doing exactly this with scipy; but this looks very straightforward in comparison. We can omit lines per segment for computational savings, but we have to be careful with Drude's, which stretch a long way to form a "pseudo-continuuum" under the PAHs themselves. We can also of course do a "final optimizer pass" in which we perform a Joint Fit with all features enabled, just to "shake things out".

One other useful idea to go with this is to have the native "spectrum" type in PAHFIT be a custom PAHFITSpectrum object, which includes a meta-data tagged set of Spectrum1D's. The tagging (like inst: "JWST/MIRI", segment: "1A", wavelength_units: UnitObj, flux_unit: UnitObj) would then key directly into the instrument pack, to tell PAHFIT what unresolved line spread resolution to use, what "cutoff" boundaries to use, etc. And of course also tell astropy.modeling what units to use for the spectral segments, which could in principle differ.

For people who still want to fit with untagged, pre-blended spectra, they are "on their own" and can specify a table (or fitting formula) for resolution. Note that the "resolution table" for SL/LL was "baked in" to PAHFIT-IDL.

I'm pulling some notes together on a YAML format for both Inst + Sci packs. The beauty of this "collection of segments" is you can pack together in one PAHFITSpectrum say Spitzer + NIRCAM + Akari (or whatever) and let it do its thing. Whether or not we turn on a very strongly worded "segment scalar" parameter is TBD. You could see it being very powerful in places, but not really PAHFIT's domain.

In terms of user-facing config, if they want to trim more from segment ends, they can copy the instrument pack and trivially customize their own. But also, at R~3000, the assumption of unresolved lines starts to be bad for many regions within galaxies. So we should probably expose a high-level doppler velocity cap, possibly per line (or class of lines), and use the convolution theorem on convolving two gaussians to predict the "broadened" gaussian line width at each wavelength/in each segment. Maybe the sci-pack is simple enough to configure this, but we could also add a top-level keyword (set).

As an aside, scipy.optimize has _much_more functionality than astropy.modeling, including trust region constraint optimizers. We may be able to "fake" that with composite fixed variables tied to each-other via a function, but we could also go "closer to the metal" and use scipy directly.

[karllark]

All great thoughts.

One comment. Getting closer to the metal via scipy.optimize is a possibility, but it would reduce the flexibility to easily swap in different optimizers (aka minmizers) and samplers (aka MCMC). One possible option would be to enable the functionality desired from scipy in astropy models.

[karllark]

How to use the JointFitter using excerpts from the notebook linked above is now in the astropy.modeling docs. https://docs.astropy.org/en/latest/modeling/jointfitter.html

[jdtsmith]

Very cool. BTW, how does the JointFitter perform for you? Is this example with two segments ~2x the time of fitting just one of them?

[karllark]

Have not checked the speed issue yet.