PAHFIT 2022

Towards an updated PAHFIT Model/Interface

Based on discussions of the April 4th PAHFIT all-hands workshop. Updated June, 28, 2022

Instrument/Science Packs

Separating the concerns of instrumental artifacts and scientific model details is important. Hence we have proposed separate “instrument” and “science” packs. Internally, information from these two types of packs is combined into a pahfit.features table, which is used both for setting up the models, inspecting model results, or altering details of the fit.

Some definitions

Feature(s)

A feature is one component of the PAHFIT model, such as a single emission line like [NeII], or a single dust feature.

Feature Group(s)

A group of features that are conceptually related. For example a group of sub-features for a composite dust feature (e.g. PAH77).

Parameter

Some value used to define a feature, e.g. wavelength for a line.

Attribute

Information about a parameter, such as its value or bounds.

Instrument Packs

The primary information an instrument pack can convey to the model is:

• The wavelength range of validity of a given spectral segment/channel/order.
• The Gaussian full-width of the line spread function at any valid rest wavelength, and its uncertainty (if any).

Other instrumental details could be implemented, which are “uninteresting” scientifically. Since instrument packs are “fixed”, and implicitly determine which instruments PAHFIT supports, and since different instruments may have different degrees of complexity in their underlying calculations of line widths, a simple table with a row for each fully distinct instrument or instrument configuration will suffice. Example:

spitzer.irs.sl1
spitzer.irs.lh
akari.irc.nirgrism
iso.sws
jwst.nirspec.R2700
jwst.miri.mrs.ch1
jwst.miri.lrs

Each row defines a polynomial expression for calculating FWHM at any applicable (rest) wavelength, and the wavelength cutoffs in microns.

Note that since the instrument pack is applied only once redshift is known, it is applied just before fitting.

Science Packs

Science packs contain all the scientific model configuration details for building a PAHFIT model. These include the list of features present, and the range of values their underlying parameters take. PAHFIT will distribute a few carefully vetted science packs, but anticipate users may make small alterations to these, or create their own. The heart of PAHFIT is its science packs. To offers a flexible format for defining packs, including grouping and identifying components for easier use, PAHFIT employs a simple YAML heirarchical plain-text format, similar to many configuration file formats. Note that internally, science packs are converted into a pahfit.features tabular format.

An outline of the syntax of a PAHFIT’s YAML science packs:

• Top level: a named set of individual features (like a single line), or groups of features (like a set of lines, continuum components, etc.)
• Elements: Each group or feature element consists of a name key, and associated attribute info. The name is just a handy means of selecting and working with a feature or a group, and is never interpreted. All top-level names in a model must be distinct, but otherwise, can be anything.
• Keys: Each feature or group contains a dict, which can have several associated keys:
  • kind (required): The kind of a feature. Current kinds, with description and the feature’s underlying model type:

    kind	Description	Default Underlying Model
    starlight_continuum	continuum from starlight emission	Blackbody
    dust_continuum	underlying dust emission continuum	Modified Blackbody
    line	individual molecular/ionic emission lines	Gaussian
    dust_feature	resolved dust emission features	Drude (possibly modified)
    attenuation	bulk attenuation, including silicate bands	named model
    absorption	individual absorption bands, such as H20	Drude

  • Other attributes, specific to each feature kind (see below).
• In addition, groups can have the following keys:
  • kind: When set on a group, kind is inherited by all features in the group. Every group or feature has a kind.
  • features: A named dictionary of features in a group, each with their own attributes.
  • attribute: A list or feature-named dict of a particular attribute for a group (e.g. wavelength). If a list, i.e. names are not included, they are generated from the group name. If multiple attribute lists are included in a single group, they must all have the same length, and the feature names (if any) are taken from the first named set encountered.

Parameters, by kind

The parameters available depend on the feature or group’s kind. All features have an implicit name (from the key atop the feature). All features have a kind, taken from a fixed set (see above), or inherited from their group.

• Continua (starlight and dust_continuum) have:
  • temperature - BB temperature, in Kelvin degrees
  • tau - the emission strength of the continuum component (f_ν units, normally unspecified)
• Resolved features (line, dust_feature) have:
  • wavelength - wavelength, in microns
  • fwhm - full-width at half-maximum, in microns.
  • power - total integrated power or intensity (normally unspecified). Unit depends on units of input spectra (Hz * [unit])

  Note: fwhm for line features is a detail of the Instrument Pack, and is not supported as input in the science pack by default, as it can vary by spectral segment or be multi-valued.

• Attenuation can have:
  • tau (normally unspecified)
  • model - a name of the attenuation model to use
  • geometry - name, currently mixed or screen, of the attenuation geometry to use. If unspecified, mixed is assumed.
• Absorption features have:
  • tau - (unitless, maximum optical depth at feature center)
  • wavelength
  • fwhm
  • geometry

A proposed but not yet implemented parameter for dust features:

• mod_asym - modified Drude asymmetry parameter (unitless, proposed)

Feature parameter format

Feature attributes can be included in one of several related formats, either individually or in a group. Here feature attributes are illustrated with the “wavelength” parameter:

1. A simple value (fixed, without bounds):

   wavelength: 9.7
       

2. A simple bounded value:

   wavelength:
       value: 10.2
       bounds: [-1%, 0.1#]
       

3. A list of simple fixed values, in a group, with implicit naming of the group features:

   cool_lines:
       wavelength: [9.7, 10.2, 15.8]
       

4. A named set of values, inside a group:

   cool_lines:
       wavelength:
           cool_line1: 9.7
           other_line: 10.2
       

5. With group-wide single parameter bounds:

   cool_group:
       bounds: [-0.1#, 1%]
       wavelength: [9.7, 10.2, 15.8]
       

6. Multiple value sets in a group, with named bounds (equivalent to the above, but useful in case more parameters than wavelength are specified for the group):

   cool_group:
       bounds:
           wavelength: [-5%, 10%]
       wavelength:
           cool_line: 9.7
           other_line: 10.2
       

7. With explicit parameter attributes in a features array:

   explicit_group:
       features:
           KrVII:
               wavelength:
                   value:  12.813
                   bounds: [-10%, 1%]
           PuVII:
               wavelength: 1.44444
       

   The latter is the only means of setting all attributes for a parameter. See below for more information on available parameter attributes.

Parameter attributes

Each feature parameter can be set to a numerical value directly (e.g. wavelength: 2.1212), or, if bounds are desired, to a dictionary with associated attributes:

• value - the (single) numerical value of the parameter.
• bounds: a 2-element list of bounds: [low, high]. Some details on bounds:
  • each of low or high can be:
    • a number, representing an absolute bound (including infinite, using YAML-specific .inf or -.inf).
    • null for no bound (same as positive/negative infinity for high/low bounds).
    • A relative offset, created by appending an #. E.g. [-0.1#, 0.1#] specifies a fixed offset of plus or minus 0.1 relative to the value.
    • A percentage offset, created by appending %. E.g. [-10%, 10%] for 10% range on either side of value, or [-5%, 0%] for a lower bound up to the value itself.
  • By default, if bounds are not included, the parameter is fixed at its given value.
  • Note that any unmentioned parameter in a given feature (except fwhm, see below) implicitly has semi-open bounds [0.0, .inf] assigned (see Default values).
  • As a special case, bounds can also be set at the group level:
    • For groups specifying only a single parameter, bounds can be a simple 2-element array as above. This sets the common bounds for that single parameter for all features in the group (typically using % or # relative offsets).
    • For groups specifying multiple parameters, bounds can also be a dict of one or more parameter: bounds pairs (e.g. wavelength: [-1%,1%]). Any parameter which is not included in the dict has no bounds set.

Proposed, but as yet unimplemented parameter attribute information:

• tied - (proposed, unimplemented) The name of another feature or group of the same kind to which this parameter is tied, or a dict of {feature: name, ratio: val_or_bounds}, naming the tied feature, where the val_or_bounds is specified as either a single value, or as a {value: val, bounds: [low, high]} dict specifying the ratio of this parameter’s value, as the numerator, to the matching tied parameter value. An omitted ratio implies a “trivial” tie, with ratio 1.0.

  Use this to tie two lines or features together with a fixed or varying ratio, typically of power or wavelength. Note this introduces a new psuedo-parameter (either fixed or independent, with bounds), and that the target tied feature can also be a group of features. See Tied constraints.

And similarly, for group-level band-summed ties:

• tied - (proposed, unimplemented) A dict of one or more name: {param: parameter, ratio: val_or_bounds} specifying another feature or group name to tie this group’s parameter to, via a sum. Typically this would be useful only for amplitudes, like power or tau. ratio is as above, except it pertains to the sum of the associated attribute values between this group and the target tied group.

For more information on tied parameter calculation for both individual features or groups, see Tied constraints.

Default values

Any parameter (except fwhm) which is required for a model, but left entirely unspecified in some feature in the science pack (e.g power for a line) is by default set to vary with lower bound of 0, and infinite upper bound ([0, .inf]). This would typically include line and feature power, continuum and attenuation τ, etc. (though these of course can be limited). As a special exception, because it will be overridden by the instrument pack, line fwhm is fixed if unmentioned.

If a parameter value is included without bounds, it is intrepreted as a fixed parameter.

Example science pack

An example of the science pack format, with several groups, and all parameter kinds:

# Starlight continuum [parameters: tau, temperature]
starlight:  # Blackbody
    kind: starlight_continuum
    temperature: 5000

# Dust continuum [parameters: tau, temperature]
dust_cont: # Modified Blackbodies
    kind: dust_continuum
    # shorthand for a group of unnamed dust_continuum features.  Naming is automatic.
    temperature: [300, 200, 135, 90, 65, 50, 40, 35]

# Unresolved Lines [parameters: wavelength, power]
H2:  # simple line group, with common bounds for all lines
    kind: line # Gaussians
    bounds: [-0.1%, 0%] # bounds can be set for the entire group, here leftward only
    # Shorthand for a group of named lines.  
    wavelength: 
        H2_S(7): 5.5115 
        H2_S(6): 6.1088

# An explicit group of unresolved lines
ionic_lines: 
    kind: line 
    # Full form of lines [parameters: wavelength, power]
    features:
        NeII:
            wavelength:
                value:  12.813
                bounds: [-10%, 1%] # fractional bounds for wavelength
        NeIII:
            wavelength: 15.555
            power:
                tied: # experimental/unimplemented
                    feature: NeII         
                    ratio:
                        value: 4.
                        bounds: [-1.5#, 0.75#]

# A group of dust (sub-)features [parameters: wavelength, power, fwhm, mod_asym]
PAH77:
    kind: dust_feature # Drude profiles
    tied:  # group-level parameter tie for the name parameter
        PAH113:
            param: power
            ratio:
                value: 3.
                bounds: [2.5,4.5]
    features:
        PAH77a:
            wavelength: 7.42
            fwhm: 0.935
        PAH77b:
            wavelength: 7.60
            fwhm: 0.3344
        PAH77c:
            wavelength: 7.85
            fwhm:
                value: 0.416
                bounds: [-5%, 15%]
            mod_asym: 0.25 # unimplemented: indicates a modified/enhanced-asymmetry Drude profile.

PAH113:
    kind: dust_feature
    wavelength:
        value: 11.322
        bounds: [-2%, 5%]
    power: 
        tied: NeII # single parameter tie, implicit ratio of 1.0

# Attenuation [parameters: tau, model]
silicate:
    kind: attenuation
    geometry: mixed
    model: kvt

# Absorption [parameters: tau, wavelength, fwhm, geometry]
absorption_features:
    kind: absorption
    features:
        H2O_3:
            wavelength: 3.05
            fwhm:
                value: .406
                bounds: [-1%, 10%]

Example Usage

Some example usage for preparing and interacting with a model in this UI.

Use Spectrum1D for inputs:

PAHFIT can operate on raw Spectrum1D objects (or iterables of such objects), as long as each is tagged with meta information identifying the instrument, e.g. instrument.jwst.nirspec.R1000. Such instrument names key directly into pahfit.instruments, which provide relevant instrument information to include in the model, such as wavelength cutoff and fwhm(lambda).

A single model can be created and updated with multiple Spectrum1D objects, from different instruments (see advanced examples, below).

An example of creating a simple random Spectrum1D with tagged instrument:

from specutils import Spectrum1D
from astropy.nddata import StdDevUncertainty
import astropy.units as u
import pahfit

# Get yourself a Spectrum1D object or two..
N=100
redshift = 0.05
flux = np.random.rand(N)*10 * u.mJy
flux_unc = StdDevUncertainty(flux/(5+10*np.random.rand(100)))
wav = np.linspace(5,28,N)*u.micron
spec = Spectrum1D(spectral_axis=wav, flux=flux, uncertainty=flux_unc,
                  meta={"instrument": "jwst.miri.MRS.segment1A"))
spec2 = Spectrum1D(spectral_axis=wav2, flux=flux2, uncertainty=flux_unc2,
                   meta={"instrument": "jwst.miri.MRS.segment1B"))

Spectrum1D redshifts are ignored.

Building a model:

A model is created by combining one or more instrument packs and a science pack. The instrument pack(s) are inferred directly from the input data at the time of the fit. Note that redshift, which is also input on model creation, is important because a given line’s width depends on which instrument segment and observed wavelength it appears at.

Input data: a Spectrum1D (or list of Spectrum1D’s), whose meta dict indicates what instrument/mode the spectrum pertains to.

# Pick desired sci-pack. Creates an underlying params pahfit.features.Features table.
scipack  = 'extended_exgal.yaml'

# Create a model, from the instrument-tagged spectra and science packs.
# Optionally pass instruments (if not tagged), and a redshift if needed.
mymodel  = pahfit.model(spec, scipack, [redshift= , instruments=('jwst.mrs.ch4')])

Estimating parameter starting guess

Once a model is created, its starting point can be estimated (guessed). An advantage of separating this step from the fit itself is it’s easy to override, e.g. for cubes (where “borrowing” an average starting position will speed convergence).

mymodel.guess(spec)
mymodel.plot(spec) # Plot the "guess", with the data

Performing the fit

Optionally, provide updated model result plotting during the fit.

mymodel.fit(spec, plot_progress=True)

Inspecting/altering parameter information

As noted above, creating a PAHFITModel internally creates an astropy.features object, which extends astropy.Table. After any fit, it also includes the model-estimated covariance information (if any), and can combine arbitrary combinations of features (with uncertainties, of the covariance is available). Since it is also an astropy.Table, all the group_by etc. functionality of that class can be employed. It also offers its own simple convenience selection functionality for group and kind.

# Get feature Table
fts = om.features # all features as a table
om.get_features('NeIII') # By name (same as fts[fts['name'] == 'NeIII'])
om.get_features(group='H2') # same as fts[fts['group'] == H2']

pah77 = om.get_features(group='PAH77').power # integrated intensity and estimated uncertainty, from covariance
total_lines = om.get_features(kind='line').power 

Plotting models

Any model can plot itself, with or without underlying data. Each model knows the (possibly disjoint) range of wavelengths to which it applies, and can sample the model with samples points over that range (for smoother plots)_.

The plot will contain a simple widget for toggling on/off by kind and group.

nsamp = 2000
mymodel.plot(samples = nsamp)  # plot, without data!
mymodel.plot(spec, samples = None) # Plot with data at same wavelength positions

mymodel.plot(include = ("continuum","lines"))

Saving/restoring models from disk

Models can be saved in any format astropy.Table supports, though the covariance can only be saved to formats that support metadata (like asdf and ecsv).

# Save to an ASDF (aka FITS TNG) file
mymodel.save(file='my_PAHFIT_model.asdf')

# Or save the table representation to a plain text file
mymodel.save(file='my_PAHFIT_model.ecsv')

# Restore a model from disk and plot it
om = pahfit.model.restore('another_older_model.ipac')
om.plot()

Inspecting model results

Separating out components of the model spectrum is easy. This can be done at any time, before or after the model has initial starting feature parameters estimated, or is fit to data.

Note that all returned model spectra are masked arrays, over relevant parameter bounds only. These model results could be Spectrum1D objects (or a lightly overloaded version that uses microns for the wavelength by default!).

Same samples keyword as for plotting.

# Model "spectrum"
model_fit = om.model # full model spectrum
pahh11 = om.get_model('PAH11.0', samples = nsamp) # Model of individual feature, by name

# Grouped features are specified in the sci-pack.  
pah7 = om.get_fit(group='PAH77') # w/ default sampling
all_lines = om.get_fit(kind='line')
h2lines = om.get_fit(group='H2') # All H2 lines (by group)
neIII = om.get_fit('NeIII', samples=nsamp) # Just one line, by name

# Refer to individual components by name
cont = om.get_fit(spec, kind='continuum') # at same sampling as the spectrum
starlight = om.get_fit('starlight') # default sampling

# Directly use the lower-level plot engine for plotting any arbitrary model combinations
pahfit.plot(spec, cont+starlight+neIII)

More advanced examples

Fitting spectra from multiple instruments at the same time (may include multiple segments/channels of the same instrument).

fts= mymodel.features
print('Model table: ', fts.info())
fts.show_in_notebook()

# Advanced: do stuff to the underlying astropy.Table if you are
# feeling cheeky
long_lines = fts['kind'] == 'line' & fts['lam_r'] > 15.
fts['amp'][long_lines] *= 1.1

# Create a model to simultaneous fit with multiple different instruments
my_combo_model = pahfit.model((specNS, specMRS), scipack)