Fitting Images
One of the key features of Octofitter.jl is the ability to search for planets directly from images of the system. Sampling from images is much more computationally demanding than sampling from astrometry, but it allows for a few very powerful results:
- You can search for a planet that is not well detected in a single image
By this, we mean you can feed in images of a system with no clear detections, and see if a planet is hiding in the noise based off of its Kepelerian motion.
- Not detecting a planet in a given image can be almost as useful as a detection for constraining its orbit.
If you have a clear detection in one epoch, but no detection in another, Octofitter can use the image from the second epoch to rule out large swathes of possible orbits.
Sampling from images can be freely combined with any known astrometry points, as well as astrometric acceleration. See advanced models for more details.
Image modelling is supported in Octofitter via the extension package OctofitterImages. To install it, run pkg> add http://github.com/sefffal/Octofitter.jl:OctofitterImages
Preparing images
The first step will be to load your images. For this, we will use our AstroImages.jl package.
Start by loading your images:
using Octofitter
using OctofitterImages
using Distributions
using Pigeons
using AstroImages
using CairoMakie
# Load individual iamges
# image1 = load("image1.fits")
# image2 = load("image2.fits")
# Or slices from a cube:
# cube = load("cube1.fits")
# image1 = cube[:,:,1]
# Download sample images from GitHub
download(
"https://github.com/sefffal/Octofitter.jl/raw/main/docs/image-examples-1.fits",
"image-examples-1.fits"
)
# Or multi-extension FITS (this example)
images = AstroImages.load("image-examples-1.fits",:)5-element Vector{AstroImages.AstroImageMat{Float64, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Tuple{}, Matrix{Float64}, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}}}:
[-0.340230578841865 -0.31239887526186355 … -0.2903753168231458 -0.34686827004814064; -0.29826568928566355 -0.27872689920999305 … -0.246452640103626 -0.2873140658333837; … ; 0.2673577868138141 0.2674003307584391 … -0.35696041486602914 -0.4320582965232348; 0.34147717271077094 0.3420908537177248 … -0.3358061981896538 -0.4139117427343716]
[-1.0477769487233675 -0.903332487139985 … 1.4299602191511185 1.5070937596027922; -0.9888724404777535 -0.8589483767137767 … 1.1201264194935676 1.1803510698319304; … ; -0.43042419506343155 -0.46545013356425746 … -0.7066171674807686 -0.8619933383603304; -0.5837006039738425 -0.6021768042432977 … -0.8745923304009366 -1.0622653729194567]
[-1.3673176596472436 -1.1985348376204175 … 1.1486692344915643 1.2929199629807906; -1.1143572848991925 -0.9810874969329584 … 0.9739105840286464 1.079311141942771; … ; 1.0292715351056276 0.8022139684685823 … 0.674123592556085 0.8519386886571362; 1.2170463982677187 0.9582307659215409 … 0.7602550494569255 0.959778160219427]
[-0.16664741987337547 -0.13437529569419648 … -1.2963013699793633 -1.413565369547865; -0.2149990713151032 -0.1847368741065949 … -1.0549153738841368 -1.1490484959706464; … ; -0.5872462145987273 -0.5703179209210852 … 0.19469039731388912 0.31932314143450863; -0.7173113967617594 -0.6848676990867931 … 0.19122281859925802 0.3434506497145075]
[0.47470453920043826 0.4747486393369815 … -0.16662161045865695 -0.16297895606178783; 0.4078627130832368 0.4041383409077344 … -0.18122203442704718 -0.17883512905931084; … ; 0.07003384675445126 0.19916060264117244 … -0.44137489323835993 -0.5253547777211444; 0.023766763613137065 0.1809119688664255 … -0.5699244785815155 -0.6776496918698217]You can preview the image using imshow2 from AstroImages:
# imshow2(image1, cmap=:magma) # for a single image
hcat(imview.(images, clims=(-1.0, 4.0))...)
Your images should either be convolved with a gaussian of diameter one λ/D, or be matched filtered. This is so that the values of the pixels in the image represent the photometry at that location.
If you want to perform the convolution in Julia, see ImageFiltering.jl.
Build the model
First, we create a table of our image data that will be attached to the Planet:
image_data = ImageLikelihood(
(band=:H, image=AstroImages.recenter(images[1]), platescale=10.0, epoch=1238.6),
(band=:H, image=AstroImages.recenter(images[2]), platescale=10.0, epoch=1584.7),
(band=:H, image=AstroImages.recenter(images[3]), platescale=10.0, epoch=3220.0),
(band=:H, image=AstroImages.recenter(images[4]), platescale=10.0, epoch=7495.9),
(band=:H, image=AstroImages.recenter(images[5]), platescale=10.0, epoch=7610.4),
)OctofitterImages.ImageLikelihood Table with 6 columns and 5 rows:
band image platescale epoch contrast ⋯
┌────────────────────────────────────────────────────────────────────────
1 │ H [-0.340231 -0.31239… 10.0 1238.6 65-element extrapol… ⋯
2 │ H [-1.04778 -0.903332… 10.0 1584.7 65-element extrapol… ⋯
3 │ H [-1.36732 -1.19853 … 10.0 3220.0 65-element extrapol… ⋯
4 │ H [-0.166647 -0.13437… 10.0 7495.9 65-element extrapol… ⋯
5 │ H [0.474705 0.474749 … 10.0 7610.4 65-element extrapol… ⋯Provide one entry for each image you want to sample from. Ensure that each image has been re-centered so that index [0,0] is the position of the star. Areas of the image where there is no data should be filled with NaN and will not contribute to the likelihood of your model. platescale should be the pixel scale of your images, in milliarseconds / pixel. epoch should be the Modified Julian Day (MJD) that your image was taken. You can use the mjd("2021-09-09") function to calculate this for you. band should be a symbol that matches the name you supplied when you created the Planet.
By default, the contrast of the images is calculated automatically, but you can supply your own contrast curve as well by also passing contrast=OctofitterImages.contrast_interp(AstroImages.recenter(my_image)).
You can freely mix and match images from different instruments as long as you specify the correct platescale. You can also provide images from multiple bands and they will be sampled independently. If you wish to tie them together, see Connecting Mass with Photometry.
Now specify the planet:
@planet b Visual{KepOrbit} begin
H ~ Normal(3.8, 0.5)
a ~ truncated(Normal(13, 4), lower=0, upper=100)
e ~ Uniform(0.0, 0.5)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,1238.6)
end image_dataPlanet b
Derived:
ω, Ω, θ, tp,
Priors:
H Distributions.Normal{Float64}(μ=3.8, σ=0.5)
a Truncated(Distributions.Normal{Float64}(μ=13.0, σ=4.0); lower=0.0, upper=100.0)
e Distributions.Uniform{Float64}(a=0.0, b=0.5)
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)
θ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)
OctofitterImages.ImageLikelihood Table with 6 columns and 5 rows:
band image platescale epoch contrast ⋯
┌────────────────────────────────────────────────────────────────────────
1 │ H [-0.340231 -0.31239… 10.0 1238.6 65-element extrapol… ⋯
2 │ H [-1.04778 -0.903332… 10.0 1584.7 65-element extrapol… ⋯
3 │ H [-1.36732 -1.19853 … 10.0 3220.0 65-element extrapol… ⋯
4 │ H [-0.166647 -0.13437… 10.0 7495.9 65-element extrapol… ⋯
5 │ H [0.474705 0.474749 … 10.0 7610.4 65-element extrapol… ⋯
Note how we also provided a prior on the photometry called H. We can put any name we want here, as long as it's used consistently throughout the model specification.
See Fit PlanetRelAstromLikelihood for a description of the different orbital parameters, and conventions used.
Finally, create the system and pass in the planet.
@system HD82134 begin
M ~ truncated(Normal(2.0, 0.1),lower=0)
plx ~ truncated(Normal(45., 0.02),lower=0)
end bSystem model HD82134
Derived:
Priors:
M Truncated(Distributions.Normal{Float64}(μ=2.0, σ=0.1); lower=0.0)
plx Truncated(Distributions.Normal{Float64}(μ=45.0, σ=0.02); lower=0.0)
Planet b
Derived:
ω, Ω, θ, tp,
Priors:
H Distributions.Normal{Float64}(μ=3.8, σ=0.5)
a Truncated(Distributions.Normal{Float64}(μ=13.0, σ=4.0); lower=0.0, upper=100.0)
e Distributions.Uniform{Float64}(a=0.0, b=0.5)
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)
θ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)
OctofitterImages.ImageLikelihood Table with 6 columns and 5 rows:
band image platescale epoch contrast ⋯
┌────────────────────────────────────────────────────────────────────────
1 │ H [-0.340231 -0.31239… 10.0 1238.6 65-element extrapol… ⋯
2 │ H [-1.04778 -0.903332… 10.0 1584.7 65-element extrapol… ⋯
3 │ H [-1.36732 -1.19853 … 10.0 3220.0 65-element extrapol… ⋯
4 │ H [-0.166647 -0.13437… 10.0 7495.9 65-element extrapol… ⋯
5 │ H [0.474705 0.474749 … 10.0 7610.4 65-element extrapol… ⋯
If you want to search for two or more planets in the same images, just create multiple Planets and pass the same images to each. You'll need to adjust the priors in some way to prevent overlap.
You can also do some very clever things like searching for planets that are co-planar and/or have a specific resonance between their periods. To do this, put the planet of the system or base period as variables of the system and derive the planet variables from those values of the system.
Sampling
Sampling from images is much more challenging than relative astrometry or proper motion anomaly, so the fitting process tends to take longer.
This is because the posterior is much "bumpier" with images. One way this manifests is very high tree depths. You might see a sampling report that says max_tree_depth_frac = 0.9 or even 1.0. To encourage the sampler to take larger steps and explore the images, it's recommended to lower the target acceptance ratio to around 0.5±0.2 and also increase the number of adapataion steps.
model = Octofitter.LogDensityModel(HD82134)
chain, pt = octofit_pigeons(model, n_rounds=10)
display(chain)[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (49.78181843537236, 55.319460409549635)
┌ Warning: Invalid log likelihood encountered. (maxlog=1)
│ θ = (M = 6.0578691279268965e28, plx = 0.0, planets = (b = (H = 0.9439541744647544, a = 1.8893994068881347e-99, e = 0.0011521166373721337, i = 3.1386143749216426, ωy = 0.2252495237299006, ωx = -0.41933411167820633, Ωy = 60.911133226193314, Ωx = -18.016295245498792, θy = 123.80148706405392, θx = -53.88100344060731, ω = -1.0778644740780776, Ω = -0.287580756257694, θ = -0.4104959520029153, tp = 1238.6),))
│ llike = NaN
│ θ_transformed =
│ 12-element Vector{Float64}:
│ 66.27374071302354
│ ⋮
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/logdensitymodel.jl:105
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 2.45 4.59 5.74 3.55e+08 31.8 3.79e-10 0.773 0.945 0.965
4 0 4.54 4.45 0.614 9.11e+06 43.7 0.00288 0.71 0.9 0.924
8 0 4.64 4.51 0.753 1.04e+06 49 0.0997 0.705 0.902 0.922
16 0 4.98 5.43 1.62 1.93e+06 51.8 0.0842 0.664 0.91 0.925
32 0 5.38 5.45 3.17 3.7e+06 58.7 0.0814 0.651 0.894 0.914
64 2 5.77 3.18 7.41 2.42e+08 64.3 0.0317 0.711 0.904 0.919
128 8 5.35 3.15 13.3 4.38e+08 64.7 0.0838 0.726 0.901 0.917
256 15 5.11 3.96 27.1 8.93e+08 64.7 0.131 0.707 0.9 0.918
512 36 5.16 4.05 53 1.75e+09 64.5 0.257 0.703 0.905 0.923
1.02e+03 56 5.28 4.62 105 3.48e+09 65.2 0.538 0.681 0.905 0.921
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────octofit_pigeons scales very well across multiple cores. Start julia with julia --threads=auto to make sure you have multiple threads available for sampling.
Diagnostics
The first thing you should do with your results is check a few diagnostics to make sure the sampler converged as intended.
The acceptance rate should be somewhat lower than when fitting just astrometry, e.g. around the 0.6 target.
You can make a trace plot:
lines(
chain["b_a"][:],
axis=(;
xlabel="iteration",
ylabel="semi-major axis (aU)"
)
)
And an auto-correlation plot:
using StatsBase
lines(
autocor(chain["b_e"][:], 1:500),
axis=(;
xlabel="lag",
ylabel="autocorrelation",
)
)
For this model, there is somewhat higher correlation between samples. Some thinning to remove this correlation is recommended.
Analysis
We can now view the orbit fit:
fig = octoplot(model, chain)
With a bit of work, we can plot one of the images under the orbit.
fig = octoplot(model, chain)
ax = fig.content[1] # grap first axis in the figure
# We have to do some annoying work to get the image orientated correctly
# since we want the RA axis increasing to the left.
image_idx = 2
platescale = image_data.table.platescale[image_idx]
img = AstroImages.recenter(AstroImage(collect(image_data.table.image[image_idx])[end:-1:begin,:]))
imgax1 = dims(img,1) .* platescale
imgax2 = dims(img,2) .* platescale
h = heatmap!(ax, imgax1, imgax2, collect(img), colormap=:greys)
Makie.translate!(h, 0,0,-1) # Send heatmap to back of the plot
# Add colorbar for image
Colorbar(fig[1,2], h, label="image flux")
Makie.resize_to_layout!(fig)
fig
Another useful view would be the orbits over a stack of the maximum pixel values of all images.
fig = octoplot(model, chain)
ax = fig.content[1] # grap first axis in the figure
# We have to do some annoying work to get the image orientated correctly
# since we want the RA axis increasing to the left.
platescale = image_data.table.platescale[image_idx]
imgs = maximum(stack(image_data.table.image),dims=3)[:,:]
img = AstroImages.recenter(AstroImage(imgs[end:-1:begin,:]))
imgax1 = dims(img,1) .* platescale
imgax2 = dims(img,2) .* platescale
h = heatmap!(ax, imgax1, imgax2, collect(img), colormap=:greys)
Makie.translate!(h, 0,0,-1) # Send heatmap to back of the plot
Makie.resize_to_layout!(fig)
fig
Pair Plot
We can show the relationships between variables on a pair plot (aka corner plot):
using CairoMakie, PairPlots
octocorner(model, chain, small=true)
Note that this time, we also show the recovered photometry in the corner plot.
Assessing Detections
To assess a detection, we can treat all the orbital variables as nuisance parameters. We start by plotting the marginal distribution of the flux parameter, H:
hist(chain["b_H"][:], axis=(xlabel="H", ylabel="counts"))
We can calculate an analog of the traditional signal to noise ratio (SNR) using that same histogram:
flux = chain["b_H"]
snr = mean(flux)/std(flux)14.132725219574658It might be better to consider a related measure, like the median flux over the interquartile distance. This will depend on your application.