Model fitting

OITOOLS supports parametric model fitting using a flat-dictionary interface compatible with PMOIRED.

Workflow overview

Model fitting in OITOOLS follows three steps:

1. Define a model dictionary (model_dict) — a Dict{String,Any} listing all parameters for every component. Values are either numbers (free or fixed) or expression strings with \$-references for derived quantities.
2. Choose which parameters to fit (list_free_params) — a Vector{String} naming the subset of keys in model_dict that the optimizer is allowed to vary. Optionally define lower/upper bounds (lb, ub) for each.
3. Fit — calling fit_model (or fit_model_lsqfit, fit_model_nested) compiles the dictionary into an efficient FlatModel internally, then optimizes the free parameters against the data.

The result contains the best-fit values (result.x_opt) and the compiled FlatModel (result.model), which you can pass to model_to_obs, model_to_image, etc.

To compile a model dictionary into a FlatModel:

model = dict_to_model(model_dict, list_free_params)
x0 = Float64[model_dict[p] for p in list_free_params]

Then pass model and x0 to fitting and evaluation functions.

Defining a model dictionary

Parameters use flat keys of the form "component,parameter":

using OITOOLS, PythonPlot
data = readoifits("data/AlphaCenA.oifits")[1]

model_dict = Dict{String,Any}(
    "star,ud"  => 8.0,    # uniform disk diameter (mas)
    "star,f"   => 1.0,    # flux fraction
)
list_free_params = ["star,ud"]  # parameters to optimize

display_model(model_dict, list_free_params)

The component name (e.g. "star", "disk", "ring") is arbitrary — you choose it. The component type is determined automatically from which geometry key is present.

Component types

Uniform disk

A single key ud sets the diameter in mas:

model_dict = Dict{String,Any}("star,ud" => 1.0, "star,f" => 1.0)

Gaussian

A single key fwhm sets the full-width at half-maximum in mas:

model_dict = Dict{String,Any}("g,fwhm" => 0.5, "g,f" => 1.0)

Limb-darkened disks

Five limb-darkening laws are available. The geometry key gives the diameter in mas, and names the law:

# Linear: I(μ) = 1 - u(1-μ)
model_dict = Dict{String,Any}("star,ldlin" => 1.0, "star,u" => 0.3, "star,f" => 1.0)

# Quadratic: I(μ) = 1 - u(1-μ) - w(1-μ)²
model_dict = Dict{String,Any}("star,ldquad" => 1.0, "star,u" => 0.2, "star,w" => 0.1, "star,f" => 1.0)

# Square root: I(μ) = 1 - u(1-μ) - w(1-√μ)
model_dict = Dict{String,Any}("star,ldsqrt" => 1.0, "star,u" => 0.3, "star,w" => 0.2, "star,f" => 1.0)

# Power-law: I(μ) = μ^α   (Hestroffer)
model_dict = Dict{String,Any}("star,ldpow" => 1.0, "star,alpha" => 0.5, "star,f" => 1.0)

# Four-parameter: I(μ) = 1 - Σₖ cₖ(1 - μ^{k/2})   (Claret)
model_dict = Dict{String,Any}("star,ldclaret4" => 1.0, "star,c1" => 0.4, "star,c2" => -0.2,
                              "star,c3" => 0.5, "star,c4" => -0.15, "star,f" => 1.0)

Uniform ring

Defined by inner and outer diameters in mas:

model_dict = Dict{String,Any}("ring,diamin" => 0.6, "ring,diamout" => 1.0, "ring,f" => 1.0)

A convenience sugar is available using diam (outer diameter) and thick (fractional thickness 0–1):

# Equivalent to diamin = 2.0*(1-0.3) = 1.4, diamout = 2.0
model_dict = Dict{String,Any}("ring,diam" => 2.0, "ring,thick" => 0.3, "ring,f" => 1.0)

Gaussian ring

A bi-Gaussian ring defined by inner and outer FWHM in mas:

model_dict = Dict{String,Any}("gr,fwhmin" => 0.2, "gr,fwhmout" => 0.5, "gr,f" => 1.0)

Crescent

Two offset uniform disks producing a crescent shape:

model_dict = Dict{String,Any}(
    "cr,crin"      => 0.8,    # inner disk diameter (mas)
    "cr,crout"     => 1.0,    # outer disk diameter (mas)
    "cr,croff"     => 0.8,    # offset factor (0 = concentric ring, 1 = max offset)
    "cr,crprojang" => 120.0,  # PA of thinnest part (degrees)
    "cr,f"         => 1.0,
)

Point source

A component with only f (no geometry key) is an unresolved point source:

model_dict = Dict{String,Any}("companion,f" => 0.05, "companion,x" => 2.0, "companion,y" => -1.0)

Fully resolved background

Sets V = 0 at all non-zero baselines (fully resolved flux):

model_dict = Dict{String,Any}("bg,resolved" => true, "bg,f" => 0.1)

Summary table

V(B) = exp(-π² · FWHM² · B² / (4 ln 2))

Common parameters

Every component supports these optional keys:

Geometric transformations

Inclination and position angle project the component on sky:

# Inclined uniform disk
model_dict = Dict{String,Any}(
    "star,ud"      => 1.0,
    "star,incl"    => 60.0,    # degrees from face-on
    "star,projang" => 30.0,    # PA in degrees
    "star,f"       => 1.0,
)

# Offset companion
model_dict = Dict{String,Any}(
    "star,ud" => 1.0, "star,f" => 0.9,
    "comp,f"  => 0.1, "comp,x" => 3.0, "comp,y" => -1.5,
)

Expression strings

Parameter values can be strings referencing other parameters with \$. References are resolved in topological order, so forward and backward references both work:

model_dict = Dict{String,Any}(
    "star,ud"      => 3.0,
    "star,f"       => 0.7,
    "disk,f"       => "1 - \$star,f",            # complement of star flux
    "disk,diamout" => "\$star,ud * 8",           # 8× the stellar diameter
    "disk,diamin"  => "\$disk,diamout * 0.5",    # half the outer diameter
)

Implicit variables

These implicit variables are available in expressions:

data = readoifits("target.oifits"; T = Float64)[1, 1]

Referencing \$WL or \$MJD in any derived expression makes the resolver broadcast all derived expressions, so every parameter becomes a per-uv-point vector — chromatic diameters and position angles work exactly like chromatic fluxes.

These enable chromatic and time-variable models:

model_dict = Dict{String,Any}(
    "star,f"        => 1.0,
    "star,spectrum"  => "(\$WL/2.2e-6)^(-4)",    # Rayleigh–Jeans spectrum
    "disk,f"        => 0.3,
    "disk,spectrum"  => "(\$WL/2.2e-6)^(-2)",    # greyer spectrum
    "disk,diamout"  => 10.0,
    "disk,profile"  => "exp(-(\$R/3.0)^2)",
)

They are not restricted to fluxes: any parameter may be an expression, including the geometric ones, in which case it takes a different value at every uv point:

model_dict = Dict{String,Any}(
    # a diameter that grows linearly with wavelength
    "star,d0"    => 3.0,
    "star,slope" => 0.05,
    "star,ud"    => "\$star,d0 * (1 + \$star,slope * (\$WL/1.6e-6 - 1))",
    "star,f"     => 1.0,
)
list_free_params = ["star,d0", "star,slope"]   # both are fitted as usual

This works for every analytic component — ud, ldlin, ldquad, ldsqrt, ldclaret4, ldpow, fwhm, the rings and the crescent — and for every one of their parameters, not only the size: a limb-darkening coefficient varying with time ("star,u" => "0.3 + 0.01*(\$MJD - 60000)") is equally valid. Gradients follow automatically, so such models can be fitted with fit_model or fit_model_lsqfit without any special handling.

Shared parameters across components

Bare keys (without a component, prefix) act as global variables that multiple components can reference:

model_dict = Dict{String,Any}(
    "PA"              => 60.0,        # global position angle
    "INC"             => 45.0,        # global inclination
    "inner,fwhm"      => 1.0,
    "inner,projang"   => "\$PA",
    "inner,incl"      => "\$INC",
    "outer,diamin"    => "2 * \$inner,fwhm",
    "outer,diamout"   => "3 * \$outer,diamin",
    "outer,projang"   => "\$PA",
    "outer,incl"      => "\$INC",
    "outer,f"         => 0.5,
)

Combining multiple components

Models with multiple components are built by giving each component a different name prefix. The total model visibility is the flux-weighted sum of all component visibilities:

model_dict = Dict{String,Any}(
    # Compact star
    "star,fwhm"      => 0.1,
    "star,spectrum"   => "\$WL^(-3)",
    # Inclined disk with profile and azimuthal modulation
    "disk,diamin"    => 0.5,
    "disk,diamout"   => 1.0,
    "disk,profile"   => "\$R^(-2)",
    "disk,az amp1"   => 1.0,
    "disk,az projang1" => 60.0,
    "disk,projang"   => 45.0,
    "disk,incl"      => -30.0,
    "disk,x"         => -0.05,
    "disk,y"         => 0.05,
    "disk,spectrum"   => "5 * \$WL^(-2)",
)

Visualise the model image and SED:

params = dict_to_model(model_dict, String[])
img = model_to_image(params, Float64[]; nx=128, pixsize=0.02, wl=1.65e-6)
imdisp(img; pixsize=0.02)

wl_grid = collect(range(1.0e-6, 2.5e-6; length=200))
f_total, f_comps = model_to_sed(params, Float64[], wl_grid)

See example_model_fitting_pmoired_models.jl for a full gallery of all component types and combinations.

Hankel models (radial profiles)

For components that cannot be described by simple analytic visibility functions, OITOOLS uses a numerical Hankel transform of an arbitrary radial profile. This is ideal for circumstellar disks, envelopes, and limb-darkened models with custom intensity laws.

Profile expressions

Define a profile expression string using \$R (radius in mas) and \$MU (cosine of the zenith angle), plus any scalar parameters:

# Disk with intensity ∝ μ^0.5 (limb-darkening)
model_dict = Dict{String,Any}(
    "star,diam"    => 1.0,
    "star,profile" => "\$MU^0.5",
    "star,f"       => 1.0,
)

# Ring with power-law profile
model_dict = Dict{String,Any}(
    "ring,udout"   => 1.0,
    "ring,profile" => "(\$R > 0.25) * \$R^(-0.5)",
    "ring,f"       => 1.0,
)

# Parabolic ring profile (YSO disk)
model_dict = Dict{String,Any}(
    "disk,profile"  => "max(0, 1 - (2*(\$R - \$Rmid)/\$width)^2)",
    "disk,diamout"  => 10.0,     # outer diameter (mas) — sets the r-grid extent
    "disk,nr"       => 200,      # radial grid points (default: 100)
    "disk,Rmid"     => 3.0,      # peak radius (custom parameter)
    "disk,width"    => 2.0,      # ring width (custom parameter)
    "disk,f"        => 0.5,
    "disk,incl"     => 45.0,
    "disk,pa"       => 120.0,
    "star,f"        => "1 - \$disk,f",
)

Hankel-specific parameters

Azimuthal modulations

Azimuthal variations are added with harmonic coefficients az ampN and az projangN for harmonic order N = 1, 2, 3, ...:

model_dict = Dict{String,Any}(
    "disk,diamin"       => 1.0,
    "disk,diamout"      => 3.0,
    "disk,profile"      => "1",
    "disk,projang"      => -20.0,
    "disk,incl"         => 60.0,
    "disk,f"            => 1.0,
    "disk,az amp1"      => 0.3,     # m=1 harmonic amplitude
    "disk,az projang1"  => 45.0,    # m=1 orientation (degrees)
    "disk,az amp2"      => 0.1,     # m=2 harmonic amplitude
    "disk,az projang2"  => 90.0,    # m=2 orientation (degrees)
)

Multiple harmonics combine additively to produce asymmetric features such as one-armed spirals or brightness asymmetries in disks.

Spiral patterns

By building N concentric rings whose az projang1 rotates with radius, you can create spiral patterns. Each ring is a thin annulus with an m=1 azimuthal modulation; the progressive phase offset produces a spiral arm. Add a spatial_kernel to smooth the result.

See example_model_fitting_pmoired_models.jl (section 5b) for a complete spiral example with 5 rings using shared global parameters (Din, Dout, INCL, PROJANG, pitch, PAin).

YSO disk fitting example

example_model_fitting_v1295aql.jl fits three YSO disk models (parabolic ring, double sigmoid, α-sigmoid) to V1295 Aql (HD 190073) H-band data, reproducing Ibrahim et al. 2023 (ApJ, 947, 68). It demonstrates chromatic models with wavelength-dependent spectra.

Spatial smoothing

The global parameter spatial_kernel applies Gaussian smoothing (in visibility space) to the model. This is useful for softening sharp edges:

model_dict = Dict{String,Any}(
    "cr,crin" => 0.5, "cr,crout" => 1.0,
    "cr,croff" => 0.8, "cr,crprojang" => 120.0,
    "cr,incl" => 45.0, "cr,projang" => 30.0,
    "cr,f" => 1.0,
    "spatial_kernel" => 0.2,    # Gaussian FWHM in mas
)

PMOIRED compatibility

The parameter dictionary format is designed to be directly compatible with PMOIRED. Models written for PMOIRED can be ported to OITOOLS with minimal changes.

Converting Python dicts to Julia

pmoired_to_dict() converts a PMOIRED Python dict literal string directly to a Julia Dict:

model_dict = pmoired_to_dict("{'star,ud': 3.2, 'ring,f': '1 - \$star,f'}")

The lower-level pmoired_to_julia() returns the Julia source string instead, if you need to inspect or edit it before evaluation.

It handles:

Converting entire scripts

pmoired_to_julia_file() converts a Python/PMOIRED script file line-by-line:

pmoired_to_julia_file("pmoired_model.py", "julia_model.jl")

Full conversion example

Given a PMOIRED model:

# PMOIRED (Python)
param = {
    'star,ud':    3.2,
    'star,f':     0.7,
    'ring,udout': '\$star,ud * 8',
    'ring,f':     '1 - \$star,f',
    'ring,incl':  30.0,
}
fitOnly = ['star,ud', 'star,f', 'ring,udout', 'ring,f']

Convert and use in Julia:

using OITOOLS, PythonPlot

# Option 1: inline conversion
model_dict = pmoired_to_dict("{'star,ud': 3.2, 'star,f': 0.7, 'ring,udout': '\$star,ud * 8', 'ring,f': '1 - \$star,f', 'ring,incl': 30.0}")

# Option 2: file conversion
pmoired_to_julia_file("pmoired_model.py", "julia_model.jl")
include("julia_model.jl")

# Option 3: write it directly in Julia (recommended)
model_dict = Dict{String,Any}(
    "star,ud"    => 3.2,
    "star,f"     => 0.7,
    "ring,udout" => raw"$star,ud * 8",
    "ring,f"     => raw"1 - $star,f",
    "ring,incl"  => 30.0,
)
list_free_params = ["star,ud", "star,f", "ring,incl"]

See example_model_fitting_pmoired_conversion.jl for more conversion examples, and example_model_fitting_pmoired_models.jl for a gallery of all PMOIRED model types reproduced in OITOOLS.

Fitting with NLopt (gradient-based)

fit_model uses NLopt for optimisation. Gradient-based methods (default :LD_LBFGS) use the analytic Wirtinger-derivative chain rule; gradient-free methods like :LN_NELDERMEAD are available for non-smooth profiles.

result = fit_model(model_dict, list_free_params, data;
    weights = [1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],  # V² + T3φ
    method  = :LD_LBFGS,
    maxeval = 500,
    verb    = true)

println("Best χ²/ν = ", result.chi2r)
println("Parameters: ", result.x_opt)

The 7-element weights vector controls per-observable weighting: [V², T3amp, T3φ, visamp, visφ, flux, diffφ].

Bounds and priors

lb = Dict("star,ud" => 0.1)   # lower bounds
ub = Dict("star,ud" => 20.0)  # upper bounds

# Gaussian priors: a Vector of (expression, target, sigma) tuples.
# The first element is an *expression* over model parameters, not just a parameter name,
# so a prior can constrain a derived quantity:
priors = [("star,ud",           8.5, 0.5),      # prior on one parameter
          ("star,f + disk,f",   1.0, 0.01)]     # prior on their sum

result = fit_model(model_dict, list_free_params, data; lb, ub, priors)

Each prior adds (value - target)^2 / sigma^2 to the objective. Priors are supported by the Dict-based fit_model method only — the FlatModel method rejects them, since the expressions have to be compiled into the model.

default_bounds(model_dict, list_free_params) suggests bounds for every free parameter if you do not want to write them out by hand; display_model then validates values against them. Pass data and the angular-size ceiling comes from the coverage itself — 2 λ/B_min, the largest scale the shortest baseline senses — rather than from a constant:

lb, ub = default_bounds(model_dict, list_free_params; data)

Constraints

A bound is a box, one parameter at a time. A relation between parameters is not, and needs a ModelConstraint:

constraints = [ModelConstraint("ring,fwhmout", ">", "ring,fwhmin"),   # a ring has an outside
               ModelConstraint("star,f + disk,f", "=", 1.0)]          # fluxes sum to one

result = fit_model(model_dict, list_free_params, data; lb, ub, constraints)

op is one of <, <=, >, >=, =. Either side may be an expression, and the right may be a number. tol (default 1e-3) is how much violation counts as none.

fit_model hands these to NLopt as real nonlinear constraints, so they hold at the optimum; an algorithm that cannot take them — including the default :LD_LBFGS — is wrapped in :AUGLAG rather than replaced. fit_model_lsqfit and fit_model_nested have no such machinery and use a one-sided quadratic penalty instead, which is soft: a steep enough χ² surface can overrule it. Use fit_model when a constraint must hold.

check_constraints(constraints, model_dict) says which ones the starting model already satisfies, and warns about the rest — a constraint written backwards looks exactly like a bad starting guess otherwise.

Saving a model and its fit settings

A model dict does not describe a fit on its own: the free list, the bounds, the constraints and the priors all change the answer. A TOML model file carries all five.

write_model_file("binary.toml", model_dict; free = list_free_params, lb, ub, constraints, priors)

m = read_model_file("binary.toml")
result = fit_model(m.model, m.free, data; m.lb, m.ub, m.constraints, m.priors)

free = ["star,ud", "disk,pa"]

[model]
"star,ud"   = 6.5
"disk,fwhm" = "$star,ud * 3"

[bounds]
"star,ud" = [0.0, 20.0]

[[constraints]]
param = "disk,fwhm"
op    = ">"
value = "star,ud"
tol   = 0.001

[[priors]]
expr   = "star,ud"
target = 6.0
sigma  = 0.5

Keys are written sorted, so the same model always produces the same bytes and two model files can be diffed. To save a fitted model, merge the result back in first:

fitted = merge(model_dict, Dict(zip(result.list_free_params, result.x_opt)))
write_model_file("binary_fitted.toml", fitted; free = result.list_free_params, lb, ub)

Fitting with LsqFit (Levenberg-Marquardt)

fit_model_lsqfit uses the Levenberg-Marquardt algorithm from LsqFit.jl. It provides parameter covariance and 1σ error bars from the Jacobian:

result = fit_model_lsqfit(model_dict, list_free_params, data;
    weights = [1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
    lb = Dict("star,ud" => 0.1),
    ub = Dict("star,ud" => 20.0),
    maxIter = 200)

println("Best χ²/ν = ", result.chi2r)
println("Parameters: ", result.x_opt)
println("1σ errors:  ", result.stderror)
println("Covariance:\n", result.covar)
println("Converged:  ", result.converged)

The analytic Jacobian is computed via the Wirtinger chain rule through the full observable pipeline (V² → T3 → residuals), so no finite differences are needed.

Bayesian inference (nested sampling)

fit_model_nested returns a posterior and the Bayesian log-evidence, which is what lets two models be compared rather than merely fitted. Bounds are required for every free parameter: they are the prior.

Two samplers implement it, and neither is a dependency of the package — loading one activates its extension:

using OITOOLS, Nautilus

result = fit_model_nested(model_dict, list_free_params, data;
    lb = Dict("star,ud" => 0.1),
    ub = Dict("star,ud" => 20.0),    # bounds required for all list_free_params
    nactive = 400,
    cornerplot = true)

println("log(Z) = ", result.logz, " ± ", result.logzerr)
println("Best χ²/ν = ", result.chi2r)
println("sampler  = ", result.backend)

The posterior samples are available as result.posterior, a nsamples × nparams matrix of equally weighted samples whichever backend produced them.

nested_backend() reports which sampler will run and set_nested_backend! chooses between them when both are loaded. With both available it is worth running each once: logz from the two should agree within logzerr, and if it does not, one of the runs has not converged.

using OITOOLS, Nautilus, PythonCall           # both backends

rj = fit_model_nested(model_dict, list_free_params, data; lb, ub, backend = :nautilus)
ru = fit_model_nested(model_dict, list_free_params, data; lb, ub, backend = :ultranest)
abs(rj.logz - ru.logz) <= 3 * sqrt(rj.logzerr^2 + ru.logzerr^2)   # they should agree

fit_model_ultranest remains as an alias pinning backend = :ultranest.

Inspecting results

Model observables

obs = model_to_obs(result.model, result.x_opt, data)
# obs.v2, obs.t3amp, obs.t3phi, obs.visamp, obs.visphi

Synthetic images

img = model_to_image(result.model, result.x_opt;
    nx=256, pixsize=0.1, wl=1.65e-6)
imdisp(img; pixsize=0.1)

Spectral energy distribution

wl_grid = range(1.5e-6, 2.5e-6, length=100)
total_flux, component_fluxes = model_to_sed(result.model, result.x_opt, wl_grid)

Uncertainty estimation

Two estimates of the parameter uncertainties are available, and they do not measure the same thing:

Both are calibrated when the error bars are right; only the bootstrap survives correlated calibration errors, where the analytic covariance is 3–5× too small even after its χ²ᵣ rescaling. The Uncertainties page measures both against simulated truth and compares them with PMOIRED.

Block bootstrap

boot = bootstrap_fit(model_dict, list_free_params, data;
    lb=lb, ub=ub, weights=[1.0, 0.0, 0.0], nboot=500, seed=42)

boot.median        # median of the bootstrap distribution
boot.sigma         # (84th − 16th percentile) / 2
boot.sigma_plus    # asymmetric error bars
boot.samples       # nboot × npar matrix, for corner plots or derived quantities
boot.covar         # parameter covariance

The resampling unit is set by granularity:

• :config (default) — one block per (MJD, baseline / triangle / telescope); all wavelength channels of a block are kept or dropped together. This is PMOIRED's "spectral vector".
• :epoch — one block per MJD. More conservative, but it needs many epochs: on a single night it becomes erratic.
• :point — one block per data point. Destroys the correlation structure and will underestimate the uncertainties on real data; provided for comparison.

mode selects how block multiplicities are drawn:

"Measured ratio" is the quoted σ divided by the true parameter scatter over 120 simulated realisations (1.00 = calibrated); see the Uncertainties page for the study, the wall times and the comparison with PMOIRED. :halfsample is the best-calibrated option when blocks are numerous, and the cheapest; with few blocks it becomes conservative. :replacement is the most uniform across regimes, which is why it is the default. :pmoired is not an estimator to quote from: it reproduces PMOIRED's construction so that its results can be cross-checked, and it draws multiplicities with variance ½ instead of 1, which makes its error bars about √2 too small. bootstrap_fit warns the first time it is used.

Derived quantities are obtained from the raw samples, which propagates the full covariance:

sep = sqrt.(boot.samples[boot.mask, 1].^2 .+ boot.samples[boot.mask, 2].^2)
println("separation = ", median(sep), " ± ", 0.5*(quantile(sep, 0.84) - quantile(sep, 0.16)))

What about adding noise to the data?

perturb_data(data) displaces every observable by a Gaussian draw from its error bar. It is useful for building simulated datasets, but refitting such replicates is not a bootstrap: it takes the quoted errors at face value, so it can only reproduce the analytic covariance — and in the validation study it understated the true scatter by a factor 11–25 once calibration systematics were present, and by exactly the factor by which the error bars were underestimated. resample_data is a deprecated alias for it.

See example_bootstrap_fit.jl for a side-by-side comparison of all schemes.