From 11daccc0efd74b6b75ff89632720d3d836555e02 Mon Sep 17 00:00:00 2001 From: Paul Tiede Date: Tue, 31 Oct 2023 17:09:18 -0400 Subject: [PATCH] Update examples --- examples/hybrid_imaging.jl | 6 +++--- examples/imaging_closures.jl | 4 ++-- examples/imaging_pol.jl | 19 ++++++++++--------- examples/imaging_vis.jl | 8 ++++---- 4 files changed, 19 insertions(+), 18 deletions(-) diff --git a/examples/hybrid_imaging.jl b/examples/hybrid_imaging.jl index 173e632b..954e2a5e 100644 --- a/examples/hybrid_imaging.jl +++ b/examples/hybrid_imaging.jl @@ -44,7 +44,7 @@ obs = ehtim.obsdata.load_uvfits(joinpath(dirname(pathof(Comrade)), "..", "exampl # Now we do some minor preprocessing: # - Scan average the data since the data have been preprocessed so that the gain phases # coherent. -obs = scan_average(obs).add_fractional_noise(0.01) +obs = scan_average(obs).add_fractional_noise(0.02) # For this tutorial we will once again fit complex visibilities since they # provide the most information once the telescope/instrument model are taken @@ -109,7 +109,7 @@ end # Now let's define our metadata. First we will define the cache for the image. This is # required to compute the numerical Fourier transform. -fovxy = μas2rad(150.0) +fovxy = μas2rad(250.0) npix = 32 grid = imagepixels(fovxy, fovxy, npix, npix) buffer = IntensityMap(zeros(npix,npix), grid) @@ -221,7 +221,7 @@ post = Posterior(lklhd, prior) xrand = prior_sample(rng, post) # and then plot the results import WGLMakie as CM -g = imagepixels(μas2rad(150.0), μas2rad(150.0), 128, 128) +g = imagepixels(μas2rad(250.0), μas2rad(250.0), 128, 128) imageviz(intensitymap(skymodel(post, xrand), g)) # ## Reconstructing the Image diff --git a/examples/imaging_closures.jl b/examples/imaging_closures.jl index 2b938b20..41d4cd4e 100644 --- a/examples/imaging_closures.jl +++ b/examples/imaging_closures.jl @@ -123,9 +123,9 @@ end # 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, Uniform(0.0, 1_000.0)) +cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat))) -prior = NamedDist(c = cprior, σimg = Uniform(0.0, 5.0), fg=Uniform(0.0, 1.0)) +prior = NamedDist(c = cprior, σimg = truncated(Normal(0.0, 0.1); lower = 0.0), fg=Uniform(0.0, 1.0)) lklhd = RadioLikelihood(sky, dlcamp, dcphase; skymeta = metadata) diff --git a/examples/imaging_pol.jl b/examples/imaging_pol.jl index 1768083d..55e42f81 100644 --- a/examples/imaging_pol.jl +++ b/examples/imaging_pol.jl @@ -102,7 +102,7 @@ Page(exportable=true, offline=true) # hide # For reproducibility we use a stable random number genreator using StableRNGs -rng = StableRNG(123) +rng = StableRNG(125) # Now we will load some synthetic polarized data. @@ -303,7 +303,7 @@ using OptimizationOptimJL using Zygote f = OptimizationFunction(tpost, Optimization.AutoZygote()) ℓ = logdensityof(tpost) -prob = Optimization.OptimizationProblem(f, prior_sample(tpost), nothing) +prob = Optimization.OptimizationProblem(f, prior_sample(rng, tpost), nothing) sol = solve(prob, LBFGS(), maxiters=15_000, g_tol=1e-1); # !!! warning @@ -329,17 +329,18 @@ img = intensitymap!(copy(imgtruesub), skymodel(post, xopt)) #Plotting the results gives import WGLMakie as CM -fig = CM.Figure(;resolution=(450, 200)); +fig = CM.Figure(;resolution=(450, 350)); polimage(fig[1,1], imgtruesub, axis=(xreversed=true, aspect=1, title="Truth", limits=((-20.0,20.0), (-20.0, 20.0))), - length_norm=1, plot_total=true, - pcolorrange=(-0.25, 0.25), pcolormap=CM.Reverse(:jet)) + length_norm=1, plot_total=true, pcolormap=:RdBu, + pcolorrange=(-0.25, 0.25),) polimage(fig[1,2], img, axis=(xreversed=true, aspect=1, title="Recon.", limits=((-20.0,20.0), (-20.0, 20.0))), - length_norm=1, plot_total=true, - pcolorrange=(-0.25, 0.25), pcolormap=CM.Reverse(:jet)) -CM.Colorbar(fig[1,3], colormap=CM.Reverse(:jet), colorrange=(-0.25, 0.25), label="Signed Polarization Fraction sign(V)*|p|") -CM.colgap!(fig.layout, 1) + length_norm=1, plot_total=true, pcolormap=:RdBu, + pcolorrange=(-0.25, 0.25),) +CM.Colorbar(fig[2,:], colormap=:RdBu, vertical=false, colorrange=(-0.25, 0.25), label="Signed Polarization Fraction sign(V)*|p|", flipaxis=false) +CM.colgap!(fig.layout, 3) +CM.rowgap!(fig.layout, 3) fig #- diff --git a/examples/imaging_vis.jl b/examples/imaging_vis.jl index 2fbaa998..05c09ea1 100644 --- a/examples/imaging_vis.jl +++ b/examples/imaging_vis.jl @@ -24,7 +24,7 @@ using LinearAlgebra # For reproducibility we use a stable random number genreator using StableRNGs -rng = StableRNG(42) +rng = StableRNG(123) @@ -39,7 +39,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 # 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)) +obs = scan_average(obs.add_fractional_noise(0.02)) # Now we extract our complex visibilities. dvis = extract_table(obs, ComplexVisibilities()) @@ -206,7 +206,7 @@ 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, Uniform(0.0, 1000.0))#InverseGamma(1.0, -log(0.01*rat))) +cprior = HierarchicalPrior(fmap, InverseGamma(1.0, -log(0.01*rat))) # We can now form our model parameter priors. Like our other imaging examples, we use a @@ -357,7 +357,7 @@ imgs = intensitymap.(samples, fovx, fovy, 128, 128) mimg = mean(imgs) simg = std(imgs) -fig = CM.Figure(;resolution=(800, 800)) +fig = CM.Figure(;resolution=(400, 400)) CM.image(fig[1,1], mimg, axis=(xreversed=true, aspect=1, title="Mean Image"), colormap=:afmhot)