VLBISkyRegularizers.jl

Appearance

Full Visibilities Imaging ​

Status `~/work/VLBISkyRegularizers.jl/VLBISkyRegularizers.jl/examples/Project.toml`
⌃ [13f3f980] CairoMakie v0.12.5
⌃ [99d987ce] Comrade v0.10.4
  [0b91fe84] DisplayAs v0.1.6
⌃ [31c24e10] Distributions v0.25.109
⌃ [7da242da] Enzyme v0.12.31
⌃ [7f7a1694] Optimization v3.27.0
  [42dfb2eb] OptimizationOptimisers v0.2.1
  [3d61700d] Pyehtim v0.1.3
  [860ef19b] StableRNGs v1.0.2
  [2913bbd2] StatsBase v0.34.3
⌃ [fd094767] Suppressor v0.2.7
  [272825ce] VLBISkyRegularizers v1.0.0-DEV `~/work/VLBISkyRegularizers.jl/VLBISkyRegularizers.jl`
  [b77e0a4c] InteractiveUtils
  [56ddb016] Logging
Info Packages marked with ⌃ have new versions available and may be upgradable.
nothing

Loading the Data ​

We start by loading Comrade and VLBISkyRegularizers.

julia
using Comrade
using VLBISkyRegularizers

We use Pyehtim to load the uvfits files.

julia
using Pyehtim
┌ Warning: Module SciMLBase with build ID fafbfcfd-c24d-fca5-0000-004b1aa0cd12 is missing from the cache.
│ This may mean SciMLBase [0bca4576-84f4-4d90-8ffe-ffa030f20462] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:1948
┌ Warning: Module BangBang with build ID fafbfcfd-7e12-50d4-0000-00459fb53c4a is missing from the cache.
│ This may mean BangBang [198e06fe-97b7-11e9-32a5-e1d131e6ad66] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:1948
┌ Warning: Module Comrade with build ID fafbfcfd-f82d-34a3-0000-007be2c79b20 is missing from the cache.
│ This may mean Comrade [99d987ce-9a1e-4df8-bc0b-1ea019aa547b] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:1948

We also use a stable random number generator

julia
using StableRNGs
rng = StableRNG(100)
StableRNGs.LehmerRNG(state=0x000000000000000000000000000000c9)

Load the observation data with Pyehtim and do some preprocessing

julia
file = joinpath(__DIR, "data", "L2V1_M87_2018_111_b3_hops_netcal_10s_StokesI.uvfits")
obseht = ehtim.obsdata.load_uvfits(file)
obs = Pyehtim.scan_average(obseht).add_fractional_noise(0.02)
Python: <ehtim.obsdata.Obsdata object at 0x7fdb612c1390>

Extract the visibilities from the obsdata object

julia
dvis = extract_table(obs, Visibilities())
EHTObservationTable{Comrade.EHTVisibilityDatum{:I}}
  source:      M87
  mjd:         58228
  bandwidth:   1.856e9
  sites:       [:AA, :AX, :GL, :LM, :MG, :MM, :PV, :SW]
  nsamples:    339

Build sky model ​

We build a model consisting of a raster of pixels and an extended Gaussian component. The pixel raster is convolved with a pulse to make a ContinuousImage object. The Gaussian component is used to model the zero baselines. The sky model consists of two components, θ and metadata. θ is a NamedTuple with parameters we want to fit. impixel here is the pixel raster, and fg controls the relative scaling between the raster and the Gaussian component. metadata is also a NamedTuple with constant parameters. ftot controls the total flux scaling, pulse is the pulse with which the pixel raster is convolved, and regularizers is the set of regularizers we will be using.

julia
function sky(θ, metadata)
    (; impixel, fg) = θ
    (; ftot, pulse, regs) = metadata

    # Transform image to the image simplex domain and scale
    rast = ftot*(1-fg).*transform_image(image_domain(regs), impixel)
	# Form the image
    img = IntensityMap(rast, VLBISkyRegularizers.grid(regs))
    m = ContinuousImage(img, pulse)
    # Add an extended Gaussian
    gauss = modify(Gaussian(), Stretch(μas2rad(250.0)), Renormalize(fg*ftot))
    return m + gauss
end
sky (generic function with 1 method)

Define Image and Regularizers ​

Next, we define the grid on which our image model will exist.

julia
global fov = 150
global npix = 32
grid_p = imagepixels(μas2rad(fov), μas2rad(fov), npix, npix)
RectiGrid(
executor: ComradeBase.Serial()
Dimensions: 
(↓ X Sampled{Float64} LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-3.5224744018114725e-10, 3.5224744018114725e-10, 32) ForwardOrdered Regular Points)
)

Now we define the regularizers we want to use. The available regulariers are:

julia
subtypes(AbstractRegularizer)
8-element Vector{Any}:
 AddRegularizer
 KLE
 L1
 PSDistance
 TSV
 TV
 WaveletL1
 WeightRegularizer

Each regularizer must be passed with four arguments: the hyperparameter, the image domain, the evaluation domain, and the grid. The hyperparameter defines the weight of the regularizer, the image domain defines the domain to which the image will be transformed, the evaluation domain defines the domain on which the regularizer will be evaluated, and the grid defines the grid of pixels on which the regularizer will be evaluated. The available domains are:

julia
subtypes(AbstractDomain)
7-element Vector{Any}:
 ALRDomain
 CLRDomain
 LinearDomain
 LogDomain
 ParameterDomain
 PolarizationDomain
 VLBISkyRegularizers.PSDDomain

Note

When using multiple regularizers, the image domains and grids of all regularizers must match.

We will use a TSV regularizer in the log domain and a wavelet-space L1-norm regularizer in the linear domain, and have a log image domain.

julia
r1 = TSV(1, LinearDomain(), LogDomain(), grid_p)
r2 = WaveletL1(1, LinearDomain(), ParameterDomain(), grid_p)
regs = r1 + r2
Added Regularizer:
   Regularizer 1:
      Regularizer:          TSV
      Hyperparameter:       1
      Evaluation Domain:    LogDomain()
   Regularizer 2:
      Regularizer:          Wavelet L1
      Hyperparameter:       1
      Evaluation Domain:    ParameterDomain()
      Wavelet:              Full db2
   Image Domain:   LinearDomain()
   Grid FOV:             150.0x150.0 μas
   Grid Size:            32x32

Build Posterior ​

Now we build our sky model.

julia
using Distributions
skymeta = (;ftot = 1.1, pulse = BSplinePulse{3}(), regs = regs)
skyprior = (
    impixel=regs,
    fg=Uniform(0,1)
)
skymodel = SkyModel(sky, skyprior, grid_p, metadata=skymeta)
SkyModel
  with map: sky
   on grid: ComradeBase.RectiGrid

And next our intstrument model.

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()

Now we can make our posterior.

julia
post = VLBIPosterior(skymodel, intmodel, dvis)
VLBIPosterior
ObservedSkyModel
  with map: sky
   on grid: VLBISkyModels.FourierDualDomainObservedInstrumentModel
  with Jones: SingleStokesGain
  with reference basis: PolarizedTypes.CirBasis()Data Products: Comrade.EHTVisibilityDatum

Image Reconstruction ​

We use the solve_opt function to solve for the MAP estimate. By default, this uses an Adam optimizer from the OptimizationOptimisers.jl package and automatic differentiation with Enzyme from the Optimization.jl package. This function runs the optimization ntrials number of times. Each optimization trial first runs half of maxiters number of iterations, adds some random noise, then runs maxiters iterations from that point. The output values are sorted by order of objective value.

julia
using Optimization
using OptimizationOptimisers
using Enzyme
solve_opt(post, Optimisers.Adam(), Optimization.AutoEnzyme(); ntrials=5, maxiters=10_000, verbose=false)
(@NamedTuple{sky::@NamedTuple{impixel::Matrix{Float64}, fg::Float64}, instrument::@NamedTuple{lg::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Float64}}, FillArrays.Fill{Float64, 1, Tuple{Base.OneTo{Int64}}}, Vector{Symbol}}, gp::Comrade.SiteArray{Float64, 1, Vector{Float64}, Vector{Comrade.IntegrationTime{Float64}}, FillArrays.Fill{Float64, 1, Tuple{Base.OneTo{Int64}}}, Vector{Symbol}}}}[(sky = (impixel = [0.014116881532772881 0.014777348727080645 … 0.01928854389001646 0.01669730603028104; 0.01505454278668511 0.01576975748986589 … 0.019111485626086586 0.0161713022013274; … ; 0.042272395180912525 0.04278591819471418 … 0.04124726158080567 0.034855794542158335; 0.03722135770597106 0.03816550957356247 … 0.03658272823542657 3.0729314908842665e-5], fg = 0.010141558999722994), instrument = (lg = [0.04489238204647651, -0.04606488786268801, 0.3470290050730102, -0.022038531834720758, 0.02262126731612938, 0.4036955614195387, 0.04122279280137001, -0.04140931310657406, -0.41630675517520277, 0.11141481523298999  …  0.079825067546233, 0.011147081592849025, -0.005309780300181497, 0.04158697906021826, -0.037676577401611885, -0.23679753741417786, -0.017909345620873693, 0.05189485939983817, 0.017688703602043946, -0.00875819102976772], gp = [0.0, -0.011845868156976034, 1.720140350167542, 0.0, 0.00817141599764446, 1.6938804603926274, 0.0, 0.003144109471192846, 0.5120597583927129, 1.632119407923321  …  -0.04401939041434011, 0.7423009347864573, 2.4433389604712628, 0.0, -0.007501962465781716, 8.035437017197811, -2.422351205594254, -0.06196343157542439, 0.7082278243311373, 2.4207586442959976])), (sky = (impixel = [0.013441215266534676 0.015548354342624054 … 0.02707289832058558 0.024435013204324972; 0.01539014379228806 0.017777575602651613 … 0.029731759959913267 0.028589881846785164; … ; 0.027361134001015235 0.03643422970704207 … 0.014264753666071477 0.012822436718106869; 0.022568201120293403 0.033584230917270336 … 0.012192195764841084 3.0240703147351368e-5], fg = 0.009519817709150594), instrument = (lg = [0.040097957190886296, -0.04191010537120548, 0.22721999414497276, -0.022174532935162604, 0.022287908668420172, 0.30786256830109054, 0.029700790665919602, -0.03142785178247717, -0.2745971018096853, 0.09322850382923188  …  0.1094489427583667, 0.01186223178355275, -0.003967101492655268, 0.0422502458025887, -0.03864964693818655, -0.0029971448063534035, 0.010756327462598263, 0.08238603141016752, 0.0213408370409757, -0.006407302427419483], gp = [0.0, -0.012842261475362103, -1.6870798194698418, 0.0, 0.0069189889759192266, -1.7290458709809853, 0.0, 0.0021612287947506536, -0.43213664974020793, -1.8315197598750144  …  1.3846039805665704, -3.9511150346761874, -2.2503796674263072, 0.0, -0.007338519881097632, -4.963865659122968, 0.5130051553624224, 1.4031081744930898, -3.936079943667181, -2.223500473300341])), (sky = (impixel = [0.015280592835368602 0.016888818125155476 … 0.03130646171660247 0.033017023239183756; 0.016476775633160592 0.01854480135445185 … 0.026951315322809963 0.028166828154406236; … ; 0.020949864649752275 0.019160442830889746 … 0.01917678438286594 0.01780645495968482; 0.020075141495790882 0.019211232338352868 … 0.018032600184138537 0.000853801091354461], fg = 0.0138143122730759), instrument = (lg = [0.04243677172789193, -0.043854345333596986, 0.2937452088940615, -0.020995131151342352, 0.021441367960363573, 0.3696426308886042, 0.040817057482141234, -0.04110249741496388, -0.3886218975333081, 0.13196341666322825  …  0.046641599489931976, 0.01209030941896644, -0.003308257462947773, 0.05260021238174424, -0.04669953793478552, -0.06532886905375648, -0.07021448243535963, 0.0122548387062525, 0.02222089760741121, -0.005340893750370735], gp = [0.0, -0.01241046231188238, 2.7109806031010253, 0.0, 0.008060508895743597, 2.7183913809693414, 0.0, 0.0026745293824776504, -0.09310167299963205, 2.741227800044074  …  -3.087535375622896, -2.243142585778503, -0.5421365037860262, 0.0, -0.006998817835539079, -1.162115596416734, 1.9156520395319696, -3.064317533887305, -2.235831771824322, -0.5233557971464764])), (sky = (impixel = [0.051425696812565466 0.06475872803906368 … 0.018211615955449054 0.01577979367541639; 0.07208025507861712 0.08631707090475153 … 0.019885954384702982 0.017119242192147298; … ; 0.012778017916469861 0.013668189487316655 … 0.00381376588987237 0.003451530960739808; 0.012512988972101428 0.01360827308072558 … 0.003533341927671574 -0.000297645862500775], fg = 0.015803390701168377), instrument = (lg = [0.04177405821399194, -0.04330546970351625, 0.2712813326794514, -0.020590909395077534, 0.02095700628839115, 0.3537031376989496, 0.02783137707982636, -0.02967921446286997, -0.3948666430699576, 0.13884614171973197  …  0.09966989297601339, 0.011270927563175194, -0.005302286116112987, 0.052726342122735415, -0.04673563134316933, 0.06621099342679113, -0.2801398476845709, 0.06265838172488969, 0.021891934865992613, -0.006292747159763147], gp = [0.0, -0.01288713221565279, -0.2877696472906244, 0.0, 0.007249201942704609, -0.3020219704073036, 0.0, 0.0021971830605524607, 1.1577281058769415, -0.33909576049638535  …  0.2765735051441467, 0.20213443852951674, 1.9032294829529983, 0.0, -0.006471082650920775, -6.305809937410017, -0.2506013964024031, 0.3303626274378504, 0.24240617637911366, 1.9546735013226877])), (sky = (impixel = [0.009495653929475802 0.009553918377790186 … 0.018653688241786995 0.019581683720821162; 0.009611069539436199 0.009611059909862172 … 0.020869904904270355 0.02409703537634546; … ; 0.013517867586278116 0.013014326889236278 … 0.006589827140816798 0.006556213260172392; 0.012774050682595697 0.01234779095467627 … 0.0065476242867046845 0.15207484573159447], fg = 0.024147653034550437), instrument = (lg = [0.04130538328457871, -0.04292728760016511, 0.25879590364472504, -0.0211015551315059, 0.02136814128878209, 0.34018959979702634, 0.03675607574922615, -0.03755831629585652, -0.42468848447183916, 0.11470804780734008  …  0.03434638092598969, 0.010770707615055702, -0.006614295348504435, 0.04331699029544688, -0.03922960383902749, -0.19207341478484163, -0.17995621692107724, 0.004314486591638323, 0.01967256757777385, -0.009262151128626557], gp = [0.0, -0.012542730277176761, 0.23964021510856776, 0.0, 0.008610804456515119, 0.28496446769177836, 0.0, 0.0028466089575114525, 2.2898392411853163, 0.3994724524939677  …  2.027752503179682, -1.5431312734677836, 0.1581862625431607, 0.0, -0.006529459493124263, -2.8983327387650273, 0.9877998699043422, 2.062756623132175, -1.539673093710795, 0.17271151730068462]))], [2039.042069660225, 2032.0132279102982, 2030.3593215491917, 557.1642661110371, -208.26492040794975])

Now we plot the MAP estimate.

julia
import CairoMakie as CM
grid_plot = imagepixels(μas2rad(fov), μas2rad(fov), npix*4, npix*4)
img = intensitymap(Comrade.skymodel(post, xopts[1]), grid_plot)
fig = imageviz(img);

This page was generated using Literate.jl.