update

examples/hybrid_imaging.jl

@@ -125,7 +125,7 @@ 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
 fwhmfac = 2*sqrt(2*log(2))
-mpr = modify(Gaussian(), Stretch(μas2rad(100.0)./fwhmfac))
+mpr = modify(Gaussian(), Stretch(μas2rad(60.0)./fwhmfac))
 imgpr = intensitymap(mpr, grid)
 
 # Now since we are actually modeling our image on the simplex we need to ensure that
@@ -170,7 +170,7 @@ lklhd = RadioLikelihood(sky, instrument, dvis;
 hh(x) = hypot(x...)
 beam = inv(maximum(hh.(uvpositions.(extract_table(obs, ComplexVisibilities()).data))))
 rat = (beam/(step(grid.X)))
-cprior = GaussMarkovRandomField(zero(meanpr), 0.1*rat, 1.0)
+cprior = GaussMarkovRandomField(zero(meanpr), 0.05*rat, 1.0)
 # additionlly we will fix the standard deviation of the field to unity and instead
 # use a pseudo non-centered parameterization for the field.
 # GaussMarkovRandomField(meanpr, 0.1*rat, 1.0, crcache)
@@ -205,7 +205,7 @@ 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.25); lower=1e-3),
+          σimg = truncated(Normal(0.0, 0.5); lower=1e-3),
           f  = Uniform(0.0, 1.0),
           r  = Uniform(μas2rad(10.0), μas2rad(30.0)),
           σ  = Uniform(μas2rad(0.1), μas2rad(5.0)),

examples/imaging_closures.jl

@@ -67,7 +67,7 @@ end
 # 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 = 32
+npix = 24
 fovxy = μas2rad(100.0)
 
 # Now, we can feed in the array information to form the cache
@@ -83,7 +83,7 @@ using VLBIImagePriors, Distributions, DistributionsAD
 # 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
 fwhmfac = 2*sqrt(2*log(2))
-mpr = modify(Gaussian(), Stretch(μas2rad(70.0)./fwhmfac))
+mpr = modify(Gaussian(), Stretch(μas2rad(50.0)./fwhmfac))
 imgpr = intensitymap(mpr, grid)
 
 # Now since we are actually modeling our image on the simplex we need to ensure that
@@ -122,7 +122,7 @@ end
 # 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)))
 
-prior = NamedDist(c = cprior, σimg = truncated(Normal(0.0, 1.0); lower=0.01))
+prior = NamedDist(c = cprior, σimg = truncated(Normal(0.0, 0.5); lower=0.01))
 
 lklhd = RadioLikelihood(model, dlcamp, dcphase;
                         skymeta = metadata)
@@ -196,14 +196,14 @@ imgs = center_image.(intensitymap.(msamples, μas2rad(100.0), μas2rad(100.0), 1
 mimg = mean(imgs)
 simg = std(imgs)
 p1 = plot(mimg, title="Mean Image")
-p2 = plot(simg./(mimg), title="1/SNR", clims=(0.0, 1.0))
+p2 = plot(simg./(mimg), 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)
 
 # Now let's see whether our residuals look better.
 p = plot();
-for s in sample(chain[1001:end], 10)
+for s in sample(chain[501:end], 10)
     residual!(p, vlbimodel(post, s), dlcamp)
 end
 ylabel!("Log-Closure Amplitude Res.");
@@ -212,10 +212,10 @@ p
 
 
 p = plot();
-for s in sample(chain[1001:end], 10)
+for s in sample(chain[501:end], 10)
     residual!(p, vlbimodel(post, s), dcphase)
 end
-ylabel!("Closure Phase Res.");
+ylabel!("|Closure Phase Res.|");
 p
 
 # And we see that the residuals are looking much better.
@@ -230,20 +230,3 @@ p
 #       you want to use a sampler that can deal with multi-modal posteriors. Check out the package [`Pigeons.jl`](https://github.com/Julia-Tempering/Pigeons.jl)
 #       for an **in-development** package that should easily enable this type of sampling.
 #-
-
-
-# ## Computing information
-# ```
-# Julia Version 1.8.5
-# Commit 17cfb8e65ea (2023-01-08 06:45 UTC)
-# Platform Info:
-#   OS: Linux (x86_64-linux-gnu)
-#   CPU: 32 × AMD Ryzen 9 7950X 16-Core Processor
-#   WORD_SIZE: 64
-#   LIBM: libopenlibm
-#   LLVM: libLLVM-13.0.1 (ORCJIT, znver3)
-#   Threads: 1 on 32 virtual cores
-# Environment:
-#   JULIA_EDITOR = code
-#   JULIA_NUM_THREADS = 1
-# ```

examples/imaging_pol.jl

@@ -160,8 +160,8 @@ function instrument(θ, metadata)
     jT = jonesT(tcache)
 
     ## Gain product parameters
-    gPa = exp.(lgp .+ 0im)
-    gRa = exp.(lgp .+ lgr .+ 0im)
+    gPa = exp.(lgp)
+    gRa = exp.(lgp .+ lgr)
     Gp = jonesG(gPa, gRa, scancache)
     ## Gain ratio
     gPp = exp.(1im.*(gpp))
@@ -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=15_000, g_tol=1e-2);
+sol = solve(prob, LBFGS(), maxiters=15_000, g_tol=1e-1);
 
 # !!! warning
 #     Fitting polarized images is generally much harder than Stokes I imaging. This difficulty means that

examples/imaging_vis.jl

@@ -201,15 +201,15 @@ end
 # and 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. If the data wants more complexity
 # then it will drive us away from the prior.
-cprior = HierarchicalPrior(fmap, NamedDist((;λ = truncated(Normal(0.0, 0.1/rat); lower=2/npix))))
+cprior = HierarchicalPrior(fmap, NamedDist((;λ = truncated(Normal(0.0, 0.25*inv(rat)); lower=2/npix))))
 
 
 # We can now form our model parameter priors. Like our other imaging examples, we use a
 # Dirichlet prior for our image pixels. For the log gain amplitudes, we use the `CalPrior`
 # which automatically constructs the prior for the given jones cache `gcache`.
 prior = NamedDist(
          fg = Uniform(0.0, 1.0),
-         σimg = truncated(Normal(0.0, 0.1); lower=0.01),
+         σimg = truncated(Normal(0.0, 0.5); lower=0.01),
          c = cprior,
          lgamp = CalPrior(distamp, gcache),
          gphase = CalPrior(distphase, gcachep)
@@ -251,7 +251,7 @@ using OptimizationOptimJL
 f = OptimizationFunction(tpost, Optimization.AutoZygote())
 prob = Optimization.OptimizationProblem(f, rand(ndim) .- 0.5, nothing)
 ℓ = logdensityof(tpost)
-sol = solve(prob, LBFGS(), maxiters=1_000, g_tol=1e-1);
+sol = solve(prob, LBFGS(), maxiters=1000, g_tol=1e-1);
 
 # Now transform back to parameter space
 xopt = transform(tpost, sol.u)
@@ -347,7 +347,7 @@ plot(ctable_ph, layout=(3,3), size=(600,500))
 plot(ctable_am, layout=(3,3), size=(600,500))
 
 # Finally let's construct some representative image reconstructions.
-samples = skymodel.(Ref(post), chain[begin:10:end])
+samples = skymodel.(Ref(post), chain[begin:2:end])
 imgs = center_image.(intensitymap.(samples, fovx, fovy, 128,  128))
 
 mimg = mean(imgs)