update examples

examples/hybrid_imaging.jl

@@ -55,20 +55,20 @@ 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, σ, ma1, mp1, ma2, mp2, fg) = θ
-    (;meanpr, grid, cache) = metadata
+    (;c, σimg, f, r, σ, τ, ξτ, ma1, mp1, ma2, mp2, fg) = θ
+    (;ftot, 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)
+    img = IntensityMap((ftot*f*(1-fg))*rast, grid)
     mimg = ContinuousImage(img, cache)
     ## Form the ring model
     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))
+    ring = smoothed(modify(MRing(α, β), Stretch(r, r*(1+τ)), Rotate(ξτ), Renormalize((ftot*(1-f)*(1-fg)))), σ)
+    gauss = modify(Gaussian(), Stretch(μas2rad(250.0)), Renormalize(ftot*f))
     ## We group the geometric models together for improved efficiency. This will be
     ## automated in future versions.
     return mimg + (ring + gauss)
@@ -123,9 +123,10 @@ cache  = create_cache(NFFTAlg(dvis), buffer, BSplinePulse{3}())
 # of priors and transformations that are useful for imaging.
 using VLBIImagePriors
 # Since we are using a Gaussian Markov random field prior we need to first specify our `mean`
-# image. For this work we will use a symmetric Gaussian with a FWHM of 40 μas
+# image. For this work we will use a symmetric Gaussian with a FWHM of 80 μas. This is larger
+# than the other examples since the ring will attempt to soak up the majority of the ring flux.
 fwhmfac = 2*sqrt(2*log(2))
-mpr = modify(Gaussian(), Stretch(μas2rad(60.0)./fwhmfac))
+mpr = modify(Gaussian(), Stretch(μas2rad(80.0)./fwhmfac))
 imgpr = intensitymap(mpr, grid)
 
 # Now since we are actually modeling our image on the simplex we need to ensure that
@@ -134,7 +135,7 @@ imgpr ./= flux(imgpr)
 meanpr = to_real(CenteredLR(), baseimage(imgpr))
 
 # Now we form the metadata
-skymetadata = (;meanpr, grid, cache)
+skymetadata = (;ftot=1.1, meanpr, grid, cache)
 
 
 # Second, we now construct our instrument model cache. This tells us how to map from the gains
@@ -167,8 +168,7 @@ lklhd = RadioLikelihood(sky, instrument, dvis;
 # To enforce this we will set the
 # length scale of the raster component equal to the beam size of the telescope in units of
 # pixel length, which is given by
-hh(x) = hypot(x...)
-beam = inv(maximum(hh.(uvpositions.(extract_table(obs, ComplexVisibilities()).data))))
+beam = beamsize(dvis)
 rat = (beam/(step(grid.X)))
 cprior = GaussMarkovRandomField(zero(meanpr), 0.05*rat, 1.0)
 # additionlly we will fix the standard deviation of the field to unity and instead
@@ -203,12 +203,14 @@ distphase = station_tuple(dvis, DiagonalVonMises(0.0, inv(π^2)))
 using VLBIImagePriors
 # Finally we can put form the total model prior
 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.5); lower=1e-3),
+          ## We use a strong smoothing prior since we want to limit the amount of high-frequency structure in the raster.
+          σimg = truncated(Normal(0.0, 1.0); lower=0.01),
           f  = Uniform(0.0, 1.0),
           r  = Uniform(μas2rad(10.0), μas2rad(30.0)),
-          σ  = Uniform(μas2rad(0.1), μas2rad(5.0)),
+          σ  = Uniform(μas2rad(0.1), μas2rad(10.0)),
+          τ  = truncated(Normal(0.0, 0.1); lower=0.0, upper=1.0),
+          ξτ = Uniform(-π/2, π/2),
           ma1 = Uniform(0.0, 0.5),
           mp1 = Uniform(0.0, 2π),
           ma2 = Uniform(0.0, 0.5),
@@ -291,8 +293,10 @@ msamples = skymodel.(Ref(post), chain[begin:2:end]);
 # The mean image is then given by
 imgs = intensitymap.(msamples, fovxy, fovxy, 128, 128)
 plot(mean(imgs), title="Mean Image")
+#-
 plot(std(imgs), title="Std Dev.")
-
+#-
+#
 # We can also split up the model into its components and analyze each separately
 comp = Comrade.components.(msamples)
 ring_samples = getindex.(comp, 2)

examples/imaging_closures.jl

@@ -43,7 +43,7 @@ obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "exampl
 #   - Scan average the data since the data have been preprocessed so that the gain phases
 #      are coherent.
 #   - Add 1% systematic noise to deal with calibration issues that cause 1% non-closing errors.
-obs = scan_average(obs).add_fractional_noise(0.01).flag_uvdist(uv_min=0.1e9)
+obs = scan_average(obs).add_fractional_noise(0.015)
 
 # Now, we extract our closure quantities from the EHT data set.
 dlcamp, dcphase  = extract_table(obs, LogClosureAmplitudes(;snrcut=3), ClosurePhases(;snrcut=3))
@@ -53,22 +53,26 @@ dlcamp, dcphase  = extract_table(obs, LogClosureAmplitudes(;snrcut=3), ClosurePh
 # convolved with some pulse or kernel to make a `ContinuousImage` object with it `Comrade's.`
 # generic image model. Note that `ContinuousImage(img, cache)` actually creates a [Comrade.modelimage](@ref)
 # object that allows `Comrade` to numerically compute the Fourier transform of the image.
-function model(θ, metadata)
-    (;c, σimg) = θ
-    (;meanpr, grid, cache) = metadata
-    ## Construct the image model
-    c = to_simplex(CenteredLR(), meanpr .+ σimg*c.params)
-    img = IntensityMap(c, grid)
-    return  ContinuousImage(img, cache)
+function sky(θ, metadata)
+    (;fg, c, σimg) = θ
+    (;K, meanpr, grid, cache) = metadata
+    ## Construct the image model we fix the flux to 0.6 Jy in this case
+    cp = meanpr .+ σimg.*c.params
+    rast = ((1-fg))*K(to_simplex(CenteredLR(), cp))
+    img = IntensityMap(rast, grid)
+    m = ContinuousImage(img, cache)
+    ## Add a large-scale gaussian to deal with the over-resolved mas flux
+    g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(fg))
+    return m + g
 end
 
 
 # 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 views; 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.
-npix = 24
-fovxy = μas2rad(100.0)
+npix = 32
+fovxy = μas2rad(150.0)
 
 # Now, we can feed in the array information to form the cache
 grid = imagepixels(fovxy, fovxy, npix, npix)
@@ -91,13 +95,12 @@ imgpr = intensitymap(mpr, grid)
 imgpr ./= flux(imgpr)
 
 meanpr = to_real(CenteredLR(), Comrade.baseimage(imgpr))
-metadata = (;meanpr, grid, cache)
+metadata = (;meanpr,K=CenterImage(imgpr), grid, cache)
 
 # In addition we want a reasonable guess for what the resolution of our image should be.
 # For radio astronomy this is given by roughly the longest baseline in the image. To put this
 # into pixel space we then divide by the pixel size.
-hh(x) = hypot(x...)
-beam = inv(maximum(hh.(uvpositions.(extract_table(obs, ComplexVisibilities()).data))))
+beam = beamsize(dlcamp)
 rat = (beam/(step(grid.X)))
 
 # To make the Gaussian Markov random field efficient we first precompute a bunch of quantities
@@ -120,11 +123,11 @@ end
 # to prefer "simple" structures. Namely, Gaussian Markov random fields are extremly flexible models.
 # To prevent overfitting it is common to use priors that penalize complexity. Therefore, we
 # want to use priors that enforce similarity to our mean image, and prefer smoothness.
-cprior = HierarchicalPrior(fmap, NamedDist(λ = truncated(Normal(0.0, 0.25*inv(rat)); lower=2/npix)))
+cprior = HierarchicalPrior(fmap, NamedDist(λ = truncated(Normal(0.0, 0.1*inv(rat)); lower=2/npix)))
 
-prior = NamedDist(c = cprior, σimg = truncated(Normal(0.0, 0.5); lower=0.01))
+prior = NamedDist(c = cprior, σimg = truncated(Normal(0.0, 1.0); lower=0.01), fg=Uniform(0.0, 1.0))
 
-lklhd = RadioLikelihood(model, dlcamp, dcphase;
+lklhd = RadioLikelihood(sky, dlcamp, dcphase;
                         skymeta = metadata)
 post = Posterior(lklhd, prior)
 
@@ -166,7 +169,7 @@ residual(skymodel(post, xopt), dcphase, ylabel="|Closure Phase Res.|")
 # Now these residuals look a bit high. However, it turns out this is because the MAP is typically
 # not a great estimator and will not provide very predictive measurements of the data. We
 # will show this below after sampling from the posterior.
-img = intensitymap(skymodel(post, xopt), μas2rad(100.0), μas2rad(100.0), 100, 100)
+img = intensitymap(skymodel(post, xopt), μas2rad(150.0), μas2rad(150.0), 100, 100)
 plot(img, title="MAP Image")
 
 # To sample from the posterior we will use HMC and more specifically the NUTS algorithm. For information about NUTS
@@ -192,11 +195,11 @@ msamples = skymodel.(Ref(post), chain[501:2:end]);
 
 # The mean image is then given by
 using StatsBase
-imgs = center_image.(intensitymap.(msamples, μas2rad(100.0), μas2rad(100.0), 128, 128))
+imgs = intensitymap.(msamples, μas2rad(150.0), μas2rad(150.0), 128, 128)
 mimg = mean(imgs)
 simg = std(imgs)
 p1 = plot(mimg, title="Mean Image")
-p2 = plot(simg./(mimg), title="1/SNR", clims=(0.0, 2.0))
+p2 = plot(simg./(max.(mimg, 1e-5)), title="1/SNR", clims=(0.0, 2.0))
 p3 = plot(imgs[1], title="Draw 1")
 p4 = plot(imgs[end], title="Draw 2")
 plot(p1, p2, p3, p4, size=(800,800), colorbar=:none)

examples/imaging_vis.jl

@@ -51,11 +51,13 @@ dvis = extract_table(obs, ComplexVisibilities())
 function sky(θ, metadata)
     (;fg, c, σimg) = θ
     (;ftot, K, meanpr, grid, cache) = metadata
-    ## Construct the image model we fix the flux to 0.6 Jy in this case
+    ## Transform to the log-ratio pixel fluxes
     cp = meanpr .+ σimg.*c.params
+    ## Transform to image space
     rast = (ftot*(1-fg))*K(to_simplex(CenteredLR(), cp))
     img = IntensityMap(rast, grid)
     m = ContinuousImage(img, cache)
+    ## Add a large-scale gaussian to deal with the over-resolved mas flux
     g = modify(Gaussian(), Stretch(μas2rad(250.0), μas2rad(250.0)), Renormalize(ftot*fg))
     return m + g
 end