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/DimensionalData/KlTtO/CondaPkg.toml
CondaPkg Found dependencies: /home/runner/.julia/packages/PythonCall/L4cjh/CondaPkg.toml
CondaPkg Found dependencies: /home/runner/.julia/packages/Pyehtim/Fm109/CondaPkg.toml
CondaPkg Resolving changes
+ ehtim (pip)
+ libstdcxx-ng
+ numpy
+ numpy (pip)
+ openssl
+ pandas
+ python
+ setuptools (pip)
+ uv
+ xarray
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, >=3, <3.6"
│ uv = ">=0.4"
│ libstdcxx-ng = ">=3.4,<15.0"
│ pandas = "<2"
│ xarray = "*"
│ numpy = ">=1.24, <2.0"
│
│ [dependencies.python]
│ channel = "conda-forge"
│ build = "*cpython*"
│ version = ">=3.8,<4, >=3.6,<=3.10"
│
│ [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.
/home/runner/work/Comrade.jl/Comrade.jl/examples/intermediate/StokesIImaging/.CondaPkg/.pixi/envs/default/lib/python3.10/site-packages/ehtim/__init__.py:58: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resourcesFor reproducibility we use a stable random number genreator
julia
using StableRNGs
rng = StableRNG(12)StableRNGs.LehmerRNG(state=0x00000000000000000000000000000019)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 0x7fc5b33163e0>obs = ehtim.obsdata.load_uvfits("~/Dropbox (Smithsonian External)/M872021Project/Data/2021/CASA/e21e18/V4/M87_calibrated_b3.uvf+EVPA_rotation+netcal_10savg+flag.uvfits") 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 0x7fc5d3f2d270>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.
julia
function sky(θ, metadata)
(; fg, c, σimg) = θ
(; ftot, mimg) = metadata
# 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}())
x0, y0 = centroid(m)
# 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))
return shifted(m, -x0, -y0) + 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 = 32
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.69663253574863e-10, 4.69663253574863e-10, 32) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.69663253574863e-10, 4.69663253574863e-10, 32) ForwardOrdered Regular Points)
)Now we need to specify our image prior. For this work we will use a Gaussian Markov Random field prior 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
using VLBIImagePriors
using Distributions
fwhmfac = 2 * sqrt(2 * log(2))
mpr = modify(Gaussian(), Stretch(μas2rad(50.0) ./ fwhmfac))
mimg = intensitymap(mpr, grid)┌ 32×32 IntensityMap{Float64, 2} ┐
├────────────────────────────────┴─────────────────────────────────────── dims ┐
↓ X Sampled{Float64} LinRange{Float64}(-4.69663253574863e-10, 4.69663253574863e-10, 32) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.69663253574863e-10, 4.69663253574863e-10, 32) ForwardOrdered Regular Points
└──────────────────────────────────────────────────────────────────────────────┘
↓ → -4.69663e-10 -4.39362e-10 -4.09062e-10 -3.78761e-10 -3.4846e-10 -3.18159e-10 -2.87858e-10 -2.57557e-10 -2.27256e-10 -1.96956e-10 -1.66655e-10 -1.36354e-10 -1.06053e-10 -7.57521e-11 -4.54513e-11 -1.51504e-11 1.51504e-11 4.54513e-11 7.57521e-11 1.06053e-10 1.36354e-10 1.66655e-10 1.96956e-10 2.27256e-10 2.57557e-10 2.87858e-10 3.18159e-10 3.4846e-10 3.78761e-10 4.09062e-10 4.39362e-10 4.69663e-10
-4.69663e-10 1.25675e-11 4.60978e-11 1.55054e-10 4.78251e-10 1.3527e-9 3.50847e-9 8.34459e-9 1.81997e-8 3.63993e-8 6.67567e-8 1.12271e-7 1.73145e-7 2.44865e-7 3.1755e-7 3.77633e-7 4.11812e-7 4.11812e-7 3.77633e-7 3.1755e-7 2.44865e-7 1.73145e-7 1.12271e-7 6.67567e-8 3.63993e-8 1.81997e-8 8.34459e-9 3.50847e-9 1.3527e-9 4.78251e-10 1.55054e-10 4.60978e-11 1.25675e-11
-4.39362e-10 4.60978e-11 1.69087e-10 5.6874e-10 1.75423e-9 4.96172e-9 1.28691e-8 3.06081e-8 6.67567e-8 1.33513e-7 2.44865e-7 4.11812e-7 6.351e-7 8.98167e-7 1.16478e-6 1.38516e-6 1.51053e-6 1.51053e-6 1.38516e-6 1.16478e-6 8.98167e-7 6.351e-7 4.11812e-7 2.44865e-7 1.33513e-7 6.67567e-8 3.06081e-8 1.28691e-8 4.96172e-9 1.75423e-9 5.6874e-10 1.69087e-10 4.60978e-11
-4.09062e-10 1.55054e-10 5.6874e-10 1.913e-9 5.90051e-9 1.66892e-8 4.32864e-8 1.02953e-7 2.24542e-7 4.49084e-7 8.23623e-7 1.38516e-6 2.13621e-6 3.02106e-6 3.91783e-6 4.65912e-6 5.0808e-6 5.0808e-6 4.65912e-6 3.91783e-6 3.02106e-6 2.13621e-6 1.38516e-6 8.23623e-7 4.49084e-7 2.24542e-7 1.02953e-7 4.32864e-8 1.66892e-8 5.90051e-9 1.913e-9 5.6874e-10 1.55054e-10
-3.78761e-10 4.78251e-10 1.75423e-9 5.90051e-9 1.81997e-8 5.14764e-8 1.33513e-7 3.1755e-7 6.92582e-7 1.38516e-6 2.5404e-6 4.27243e-6 6.58898e-6 9.31823e-6 1.20843e-5 1.43707e-5 1.56713e-5 1.56713e-5 1.43707e-5 1.20843e-5 9.31823e-6 6.58898e-6 4.27243e-6 2.5404e-6 1.38516e-6 6.92582e-7 3.1755e-7 1.33513e-7 5.14764e-8 1.81997e-8 5.90051e-9 1.75423e-9 4.78251e-10
-3.4846e-10 1.3527e-9 4.96172e-9 1.66892e-8 5.14764e-8 1.45597e-7 3.77633e-7 8.98167e-7 1.95892e-6 3.91783e-6 7.18534e-6 1.20843e-5 1.86365e-5 2.63559e-5 3.41794e-5 4.06464e-5 4.43252e-5 4.43252e-5 4.06464e-5 3.41794e-5 2.63559e-5 1.86365e-5 1.20843e-5 7.18534e-6 3.91783e-6 1.95892e-6 8.98167e-7 3.77633e-7 1.45597e-7 5.14764e-8 1.66892e-8 4.96172e-9 1.3527e-9
-3.18159e-10 3.50847e-9 1.28691e-8 4.32864e-8 1.33513e-7 3.77633e-7 9.79458e-7 2.32956e-6 5.0808e-6 1.01616e-5 1.86365e-5 3.13427e-5 4.8337e-5 6.83588e-5 8.86505e-5 0.000105424 0.000114965 0.000114965 0.000105424 8.86505e-5 6.83588e-5 4.8337e-5 3.13427e-5 1.86365e-5 1.01616e-5 5.0808e-6 2.32956e-6 9.79458e-7 3.77633e-7 1.33513e-7 4.32864e-8 1.28691e-8 3.50847e-9
-2.87858e-10 8.34459e-9 3.06081e-8 1.02953e-7 3.1755e-7 8.98167e-7 2.32956e-6 5.54065e-6 1.20843e-5 2.41685e-5 4.43252e-5 7.45458e-5 0.000114965 0.000162586 0.000210847 0.000250741 0.000273435 0.000273435 0.000250741 0.000210847 0.000162586 0.000114965 7.45458e-5 4.43252e-5 2.41685e-5 1.20843e-5 5.54065e-6 2.32956e-6 8.98167e-7 3.1755e-7 1.02953e-7 3.06081e-8 8.34459e-9
-2.57557e-10 1.81997e-8 6.67567e-8 2.24542e-7 6.92582e-7 1.95892e-6 5.0808e-6 1.20843e-5 2.63559e-5 5.27119e-5 9.6674e-5 0.000162586 0.000250741 0.000354602 0.000459862 0.000546871 0.000596367 0.000596367 0.000546871 0.000459862 0.000354602 0.000250741 0.000162586 9.6674e-5 5.27119e-5 2.63559e-5 1.20843e-5 5.0808e-6 1.95892e-6 6.92582e-7 2.24542e-7 6.67567e-8 1.81997e-8
-2.27256e-10 3.63993e-8 1.33513e-7 4.49084e-7 1.38516e-6 3.91783e-6 1.01616e-5 2.41685e-5 5.27119e-5 0.000105424 0.000193348 0.000325171 0.000501483 0.000709204 0.000919723 0.00109374 0.00119273 0.00119273 0.00109374 0.000919723 0.000709204 0.000501483 0.000325171 0.000193348 0.000105424 5.27119e-5 2.41685e-5 1.01616e-5 3.91783e-6 1.38516e-6 4.49084e-7 1.33513e-7 3.63993e-8
-1.96956e-10 6.67567e-8 2.44865e-7 8.23623e-7 2.5404e-6 7.18534e-6 1.86365e-5 4.43252e-5 9.6674e-5 0.000193348 0.000354602 0.000596367 0.000919723 0.00130069 0.00168678 0.00200593 0.00218748 0.00218748 0.00200593 0.00168678 0.00130069 0.000919723 0.000596367 0.000354602 0.000193348 9.6674e-5 4.43252e-5 1.86365e-5 7.18534e-6 2.5404e-6 8.23623e-7 2.44865e-7 6.67567e-8
-1.66655e-10 1.12271e-7 4.11812e-7 1.38516e-6 4.27243e-6 1.20843e-5 3.13427e-5 7.45458e-5 0.000162586 0.000325171 0.000596367 0.00100297 0.00154678 0.00218748 0.00283681 0.00337356 0.00367889 0.00367889 0.00337356 0.00283681 0.00218748 0.00154678 0.00100297 0.000596367 0.000325171 0.000162586 7.45458e-5 3.13427e-5 1.20843e-5 4.27243e-6 1.38516e-6 4.11812e-7 1.12271e-7
-1.36354e-10 1.73145e-7 6.351e-7 2.13621e-6 6.58898e-6 1.86365e-5 4.8337e-5 0.000114965 0.000250741 0.000501483 0.000919723 0.00154678 0.00238547 0.00337356 0.00437497 0.00520274 0.00567363 0.00567363 0.00520274 0.00437497 0.00337356 0.00238547 0.00154678 0.000919723 0.000501483 0.000250741 0.000114965 4.8337e-5 1.86365e-5 6.58898e-6 2.13621e-6 6.351e-7 1.73145e-7
-1.06053e-10 2.44865e-7 8.98167e-7 3.02106e-6 9.31823e-6 2.63559e-5 6.83588e-5 0.000162586 0.000354602 0.000709204 0.00130069 0.00218748 0.00337356 0.00477093 0.00618714 0.00735779 0.00802372 0.00802372 0.00735779 0.00618714 0.00477093 0.00337356 0.00218748 0.00130069 0.000709204 0.000354602 0.000162586 6.83588e-5 2.63559e-5 9.31823e-6 3.02106e-6 8.98167e-7 2.44865e-7
-7.57521e-11 3.1755e-7 1.16478e-6 3.91783e-6 1.20843e-5 3.41794e-5 8.86505e-5 0.000210847 0.000459862 0.000919723 0.00168678 0.00283681 0.00437497 0.00618714 0.00802372 0.00954187 0.0104055 0.0104055 0.00954187 0.00802372 0.00618714 0.00437497 0.00283681 0.00168678 0.000919723 0.000459862 0.000210847 8.86505e-5 3.41794e-5 1.20843e-5 3.91783e-6 1.16478e-6 3.1755e-7
-4.54513e-11 3.77633e-7 1.38516e-6 4.65912e-6 1.43707e-5 4.06464e-5 0.000105424 0.000250741 0.000546871 0.00109374 0.00200593 0.00337356 0.00520274 0.00735779 0.00954187 0.0113473 0.0123743 0.0123743 0.0113473 0.00954187 0.00735779 0.00520274 0.00337356 0.00200593 0.00109374 0.000546871 0.000250741 0.000105424 4.06464e-5 1.43707e-5 4.65912e-6 1.38516e-6 3.77633e-7
-1.51504e-11 4.11812e-7 1.51053e-6 5.0808e-6 1.56713e-5 4.43252e-5 0.000114965 0.000273435 0.000596367 0.00119273 0.00218748 0.00367889 0.00567363 0.00802372 0.0104055 0.0123743 0.0134942 0.0134942 0.0123743 0.0104055 0.00802372 0.00567363 0.00367889 0.00218748 0.00119273 0.000596367 0.000273435 0.000114965 4.43252e-5 1.56713e-5 5.0808e-6 1.51053e-6 4.11812e-7
1.51504e-11 4.11812e-7 1.51053e-6 5.0808e-6 1.56713e-5 4.43252e-5 0.000114965 0.000273435 0.000596367 0.00119273 0.00218748 0.00367889 0.00567363 0.00802372 0.0104055 0.0123743 0.0134942 0.0134942 0.0123743 0.0104055 0.00802372 0.00567363 0.00367889 0.00218748 0.00119273 0.000596367 0.000273435 0.000114965 4.43252e-5 1.56713e-5 5.0808e-6 1.51053e-6 4.11812e-7
4.54513e-11 3.77633e-7 1.38516e-6 4.65912e-6 1.43707e-5 4.06464e-5 0.000105424 0.000250741 0.000546871 0.00109374 0.00200593 0.00337356 0.00520274 0.00735779 0.00954187 0.0113473 0.0123743 0.0123743 0.0113473 0.00954187 0.00735779 0.00520274 0.00337356 0.00200593 0.00109374 0.000546871 0.000250741 0.000105424 4.06464e-5 1.43707e-5 4.65912e-6 1.38516e-6 3.77633e-7
7.57521e-11 3.1755e-7 1.16478e-6 3.91783e-6 1.20843e-5 3.41794e-5 8.86505e-5 0.000210847 0.000459862 0.000919723 0.00168678 0.00283681 0.00437497 0.00618714 0.00802372 0.00954187 0.0104055 0.0104055 0.00954187 0.00802372 0.00618714 0.00437497 0.00283681 0.00168678 0.000919723 0.000459862 0.000210847 8.86505e-5 3.41794e-5 1.20843e-5 3.91783e-6 1.16478e-6 3.1755e-7
1.06053e-10 2.44865e-7 8.98167e-7 3.02106e-6 9.31823e-6 2.63559e-5 6.83588e-5 0.000162586 0.000354602 0.000709204 0.00130069 0.00218748 0.00337356 0.00477093 0.00618714 0.00735779 0.00802372 0.00802372 0.00735779 0.00618714 0.00477093 0.00337356 0.00218748 0.00130069 0.000709204 0.000354602 0.000162586 6.83588e-5 2.63559e-5 9.31823e-6 3.02106e-6 8.98167e-7 2.44865e-7
1.36354e-10 1.73145e-7 6.351e-7 2.13621e-6 6.58898e-6 1.86365e-5 4.8337e-5 0.000114965 0.000250741 0.000501483 0.000919723 0.00154678 0.00238547 0.00337356 0.00437497 0.00520274 0.00567363 0.00567363 0.00520274 0.00437497 0.00337356 0.00238547 0.00154678 0.000919723 0.000501483 0.000250741 0.000114965 4.8337e-5 1.86365e-5 6.58898e-6 2.13621e-6 6.351e-7 1.73145e-7
1.66655e-10 1.12271e-7 4.11812e-7 1.38516e-6 4.27243e-6 1.20843e-5 3.13427e-5 7.45458e-5 0.000162586 0.000325171 0.000596367 0.00100297 0.00154678 0.00218748 0.00283681 0.00337356 0.00367889 0.00367889 0.00337356 0.00283681 0.00218748 0.00154678 0.00100297 0.000596367 0.000325171 0.000162586 7.45458e-5 3.13427e-5 1.20843e-5 4.27243e-6 1.38516e-6 4.11812e-7 1.12271e-7
1.96956e-10 6.67567e-8 2.44865e-7 8.23623e-7 2.5404e-6 7.18534e-6 1.86365e-5 4.43252e-5 9.6674e-5 0.000193348 0.000354602 0.000596367 0.000919723 0.00130069 0.00168678 0.00200593 0.00218748 0.00218748 0.00200593 0.00168678 0.00130069 0.000919723 0.000596367 0.000354602 0.000193348 9.6674e-5 4.43252e-5 1.86365e-5 7.18534e-6 2.5404e-6 8.23623e-7 2.44865e-7 6.67567e-8
2.27256e-10 3.63993e-8 1.33513e-7 4.49084e-7 1.38516e-6 3.91783e-6 1.01616e-5 2.41685e-5 5.27119e-5 0.000105424 0.000193348 0.000325171 0.000501483 0.000709204 0.000919723 0.00109374 0.00119273 0.00119273 0.00109374 0.000919723 0.000709204 0.000501483 0.000325171 0.000193348 0.000105424 5.27119e-5 2.41685e-5 1.01616e-5 3.91783e-6 1.38516e-6 4.49084e-7 1.33513e-7 3.63993e-8
2.57557e-10 1.81997e-8 6.67567e-8 2.24542e-7 6.92582e-7 1.95892e-6 5.0808e-6 1.20843e-5 2.63559e-5 5.27119e-5 9.6674e-5 0.000162586 0.000250741 0.000354602 0.000459862 0.000546871 0.000596367 0.000596367 0.000546871 0.000459862 0.000354602 0.000250741 0.000162586 9.6674e-5 5.27119e-5 2.63559e-5 1.20843e-5 5.0808e-6 1.95892e-6 6.92582e-7 2.24542e-7 6.67567e-8 1.81997e-8
2.87858e-10 8.34459e-9 3.06081e-8 1.02953e-7 3.1755e-7 8.98167e-7 2.32956e-6 5.54065e-6 1.20843e-5 2.41685e-5 4.43252e-5 7.45458e-5 0.000114965 0.000162586 0.000210847 0.000250741 0.000273435 0.000273435 0.000250741 0.000210847 0.000162586 0.000114965 7.45458e-5 4.43252e-5 2.41685e-5 1.20843e-5 5.54065e-6 2.32956e-6 8.98167e-7 3.1755e-7 1.02953e-7 3.06081e-8 8.34459e-9
3.18159e-10 3.50847e-9 1.28691e-8 4.32864e-8 1.33513e-7 3.77633e-7 9.79458e-7 2.32956e-6 5.0808e-6 1.01616e-5 1.86365e-5 3.13427e-5 4.8337e-5 6.83588e-5 8.86505e-5 0.000105424 0.000114965 0.000114965 0.000105424 8.86505e-5 6.83588e-5 4.8337e-5 3.13427e-5 1.86365e-5 1.01616e-5 5.0808e-6 2.32956e-6 9.79458e-7 3.77633e-7 1.33513e-7 4.32864e-8 1.28691e-8 3.50847e-9
3.4846e-10 1.3527e-9 4.96172e-9 1.66892e-8 5.14764e-8 1.45597e-7 3.77633e-7 8.98167e-7 1.95892e-6 3.91783e-6 7.18534e-6 1.20843e-5 1.86365e-5 2.63559e-5 3.41794e-5 4.06464e-5 4.43252e-5 4.43252e-5 4.06464e-5 3.41794e-5 2.63559e-5 1.86365e-5 1.20843e-5 7.18534e-6 3.91783e-6 1.95892e-6 8.98167e-7 3.77633e-7 1.45597e-7 5.14764e-8 1.66892e-8 4.96172e-9 1.3527e-9
3.78761e-10 4.78251e-10 1.75423e-9 5.90051e-9 1.81997e-8 5.14764e-8 1.33513e-7 3.1755e-7 6.92582e-7 1.38516e-6 2.5404e-6 4.27243e-6 6.58898e-6 9.31823e-6 1.20843e-5 1.43707e-5 1.56713e-5 1.56713e-5 1.43707e-5 1.20843e-5 9.31823e-6 6.58898e-6 4.27243e-6 2.5404e-6 1.38516e-6 6.92582e-7 3.1755e-7 1.33513e-7 5.14764e-8 1.81997e-8 5.90051e-9 1.75423e-9 4.78251e-10
4.09062e-10 1.55054e-10 5.6874e-10 1.913e-9 5.90051e-9 1.66892e-8 4.32864e-8 1.02953e-7 2.24542e-7 4.49084e-7 8.23623e-7 1.38516e-6 2.13621e-6 3.02106e-6 3.91783e-6 4.65912e-6 5.0808e-6 5.0808e-6 4.65912e-6 3.91783e-6 3.02106e-6 2.13621e-6 1.38516e-6 8.23623e-7 4.49084e-7 2.24542e-7 1.02953e-7 4.32864e-8 1.66892e-8 5.90051e-9 1.913e-9 5.6874e-10 1.55054e-10
4.39362e-10 4.60978e-11 1.69087e-10 5.6874e-10 1.75423e-9 4.96172e-9 1.28691e-8 3.06081e-8 6.67567e-8 1.33513e-7 2.44865e-7 4.11812e-7 6.351e-7 8.98167e-7 1.16478e-6 1.38516e-6 1.51053e-6 1.51053e-6 1.38516e-6 1.16478e-6 8.98167e-7 6.351e-7 4.11812e-7 2.44865e-7 1.33513e-7 6.67567e-8 3.06081e-8 1.28691e-8 4.96172e-9 1.75423e-9 5.6874e-10 1.69087e-10 4.60978e-11
4.69663e-10 1.25675e-11 4.60978e-11 1.55054e-10 4.78251e-10 1.3527e-9 3.50847e-9 8.34459e-9 1.81997e-8 3.63993e-8 6.67567e-8 1.12271e-7 1.73145e-7 2.44865e-7 3.1755e-7 3.77633e-7 4.11812e-7 4.11812e-7 3.77633e-7 3.1755e-7 2.44865e-7 1.73145e-7 1.12271e-7 6.67567e-8 3.63993e-8 1.81997e-8 8.34459e-9 3.50847e-9 1.3527e-9 4.78251e-10 1.55054e-10 4.60978e-11 1.25675e-11Now we can form our metadata we need to fully define our model. We will also fix the total flux to be the observed value 1.1. This is because total flux is degenerate with a global shift in the gain amplitudes making the problem degenerate. To fix this we use the observed total flux as our value.
julia
skymeta = (; ftot = 1.1, mimg = mimg ./ flux(mimg))(ftot = 1.1, mimg = [1.2567547416996487e-11 4.60979671782088e-11 1.5505446140255663e-10 4.782533378056425e-10 1.3527047131498763e-9 3.5084819555512232e-9 8.344623408800353e-9 1.8199752706952592e-8 3.639950541390525e-8 6.675698727040275e-8 1.1227142257763898e-7 1.7314620328318453e-7 2.448657089564883e-7 3.175515369524363e-7 3.7763454712388683e-7 4.1181339376014135e-7 4.1181339376014135e-7 3.7763454712388683e-7 3.175515369524363e-7 2.448657089564883e-7 1.7314620328318453e-7 1.1227142257763922e-7 6.675698727040275e-8 3.639950541390525e-8 1.8199752706952592e-8 8.344623408800353e-9 3.5084819555512232e-9 1.3527047131498763e-9 4.782533378056425e-10 1.5505446140255717e-10 4.60979671782088e-11 1.2567547416996487e-11; 4.60979671782088e-11 1.690880891437348e-10 5.687422720922681e-10 1.7542409777756147e-9 4.9617427648818055e-9 1.2869168555004412e-8 3.060821361956108e-8 6.675698727040275e-8 1.335139745408057e-7 2.448657089564883e-7 4.1181339376014066e-7 6.351030739048726e-7 8.981713806211123e-7 1.1647841732449665e-6 1.3851696262654769e-6 1.510538188495548e-6 1.510538188495548e-6 1.3851696262654769e-6 1.1647841732449665e-6 8.981713806211123e-7 6.351030739048726e-7 4.1181339376014135e-7 2.448657089564883e-7 1.335139745408057e-7 6.675698727040275e-8 3.060821361956108e-8 1.2869168555004412e-8 4.9617427648818055e-9 1.7542409777756147e-9 5.687422720922701e-10 1.690880891437348e-10 4.60979671782088e-11; 1.5505446140255663e-10 5.687422720922681e-10 1.913013351222567e-9 5.900539798810739e-9 1.6689246817600684e-8 4.328655082079613e-8 1.0295334844003551e-7 2.2454284515527802e-7 4.4908569031055694e-7 8.23626787520283e-7 1.3851696262654769e-6 2.136223592652892e-6 3.021076376991096e-6 3.917851343303821e-6 4.659136692979867e-6 5.0808245912389815e-6 5.0808245912389815e-6 4.659136692979867e-6 3.917851343303814e-6 3.021076376991096e-6 2.136223592652892e-6 1.3851696262654794e-6 8.236267875202814e-7 4.4908569031055694e-7 2.2454284515527802e-7 1.0295334844003551e-7 4.328655082079613e-8 1.6689246817600684e-8 5.900539798810739e-9 1.9130133512225742e-9 5.687422720922681e-10 1.5505446140255663e-10; 4.782533378056425e-10 1.7542409777756147e-9 5.900539798810739e-9 1.819975270695265e-8 5.1476674220017655e-8 1.3351397454080597e-7 3.175515369524369e-7 6.925848131327384e-7 1.3851696262654794e-6 2.5404122956194908e-6 4.272447185305786e-6 6.589014300162271e-6 9.318273385959749e-6 1.2084305507964394e-5 1.4370742089937826e-5 1.5671405373215287e-5 1.5671405373215287e-5 1.4370742089937826e-5 1.2084305507964394e-5 9.318273385959749e-6 6.589014300162271e-6 4.272447185305794e-6 2.5404122956194908e-6 1.3851696262654794e-6 6.925848131327384e-7 3.175515369524369e-7 1.3351397454080597e-7 5.1476674220017655e-8 1.819975270695265e-8 5.90053979881076e-9 1.7542409777756147e-9 4.782533378056425e-10; 1.3527047131498763e-9 4.9617427648818055e-9 1.6689246817600684e-8 5.1476674220017655e-8 1.455980216556205e-7 3.7763454712388614e-7 8.981713806211139e-7 1.9589256716519037e-6 3.917851343303814e-6 7.185371044968894e-6 1.2084305507964365e-5 1.8636546771919454e-5 2.6356057200649016e-5 3.41795774824462e-5 4.0646596729911805e-5 4.432542804049519e-5 4.432542804049519e-5 4.0646596729911805e-5 3.41795774824462e-5 2.6356057200649016e-5 1.8636546771919454e-5 1.2084305507964374e-5 7.185371044968885e-6 3.917851343303814e-6 1.9589256716519037e-6 8.981713806211139e-7 3.7763454712388614e-7 1.455980216556205e-7 5.1476674220017655e-8 1.6689246817600743e-8 4.9617427648818055e-9 1.3527047131498763e-9; 3.5084819555512232e-9 1.2869168555004412e-8 4.328655082079613e-8 1.3351397454080597e-7 3.7763454712388614e-7 9.79462835825955e-7 2.329568346489937e-6 5.0808245912389815e-6 1.0161649182477983e-5 1.863654677191947e-5 3.134281074643055e-5 4.833722203185754e-5 6.835915496489247e-5 8.865085608099055e-5 0.00010542422880259616 0.00011496593671950256 0.00011496593671950256 0.00010542422880259628 8.865085608099055e-5 6.835915496489247e-5 4.833722203185754e-5 3.1342810746430575e-5 1.863654677191947e-5 1.0161649182477983e-5 5.0808245912389815e-6 2.329568346489937e-6 9.79462835825955e-7 3.7763454712388614e-7 1.3351397454080597e-7 4.3286550820796285e-8 1.2869168555004412e-8 3.5084819555512232e-9; 8.344623408800353e-9 3.060821361956108e-8 1.0295334844003551e-7 3.175515369524369e-7 8.981713806211139e-7 2.329568346489937e-6 5.540678505061918e-6 1.2084305507964385e-5 2.4168611015928815e-5 4.4325428040495274e-5 7.454618708767795e-5 0.00011496593671950265 0.0001625863869196475 0.0002108484576051926 0.0002507424859714445 0.00027343661985957014 0.00027343661985957014 0.0002507424859714447 0.0002108484576051926 0.0001625863869196475 0.00011496593671950265 7.454618708767802e-5 4.4325428040495274e-5 2.4168611015928815e-5 1.2084305507964385e-5 5.540678505061918e-6 2.329568346489937e-6 8.981713806211139e-7 3.175515369524369e-7 1.0295334844003588e-7 3.060821361956108e-8 8.344623408800353e-9; 1.8199752706952592e-8 6.675698727040275e-8 2.2454284515527802e-7 6.925848131327384e-7 1.9589256716519037e-6 5.0808245912389815e-6 1.2084305507964385e-5 2.6356057200648992e-5 5.271211440129808e-5 9.667444406371495e-5 0.00016258638691964722 0.00025074248597144427 0.0003546034243239617 0.0004598637468780097 0.0005468732397191394 0.0005963694967014225 0.0005963694967014225 0.0005468732397191394 0.0004598637468780094 0.0003546034243239617 0.00025074248597144427 0.0001625863869196475 9.667444406371485e-5 5.271211440129808e-5 2.6356057200648992e-5 1.2084305507964385e-5 5.0808245912389815e-6 1.9589256716519037e-6 6.925848131327384e-7 2.2454284515527882e-7 6.675698727040275e-8 1.8199752706952592e-8; 3.639950541390525e-8 1.335139745408057e-7 4.4908569031055694e-7 1.3851696262654794e-6 3.917851343303814e-6 1.0161649182477983e-5 2.4168611015928815e-5 5.271211440129808e-5 0.00010542422880259636 0.00019334888812743022 0.000325172773839295 0.0005014849719428892 0.0007092068486479246 0.0009197274937560213 0.0010937464794382808 0.0011927389934028472 0.0011927389934028472 0.0010937464794382808 0.0009197274937560206 0.0007092068486479246 0.0005014849719428892 0.0003251727738392957 0.00019334888812743006 0.00010542422880259636 5.271211440129808e-5 2.4168611015928815e-5 1.0161649182477983e-5 3.917851343303814e-6 1.3851696262654794e-6 4.490856903105585e-7 1.335139745408057e-7 3.639950541390525e-8; 6.675698727040275e-8 2.448657089564883e-7 8.23626787520283e-7 2.5404122956194908e-6 7.185371044968894e-6 1.863654677191947e-5 4.4325428040495274e-5 9.667444406371495e-5 0.00019334888812743022 0.0003546034243239617 0.0005963694967014229 0.0009197274937560198 0.0013006910953571786 0.001686787660841538 0.0020059398877715537 0.002187492958876558 0.002187492958876558 0.0020059398877715537 0.0016867876608415372 0.0013006910953571786 0.0009197274937560198 0.0005963694967014233 0.00035460342432396154 0.00019334888812743022 9.667444406371495e-5 4.4325428040495274e-5 1.863654677191947e-5 7.185371044968894e-6 2.5404122956194908e-6 8.23626787520286e-7 2.448657089564883e-7 6.675698727040275e-8; 1.1227142257763898e-7 4.1181339376014066e-7 1.3851696262654769e-6 4.272447185305786e-6 1.2084305507964365e-5 3.134281074643055e-5 7.454618708767795e-5 0.00016258638691964722 0.000325172773839295 0.0005963694967014229 0.0010029699438857767 0.0015467911050194411 0.002187492958876558 0.0028368273945916967 0.0033735753216830766 0.003678909975024079 0.003678909975024079 0.0033735753216830766 0.0028368273945916954 0.002187492958876558 0.0015467911050194411 0.001002969943885778 0.0005963694967014221 0.000325172773839295 0.00016258638691964722 7.454618708767795e-5 3.134281074643055e-5 1.2084305507964365e-5 4.272447185305786e-6 1.3851696262654818e-6 4.1181339376014066e-7 1.1227142257763898e-7; 1.7314620328318453e-7 6.351030739048726e-7 2.136223592652892e-6 6.589014300162271e-6 1.8636546771919454e-5 4.833722203185754e-5 0.00011496593671950265 0.00025074248597144427 0.0005014849719428892 0.0009197274937560198 0.0015467911050194411 0.002385477986805694 0.003373575321683077 0.004374985917753121 0.00520276438142872 0.005673654789183394 0.005673654789183394 0.00520276438142872 0.004374985917753119 0.003373575321683077 0.002385477986805694 0.0015467911050194428 0.0009197274937560193 0.0005014849719428892 0.00025074248597144427 0.00011496593671950265 4.833722203185754e-5 1.8636546771919454e-5 6.589014300162271e-6 2.1362235926529e-6 6.351030739048726e-7 1.7314620328318453e-7; 2.448657089564883e-7 8.981713806211123e-7 3.021076376991096e-6 9.318273385959749e-6 2.6356057200649016e-5 6.835915496489247e-5 0.0001625863869196475 0.0003546034243239617 0.0007092068486479246 0.0013006910953571786 0.002187492958876558 0.003373575321683077 0.004770955973611384 0.006187164420077763 0.007357819950048159 0.008023759551086217 0.008023759551086217 0.007357819950048159 0.006187164420077759 0.004770955973611384 0.003373575321683077 0.002187492958876561 0.001300691095357178 0.0007092068486479246 0.0003546034243239617 0.0001625863869196475 6.835915496489247e-5 2.6356057200649016e-5 9.318273385959749e-6 3.0210763769911065e-6 8.981713806211123e-7 2.448657089564883e-7; 3.175515369524363e-7 1.1647841732449665e-6 3.917851343303821e-6 1.2084305507964394e-5 3.41795774824462e-5 8.865085608099055e-5 0.0002108484576051926 0.0004598637468780097 0.0009197274937560213 0.001686787660841538 0.0028368273945916967 0.004374985917753121 0.006187164420077763 0.008023759551086224 0.009541911947222773 0.010405528762857438 0.010405528762857438 0.009541911947222773 0.008023759551086218 0.006187164420077763 0.004374985917753121 0.002836827394591701 0.001686787660841538 0.0009197274937560213 0.0004598637468780097 0.0002108484576051926 8.865085608099055e-5 3.41795774824462e-5 1.2084305507964394e-5 3.917851343303835e-6 1.1647841732449665e-6 3.175515369524363e-7; 3.7763454712388683e-7 1.3851696262654769e-6 4.659136692979867e-6 1.4370742089937826e-5 4.0646596729911805e-5 0.00010542422880259616 0.0002507424859714445 0.0005468732397191394 0.0010937464794382808 0.0020059398877715537 0.0033735753216830766 0.00520276438142872 0.007357819950048159 0.009541911947222773 0.011347309578366787 0.012374328840155524 0.012374328840155524 0.011347309578366787 0.009541911947222767 0.007357819950048159 0.00520276438142872 0.0033735753216830814 0.0020059398877715524 0.0010937464794382808 0.0005468732397191394 0.0002507424859714445 0.00010542422880259616 4.0646596729911805e-5 1.4370742089937826e-5 4.659136692979882e-6 1.3851696262654769e-6 3.7763454712388683e-7; 4.1181339376014135e-7 1.510538188495548e-6 5.0808245912389815e-6 1.5671405373215287e-5 4.432542804049519e-5 0.00011496593671950256 0.00027343661985957014 0.0005963694967014225 0.0011927389934028472 0.002187492958876558 0.003678909975024079 0.005673654789183394 0.008023759551086217 0.010405528762857438 0.012374328840155524 0.013494301286732307 0.013494301286732307 0.012374328840155525 0.010405528762857434 0.008023759551086217 0.005673654789183394 0.0036789099750240833 0.0021874929588765572 0.0011927389934028472 0.0005963694967014225 0.00027343661985957014 0.00011496593671950256 4.432542804049519e-5 1.5671405373215287e-5 5.080824591239e-6 1.510538188495548e-6 4.1181339376014135e-7; 4.1181339376014135e-7 1.510538188495548e-6 5.0808245912389815e-6 1.5671405373215287e-5 4.432542804049519e-5 0.00011496593671950256 0.00027343661985957014 0.0005963694967014225 0.0011927389934028472 0.002187492958876558 0.003678909975024079 0.005673654789183394 0.008023759551086217 0.010405528762857438 0.012374328840155524 0.013494301286732307 0.013494301286732307 0.012374328840155525 0.010405528762857434 0.008023759551086217 0.005673654789183394 0.0036789099750240833 0.0021874929588765572 0.0011927389934028472 0.0005963694967014225 0.00027343661985957014 0.00011496593671950256 4.432542804049519e-5 1.5671405373215287e-5 5.080824591239e-6 1.510538188495548e-6 4.1181339376014135e-7; 3.7763454712388683e-7 1.3851696262654769e-6 4.659136692979867e-6 1.4370742089937826e-5 4.0646596729911805e-5 0.00010542422880259628 0.0002507424859714447 0.0005468732397191394 0.0010937464794382808 0.0020059398877715537 0.0033735753216830766 0.00520276438142872 0.007357819950048159 0.009541911947222773 0.011347309578366787 0.012374328840155525 0.012374328840155525 0.011347309578366788 0.009541911947222769 0.007357819950048159 0.00520276438142872 0.0033735753216830814 0.0020059398877715524 0.0010937464794382808 0.0005468732397191394 0.0002507424859714447 0.00010542422880259628 4.0646596729911805e-5 1.4370742089937826e-5 4.659136692979882e-6 1.3851696262654769e-6 3.7763454712388683e-7; 3.175515369524363e-7 1.1647841732449665e-6 3.917851343303814e-6 1.2084305507964394e-5 3.41795774824462e-5 8.865085608099055e-5 0.0002108484576051926 0.0004598637468780094 0.0009197274937560206 0.0016867876608415372 0.0028368273945916954 0.004374985917753119 0.006187164420077759 0.008023759551086218 0.009541911947222767 0.010405528762857434 0.010405528762857434 0.009541911947222769 0.008023759551086213 0.006187164420077759 0.004374985917753119 0.0028368273945916997 0.0016867876608415364 0.0009197274937560206 0.0004598637468780094 0.0002108484576051926 8.865085608099055e-5 3.41795774824462e-5 1.2084305507964394e-5 3.917851343303829e-6 1.1647841732449665e-6 3.175515369524363e-7; 2.448657089564883e-7 8.981713806211123e-7 3.021076376991096e-6 9.318273385959749e-6 2.6356057200649016e-5 6.835915496489247e-5 0.0001625863869196475 0.0003546034243239617 0.0007092068486479246 0.0013006910953571786 0.002187492958876558 0.003373575321683077 0.004770955973611384 0.006187164420077763 0.007357819950048159 0.008023759551086217 0.008023759551086217 0.007357819950048159 0.006187164420077759 0.004770955973611384 0.003373575321683077 0.002187492958876561 0.001300691095357178 0.0007092068486479246 0.0003546034243239617 0.0001625863869196475 6.835915496489247e-5 2.6356057200649016e-5 9.318273385959749e-6 3.0210763769911065e-6 8.981713806211123e-7 2.448657089564883e-7; 1.7314620328318453e-7 6.351030739048726e-7 2.136223592652892e-6 6.589014300162271e-6 1.8636546771919454e-5 4.833722203185754e-5 0.00011496593671950265 0.00025074248597144427 0.0005014849719428892 0.0009197274937560198 0.0015467911050194411 0.002385477986805694 0.003373575321683077 0.004374985917753121 0.00520276438142872 0.005673654789183394 0.005673654789183394 0.00520276438142872 0.004374985917753119 0.003373575321683077 0.002385477986805694 0.0015467911050194428 0.0009197274937560193 0.0005014849719428892 0.00025074248597144427 0.00011496593671950265 4.833722203185754e-5 1.8636546771919454e-5 6.589014300162271e-6 2.1362235926529e-6 6.351030739048726e-7 1.7314620328318453e-7; 1.1227142257763922e-7 4.1181339376014135e-7 1.3851696262654794e-6 4.272447185305794e-6 1.2084305507964374e-5 3.1342810746430575e-5 7.454618708767802e-5 0.0001625863869196475 0.0003251727738392957 0.0005963694967014233 0.001002969943885778 0.0015467911050194428 0.002187492958876561 0.002836827394591701 0.0033735753216830814 0.0036789099750240833 0.0036789099750240833 0.0033735753216830814 0.0028368273945916997 0.002187492958876561 0.0015467911050194428 0.0010029699438857793 0.0005963694967014233 0.0003251727738392957 0.0001625863869196475 7.454618708767802e-5 3.1342810746430575e-5 1.2084305507964374e-5 4.272447185305794e-6 1.3851696262654843e-6 4.1181339376014135e-7 1.1227142257763922e-7; 6.675698727040275e-8 2.448657089564883e-7 8.236267875202814e-7 2.5404122956194908e-6 7.185371044968885e-6 1.863654677191947e-5 4.4325428040495274e-5 9.667444406371485e-5 0.00019334888812743006 0.00035460342432396154 0.0005963694967014221 0.0009197274937560193 0.001300691095357178 0.001686787660841538 0.0020059398877715524 0.0021874929588765572 0.0021874929588765572 0.0020059398877715524 0.0016867876608415364 0.001300691095357178 0.0009197274937560193 0.0005963694967014233 0.0003546034243239614 0.00019334888812743006 9.667444406371485e-5 4.4325428040495274e-5 1.863654677191947e-5 7.185371044968885e-6 2.5404122956194908e-6 8.236267875202844e-7 2.448657089564883e-7 6.675698727040275e-8; 3.639950541390525e-8 1.335139745408057e-7 4.4908569031055694e-7 1.3851696262654794e-6 3.917851343303814e-6 1.0161649182477983e-5 2.4168611015928815e-5 5.271211440129808e-5 0.00010542422880259636 0.00019334888812743022 0.000325172773839295 0.0005014849719428892 0.0007092068486479246 0.0009197274937560213 0.0010937464794382808 0.0011927389934028472 0.0011927389934028472 0.0010937464794382808 0.0009197274937560206 0.0007092068486479246 0.0005014849719428892 0.0003251727738392957 0.00019334888812743006 0.00010542422880259636 5.271211440129808e-5 2.4168611015928815e-5 1.0161649182477983e-5 3.917851343303814e-6 1.3851696262654794e-6 4.490856903105585e-7 1.335139745408057e-7 3.639950541390525e-8; 1.8199752706952592e-8 6.675698727040275e-8 2.2454284515527802e-7 6.925848131327384e-7 1.9589256716519037e-6 5.0808245912389815e-6 1.2084305507964385e-5 2.6356057200648992e-5 5.271211440129808e-5 9.667444406371495e-5 0.00016258638691964722 0.00025074248597144427 0.0003546034243239617 0.0004598637468780097 0.0005468732397191394 0.0005963694967014225 0.0005963694967014225 0.0005468732397191394 0.0004598637468780094 0.0003546034243239617 0.00025074248597144427 0.0001625863869196475 9.667444406371485e-5 5.271211440129808e-5 2.6356057200648992e-5 1.2084305507964385e-5 5.0808245912389815e-6 1.9589256716519037e-6 6.925848131327384e-7 2.2454284515527882e-7 6.675698727040275e-8 1.8199752706952592e-8; 8.344623408800353e-9 3.060821361956108e-8 1.0295334844003551e-7 3.175515369524369e-7 8.981713806211139e-7 2.329568346489937e-6 5.540678505061918e-6 1.2084305507964385e-5 2.4168611015928815e-5 4.4325428040495274e-5 7.454618708767795e-5 0.00011496593671950265 0.0001625863869196475 0.0002108484576051926 0.0002507424859714445 0.00027343661985957014 0.00027343661985957014 0.0002507424859714447 0.0002108484576051926 0.0001625863869196475 0.00011496593671950265 7.454618708767802e-5 4.4325428040495274e-5 2.4168611015928815e-5 1.2084305507964385e-5 5.540678505061918e-6 2.329568346489937e-6 8.981713806211139e-7 3.175515369524369e-7 1.0295334844003588e-7 3.060821361956108e-8 8.344623408800353e-9; 3.5084819555512232e-9 1.2869168555004412e-8 4.328655082079613e-8 1.3351397454080597e-7 3.7763454712388614e-7 9.79462835825955e-7 2.329568346489937e-6 5.0808245912389815e-6 1.0161649182477983e-5 1.863654677191947e-5 3.134281074643055e-5 4.833722203185754e-5 6.835915496489247e-5 8.865085608099055e-5 0.00010542422880259616 0.00011496593671950256 0.00011496593671950256 0.00010542422880259628 8.865085608099055e-5 6.835915496489247e-5 4.833722203185754e-5 3.1342810746430575e-5 1.863654677191947e-5 1.0161649182477983e-5 5.0808245912389815e-6 2.329568346489937e-6 9.79462835825955e-7 3.7763454712388614e-7 1.3351397454080597e-7 4.3286550820796285e-8 1.2869168555004412e-8 3.5084819555512232e-9; 1.3527047131498763e-9 4.9617427648818055e-9 1.6689246817600684e-8 5.1476674220017655e-8 1.455980216556205e-7 3.7763454712388614e-7 8.981713806211139e-7 1.9589256716519037e-6 3.917851343303814e-6 7.185371044968894e-6 1.2084305507964365e-5 1.8636546771919454e-5 2.6356057200649016e-5 3.41795774824462e-5 4.0646596729911805e-5 4.432542804049519e-5 4.432542804049519e-5 4.0646596729911805e-5 3.41795774824462e-5 2.6356057200649016e-5 1.8636546771919454e-5 1.2084305507964374e-5 7.185371044968885e-6 3.917851343303814e-6 1.9589256716519037e-6 8.981713806211139e-7 3.7763454712388614e-7 1.455980216556205e-7 5.1476674220017655e-8 1.6689246817600743e-8 4.9617427648818055e-9 1.3527047131498763e-9; 4.782533378056425e-10 1.7542409777756147e-9 5.900539798810739e-9 1.819975270695265e-8 5.1476674220017655e-8 1.3351397454080597e-7 3.175515369524369e-7 6.925848131327384e-7 1.3851696262654794e-6 2.5404122956194908e-6 4.272447185305786e-6 6.589014300162271e-6 9.318273385959749e-6 1.2084305507964394e-5 1.4370742089937826e-5 1.5671405373215287e-5 1.5671405373215287e-5 1.4370742089937826e-5 1.2084305507964394e-5 9.318273385959749e-6 6.589014300162271e-6 4.272447185305794e-6 2.5404122956194908e-6 1.3851696262654794e-6 6.925848131327384e-7 3.175515369524369e-7 1.3351397454080597e-7 5.1476674220017655e-8 1.819975270695265e-8 5.90053979881076e-9 1.7542409777756147e-9 4.782533378056425e-10; 1.5505446140255717e-10 5.687422720922701e-10 1.9130133512225742e-9 5.90053979881076e-9 1.6689246817600743e-8 4.3286550820796285e-8 1.0295334844003588e-7 2.2454284515527882e-7 4.490856903105585e-7 8.23626787520286e-7 1.3851696262654818e-6 2.1362235926529e-6 3.0210763769911065e-6 3.917851343303835e-6 4.659136692979882e-6 5.080824591239e-6 5.080824591239e-6 4.659136692979882e-6 3.917851343303829e-6 3.0210763769911065e-6 2.1362235926529e-6 1.3851696262654843e-6 8.236267875202844e-7 4.490856903105585e-7 2.2454284515527882e-7 1.0295334844003588e-7 4.3286550820796285e-8 1.6689246817600743e-8 5.90053979881076e-9 1.913013351222581e-9 5.687422720922701e-10 1.5505446140255717e-10; 4.60979671782088e-11 1.690880891437348e-10 5.687422720922681e-10 1.7542409777756147e-9 4.9617427648818055e-9 1.2869168555004412e-8 3.060821361956108e-8 6.675698727040275e-8 1.335139745408057e-7 2.448657089564883e-7 4.1181339376014066e-7 6.351030739048726e-7 8.981713806211123e-7 1.1647841732449665e-6 1.3851696262654769e-6 1.510538188495548e-6 1.510538188495548e-6 1.3851696262654769e-6 1.1647841732449665e-6 8.981713806211123e-7 6.351030739048726e-7 4.1181339376014135e-7 2.448657089564883e-7 1.335139745408057e-7 6.675698727040275e-8 3.060821361956108e-8 1.2869168555004412e-8 4.9617427648818055e-9 1.7542409777756147e-9 5.687422720922701e-10 1.690880891437348e-10 4.60979671782088e-11; 1.2567547416996487e-11 4.60979671782088e-11 1.5505446140255663e-10 4.782533378056425e-10 1.3527047131498763e-9 3.5084819555512232e-9 8.344623408800353e-9 1.8199752706952592e-8 3.639950541390525e-8 6.675698727040275e-8 1.1227142257763898e-7 1.7314620328318453e-7 2.448657089564883e-7 3.175515369524363e-7 3.7763454712388683e-7 4.1181339376014135e-7 4.1181339376014135e-7 3.7763454712388683e-7 3.175515369524363e-7 2.448657089564883e-7 1.7314620328318453e-7 1.1227142257763922e-7 6.675698727040275e-8 3.639950541390525e-8 1.8199752706952592e-8 8.344623408800353e-9 3.5084819555512232e-9 1.3527047131498763e-9 4.782533378056425e-10 1.5505446140255717e-10 4.60979671782088e-11 1.2567547416996487e-11])To 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: (32, 32)
)
) hyper prior:
Truncated(Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.054358121593522137)
θ: 18.396515013483654
)
; lower=1.0, upper=64.0)
)Putting everything together the total prior is then our image prior, a prior on the standard deviation of the MRF, and a prior on the fractional flux of the Gaussian component.
julia
prior = (
c = cprior,
σimg = truncated(Normal(0.0, 0.5); lower = 0.0),
fg = Uniform(0.0, 1.0),
)(c = HierarchicalPrior(
map:
ConditionalMarkov(
Random Field: VLBIImagePriors.GaussMarkovRandomField
Graph: MarkovRandomFieldGraph{1}(
dims: (32, 32)
)
) hyper prior:
Truncated(Distributions.InverseGamma{Float64}(
invd: Distributions.Gamma{Float64}(α=1.0, θ=0.054358121593522137)
θ: 18.396515013483654
)
; lower=1.0, upper=64.0)
)
, σimg = Truncated(Distributions.Normal{Float64}(μ=0.0, σ=0.5); lower=0.0), fg = Distributions.Uniform{Float64}(a=0.0, b=1.0))Now we can construct our sky model.
julia
skym = SkyModel(sky, prior, grid; metadata = skymeta)SkyModel
with map: sky
on grid:
RectiGrid(
executor: ComradeBase.Serial()
Dimensions:
(↓ X Sampled{Float64} LinRange{Float64}(-4.69663253574863e-10, 4.69663253574863e-10, 32) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.69663253574863e-10, 4.69663253574863e-10, 32) 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
G = SingleStokesGain() do x
lg = x.lg
gp = x.gp
return exp(lg + 1im * gp)
end
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. Additionally, since we want to use gradients we need to specify the AD mode. Essentially for all modes we recommend using Enzyme.set_runtime_activity(Enzyme.Reverse). Eventually as Comrade and Enzyme matures we will no need set_runtime_activity.
julia
using Enzyme
post = VLBIPosterior(skym, intmodel, dvis; admode = set_runtime_activity(Enzyme.Reverse))VLBIPosterior
ObservedSkyModel
with map: 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.69663253574863e-10, 4.69663253574863e-10, 32) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.69663253574863e-10, 4.69663253574863e-10, 32) ForwardOrdered Regular Points)
)
Visibility Domain: UnstructredDomain(
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.EHTVisibilityDatumdone using the asflat function.
julia
tpost = asflat(post)
ndim = dimension(tpost)1355We can now also find the dimension of our posterior or the number of parameters we are going to sample.
Warning
This can often be different from what you would expect. This is especially true when using angular variables where we often artificially increase the dimension of the parameter space to make sampling easier.
To initialize our sampler we will use optimize using Adam
julia
using Optimization
using OptimizationOptimisers
xopt, sol = comrade_opt(post, Optimisers.Adam(); initial_params = prior_sample(rng, post), maxiters = 20_000, 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. To do this, Comrade provides a caltable function that converts the flattened gain parameters to a tabular format based on the time and its segmentation.
julia
intopt = instrumentmodel(post, xopt)
gt = Comrade.caltable(angle.(intopt))
plotcaltable(gt) |> DisplayAs.PNG |> DisplayAs.Text
The gain phases are pretty random, although much of this is due to us picking a random reference sites for each scan.
Moving onto the gain amplitudes, we see that most of the gain variation is within 10% as expected except LMT, which has massive variations.
julia
gt = Comrade.caltable(abs.(intopt))
plotcaltable(gt) |> DisplayAs.PNG |> DisplayAs.Text
To sample from the posterior, we will use HMC, specifically the NUTS algorithm. For information about NUTS, see Michael Betancourt's notes. However, due to the need to sample a large number of gain parameters, constructing the posterior is rather time-consuming. Therefore, for this tutorial, we will only do a quick preliminary run
julia
using AdvancedHMC
chain = sample(rng, post, NUTS(0.8), 1_000; n_adapts = 500, progress = false, initial_params = xopt)PosteriorSamples
Samples size: (1000,)
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.00301803 0.0251345 … 0.0340956 0.0432541; 0.0197188 0.0106893 … 0.0608393 0.0580864; … ; -0.045135 -0.0685784 … 0.030086 0.0209466; -0.00239971 -0.0250522 … 0.0255903 -0.00555539], hyperparams = 25.1021), σimg = 2.14328, fg = 0.41268) │ (lg = [0.0232339, 0.0295758, 0.00744882, 0.0270675, -0.330611, 0.145062, 0.00380032, 0.0290773, -0.150635, 0.104457 … -0.0497594, 0.0140417, -0.824617, 0.0148191, 0.00853384, 0.0275933, -0.115967, 0.0224276, -0.870169, 0.011405], gp = [0.0, -0.669657, 0.0, -2.19186, 1.10084, -0.830032, 0.0, -2.24137, 1.16631, -1.03397 … -2.18596, 0.0800381, -1.16817, -0.415858, 0.0, -1.83109, -2.31458, 0.104126, 0.467074, -0.551113]) │
└──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
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.557979 0.647511 … 0.615322 0.557742; 0.556651 0.641241 … 0.62868 0.5629; … ; 0.58708 0.668228 … 0.662141 0.574723; 0.535372 0.575473 … 0.561484 0.511736], hyperparams = 15.2474), σimg = 0.283728, fg = 0.0671461) │ (lg = [0.147382, 0.145279, 0.0200271, 0.0217664, 0.094113, 0.0627787, 0.0201278, 0.0205804, 0.0857544, 0.0621377 … 0.058894, 0.0182279, 0.0581472, 0.0168254, 0.0170794, 0.0208763, 0.059951, 0.0186828, 0.0594319, 0.0160731], gp = [0.0, 0.446905, 0.0, 0.0182506, 0.245196, 0.450281, 0.0, 0.0177678, 0.247418, 0.455295 … 0.26575, 0.313978, 2.70395, 2.8514, 0.0, 0.0188364, 0.263533, 0.310172, 2.93748, 2.82886]) │
└───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘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 fix that you can include the keyword argument saveto = DiskStore() 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[501:end]PosteriorSamples
Samples size: (500,)
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.0318392 -0.017169 … 0.00353679 0.0380453; -0.0263066 -0.0720655 … 0.00283436 0.00986243; … ; -0.0414745 -0.0748675 … -0.000838874 0.0262174; -0.0108099 -0.0520351 … 0.00359707 -0.00221917], hyperparams = 24.3082), σimg = 2.02729, fg = 0.394152) │ (lg = [0.0202753, 0.0334072, 0.00709598, 0.0265342, -0.351078, 0.147696, 0.00290538, 0.0286557, -0.169272, 0.108936 … -0.0630879, 0.0136652, -0.842809, 0.0155061, 0.00917717, 0.0269159, -0.129585, 0.021719, -0.888442, 0.012459], gp = [0.0, -0.571426, 0.0, -2.19218, 1.07015, -0.733226, 0.0, -2.24108, 1.13101, -0.937959 … -2.23417, -0.06021, -0.938463, 0.37406, 0.0, -1.83106, -2.36145, -0.034974, 0.762532, 0.120358]) │
└─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
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.58302 0.641178 … 0.608284 0.556362; 0.559814 0.646989 … 0.648767 0.586031; … ; 0.579242 0.669383 … 0.650686 0.601027; 0.544944 0.592417 … 0.570602 0.511895], hyperparams = 15.2818), σimg = 0.256083, fg = 0.066755) │ (lg = [0.140254, 0.146408, 0.0190825, 0.0205538, 0.0803843, 0.0620398, 0.0194138, 0.0202098, 0.0765317, 0.0616648 … 0.0557461, 0.0183508, 0.0571468, 0.0166556, 0.0176777, 0.0204224, 0.056692, 0.0192449, 0.0578287, 0.0167508], gp = [0.0, 0.332873, 0.0, 0.0188131, 0.177673, 0.336652, 0.0, 0.0174729, 0.178169, 0.340937 … 0.195439, 0.19986, 2.83573, 2.95827, 0.0, 0.0178702, 0.191137, 0.197573, 2.89838, 2.98471]) │
└─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘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.57±0.33 missing
1.22 hr │ 227.070703125 GHz │ 0.0±0.0 -2.192±0.019 missing missing 1.07±0.18 -0.73±0.34 missing
1.52 hr │ 227.070703125 GHz │ 0.0±0.0 -2.241±0.017 missing missing 1.13±0.18 -0.94±0.34 missing
1.82 hr │ 227.070703125 GHz │ 0.0±0.0 -2.275±0.016 missing missing 1.16±0.18 -1.12±0.35 missing
2.12 hr │ 227.070703125 GHz │ 0.0±0.0 -2.328±0.016 missing missing 1.22±0.18 -1.31±0.35 missing
2.45 hr │ 227.070703125 GHz │ missing 0.0±0.0 missing missing 1.73±0.19 0.88±0.36 missing
2.75 hr │ 227.070703125 GHz │ 0.0±0.0 -2.401±0.019 missing missing 1.23±0.18 -1.68±0.35 missing
3.05 hr │ 227.070703125 GHz │ 0.0±0.0 -2.435±0.018 missing missing 1.23±0.18 -1.86±0.36 missing
3.35 hr │ 227.070703125 GHz │ 0.0±0.0 -2.465±0.017 missing missing 1.21±0.18 -2.06±0.36 missing
3.68 hr │ 227.070703125 GHz │ 0.0±0.0 -2.469±0.018 0.86±0.2 missing 1.17±0.18 -2.24±0.54 missing
3.98 hr │ 227.070703125 GHz │ 0.0±0.0 -2.438±0.019 0.72±0.2 missing 1.09±0.18 -2.4±0.73 missing
4.28 hr │ 227.070703125 GHz │ 0.0±0.0 -2.455±0.02 0.55±0.2 missing 0.97±0.18 -2.2±1.5 missing
4.58 hr │ 227.070703125 GHz │ 0.0±0.0 -2.373±0.021 0.38±0.2 -1.01±0.2 0.89±0.18 -0.8±2.7 missing
4.92 hr │ 227.070703125 GHz │ 0.0±0.0 -2.424±0.019 0.17±0.2 -0.85±0.2 0.73±0.19 0.3±2.9 -2.1±0.2
5.18 hr │ 227.070703125 GHz │ 0.0±0.0 -2.385±0.019 -0.009±0.2 -0.73±0.21 0.54±0.19 1.4±2.4 -2.32±0.21
5.45 hr │ 227.070703125 GHz │ 0.0±0.0 -2.302±0.017 -0.17±0.2 -0.59±0.21 0.32±0.19 1.9±1.8 -2.55±0.21
5.72 hr │ 227.070703125 GHz │ 0.0±0.0 missing -0.4±0.2 -0.48±0.21 -0.036±0.19 missing -2.5±1.2
6.05 hr │ 227.070703125 GHz │ 0.0±0.0 -2.224±0.02 -0.69±0.2 -0.34±0.22 -0.8±0.19 missing missing
6.32 hr │ 227.070703125 GHz │ 0.0±0.0 -2.163±0.019 -0.9±0.2 -0.26±0.22 -1.49±0.19 missing -1.4±2.6
6.58 hr │ 227.070703125 GHz │ 0.0±0.0 missing -1.17±0.2 -0.19±0.22 -1.93±0.19 missing -0.81±2.8
6.85 hr │ 227.070703125 GHz │ 0.0±0.0 -2.109±0.022 -1.47±0.2 -0.13±0.21 -2.27±0.19 missing 0.36±2.9
7.18 hr │ 227.070703125 GHz │ 0.0±0.0 -2.009±0.018 -1.74±0.2 -0.083±0.21 -2.61±0.18 missing 0.52±2.9
7.45 hr │ 227.070703125 GHz │ 0.0±0.0 -1.969±0.018 -2.0±0.2 -0.06±0.21 -2.79±0.49 missing 0.8±2.9
7.72 hr │ 227.070703125 GHz │ 0.0±0.0 -1.893±0.019 -2.23±0.2 -0.06±0.2 -0.94±2.8 missing 0.37±3.0
7.98 hr │ 227.070703125 GHz │ 0.0±0.0 -1.831±0.018 -2.36±0.19 -0.035±0.2 0.76±2.9 missing 0.12±3.0
─────────┴───────────────────┴─────────────────────────────────────────────────────────────────────────────────────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(; resolution = (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 = "Std"
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.