New internal API: agnostic functional model setup / results retrieval

[jdtsmith]

To allow for different model-fitting cores to be optionally employed within PAHFIT, we need to specify an internal API: Fitter. Any compliant fitter can then be utilized to update the fit, and internal fitting details (including: internal units employed, the fitting framework setup and teardown, etc.) are all hidden behind the Fitter API. Many aspect of PAHFIT will then sit outside of Fitter, and interact as needed with it via the API. These include:

• Features reading and manipulation (already in hand).
• Plotting models
• "Guessing" Features parameters; Fitter-agnostic.
• Evaluating specific component models (Blackbody, Drude, etc.) and specified wavelengths.
• tabulate (mostly using the above with a nicer interface).
• Saving & re-reading models from disk (provided by Model + Features)
• Updating the Features table by performing the fit, translating as needed to the working features units.
• Instrument pack updates to the model parameters. Here the Fitter can specific address a new Instrument class at the last minute to update e.g. FWHM (including, down the road, for resolved features).

Some thoughts/details we'll need to settle on:

• When and how the Instrument pack is married to the Features input model information. This mostly (for now) involves updating the FWHM in the Features table, and disabling some features.
• The logic governing user-defined vs. instrument-pack-derived input parameters should happen outside Fitter.
• Perhaps Fitter should not know anything about the Features table, but accept simpler parameter inputs/bounds/etc.
• Perhaps Fitter should similarly return simplified outputs:
  • Best-fit parameter values for the varying params
  • Covariance matrix among varying parameters.
• It is the job of Model to adapt Features, marrying in the Instrument information, to Fitter inputs, and to translate Fitter outputs back into the Features table.

[drvdputt]

Just a reminder that doing this will help us to finally remove base.py. The old features-to-model converter can be moved to a Fitter subclass, and the plotting code can be moved into Model.

[drvdputt]

To start, we can make a list of functions / quantities we need access to. We can then think about what the ideal function/class structure should be.

Here’s a start:

• Choose a Fitter backend (e.g. choose Fitter subclass)
• Pass current numbers from features table to set up model. Ideally with no conversion needed (e.g. Gaussian - - components take wavelength, FWHM, and power, instead of mean, stddev and amplitude).
• The fitting itself. Spectrum data will be passed.
• Getting the fit results in a format that can easily be written to the features table, by just providing the feature name. The current “parse_astropy_result” will be hidden behind the API.
  • Getting uncertainties if available, and fit status / quality.
• Evaluating the model, or single components. Needed for plotting and tabulating. Either by passing feature name, or passing a row/subset of the feature table.

Things that don’t go behind the API

• Plotting
• Saving
• Guessing (no model is used, numbers filled directly in to the feature table.)
• Unit conversions (everything needs to be converted first. API can assume internal units)
• Redshift handling (data is shifted before handing to API)
• Instrument model (But API will need ‘fixed’ options to enable our fixed FWHM of course)

[jdtsmith]

This is the highest priority action. The one thing I'm questioning is whether the fitter cores (Fitter subclasses, e.g. APFitter and SPFitter for astropy and scipy fitters) should "know about" and directly interact with Features tables.

The advantage of having them do so is it's less work: you can pass a Features table directly in to e.g. fit, and it will do the right thing and update it in place (and attach the covariance matrix if available).

The disadvantage is now you are quite strongly coupled, so any future changes to Features would need to preserve backward compatibility.

One thing is clear is that Fitter's should not need to know anything about instruments.

Maybe it helps to think about what kind of thing Fitter.SomeSubClass.fit would accept and return, if not a Features table? The IDL MPFIT equivalent is a list of structures with parameter name, starting value, boundedness, and bounds; the so called parinfo structure. Somebody actually ported this to python; here's their parinfo details.

Another detail, to save room for future improvements:

Soon we'd like to fit collections of Spectrum1D objects. My fitter's approach is to concat these into one big "franken-spectrum" with non-monotonic wavelengths, and then specify segment constraints. So e.g. a narrow line will apply to one segment, and a wider (lower-resolution) copy will apply to another. Once we dilate the width of Drude's too, this will happen commonly.

It would be easy enough to just hand fit a list of Spectrum1D's (suitably tagged with instrument name), and let it set up its fit internally for that. I don't know if AstroPy will get there, but we should save room for it to do so via the JointFitter. In some ways this argues for letting each Fitter sub-class know about and read/interact with Features tables directly.

Thoughts?

[drvdputt]

So, I started writing down my thoughts, and it became a bit long. See a little summary at the end.

This is the highest priority action. The one thing I'm questioning is whether the fitter cores (Fitter subclasses, e.g. APFitter and SPFitter for astropy and scipy fitters) should "know about" and directly interact with Features tables.

The advantage of having them do so is it's less work: you can pass a Features table directly in to e.g. fit, and it will do the right thing and update it in place (and attach the covariance matrix if available).

The disadvantage is now you are quite strongly coupled, so any future changes to Features would need to preserve backward compatibility.

I think there will aways be a coupling here. The "fitters" need to build up a model, based on the info in the features table. When more different shapes for a feature are supported (e.g. allow for asymmetry), then the complexity of all "fitter" subclasses will increase.

Changes to Features, will usually be to support a new kind of fit component or parameter. But I do agree that in the future, the column names or the way the bounds/fixed are stored might change. Some backward compatibility would then be nice. I've made some suggestions below.

One thing is clear is that Fitter's should not need to know anything about instruments.

This is reasonable. To put it into strong words: The "Model" class should combine the information in Features and instrument to tell "Fitter" what to do.

Maybe it helps to think about what kind of thing Fitter.SomeSubClass.fit would accept and return, if not a Features table? The IDL MPFIT equivalent is a list of structures with parameter name, starting value, boundedness, and bounds; the so called parinfo structure. Somebody actually ported this to python; here's their parinfo details.

I feel like there would be too much redundancy, since the info you mentioned above is already contained in the Features table. In fact, there is already param_info in the "old" python pahfit code, which I was trying to phase out. I found this type of object difficult to maintain, since it needed to be unpacked again when the Astropy model is actually constructed, and any conceptual changes would need to be propagated from Features -> param_info -> astropy or other model construction.

If we really want backwards compatibility, perhaps a series of functions, instead of a list of dicts, could be provided. E.g. "add_gaussian_component(fwhm, wavelength, ...)". The param_info dicts would likely (and annoyingly) be slightly different per type of component. A collection of functions would be more transparent, since the available functions define what types of components are possible, and the function signatures enforce what input is needed.

I once coded up an astropy model constructor which built the model directly from Features, and it turned out to be quite straightforward. Just a loop over the rows, and depending on the feature type, components of the right class were constructed and added to the astropy composite model. The initial, range, and "fixed", values were set according to what was in the Features table + any overrides. With the "add component" functions, this loop would just make the appropriate function calls. This would also deduplicate the loop over the Features.

As for the return values of the fitters, I have one practical consideration which might provide some insight: plotting individual components and the tabulate function. In both of these, I am still dealing with Astropy model objects. Since the underlying objects (Astropy model vs your scipy model class) will be different in terms of how they are evaluated (at a list of wavelengths) or indexed (acessing individual components), I feel like the Fitter API should provide a uniform way to evaluate the total model or its subcomponents.

Besides that, there should be a uniform way to extract the fit results and fill them into Features. Either by returning the parameters for individual components through a function, in a common format, or by simply having some sort of "export to Features" option.

Another detail, to save room for future improvements:

Soon we'd like to fit collections of Spectrum1D objects. My fitter's approach is to concat these into one big "franken-spectrum" with non-monotonic wavelengths, and then specify segment constraints. So e.g. a narrow line will apply to one segment, and a wider (lower-resolution) copy will apply to another. Once we dilate the width of Drude's too, this will happen commonly.

It would be easy enough to just hand fit a list of Spectrum1D's (suitably tagged with instrument name), and let it set up its fit internally for that. I don't know if AstroPy will get there, but we should save room for it to do so via the JointFitter. In some ways this argues for letting each Fitter sub-class know about and read/interact with Features tables directly.

Yes, dealing with multi-segment input data is a complexity which I expect to depend entirely on the fitter (whether JointFitter is used for the Astropy backend, or something I could home-cook). We could pass the list of Spectrum1D, but this has the problem that we are tagging the Spectrum1D with instrument names. This has no real practical implications, but conceptually, it does mean that some instrument info will leak into the Fitter.

Summarizing:

We can avoid duplicating any loops over the Features table by keeping the features table out of the "Fitter" concept. But at the same time, I don't want something like param_info, for its lack of explicitness (it doesn't advertise which components are possible, and which parameters they require). The alternative suggestion, is that he API has one function for each component that can be added to the model, which advertises the data needed for the setup of the respective component.

The Fitter workflow would then be

• Input: functions that set up components
• do the fit (pass spectra in a common format)
• Output: functions that return the results in a common format (also per component type?)

The Model implementation should drive the above functions, according to the contents of Features, the instrument, and any other options.

[jdtsmith]

Thanks for these good thoughts.

I agree, maintaining two roughly-similar containers for parameter information seems like an unnecessary hassle. There may still be reasons to do so, to make a "rawer" intermediate representation of the parameter set (IRPS).

RE: tabulate: absolutely, the Fitter API will specify a set of required methods per component type (MBB, BB, Drude, Gaussian, ...), which just take the "basic model component parameters" like central wavelength, FWHM, etc. in fixed, canonical units (#259, which we should agree on), and do the calculation then in those units however they want. Otherwise each fitter has to worry about/convert the units in a Features table itself.

So standardizing the units is one argument for a raw IRPS. Different fitters may do this differently (and verifying they produce same/similar results is low-hanging test fruit). Then Model can be the one to tabulate up different subsets of the model ("all the lines", "continuum-only", "now with extinction/absorption", "just this set of PAH bands"). Fitter does not get any say whatsoever on the meaning of these features and their parameters: its job is to evaluate them faithfully, optimize them based on the input spectrum/a, and estimate their uncertainty and covariance.

On multi-segment fitting, this is complex. You are right, this does constitute an "instrument leakage" into Fitter. Right now Model jiggers the Gaussian fwhm directly in the Features table, right? So either Fitter will need some way to "fork" a single model component into multiple power-tied components (main line, ghost line #1, ghost line #2, ....), or this will need to be done from the outside, i.e. the IRPS will need to include feature ties, segment maps, etc.

This was one specific advantage of a more verbose "intermediate representation" for parameter info: Model can set parameters up to include segment maps (these parameters here apply only to this range of the spectrum), etc. Some of this I guess depends on what ambitions we have for the Astropy Modeling based fitter: multi-segment JointFitter? Parameter ties?

To help with all this, we can provide some tooling over on the Instrument side (maybe this deserves a new class) to ease the situation with multi-segment fits, so fitters don't have to do this all themselves. For now your astropy one will probably just raise and error if passed multiple segments anyway ;).

I'll work up a PR Fitter abc in a separate branch and we can discuss specifics there.

[drvdputt]

I have made some progress on this. I have short term goal to get the asymmetric features ready, as I really need those profiles for my MIRI MRS data. However, I did not think it was a good idea to implement those using the old param_info framework, since it would have needed a lot of boilerplate to add a new type of feature to the current Astropy model construction.

Instead, I implemented a draft for the AstropyFitter subclass of our hypothetical Fitter API. I could then straightforwardly rewrite the model builder, as a loop over the features table in Model. It seems to work with my existing notebooks, and with the test cases for Model I wrote in the past I think I squashed most bugs already.

If you're interested in my draft, you can take a look at this branch. https://github.com/drvdputt/pahfit/tree/local-changes%2Bmy-model-builder. See mainly fitter.py and model.py, as well as the fact that base.py was removed.

Caveat: my branch also includes other changes, such as a re-write of the default plot, and fixes for some recent known issues. So we won't be able to merge it all in one go. But if you like this version of the API, I would be happy to copy the right parts into a fresh branch.

A couple of pointers as to how it works:

1. Setting up the components: Model loops over its Features, and passes the appropriate numbers using calls such as AstropyFitter.register_starlight('name', tau, temperature), and AstropyFitter.register_line('name', power, wavelength, fwhm). Each of the arguments is a 3-tuple from Features, and the bounds and fixed settings for the Astropy model are set internally based on the values in that tuple. At the end of this loop, AstropyFitter.finalize_model() combines the registered components, and is then ready for fitting.
2. For the fitting, the arguments passed are Astropy.Fitter.fit(x, y, unc), where the wavelength x, the flux y and the uncertainty unc have to be in internal units. The Fitter class only deals with the numbers, without thinking about redshift or instruments. That responsibility rests with Model, where x, y, and unc are adjusted for redshift before passing them.

3. To retrieve the parameters, I went with something simple: AstropyFitter.get_result(name) will just return a dictionary with the fitted parameter names and values using standard PAHFIT nomenclature (column names of the Features table) and units.

   Things that are ignored for now: multi-segment fitting, uncertainties.

[jdtsmith]

Sounds like a step in the right direction. Do you want to spin up a simple branch that contains just an empty Fitter class (potentially derived from abc.ABC) and we can hack on it there?

One issue with AsymDrude's I thought of a while back: their power and amplitude representations are non-trivially related (unlike Gaussians and Drudes).

[jdtsmith]

@drvdputt where does this stand now in the context of #280?

[drvdputt]

This is now the necessary next step for #280, as the current state of the dev branch was broken by pulling in https://github.com/PAHFIT/pahfit/pull/283. Starting work on this right now. Will create a new branch that supersedes #270.

[jdtsmith]

Closing in favor of #289.