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.

Note

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-8b2a-b38e-0000-011b7844e76c 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:2018

Download 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"), spectrum_var=:contrast_F480M, use_vis2=false),
    (; filename="Sim_data_2023_2_.oifits", epoch=mjd("2023-08-15"), spectrum_var=:contrast_F480M, use_vis2=false),
    (; filename="Sim_data_2024_1_.oifits", epoch=mjd("2024-06-01"), spectrum_var=:contrast_F480M, use_vis2=false),
)
OctofitterInterferometry.InterferometryLikelihood Table with 14 columns and 3 rows:
     filename              epoch    spectrum_var    use_vis2  ⋯
   ┌───────────────────────────────────────────────────────────
 1 │ Sim_data_2023_1_.oi…  60096.0  contrast_F480M  false     ⋯
 2 │ Sim_data_2023_2_.oi…  60171.0  contrast_F480M  false     ⋯
 3 │ Sim_data_2024_1_.oi…  60462.0  contrast_F480M  false     ⋯

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
Example block output
@planet b Visual{KepOrbit} begin
    a ~ truncated(Normal(2,0.1), lower=0.1)
    e ~ truncated(Normal(0, 0.05),lower=0, upper=0.90)
    i ~ Sine()
    ω ~ UniformCircular()
    Ω ~ UniformCircular()

    # Our prior on the planet's photometry
    # 0 +- 10% of stars brightness (assuming this is unit of data files)
    contrast_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.1)
    plx ~ truncated(Normal(100., 0.1), lower=0.1)
end vis_like b
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:
  ω, Ω, θ, 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()
  ω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)
  contrast_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    spectrum_var    use_vis2  ⋯
   ┌───────────────────────────────────────────────────────────
 1 │ Sim_data_2023_1_.oi…  60096.0  contrast_F480M  false     ⋯
 2 │ Sim_data_2023_2_.oi…  60171.0  contrast_F480M  false     ⋯
 3 │ Sim_data_2024_1_.oi…  60462.0  contrast_F480M  false     ⋯

Create the model object and run octofit_pigeons:

model = Octofitter.LogDensityModel(Tutoria)

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.000029 seconds (15 allocations: 3.984 KiB)
∇ℓπcallback(θ): 0.000068 seconds (16 allocations: 26.516 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 parameters held fixed).
┌ Info: Found sample of initial positions
│   logpost_range = (185.03141950155475, 195.90911201003175)
└   mean_logpost = 192.56230385389688
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
        2          0       3.61       4.07      0.368   8.83e+07  -1.77e+04          0      0.752          1          1
        4          0       4.83        3.9       0.49   1.63e+08        164          0      0.718          1          1
        8          0       4.82       4.67      0.852   3.17e+08        180   4.62e-82      0.694          1          1
       16          0       5.12       6.14        1.8   6.37e+08        178      0.125      0.637          1          1
       32          0       5.08       6.21       3.57   1.27e+09        179     0.0625      0.636          1          1
       64          5       5.45       2.59       7.34   2.61e+09        176     0.0768      0.741          1          1
      128         10        5.8       2.61       14.4   5.27e+09        176      0.114      0.729          1          1
      256         33       5.91       2.67       28.7   1.05e+10        176      0.295      0.723          1          1
      512         74       6.06       2.57       57.2    2.1e+10        176      0.503      0.722          1          1
 1.02e+03        158       6.12       2.68        115   4.19e+10        176      0.542      0.716          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[:b_contrast_F480M][:], axis=(;xlabel="F480M"))
Example block output

Determine the significance of the detection:

using Statistics
phot = results[:b_contrast_F480M][:]
snr = mean(phot)/std(phot)
7.3506780473875715

Plot the resulting orbit:

octoplot(model, results)
Example block output

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
Example block output

Finally we can examine the joint photometry and orbit posterior as a corner plot:

using PairPlots
using CairoMakie: Makie
octocorner(model, results)
Example block output