Fitting From a Configuration File#

galspec.pipeline fits one spectrum, or a whole sample with one model, from a TOML file instead of a Python script:

./example/fit_iras07546_config.sh            # the published IRAS 07546 fit
./example/fit_iras07546_config.sh --dry-run  # what would be fitted, no sampling

The pipeline is a consumer of the models: it builds the same JaxCompoundModel a script would assemble, and runs the same JNesty_Fit. It adds no model code and changes none, so a fit produced from a file and a fit produced by hand of the same model and data are the same fit. galspec/__init__.py does not import it; nothing changes for code that does not use it.

The running example in this chapter is example/fit_iras07546_config.toml: the seventeen-component, 27-parameter model of example/fit_iras07546_jax.py, written as a file and checked against the script parameter by script parameter (develop/dev_config/test_iras07546_config.py). Read it next to this page.

A file can be much shorter than that one. Nested sampling draws from the prior, so neither the starting values nor the priors have to be written down – the IRAS file below writes no starting value at all, only what the model is (line centres, a fixed redshift, the ties) – and the narrow-line profile can travel inside the spectrum’s own FITS file. A complete configuration for one emission line is:

[data]
file = "spectrum.fits"
narrow_template_key = "NARROW"      # optional; a Gaussian if omitted

[model]
[[model.component]]
type = "Line_template"
wavec = "@Halpha"

That file has two free parameters, the profile’s amplitude and its velocity, and a continuum component would be the natural next line. --dry-run prints the prior every omitted parameter ended up with, and the same list goes into run_report.json.

Starting from a template#

The file format has a lot of optional keys, so the command will write you a commented starter:

$ python -m galspec.pipeline --template config.toml
wrote config.toml
  edit [data] file, then
    python -m galspec.pipeline config.toml --dry-run    # what would be fitted
    python -m galspec.pipeline config.toml              # the fit

--template alone writes config.toml in the current directory; give it a path to write elsewhere, - to print to standard output, and --force to overwrite a file that is already there (without it an existing file is an error, so the command cannot eat a configuration you have written). It exits without touching any data.

The starter is a working configuration rather than a sheet of placeholders: in a checkout it points at the two-line mock spectrum in example/data, so the first --dry-run reports one source and eight free parameters before you edit anything, and the fit it describes takes about forty seconds on a CPU. From there the edits are: [data] file (or glob/catalog for a sample), the components in [model], and the priors. Every section it writes is explained below.

Why a file#

A sample is fitted with the same model, the same priors and the same sampler settings for every source, and only the data and a few flux scales differ. Writing that down once – instead of copy-pasting a script per object – is what makes a sample fit reproducible: the file is hashed into the output, and the run report records every setting that produced each result.

Requirements on the data#

One FITS file per source. The pipeline reads one file at a time; a “sample” is a list of files, not a list of rows inside one file.
Three arrays: wavelength, flux and error. All three are mandatory. A missing error array is an error, not a guess: an uncertainty invented from std(flux) is not a measurement, so the pipeline refuses to make one.

The arrays are found either as columns of a binary table or as image HDUs whose EXTNAME matches the key:
```
[data]
wave_key = "WAVELENGTH"     # default; matched case-insensitively
flux_key = "FLUX"
error_key = "ERROR"
```
A column name is matched first; an image HDU is looked up by extension name. Use wave_hdu/flux_hdu/error_hdu when a file has more than one candidate (they take an HDU index or an extension name).
Wavelength is converted to Angstrom, flux is not converted. The unit is taken from TUNITn (a table column) or BUNIT/CUNIT1 (an image), or from [data] wave_unit/flux_unit. A wavelength in nm or um is converted; an axis that is not a length (velocity, frequency, pixels) is rejected, because the models are evaluated on a wavelength grid. The flux unit is recorded in the output header and the report and otherwise left alone: no f_lambda/f_nu conversion is attempted.

The header supplies the two corrections. A spectrum is fitted in the rest frame, with the Galaxy’s reddening removed, and both numbers can be read from the file rather than written into the configuration. One keyword each:

quantity	header keyword (`[data]`)	what the pipeline does
redshift	`z_key`, default `"ZRED"`	divides the wavelength by `1 + z`; the model’s line centres are rest wavelengths (`"@Halpha"` or `6562.819`)
Galactic reddening	`e_bv_key`, default `"EBV"`	de-reddens the flux and the error with `extinction.ccm89` at `extinction_rv`, on the observed grid before the wavelength is shifted. The keyword’s name says how to read it: a colour excess is taken as is, an `A_V` is divided by `R_V` first, so `e_bv_key = "A_V"` is a complete statement.

A single keyword, not a list of spellings, because a file either has the keyword or it does not: a list turns “this file has no redshift in it” into a quiet fallback onto whichever of five keys happened to be there for an unrelated reason. If the file stores it under another name, name it once (z_key = "MYZ"). extinction_rv is a configuration value only – no header is trusted for it, since a survey that stores A_V leans on R_V and a wrong R_V would propagate into every flux without leaving a trace. It defaults to 3.1, the Milky Way value.

The precedence is: the catalog’s z column, the header, z_default. A value that is not in the file is assumed to be zero and the run warns: the report records z_source and e_bv_source per source, so a correction that was assumed is visible next to one that was applied. Writing the key in the configuration (z_default = 0.0, or extinction = false for a synthetic spectrum) counts as stating it, and the run is then silent.

Note

The pipeline does not mask sky lines and does not synthesize errors: if your spectra need that, prepare a FITS file with the three arrays – that file is the input, and it keeps the processing visible and repeatable. See Preparing a file for the IRAS spectrum, which needs both.

Preparing a file#

Two files in example/data/ show why the three arrays are a requirement rather than a convenience, and what “prepare it” means in practice. A template is the one extra thing a file may carry, and it is deliberately optional: a spectrum without one is fitted with Gaussian lines.

spec-0651-52141-0434.fits is a BOSS coadd: it has flux, loglam and ivar, so there is a wavelength grid and a noise estimate, but no error array and no length unit. One short script turns it into pipeline input – a table with WAVELENGTH = 10**loglam, FLUX = flux and ERROR = 1/sqrt(ivar), with the masked pixels dropped – and the choices it makes (which pixels to keep, what to do with zero ivar) stay in the file instead of being buried in a fitting script.

ir07546sed.fits is a single image with no error and an observed-frame logarithmic axis; the published fit (example.iras07546_model) dereddens it with extinction.ccm89, divides by 1 + z and derives the error from a running median. Only the corrections can move into the pipeline, because the header can carry them; the error cannot, because it is not a number in the file. develop/dev_config/convert_iras07546.py writes the input therefore:

python develop/dev_config/convert_iras07546.py     # -> example/data/iras07546_pipeline.fits

It stores the observed-frame WAVELENGTH, FLUX and ERROR columns, puts REDSHIFT and EBV (and EBVRV) in the header, adds the measured narrow-line profile as a fourth extension (NARROW: velocity in km/s, normalised flux – what [data] narrow_template_key reads, so the profile travels with the spectrum instead of beside it), and then checks itself: it applies the pipeline’s own de-reddening and de-shifting to what it wrote and asserts that the result is iras07546_model.load_spectrum() to 1e-12 relative. Wavelengths round trip exactly (0.0), flux and error at 1e-16. That check is what makes the configuration a reproduction of the script rather than an approximation of it. python develop/dev_config/convert_iras07546.py --check-only runs the check without writing.

A configuration is then a file of declarations, and the section-by-section walkthrough below reads example/fit_iras07546_config.toml in order.

The run section#

[run]
output_dir = "out/iras07546"
workers = 1

[run] says where the outputs go and how many processes do the work:

key	default	meaning
`output_dir`	`"out"`	resolved relative to the configuration file, not the working directory.
`workers`	`1`	`1` fits in this process; `>1` spawns that many processes, each fitting its own sources.
`force`	`false`	refit even when a source’s outputs exist.
`progress`	unset	the sampler’s progress bar: unset means “on a terminal, with `workers = 1`”; `true` forces it, `false` never draws it.

The data section#

[data]
file = "data/iras07546_pipeline.fits"
narrow_template_key = "NARROW"     # the measured line profile, inside the file

z_key = "REDSHIFT"                 # one keyword; the default is "ZRED"
e_bv_key = "EBV"                   # ... and "EBV"
extinction = true
extinction_rv = 3.1                # R_V of the law; the default, config-only

wave_min = 4200.0
wave_max = 6800.0
continuum_percentile = 15.0
min_pixels = 50000

key	default	meaning
`file`, `files`, `glob`, `catalog`	–	exactly one of them: one spectrum, a list of spectra, a shell pattern (relative to the config file, or `**` for recursion), or a CSV/TSV with one row per source.
`catalog_id_column`	`id`	catalog column with the source identifier (falls back to the file stem).
`catalog_file_column`	`file`	catalog column with the spectrum path, relative to the catalog.
`catalog_z_column`	–	catalog column with the redshift, if it has one.
`wave_key`, `flux_key`, `error_key`	`WAVELENGTH`, `FLUX`, `ERROR`	array names, matched case-insensitively.
`wave_hdu`, `flux_hdu`, `error_hdu`	–	restrict a lookup to one HDU (index or extension name).
`wave_unit`, `flux_unit`	from the header	override the header’s unit.
`narrow_template_key`	–	name (or index) of an HDU holding the narrow-line profile as a two-column table: velocity in km/s, then normalised flux. Each spectrum then carries its own profile, which is what a sample usually needs – the lines of different galaxies are not the same shape. A component’s `template_file` overrides it; components whose model takes no profile ignore it, and every other array lookup is unaffected. The columns are read by position, so their names do not matter; a velocity column with a unit is converted to km/s, and one without is taken to be km/s already.
`z_key`	`"ZRED"`	the header keyword holding the object’s redshift. One keyword, not a list: name the one your files use (`"MYZ"`); a file without it falls back to `z_default` and says so.
`z_default`	`0.0`	redshift when neither the catalog nor the header says. Declaring it explicitly is what silences the “assumed” warning.
`e_bv_key`	`"EBV"`	the header keyword holding the Galactic reddening. A colour-excess name (`EBV`, `E(B-V)`) is read as is; an extinction name (`A_V`) is divided by `extinction_rv`.
`extinction`	`true`	de-redden with that keyword; `false` declares that this spectrum needs no correction (and suppresses the warning).
`extinction_rv`	`3.1`	`R_V` of the reddening law. Configuration only – no header is consulted for it, because a wrong `R_V` would scale every flux silently.
`wave_min`, `wave_max`	–	rest-frame range to fit.
`mask`	–	list of `[low, high]` rest-frame windows to drop, e.g. residuals of a telluric band.
`min_pixels`	`50`	refuse a source with fewer pixels than this after the cuts. The IRAS window keeps 56956; the value in the file fails loudly if a wrong window is used by accident.
`continuum_percentile`	`50`	percentile behind the `continuum_level` expression token. The example script uses `np.percentile(flux, 15.0)`, so the file says `15`.

A catalog looks like:

id,file,z
J0925,spectra/J0925.fits,0.0312
J0931,spectra/J0931.fits,0.0455

The fit section#

[fit]
nlive = 400
dlogz = 0.01
max_iterations = 200000
share_compilation = true
device = "auto"

[fit] passes sampler settings to JNesty_Fit; the values above are the ones example/fit_iras07546_jax.py uses (NLIVE = 400, DLOGZ = 0.01).

key	default	meaning
`nlive`	`500`	live points. The main cost knob; 500 to 1000 is usual for a dozen parameters, 400 is what the IRAS model was published with.
`dlogz`	`0.01`	stop when the remaining evidence is below this (nats).
`max_iterations`	`200000`	safety cap.
`bound`	`"multi"`	`none`/`single`/`multi` bounding.
`device`	`"auto"`	`auto`/`cpu`/`gpu`.
`queue_size`, `batch_size`, `rwalk_K`, `rwalk_step_scale`	automatic	forwarded to JNesty; leave unset unless you are tuning.
`max_ellipsoids`, `memory_frac`, `min_eff`, `min_ncall`	`20`, `0.9`, `10`, –	forwarded to JNesty.
`share_compilation`	`true`	compile the likelihood once for the whole sample (see Compiling once for a sample).
`compilation_cache_dir`	–	a directory for JAX’s on-disk compilation cache.

The output section#

[output]
figures = "pdf"                    # one <id>_fit.pdf; "png" or "both" also work
fit_plot = true
corner = true
corner_params = [
    "Continuum.amplitude", "Continuum.alpha", "FeII.amplitude",
    "b_halpha.amp_c", "b_halpha.amp_w0", "b_halpha.dv_c", "b_halpha.sigma_c",
    # ... five more
]

histogram = true
hist_params = ["Continuum.alpha", "FeII.amplitude", "b_halpha.amp_c", # ...]

[output] controls the two files. figures (default "pdf") chooses the container of the figures – one <id>_fit.pdf, one PNG per page ("png", for a talk), or both ("both"); the pages themselves are the same either way. thin (default 1) keeps every n-th posterior sample in the FITS – the FITS is the only copy of the posterior, so a value above 1 throws samples away and is meant for a run whose file size, not whose analysis, is the problem. fit_plot (default true) draws the fit page; corner and histogram (default false) add those pages. The summary page is always there.

corner_params and hist_params select the parameters to draw, and both spellings are accepted: the qualified name (b_halpha.amp_c), which does not depend on the summation order, or the flat one (amp_c_3), which is what the FITS column order stores. A name that is not a free parameter of the model is an error listing the free ones. A 27-parameter corner plot is a 27x27 square nobody reads, which is why the IRAS file names thirteen.

The figures are pages of one PDF, written in the order [output] lists them. They are the library’s, not the pipeline’s own drawings: the fit plot is galspec.plot.plot_fit_jax() – data, total model, every component and the residuals in units of the error, the same figure example/fit_iras07546_jax.py writes – the histograms are galspec.plot.plot_posterior_hist(), and the corner is galspec.plot.plot_corner(). So a pipeline figure and the figure a hand-written script makes of the same fit agree, and the qualified LaTeX labels and prior annotations are the ones the rest of GalSpec uses. One consequence worth knowing: a corner plot needs a posterior that moved, so a run stopped before the sampler explored skips that one figure with a warning instead of failing the fit.

The model section#

[model] holds require_priors (default false) and one [[model.component]] per additive term, in the order they are summed:

[model]
# require_priors = true        # every prior must be written out to fit

[[model.component]]
type = "PowerLaw1D"            # any kind galspec.jax_models knows
name = "Continuum"             # used by ties, priors and plots
amplitude = 0.55               # optional: the model has a default too
x_0 = 5500.0                   # not optional: it is where the power law sits
alpha = 1.0
fixed = { x_0 = true }

[model.component.priors.amplitude]
kind = "loguniform"            # optional: the default is loguniform over
lower = "0.2*continuum_level"  # 1e-3 .. 3*max_flux, which is wider
upper = "4*continuum_level"

The order matters twice: it is the order of the sum, and the flat parameter names carry the index (amp_c_3 is the fourth component’s amp_c). Reordering the file therefore renames parameters; the qualified spelling is the stable one, which is why ties, priors and plots use it.

Keys of a component are one of:

type and name;
a parameter value – anything the model names as a parameter (amplitude, dv, sigma, amp_c0, …), which may be a number or an expression. Optional: an omitted value keeps the model class’s default, which nested sampling ignores anyway. Three values are not optional in practice – wavec/wavec0/wavec1 and x_0 – because their model default (5000 A, 1 A) does not describe the line or the pivot you mean: a line component has to name its centre, with a number or with "@Halpha", or pin it with [model.component.bounds.wavec];
a constructor argument (n_components, par_w, degree, template_name, Star_types, R, template_file) – these are passed to the model’s constructor and must not collide with a parameter name. A key ending in _file is a path: it is resolved relative to the configuration file (an absolute path is used as written) and never parsed as an expression;
the tables fixed, bounds, params, ties and priors.

Anything else is rejected with its TOML path, so a typo cannot become an unused keyword argument that silently does nothing.

Expressions#

Any number in [model] and [data] may be an expression over the statistics of the source being fitted, which is what lets one file serve a sample whose flux scale differs from spectrum to spectrum:

token	value for the current source
`peak_flux`	`max(\|flux\|)`, after despiking (see below)
`max_flux`, `min_flux`	`max(flux)`, `min(flux)`, after despiking
`median_flux`, `continuum_level`	the median, and the `continuum_percentile` percentile, of the raw pixels
`rms`	standard deviation of the flux
`n_pixels`	pixels being fitted
`z`	the redshift used

The language is numbers, those names, + - * / ** and parentheses. Everything else – a call, an attribute, a subscript – is rejected when the file is read, before any spectrum is touched. "3*max_flux" sets a starting value; [model.component.bounds] amplitude = ["0.01*peak_flux", "10*peak_flux"] sets the box in the same units, which is what makes a magnitude-limited sample fit with one file.

How the statistics are measured#

The tokens are computed once per source, from the pixels that are actually fitted – after the rest-frame window, the masks and the de-reddening, but in the flux unit the file carries, so a token inherits that unit and no conversion is attempted.

continuum_level and median_flux come straight from the flux array: np.percentile(flux, continuum_percentile) and np.median(flux). A percentile is what a continuum level should be – it is a statement about where the spectrum sits, not about its brightest feature – and it is already immune to a few bad pixels, so nothing is smoothed. continuum_percentile therefore has to be chosen for the spectrum: a spectrum with 60% of its window filled by emission lines wants something like 15 or 20, one that is mostly continuum wants the default 50.

max_flux, min_flux and peak_flux set the scale of every amplitude default and are measured on a despiked copy of the flux, because a single hot pixel would otherwise decide them:

a running median of five pixels, which outvotes any outlier up to two pixels wide – a cosmic ray, a bad column, a sky residual left in one channel;
a sigma clip of the residual against that running median: a pixel that stands more than three residual spreads away and agrees with neither neighbour is replaced by the median. The neighbour test is what keeps a real line, whose flanks do not look like that;
the extremes of the result.

A real line is wide compared with five pixels at any sensible dispersion, so its peak survives: the despiked maximum of the clean IRAS spectrum differs from max(flux) by 4e-4 in relative terms – while a one-pixel spike of any height does not. A narrow line loses the tip of its profile – a one-pixel line keeps 80% of its peak, a two-pixel line 93%, five pixels or more 98% – which is the right trade for a number whose job is to set a prior scale. rms keeps the raw array, since it is a noise estimate and despiking it would remove exactly what it measures.

The IRAS file uses max_flux rather than peak_flux on purpose: the script’s set_priors() uses np.max(flux), which is not max(|flux|) when the noise dips below zero, and a reproduction should not depend on that difference.

Arithmetic between constants is allowed and is used for the atomic ratios, so a starting value is the ratio itself rather than a rounded copy:

amp_c0 = 1.8
amp_c1 = "1.8/2.98"            # the [O III] 5007/4959 ratio

Line names such as "@Halpha" resolve to the rest wavelength from galspec.utils.line_wave_dict (Halpha, Hbeta, OIII_5007, NII_6583, SII_6716, HeI_5876, …); other values are written in Angstrom.

Fixed, bounds, params#

To hold a parameter at its value:

[model.component]
type = "Line_Gaussian"
name = "Halpha"
fixed = { wavec = true }       # Line_Gaussian already fixes it by default

To change a box:

[model.component.bounds]
amplitude = ["0.01*peak_flux", "5*peak_flux"]
sigma = [50.0, 1200.0]

Each parameter can also be written as its own table, which is convenient when it needs several settings:

[model.component.params.sigma]
value = 300.0
bounds = [50.0, 1200.0]
prior = { kind = "uniform", lower = 50.0, upper = 1200.0 }

Ties#

A tie makes one parameter follow another. The plain form is a reference:

[model.component.ties]
dv = "Halpha.dv"                 # this line moves with Halpha
sigma = { ref = "Halpha.sigma", scale = 1.5 }   # 1.5 * Halpha.sigma

The named kinds are the GalSpec conventions, written as reference * scale + offset so that they survive serialisation as pure data:

`kind`	reference parameter	value
`ref`, `follow`, `copy`	`param` (required)	`reference`
`doublet_ratio`	`amp_c0`	`amp_c0 / ratio` (e.g. 2.98 for [O III])
`multigauss_amp_c`	`amp_c`	`amp_c / ratio`
`template_amplitude`	`amplitude`	`amplitude / ratio`
`logtau0_ratio`	`logtau0`	`logtau0 - log10(ratio)`
`template_dv`, `absorption_dv`, `absorption_sigma`, `absorption_cf`, `starspectrum_sigma`, `starspectrum_deltaz`	that parameter	`reference`

The IRAS model is held together by 29 such ties, and they are the whole physical content of the model beyond the line centres and the fixed values: one narrow-line velocity, one broad-line kinematics, and three atomic ratios. Written as they appear in the file, for the four kinds of reader – the profile comes from the data file ([data] narrow_template_key), so no component repeats a path:

# long-wavelength narrow lines share the [O III] velocity
[[model.component]]
type = "Line_template_doublet"
name = "NII"
amp_c0 = 1.0
amp_c1 = "1.0/2.96"
dv_c = 0.0
wavec0 = "@NII_6583"
wavec1 = "@NII_6548"

[model.component.ties]
amp_c1 = { kind = "doublet_ratio", ref = "NII", ratio = 2.96 }
dv_c = "OIII.dv_c"                       # ... and [O III] owns the velocity

A component that needs a different profile from the rest of the spectrum writes template_file = "profiles/nii.txt" (two columns, velocity then flux, resolved relative to the configuration) and overrides narrow_template_key for itself.

# broad H-beta keeps its own peak and its own wind amplitude, and copies the # shape of broad H-alpha, which is the only line whose kinematics are free [[model.component]] type = “Line_MultiGauss” name = “b_hbeta” n_components = 2 amp_c = 0.22 dv_c = 0.0 sigma_c = 1500.0 wavec = “@Hbeta” amp_w0 = 0.20 dv_w0 = 0.0 sigma_w0 = 6000.0

[model.component.ties] dv_c = “b_halpha.dv_c” sigma_c = “b_halpha.sigma_c” dv_w0 = “b_halpha.dv_w0” sigma_w0 = “b_halpha.sigma_w0”

The iron pseudo-continuum is tied the same way (five kernel parameters follow broad H-alpha), so the iron costs the sampler exactly one number, its amplitude.

Ties are validated when the model is built: a reference to a component or parameter that does not exist is an error naming the alternatives. Tying a parameter that is fixed is allowed (the tie then simply never matters); putting a prior on a tied parameter is rejected, because the prior would never be used.

Priors#

A parameter with no prior is uniform over its bounds. Every prior takes the arguments of its class in galspec.priors, plus expressions where a number is expected:

[model.component.priors.amp_c0]
kind = "loguniform"            # uniform, tophat, loguniform, spikeandslab,
lower = "0.01*peak_flux"       # normal, clippednormal
upper = "5*peak_flux"

[model.component.priors.dv]
kind = "normal"
mean = 0.0
sigma = 150.0

[model.component.priors.sigma_w0]
kind = "normal"
mean = "Halpha.sigma_c"        # centred on another parameter
sigma = 0.4
relative = true
nsigma = 3

Default priors#

Nested sampling draws from the prior, so a configuration does not have to say where to start: an omitted parameter value falls back to the model class’s own default (amplitude = 1, sigma = 200 km/s, …), which never biases the answer – it only changes the point the sampler prints while it warms up.

A parameter the file says nothing about also gets a default prior, chosen from its name and resolved against the source’s own statistics:

parameter	default prior	why
`amplitude`, `amp_c`, `amp_c0`, `amp_c1`, …	`loguniform(1e-3continuum_level, 3max_flux)`	a line is between a thousandth of the continuum and three times the brightest pixel, and “twice as faint” and “twice as bright” are equally probable. Both ends are scaled by the spectrum, so one file serves a sample whose flux spans decades; an absolute `1e-3` would be a different physical limit for every source. If the resolved interval is empty, or its lower end is not positive (a continuum level of zero or below), it falls back to `1e-3*max_flux` and `--dry-run` says so.
`amp_w0`, `amp_w1`, …	`uniform(0, 1)`	a wind amplitude is a fraction of the line.
`dv`, `dv_c`	`uniform(-3000, 3000)` km/s	systematics and a few hundred km/s of gravitational slop.
`dv_w0`, `dv_w1`, …	`uniform(-6000, 6000)` km/s	a wind is wider and can be more offset.
`sigma`, `sigma_c`, `stddev`	`uniform(10, 8000)` km/s	any width a spectrograph can resolve, from a stellar line in a high-resolution spectrum to a broad-line-region profile. A broad line wants a floor near 300 km/s, which is a statement about the object rather than about units, so write it out – the IRAS file does.
`sigma_w0`, `sigma_w1`, …	`uniform(10, 8000)` km/s	the wing, i.e. an outflow.
`alpha`	`uniform(-5, 5)`	a power-law index, both signs.
`h3`, `h4`	`uniform(-0.5, 0.5)`	Hermite moments.
`logtau`	`uniform(-3, 3)`	optical depth of an absorption model.
`cf`	`uniform(0, 1)`	covering fraction.
`z`, `redshift`	`uniform(-0.1, 0.1)`	a redshift correction.
`velscale`	`uniform(10, 1000)` km/s	a broadening kernel.

A prior or a bound written in the file always wins, so the defaults are what is in force only where the file is silent. --dry-run lists every parameter that took one, with the numbers it resolved to, and the same list is in run_report.json under model.default_priors: nothing about the prior is invisible.

Two kinds of parameter are deliberately not defaulted, because the right range depends on the data in a way a name cannot express – a linear continuum’s slope and intercept, and a free wavec: a wrong box there changes the evidence silently instead of failing. Those still need

[model.component.bounds.slope]
...  # = [low, high]

require_priors = true demands a prior written in the file on every free parameter – defaults do not count – so it is the way to make a configuration state its priors even though it no longer has to. Two checks run before the sampler starts, so a mistake costs a second instead of a night:

every free parameter needs a finite sampling range – both bounds, or a prior with finite limits (a centred prior such as normal additionally needs the bounds, because they truncate it);
no free parameter means nothing to fit, which is also an error.

The IRAS file leaves twenty-one of its twenty-seven priors to the defaults: every line amplitude (loguniform(1e-3*continuum_level, 3*max_flux), the reference script’s loguniform(1e-3, 3*max_flux) with a floor that scales instead of one absolute number), alpha, the three wind amplitudes and one wind velocity. Six are written out, and each is worth knowing about when writing a new file:

the broad lines’ widths: 300-8000 km/s, while the default floor is 10 km/s, which a narrow line needs. A floor of 300 is a statement about the object, not about units, so it cannot be a default;
the continuum’s 0.2-4*continuum_level and the iron template’s 3*continuum_level ceiling – a pseudo-continuum is scaled by the continuum, not by the brightest emission line;
the two narrow-line velocity windows, broad H-alpha at 2500 km/s and [O III] at 1000 km/s, which the systemic reference needs.

The scaled amplitude floor is not cosmetic for this object: 1e-3 of a continuum level of 0.4991 is 5e-4, half the script’s absolute floor, so a log-uniform prior with the same ceiling is wider by ln(12.99/5e-4) / ln(12.99/1e-3) = 1.074 in its log-range, and a log-uniform density scales as 1/ln(b/a). Each of the sixteen amplitudes therefore lowers the evidence by ln(1.074) = 0.071, and the total is about -1.1: the published 130677.27 belongs to the script’s prior, and a file that uses the scaled floor would report about 130676.1.

Only the number moves. The two priors differ by a constant factor wherever both are non-zero, so the posterior is the same up to the renormalisation the added faint tail implies – and that tail (amplitudes between 5e-4 and 1e-3, where no line of this spectrum lives) can only bring the shift back toward zero. The offset is also worth putting next to the ±0.58 sampling error: it is a deliberate difference in the family of models, not a defect in the fit. A short run cannot measure it either – a nested-sampling run truncated at a few hundred iterations is dominated by its un-converged initial state, and the numbers it reports differ by far more than one nat – which is why the value above is the exact prior-volume ratio rather than a measured difference. The mapping is one line of the script per table of the file, which is the property to aim for when writing a new configuration:

in `fit_iras07546_jax.py`	in `fit_iras07546_config.toml`
`build_continuum()`	the `Continuum` component
`build_iron_template()` + `tie_iron_template()`	the `FeII` component and its five ties
`build_narrow_lines()` + `tie_narrow_lines()`	eight `Line_template`/`Line_template_doublet` components, their `dv` ties and three ratio ties
`build_broad_lines()` + `tie_broad_lines()`	five `Line_MultiGauss` components and their ties
`set_priors(model, flux)`	the six `[model.component.priors.*]` tables, plus the pipeline’s default prior for the other twenty-one parameters
`JNesty_Fit(..., nlive=NLIVE, dlogz=DLOGZ, ...)`	`[fit]`

Output#

For a source <id>:

<output_dir>/<id>/<id>_fit.fits    # the model, the data and the whole posterior
<output_dir>/<id>/<id>_fit.pdf     # a summary page, then the figures
<output_dir>/run_report.json       # what was fitted, with what, how

Two files per source, and nothing else. There is no <id>_posterior.npz and no directory of .png figures: the FITS already holds every array such an archive would hold, and a second copy of a posterior is a second thing to keep in step with the first. (A run made before this changed wrote the archive; its arrays are in the FITS under the names below, which is what a script reading np.load should read instead.)

<id>_fit.fits is the file galspec.jnesty_fit_io.save_jnesty_fit_to_fits() writes. The primary header carries the provenance – source identifier, input file, pixels, redshift and where it came from, the extinction removed (EBV, EBVRV, EXTINCT), the configuration’s name and the first 16 characters of its SHA-256, the fit date – and the sampler’s verdict: LOGZ/LOGZERR, NITER (iterations, i.e. dead points), NSAMP (rows in the sample arrays – NITER + NLIVE before thinning), NCALL, NLIVE, NDIM, BOUND, THIN, DEVICE, RUNTIME, and the diagnostics INFORM (the Bayesian information H), ACCRATE, EFF, CONVERGED, DLZ, DLZTHR.

The extensions hold the run:

HDU	what it holds
`SAMPLES`, `SAMPLES_EQ`	every sample (`NSAMP` rows – `NITER` dead points plus the live points, or fewer under `thin`; one column per free parameter) and the equal-weight resample the percentiles are computed from
`LOGL`, `LOGWT`, `LOGVOL`, `SAMPLES_U`	one entry per sample: log-likelihood, log-weight, log-volume, and the same point in the unit cube
`LOGZ`, `LOGZERR`, `DLZ_TRAJ`, `SCALE_TR`	the run’s trajectories, which is what the dynesty-style `runplot` / `traceplot` diagnostics are drawn from
`PARAM_NAMES`, `FULL_PARAM`, `PARAM_BOUNDS`, `LOG_SCALE_PARAMS`	the column order of `SAMPLES`, the qualified names, the prior limits, and which parameters were sampled in log space
`THETA_INIT`, `THETA_BEST`	the point the fit started from and the maximum-likelihood sample
`WAVE_USE`, `FLUX_USE`, `FERR`	the arrays that were fitted, after the cuts and the two corrections
`MODEL_STR`, `MODEL_ARR`, `MODEL_PARAM`, `SUBMODELS`, `CUSTOM_PRIORS`	the model itself, so the fit can be redrawn without the configuration file

That is enough to re-derive every number in this chapter from the file alone:

from astropy.io import fits
hdul = fits.open("out/<id>/<id>_fit.fits")
samples = hdul["SAMPLES"].data             # (niter, ndim), in PARAM_NAMES order
weights = hdul["LOGWT"].data               # ln(weight) of each sample
equal = hdul["SAMPLES_EQ"].data            # already resampled
print(hdul[0].header["LOGZ"], hdul[0].header["INFORM"])

With figures = "png" the same pages are written beside the FITS as <id>_summary.png, <id>_fit.png, <id>_corner.png, <id>_histogram.png, and figures = "both" writes both sets; nothing else about the run changes, so one can be re-rendered from the other without re-sampling – except that only the PDF can be re-rendered from anything, because the FITS carries every number the pages are drawn from.

<id>_fit.pdf is one document. Page one is text – the evidence and the run’s diagnostics, the data and the two corrections the pipeline applied, and a table of every free parameter with its median, 16th and 84th percentile and maximum-likelihood value, marking with * a parameter on the pipeline’s default prior and with ! one pressed against a prior limit (a limit, not a measurement). Then come the pages [output] asks for: the fit, the corner plot, the marginal posteriors. The pages are vector drawings, so a 57000-pixel spectrum stays sharp at any zoom – which is also why losing the separate PNGs costs nothing.

run_report.json records the settings, library versions, and per source: the status, evidence, iterations, timing (compile and sampling separately), a parameters table with the median, the 16th and 84th percentile and the maximum-likelihood value of every free parameter, which parameters are pressed against a prior limit (prior_dominated), which parameters took a default prior (model.default_priors, with the numbers it resolved to for this source), the flux tokens that were used, the redshift and extinction provenance, any warning raised while reading, and the outputs. A source that fails is recorded with its error and the run continues; the report is written even when every source fails.

Resuming and parallelism#

A source whose output files all exist is skipped unless force = true (or --force), so an interrupted run is continued by running the same command again. workers > 1 fits sources in spawned processes; each process writes its own sources, and the report is assembled from their results. Failures are isolated per source in both modes. Set GALSPEC_PIPELINE_TRACE=1 to keep the traceback in the report.

Compiling once for a sample#

Every JAX program has to be compiled before it runs, and compiling the likelihood and the prior of a full model takes a fraction of a second to a second (0.6 s for the 39-parameter example in develop/dev_jax/bench_shared_compilation.py). With share_compilation = true – the pipeline’s default – the spectrum travels to the compiled likelihood as an argument, so the compiled program does not contain the data and is reused for the next source. In a four-worker run of develop/dev_config/mock_sample.toml the report shows 0.26-0.38 s of compiling for the first source each worker fitted and 0.01 s for the second.

JNesty’s own loops accept the same context, so its kernels are reused as well, provided the model structure (not the data) stays the same. The report records compile_seconds and sample_seconds separately, so the saving – or its absence – is visible per source rather than argued about.

The one thing that defeats the reuse is a model that carries per-source numbers of its own: a narrow-line profile read out of each spectrum ([data] narrow_template_key) is a constant inside the traced program, not an argument, so two sources with different profiles trace two different programs. The cost is one recompile of the whole likelihood kernel – there is no “template part” to compile separately, because the profile is arithmetic inside the same expression as the model sum and the residual. Measured on the 27-parameter IRAS model (56956 pixels, CPU): 0.5-0.7 s for a new profile, 0.07 s when the profile repeats (the trace plus a registry hit). The prior transform is a separate program with no model in it, so it is reused; JNesty’s own sampling kernels never touch the model and are not recompiled at all.

Avoiding even that would mean passing the profile as a runtime array argument, which needs the model to accept one – a change to galspec.jax_models, and the pipeline deliberately changes no model. Half a second per source against minutes of sampling is the trade the file takes; the report records it per source as compile_seconds.

The registry hashes the program including its constants, which matters here: JAX prints a large array constant in abbreviated form, so two models that differ only in their profile trace to the same text, and a key made from the text alone would hand the second source the first source’s kernel – a fit with the wrong line shape, and nothing downstream would show it. tests/test_shared_compilation.py checks that a model differing only in a constant compiles again and gives its own likelihood.

This is a compilation saving, not a speed-up of the sampling: the compiled arithmetic is the same, and per-iteration cost is unchanged within run-to-run noise. The numbers do not change either – the same computation is traced, and tests/test_shared_compilation.py checks that the likelihood, the prior transform and a short run agree with the per-source path. compilation_cache_dir additionally keeps JAX’s on-disk cache, which is what helps across processes (workers > 1) and across runs.

Command line#

python -m galspec.pipeline CONFIG.toml [options]

--template [PATH]  write a starter configuration (PATH defaults to config.toml;
                   "-" prints it) and exit; --force overwrites
--dry-run          read every spectrum and report what would be fitted
--list-sources     print the source ids and files the config selects
--workers N        fit in N processes
--force            refit even when the outputs exist
--output DIR       write here instead of [run] output_dir
--device CPU|GPU   override [fit] device

--quiet            only failures and the final summary

A fit prints nothing between the fitting line and its result, which for a 27-parameter model on a CPU means hours: nested sampling converges by itself and has nothing to report in between. In a terminal the sampler draws a progress bar, so silence there means something else; when the output is redirected the bar is suppressed (it is written with carriage returns, which is noise in a log) and the run says so before it starts. [run] progress forces either behaviour, and --quiet turns both the log and the bar off.

--dry-run is the one to run first on a new sample: it reads every file, resolves every expression against that source’s flux, builds the model, and reports the pixels, the rest-frame range, the redshift and its source, the tokens and the number of free parameters – without sampling. For the IRAS file:

$ ./example/fit_iras07546_config.sh --dry-run
[1/1] iras07546_pipeline: 56956 pixels, 4200.0-6800.0 A rest, z=0.0953 (header:REDSHIFT), peak_flux=4.332, 27 free parameters

example/fit_iras07546_config.sh is a two-line wrapper around the command; it passes its arguments through, so --device cpu or --output work with it too. The same thing from Python:

from galspec.pipeline import load_config, run_pipeline

config = load_config("sample.toml")
report = run_pipeline(config)            # or workers=8, force=True
print(report["counts"], report["sources"][0]["log_evidence"])

The IRAS reproduction#

example/fit_iras07546_config.toml is the model of example/fit_iras07546_jax.py – seventeen additive components, 27 free parameters, 56956 pixels – and it reproduces the published run:

./example/fit_iras07546_config.sh                # ~12 min on a GPU
./example/fit_iras07546_config.sh --device cpu   # hours

The claim is checkable and is checked by develop/dev_config/test_iras07546_config.py, which runs in a few minutes:

the spectrum equals iras07546_model.load_spectrum() to 1.3e-16 relative – the best possible, since both are the same float64 numbers;
the model is the same: same 17 components in the same order, same 27 free parameters in the same order, same starting values, bounds, priors (kind and numbers) and the same 29 ties;
prior_transform and the log-likelihood agree point by point, so the sampling problem is identical;
a capped fit (50 live points, 300 iterations) drives the sampler end to end and writes every output.

Published numbers to compare against (example/iras07546_model.REFERENCE):

Gaussian narrow lines          30 free   log Z = 130889.60 +/- 0.62   (10.2 min)
measured narrow profile        27 free   log Z = 130677.27 +/- 0.58   (11.7 min)   <- this file

Repeat runs of one configuration have given 130656.22, 130656.34 and 130658.06 – about two nats of scatter, from the last bits of a GPU reduction reaching the sampler’s ordering – and re-extracting the narrow-line profile moved the evidence by about twenty. Quote an evidence from this family to a nat or two; the differences the model is judged on are two hundred. What the file reproduces exactly is the problem – the data, the model, the priors – not the floating-point trajectory of the sampler.

develop/dev_config/index.html (written by develop/dev_config/make_report.py) is the illustrated version of this walkthrough: it draws the spectrum the converter writes, the model at its starting values split into continuum, iron, narrow and broad, the prior boxes, and the parity figures.

A small example to copy from#

develop/dev_config/mock_two_line.toml fits the two-line mock spectrum that ships with GalSpec – a linear continuum plus two Gaussians, one tied velocity, token priors:

python -m galspec.pipeline develop/dev_config/mock_two_line.toml --dry-run
python -m galspec.pipeline develop/dev_config/mock_two_line.toml

The fit recovers the line centres to a fraction of a pixel and the two very different widths of the mock (about 1300 and 160 km/s), with residuals consistent with the 1% noise. develop/dev_config/mock_sample.toml is the same model on six spectra of different brightness: python develop/dev_config/make_mock_sample.py writes them, then --workers 4 fits them in about 45 s and the report shows both widths recovered in every source while the amplitudes follow each source’s peak_flux.

These two live in develop/ because they are tests of the notation, not science; example/ holds the two examples a reader should start from (fit_iras07546_jax.py in Python and fit_iras07546_config.toml in TOML).

Limitations#

One FITS file per source; no row-slicing of a multispec file.
No error synthesis, no sky-line masking and no flux-unit conversion: prepare the file you want fitted. The two corrections the pipeline does apply, the redshift and the Galactic extinction, are read from the header and recorded.
The model is a sum of components; a product (e.g. a convolution wrapper) is not expressible in the file yet.
A prior cannot be nested inside another prior (spikeandslab uses its default slab).
dynamic.py’s DynamicNestedSampler is not covered by the shared compilation path; JNesty_Fit uses NestedSampler.