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 ffffffff-ffff-ffff-0000-004287710a13 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:1948Download 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:
vis_like = InterferometryLikelihood(
(; filename="Sim_data_2023_1_.oifits", epoch=mjd("2023-06-01"), band=:F480M, use_vis2=false),
(; filename="Sim_data_2023_2_.oifits", epoch=mjd("2023-08-15"), band=:F480M, use_vis2=false),
(; filename="Sim_data_2024_1_.oifits", epoch=mjd("2024-06-01"), band=:F480M, use_vis2=false),
)OctofitterInterferometry.InterferometryLikelihood Table with 14 columns and 3 rows:
filename epoch band use_vis2 u ⋯
┌────────────────────────────────────────────────────────────────────────
1 │ Sim_data_2023_1_.oi… 60096.0 F480M false [-3.46641e5; 6.7596… ⋯
2 │ Sim_data_2023_2_.oi… 60171.0 F480M false [-3.46641e5; 6.7596… ⋯
3 │ Sim_data_2024_1_.oi… 60462.0 F480M false [-3.46641e5; 6.7596… ⋯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 Visual{KepOrbit} begin
a ~ truncated(Normal(2,0.1), lower=0)
e ~ truncated(Normal(0, 0.05),lower=0, upper=1.0)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
# Our prior on the planet's photometry
# 0 +- 10% of stars brightness (assuming this is unit of data files)
F480M ~ truncated(Normal(0, 0.1),lower=0)
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,60171) # reference epoch for θ. Choose an MJD date near your data.
end
@system Tutoria begin
M ~ truncated(Normal(1.5, 0.01), lower=0)
plx ~ truncated(Normal(100., 0.1), lower=0)
end vis_like bSystem model Tutoria
Derived:
Priors:
M Truncated(Distributions.Normal{Float64}(μ=1.5, σ=0.01); lower=0.0)
plx Truncated(Distributions.Normal{Float64}(μ=100.0, σ=0.1); lower=0.0)
Planet b
Derived:
ω, Ω, θ, tp,
Priors:
a Truncated(Distributions.Normal{Float64}(μ=2.0, σ=0.1); lower=0.0)
e Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.05); lower=0.0, upper=1.0)
i Sine()
ω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)
F480M Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.1); lower=0.0)
θy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
θx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Octofitter.UnitLengthPrior{:ωy, :ωx}: √(ωy^2+ωx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:Ωy, :Ωx}: √(Ωy^2+Ωx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:θy, :θx}: √(θy^2+θx^2) ~ LogNormal(log(1), 0.02)
OctofitterInterferometry.InterferometryLikelihood Table with 14 columns and 3 rows:
filename epoch band use_vis2 u ⋯
┌────────────────────────────────────────────────────────────────────────
1 │ Sim_data_2023_1_.oi… 60096.0 F480M false [-3.46641e5; 6.7596… ⋯
2 │ Sim_data_2023_2_.oi… 60171.0 F480M false [-3.46641e5; 6.7596… ⋯
3 │ Sim_data_2024_1_.oi… 60462.0 F480M false [-3.46641e5; 6.7596… ⋯
Create the model object and run octofit_pigeons:
model = Octofitter.LogDensityModel(Tutoria)
using Pigeons
results,pt = octofit_pigeons(model, n_rounds=10);[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (185.17099342738285, 193.07781517261017)
┌ Warning: Invalid log likelihood encountered. (maxlog=1)
│ θ = (M = 5.059927634544915e24, plx = Inf, planets = (b = (a = 9.5267808190834e6, e = 0.08077808550150817, i = 2.1757296766378484, ωy = -1.896877841334627, ωx = -1.3190299232134697, Ωy = -1.5082081702986658, Ωx = 0.15697632703720482, F480M = 0.00016486720032773768, θy = 0.7712755317770947, θx = 1.3411209545116216, ω = -2.5339815537799977, Ω = 3.0378847248402727, θ = 1.0488886509192878, tp = 60167.57544379029),))
│ llike = NaN
│ θ_transformed =
│ 12-element Vector{Float64}:
│ 56.88339441358087
│ ⋮
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/logdensitymodel.jl:105
┌ Warning: Negative semi-major is required for hyperbolic (e>1) orbits. Flipping sign (maxlog=1).
└ @ PlanetOrbits ~/.julia/packages/PlanetOrbits/XOx9y/src/orbit-keplerian.jl:61
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.3 5.34 7.68 1.35e+09 -8.24e+04 0 0.721 0.937 0.969
4 0 4.25 5.15 1.55 1.48e+09 -1.59e+04 0 0.697 0.887 0.933
8 0 4.63 4.54 2.35 2.4e+09 168 2.51e-12 0.704 0.904 0.929
16 0 6.53 4.85 4.83 4.98e+09 177 0.125 0.633 0.901 0.921
32 0 6.77 6.13 9.45 9.71e+09 171 2.88e-07 0.584 0.901 0.919
64 4 6.99 2.7 20.1 1.97e+10 176 0.371 0.688 0.905 0.919
128 9 6.66 2.82 38.8 3.96e+10 176 0.406 0.694 0.91 0.922
256 29 6.95 2.6 77.6 7.91e+10 176 0.474 0.692 0.907 0.922
512 75 6.93 2.74 155 1.58e+11 176 0.483 0.688 0.912 0.921
1.02e+03 151 6.83 2.64 310 3.16e+11 176 0.481 0.694 0.91 0.921
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────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[:b_F480M][:], axis=(;xlabel="F480M"))
Determine the significance of the detection:
using Statistics
phot = results[:b_F480M][:]
snr = mean(phot)/std(phot)7.432261600000638Plot the resulting orbit:
octoplot(model, results)
Plot position at each epoch:
using PlanetOrbits
els = Octofitter.construct_elements(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
We can use PairPlots.jl to create a contour plot of positions at all three epochs:
els = Octofitter.construct_elements(results,:b,:);
fig = pairplot(
[
(;
ra=raoff.(els, epoch)[:],
dec=decoff.(els, epoch)[:],
)=>(PairPlots.Contourf(),)
for epoch in vis_like.table.epoch
]...,
bodyaxis=(;width=400,height=400),
axis=(;
ra=(;reversed=true, lims=(;low=250,high=-250,)),
dec=(;lims=(;low=-250,high=250,)),
),
labels=Dict(:ra=>"ra offset [mas]", :dec=>"dec offset [mas]"),
)
Makie.scatter!(fig.content[1], [0],[0],marker='⭐', markersize=30, color=:black)
fig
Finally we can examine the joint photometry and orbit posterior as a corner plot:
using PairPlots
using CairoMakie: Makie
octocorner(model, results)