Skip to content

Commit

Permalink
Switch to WGLMakie
Browse files Browse the repository at this point in the history
  • Loading branch information
ptiede committed Oct 30, 2023
1 parent 5fa3b6f commit 5b3b5cf
Show file tree
Hide file tree
Showing 6 changed files with 23 additions and 22 deletions.
25 changes: 13 additions & 12 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,18 +7,6 @@ Pkg.activate(script_dir)
parent_dir = dirname(script_dir)
Pkg.develop(PackageSpec(path=parent_dir))


using Documenter
using Pyehtim
using Zygote
using Comrade
using ComradeBase
using VLBISkyModels
using InteractiveUtils

using Literate
using Pkg

function dev_subpkg(subpkg)
subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
Pkg.develop(PackageSpec(path=subpkg_path))
Expand All @@ -32,6 +20,19 @@ dev_subpkg("ComradeNested")
dev_subpkg("ComradeDynesty")
dev_subpkg("ComradeAdaptMCMC")


using Documenter
using Pyehtim
using Zygote
using Comrade
using ComradeBase
using VLBISkyModels
using InteractiveUtils

using Literate
using Pkg


using ComradeAHMC
using ComradeOptimization
using ComradeNested
Expand Down
2 changes: 1 addition & 1 deletion examples/geometric_modeling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ xopt = transform(fpost, sol)

# Given this we can now plot the optimal image or the *maximum a posteriori* (MAP) image.

import CairoMakie as CM
import WGLMakie as CM
g = imagepixels(μas2rad(200.0), μas2rad(200.0), 256, 256)
fig, ax, plt = CM.image(g, model(xopt); axis=(xreversed=true, aspect=1, xlabel="RA (μas)", ylabel="Dec (μas)"), figure=(;resolution=(650,500),) ,colormap=:afmhot)

Expand Down
8 changes: 4 additions & 4 deletions examples/hybrid_imaging.jl
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@ post = Posterior(lklhd, prior)
# To sample from our prior we can do
xrand = prior_sample(rng, post)
# and then plot the results
import CairoMakie as CM
import WGLMakie as CM
g = imagepixels(μas2rad(150.0), μas2rad(150.0), 128, 128)
imageviz(intensitymap(skymodel(post, xrand), g))

Expand Down Expand Up @@ -256,7 +256,7 @@ residual(vlbimodel(post, xopt), dvis, ylabel="Correlated Flux Residual")
# 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.
CM.image(g, skymodel(post, xopt), axis=(aspect=1, xreversed=true, title="MAP"), colormap=:afmhot)
CM.image(g, skymodel(post, xopt), axis=(aspect=1, xreversed=true, title="MAP"), colormap=:afmhot, figure=(;resolution=(400, 400),))


# We will now move directly to sampling at this point.
Expand Down Expand Up @@ -295,7 +295,7 @@ rast_imgs = intensitymap.(rast_samples, fovxy, fovxy, 128, 128)
ring_mean, ring_std = mean_and_std(ring_imgs)
rast_mean, rast_std = mean_and_std(rast_imgs)

fig = CM.Figure(; resolution=(800, 800))
fig = CM.Figure(; resolution=(400, 400))
axes = [CM.Axis(fig[i, j], xreversed=true, aspect=CM.DataAspect()) for i in 1:2, j in 1:2]
CM.image!(axes[1,1], ring_mean, colormap=:afmhot); axes[1,1].title = "Ring Mean"
CM.image!(axes[1,2], ring_std, colormap=:afmhot); axes[1,2].title = "Ring Std. Dev."
Expand All @@ -304,7 +304,7 @@ CM.image!(axes[2,2], rast_std, colormap=:afmhot); axes[2,2].title = "Rast Std. D
fig

# Finally, let's take a look at some of the ring parameters
figd = CM.Figure(;resolution=(900, 600))
figd = CM.Figure(;resolution=(600, 400))
p1 = CM.density(figd[1,1], rad2μas(chain.r)*2, axis=(xlabel="Ring Diameter (μas)",))
p2 = CM.density(figd[1,2], rad2μas(chain.σ)*2*sqrt(2*log(2)), axis=(xlabel="Ring FWHM (μas)",))
p3 = CM.density(figd[1,3], -rad2deg.(chain.mp1) .+ 360.0, axis=(xlabel = "Ring PA (deg) E of N",))
Expand Down
6 changes: 3 additions & 3 deletions examples/imaging_closures.jl
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,9 @@ residual(skymodel(post, xopt), dlcamp, ylabel="Log Closure Amplitude Res.")
residual(skymodel(post, xopt), dcphase, ylabel="|Closure Phase Res.|")

# Now let's plot the MAP estimate.
import CairoMakie as CM
import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), μas2rad(150.0), μas2rad(150.0), 100, 100)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot, figure=(;resolution=(400, 400),))

# To sample from the posterior we will use HMC and more specifically the NUTS algorithm. For information about NUTS
# see Michael Betancourt's [notes](https://arxiv.org/abs/1701.02434).
Expand Down Expand Up @@ -194,7 +194,7 @@ using StatsBase
imgs = intensitymap.(msamples, μas2rad(150.0), μas2rad(150.0), 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)
Expand Down
2 changes: 1 addition & 1 deletion examples/imaging_pol.jl
Original file line number Diff line number Diff line change
Expand Up @@ -323,7 +323,7 @@ imgtrue = Comrade.load(joinpath(dirname(pathof(Comrade)), "..", "examples", "Pol
imgtruesub = imgtrue(Interval(-fovx/2, fovx/2), Interval(-fovy/2, fovy/2))
img = intensitymap!(copy(imgtruesub), skymodel(post, xopt))
import CairoMakie as CM
fig = CM.Figure(;resolution=(900, 400))
fig = CM.Figure(;resolution=(450, 200));
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,
Expand Down
2 changes: 1 addition & 1 deletion examples/imaging_vis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -269,7 +269,7 @@ residual(vlbimodel(post, xopt), dvis)
# improved in a few ways, but that is beyond the goal of this quick tutorial.
# Plotting the image, we see that we have a much cleaner version of the closure-only image from
# [Imaging a Black Hole using only Closure Quantities](@ref).
import CairoMakie as CM
import WGLMakie as CM
img = intensitymap(skymodel(post, xopt), fovx, fovy, 128, 128)
CM.image(img, axis=(xreversed=true, aspect=1, title="MAP Image"), colormap=:afmhot)

Expand Down

0 comments on commit 5b3b5cf

Please sign in to comment.