Fitting Interferometric Observables
In this tutorial, we fit a planet & orbit model to a sequence of interferometric observations. Closure phases and squared visibilities are supported.
We load the observations in OI-FITS format and model them as a point source orbiting a star.
Interferometer modelling is supported in Octofitter via the extension package OctofitterInterferometry. To install it, run pkg> add http://github.com/sefffal/Octofitter.jl:OctofitterInterferometry
using Octofitter
using OctofitterInterferometry
using Distributions
using CairoMakie
using PairPlots┌ Warning: Module Octofitter with build ID fafbfcfd-73cd-724a-0000-001de4318a31 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2018Download simulated JWST AMI observations from our examples folder on GitHub:
download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2023_1_.oifits", "Sim_data_2023_1_.oifits")
download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2023_2_.oifits", "Sim_data_2023_2_.oifits")
download("https://github.com/sefffal/Octofitter.jl/raw/main/examples/AMI_data/Sim_data_2024_1_.oifits", "Sim_data_2024_1_.oifits")"Sim_data_2024_1_.oifits"Create the likelihood object:
data = Table([
(; filename="Sim_data_2023_1_.oifits", epoch=mjd("2023-06-01"), use_vis2=false),
(; filename="Sim_data_2023_2_.oifits", epoch=mjd("2023-08-15"), use_vis2=false),
(; filename="Sim_data_2024_1_.oifits", epoch=mjd("2024-06-01"), use_vis2=false),
])
vis_like = InterferometryLikelihood(
data,
name="NIRISS-AMI",
variables=@variables begin
# For single planet:
flux ~ truncated(Normal(0, 0.1), lower=0) # Planet flux/contrast (array with one element)
# For multiple planets (array - one per planet):
# flux ~ Product([truncated(Normal(0, 0.1), lower=0), truncated(Normal(0, 0.1), lower=0)])
# Optional calibration parameters:
platescale = 1.0 # Platescale multiplier [could use: platescale ~ truncated(Normal(1, 0.01), lower=0)]
northangle = 0.0 # North angle offset in radians [could use: northangle ~ Normal(0, deg2rad(1))]
σ_cp_jitter = 0.0 # closure phase jitter [could use ~ LogUniform(0.1, 100))
end
)OctofitterInterferometry.InterferometryLikelihood Table with 13 columns and 3 rows:
filename epoch use_vis2 u ⋯
┌─────────────────────────────────────────────────────────────────
1 │ Sim_data_2023_1_.oi… 60096.0 false [-3.46641e5; 6.7596… ⋯
2 │ Sim_data_2023_2_.oi… 60171.0 false [-3.46641e5; 6.7596… ⋯
3 │ Sim_data_2024_1_.oi… 60462.0 false [-3.46641e5; 6.7596… ⋯If you want to include multiple bands, group these into different InterferometryLikelihood objects with different instrument names (i.e. include the band in the name for the sake of bookkeeping)
Plot the closure phases:
fig = Makie.Figure()
ax = Axis(
fig[1,1],
xlabel="index",
ylabel="closure phase",
)
Makie.stem!(
vis_like.table.cps_data[1][:],
label="epoch 1",
)
Makie.stem!(
vis_like.table.cps_data[2][:],
label="epoch 2"
)
Makie.stem!(
vis_like.table.cps_data[3][:],
label="epoch 3"
)
Makie.Legend(fig[1,2], ax)
fig
planet_b = Planet(
name="b",
basis=Visual{KepOrbit},
likelihoods=[],
variables=@variables begin
M = system.M
a ~ truncated(Normal(2,0.1), lower=0.1)
e ~ truncated(Normal(0, 0.05),lower=0, upper=0.90)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(θ, 60171; M, e, a, i, ω, Ω) # reference epoch for θ. Choose an MJD date near your data.
end
)
sys = System(
name="Tutoria",
companions=[planet_b],
likelihoods=[vis_like],
variables=@variables begin
M ~ truncated(Normal(1.5, 0.01), lower=0.1)
plx ~ truncated(Normal(100., 0.1), lower=0.1)
end
)System model Tutoria
Derived:
Priors:
M Truncated(Distributions.Normal{Float64}(μ=1.5, σ=0.01); lower=0.1)
plx Truncated(Distributions.Normal{Float64}(μ=100.0, σ=0.1); lower=0.1)
Planet b
Derived:
ω, Ω, θ, M, tp,
Priors:
a Truncated(Distributions.Normal{Float64}(μ=2.0, σ=0.1); lower=0.1)
e Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.05); lower=0.0, upper=0.9)
i Sine()
ωx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
ωy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Ωx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Ωy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
θx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
θy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Octofitter.UnitLengthPrior{:ωx, :ωy}: √(ωx^2+ωy^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:Ωx, :Ωy}: √(Ωx^2+Ωy^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:θx, :θy}: √(θx^2+θy^2) ~ LogNormal(log(1), 0.02)
OctofitterInterferometry.InterferometryLikelihood Table with 13 columns and 3 rows:
filename epoch use_vis2 u ⋯
┌─────────────────────────────────────────────────────────────────
1 │ Sim_data_2023_1_.oi… 60096.0 false [-3.46641e5; 6.7596… ⋯
2 │ Sim_data_2023_2_.oi… 60171.0 false [-3.46641e5; 6.7596… ⋯
3 │ Sim_data_2024_1_.oi… 60462.0 false [-3.46641e5; 6.7596… ⋯
Create the model object and run octofit_pigeons:
model = Octofitter.LogDensityModel(sys)
using Pigeons
results,pt = octofit_pigeons(model, n_rounds=10);[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 12
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000031 seconds (12 allocations: 3.797 KiB)
∇ℓπcallback(θ): 0.000081 seconds (13 allocations: 26.047 KiB)
┌ Warning: This model has priors that cannot be sampled IID.
└ @ OctofitterPigeonsExt ~/work/Octofitter.jl/Octofitter.jl/ext/OctofitterPigeonsExt/OctofitterPigeonsExt.jl:95
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: Performing global optimization with 12 parameters (0 initial parameter value provided).
┌ Warning: Unrecognized stop reason: Too many steps (101) without any function evaluations (probably search has converged). Defaulting to ReturnCode.Default.
└ @ Optimization ~/.julia/packages/Optimization/hfEs4/src/utils.jl:93
┌ Info: Found sample of initial positions
│ logpost_range = (184.8741324911994, 195.99828422111315)
└ mean_logpost = 192.66587769078745
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.06 3.9 0.444 8.57e+07 -2.69e+03 0 0.776 1 1
4 0 5.42 4.01 0.396 1.48e+08 -2.89e+04 0 0.696 1 1
8 0 7.04 5.22 0.858 2.95e+08 28.7 4.49e-132 0.605 1 1
16 0 6.26 6.13 1.66 5.79e+08 178 0.0067 0.6 1 1
32 0 6.67 6.27 3.31 1.15e+09 176 0.00439 0.583 1 1
64 4 6.64 2.93 7.02 2.43e+09 176 0.21 0.691 1 1
128 10 6.79 2.43 13.8 4.91e+09 176 0.318 0.703 1 1
256 34 6.84 2.7 27.5 9.79e+09 176 0.418 0.692 1 1
512 74 6.79 2.84 54.6 1.95e+10 176 0.435 0.689 1 1
1.02e+03 157 6.79 2.62 109 3.9e+10 176 0.478 0.696 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Note that we use Pigeons paralell tempered sampling (octofit_pigeons) instead of HMC (octofit) because interferometry data is almost always multi-modal (or more precisely non-convex, there is often still a single mode that dominates).
Examine the recovered photometry posterior:
hist(results[:NIRISS_AMI_flux][:], axis=(;xlabel="flux"))
Determine the significance of the detection:
using Statistics
phot = results[:NIRISS_AMI_flux][:]
snr = mean(phot)/std(phot)7.306387749146333Plot the resulting orbit:
octoplot(model, results)
Plot only the position at each epoch:
using PlanetOrbits
els = Octofitter.construct_elements(model, results,:b,:);
fig = Makie.Figure()
ax = Makie.Axis(
fig[1,1],
autolimitaspect = 1,
xreversed=true,
xlabel="ΔR.A. (mas)",
ylabel="ΔDec. (mas)",
)
for epoch in vis_like.table.epoch
Makie.scatter!(
ax,
raoff.(els, epoch)[:],
decoff.(els, epoch)[:],
label=string(mjd2date(epoch)),
markersize=1.5,
)
end
Makie.Legend(fig[1,2], ax, "date")
fig
Finally we can examine the joint photometry and orbit posterior as a corner plot:
using PairPlots
using CairoMakie: Makie
octocorner(model, results)