Stokes I Simultaneous Image and Instrument Modeling
In this tutorial, we will create a preliminary reconstruction of the 2017 M87 data on April 6 by simultaneously creating an image and model for the instrument. By instrument model, we mean something akin to self-calibration in traditional VLBI imaging terminology. However, unlike traditional self-cal, we will solve for the gains each time we update the image self-consistently. This allows us to model the correlations between gains and the image.
To get started we load Comrade.
julia
using Comrade
using Pyehtim
using LinearAlgebra CondaPkg Found dependencies: /home/runner/.julia/packages/CondaPkg/8GjrP/CondaPkg.toml
CondaPkg Found dependencies: /home/runner/.julia/packages/PythonCall/83z4q/CondaPkg.toml
CondaPkg Found dependencies: /home/runner/.julia/packages/Pyehtim/bQtHC/CondaPkg.toml
CondaPkg Resolving changes
+ ehtim (pip)
+ libstdcxx
+ libstdcxx-ng
+ numpy
+ numpy (pip)
+ openssl
+ pandas
+ python
+ setuptools (pip)
+ uv
CondaPkg Initialising pixi
│ /home/runner/.julia/artifacts/cefba4912c2b400756d043a2563ef77a0088866b/bin/pixi
│ init
│ --format pixi
└ /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg
✔ Created /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/pixi.toml
CondaPkg Wrote /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/pixi.toml
│ [dependencies]
│ openssl = ">=3, <3.6"
│ libstdcxx = ">=3.4,<15.0"
│ uv = ">=0.4"
│ libstdcxx-ng = ">=3.4,<15.0"
│ pandas = "<2"
│ numpy = ">=1.24, <2.0"
│
│ [dependencies.python]
│ channel = "conda-forge"
│ build = "*cp*"
│ version = ">=3.10,!=3.14.0,!=3.14.1,<4, >=3.6,<=3.12"
│
│ [project]
│ name = ".CondaPkg"
│ platforms = ["linux-64"]
│ channels = ["conda-forge"]
│ channel-priority = "strict"
│ description = "automatically generated by CondaPkg.jl"
│
│ [pypi-dependencies]
│ ehtim = ">=1.2.10, <2.0"
│ numpy = ">=1.24, <2.0"
└ setuptools = "*"
CondaPkg Installing packages
│ /home/runner/.julia/artifacts/cefba4912c2b400756d043a2563ef77a0088866b/bin/pixi
│ install
└ --manifest-path /home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/pixi.toml
✔ The default environment has been installed.For reproducibility we use a stable random number genreator
julia
using StableRNGs
rng = StableRNG(42)StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)Load the Data
To download the data visit https://doi.org/10.25739/g85n-f134 First we will load our data:
julia
obs = ehtim.obsdata.load_uvfits(joinpath(__DIR, "..", "..", "Data", "SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))Python: <ehtim.obsdata.Obsdata object at 0x7f359acdf9d0>Now we do some minor preprocessing:
Scan average the data since the data have been preprocessed so that the gain phases coherent.
Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.
julia
obs = scan_average(obs).add_fractional_noise(0.02)Python: <ehtim.obsdata.Obsdata object at 0x7f35a5bf8150>Now we extract our complex visibilities.
julia
dvis = extract_table(obs, Visibilities())EHTObservationTable{Comrade.EHTVisibilityDatum{:I}}
source: M87
mjd: 57849
bandwidth: 1.856e9
sites: [:AA, :AP, :AZ, :JC, :LM, :PV, :SM]
nsamples: 274##Building the Model/Posterior
Now, we must build our intensity/visibility model. That is, the model that takes in a named tuple of parameters and perhaps some metadata required to construct the model. For our model, we will use a raster or ContinuousImage for our image model. The model is given below:
The model construction is very similar to Imaging a Black Hole using only Closure Quantities, except we include a large scale gaussian since we want to model the zero baselines. For more information about the image model please read the closure-only example. We use the @sky macro to bundle the model and its priors in one block. Each name ~ dist line adds an entry to the prior; everything else is the model body.
julia
using VLBIImagePriors
using Distributions
@sky function sky(grid; ftot, mimg, cprior)
c ~ cprior
σimg ~ truncated(Normal(0.0, 0.5); lower = 0.0)
fg ~ Uniform(0.0, 1.0)
# Apply the GMRF fluctuations to the image
rast = apply_fluctuations(CenteredLR(), mimg, σimg .* c.params)
pimg = parent(rast)
@. pimg = (ftot * (1 - fg)) * pimg
m = ContinuousImage(rast, BSplinePulse{3}())
# Add a large-scale gaussian to deal with the over-resolved mas flux
g = modify(Gaussian(), Stretch(μas2rad(500.0), μas2rad(500.0)), Renormalize(ftot * fg))
x, y = centroid(m)
return shifted(m, -x, -y) + g
endsky (generic function with 1 method)Now, let's set up our image model. The EHT's nominal resolution is 20-25 μas. Additionally, the EHT is not very sensitive to a larger field of view. Typically 60-80 μas is enough to describe the compact flux of M87. Given this, we only need to use a small number of pixels to describe our image.
julia
npix = 48
fovx = μas2rad(200.0)
fovy = μas2rad(200.0)9.69627362219072e-10Now let's form our cache's. First, we have our usual image cache which is needed to numerically compute the visibilities.
julia
grid = imagepixels(fovx, fovy, npix, npix)RectiGrid(
executor: ComradeBase.Serial()
Dimensions:
(↓ X Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points)
)Since we are using a Gaussian Markov random field prior we need to first specify our mean image. This behaves somewhat similary to a entropy regularizer in that it will start with an initial guess for the image structure. For this tutorial we will use a a symmetric Gaussian with a FWHM of 50 μas
julia
fwhmfac = 2 * sqrt(2 * log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(60.0) ./ fwhmfac))
mimg_raw = intensitymap(mpr, grid)
mimg = mimg_raw ./ flux(mimg_raw)┌ 48×48 IntensityMap{Float64, 2} ┐
├────────────────────────────────┴─────────────────────────────────────── dims ┐
↓ X Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points
└──────────────────────────────────────────────────────────────────────────────┘
↓ → -4.74713e-10 -4.54513e-10 -4.34312e-10 -4.14112e-10 -3.93911e-10 -3.73711e-10 -3.5351e-10 -3.33309e-10 -3.13109e-10 -2.92908e-10 -2.72708e-10 -2.52507e-10 -2.32307e-10 -2.12106e-10 -1.91905e-10 -1.71705e-10 -1.51504e-10 -1.31304e-10 -1.11103e-10 -9.09026e-11 -7.0702e-11 -5.05014e-11 -3.03009e-11 -1.01003e-11 1.01003e-11 3.03009e-11 5.05014e-11 7.0702e-11 9.09026e-11 1.11103e-10 1.31304e-10 1.51504e-10 1.71705e-10 1.91905e-10 2.12106e-10 2.32307e-10 2.52507e-10 2.72708e-10 2.92908e-10 3.13109e-10 3.33309e-10 3.5351e-10 3.73711e-10 3.93911e-10 4.14112e-10 4.34312e-10 4.54513e-10 4.74713e-10
-4.74713e-10 1.64194e-9 3.03721e-9 5.46989e-9 9.59111e-9 1.63736e-8 2.72149e-8 4.40409e-8 6.9389e-8 1.06442e-7 1.58971e-7 2.3116e-7 3.27259e-7 4.51084e-7 6.05354e-7 7.90948e-7 1.00617e-6 1.24618e-6 1.50272e-6 1.76426e-6 2.01665e-6 2.24433e-6 2.4318e-6 2.5654e-6 2.63493e-6 2.63493e-6 2.5654e-6 2.4318e-6 2.24433e-6 2.01665e-6 1.76426e-6 1.50272e-6 1.24618e-6 1.00617e-6 7.90948e-7 6.05354e-7 4.51084e-7 3.27259e-7 2.3116e-7 1.58971e-7 1.06442e-7 6.9389e-8 4.40409e-8 2.72149e-8 1.63736e-8 9.59111e-9 5.46989e-9 3.03721e-9 1.64194e-9
-4.54513e-10 3.03721e-9 5.61814e-9 1.0118e-8 1.77413e-8 3.02874e-8 5.03414e-8 8.14655e-8 1.28354e-7 1.96893e-7 2.9406e-7 4.27592e-7 6.05354e-7 8.34402e-7 1.11977e-6 1.46307e-6 1.86119e-6 2.30516e-6 2.77969e-6 3.26347e-6 3.73034e-6 4.15149e-6 4.49827e-6 4.7454e-6 4.87402e-6 4.87402e-6 4.7454e-6 4.49827e-6 4.15149e-6 3.73034e-6 3.26347e-6 2.77969e-6 2.30516e-6 1.86119e-6 1.46307e-6 1.11977e-6 8.34402e-7 6.05354e-7 4.27592e-7 2.9406e-7 1.96893e-7 1.28354e-7 8.14655e-8 5.03414e-8 3.02874e-8 1.77413e-8 1.0118e-8 5.61814e-9 3.03721e-9
-4.34312e-10 5.46989e-9 1.0118e-8 1.82222e-8 3.19514e-8 5.45465e-8 9.06628e-8 1.46716e-7 2.3116e-7 3.54596e-7 5.29591e-7 7.70077e-7 1.09022e-6 1.50272e-6 2.01665e-6 2.63493e-6 3.35192e-6 4.15149e-6 5.00611e-6 5.87738e-6 6.7182e-6 7.47666e-6 8.1012e-6 8.54628e-6 8.7779e-6 8.7779e-6 8.54628e-6 8.1012e-6 7.47666e-6 6.7182e-6 5.87738e-6 5.00611e-6 4.15149e-6 3.35192e-6 2.63493e-6 2.01665e-6 1.50272e-6 1.09022e-6 7.70077e-7 5.29591e-7 3.54596e-7 2.3116e-7 1.46716e-7 9.06628e-8 5.45465e-8 3.19514e-8 1.82222e-8 1.0118e-8 5.46989e-9
-4.14112e-10 9.59111e-9 1.77413e-8 3.19514e-8 5.60248e-8 9.56437e-8 1.58971e-7 2.57257e-7 4.05324e-7 6.21761e-7 9.28604e-7 1.35028e-6 1.91163e-6 2.63493e-6 3.53607e-6 4.62018e-6 5.87738e-6 7.27937e-6 8.7779e-6 1.03056e-5 1.17799e-5 1.31099e-5 1.42049e-5 1.49854e-5 1.53915e-5 1.53915e-5 1.49854e-5 1.42049e-5 1.31099e-5 1.17799e-5 1.03056e-5 8.7779e-6 7.27937e-6 5.87738e-6 4.62018e-6 3.53607e-6 2.63493e-6 1.91163e-6 1.35028e-6 9.28604e-7 6.21761e-7 4.05324e-7 2.57257e-7 1.58971e-7 9.56437e-8 5.60248e-8 3.19514e-8 1.77413e-8 9.59111e-9
-3.93911e-10 1.63736e-8 3.02874e-8 5.45465e-8 9.56437e-8 1.6328e-7 2.71391e-7 4.39181e-7 6.91956e-7 1.06145e-6 1.58528e-6 2.30516e-6 3.26347e-6 4.49827e-6 6.03667e-6 7.88743e-6 1.00337e-5 1.24271e-5 1.49854e-5 1.75934e-5 2.01103e-5 2.23807e-5 2.42502e-5 2.55825e-5 2.62759e-5 2.62759e-5 2.55825e-5 2.42502e-5 2.23807e-5 2.01103e-5 1.75934e-5 1.49854e-5 1.24271e-5 1.00337e-5 7.88743e-6 6.03667e-6 4.49827e-6 3.26347e-6 2.30516e-6 1.58528e-6 1.06145e-6 6.91956e-7 4.39181e-7 2.71391e-7 1.6328e-7 9.56437e-8 5.45465e-8 3.02874e-8 1.63736e-8
-3.73711e-10 2.72149e-8 5.03414e-8 9.06628e-8 1.58971e-7 2.71391e-7 4.51084e-7 7.29972e-7 1.15011e-6 1.76426e-6 2.63493e-6 3.83144e-6 5.42428e-6 7.47666e-6 1.00337e-5 1.31099e-5 1.66772e-5 2.06554e-5 2.49075e-5 2.92423e-5 3.34258e-5 3.71995e-5 4.03068e-5 4.25212e-5 4.36736e-5 4.36736e-5 4.25212e-5 4.03068e-5 3.71995e-5 3.34258e-5 2.92423e-5 2.49075e-5 2.06554e-5 1.66772e-5 1.31099e-5 1.00337e-5 7.47666e-6 5.42428e-6 3.83144e-6 2.63493e-6 1.76426e-6 1.15011e-6 7.29972e-7 4.51084e-7 2.71391e-7 1.58971e-7 9.06628e-8 5.03414e-8 2.72149e-8
-3.5351e-10 4.40409e-8 8.14655e-8 1.46716e-7 2.57257e-7 4.39181e-7 7.29972e-7 1.18129e-6 1.86119e-6 2.85503e-6 4.26401e-6 6.20028e-6 8.7779e-6 1.20992e-5 1.62371e-5 2.12152e-5 2.6988e-5 3.34258e-5 4.03068e-5 4.73218e-5 5.40916e-5 6.01984e-5 6.52269e-5 6.88104e-5 7.06754e-5 7.06754e-5 6.88104e-5 6.52269e-5 6.01984e-5 5.40916e-5 4.73218e-5 4.03068e-5 3.34258e-5 2.6988e-5 2.12152e-5 1.62371e-5 1.20992e-5 8.7779e-6 6.20028e-6 4.26401e-6 2.85503e-6 1.86119e-6 1.18129e-6 7.29972e-7 4.39181e-7 2.57257e-7 1.46716e-7 8.14655e-8 4.40409e-8
-3.33309e-10 6.9389e-8 1.28354e-7 2.3116e-7 4.05324e-7 6.91956e-7 1.15011e-6 1.86119e-6 2.93241e-6 4.49827e-6 6.7182e-6 9.76891e-6 1.38301e-5 1.9063e-5 2.55825e-5 3.34258e-5 4.25212e-5 5.26643e-5 6.35057e-5 7.45582e-5 8.52246e-5 9.48462e-5 0.000102769 0.000108415 0.000111353 0.000111353 0.000108415 0.000102769 9.48462e-5 8.52246e-5 7.45582e-5 6.35057e-5 5.26643e-5 4.25212e-5 3.34258e-5 2.55825e-5 1.9063e-5 1.38301e-5 9.76891e-6 6.7182e-6 4.49827e-6 2.93241e-6 1.86119e-6 1.15011e-6 6.91956e-7 4.05324e-7 2.3116e-7 1.28354e-7 6.9389e-8
-3.13109e-10 1.06442e-7 1.96893e-7 3.54596e-7 6.21761e-7 1.06145e-6 1.76426e-6 2.85503e-6 4.49827e-6 6.90028e-6 1.03056e-5 1.49854e-5 2.12152e-5 2.92423e-5 3.92432e-5 5.12746e-5 6.52269e-5 8.07862e-5 9.74168e-5 0.000114371 0.000130733 0.000145493 0.000157646 0.000166307 0.000170814 0.000170814 0.000166307 0.000157646 0.000145493 0.000130733 0.000114371 9.74168e-5 8.07862e-5 6.52269e-5 5.12746e-5 3.92432e-5 2.92423e-5 2.12152e-5 1.49854e-5 1.03056e-5 6.90028e-6 4.49827e-6 2.85503e-6 1.76426e-6 1.06145e-6 6.21761e-7 3.54596e-7 1.96893e-7 1.06442e-7
-2.92908e-10 1.58971e-7 2.9406e-7 5.29591e-7 9.28604e-7 1.58528e-6 2.63493e-6 4.26401e-6 6.7182e-6 1.03056e-5 1.53915e-5 2.23807e-5 3.1685e-5 4.36736e-5 5.861e-5 7.6579e-5 9.74168e-5 0.000120655 0.000145493 0.000170814 0.000195251 0.000217294 0.000235445 0.00024838 0.000255112 0.000255112 0.00024838 0.000235445 0.000217294 0.000195251 0.000170814 0.000145493 0.000120655 9.74168e-5 7.6579e-5 5.861e-5 4.36736e-5 3.1685e-5 2.23807e-5 1.53915e-5 1.03056e-5 6.7182e-6 4.26401e-6 2.63493e-6 1.58528e-6 9.28604e-7 5.29591e-7 2.9406e-7 1.58971e-7
-2.72708e-10 2.3116e-7 4.27592e-7 7.70077e-7 1.35028e-6 2.30516e-6 3.83144e-6 6.20028e-6 9.76891e-6 1.49854e-5 2.23807e-5 3.25437e-5 4.60731e-5 6.35057e-5 8.52246e-5 0.000111353 0.000141653 0.000175444 0.00021156 0.00024838 0.000283914 0.000315967 0.00034236 0.000361169 0.000370958 0.000370958 0.000361169 0.00034236 0.000315967 0.000283914 0.00024838 0.00021156 0.000175444 0.000141653 0.000111353 8.52246e-5 6.35057e-5 4.60731e-5 3.25437e-5 2.23807e-5 1.49854e-5 9.76891e-6 6.20028e-6 3.83144e-6 2.30516e-6 1.35028e-6 7.70077e-7 4.27592e-7 2.3116e-7
-2.52507e-10 3.27259e-7 6.05354e-7 1.09022e-6 1.91163e-6 3.26347e-6 5.42428e-6 8.7779e-6 1.38301e-5 2.12152e-5 3.1685e-5 4.60731e-5 6.52269e-5 8.99068e-5 0.000120655 0.000157646 0.000200543 0.00024838 0.000299512 0.000351639 0.000401944 0.000447323 0.000484688 0.000511317 0.000525175 0.000525175 0.000511317 0.000484688 0.000447323 0.000401944 0.000351639 0.000299512 0.00024838 0.000200543 0.000157646 0.000120655 8.99068e-5 6.52269e-5 4.60731e-5 3.1685e-5 2.12152e-5 1.38301e-5 8.7779e-6 5.42428e-6 3.26347e-6 1.91163e-6 1.09022e-6 6.05354e-7 3.27259e-7
-2.32307e-10 4.51084e-7 8.34402e-7 1.50272e-6 2.63493e-6 4.49827e-6 7.47666e-6 1.20992e-5 1.9063e-5 2.92423e-5 4.36736e-5 6.35057e-5 8.99068e-5 0.000123925 0.000166307 0.000217294 0.000276422 0.00034236 0.000412838 0.000484688 0.000554028 0.000616576 0.00066808 0.000704784 0.000723885 0.000723885 0.000704784 0.00066808 0.000616576 0.000554028 0.000484688 0.000412838 0.00034236 0.000276422 0.000217294 0.000166307 0.000123925 8.99068e-5 6.35057e-5 4.36736e-5 2.92423e-5 1.9063e-5 1.20992e-5 7.47666e-6 4.49827e-6 2.63493e-6 1.50272e-6 8.34402e-7 4.51084e-7
-2.12106e-10 6.05354e-7 1.11977e-6 2.01665e-6 3.53607e-6 6.03667e-6 1.00337e-5 1.62371e-5 2.55825e-5 3.92432e-5 5.861e-5 8.52246e-5 0.000120655 0.000166307 0.000223183 0.000291608 0.000370958 0.000459446 0.000554028 0.000650451 0.000743504 0.000827444 0.000896562 0.000945818 0.000971453 0.000971453 0.000945818 0.000896562 0.000827444 0.000743504 0.000650451 0.000554028 0.000459446 0.000370958 0.000291608 0.000223183 0.000166307 0.000120655 8.52246e-5 5.861e-5 3.92432e-5 2.55825e-5 1.62371e-5 1.00337e-5 6.03667e-6 3.53607e-6 2.01665e-6 1.11977e-6 6.05354e-7
-1.91905e-10 7.90948e-7 1.46307e-6 2.63493e-6 4.62018e-6 7.88743e-6 1.31099e-5 2.12152e-5 3.34258e-5 5.12746e-5 7.6579e-5 0.000111353 0.000157646 0.000217294 0.000291608 0.000381012 0.000484688 0.000600306 0.000723885 0.00084987 0.000971453 0.00108113 0.00117144 0.00123579 0.00126929 0.00126929 0.00123579 0.00117144 0.00108113 0.000971453 0.00084987 0.000723885 0.000600306 0.000484688 0.000381012 0.000291608 0.000217294 0.000157646 0.000111353 7.6579e-5 5.12746e-5 3.34258e-5 2.12152e-5 1.31099e-5 7.88743e-6 4.62018e-6 2.63493e-6 1.46307e-6 7.90948e-7
-1.71705e-10 1.00617e-6 1.86119e-6 3.35192e-6 5.87738e-6 1.00337e-5 1.66772e-5 2.6988e-5 4.25212e-5 6.52269e-5 9.74168e-5 0.000141653 0.000200543 0.000276422 0.000370958 0.000484688 0.000616576 0.000763655 0.000920861 0.00108113 0.00123579 0.00137531 0.00149019 0.00157206 0.00161467 0.00161467 0.00157206 0.00149019 0.00137531 0.00123579 0.00108113 0.000920861 0.000763655 0.000616576 0.000484688 0.000370958 0.000276422 0.000200543 0.000141653 9.74168e-5 6.52269e-5 4.25212e-5 2.6988e-5 1.66772e-5 1.00337e-5 5.87738e-6 3.35192e-6 1.86119e-6 1.00617e-6
-1.51504e-10 1.24618e-6 2.30516e-6 4.15149e-6 7.27937e-6 1.24271e-5 2.06554e-5 3.34258e-5 5.26643e-5 8.07862e-5 0.000120655 0.000175444 0.00024838 0.00034236 0.000459446 0.000600306 0.000763655 0.000945818 0.00114052 0.00133902 0.00153058 0.00170338 0.00184567 0.00194707 0.00199984 0.00199984 0.00194707 0.00184567 0.00170338 0.00153058 0.00133902 0.00114052 0.000945818 0.000763655 0.000600306 0.000459446 0.00034236 0.00024838 0.000175444 0.000120655 8.07862e-5 5.26643e-5 3.34258e-5 2.06554e-5 1.24271e-5 7.27937e-6 4.15149e-6 2.30516e-6 1.24618e-6
-1.31304e-10 1.50272e-6 2.77969e-6 5.00611e-6 8.7779e-6 1.49854e-5 2.49075e-5 4.03068e-5 6.35057e-5 9.74168e-5 0.000145493 0.00021156 0.000299512 0.000412838 0.000554028 0.000723885 0.000920861 0.00114052 0.00137531 0.00161467 0.00184567 0.00205404 0.00222561 0.00234789 0.00241152 0.00241152 0.00234789 0.00222561 0.00205404 0.00184567 0.00161467 0.00137531 0.00114052 0.000920861 0.000723885 0.000554028 0.000412838 0.000299512 0.00021156 0.000145493 9.74168e-5 6.35057e-5 4.03068e-5 2.49075e-5 1.49854e-5 8.7779e-6 5.00611e-6 2.77969e-6 1.50272e-6
-1.11103e-10 1.76426e-6 3.26347e-6 5.87738e-6 1.03056e-5 1.75934e-5 2.92423e-5 4.73218e-5 7.45582e-5 0.000114371 0.000170814 0.00024838 0.000351639 0.000484688 0.000650451 0.00084987 0.00108113 0.00133902 0.00161467 0.00189569 0.00216689 0.00241152 0.00261296 0.00275651 0.00283122 0.00283122 0.00275651 0.00261296 0.00241152 0.00216689 0.00189569 0.00161467 0.00133902 0.00108113 0.00084987 0.000650451 0.000484688 0.000351639 0.00024838 0.000170814 0.000114371 7.45582e-5 4.73218e-5 2.92423e-5 1.75934e-5 1.03056e-5 5.87738e-6 3.26347e-6 1.76426e-6
-9.09026e-11 2.01665e-6 3.73034e-6 6.7182e-6 1.17799e-5 2.01103e-5 3.34258e-5 5.40916e-5 8.52246e-5 0.000130733 0.000195251 0.000283914 0.000401944 0.000554028 0.000743504 0.000971453 0.00123579 0.00153058 0.00184567 0.00216689 0.00247688 0.00275651 0.00298677 0.00315086 0.00323626 0.00323626 0.00315086 0.00298677 0.00275651 0.00247688 0.00216689 0.00184567 0.00153058 0.00123579 0.000971453 0.000743504 0.000554028 0.000401944 0.000283914 0.000195251 0.000130733 8.52246e-5 5.40916e-5 3.34258e-5 2.01103e-5 1.17799e-5 6.7182e-6 3.73034e-6 2.01665e-6
-7.0702e-11 2.24433e-6 4.15149e-6 7.47666e-6 1.31099e-5 2.23807e-5 3.71995e-5 6.01984e-5 9.48462e-5 0.000145493 0.000217294 0.000315967 0.000447323 0.000616576 0.000827444 0.00108113 0.00137531 0.00170338 0.00205404 0.00241152 0.00275651 0.00306772 0.00332397 0.00350659 0.00360162 0.00360162 0.00350659 0.00332397 0.00306772 0.00275651 0.00241152 0.00205404 0.00170338 0.00137531 0.00108113 0.000827444 0.000616576 0.000447323 0.000315967 0.000217294 0.000145493 9.48462e-5 6.01984e-5 3.71995e-5 2.23807e-5 1.31099e-5 7.47666e-6 4.15149e-6 2.24433e-6
-5.05014e-11 2.4318e-6 4.49827e-6 8.1012e-6 1.42049e-5 2.42502e-5 4.03068e-5 6.52269e-5 0.000102769 0.000157646 0.000235445 0.00034236 0.000484688 0.00066808 0.000896562 0.00117144 0.00149019 0.00184567 0.00222561 0.00261296 0.00298677 0.00332397 0.00360162 0.0037995 0.00390247 0.00390247 0.0037995 0.00360162 0.00332397 0.00298677 0.00261296 0.00222561 0.00184567 0.00149019 0.00117144 0.000896562 0.00066808 0.000484688 0.00034236 0.000235445 0.000157646 0.000102769 6.52269e-5 4.03068e-5 2.42502e-5 1.42049e-5 8.1012e-6 4.49827e-6 2.4318e-6
-3.03009e-11 2.5654e-6 4.7454e-6 8.54628e-6 1.49854e-5 2.55825e-5 4.25212e-5 6.88104e-5 0.000108415 0.000166307 0.00024838 0.000361169 0.000511317 0.000704784 0.000945818 0.00123579 0.00157206 0.00194707 0.00234789 0.00275651 0.00315086 0.00350659 0.0037995 0.00400824 0.00411687 0.00411687 0.00400824 0.0037995 0.00350659 0.00315086 0.00275651 0.00234789 0.00194707 0.00157206 0.00123579 0.000945818 0.000704784 0.000511317 0.000361169 0.00024838 0.000166307 0.000108415 6.88104e-5 4.25212e-5 2.55825e-5 1.49854e-5 8.54628e-6 4.7454e-6 2.5654e-6
-1.01003e-11 2.63493e-6 4.87402e-6 8.7779e-6 1.53915e-5 2.62759e-5 4.36736e-5 7.06754e-5 0.000111353 0.000170814 0.000255112 0.000370958 0.000525175 0.000723885 0.000971453 0.00126929 0.00161467 0.00199984 0.00241152 0.00283122 0.00323626 0.00360162 0.00390247 0.00411687 0.00422845 0.00422845 0.00411687 0.00390247 0.00360162 0.00323626 0.00283122 0.00241152 0.00199984 0.00161467 0.00126929 0.000971453 0.000723885 0.000525175 0.000370958 0.000255112 0.000170814 0.000111353 7.06754e-5 4.36736e-5 2.62759e-5 1.53915e-5 8.7779e-6 4.87402e-6 2.63493e-6
1.01003e-11 2.63493e-6 4.87402e-6 8.7779e-6 1.53915e-5 2.62759e-5 4.36736e-5 7.06754e-5 0.000111353 0.000170814 0.000255112 0.000370958 0.000525175 0.000723885 0.000971453 0.00126929 0.00161467 0.00199984 0.00241152 0.00283122 0.00323626 0.00360162 0.00390247 0.00411687 0.00422845 0.00422845 0.00411687 0.00390247 0.00360162 0.00323626 0.00283122 0.00241152 0.00199984 0.00161467 0.00126929 0.000971453 0.000723885 0.000525175 0.000370958 0.000255112 0.000170814 0.000111353 7.06754e-5 4.36736e-5 2.62759e-5 1.53915e-5 8.7779e-6 4.87402e-6 2.63493e-6
3.03009e-11 2.5654e-6 4.7454e-6 8.54628e-6 1.49854e-5 2.55825e-5 4.25212e-5 6.88104e-5 0.000108415 0.000166307 0.00024838 0.000361169 0.000511317 0.000704784 0.000945818 0.00123579 0.00157206 0.00194707 0.00234789 0.00275651 0.00315086 0.00350659 0.0037995 0.00400824 0.00411687 0.00411687 0.00400824 0.0037995 0.00350659 0.00315086 0.00275651 0.00234789 0.00194707 0.00157206 0.00123579 0.000945818 0.000704784 0.000511317 0.000361169 0.00024838 0.000166307 0.000108415 6.88104e-5 4.25212e-5 2.55825e-5 1.49854e-5 8.54628e-6 4.7454e-6 2.5654e-6
5.05014e-11 2.4318e-6 4.49827e-6 8.1012e-6 1.42049e-5 2.42502e-5 4.03068e-5 6.52269e-5 0.000102769 0.000157646 0.000235445 0.00034236 0.000484688 0.00066808 0.000896562 0.00117144 0.00149019 0.00184567 0.00222561 0.00261296 0.00298677 0.00332397 0.00360162 0.0037995 0.00390247 0.00390247 0.0037995 0.00360162 0.00332397 0.00298677 0.00261296 0.00222561 0.00184567 0.00149019 0.00117144 0.000896562 0.00066808 0.000484688 0.00034236 0.000235445 0.000157646 0.000102769 6.52269e-5 4.03068e-5 2.42502e-5 1.42049e-5 8.1012e-6 4.49827e-6 2.4318e-6
7.0702e-11 2.24433e-6 4.15149e-6 7.47666e-6 1.31099e-5 2.23807e-5 3.71995e-5 6.01984e-5 9.48462e-5 0.000145493 0.000217294 0.000315967 0.000447323 0.000616576 0.000827444 0.00108113 0.00137531 0.00170338 0.00205404 0.00241152 0.00275651 0.00306772 0.00332397 0.00350659 0.00360162 0.00360162 0.00350659 0.00332397 0.00306772 0.00275651 0.00241152 0.00205404 0.00170338 0.00137531 0.00108113 0.000827444 0.000616576 0.000447323 0.000315967 0.000217294 0.000145493 9.48462e-5 6.01984e-5 3.71995e-5 2.23807e-5 1.31099e-5 7.47666e-6 4.15149e-6 2.24433e-6
9.09026e-11 2.01665e-6 3.73034e-6 6.7182e-6 1.17799e-5 2.01103e-5 3.34258e-5 5.40916e-5 8.52246e-5 0.000130733 0.000195251 0.000283914 0.000401944 0.000554028 0.000743504 0.000971453 0.00123579 0.00153058 0.00184567 0.00216689 0.00247688 0.00275651 0.00298677 0.00315086 0.00323626 0.00323626 0.00315086 0.00298677 0.00275651 0.00247688 0.00216689 0.00184567 0.00153058 0.00123579 0.000971453 0.000743504 0.000554028 0.000401944 0.000283914 0.000195251 0.000130733 8.52246e-5 5.40916e-5 3.34258e-5 2.01103e-5 1.17799e-5 6.7182e-6 3.73034e-6 2.01665e-6
1.11103e-10 1.76426e-6 3.26347e-6 5.87738e-6 1.03056e-5 1.75934e-5 2.92423e-5 4.73218e-5 7.45582e-5 0.000114371 0.000170814 0.00024838 0.000351639 0.000484688 0.000650451 0.00084987 0.00108113 0.00133902 0.00161467 0.00189569 0.00216689 0.00241152 0.00261296 0.00275651 0.00283122 0.00283122 0.00275651 0.00261296 0.00241152 0.00216689 0.00189569 0.00161467 0.00133902 0.00108113 0.00084987 0.000650451 0.000484688 0.000351639 0.00024838 0.000170814 0.000114371 7.45582e-5 4.73218e-5 2.92423e-5 1.75934e-5 1.03056e-5 5.87738e-6 3.26347e-6 1.76426e-6
1.31304e-10 1.50272e-6 2.77969e-6 5.00611e-6 8.7779e-6 1.49854e-5 2.49075e-5 4.03068e-5 6.35057e-5 9.74168e-5 0.000145493 0.00021156 0.000299512 0.000412838 0.000554028 0.000723885 0.000920861 0.00114052 0.00137531 0.00161467 0.00184567 0.00205404 0.00222561 0.00234789 0.00241152 0.00241152 0.00234789 0.00222561 0.00205404 0.00184567 0.00161467 0.00137531 0.00114052 0.000920861 0.000723885 0.000554028 0.000412838 0.000299512 0.00021156 0.000145493 9.74168e-5 6.35057e-5 4.03068e-5 2.49075e-5 1.49854e-5 8.7779e-6 5.00611e-6 2.77969e-6 1.50272e-6
1.51504e-10 1.24618e-6 2.30516e-6 4.15149e-6 7.27937e-6 1.24271e-5 2.06554e-5 3.34258e-5 5.26643e-5 8.07862e-5 0.000120655 0.000175444 0.00024838 0.00034236 0.000459446 0.000600306 0.000763655 0.000945818 0.00114052 0.00133902 0.00153058 0.00170338 0.00184567 0.00194707 0.00199984 0.00199984 0.00194707 0.00184567 0.00170338 0.00153058 0.00133902 0.00114052 0.000945818 0.000763655 0.000600306 0.000459446 0.00034236 0.00024838 0.000175444 0.000120655 8.07862e-5 5.26643e-5 3.34258e-5 2.06554e-5 1.24271e-5 7.27937e-6 4.15149e-6 2.30516e-6 1.24618e-6
1.71705e-10 1.00617e-6 1.86119e-6 3.35192e-6 5.87738e-6 1.00337e-5 1.66772e-5 2.6988e-5 4.25212e-5 6.52269e-5 9.74168e-5 0.000141653 0.000200543 0.000276422 0.000370958 0.000484688 0.000616576 0.000763655 0.000920861 0.00108113 0.00123579 0.00137531 0.00149019 0.00157206 0.00161467 0.00161467 0.00157206 0.00149019 0.00137531 0.00123579 0.00108113 0.000920861 0.000763655 0.000616576 0.000484688 0.000370958 0.000276422 0.000200543 0.000141653 9.74168e-5 6.52269e-5 4.25212e-5 2.6988e-5 1.66772e-5 1.00337e-5 5.87738e-6 3.35192e-6 1.86119e-6 1.00617e-6
1.91905e-10 7.90948e-7 1.46307e-6 2.63493e-6 4.62018e-6 7.88743e-6 1.31099e-5 2.12152e-5 3.34258e-5 5.12746e-5 7.6579e-5 0.000111353 0.000157646 0.000217294 0.000291608 0.000381012 0.000484688 0.000600306 0.000723885 0.00084987 0.000971453 0.00108113 0.00117144 0.00123579 0.00126929 0.00126929 0.00123579 0.00117144 0.00108113 0.000971453 0.00084987 0.000723885 0.000600306 0.000484688 0.000381012 0.000291608 0.000217294 0.000157646 0.000111353 7.6579e-5 5.12746e-5 3.34258e-5 2.12152e-5 1.31099e-5 7.88743e-6 4.62018e-6 2.63493e-6 1.46307e-6 7.90948e-7
2.12106e-10 6.05354e-7 1.11977e-6 2.01665e-6 3.53607e-6 6.03667e-6 1.00337e-5 1.62371e-5 2.55825e-5 3.92432e-5 5.861e-5 8.52246e-5 0.000120655 0.000166307 0.000223183 0.000291608 0.000370958 0.000459446 0.000554028 0.000650451 0.000743504 0.000827444 0.000896562 0.000945818 0.000971453 0.000971453 0.000945818 0.000896562 0.000827444 0.000743504 0.000650451 0.000554028 0.000459446 0.000370958 0.000291608 0.000223183 0.000166307 0.000120655 8.52246e-5 5.861e-5 3.92432e-5 2.55825e-5 1.62371e-5 1.00337e-5 6.03667e-6 3.53607e-6 2.01665e-6 1.11977e-6 6.05354e-7
2.32307e-10 4.51084e-7 8.34402e-7 1.50272e-6 2.63493e-6 4.49827e-6 7.47666e-6 1.20992e-5 1.9063e-5 2.92423e-5 4.36736e-5 6.35057e-5 8.99068e-5 0.000123925 0.000166307 0.000217294 0.000276422 0.00034236 0.000412838 0.000484688 0.000554028 0.000616576 0.00066808 0.000704784 0.000723885 0.000723885 0.000704784 0.00066808 0.000616576 0.000554028 0.000484688 0.000412838 0.00034236 0.000276422 0.000217294 0.000166307 0.000123925 8.99068e-5 6.35057e-5 4.36736e-5 2.92423e-5 1.9063e-5 1.20992e-5 7.47666e-6 4.49827e-6 2.63493e-6 1.50272e-6 8.34402e-7 4.51084e-7
2.52507e-10 3.27259e-7 6.05354e-7 1.09022e-6 1.91163e-6 3.26347e-6 5.42428e-6 8.7779e-6 1.38301e-5 2.12152e-5 3.1685e-5 4.60731e-5 6.52269e-5 8.99068e-5 0.000120655 0.000157646 0.000200543 0.00024838 0.000299512 0.000351639 0.000401944 0.000447323 0.000484688 0.000511317 0.000525175 0.000525175 0.000511317 0.000484688 0.000447323 0.000401944 0.000351639 0.000299512 0.00024838 0.000200543 0.000157646 0.000120655 8.99068e-5 6.52269e-5 4.60731e-5 3.1685e-5 2.12152e-5 1.38301e-5 8.7779e-6 5.42428e-6 3.26347e-6 1.91163e-6 1.09022e-6 6.05354e-7 3.27259e-7
2.72708e-10 2.3116e-7 4.27592e-7 7.70077e-7 1.35028e-6 2.30516e-6 3.83144e-6 6.20028e-6 9.76891e-6 1.49854e-5 2.23807e-5 3.25437e-5 4.60731e-5 6.35057e-5 8.52246e-5 0.000111353 0.000141653 0.000175444 0.00021156 0.00024838 0.000283914 0.000315967 0.00034236 0.000361169 0.000370958 0.000370958 0.000361169 0.00034236 0.000315967 0.000283914 0.00024838 0.00021156 0.000175444 0.000141653 0.000111353 8.52246e-5 6.35057e-5 4.60731e-5 3.25437e-5 2.23807e-5 1.49854e-5 9.76891e-6 6.20028e-6 3.83144e-6 2.30516e-6 1.35028e-6 7.70077e-7 4.27592e-7 2.3116e-7
2.92908e-10 1.58971e-7 2.9406e-7 5.29591e-7 9.28604e-7 1.58528e-6 2.63493e-6 4.26401e-6 6.7182e-6 1.03056e-5 1.53915e-5 2.23807e-5 3.1685e-5 4.36736e-5 5.861e-5 7.6579e-5 9.74168e-5 0.000120655 0.000145493 0.000170814 0.000195251 0.000217294 0.000235445 0.00024838 0.000255112 0.000255112 0.00024838 0.000235445 0.000217294 0.000195251 0.000170814 0.000145493 0.000120655 9.74168e-5 7.6579e-5 5.861e-5 4.36736e-5 3.1685e-5 2.23807e-5 1.53915e-5 1.03056e-5 6.7182e-6 4.26401e-6 2.63493e-6 1.58528e-6 9.28604e-7 5.29591e-7 2.9406e-7 1.58971e-7
3.13109e-10 1.06442e-7 1.96893e-7 3.54596e-7 6.21761e-7 1.06145e-6 1.76426e-6 2.85503e-6 4.49827e-6 6.90028e-6 1.03056e-5 1.49854e-5 2.12152e-5 2.92423e-5 3.92432e-5 5.12746e-5 6.52269e-5 8.07862e-5 9.74168e-5 0.000114371 0.000130733 0.000145493 0.000157646 0.000166307 0.000170814 0.000170814 0.000166307 0.000157646 0.000145493 0.000130733 0.000114371 9.74168e-5 8.07862e-5 6.52269e-5 5.12746e-5 3.92432e-5 2.92423e-5 2.12152e-5 1.49854e-5 1.03056e-5 6.90028e-6 4.49827e-6 2.85503e-6 1.76426e-6 1.06145e-6 6.21761e-7 3.54596e-7 1.96893e-7 1.06442e-7
3.33309e-10 6.9389e-8 1.28354e-7 2.3116e-7 4.05324e-7 6.91956e-7 1.15011e-6 1.86119e-6 2.93241e-6 4.49827e-6 6.7182e-6 9.76891e-6 1.38301e-5 1.9063e-5 2.55825e-5 3.34258e-5 4.25212e-5 5.26643e-5 6.35057e-5 7.45582e-5 8.52246e-5 9.48462e-5 0.000102769 0.000108415 0.000111353 0.000111353 0.000108415 0.000102769 9.48462e-5 8.52246e-5 7.45582e-5 6.35057e-5 5.26643e-5 4.25212e-5 3.34258e-5 2.55825e-5 1.9063e-5 1.38301e-5 9.76891e-6 6.7182e-6 4.49827e-6 2.93241e-6 1.86119e-6 1.15011e-6 6.91956e-7 4.05324e-7 2.3116e-7 1.28354e-7 6.9389e-8
3.5351e-10 4.40409e-8 8.14655e-8 1.46716e-7 2.57257e-7 4.39181e-7 7.29972e-7 1.18129e-6 1.86119e-6 2.85503e-6 4.26401e-6 6.20028e-6 8.7779e-6 1.20992e-5 1.62371e-5 2.12152e-5 2.6988e-5 3.34258e-5 4.03068e-5 4.73218e-5 5.40916e-5 6.01984e-5 6.52269e-5 6.88104e-5 7.06754e-5 7.06754e-5 6.88104e-5 6.52269e-5 6.01984e-5 5.40916e-5 4.73218e-5 4.03068e-5 3.34258e-5 2.6988e-5 2.12152e-5 1.62371e-5 1.20992e-5 8.7779e-6 6.20028e-6 4.26401e-6 2.85503e-6 1.86119e-6 1.18129e-6 7.29972e-7 4.39181e-7 2.57257e-7 1.46716e-7 8.14655e-8 4.40409e-8
3.73711e-10 2.72149e-8 5.03414e-8 9.06628e-8 1.58971e-7 2.71391e-7 4.51084e-7 7.29972e-7 1.15011e-6 1.76426e-6 2.63493e-6 3.83144e-6 5.42428e-6 7.47666e-6 1.00337e-5 1.31099e-5 1.66772e-5 2.06554e-5 2.49075e-5 2.92423e-5 3.34258e-5 3.71995e-5 4.03068e-5 4.25212e-5 4.36736e-5 4.36736e-5 4.25212e-5 4.03068e-5 3.71995e-5 3.34258e-5 2.92423e-5 2.49075e-5 2.06554e-5 1.66772e-5 1.31099e-5 1.00337e-5 7.47666e-6 5.42428e-6 3.83144e-6 2.63493e-6 1.76426e-6 1.15011e-6 7.29972e-7 4.51084e-7 2.71391e-7 1.58971e-7 9.06628e-8 5.03414e-8 2.72149e-8
3.93911e-10 1.63736e-8 3.02874e-8 5.45465e-8 9.56437e-8 1.6328e-7 2.71391e-7 4.39181e-7 6.91956e-7 1.06145e-6 1.58528e-6 2.30516e-6 3.26347e-6 4.49827e-6 6.03667e-6 7.88743e-6 1.00337e-5 1.24271e-5 1.49854e-5 1.75934e-5 2.01103e-5 2.23807e-5 2.42502e-5 2.55825e-5 2.62759e-5 2.62759e-5 2.55825e-5 2.42502e-5 2.23807e-5 2.01103e-5 1.75934e-5 1.49854e-5 1.24271e-5 1.00337e-5 7.88743e-6 6.03667e-6 4.49827e-6 3.26347e-6 2.30516e-6 1.58528e-6 1.06145e-6 6.91956e-7 4.39181e-7 2.71391e-7 1.6328e-7 9.56437e-8 5.45465e-8 3.02874e-8 1.63736e-8
4.14112e-10 9.59111e-9 1.77413e-8 3.19514e-8 5.60248e-8 9.56437e-8 1.58971e-7 2.57257e-7 4.05324e-7 6.21761e-7 9.28604e-7 1.35028e-6 1.91163e-6 2.63493e-6 3.53607e-6 4.62018e-6 5.87738e-6 7.27937e-6 8.7779e-6 1.03056e-5 1.17799e-5 1.31099e-5 1.42049e-5 1.49854e-5 1.53915e-5 1.53915e-5 1.49854e-5 1.42049e-5 1.31099e-5 1.17799e-5 1.03056e-5 8.7779e-6 7.27937e-6 5.87738e-6 4.62018e-6 3.53607e-6 2.63493e-6 1.91163e-6 1.35028e-6 9.28604e-7 6.21761e-7 4.05324e-7 2.57257e-7 1.58971e-7 9.56437e-8 5.60248e-8 3.19514e-8 1.77413e-8 9.59111e-9
4.34312e-10 5.46989e-9 1.0118e-8 1.82222e-8 3.19514e-8 5.45465e-8 9.06628e-8 1.46716e-7 2.3116e-7 3.54596e-7 5.29591e-7 7.70077e-7 1.09022e-6 1.50272e-6 2.01665e-6 2.63493e-6 3.35192e-6 4.15149e-6 5.00611e-6 5.87738e-6 6.7182e-6 7.47666e-6 8.1012e-6 8.54628e-6 8.7779e-6 8.7779e-6 8.54628e-6 8.1012e-6 7.47666e-6 6.7182e-6 5.87738e-6 5.00611e-6 4.15149e-6 3.35192e-6 2.63493e-6 2.01665e-6 1.50272e-6 1.09022e-6 7.70077e-7 5.29591e-7 3.54596e-7 2.3116e-7 1.46716e-7 9.06628e-8 5.45465e-8 3.19514e-8 1.82222e-8 1.0118e-8 5.46989e-9
4.54513e-10 3.03721e-9 5.61814e-9 1.0118e-8 1.77413e-8 3.02874e-8 5.03414e-8 8.14655e-8 1.28354e-7 1.96893e-7 2.9406e-7 4.27592e-7 6.05354e-7 8.34402e-7 1.11977e-6 1.46307e-6 1.86119e-6 2.30516e-6 2.77969e-6 3.26347e-6 3.73034e-6 4.15149e-6 4.49827e-6 4.7454e-6 4.87402e-6 4.87402e-6 4.7454e-6 4.49827e-6 4.15149e-6 3.73034e-6 3.26347e-6 2.77969e-6 2.30516e-6 1.86119e-6 1.46307e-6 1.11977e-6 8.34402e-7 6.05354e-7 4.27592e-7 2.9406e-7 1.96893e-7 1.28354e-7 8.14655e-8 5.03414e-8 3.02874e-8 1.77413e-8 1.0118e-8 5.61814e-9 3.03721e-9
4.74713e-10 1.64194e-9 3.03721e-9 5.46989e-9 9.59111e-9 1.63736e-8 2.72149e-8 4.40409e-8 6.9389e-8 1.06442e-7 1.58971e-7 2.3116e-7 3.27259e-7 4.51084e-7 6.05354e-7 7.90948e-7 1.00617e-6 1.24618e-6 1.50272e-6 1.76426e-6 2.01665e-6 2.24433e-6 2.4318e-6 2.5654e-6 2.63493e-6 2.63493e-6 2.5654e-6 2.4318e-6 2.24433e-6 2.01665e-6 1.76426e-6 1.50272e-6 1.24618e-6 1.00617e-6 7.90948e-7 6.05354e-7 4.51084e-7 3.27259e-7 2.3116e-7 1.58971e-7 1.06442e-7 6.9389e-8 4.40409e-8 2.72149e-8 1.63736e-8 9.59111e-9 5.46989e-9 3.03721e-9 1.64194e-9To make the Gaussian Markov random field efficient we first precompute a bunch of quantities that allow us to scale things linearly with the number of image pixels. The returns a functional that accepts a single argument related to the correlation length of the field. The second argument defines the underlying random field of the Markov process. Here we are using a zero mean and unit variance Gaussian Markov random field. For this tutorial we will use the first order random field
julia
cprior = corr_image_prior(grid, dvis)HierarchicalPrior(
map:
ConditionalMarkov(
Random Field: VLBIImagePriors.GaussMarkovRandomField
Graph: MarkovRandomFieldGraph{1}(
dims: (48, 48)
)
) hyper prior:
Truncated(Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.03623874772901475)
θ: 27.594772520225487
)
; lower=1.0, upper=96.0)
)Now we can construct our sky model. We fix the total flux to the observed value 1.1 because total flux is degenerate with a global shift in the gain amplitudes.
julia
skym = sky(grid; ftot = 1.1, mimg, cprior)SkyModel
with map: _sky_sky
on grid:
RectiGrid(
executor: ComradeBase.Serial()
Dimensions:
(↓ X Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points)
)
)Unlike other imaging examples (e.g., Imaging a Black Hole using only Closure Quantities) we also need to include a model for the instrument, i.e., gains. The gains will be broken into two components
Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
Gain phases which are more difficult to constrain and can shift rapidly.
julia
gain(x) = exp(x.lg + 1im * x.gp)
G = SingleStokesGain(gain)
intpr = (
lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.2)); LM = IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
gp = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
)
intmodel = InstrumentModel(G, intpr)InstrumentModel
with Jones: SingleStokesGain
with reference basis: PolarizedTypes.CirBasis()To form the posterior we just combine the skymodel, instrument model and the data. To utilize gradients of the posterior we also need to load Enzyme.jl. Under the hood, Comrade will use Enzyme to compute the gradients of the posterior.
julia
using Enzyme
post = VLBIPosterior(skym, intmodel, dvis)VLBIPosterior
ObservedSkyModel
with map: _sky_sky
on grid:
FourierDualDomain(
Algorithm: VLBISkyModels.NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 1.0e-9, AbstractNFFTs.TENSOR, 0x00000000)
Image Domain: RectiGrid(
executor: ComradeBase.Serial()
Dimensions:
(↓ X Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.747133960864206e-10, 4.747133960864206e-10, 48) ForwardOrdered Regular Points)
)
Visibility Domain: UnstructuredDomain(
executor: ComradeBase.Serial()
Dimensions:
274-element StructArray(::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}) with eltype @NamedTuple{U::Float64, V::Float64, Ti::Float64, Fr::Float64}:
(U = -4.405690154666661e9, V = -4.523017159111106e9, Ti = 0.9166666567325592, Fr = 2.27070703125e11)
(U = 787577.6145833326, V = -1.6838098888888871e6, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
(U = -4.444299918222218e9, V = -4.597825294222218e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
(U = 1.337045162666665e9, V = -3.765300401777774e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
(U = -1.336260540444443e9, V = 3.763616127999996e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
(U = 4.445088654222218e9, V = 4.596145080888884e9, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
(U = 5.781345607111105e9, V = 8.325259893333325e8, Ti = 1.2166666388511658, Fr = 2.27070703125e11)
(U = 757554.6649305547, V = -1.6707483020833314e6, Ti = 1.516666665673256, Fr = 2.27070703125e11)
(U = 1.4806382151111097e9, V = -3.741479615999996e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
(U = -4.455366328888884e9, V = -4.673060451555551e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
(U = -1.4798758791111097e9, V = 3.739809735111107e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
(U = 4.456123861333328e9, V = 4.671391715555551e9, Ti = 1.516666665673256, Fr = 2.27070703125e11)
(U = 5.936013027555549e9, V = 9.315912497777768e8, Ti = 1.516666665673256, Fr = 2.27070703125e11)
(U = 722830.7065972214, V = -1.6582321493055536e6, Ti = 1.816666603088379, Fr = 2.27070703125e11)
(U = -4.438811278222218e9, V = -4.748261176888884e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
(U = 1.615060401777776e9, V = -3.715306943999996e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
(U = -1.6143345528888872e9, V = 3.713649343999996e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
(U = 4.439536184888884e9, V = 4.746598001777773e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
(U = 6.053865329777771e9, V = 1.0329465084444433e9, Ti = 1.816666603088379, Fr = 2.27070703125e11)
(U = 683620.5937499993, V = -1.6463396631944429e6, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
(U = -4.394744988444439e9, V = -4.822934485333328e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
(U = 1.7394626631111093e9, V = -3.686947740444441e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
(U = -1.7387777884444425e9, V = 3.6853017173333297e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
(U = 4.395427669333329e9, V = 4.821289443555551e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
(U = 6.134203007999993e9, V = 1.1359907128888876e9, Ti = 2.1166666746139526, Fr = 2.27070703125e11)
(U = -1.8643712995555537e9, V = 3.651452543999996e9, Ti = 2.449999988079071, Fr = 2.27070703125e11)
(U = 6.178850375111104e9, V = 1.2516714417777765e9, Ti = 2.449999988079071, Fr = 2.27070703125e11)
(U = 4.314477112888885e9, V = 4.903116785777773e9, Ti = 2.449999988079071, Fr = 2.27070703125e11)
(U = 587276.6180555549, V = -1.6236182569444429e6, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
(U = 1.9658264106666646e9, V = -3.6206990791111073e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
(U = -4.212778346666662e9, V = -4.976839239111106e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
(U = -1.965238759111109e9, V = 3.6190756053333297e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
(U = 6.178607487999993e9, V = 1.356139594666665e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
(U = 4.213364039111107e9, V = 4.975216568888884e9, Ti = 2.7500000596046448, Fr = 2.27070703125e11)
(U = 535806.2881944438, V = -1.6141240694444429e6, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
(U = 2.0544588266666646e9, V = -3.586710200888885e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
(U = -4.085601863111107e9, V = -5.046995640888884e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
(U = -2.0539194026666646e9, V = 3.5850976213333297e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
(U = 4.0861423288888845e9, V = 5.045379299555551e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
(U = 6.140063729777771e9, V = 1.4602817422222207e9, Ti = 3.0500001311302185, Fr = 2.27070703125e11)
(U = 481010.2230902773, V = -1.6055275520833316e6, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
(U = 2.130349845333331e9, V = -3.5513328995555515e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
(U = -3.933104519111107e9, V = -5.114784767999995e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
(U = -2.1298682808888867e9, V = 3.5497276302222185e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
(U = 6.063457521777771e9, V = 1.5634511359999983e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
(U = 3.9335861404444404e9, V = 5.113178993777773e9, Ti = 3.3499998450279236, Fr = 2.27070703125e11)
(U = 416648.93055555515, V = -1.5970948472222206e6, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = 2.578938929777775e9, V = -4.734788238222218e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = -3.735105663999996e9, V = -5.186825713777772e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = 2.1991734328888865e9, V = -3.5106569031111073e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = -2.5785191964444413e9, V = 4.733192191999995e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = -2.198756124444442e9, V = 3.509060273777774e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = 3.797641537777774e8, V = -1.2241308906666653e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = 3.735526620444441e9, V = 5.185227192888883e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = 6.314048611555549e9, V = 4.520333351111106e8, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = 5.934280931555549e9, V = 1.6761668835555537e9, Ti = 3.6833333373069763, Fr = 2.27070703125e11)
(U = 364284.04464285675, V = -1.591376535714284e6, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = -3.5324296035555515e9, V = -5.248264732444439e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = 2.246733795555553e9, V = -3.4730721351111073e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = 2.7040352142222195e9, V = -4.690125752888884e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = -2.687901513142854e9, V = 4.694663167999995e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = -2.2408562468571405e9, V = 3.476569380571425e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = 4.573011991111106e8, V = -1.2170540408888876e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = 3.561408036571425e9, V = 5.238641298285708e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = 5.779150734222217e9, V = 1.7751998862222204e9, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = 6.236464440888882e9, V = 5.581377528888882e8, Ti = 3.98333340883255, Fr = 2.27070703125e11)
(U = 292810.0842803027, V = -1.5850529772727257e6, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = 2.81236712533333e9, V = -4.643490716444439e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = 2.280368668444442e9, V = -3.43479938133333e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = -3.307854179555552e9, V = -5.306092316444439e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = -2.812399111529409e9, V = 4.641748931764702e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = -2.280159088941174e9, V = 3.43308967905882e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = 5.3200019644444394e8, V = -1.2086904924444432e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = 3.307383679999996e9, V = 5.304687194352936e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = 6.120220771555549e9, V = 6.626024835555549e8, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = 5.588211185777772e9, V = 1.8712975182222202e9, Ti = 4.283333480358124, Fr = 2.27070703125e11)
(U = 238138.85156249977, V = -1.5812476874999984e6, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 5.092881095111105e9, V = -4.199598421333329e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 2.903268209777775e9, V = -4.595170062222218e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 2.299867178666664e9, V = -3.396077226666663e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = -3.062767004444441e9, V = -5.359950862222217e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 5.050181503999994e9, V = -4.210810559999995e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = -2.8907034879999967e9, V = 4.600886271999995e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 2.189606492444442e9, V = 3.955690115555551e8, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = -2.2976930559999976e9, V = 3.400280575999996e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 2.79301371733333e9, V = -8.035218719999992e8, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 6.034046399999994e8, V = -1.199091349333332e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 3.1008046079999967e9, V = 5.350639359999994e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 8.155639423999991e9, V = 1.1603480675555544e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 5.966032455111105e9, V = 7.64783736888888e8, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 5.36262916266666e9, V = 1.963875015111109e9, Ti = 4.583333194255829, Fr = 2.27070703125e11)
(U = 163232.87499999983, V = -1.5773965972222206e6, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 5.388079544888883e9, V = -4.101135303111107e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 2.3048084479999976e9, V = -3.3528182257777743e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = -2.768253631999997e9, V = -5.414730879999994e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 5.388055679999994e9, V = -4.101249678222218e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 2.9831417315555525e9, V = -4.539868501333328e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 5.356877311999994e9, V = -4.110922695111107e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = -2.975202047999997e9, V = 4.544589937777773e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 2.4049356302222195e9, V = 4.387324826666662e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 2.4049011413333306e9, V = 4.3861474933333284e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = -2.305032220444442e9, V = 3.356104334222219e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 3.083277411555552e9, V = -7.48315953777777e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 3.0832264319999967e9, V = -7.484370453333325e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 6.783317208888881e8, V = -1.187050453333332e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 2.8027259164444413e9, V = 5.407289457777772e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 8.156328149333324e9, V = 1.3135985066666653e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 8.156297400888881e9, V = 1.3134812622222207e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 5.751382712888883e9, V = 8.74867279999999e8, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 5.073054606222218e9, V = 2.0619144959999979e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = -5.356795278222217e9, V = 4.1110577777777734e9, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 30750.899088541635, V = 116905.95138888876, Ti = 4.916666686534882, Fr = 2.27070703125e11)
(U = 99881.10751488084, V = -1.575297547619046e6, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 5.59465437866666e9, V = -4.018611114666662e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 3.030624810666663e9, V = -4.494683207111106e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 2.2960602239999976e9, V = -3.3182473599999967e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 5.594626830222217e9, V = -4.0187286115555515e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = -2.5172720142222195e9, V = -5.45444669155555e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 5.579493814857137e9, V = -4.023605065142853e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = -3.0274068479999967e9, V = 4.496662674285709e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 2.5640353279999976e9, V = 4.7607165511111057e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 2.564000689777775e9, V = 4.7595243111111057e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = -2.2970926933333306e9, V = 3.3193651078095202e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 3.2985927679999967e9, V = -7.003643235555549e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 7.345665155555549e8, V = -1.1764336924444432e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 3.2985603128888855e9, V = -7.004829137777771e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 2.5374800944761877e9, V = 5.44992148723809e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 8.111924807111103e9, V = 1.4358354488888874e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 4.813340984888884e9, V = 2.136196949333331e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 5.547893447111105e9, V = 9.597677279999989e8, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 8.111898709333324e9, V = 1.4357162488888874e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = -5.579455049142851e9, V = 4.0237272990476146e9, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 27098.63357204858, V = 117340.62348090265, Ti = 5.183333337306976, Fr = 2.27070703125e11)
(U = 45690.53906249995, V = -1.5743099999999984e6, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 5.773825464888883e9, V = -3.933187214222218e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 2.2760630257777753e9, V = -3.283890659555552e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 3.0632664319999967e9, V = -4.448890581333328e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 5.773791857777772e9, V = -3.9333102648888845e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = -2.253957283555553e9, V = -5.490298367999994e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 5.744240127999994e9, V = -3.946954495999996e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = -3.0586252799999967e9, V = 4.455449599999995e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 2.710558015999997e9, V = 5.1570311111111057e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 2.7105363982222195e9, V = 5.1558591999999946e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = -2.2804378879999976e9, V = 3.2883723519999967e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 3.4977656675555515e9, V = -6.492943324444438e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 7.872002471111102e8, V = -1.164999438222221e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 3.497740145777774e9, V = -6.494127164444437e8, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 2.3015257599999976e9, V = 5.482690303999994e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 8.027784391111103e9, V = 1.5571086684444427e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 5.317234005333327e9, V = 1.0414041991111101e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 8.027758634666658e9, V = 1.556993443555554e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 4.530021788444439e9, V = 2.206406378666664e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = -5.744238079999993e9, V = 3.9470612479999957e9, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 23313.64084201386, V = 117719.41384548598, Ti = 5.449999988079071, Fr = 2.27070703125e11)
(U = 5.924709319111105e9, V = -3.8452848639999957e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 5.924688497777772e9, V = -3.845403854222218e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 3.08089890133333e9, V = -4.402725802666661e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 2.244918698666664e9, V = -3.2499217208888855e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 2.8438056746666636e9, V = 5.57437258666666e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 2.8437909475555525e9, V = 5.573208195555549e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 3.6797841706666627e9, V = -5.953657653333327e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 3.679770524444441e9, V = -5.954820355555549e8, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 8.35980113777777e8, V = -1.1528026026666653e9, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = 19414.168077256923, V = 118040.4856770832, Ti = 5.716666638851166, Fr = 2.27070703125e11)
(U = -92093.64322916657, V = -1.5751088611111094e6, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = 3.081715143111108e9, V = -4.344833365333329e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = 6.07242564266666e9, V = -3.732578396444441e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = 2.190544647111109e9, V = -3.208251676444441e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = 6.05836305066666e9, V = -3.7438256071111073e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = -3.0829770239999967e9, V = 4.34976796444444e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = 2.99071138133333e9, V = 6.122509119999994e8, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = -2.1976609848888865e9, V = 3.211308942222219e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = 3.881882190222218e9, V = -5.2432489066666615e8, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = 8.911678026666657e8, V = -1.1365769386666653e9, Ti = 6.049999952316284, Fr = 2.27070703125e11)
(U = -150875.45052083317, V = -1.5769051249999984e6, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 6.157184056888882e9, V = -3.6406887111111073e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 6.157173617777771e9, V = -3.6408075306666627e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 2.1349524657777758e9, V = -3.175750087111108e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 3.065380067555552e9, V = -4.298637105777774e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 6.14621098666666e9, V = -3.6535959466666627e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = -3.0691708159999967e9, V = 4.304277162666662e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 3.091809799111108e9, V = 6.579570684444437e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 3.091795150222219e9, V = 6.578363697777771e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = -2.1445705386666644e9, V = 3.1792032426666636e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 4.0222341617777734e9, V = -4.6493393599999946e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 4.0222294044444404e9, V = -4.650499119999995e8, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 9.304277457777768e8, V = -1.1228893119999988e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = -6.146197247999993e9, V = 3.653717759999996e9, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 10322.971652560753, V = 118544.06944444432, Ti = 6.316666603088379, Fr = 2.27070703125e11)
(U = 6.211784632888882e9, V = -3.5477483875555515e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 6.211780195555549e9, V = -3.547863480888885e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 2.068902503111109e9, V = -3.1441626808888855e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 3.034031331555552e9, V = -4.2528075377777734e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 3.177752227555552e9, V = 7.050634364444437e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 3.177752810666663e9, V = 7.049488924444437e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 4.1428866915555515e9, V = -4.035797448888885e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 9.651305439999989e8, V = -1.108645550222221e9, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 4.1428811519999957e9, V = -4.036979662222218e8, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = 6186.6973334418335, V = 118668.12304687487, Ti = 6.583333253860474, Fr = 2.27070703125e11)
(U = -266936.3645833331, V = -1.583141166666665e6, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 6.235954645333326e9, V = -3.4543274026666627e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 2.987818652444441e9, V = -4.207558826666662e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 6.235956423111104e9, V = -3.454211199999996e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 1.9927185599999979e9, V = -3.1136442239999967e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 6.23440708266666e9, V = -3.468812543999996e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = -2.9972223146666636e9, V = 4.213749674666662e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 3.2481352177777743e9, V = 7.533482702222215e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 3.2481338097777743e9, V = 7.532299875555547e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = -2.006915669333331e9, V = 3.1172445866666636e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 4.243231630222218e9, V = -3.4069073333333296e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 4.243239466666662e9, V = -3.4056780622222185e8, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 9.951025973333323e8, V = -1.0939171839999988e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = -6.23440452266666e9, V = 3.4689294506666627e9, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = 2020.1427137586786, V = 118729.78884548598, Ti = 6.849999904632568, Fr = 2.27070703125e11)
(U = -343950.03624999966, V = -1.589413219999998e6, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 6.223212003555549e9, V = -3.337163527111108e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 6.223215416888882e9, V = -3.337282766222219e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 2.909506687999997e9, V = -4.152156949333329e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 1.8838133048888867e9, V = -3.07722722133333e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 6.225276682239994e9, V = -3.340931624959996e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = -2.913978071039997e9, V = 4.1530716876799955e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 3.313704135111108e9, V = 8.14994561777777e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 3.3137085368888855e9, V = 8.148773351111102e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = -1.889491799039998e9, V = 3.077254901759997e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 4.339403889777773e9, V = -2.5993809599999973e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 4.339408696888885e9, V = -2.6005522088888863e8, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = 1.0256997831111101e9, V = -1.0749325617777767e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = -6.225278402559994e9, V = 3.3410462617599964e9, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = -3199.665127224389, V = 118718.70724826376, Ti = 7.183333396911621, Fr = 2.27070703125e11)
(U = -401459.5972222218, V = -1.5953378194444429e6, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 6.178709461333326e9, V = -3.243972359111108e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 2.8307861475555525e9, V = -4.1090220657777734e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 6.178720483555549e9, V = -3.244096682666663e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 1.786251743999998e9, V = -3.049647466666663e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 6.179113016888882e9, V = -3.242380138666663e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 3.347918890666663e9, V = 8.650494186666657e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = -2.8311915662222195e9, V = 4.107428693333329e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 3.347926620444441e9, V = 8.649315555555546e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = -1.7866596017777758e9, V = 3.048053795555552e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 4.392453788444439e9, V = -1.9432878977777755e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 4.392462890666661e9, V = -1.944448759999998e8, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = 1.0445360586666656e9, V = -1.0593764586666656e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = -6.179119544888882e9, V = 3.2424975217777743e9, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = -7359.556450737839, V = 118639.36892361098, Ti = 7.450000047683716, Fr = 2.27070703125e11)
(U = -453840.1640624995, V = -1.6017645034722206e6, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 2.7382116835555525e9, V = -4.0671804444444404e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 6.103944163555549e9, V = -3.151683946666663e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 1.6799544319999983e9, V = -3.023603527111108e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 6.103956750222216e9, V = -3.1518036479999967e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 6.10439825066666e9, V = -3.150082410666663e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 3.3657344497777743e9, V = 9.154982524444435e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = -2.7386646755555525e9, V = 4.065578346666662e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 3.365744974222219e9, V = 9.153749564444435e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = -1.6804040142222204e9, V = 3.02200078933333e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 4.423988807111106e9, V = -1.2808080399999987e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 4.424001663999995e9, V = -1.2819521311111099e8, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = 1.0582553617777767e9, V = -1.0435755999999989e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = -6.104410581333326e9, V = 3.150201891555552e9, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = -11483.574490017349, V = 118497.78168402765, Ti = 7.7166666984558105, Fr = 2.27070703125e11)
(U = -504000.3645833328, V = -1.6089620833333316e6, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 5.999277468444439e9, V = -3.060741980444441e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 2.6322108302222195e9, V = -4.0268253795555515e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 5.999288348444438e9, V = -3.060856782222219e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 1.565423818666665e9, V = -2.999217870222219e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 5.999780252444438e9, V = -3.0591319537777743e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = -2.6327225742222195e9, V = 4.025219150222218e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 3.3670615039999967e9, V = 9.660856764444433e8, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 3.367077162666663e9, V = 9.659668408888879e8, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 4.433867007999995e9, V = -6.164012344444438e7, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = -1.5659301759999983e9, V = 2.99760936533333e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 4.433851406222218e9, V = -6.1523151111111045e7, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = 1.0667900266666656e9, V = -1.0276084337777767e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = -5.999801315555549e9, V = 3.059254300444441e9, Ti = 7.983333349227905, Fr = 2.27070703125e11)
(U = -15551.297851562484, V = 118294.64453124987, Ti = 7.983333349227905, Fr = 2.27070703125e11)
)
)
ObservedInstrumentModel
with Jones: SingleStokesGain
with reference basis: PolarizedTypes.CirBasis()Data Products: Comrade.EHTVisibilityDatumOptimization and Sampling
Now we need to actually compute our image. For this we will first follow standard approaches in VLBI and find the maximum a posteriori (MAP) estimate of the image and instrument model. This is done using the comrade_opt function which accepts the posterior, an optimization algorithm, and some keyword arguments. For this tutorial we will use the L-BFGS algorithm. For more information about the optimization algorithms available see the Optimization.jl docs.
julia
using Optimization, OptimizationLBFGSB
xopt, sol = comrade_opt(
post, LBFGSB(); initial_params = prior_sample(rng, post),
maxiters = 2000, g_tol = 1.0e-1
);Warning
Fitting gains tends to be very difficult, meaning that optimization can take a lot longer. The upside is that we usually get nicer images.
First we will evaluate our fit by plotting the residuals
julia
using CairoMakie
using DisplayAs
res = residuals(post, xopt)
plotfields(res[1], :uvdist, :res) |> DisplayAs.PNG |> DisplayAs.Text
These look reasonable, although there may be some minor overfitting. This could be improved in a few ways, but that is beyond the goal of this quick tutorial. Plotting the image, we see that we have a much cleaner version of the closure-only image from Imaging a Black Hole using only Closure Quantities.
julia
g = imagepixels(fovx, fovy, 128, 128)
img = intensitymap(skymodel(post, xopt), g)
imageviz(img, size = (500, 400)) |> DisplayAs.PNG |> DisplayAs.Text
Because we also fit the instrument model, we can inspect their parameters. First, let's query the posterior object with the optimal parameters to get the instrument model.
julia
intopt = instrumentmodel(post, xopt)126-element Comrade.SiteArray{ComplexF64, 1, Vector{ComplexF64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}:
0.9998314993460414 + 0.0im
-0.6161752412174768 + 0.7871195849877409im
1.0101103009739858 + 0.0im
-0.5974342447980401 - 0.8350839774074871im
-0.4054936644642261 + 0.7245001322504171im
-0.5197965486303293 + 0.9503602534651909im
1.0051934067466768 + 0.0im
-0.6398520928971607 - 0.8069274411567765im
-0.5458329976777295 + 0.8180925699561006im
-0.30200027969810345 + 0.9959470952170367im
1.019427727067682 + 0.0im
-0.660081619241941 - 0.7766677449545305im
-0.5449090864349937 + 0.7395375167806303im
-0.11437022489889351 + 1.0088028477798225im
0.9964440832070759 + 0.0im
-0.7127821336830773 - 0.7537340905292965im
-0.60095385150104 + 0.7171530600780327im
0.07382337837844767 + 0.9581373870298888im
1.0359324333401936 + 0.0im
-0.6338881028419866 + 0.24667133394810792im
-0.7570051843437985 - 0.464808147047603im
1.0207207086852752 + 0.0im
-0.7535455357325787 - 0.6869864799205821im
-0.40493393529622795 + 0.479016000819441im
0.43161983025261 + 0.8872926365667612im
1.0090394957058721 + 0.0im
-0.7752409572963028 - 0.6621729245226519im
-0.37634878523589665 + 0.4493905823953936im
0.6102865368899363 + 0.8083483617128656im
1.0209362111260323 + 0.0im
-0.7944181312460461 - 0.6376044406727837im
-0.36334055811338867 + 0.4661593393478189im
0.7642609448278811 + 0.6633916589477035im
1.0041187561573797 + 0.0im
-0.8040610723348002 - 0.6404393193511961im
-0.2356441110038073 + 0.9763812673846272im
-0.30220636573457077 + 0.44936828978543625im
0.9282146450693155 + 0.471480310660943im
0.9875335596393341 + 0.0im
-0.7995211138738271 - 0.6782648192532074im
-0.07565206990931189 + 0.9553034929571882im
-0.28592140912536623 + 0.5335857988848399im
0.9780340336940209 + 0.22241123384400385im
1.0505799411129848 + 0.0im
-0.771251666470603 - 0.6301595660630057im
0.11411846953473051 + 0.9716616594332879im
-0.24813488986305218 + 0.7001101325062572im
1.0617107515958 - 0.024839156360999558im
0.9849180549364352 + 0.0im
-0.7553076135672728 - 0.7313288132089844im
0.31286402593375495 + 0.9302953459439481im
0.7781799678471721 - 0.6680119176823627im
-0.15279272047853576 + 0.6589527863405316im
1.0496488852023131 - 0.327139570360232im
1.036820623387895 + 0.0im
-0.7523359820810424 - 0.655486485194263im
0.5279044881208667 + 0.8370111674347079im
0.8738583602547666 - 0.5241771422827712im
-0.031300634373475016 + 0.670591089765447im
0.9235329122469286 - 0.6295609875431064im
-0.21813973112911655 - 0.9873882129225111im
0.9663878267147861 + 0.0im
-0.7715959484247773 - 0.7282342048043099im
0.6850892781485982 + 0.7193286771143406im
0.910300031748818 - 0.4260129861046372im
0.09935399479971224 + 0.6266489793179821im
0.882918466641291 - 0.9251119520718835im
-0.4443735647643654 - 0.911715543885358im
0.8869710043080896 + 0.0im
-0.7784268861490463 - 0.8649234521873996im
0.790391509240728 + 0.586274777390951im
0.9625134838098721 - 0.3381835943530614im
0.24093244898466096 + 0.598356216707624im
0.5240837179610733 - 0.8464238035515451im
-0.6825675556938772 - 0.7697857918316992im
1.030065041242358 + 0.0im
0.9174190783231865 + 0.3878726909523574im
0.9770434191429721 - 0.27729169316004876im
0.4336601251714223 + 0.4767895888978542im
-0.8549574966739072 - 0.5556997714405668im
1.0708596009843294 + 0.0im
-0.5931559467404359 - 0.7735887378744822im
1.0071593527434357 + 0.10413003395308772im
1.012869063087588 - 0.2314837919641679im
0.6790724689194921 + 0.05501403204369156im
1.0214520461783725 + 0.0im
-0.5666873919013861 - 0.8424827658026268im
1.0015969295502896 - 0.09672670304067059im
0.9906654239025717 - 0.22264509415575526im
0.4815796697887182 - 0.31852369074869125im
-0.9895873568492877 - 0.20553875645690214im
1.0866699024533295 + 0.0im
0.9680417973677655 - 0.3615391253463441im
0.9938354940091062 - 0.21598996861857364im
0.33127899430826774 - 0.5187741030253592im
-1.015262506675966 - 0.03327955473924759im
1.1090431729832717 + 0.0im
-0.4790977928407483 - 0.8032167128387977im
0.8177856041683886 - 0.6033134228483182im
0.9893743778206388 - 0.2281329745195393im
0.1355211760331925 - 0.5035675535701557im
-1.0026945430679037 + 0.14108623282797936im
1.0292048121437136 + 0.0im
-0.4277412489460898 - 0.9122888131468191im
0.6693726239734855 - 0.7769127561420285im
0.9825011731506846 - 0.25936292387498794im
-0.014639718735743679 - 0.5654935300351742im
-0.9929498650570969 + 0.23521318892335016im
1.0044120123110625 + 0.0im
-0.40131832282747265 - 0.9547951088390516im
0.4744612074698839 - 0.894981542175112im
0.9731262807402614 - 0.29033156225238005im
-0.14656626505859038 - 0.7005519749689655im
-0.9727766963553338 + 0.3096098821916213im
0.9809512227772282 + 0.0im
-0.32905054344643586 - 0.9924096250866589im
0.2784736328203881 - 0.9514222104556507im
0.9580392077187571 - 0.333718483120392im
-0.25602474497601013 - 0.651964493719375im
-0.9625288963326686 + 0.314014684968358im
1.010340117249841 + 0.0im
-0.2639439737764915 - 0.9940988139488873im
0.17683230008901474 - 0.9012971926186393im
0.9620964298074085 - 0.34508312478043823im
-0.32612268864747335 - 0.5856660889390896im
-0.9577102553010632 + 0.32540891399578137imThis returns a SiteArray object which contains the gains as a flat vector with metadata about the sites, time, and frequency. To visualize the gains we can use the plotcaltable function which automatically plots the gains. Since the gains are complex we will first plot the phases.
julia
plotcaltable(angle.(intopt)) |> DisplayAs.PNG |> DisplayAs.Text
Due to the a priori calibration of the data, the gain phases are quite stable and just drift over time. Note that the one outlier around 2.5 UT is when ALMA leaves the array. As a result we shift the reference antenna to APEX on that scan causing the gain phases to jump.
We can also plot the gain amplitudes
julia
plotcaltable(abs.(intopt)) |> DisplayAs.PNG |> DisplayAs.Text
Here we find relatively stable gains for most stations. The exception is LMT which has a large offset after 2.5 UT. This is a known issue with the LMT in 2017 and is due to pointing issues. However, we can see the power of simultaneous imaging and instrument modeling as we are able to solve for the gain amplitudes and get a reasonable image.
One problem with the MAP estimate is that it does not provide uncertainties on the image. That is we are unable to statistically assess which components of the image are certain. Comrade is really a Bayesian imaging and calibration package for VLBI. Therefore, our goal is to sample from the posterior distribution of the image and instrument model. This is a very high-dimensional distribution with typically 1,000 - 100,000 parameters. To sample from this very high dimensional distribution, Comrade has an array of samplers that can be used. However, note that Comrade also satisfies the LogDensityProblems.jl interface. Therefore, you can use any package that supports LogDensityProblems.jl if you have your own fancy sampler.
For this example, we will use HMC, specifically the NUTS algorithm. For However, due to the need to sample a large number of gain parameters, constructing the posterior can take a few minutes. Therefore, for this tutorial, we will only do a quick preliminary run.
julia
using AdvancedHMC
chain = sample(rng, post, NUTS(0.8), 700; n_adapts = 500, initial_params = xopt)PosteriorSamples
Samples size: (700,)
sampler used: AHMC
Mean
┌──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky │ instrument │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64} │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}} │
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [-0.0249574 -0.0251727 … 0.0243484 0.0362277; -0.0301597 -0.0389172 … 0.0782692 0.0141458; … ; 0.0304332 0.074587 … -0.0109702 0.0304977; -0.0119619 0.0624222 … -0.0250868 0.0014627], hyperparams = 40.3174), σimg = 1.85952, fg = 0.32527) │ (lg = [0.0539112, 0.0367574, 0.00575793, 0.0272572, -0.273768, 0.160172, 0.00164779, 0.0294279, -0.0839277, 0.101536 … -0.0107824, 0.0133733, -0.658941, 0.0130774, 0.00859295, 0.0278123, -0.0753224, 0.0213296, -0.711777, 0.0118118], gp = [0.0, 1.19292, 0.0, -2.193, 2.05583, 1.21722, 0.0, -2.24218, 2.05209, 1.20651 … -1.10544, -0.185126, -1.99453, 1.48727, 0.0, -1.82959, -1.20862, -0.168047, -2.12359, 1.38308]) │
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
Std. Dev.
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky │ instrument │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64} │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}} │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [0.542464 0.545913 … 0.587245 0.58272; 0.577962 0.661412 … 0.645884 0.612747; … ; 0.592717 0.634832 … 0.62159 0.549095; 0.530158 0.606667 … 0.561158 0.534056], hyperparams = 23.0264), σimg = 0.206788, fg = 0.0880031) │ (lg = [0.15271, 0.150921, 0.0206144, 0.0224107, 0.0998478, 0.0707721, 0.0188178, 0.020438, 0.0933174, 0.0655478 … 0.0543661, 0.0194477, 0.105723, 0.0187268, 0.0182161, 0.0212383, 0.0568516, 0.0194572, 0.107184, 0.0173747], gp = [0.0, 1.78008, 0.0, 0.0169901, 1.02049, 1.67079, 0.0, 0.0177288, 1.15894, 1.54321 … 0.352495, 0.385118, 0.360936, 2.37464, 0.0, 0.0197237, 0.366945, 0.383723, 0.419036, 2.44346]) │
└─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘Note
The above sampler will store the samples in memory, i.e. RAM. For large models this can lead to out-of-memory issues. To avoid this we recommend using the saveto = DiskStore() kwargs which periodically saves the samples to disk limiting memory useage. You can load the chain using load_samples(diskout) where diskout is the object returned from sample.
Now we prune the adaptation phase
julia
chain = chain[(begin + 500):end]PosteriorSamples
Samples size: (200,)
sampler used: AHMC
Mean
┌────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky │ instrument │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64} │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}} │
├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [-0.0218183 0.0202816 … 0.0748636 0.0642358; 0.0237347 0.0229201 … 0.0887402 0.0307705; … ; 0.0144606 0.00232883 … -0.0923961 0.0158294; 0.0064829 -0.000490427 … -0.111853 -0.0123247], hyperparams = 47.5333), σimg = 1.95384, fg = 0.384945) │ (lg = [0.0482547, 0.0346024, 0.00450664, 0.0269286, -0.295181, 0.16846, 0.00017704, 0.0315785, -0.10584, 0.11938 … -0.0038557, 0.0139322, -0.732362, 0.0131829, 0.00967142, 0.0256086, -0.0660486, 0.0209722, -0.786968, 0.011531], gp = [0.0, 0.885165, 0.0, -2.19329, 1.87278, 0.747851, 0.0, -2.24237, 1.92797, 0.558443 … -1.45618, 0.192231, -2.37505, -0.383558, 0.0, -1.83002, -1.57635, 0.201507, -2.49644, -0.504677]) │
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
Std. Dev.
┌──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ sky │ instrument │
│ @NamedTuple{c::@NamedTuple{params::Matrix{Float64}, hyperparams::Float64}, σimg::Float64, fg::Float64} │ @NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Int64, Float64}}, Vector{Comrade.FrequencyChannel{Float64, Int64}}, Vector{Symbol}}} │
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ (c = (params = [0.505213 0.543242 … 0.607785 0.535246; 0.576354 0.592984 … 0.590116 0.622024; … ; 0.593724 0.629401 … 0.70789 0.556253; 0.521556 0.577515 … 0.563045 0.494608], hyperparams = 24.1244), σimg = 0.194952, fg = 0.0591083) │ (lg = [0.135707, 0.141914, 0.0190027, 0.0211961, 0.0808575, 0.0626726, 0.0196254, 0.0219053, 0.0691894, 0.062647 … 0.0463151, 0.0196659, 0.0666861, 0.0188531, 0.0185976, 0.0226825, 0.0517853, 0.018102, 0.0687197, 0.0178776], gp = [0.0, 0.541519, 0.0, 0.0146162, 0.317102, 0.548648, 0.0, 0.0171193, 0.313222, 0.557523 … 0.287571, 0.41182, 0.248492, 2.74225, 0.0, 0.0231788, 0.290906, 0.405344, 0.465701, 2.72856]) │
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘Warning
This should be run for likely an order of magnitude more steps to properly estimate expectations of the posterior
Now that we have our posterior, we can put error bars on all of our plots above. Let's start by finding the mean and standard deviation of the gain phases
julia
mchain = Comrade.rmap(mean, chain);
schain = Comrade.rmap(std, chain);Now we can use the measurements package to automatically plot everything with error bars. First we create a caltable the same way but making sure all of our variables have errors attached to them.
julia
using Measurements
gmeas = instrumentmodel(post, (; instrument = map((x, y) -> Measurements.measurement.(x, y), mchain.instrument, schain.instrument)))
ctable_am = caltable(abs.(gmeas))
ctable_ph = caltable(angle.(gmeas))─────────┬───────────────────┬────────────────────────────────────────────────────────────────────────────────────
Ti │ Fr │ AA AP AZ JC LM PV SM
─────────┼───────────────────┼────────────────────────────────────────────────────────────────────────────────────
0.92 hr │ 227.070703125 GHz │ 0.0±0.0 missing missing missing missing 0.88±0.54 missing
1.22 hr │ 227.070703125 GHz │ 0.0±0.0 -2.193±0.015 missing missing 1.87±0.32 0.75±0.55 missing
1.52 hr │ 227.070703125 GHz │ 0.0±0.0 -2.242±0.017 missing missing 1.93±0.31 0.56±0.56 missing
1.82 hr │ 227.070703125 GHz │ 0.0±0.0 -2.274±0.021 missing missing 1.96±0.31 0.4±0.56 missing
2.12 hr │ 227.070703125 GHz │ 0.0±0.0 -2.327±0.017 missing missing 2.01±0.31 0.22±0.57 missing
2.45 hr │ 227.070703125 GHz │ missing 0.0±0.0 missing missing 2.51±0.3 1.6±1.7 missing
2.75 hr │ 227.070703125 GHz │ 0.0±0.0 -2.403±0.018 missing missing 2.0±0.3 -0.13±0.57 missing
3.05 hr │ 227.070703125 GHz │ 0.0±0.0 -2.435±0.02 missing missing 1.99±0.29 -0.31±0.57 missing
3.35 hr │ 227.070703125 GHz │ 0.0±0.0 -2.464±0.018 missing missing 1.95±0.28 -0.51±0.57 missing
3.68 hr │ 227.070703125 GHz │ 0.0±0.0 -2.468±0.015 1.73±0.32 missing 1.89±0.28 -0.73±0.56 missing
3.98 hr │ 227.070703125 GHz │ 0.0±0.0 -2.437±0.018 1.58±0.32 missing 1.8±0.27 -0.96±0.55 missing
4.28 hr │ 227.070703125 GHz │ 0.0±0.0 -2.456±0.018 1.4±0.31 missing 1.67±0.27 -1.18±0.54 missing
4.58 hr │ 227.070703125 GHz │ 0.0±0.0 -2.371±0.02 1.22±0.3 -0.44±0.32 1.57±0.26 -1.43±0.53 missing
4.92 hr │ 227.070703125 GHz │ 0.0±0.0 -2.425±0.02 1.0±0.3 -0.32±0.33 1.4±0.25 -1.7±0.53 -1.57±0.33
5.18 hr │ 227.070703125 GHz │ 0.0±0.0 -2.384±0.02 0.8±0.3 -0.23±0.34 1.2±0.25 -1.9±0.52 -1.82±0.34
5.45 hr │ 227.070703125 GHz │ 0.0±0.0 -2.303±0.02 0.63±0.29 -0.13±0.35 0.96±0.25 -2.01±0.81 -2.09±0.35
5.72 hr │ 227.070703125 GHz │ 0.0±0.0 missing 0.4±0.29 -0.047±0.37 0.6±0.25 missing -2.34±0.37
6.05 hr │ 227.070703125 GHz │ 0.0±0.0 -2.224±0.016 0.093±0.28 0.052±0.38 -0.17±0.25 missing missing
6.32 hr │ 227.070703125 GHz │ 0.0±0.0 -2.164±0.018 -0.12±0.28 0.1±0.4 -0.86±0.25 missing -2.2±1.5
6.58 hr │ 227.070703125 GHz │ 0.0±0.0 missing -0.4±0.28 0.15±0.41 -1.3±0.25 missing -1.2±2.4
6.85 hr │ 227.070703125 GHz │ 0.0±0.0 -2.108±0.018 -0.7±0.28 0.19±0.42 -1.64±0.25 missing -0.67±2.7
7.18 hr │ 227.070703125 GHz │ 0.0±0.0 -2.01±0.017 -0.97±0.28 0.21±0.42 -1.97±0.25 missing -0.36±2.7
7.45 hr │ 227.070703125 GHz │ 0.0±0.0 -1.968±0.018 -1.22±0.29 0.21±0.42 -2.18±0.25 missing -0.34±2.7
7.72 hr │ 227.070703125 GHz │ 0.0±0.0 -1.889±0.017 -1.46±0.29 0.19±0.41 -2.38±0.25 missing -0.38±2.7
7.98 hr │ 227.070703125 GHz │ 0.0±0.0 -1.83±0.023 -1.58±0.29 0.2±0.4 -2.5±0.47 missing -0.5±2.7
─────────┴───────────────────┴────────────────────────────────────────────────────────────────────────────────────Now let's plot the phase curves
julia
plotcaltable(ctable_ph) |> DisplayAs.PNG |> DisplayAs.Text
and now the amplitude curves
julia
plotcaltable(ctable_am) |> DisplayAs.PNG |> DisplayAs.Text
Finally let's construct some representative image reconstructions.
julia
samples = skymodel.(Ref(post), chain[begin:5:end])
imgs = intensitymap.(samples, Ref(g))
mimg = mean(imgs)
simg = std(imgs)
fig = Figure(; size = (700, 700));
axs = [Axis(fig[i, j], xreversed = true, aspect = 1) for i in 1:2, j in 1:2]
image!(axs[1, 1], mimg, colormap = :afmhot); axs[1, 1].title = "Mean"
image!(axs[1, 2], simg ./ (max.(mimg, 1.0e-8)), colorrange = (0.0, 2.0), colormap = :afmhot);axs[1, 2].title = "Frac. Uncer."
image!(axs[2, 1], imgs[1], colormap = :afmhot);
image!(axs[2, 2], imgs[end], colormap = :afmhot);
hidedecorations!.(axs)
fig |> DisplayAs.PNG |> DisplayAs.Text
And viola, you have just finished making a preliminary image and instrument model reconstruction. In reality, you should run the sample step for many more MCMC steps to get a reliable estimate for the reconstructed image and instrument model parameters.
This page was generated using Literate.jl.