update examples

examples/geometric_modeling.jl

@@ -93,10 +93,10 @@ prior = NamedDist(
 
 post = Posterior(lklhd, prior)
 
-# !!! notes
-# As of Comrade 0.9 we have switched to the proper covariant closure likelihood.
-# This is slower than the naieve diagonal liklelihood, but takes into account the
-# correlations between closures that share the same baselines.
+# !!!warn
+#    As of Comrade 0.9 we have switched to the proper covariant closure likelihood.
+#    This is slower than the naieve diagonal liklelihood, but takes into account the
+#    correlations between closures that share the same baselines.
 
 # This constructs a posterior density that can be evaluated by calling `logdensityof`.
 # For example,

examples/hybrid_imaging.jl

@@ -55,18 +55,19 @@ dvis  = extract_table(obs, ComplexVisibilities())
 # and a large asymmetric Gaussian component to model the unresolved short-baseline flux.
 
 function sky(θ, metadata)
-    (;c, σimg, f, r, σ, ma, mp, fg) = θ
+    (;c, σimg, f, r, σ, ma1, mp1, ma2, mp2, fg) = θ
     (;meanpr, grid, cache) = metadata
     ## Form the image model
     ## First transform to simplex space first applying the non-centered transform
     rast = to_simplex(CenteredLR(), meanpr .+ σimg.*c)
     img = IntensityMap((f*(1-fg))*rast, grid)
     mimg = ContinuousImage(img, cache)
     ## Form the ring model
-    s,c = sincos(mp)
-    α = ma*c
-    β = ma*s
-    ring = ((1-f)*(1-fg))*smoothed(stretched(MRing(α, β), r, r),σ)
+    s1,c1 = sincos(mp1)
+    s2,c2 = sincos(mp2)
+    α = (ma1*c1, ma2*c2)
+    β = (ma1*s1, ma2*s2)
+    ring = ((1-f)*(1-fg))*smoothed(stretched(MRing(α, β), r, r), σ)
     gauss = fg*stretched(Gaussian(), μas2rad(200.0), μas2rad(200.0))
     ## We group the geometric models together for improved efficiency. This will be
     ## automated in future versions.
@@ -183,7 +184,7 @@ cprior = GaussMarkovRandomField(zero(meanpr), 0.1*rat, 1.0)
 # we let the prior expand to 100% due to the known pointing issues LMT had in 2017.
 using Distributions
 using DistributionsAD
-distamp = station_tuple(dvis, Normal(0.0, 0.1); LM = Normal(1.0))
+distamp = station_tuple(dvis, Normal(0.0, 0.1); LM = Normal(0.0, 1.0))
 
 # For the phases, as mentioned above, we will use a segmented gain prior.
 # This means that rather than the parameters
@@ -204,12 +205,14 @@ using VLBIImagePriors
 prior = NamedDist(
           ## We use a strong smoothing prior since we want to limit the amount of high-frequency structure in the raster.
           c  = cprior,
-          σimg = truncated(Normal(0.0, 0.2); lower=1e-3),
+          σimg = truncated(Normal(0.0, 0.25); lower=1e-3),
           f  = Uniform(0.0, 1.0),
           r  = Uniform(μas2rad(10.0), μas2rad(30.0)),
-          σ  = Uniform(μas2rad(0.1), μas2rad(20.0)),
-          ma = Uniform(0.0, 0.5),
-          mp = Uniform(0.0, 2π),
+          σ  = Uniform(μas2rad(0.1), μas2rad(5.0)),
+          ma1 = Uniform(0.0, 0.5),
+          mp1 = Uniform(0.0, 2π),
+          ma2 = Uniform(0.0, 0.5),
+          mp2 = Uniform(0.0, 2π),
           fg = Uniform(0.0, 1.0),
           lgamp = CalPrior(distamp, gcache),
           gphase = CalPrior(distphase, gcachep),
@@ -246,7 +249,7 @@ using ComradeOptimization
 using OptimizationOptimJL
 using Zygote
 f = OptimizationFunction(tpost, Optimization.AutoZygote())
-prob = Optimization.OptimizationProblem(f, rand(ndim) .- 0.5, nothing)
+prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing)
 sol = solve(prob, LBFGS(); maxiters=5_000);
 
 
@@ -311,8 +314,8 @@ plot(p1,p2,p3,p4, layout=(2,2), size=(650, 650))
 using StatsPlots
 p1 = density(rad2μas(chain.r)*2, xlabel="Ring Diameter (μas)")
 p2 = density(rad2μas(chain.σ)*2*sqrt(2*log(2)), xlabel="Ring FWHM (μas)")
-p3 = density(-rad2deg.(chain.mp) .+ 360.0, xlabel = "Ring PA (deg) E of N")
-p4 = density(2*chain.ma, xlabel="Brightness asymmetry")
+p3 = density(-rad2deg.(chain.mp1) .+ 360.0, xlabel = "Ring PA (deg) E of N")
+p4 = density(2*chain.ma1, xlabel="Brightness asymmetry")
 p5 = density(1 .- chain.f, xlabel="Ring flux fraction")
 plot(p1, p2, p3, p4, p5, size=(900, 600), legend=nothing)
 

examples/imaging_closures.jl

@@ -188,7 +188,7 @@ chain, stats = sample(post, AHMC(;metric, autodiff=Val(:Zygote)), 700; nadapts=5
 # unable to assess uncertainty in their reconstructions.
 #
 # To explore our posterior let's first create images from a bunch of draws from the posterior
-msamples = skymodel.(Ref(post), chain[1001:10:end]);
+msamples = skymodel.(Ref(post), chain[501:2:end]);
 
 # The mean image is then given by
 using StatsBase

examples/imaging_pol.jl

@@ -302,7 +302,7 @@ using Zygote
 f = OptimizationFunction(tpost, Optimization.AutoZygote())
 ℓ = logdensityof(tpost)
 prob = Optimization.OptimizationProblem(f, prior_sample(tpost), nothing)
-sol = solve(prob, LBFGS(), maxiters=10_000, g_tol=1e-2, callback=(x,p)->(@info f(x,p); false));
+sol = solve(prob, LBFGS(), maxiters=15_000, g_tol=1e-2);
 
 # !!! warning
 #     Fitting polarized images is generally much harder than Stokes I imaging. This difficulty means that