diff --git a/Project.toml b/Project.toml
index 4d55fbd44..f83230fb3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.11.26"
 [deps]
 AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001"
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
 AstroLib = "c7932e45-9af1-51e7-9da9-f004cd3a462b"
 AstroTime = "c61b5328-d09d-5e37-a9a8-0eb41c39009c"
@@ -30,6 +31,7 @@ PolarizedTypes = "d3c5d4cd-a8ee-40d6-aac7-e34df5a20044"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReactantCore = "a3311ec8-5e00-46d5-b541-4f83e724a433"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -53,6 +55,7 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 Pigeons = "0eb8d820-af6a-4919-95ae-11206f830c31"
 Pyehtim = "3d61700d-6e5b-419a-8e22-9c066cf00468"
+Reactant = "3c362404-f566-11ee-1572-e11a4b42c853"
 VLBIFiles = "c1ebf4c8-f9d4-409a-8daf-7009448f4e6e"
 
 [extensions]
@@ -63,9 +66,11 @@ ComradeMakieExt = "Makie"
 ComradeOptimizationExt = "Optimization"
 ComradePigeonsExt = "Pigeons"
 ComradePyehtimExt = "Pyehtim"
+ComradeReactantExt = "Reactant"
 ComradeVLBIFilesExt = "VLBIFiles"
 
 [compat]
+Adapt = "4"
 AbstractMCMC = "3, 4, 5"
 Accessors = "0.1"
 AdvancedHMC = "0.6, 0.7, 0.8"
@@ -76,7 +81,7 @@ ChainRulesCore = "1"
 ComradeBase = "0.9.13"
 DelimitedFiles = "1"
 DensityInterface = "0.4"
-DimensionalData = "0.29 - 0.29.24, ^0.29.26, 0.30"
+DimensionalData = "0.30"
 Distributions = "0.25"
 DocStringExtensions = "0.8, 0.9"
 Dynesty = "0.4"
@@ -92,9 +97,11 @@ Optimization = "5"
 PaddedViews = "0.5"
 ParameterHandling = "0.4, 0.5"
 Pigeons = "0.4, 0.5"
-PolarizedTypes = "0.1.0 - 0.1.3"
+PolarizedTypes = "0.1"
 PrettyTables = "3"
 Pyehtim = "0.2"
+Reactant = "0.2"
+ReactantCore = "0.1"
 RecipesBase = "1"
 Reexport = "1"
 SpecialFunctions = "0.10, 1, 2"
@@ -105,8 +112,8 @@ StructArrays = "0.5,0.6, 0.7"
 Tables = "1"
 TransformVariables = "^0.8.19"
 VLBIFiles = "0.4.16"
-VLBIImagePriors = "0.10.8"
-VLBILikelihoods = "^0.2.6"
+VLBIImagePriors = "^0.10.9"
+VLBILikelihoods = "0.3.0"
 VLBISkyModels = "^0.6.24"
 julia = "1.10"
 
@@ -119,8 +126,9 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 Pigeons = "0eb8d820-af6a-4919-95ae-11206f830c31"
 Pyehtim = "3d61700d-6e5b-419a-8e22-9c066cf00468"
+Reactant = "3c362404-f566-11ee-1572-e11a4b42c853"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 VLBIFiles = "c1ebf4c8-f9d4-409a-8daf-7009448f4e6e"
 
 [targets]
-test = ["Test", "AdvancedHMC", "Dynesty", "Enzyme", "Makie", "Optimization", "Pigeons", "Pyehtim", "VLBIFiles"]
+test = ["Test", "AdvancedHMC", "Dynesty", "Enzyme", "Makie", "Optimization", "Pigeons", "Pyehtim", "VLBIFiles", "Reactant"]
diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md
index 6492d912a..ef2be3f69 100644
--- a/docs/src/benchmarks.md
+++ b/docs/src/benchmarks.md
@@ -60,15 +60,15 @@ function model(θ, p)
     return ring + g
 end
 prior = (
-          rad = Uniform(10.0, 30.0),
-          wid = Uniform(1.0, 10.0),
-          a = Uniform(-0.5, 0.5), b = Uniform(-0.5, 0.5),
-          f = Uniform(0.0, 1.0),
-          sig = Uniform((1.0), (60.0)),
-          asy = Uniform(0.0, 0.9),
-          pa = Uniform(0.0, 1π),
-          x = Uniform(-(80.0), (80.0)),
-          y = Uniform(-(80.0), (80.0))
+          rad = VLBIUniform(10.0, 30.0),
+          wid = VLBIUniform(1.0, 10.0),
+          a = VLBIUniform(-0.5, 0.5), b = VLBIUniform(-0.5, 0.5),
+          f = VLBIUniform(0.0, 1.0),
+          sig = VLBIUniform((1.0), (60.0)),
+          asy = VLBIUniform(0.0, 0.9),
+          pa = VLBIUniform(0.0, 1π),
+          x = VLBIUniform(-(80.0), (80.0)),
+          y = VLBIUniform(-(80.0), (80.0))
         )
 # Now form the posterior
 skym = SkyModel(model, prior, imagepixels(μas2rad(150.0), μas2rad(150.0), 128, 128))
diff --git a/docs/src/tutorials/index.md b/docs/src/tutorials/index.md
index 6e8ea7919..794b70bbb 100644
--- a/docs/src/tutorials/index.md
+++ b/docs/src/tutorials/index.md
@@ -6,8 +6,8 @@ const beginner = [
   {
     href: "beginner/LoadingData", 
     src: "../vis.png",
-    caption: "Loading Data with Pyehtim",
-    desc: "How to load data using standard eht-imaging in Julia." 
+    caption: "Loading Data with VLBIFiles",
+    desc: "How to load data into the Comrade format" 
   },
   {
     href: "beginner/GeometricModeling", 
diff --git a/examples/advanced/FitPS/main.jl b/examples/advanced/FitPS/main.jl
index 2c649aadd..48f57ee78 100644
--- a/examples/advanced/FitPS/main.jl
+++ b/examples/advanced/FitPS/main.jl
@@ -77,9 +77,9 @@ using Distributions
 # arguments and are in scope inside both the prior expressions and the body.
 @sky function sky(grid; mimg, pl, cprior, ρmax)
     c ~ cprior
-    ρs ~ ntuple(Returns(Uniform(0.01, ρmax)), 3)
-    σimg ~ Exponential(2.0)
-    fb ~ Uniform(0.0, 1.0)
+    ρs ~ ntuple(Returns(VLBIUniform(0.01, ρmax)), 3)
+    σimg ~ VLBIExponential(2.0)
+    fb ~ VLBIUniform(0.0, 1.0)
     ## Apply the GMRF fluctuations to the image
     x = genfield(StationaryRandomField(MarkovPS(ρs), pl), c)
     x .= σimg .* x
@@ -123,7 +123,7 @@ skym = sky(grid; mimg, pl, cprior, ρmax = max(size(grid)...))
 # Since we are fitting closures we do not need to include an instrument model, since
 # the closure likelihood is approximately independent of gains in the high SNR limit.
 using Enzyme
-post = VLBIPosterior(skym, dlcamp, dcphase)
+post = VLBIPosterior(skym, dlcamp, dcphase);
 
 # ## Reconstructing the Image
 
diff --git a/examples/advanced/Hibi/main.jl b/examples/advanced/Hibi/main.jl
index 432520988..44543fb8b 100644
--- a/examples/advanced/Hibi/main.jl
+++ b/examples/advanced/Hibi/main.jl
@@ -77,10 +77,10 @@ using VLBIImagePriors
 using Distributions
 @sky function sky(grid; ftot, cprior)
     c ~ cprior
-    σimg ~ truncated(Normal(0.0, 0.5); lower = 0.0)
-    r ~ Uniform(μas2rad(10.0), μas2rad(40.0))
-    ain ~ Uniform(1.0, 20.0)
-    aout ~ Uniform(1.0, 20.0)
+    σimg ~ VLBITruncated(VLBIGaussian(0.0, 0.5); lower = 0.0)
+    r ~ VLBIUniform(μas2rad(10.0), μas2rad(40.0))
+    ain ~ VLBIUniform(1.0, 20.0)
+    aout ~ VLBIUniform(1.0, 20.0)
     ## Form the image model
     mb = RingTemplate(RadialDblPower(ain, aout), AzimuthalUniform())
     mr = modify(mb, Stretch(r))
@@ -97,11 +97,13 @@ end
 #   - Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
 #   - Gain phases which are more difficult to constrain and can shift rapidly.
 
-fgain(x) = exp(x.lg + 1im * x.gp)
+using VLBIImagePriors
+using Distributions
+fgain(x) = exp(complex(x.lg, x.gp))
 G = SingleStokesGain(fgain)
 
 intpr = (
-    lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.2)); LM = IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
+    lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.2)); LM = IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0))),
     gp = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))); refant = SEFDReference(0.0)),
 )
 intmodel = InstrumentModel(G, intpr)
@@ -123,7 +125,7 @@ skym = sky(g; ftot = 1.1, cprior)
 # This is everything we need to specify our posterior distribution, which our is the main
 # object of interest in image reconstructions when using Bayesian inference.
 using Enzyme
-post = VLBIPosterior(skym, intmodel, dvis)
+post = VLBIPosterior(skym, intmodel, dvis);
 
 # We can sample from the prior to see what the model looks like
 using DisplayAs #hide
diff --git a/examples/advanced/NeuralFields/Project.toml b/examples/advanced/NeuralFields/Project.toml
new file mode 100644
index 000000000..f4c1afb13
--- /dev/null
+++ b/examples/advanced/NeuralFields/Project.toml
@@ -0,0 +1,32 @@
+[deps]
+AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Comrade = "99d987ce-9a1e-4df8-bc0b-1ea019aa547b"
+DisplayAs = "0b91fe84-8a4c-11e9-3e1d-67c38462b6d6"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
+Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
+NonuniformFFTs = "cd96f58b-6017-4a02-bb9e-f4d81626177f"
+Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Pyehtim = "3d61700d-6e5b-419a-8e22-9c066cf00468"
+Reactant = "3c362404-f566-11ee-1572-e11a4b42c853"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+VLBIFiles = "c1ebf4c8-f9d4-409a-8daf-7009448f4e6e"
+VLBIImagePriors = "b1ba175b-8447-452c-b961-7db2d6f7a029"
+VLBILikelihoods = "90db92cd-0007-4c0a-8e51-dbf0782ce592"
+VLBISkyModels = "d6343c73-7174-4e0f-bb64-562643efbeca"
+
+[compat]
+CairoMakie = "0.15"
+DisplayAs = "0.1"
+Distributions = "0.25"
+Optimisers = "0.4"
+Pkg = "1"
+Pyehtim = "0.2"
+Reactant = "0.2"
+StableRNGs = "1"
+VLBIImagePriors = "0.10"
diff --git a/examples/advanced/NeuralFields/main.jl b/examples/advanced/NeuralFields/main.jl
new file mode 100644
index 000000000..9cf4e6b9b
--- /dev/null
+++ b/examples/advanced/NeuralFields/main.jl
@@ -0,0 +1,328 @@
+import Pkg #hide
+__DIR = @__DIR__ #hide
+pkg_io = open(joinpath(__DIR, "pkg.log"), "w") #hide
+Pkg.activate(__DIR; io = pkg_io) #hide
+Pkg.develop(; path = joinpath(__DIR, "..", "..", ".."), io = pkg_io) #hide
+Pkg.instantiate(; io = pkg_io) #hide
+Pkg.precompile(; io = pkg_io) #hide
+close(pkg_io) #hide
+
+# # Neural Field Imaging with Lux in Comrade
+
+# In this tutorial, we will do closure-only modeling of the AGN DA 193 observed with the VLBA
+# at 15 GHz with the Mojave AGN project. Unlike other tutorials, we will not attempt a Bayesian
+# analysis and instead utilize a **neural field** to obtain a MAP estimate of the image.
+# This is mainly due to the challenges of characterizing posterior uncertainty with neural fields.
+# The goal here is not to provide a optimal analysis of this source, or neural fields, but rather
+# to demonstrate how Comrade's flexible modeling framework can be used to implement neural fields
+# with only minor modifications to existing code.
+
+
+# To get started, we will load Comrade
+using Comrade
+using LinearAlgebra
+using VLBIFiles
+using Lux
+
+# We will accelerate the optimization by tracing the forward+gradient through MLIR/XLA
+# with [Reactant.jl](https://github.com/EnzymeAD/Reactant.jl). Reactant precompiles a single
+# loss-and-gradient step into an XLA program, which is then re-used for every iteration of
+# our optimizer.
+using Reactant
+using VLBISkyModels: ReactantNUFFTAlg
+using Comrade: ReactantEx
+
+# For reproducibility we use a stable random number genreator
+using StableRNGs
+rng = StableRNG(123)
+
+
+# ## Load the Data
+# For this tutorial we will image publicly available VLBA data of the AGN
+# 1308+326 observed on 2021/03/19 at 43 GHz as part of the Boston University blazar monitoring program.
+file = Base.download("https://www.bu.edu/blazars/VLBA_GLAST/1308/1308+326Q.2021-03-19.UVP.gz")
+uvd = VLBIFiles.load(
+    VLBIFiles.UVData,
+    file
+)
+
+
+# For this tutorial we will only use closure quantities to reconstruct the image however, polarized
+# or complex visibilities can also be used with instrumental models following the other tutorials.
+dvis = extract_table(
+    uvd,
+    Visibilities(; time_average = VLBI.GapBasedScans()),
+)
+add_fractional_noise!(dvis, 0.005)
+
+
+# ## Build the Model/Posterior
+# For our model we will use a neural field to model the log image fluctuations. That is our model
+# will look like
+#
+#   δ(x) = NN(γ(x); θ)
+#
+# where NN is a neural network with weights and biases θ and $x$ are the spatial coordinates. γ denotes
+# a Fourier feature embedding of the spatial coordinates, which we include to potentially help with
+# learning high frequency structure in the image. This kind of model for VLBI is described in more detail in
+# the excellent paper [kine](@citet).
+#
+# Let's first define the Fourier feature embedding function.
+function fourierfeature(grid::RectiGrid; m = 4)
+    nx, ny = size(grid)
+    freqs = ntuple(n -> 2^(n - 1), m)
+    feats = zeros(Float32, 4 * m, nx, ny)
+    x = range(-1.0f0 * π, 1.0f0 * π, length = nx + 1)[begin:(end - 1)] #clip the endpoint
+    y = range(-1.0f0 * π, 1.0f0 * π, length = ny + 1)[begin:(end - 1)] #clip the endpoint
+    for i in 1:nx, j in 1:ny, (ik, k) in enumerate(1:4:(4 * m))
+        s, c = sincos(freqs[ik] * x[j])
+        feats[k, j, i] = s
+        feats[k + 1, j, i] = c
+        s, c = sincos(freqs[ik] * y[i])
+        feats[k + 2, j, i] = s
+        feats[k + 3, j, i] = c
+    end
+    return reshape(feats, 4 * m, nx * ny)
+end
+
+# For this tutorial we decided to image a very compact AGN. Thus, we will use a small FOV for a 15 GHz
+# observation. Namely, we will use a 1000 μas FOV with 64x64 pixels.
+nx = 128
+ny = 128
+fovx = μas2rad(1_000)
+fovy = fovx * ny / nx
+grid = imagepixels(fovx, fovy, nx, ny, μas2rad(150.0), -μas2rad(150.0))
+
+
+# Next we define neural network using the Lux.jl package.
+using Lux
+
+# We only use 1 Fourier feature embedding here for simplicity and because VLBI practioners typically
+# prefer very smooth images. However, more Fourier features can be used to model
+# finer scale structures.
+ff = fourierfeature(grid; m = 2)
+
+
+# For the neural field we use a very simple MLP with 3 hidden layers of 64 units each
+# and a final output layer with a single unit. We use the `tanh_fast` activation for speed
+# and to ensure smoothness in the image. Also be disable the bias in the last layer since
+# in some sense is prefers zero mean output which is desirable for fluctuation modeling.
+nnmodel = Chain(
+    Dense(size(ff, 1) => 64, Lux.tanh_fast),
+    Dense(64 => 64, Lux.tanh_fast),
+    Dense(64 => 64, Lux.tanh_fast),
+    Dense(64 => 1, use_bias = false),
+)
+
+
+# We can now define our sky model function. This function takes in the neural network parameters `nn`
+# and hyperparameters `fb, σ` and returns a `ContinuousImage` representing the sky brightness distribution.
+# This should look similar to essentially every other imaging model in the Comrade tutorials.
+function sky(θ, metadata)
+    (; nn, σ, fb, ftot) = θ
+    (; nnmodel, mimg, ff, st) = metadata
+    fbn = fb / length(mimg)
+    bimg = baseimage(mimg)
+    mb = bimg .* (1 - fb) .+ fbn
+
+    ## Compute the neural field
+    δ = reshape(first(Lux.apply(nnmodel, ff, nn, st)), size(mimg))
+    δ .*= σ
+    mδ = maximum(δ)
+    rast = @. exp(δ - mδ) * mb
+    rast .*= ftot / sum(rast)
+    return ContinuousImage(rast, axisdims(mimg), DeltaPulse{Float64}())
+end
+
+
+# For the prior on the network parameters we use the so-called *NNGP parameterization* of
+# [Neal (1996)](@cite) and [Lee et al. (2018)](@cite). For a `Dense` layer with weight
+# matrix `W ∈ ℝ^{out × in}` and bias `b ∈ ℝ^{out}` we use
+#
+#   W_{ij} ∼ 𝒩(0, σ_w² / fan_in),   b_i ∼ 𝒩(0, σ_b²)
+#
+# i.e. the weight variance is scaled by the inverse fan-in. With this scaling the pre-activation
+# variance is independent of the layer width, and as the hidden widths are taken to infinity the
+# network output converges to a well-defined Gaussian process (the NNGP). With unit-variance
+# `StdNormal` priors on every parameter the pre-activation variance grows linearly with width and
+# no proper GP limit exists, so this scaling is the principled choice for a Bayesian neural field.
+# We use `Lux.setup` to get the parameter layout of the network and then use `Comrade.rmap` to
+# walk the parameter tree, dispatching on parameter rank: rank-2 arrays are weight matrices (with
+# fan-in equal to `size(W, 2)`) and rank-1 arrays are biases.
+using Random
+ps, st = Lux.setup(Random.default_rng(), nnmodel) |> reactant_device()
+function nngp_prior(x::AbstractArray)
+    if ndims(x) == 2
+        fan_in = size(x, 2)
+        return VLBIGaussian(zero(eltype(x)), one(eltype(x)), size(x))
+    else
+        return VLBIImagePriors.StdNormal(size(x))
+    end
+end
+nnprior = Comrade.rmap(nngp_prior, ps)
+
+# Similar to the RF models we will use a so-called `mean image`. Note that in this case is isn't
+# actually clean what the mean image is, but this will help with centroid drift so we include it.
+# We will use a symmetric Gaussian with a FWHM equal to the approximate
+# beamsize of the array. This models the fact that we expect the AGN core to be compact.
+fwhmfac = 2 * sqrt(2 * log(2))
+mpr = modify(TBlob(4.0), Stretch(beamsize(dvis) / 3 / fwhmfac))
+imgpr = intensitymap(mpr, grid)
+skymeta = (; mimg = imgpr ./ sum(imgpr), nnmodel, ff, st);
+
+using Distributions
+# The only other hyperparameters we need to set priors for are the marginal standard deviation
+# and the high-frequency tapering parameter. For these we will use an Exponential prior with shape 1.0
+# for the standard deviation and a Uniform(0, 1) prior for constant background flux ratio.
+prior = (;
+    nn = nnprior,
+    fb = VLBIUniform(0.0, 1.0),
+    σ = VLBIExponential(0.5),
+    ftot = VLBIUniform(0.1, 10.0),
+)
+
+# We can then define our sky model. We pass `algorithm = ReactantAlg()` so the non-uniform
+# FFT used to evaluate visibilities runs through Reactant's XLA-traced NUFT plan rather
+# than the default CPU `NFFTAlg()`.
+skym = SkyModel(sky, prior, grid; algorithm = VLBISkyModels.ReactantNUFFTAlg(), metadata = skymeta)
+
+
+# Now we will fit gains since our data product are complex visibilities
+g(x) = exp(complex(x.lg, x.gp))
+G = SingleStokesGain(g)
+intpr = (
+    lg = ArrayPrior(
+        IIDSitePrior(IntegSeg(), VLBIGaussian(0.0, 0.5))
+    ),
+    gp = ArrayPrior(
+        IIDSitePrior(IntegSeg(), DiagonalVonMises(0.0, inv(0.5^2)));
+        refant = SEFDReference(0.0),
+        phase = true,
+    ),
+)
+intmodel = InstrumentModel(G, intpr)
+
+
+# Since we are fitting closures we do not need to include an instrument model, since
+# the closure likelihood is approximately independent of gains in the high SNR limit.
+using Enzyme
+post_cpu = VLBIPosterior(skym, intmodel, dvis)
+
+# Move every array Reactant cares about (sky-model grid + metadata, data tables, etc.)
+# onto the Reactant device. `prepare_device` walks the posterior tree and replaces the
+# concrete arrays with `Reactant.RArray`s, leaving priors and the (empty) instrument
+# model untouched.
+post = Comrade.prepare_device(post_cpu, ComradeBase.ReactantEx())
+
+# ## Reconstructing the Image
+
+# Because the optimizer needs to traverse the entire parameter space (network weights +
+# `σ`, `fb`) we work on the unconstrained-flat representation produced by `asflat`. This
+# wraps the posterior in a `TransformedVLBIPosterior` whose log-density takes a flat
+# `Vector{Float64}` and includes the log-Jacobian of the bijection.
+tpost = asflat(post)
+
+# Draw initial parameters from the prior, push them through the inverse transform to get a
+# flat vector, and visualize the corresponding initial image. `skymodel` traces through
+# Reactant because `post` lives on the device, so we `@jit` the call and then materialize
+# the underlying raster back onto the host before handing it to Makie.
+using CairoMakie, DisplayAs
+xinit_nt = prior_sample(rng, post) ## NamedTuple in the original parameter space
+init_img_r = @jit skymodel(post, Reactant.to_rarray(xinit_nt))
+init_img = Comrade.Adapt.adapt(Array, init_img_r.img)
+imageviz(init_img, size = (500, 400), colorscale = log10, colorrange = (1.0e-6, 1.0e-4)) |> DisplayAs.PNG |> DisplayAs.Text
+
+# Move the parameter vector onto the Reactant device.
+xinit_r = Reactant.to_rarray(Comrade.inverse(tpost, xinit_nt))
+
+# We use [Optimisers.jl](https://github.com/FluxML/Optimisers.jl) directly (rather than the
+# `comrade_opt` wrapper) because Reactant compiles the whole step — including the optimizer
+# update — into a single XLA program. AdamW with a learning rate of `3.0e-4` matches the
+# Enzyme/CPU version of this tutorial.
+using Optimisers
+opt = Optimisers.AdamW(3.0e-3)
+opt_state = @jit Optimisers.setup(opt, xinit_r)
+
+# A single training step: compute the negative log-density and its Enzyme reverse-mode
+# gradient, then apply one Optimisers.jl update. Reactant traces this whole function once
+# and reuses the compiled XLA program for every iteration.
+function step(tpost, x, opt_state)
+    dx = Enzyme.make_zero(x)
+    _, primal = Enzyme.autodiff(
+        Enzyme.WithPrimal(Enzyme.set_strong_zero(Enzyme.Reverse)),
+        Comrade.logdensityof, Active, Const(tpost), Duplicated(x, dx)
+    )
+    ## Optimisers expects gradients of the loss we are *minimising*; logdensityof returns
+    ## the log-posterior so we negate.
+    grads = -1 .* dx
+    new_state, new_x = Optimisers.update(opt_state, x, grads)
+    return new_state, new_x, -primal
+end
+
+# Compile the step. The first call is the (slow) MLIR/XLA compile; every subsequent call is
+# a single dispatched program.
+step_jit = @compile sync = true step(tpost, xinit_r, opt_state)
+
+# Run the optimizer. We log the loss every 250 steps so the docs build still shows
+# convergence without flooding stdout.
+xcur = xinit_r
+state_cur = opt_state
+maxiters = 10_000
+for i in 1:maxiters
+    state_cur, xcur, loss = step_jit(tpost, xcur, state_cur)
+    if i == 1 || i % 250 == 0 || i == maxiters
+        @info "iter $i  -logdensity = $(Reactant.@allowscalar Float64(loss))"
+    end
+end
+
+# Pull the optimized parameters back to the host as a NamedTuple in the original space.
+xopt = @jit Comrade.transform(tpost, xcur);
+
+using CairoMakie
+using DisplayAs #hide
+# We can plot the MAP image using `imageviz`.
+imgmap = parent(@jit skymodel(post, xopt));
+crange = (1.0e-7, 1.0e-3)
+fig = imageviz(Comrade.Adapt.adapt(Array, imgmap), colorscale = log10, colorrange = crange, size = (650, 500));
+DisplayAs.Text(DisplayAs.PNG(fig)) #hide
+
+
+# To see how well the MAP estimate fits the data we can plot the residuals.
+xopt_cpu = Comrade.transform(asflat(post_cpu), Array(xcur))
+res = @jit(Comrade.residuals(post, xopt))[1]
+res_cpu = Comrade.Adapt.adapt(Array, res)
+fig = Figure(; size = (800, 300))
+plotfields!(fig[1, 1], res_cpu, uvdist, :res);
+fig |> DisplayAs.PNG |> DisplayAs.Text
+
+
+# Overall, the image looks reasonable. Unlike other models in `Comrade`, estimating posterior uncertainty
+# with neural fields is highly challenging so we skip that step here.
+
+# We can also compare the Comrade neural field reconstruction to the CLEAN reconstruction of the same data.
+cleanf = Base.download("https://www.bu.edu/blazars/VLBA_GLAST/1308/1308+326Q.2021-03-19.IMAP.gz")
+# By default this will load the clean components with the beam defined in the FITS header.
+mcl = load_clean_components(cleanf)
+# We can also choose the load the clean components with a user-defined beam.
+mcl_25 = load_clean_components(cleanf, modify(Gaussian(), Stretch(beamsize(dvis) / 4 / fwhmfac)))
+
+# Now we can produce the CLEAN images on the same grid as our Comrade reconstruction.
+cleanimg = intensitymap(mcl, grid)
+cleanimg25 = intensitymap(mcl_25, grid)
+
+fig = Figure(; size = (900, 350));
+axs = [Axis(fig[1, j], xreversed = true, aspect = DataAspect()) for j in 1:3]
+crange = (1.0e-4, 1.0e-1)
+image!(axs[1], Comrade.Adapt.adapt(Array, imgmap), colormap = :afmhot, colorscale = log10, colorrange = crange); axs[1].title = "Comrade Neural Field"
+image!(axs[2], max.(cleanimg, 1.0e-20), colormap = :afmhot, colorscale = log10, colorrange = crange); axs[2].title = "CLEAN (Nominal beam)"
+image!(axs[3], max.(cleanimg25, 1.0e-20), colormap = :afmhot, colorscale = log10, colorrange = crange); axs[3].title = "CLEAN (25% beam)"
+hidedecorations!.(axs)
+fig |> DisplayAs.PNG |> DisplayAs.Text
+
+# From the plot you can see that the neural field reconstruction appears to capture a lot of fine detail.
+# Interpreting the reliability of these features is rather challenging without either 1 posterior uncertainty
+# estimates or 2 knowledege of what these sources look like. However, this example demonstrates the flexibility of neural fields
+# for VLBI imaging.
+#
+# A disclaimer is that we have not made any attempt to optimize the neural network architecture or hyperparameters here.
+# For a more thorough and insightful analysis a paper by Foschi et al is in preparation.
diff --git a/examples/beginner/GeometricModeling/main.jl b/examples/beginner/GeometricModeling/main.jl
index 09ab65598..f1008f7c5 100644
--- a/examples/beginner/GeometricModeling/main.jl
+++ b/examples/beginner/GeometricModeling/main.jl
@@ -42,13 +42,14 @@ uvd = VLBIFiles.load(
 )
 
 # Since the gains in this dataset are coherent across a scan, we extract scan-averaged
-# closures.
+# closures. Since uvfits aren't guaranteed to have the scan table we use a gap-based
+# scan averaging scheme.
 dlcamp, dcphase = extract_table(
     uvd,
     LogClosureAmplitudes(; time_average = VLBI.GapBasedScans()),
     ClosurePhases(; time_average = VLBI.GapBasedScans()),
 )
-# Kill the trivial 0-baseline triangles (we don't care about large-scale flux) and inflate
+# We also remove the trivial 0-baseline triangles (we don't care about large-scale flux) and inflate
 # the noise by 2% in quadrature to account for residual systematics.
 isshort(d) = any(b -> hypot(b.U, b.V) < 0.1e9, d.baseline)
 dlcamp = flag(isshort, dlcamp)
@@ -72,18 +73,18 @@ add_fractional_noise!(dcphase, 0.02)
 using Distributions, VLBIImagePriors
 
 @sky function sky(grid)
-    radius ~ Uniform(μas2rad(10.0), μas2rad(30.0))
-    width ~ Uniform(μas2rad(1.0), μas2rad(10.0))
-    ma ~ (Uniform(0.0, 0.5),)
-    mp ~ (Uniform(0, 2π),)
-    τ ~ Uniform(0.0, 1.0)
-    ξτ ~ Uniform(0.0, π)
-    f ~ Uniform(0.0, 1.0)
-    σG ~ Uniform(μas2rad(1.0), μas2rad(100.0))
-    τG ~ Exponential(1.0)
-    ξG ~ Uniform(0.0, 1π)
-    xG ~ Uniform(-μas2rad(80.0), μas2rad(80.0))
-    yG ~ Uniform(-μas2rad(80.0), μas2rad(80.0))
+    radius ~ VLBIUniform(μas2rad(10.0), μas2rad(30.0))
+    width ~ VLBIUniform(μas2rad(1.0), μas2rad(10.0))
+    ma ~ (VLBIUniform(0.0, 0.5),)
+    mp ~ (VLBIUniform(0, 2π),)
+    τ ~ VLBIUniform(0.0, 1.0)
+    ξτ ~ VLBIUniform(0.0, π)
+    f ~ VLBIUniform(0.0, 1.0)
+    σG ~ VLBIUniform(μas2rad(1.0), μas2rad(100.0))
+    τG ~ VLBIExponential(1.0)
+    ξG ~ VLBIUniform(0.0, 1π)
+    xG ~ VLBIUniform(-μas2rad(80.0), μas2rad(80.0))
+    yG ~ VLBIUniform(-μas2rad(80.0), μas2rad(80.0))
     α = ma .* cos.(mp .- ξτ)
     β = ma .* sin.(mp .- ξτ)
     ring = f * smoothed(modify(MRing(α, β), Stretch(radius, radius * (1 + τ)), Rotate(ξτ)), width)
@@ -91,19 +92,22 @@ using Distributions, VLBIImagePriors
     return ring + g
 end
 
-# Each `name ~ dist` line above contributes an entry to the prior, while every other line
-# is part of the model body. Note that for `α` and `β` we use a product distribution to
-# signify that we want to use a multivariate uniform for the mring components.
+# The syntax here follows standard probablisitic notation. Namely `~` denotes that the
+# parameter on the left is distributed according to the distribution on the right.
+# These are essentially the parameters of our model. For example, `radius` is the radius
+# of the ring and is distributed according to a uniform distribution between 10 and 30 μas.
+# Other variables like `α` and `β` are deterministic functions of the parameters on the left-hand side of
+# the `~`. The model function must return an object that implements the `VLBISkyModel` interface.
+# In this case, we return a sum of a smoothed MRing and a Gaussian.
 
-# We can now construct our Sky model, which typically takes a grid. Note that since our
-# model is analytic the grid is not directly used when computing visibilities.
+# The `@sky` macro then transforms this into a function called `sky` that we can then use to
+# initialize or construct our sky model for some given grid.
 skym = sky(imagepixels(μas2rad(200.0), μas2rad(200.0), 128, 128))
 
 
-# In this tutorial we will be using closure products as our data. As such we do not need to specify a
-# instrument model, since for stokes I imaging, the likelihood is approximately invariant to the instrument
-# model.
-post = VLBIPosterior(skym, dlcamp, dcphase)
+# Since we are fitting closures this is essentialy the entire model. In other tutorials we will
+# show how to model the instrument directly, e.g., gains.
+post = VLBIPosterior(skym, dlcamp, dcphase);
 
 # !!! note
 #     When fitting visibilities a instrument is required, and a reader can refer to
@@ -142,12 +146,12 @@ logdensityof(
 # parameters in the unit hypercube. To transform the posterior to the unit hypercube, we
 # can use the `ascube` function
 
-cpost = ascube(post)
+cpost = ascube(post);
 
 # If we want to flatten the parameter space and move from constrained parameters to (-∞, ∞)
 # support we can use the `asflat` function
 
-fpost = asflat(post)
+fpost = asflat(post);
 
 # These transformed posterior expect a vector of parameters. For example, we can draw from the
 # prior in our usual parameter space
@@ -157,7 +161,11 @@ p = prior_sample(rng, post)
 logdensityof(cpost, Comrade.inverse(cpost, p))
 logdensityof(fpost, Comrade.inverse(fpost, p))
 
-# note that the log densit is not the same since the transformation has causes a jacobian to ensure volume is preserved.
+# note that the log densit is not the same since the transformation has causes a
+# jacobian to ensure volume is preserved. Note that this is rather critical because it
+# means that the maximum a posteriori (MAP) estimate in bother of these examples will differ.
+# In general, the MAP estimate is parameterization dependent. This is not the case for estimates
+# that are derived from expectations of the posterior, e.g., the mean image.
 
 # ### Finding the Optimal Image
 
diff --git a/examples/intermediate/ClosureImaging/main.jl b/examples/intermediate/ClosureImaging/main.jl
index b2d079c46..1b913bbab 100644
--- a/examples/intermediate/ClosureImaging/main.jl
+++ b/examples/intermediate/ClosureImaging/main.jl
@@ -75,8 +75,8 @@ add_fractional_noise!(dcphase, 0.02)
 using VLBIImagePriors, Distributions
 @sky function sky(grid; mimg, cprior)
     c ~ cprior
-    σimg ~ Exponential(0.1)
-    fg ~ Uniform(0.0, 1.0)
+    σimg ~ VLBIExponential(0.1)
+    fg ~ VLBIUniform(0.0, 1.0)
     ## Apply the GMRF fluctuations to the image
     rast = apply_fluctuations(CenteredLR(), mimg, σimg .* c.params)
     m = ContinuousImage(((1 - fg)) .* rast, BSplinePulse{3}())
@@ -84,29 +84,28 @@ using VLBIImagePriors, Distributions
     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.
+# Let's explain this model a bit. In general Comrade aims to provide an extremely flexible set of
+# possible image models to consider. The basic image model is `ContinuousImage`, which is a raster of point sources
+# convolved with some kernel or pulse. The parameters of the model are the fluxes of the point sources,
+# which are given by `c.params` in the code above. This is essentially identical for every imaging model we consider.
+# The only difference between different image models is then the prior we place on the pixel fluxes.
+# Other tutorials will consider a vast array of different image priors, from Gaussian processes like Matern kernels to
+# neural fields. In this tutorial we will a very simple Gaussian Markov random field (GMRF) prior, which is a type of
+# Gaussian process that is very fast to sample from and evaluate. Note that our prior actually lives in the log-ratio space.
+# We do this to 1 ensure positivity of the image and 2 ensure that the total flux is fixed to some value. This ensures that
+# we have more directly control over a classic VLBI degeneracy the total flux of the image, which actually is not constrained
+# by closures.
+
+
+# To define the GMRF we first specify our grid on which the GMRF is defined.
+# Given that M87* is compact and the EHT is not very sensitive in 2018 we use
+# a small FOV
 npix = 32
 fovxy = μas2rad(150.0)
-
-# To define the image model we need to specify both the grid we will be using and the
-# FT algorithm we will use, in this case the NFFT which is the most efficient.
 grid = imagepixels(fovxy, fovxy, npix, npix)
 
-
-# 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 50 μas
-fwhmfac = 2 * sqrt(2 * log(2))
-mpr = modify(Gaussian(), Stretch(μas2rad(50.0) ./ fwhmfac))
-imgpr = intensitymap(mpr, grid)
-mimg = imgpr ./ Comrade.flux(imgpr);
-
-
-# Now we can finally form our image prior. For this we use a heirarchical prior where the
+# Given this grid we can now define our GMRF prior called `cprior` above.
+# For this we use a heirarchical prior where the
 # direct log-ratio image prior is a Gaussian Markov Random Field. The correlation length
 # of the GMRF is a hyperparameter that is fit during imaging. We pass the data to the prior
 # to estimate what the maximumal resolutoin of the array is and prevent the prior from allowing
@@ -114,13 +113,25 @@ mimg = imgpr ./ Comrade.flux(imgpr);
 # has unit variance. For more information on the GMRF prior see the [`corr_image_prior`](@ref) doc string.
 cprior = corr_image_prior(grid, dlcamp)
 
-# We can then define our sky model.
+# The other aspect of the image model is the mean image. By mean image we mean the image structure
+# about which the kind of fluctuations occur. In this case we can view the GMRF flucations `c` as a
+# kind of multiplicative turbulence aboue our apriori mean structure. For the EHT we will follow the
+# original publications and essentially assume that this mean structure is a Gaussian. For simplicitly in
+# in this tutorial we assume that the Gaussian is symmetric with a FWHM of 50 μas which is rougly twice
+# the beam of the EHT. In reality we could also easily fit for the size of this Gaussian.
+fwhmfac = 2 * sqrt(2 * log(2))
+mpr = modify(Gaussian(), Stretch(μas2rad(50.0) ./ fwhmfac))
+imgpr = intensitymap(mpr, grid)
+mimg = imgpr ./ Comrade.flux(imgpr);
+
+
+# Given these two ingredients we can now construct our sky model.
 skym = sky(grid; mimg, cprior)
 
-# Since we are fitting closures we do not need to include an instrument model, since
-# the closure likelihood is approximately independent of gains in the high SNR limit.
+# At this point since we are fitting closures we have essentially finished our model specification and can form
+# our VLBIPosterior.
 using Enzyme
-post = VLBIPosterior(skym, dlcamp, dcphase)
+post = VLBIPosterior(skym, dlcamp, dcphase);
 
 # ## Reconstructing the Image
 
diff --git a/examples/intermediate/PolarizedImaging/Project.toml b/examples/intermediate/PolarizedImaging/Project.toml
index db7939465..431bcc227 100644
--- a/examples/intermediate/PolarizedImaging/Project.toml
+++ b/examples/intermediate/PolarizedImaging/Project.toml
@@ -16,4 +16,5 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 VLBIFiles = "c1ebf4c8-f9d4-409a-8daf-7009448f4e6e"
 VLBIImagePriors = "b1ba175b-8447-452c-b961-7db2d6f7a029"
+VLBILikelihoods = "90db92cd-0007-4c0a-8e51-dbf0782ce592"
 VLBISkyModels = "d6343c73-7174-4e0f-bb64-562643efbeca"
diff --git a/examples/intermediate/PolarizedImaging/main.jl b/examples/intermediate/PolarizedImaging/main.jl
index 0e16ee0d2..ff6463d05 100644
--- a/examples/intermediate/PolarizedImaging/main.jl
+++ b/examples/intermediate/PolarizedImaging/main.jl
@@ -122,6 +122,9 @@ dvis = extract_table(
         time_average = VLBI.GapBasedScans(),
     )
 )
+
+# Often the EHT does not take into acount feed offsets. To load those we read in the standard
+# ehtim style array.txt files.
 Comrade.reset_mounts!(dvis, Comrade.load_array_txt(joinpath(__DIR, "../../Data/array.txt")))
 # Inflate noise by 1% and drop short (uvdist < 0.1 Gλ) baselines.
 add_fractional_noise!(dvis, 0.01)
@@ -154,28 +157,37 @@ dvis = flag(d -> hypot(d.baseline.U, d.baseline.V) < 0.1e9, dvis)
 using Distributions
 using VLBIImagePriors
 @sky function sky(grid; mimg, ftot, cprior)
-    σs ~ ntuple(Returns(truncated(Normal(0.0, 0.5); lower = 0.0)), 4)
+    σs ~ ntuple(Returns(VLBITruncated(VLBIGaussian(0.0, 0.5); lower = 0.0)), 4)
     as ~ ntuple(Returns(cprior), 4)
     ## Build the stokes I model
     δs = ntuple(Val(4)) do i
         σs[i] .* as[i].params
     end
 
-    ## Convert hyperbolic polarization basis to Stokes basis
+    ## Convert hyperbolic polarization parameterization to Stokes parameters.
     pmap = VLBISkyModels.PolExp2Map!(δs..., axisdims(mimg))
 
     ## We now add a mean image. Namely, we assume that `pmap` are multiplicative fluctuations
     ## about some mean image `mimg`. We also compute the total flux of the Stokes I image
     ## for normalization purposes below.
     bpmap = baseimage(pmap)
+    bmimg = baseimage(mimg)
     ft = zero(ftot)
-    for i in eachindex(bpmap)
-        bpmap[i] = mimg[i] * bpmap[i]
-        ft += bpmap[i].I
+    bpmapI = stokes(bpmap, :I)
+    bpmapQ = stokes(bpmap, :Q)
+    bpmapU = stokes(bpmap, :U)
+    bpmapV = stokes(bpmap, :V)
+
+    @inbounds for i in eachindex(bpmap, bmimg)
+        bpmapI[i] = bmimg[i] * bpmapI[i]
+        bpmapQ[i] = bmimg[i] * bpmapQ[i]
+        bpmapU[i] = bmimg[i] * bpmapU[i]
+        bpmapV[i] = bmimg[i] * bpmapV[i]
+        ft += bpmapI[i]
     end
 
     ## Renormlize to get the correct total flux.
-    for i in eachindex(bpmap)
+    @inbounds for i in eachindex(bpmap)
         bpmap[i] *= ftot / ft
     end
 
@@ -230,8 +242,8 @@ skym = sky(grid; mimg, ftot = 0.6, cprior)
 # The first element of the tuple is the gain for the first polarization feed (R) and the
 # second is the gain for the second polarization feed (L).
 function fgain(x)
-    gR = exp(x.lgR + 1im * x.gpR)
-    gL = gR * exp(x.lgrat + 1im * x.gprat)
+    gR = exp(complex(x.lgR, x.gpR))
+    gL = gR * exp(complex(x.lgrat, x.gprat))
     return gR, gL
 end
 G = JonesG(fgain)
@@ -259,7 +271,7 @@ R = JonesR(; add_fr = true)
 # first argument is a function that specifies how to combine each Jones matrix. In this case
 # we will use the standard decomposition J = adjoint(R)*G*D*R, where we need to apply the adjoint
 # of the feed rotaion matrix `R` because the data has feed rotation calibration.
-js(g, d, r) = adjoint(r) * g * d * r
+@inline js(g, d, r) = adjoint(r) * g * d * r
 J = JonesSandwich(js, G, D, R)
 
 # !!! note
@@ -283,14 +295,14 @@ J = JonesSandwich(js, G, D, R)
 # so we use `ScanSeg` for those quantities. The d-terms are typically stable over the track
 # so we use `TrackSeg` for those.
 intprior = (
-    lgR = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.2)); LM = IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
-    lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
+    lgR = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.2)); LM = IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0))),
+    lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
     gpR = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
     gprat = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(0.1^2))); refant = SingleReference(:AA, 0.0), phase = true),
-    dRx = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-    dRy = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-    dLx = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-    dLy = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
+    dRx = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+    dRy = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+    dLx = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+    dLy = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
 )
 
 # Finally, we can build our instrument model which takes a model for the Jones matrix `J`
@@ -301,14 +313,14 @@ intmodel = InstrumentModel(J, intprior)
 # Putting it all together, we form our likelihood and posterior objects for optimization and
 # sampling, and specifying to use Enzyme.Reverse with runtime activity for AD.
 using Enzyme
-post = VLBIPosterior(skym, intmodel, dvis)
+post = VLBIPosterior(skym, intmodel, dvis);
 
 # ## Reconstructing the Image and Instrument Effects
 
 # To sample from this posterior, it is convenient to move from our constrained parameter space
 # to an unconstrained one (i.e., the support of the transformed posterior is (-∞, ∞)). This transformation is
 # done using the `asflat` function.
-tpost = asflat(post)
+tpost = asflat(post);
 
 # We can also query the dimension of our posterior or the number of parameters we will sample.
 # !!! warning
diff --git a/examples/intermediate/StokesIImaging/main.jl b/examples/intermediate/StokesIImaging/main.jl
index b4d4dbc70..124ccf366 100644
--- a/examples/intermediate/StokesIImaging/main.jl
+++ b/examples/intermediate/StokesIImaging/main.jl
@@ -59,8 +59,8 @@ using VLBIImagePriors
 using Distributions
 @sky function sky(grid; ftot, mimg, cprior)
     c ~ cprior
-    σimg ~ truncated(Normal(0.0, 0.5); lower = 0.0)
-    fg ~ Uniform(0.0, 1.0)
+    σimg ~ VLBITruncated(VLBIGaussian(0.0, 0.5); lower = 0.0)
+    fg ~ VLBIUniform(0.0, 1.0)
     ## Apply the GMRF fluctuations to the image
     rast = apply_fluctuations(CenteredLR(), mimg, σimg .* c.params)
     pimg = parent(rast)
@@ -114,11 +114,11 @@ skym = sky(grid; ftot = 1.1, mimg, cprior)
 #   - Gain amplitudes which are typically known to 10-20%, except for LMT, which has amplitudes closer to 50-100%.
 #   - Gain phases which are more difficult to constrain and can shift rapidly.
 
-gain(x) = exp(x.lg + 1im * x.gp)
+gain(x) = exp(complex(x.lg, x.gp))
 G = SingleStokesGain(gain)
 
 intpr = (
-    lg = ArrayPrior(IIDSitePrior(IntegSeg(), Normal(0.0, 0.2)); LM = IIDSitePrior(IntegSeg(), Normal(0.0, 1.0))),
+    lg = ArrayPrior(IIDSitePrior(IntegSeg(), VLBIGaussian(0.0, 0.2)); LM = IIDSitePrior(IntegSeg(), VLBIGaussian(0.0, 1.0))),
     gp = ArrayPrior(IIDSitePrior(IntegSeg(), DiagonalVonMises(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
 )
 intmodel = InstrumentModel(G, intpr)
@@ -128,7 +128,7 @@ intmodel = InstrumentModel(G, intpr)
 # gradients of the posterior we also need to load `Enzyme.jl`. Under the hood, Comrade will use
 # Enzyme to compute the gradients of the posterior.
 using Enzyme
-post = VLBIPosterior(skym, intmodel, dvis)
+post = VLBIPosterior(skym, intmodel, dvis);
 
 # ## Optimization and Sampling
 
@@ -186,7 +186,7 @@ plotcaltable(abs.(intopt)) |> DisplayAs.PNG |> DisplayAs.Text
 # instrument modeling as we are able to solve for the gain amplitudes and get a reasonable
 # image.
 
-
+# # Why we should Sample from the Posterior
 # One problem with the MAP estimate is that it does not provide uncertainties on the image.
 # That is we are unable to statistically assess which components of the image are certain.
 # Comrade is really a Bayesian imaging and calibration package for VLBI. Therefore, our
diff --git a/ext/ComradeReactantExt.jl b/ext/ComradeReactantExt.jl
new file mode 100644
index 000000000..38c3739b5
--- /dev/null
+++ b/ext/ComradeReactantExt.jl
@@ -0,0 +1,93 @@
+module ComradeReactantExt
+using Comrade
+using ComradeBase: ReactantEx
+using Reactant
+using StructArrays
+using StaticArraysCore
+using Accessors
+using PolarizedTypes
+
+@inline function Comrade._apply_instrument!(
+        vout::Union{Reactant.AnyTracedRArray, AbstractArray{<:SMatrix{2, 2, <:Reactant.TracedRNumber}}},
+        vis,
+        J::Comrade.ObservedInstrumentModel,
+        xint
+    )
+    Reactant.@allowscalar @trace track_numbers = false for i in eachindex(vis)
+        vout[i] = Comrade.apply_jones(vis[i], i, J, xint)
+    end
+    return vout
+end
+
+function StructArrays.createinstance(::Type{<:StokesParams}, args...)
+    return StokesParams(args...)
+end
+
+@inline function Comrade.site_sum(y::Reactant.AnyTracedRArray, site_map::Comrade.SiteLookup)
+    yout = similar(y)
+    vals = values(lookup(site_map))
+    #unroll this?
+    for site in vals
+        ys = y[site]
+        yout[site] = cumsum(ys)
+    end
+    return yout
+end
+
+# @inline function Comrade.site_sum(y::Reactant.AnyTracedRArray, site_map::Comrade.SiteLookup)
+#     yout = similar(y)
+#     vals = values(lookup(site_map))
+#     #unroll this?
+#     @trace track_numbers=false for i in eachindex(vals)
+#         site = vals[i]
+#         ys = y[site]
+#         acc = zero(eltype(y))
+#         @inbounds ys = @view y[site]
+#         for i in eachindex(ys)
+#             acc += ys[i]
+#             ys[i] = acc
+#         end
+#         yout[site] = ys
+#     end
+#     return yout
+# end
+
+
+function Comrade.branchcut(x::Reactant.TracedRNumber)
+    xmod = mod(x, oftype(x, 2π))
+    return ifelse(xmod > π, xmod - oftype(x, 2π), xmod)
+end
+
+function Comrade.prepare_device(post::VLBIPosterior, ex::ReactantEx)
+    ## TODO add Sharding and other options here
+    for p in propertynames(post)
+        p == :prior && continue
+        post = Accessors.set(post, PropertyLens(p), Comrade.prepare_device(getproperty(post, p), ex))
+    end
+    return post
+end
+
+function Comrade.prepare_device(m, ex::ReactantEx)
+    return Reactant.to_rarray(m)
+end
+
+function Comrade.prepare_device(m::Comrade.ObservedSkyModel, ex::ReactantEx)
+    grid = Comrade.prepare_device(m.grid, ex)
+    return Comrade.ObservedSkyModel(m.f, grid, Reactant.to_rarray(m.metadata))
+end
+
+function Comrade.prepare_device(m::Comrade.ObservedInstrumentModel, ex::ReactantEx)
+    return Comrade.ObservedInstrumentModel(Reactant.to_rarray(m.instrument), m.refbasis, m.bsitelookup)
+end
+
+function Comrade.prepare_device(grid::VLBISkyModels.FourierDualDomain, ex::ReactantEx)
+    gimr = @jit identity(grid.imgdomain)
+    guvr = Reactant.to_rarray(grid.visdomain)
+    if grid.algorithm isa VLBISkyModels.ReactantNUFFTAlg
+        algr = grid.algorithm
+    else
+        algr = VLBISkyModels.ReactantNUFFTAlg(eltype(grid.imgdomain))
+    end
+    return VLBISkyModels.FourierDualDomain(gimr, guvr, algr)
+end
+end
diff --git a/ext/ComradeVLBIFilesExt.jl b/ext/ComradeVLBIFilesExt.jl
index 18bbab413..7430035f8 100644
--- a/ext/ComradeVLBIFilesExt.jl
+++ b/ext/ComradeVLBIFilesExt.jl
@@ -160,8 +160,8 @@ function _arrayconfig(
     Ti = Vector{Float64}(undef, n)
     Fr = Vector{Float64}(undef, n)
     sites_v = Vector{Tuple{Symbol, Symbol}}(undef, n)
-    elevation = Vector{Tuple{Float64, Float64}}(undef, n)
-    parallactic = Vector{Tuple{Float64, Float64}}(undef, n)
+    elevation = StructVector{Tuple{Float64, Float64}}(undef, n)
+    parallactic = StructVector{Tuple{Float64, Float64}}(undef, n)
 
     # caches built from antenna metadata
     site_xyz = Dict{Symbol, NTuple{3, Float64}}()
@@ -235,7 +235,7 @@ function _prep_uvtable(
         uvtbl;
         time_average = nothing, frequency_average = true
     )
-    scan_table = VLBIData.scan_intervals(time_average, uvtbl)
+    scan_table = VLBIData.scan_intervals(VLBI.GapBasedScans(), uvtbl)
     rows = uvtbl
     if frequency_average !== nothing && frequency_average !== false
         rows = VLBI.average_data(VLBI.ByFrequency(), rows)
@@ -244,7 +244,6 @@ function _prep_uvtable(
         rows = VLBI.average_data(time_average, rows)
     end
 
-
     return rows, scan_table
 end
 
diff --git a/src/Comrade.jl b/src/Comrade.jl
index aad90e5da..e7dcc15cb 100755
--- a/src/Comrade.jl
+++ b/src/Comrade.jl
@@ -1,5 +1,6 @@
 module Comrade
 
+using Adapt
 using AbstractMCMC
 using Accessors: @set
 using ArgCheck: @argcheck
@@ -29,6 +30,8 @@ using Tables
 import TransformVariables as TV
 using ComradeBase: AbstractDomain, AbstractSingleDomain, AbstractRectiGrid
 using VLBISkyModels: FourierTransform, FourierDualDomain
+using ReactantCore
+using ComradeBase: rgetindex, rsetindex!
 
 # Reexport the core libraries for Comrade
 @reexport using VLBISkyModels
diff --git a/src/instrument/instrument_transforms.jl b/src/instrument/instrument_transforms.jl
index 130fe002a..a2336b963 100644
--- a/src/instrument/instrument_transforms.jl
+++ b/src/instrument/instrument_transforms.jl
@@ -7,7 +7,8 @@ ascube(t::AbstractInstrumentTransform) = t
 function TV.transform_with(flag::TV.LogJacFlag, m::AbstractInstrumentTransform, x, index)
     y, ℓ, index = _instrument_transform_with(flag, m, x, index)
     sm = site_map(m)
-    return SiteArray(y, sm), ℓ, index
+    sa = SiteArray(y, sm)
+    return sa, ℓ, index
 end
 
 EnzymeRules.inactive_type(::Type{AbstractArray{<:IntegrationTime}}) = true
diff --git a/src/instrument/jonesmatrices.jl b/src/instrument/jonesmatrices.jl
index 5848513fd..2bf042f26 100644
--- a/src/instrument/jonesmatrices.jl
+++ b/src/instrument/jonesmatrices.jl
@@ -47,8 +47,8 @@ is the complex gain `g1` and the second element is the complex gain `g2`.
 ## Example
 ```julia
 G = JonesG() do
-    g1 = exp(x.lg1 + 1im.*x.gp1)
-    g2 = g1*exp(x.lgratio + 1im.*x.gpratio)
+    g1 = exp(complex(x.lg1, x.gp1))
+    g2 = g1*exp(complex(x.lgratio, x.gpratio))
     return g1, g2
 end
 ```
@@ -112,7 +112,16 @@ end
 struct GenericJones{F} <: AbstractJonesMatrix
     param_map::F
 end
-Base.@propagate_inbounds construct_jones(::GenericJones, x::NTuple{4, T}, index, site) where {T} = SMatrix{2, 2, T, 4}(x[1], x[2], x[3], x[4])
+@inline construct_jones(::GenericJones, x::NTuple{4, T}, index, site) where {T} = SMatrix{2, 2, T, 4}(x[1], x[2], x[3], x[4])
+
+
+struct JonesConst{M} <: AbstractJonesMatrix
+    matrices::M
+end
+JonesConst(m1, m2) = JonesConst((m1, m2))
+@inline construct_jones(J::JonesConst, x, index, ::Val{N}) where {N} = @inbounds(rgetindex(J.matrices[N], index))
+param_map(::JonesConst, x) = x
+
 
 """
     JonesF(;add_fr=true)
@@ -128,13 +137,14 @@ struct JonesF{M} <: AbstractJonesMatrix
     matrices::M
 end
 JonesF() = JonesF(nothing)
-Base.@propagate_inbounds construct_jones(J::JonesF, x, index, ::Val{M}) where {M} = J.matrices[index][M]
+@inline construct_jones(J::JonesF, x, index, ::Val{M}) where {M} = @inbounds J.matrices[index][M]
 param_map(::JonesF, x) = x
 function preallocate_jones(::JonesF, array::AbstractArrayConfiguration, ref)
     field_rotations = build_feedrotation(array)
-    return JonesF(field_rotations)
+    return JonesConst(field_rotations[1], field_rotations[2])
 end
 
+
 """
     JonesR(;add_fr=true)
 
@@ -149,7 +159,7 @@ Base.@kwdef struct JonesR{M} <: AbstractJonesMatrix
     matrices::M = nothing
     add_fr::Bool = true
 end
-Base.@propagate_inbounds construct_jones(J::JonesR, x, index, ::Val{M}) where {M} = J.matrices[M][index]
+@inline construct_jones(J::JonesR, x, index, ::Val{M}) where {M} = @inbounds rgetindex(J.matrices[M], index)
 param_map(::JonesR, x) = x
 
 function preallocate_jones(J::JonesR, array::AbstractArrayConfiguration, ref)
@@ -162,7 +172,7 @@ function preallocate_jones(J::JonesR, array::AbstractArrayConfiguration, ref)
         @. T1 .= Tcirc1 * field_rotations[1] * adjoint(Tcirc1) * T1
         @. T2 .= Tcirc2 * field_rotations[2] * adjoint(Tcirc2) * T2
     end
-    return JonesR((T1, T2), J.add_fr)
+    return JonesConst(T1, T2)
 
 end
 
diff --git a/src/instrument/model.jl b/src/instrument/model.jl
index af940bcaf..ab3f8b76a 100644
--- a/src/instrument/model.jl
+++ b/src/instrument/model.jl
@@ -106,8 +106,8 @@ visibilities.
 
 A Stokes I example is
 ```julia-repl
-julia> G = SingleStokesGain(x->exp(x.lg + 1im*x.pg))
-julia> intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
+julia> G = SingleStokesGain(x->exp(complex(x.lg, x.pg)))
+julia> intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
             pg = ArrayPrior(IIDSitePrior(ScanSeg(), DiagVonMises(0.0, inv(π^2))))
             )
 
@@ -117,8 +117,8 @@ julia> intm = InstrumentModel(G, intprior)
 A standard polarized example is
 ```julia-repl
 julia> G = JonesG() do
-        gR = exp.(x.lgr + 1im*x.gpr)
-        gL = gr*exp.(x.lgrat + 1im*x.gprat)
+        gR = exp.(complex(x.lgr, x.gpr))
+        gL = gr*exp.(complex(x.lgrat, x.gprat))
         return gR, gL
     end
 julia> D = JonesD() do
@@ -128,14 +128,14 @@ julia> D = JonesD() do
     end
 julia> R = JonesR()
 julia> J = JonesSandwich(G, D, R)
-julia> intprior = (lgr = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)),
+julia> intprior = (lgr = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)),
                     gpr = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))),
-                    lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)),
+                    lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)),
                     gprat = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))),
-                    dRre = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.1)),
-                    dRim = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.1)),
-                    dLre = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.1)),
-                    dLim = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.1))
+                    dRre = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1)),
+                    dRim = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1)),
+                    dLre = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1)),
+                    dLim = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1))
                     )
 julia> intm = InstrumentModel(J, intprior)
 ```
@@ -144,7 +144,7 @@ which construct the gain matrix from R and ratios, and D is the small leakage ma
 is the *response matrix* that controls how the site responds to the ideal visibility in the reference
 basis.
 """
-function InstrumentModel(jones::AbstractJonesMatrix, prior::NamedTuple{N, <:NTuple{M, ArrayPrior}}; refbasis = CirBasis()) where {N, M}
+function InstrumentModel(jones::AbstractJonesMatrix, prior::NamedTuple{N}; refbasis = CirBasis()) where {N}
     return InstrumentModel(jones, prior, refbasis)
 end
 
@@ -169,7 +169,7 @@ function set_array(int::IdealInstrumentModel, ::AbstractArrayConfiguration)
     return (int, NamedDist((;)))
 end
 
-struct BaselineSiteLookup{V <: AbstractArray{<:Integer}}
+struct BaselineSiteLookup{V <: AbstractArray}
     indices_1::V
     indices_2::V
 end
@@ -178,7 +178,6 @@ function _construct_baselinemap(array::EHTArrayConfiguration, x::SiteArray)
     T = array[:Ti]
     F = array[:Fr]
     bl = array[:sites]
-
     return _construct_baselinemap(T, F, bl, x)
 end
 
@@ -206,44 +205,31 @@ function _construct_baselinemap(T, F, bl, x::SiteArray)
 end
 
 
-@inline intout(vis::AbstractArray{<:StokesParams{T}}) where {T <: Real} = similar(vis, SMatrix{2, 2, Complex{T}, 4})
-@inline intout(vis::AbstractArray{T}) where {T <: Real} = similar(vis, Complex{T})
-@inline intout(vis::AbstractArray{<:CoherencyMatrix{A, B, T}}) where {A, B, T <: Real} = similar(vis, SMatrix{2, 2, Complex{T}, 4})
-
-@inline intout(vis::AbstractArray{<:StokesParams{T}}) where {T <: Complex} = similar(vis, SMatrix{2, 2, T, 4})
-@inline intout(vis::AbstractArray{T}) where {T <: Complex} = similar(vis, T)
-@inline intout(vis::AbstractArray{<:CoherencyMatrix{A, B, T}}) where {A, B, T <: Complex} = similar(vis, SMatrix{2, 2, T, 4})
-
-intout(vis::StructArray{<:StokesParams{T}}) where {T <: Complex} = StructArray{SMatrix{2, 2, T, 4}}((vis.I, vis.Q, vis.U, vis.V))
+@inline intout(vis::AbstractArray{<:StokesParams{T}}) where {T} = similar(vis, SMatrix{2, 2, complex(T), 4})
+@inline intout(vis::AbstractArray{T}) where {T} = similar(vis, complex(T))
+@inline intout(vis::AbstractArray{<:CoherencyMatrix{A, B, T}}) where {A, B, T} = similar(vis, SMatrix{2, 2, complex(T), 4})
+@inline intout(vis::StructArray{<:StokesParams{T}}) where {T} = StructArray{SMatrix{2, 2, complex(T), 4}}((similar(vis.I), similar(vis.Q), similar(vis.U), similar(vis.V)))
 
 @inline function apply_instrument(vis, J::ObservedInstrumentModel, x)
-    vout = intout(parent(vis))
-    # Grab parent arrary so that type inference works better for Enzyme Reverse pass
+    vout = parent(intout(vis))
+    # Grab parent arrary so we don't trace through SiteArray since that stuff is constant
     xint = map(parent, x.instrument)
-    @inbounds for i in eachindex(vis, vout)
-        vout[i] = @inline apply_jones(vis[i], i, J, xint)
-    end
-    # TODO this randomly segfaults when hitting the GC if we figure out why
-    # we will revert to broadcast so it works on the GPU
-    # RJ = Ref(J)
-    # Rx = Ref(xint)
-    # vout .= apply_jones.(vis, eachindex(vis), RJ, Rx)
+    _apply_instrument!(vout, parent(vis), J, xint)
     return vout
 end
 
-# function apply_instrument(vis, J::ObservedInstrumentModel, x)
-#     xint = x.instrument
-#     vout = map(Array(vis), eachindex(vis)) do v, i
-#         return apply_jones(v, i, J, xint)
-#     end
-#     # vout = apply_jones.(vis, eachindex(vis), Ref(J), Ref(x.instrument))
-#     return UnstructuredMap(StructArray(vout), axisdims(vis))
-# end
+
+@inline function _apply_instrument!(vout::AbstractArray, vis, J::ObservedInstrumentModel, xint)
+    @inbounds for i in eachindex(vout, vis)
+        vout[i] = @inline apply_jones(vis[i], i, J, xint)
+    end
+    return
+end
 
 EnzymeRules.inactive_type(::Type{<:ObservedInstrumentModel}) = true
 
 
-@inline function apply_instrument(vis, J::ObservedInstrumentModel{<:Union{JonesR, JonesF}}, x)
+@inline function apply_instrument(vis, J::ObservedInstrumentModel{<:JonesConst}, x)
     vout = intout(parent(vis))
     vout .= apply_jones.(vis, eachindex(vis), Ref(J), Ref((;)))
     return vout
@@ -260,26 +246,31 @@ end
 #     return nothing
 # end
 
-@inline get_indices(bsitemaps, index, ::Val{1}) = map(x -> getindex(x.indices_1, index), bsitemaps)
-@inline get_indices(bsitemaps, index, ::Val{2}) = map(x -> getindex(x.indices_2, index), bsitemaps)
-@inline get_params(x::NamedTuple{N}, indices::NamedTuple{N}) where {N} = NamedTuple{N}(map(getindex, values(x), values(indices)))
+@inline get_indices(bsitemaps, index, ::Val{1}) = map(Base.Fix2(rgetindex, index) ∘ Base.Fix2(getproperty, :indices_1), bsitemaps)
+@inline get_indices(bsitemaps, index, ::Val{2}) = map(Base.Fix2(rgetindex, index) ∘ Base.Fix2(getproperty, :indices_2), bsitemaps)
+@inline get_params(x::NamedTuple{N}, indices::NamedTuple{N}) where {N} = NamedTuple{N}(map(rgetindex, values(x), values(indices)))
 # @inline get_params(x::NamedTuple{N}, indices::NamedTuple{N}) where {N} = NamedTuple{N}(ntuple(i->getindex(x[i], indices[i]), Val(length(N))))
 
 # We need this because Enzyme seems to crash when generating code for this
 # TODO try to find MWE and post to Enzyme.jl
 EnzymeRules.inactive(::typeof(get_indices), args...) = nothing
 
+@inline function unrollmap(x::NTuple{N}, inds::NTuple{N}) where {N}
+    return ntuple(Val(N)) do k
+        Base.@_inline_meta
+        @inbounds getindex(x[k], inds[k])
+    end
+end
 
-Base.@propagate_inbounds function apply_jones(v, index::Int, J::ObservedInstrumentModel, x::NamedTuple{N}) where {N}
+Base.@propagate_inbounds function apply_jones(v, index, J::ObservedInstrumentModel, x::NamedTuple{N}) where {N}
     # First lhs station
-    id1 = Base.Fix2(getindex, index) ∘ Base.Fix2(getproperty, :indices_1)
-    indices1 = map(x -> getindex(x.indices_1, index), sitelookup(J)) #get_indices(sitelookup(J), index, Val(N))
-    params1 = NamedTuple{N}(map(getindex, values(x), values(indices1)))
+    indices1 = get_indices(sitelookup(J), index, Val(1))
+    params1 = NamedTuple{N}(unrollmap(values(x), values(indices1)))
     j1 = jonesmatrix(instrument(J), params1, index, Val(1))
 
     # Second RHS station
-    indices2 = map(x -> getindex(x.indices_2, index), sitelookup(J)) #get_indices(sitelookup(J), index, Val(N))
-    params2 = NamedTuple{N}(map(getindex, values(x), values(indices2)))
+    indices2 = get_indices(sitelookup(J), index, Val(2))
+    params2 = NamedTuple{N}(unrollmap(values(x), values(indices2)))
     j2 = jonesmatrix(instrument(J), params2, index, Val(2))
 
     vout = _apply_jones(v, j1, j2, refbasis(J))
diff --git a/src/instrument/priors/array_priors.jl b/src/instrument/priors/array_priors.jl
index def1714dc..f0c4c98d5 100644
--- a/src/instrument/priors/array_priors.jl
+++ b/src/instrument/priors/array_priors.jl
@@ -23,7 +23,7 @@ recommend setting it to false currently.*
 # Example
 
 ```julia
-p = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0, 0.1)); LM = IIDSitePrior(ScanSeg(), Normal(0.0, 1.0)) refant=SEFDReference())
+p = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0, 0.1)); LM = IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0)) refant=SEFDReference())
 ```
 
 means that every site has a normal prior with mean 0 and 0.1 std. dev. except LM which is mean
@@ -49,7 +49,7 @@ struct ObservedArrayPrior{D, S} <: Distributions.ContinuousMultivariateDistribut
 end
 Base.eltype(d::ObservedArrayPrior) = eltype(d.dists)
 Base.length(d::ObservedArrayPrior) = length(d.dists)
-Dists._logpdf(d::ObservedArrayPrior, x::AbstractArray{<:Real}) = Dists.logpdf(d.dists, parent(x))
+Dists.logpdf(d::ObservedArrayPrior, x::AbstractVector{<:Number}) = Dists.logpdf(d.dists, parent(x))
 Dists._rand!(rng::Random.AbstractRNG, d::ObservedArrayPrior, x::AbstractArray{<:Real}) = SiteArray(Dists._rand!(rng, d.dists, x), d.sitemap)
 function asflat(d::ObservedArrayPrior)
     d.phase && MarkovInstrumentTransform(asflat(d.dists), d.sitemap)
@@ -69,6 +69,7 @@ function build_sitemap(d::ArrayPrior, array)
     T = array[:Ti]
     F = array[:Fr]
 
+
     # Ok to so this we are going to construct the schema first over sites.
     # At the end we may re-order depending on the schema ordering we want
     # to use.
@@ -76,7 +77,7 @@ function build_sitemap(d::ArrayPrior, array)
         seg = segmentation(sites_prior[s])
         # get all the indices where this site is present
         inds_s = findall(x -> ((x[1] == s)||x[2] == s), array[:sites])
-        # Get all the unique times
+        # Get all the times and frequencies for that site
         ts = T[inds_s]
         fs = F[inds_s]
         tfs = zip(ts, fs)
@@ -179,7 +180,7 @@ Base.eltype(d::PartiallyConditionedDist) = eltype(d.dist)
 
 Distributions.sampler(d::PartiallyConditionedDist) = d
 
-function Distributions._logpdf(d::PartiallyConditionedDist, x)
+function Distributions.logpdf(d::PartiallyConditionedDist, x::AbstractVector)
     xv = @view parent(x)[d.variate_index]
     return Dists.logpdf(d.dist, xv)
 end
diff --git a/src/instrument/priors/independent.jl b/src/instrument/priors/independent.jl
index 1c7501576..a30206a0b 100644
--- a/src/instrument/priors/independent.jl
+++ b/src/instrument/priors/independent.jl
@@ -14,7 +14,7 @@ segmentation is. The `dist` argument is the distribution of the site prior.
 ## Example
 
 ```julia
-A = IIDSitePrior(ScanSeg(), Normal(0, 1))
+A = IIDSitePrior(ScanSeg(), VLBIGaussian(0, 1))
 ```
 
 creates a site prior that is constant across scans and each scan has a unit Normal prior.
diff --git a/src/instrument/priors/refant.jl b/src/instrument/priors/refant.jl
index cbd90c232..445c22f17 100644
--- a/src/instrument/priors/refant.jl
+++ b/src/instrument/priors/refant.jl
@@ -36,7 +36,7 @@ SEFDReference(val::Number) = SEFDReference(val, 0)
 reference_indices(::AbstractArrayConfiguration, st::SiteLookup, ::NoReference) = [], nothing
 function reference_indices(::AbstractArrayConfiguration, st::SiteLookup, p::SingleReference)
     inds = findall(==(p.site), st.sites)
-    return inds, Fill(p.value, length(inds))
+    return inds, fill(p.value, length(inds))
 end
 
 function reference_indices(array::AbstractArrayConfiguration, st::SiteLookup, r::SEFDReference)
@@ -55,5 +55,5 @@ function reference_indices(array::AbstractArrayConfiguration, st::SiteLookup, r:
         _, ind = findmin(map(s -> sefd[s], sites))
         push!(fixedinds, inds[ind])
     end
-    return fixedinds, Fill(r.value, length(fixedinds))
+    return fixedinds, fill(r.value, length(fixedinds))
 end
diff --git a/src/instrument/priors/segmentation.jl b/src/instrument/priors/segmentation.jl
index 3291dfee0..06fa559f2 100644
--- a/src/instrument/priors/segmentation.jl
+++ b/src/instrument/priors/segmentation.jl
@@ -69,7 +69,7 @@ function timestamps(::IntegSeg, array)
     # TODO build in the dt into the data format
     if length(ts) <= 1
         # arbritrarily set the dt to 1
-        dt = 1 / 3600
+        dt = 0.1 / 3600
     else
         dt = minimum(diff(ts))
     end
diff --git a/src/instrument/site_array.jl b/src/instrument/site_array.jl
index 909ed2506..397374f0d 100644
--- a/src/instrument/site_array.jl
+++ b/src/instrument/site_array.jl
@@ -23,6 +23,16 @@ struct SiteArray{T, N, A <: AbstractArray{T, N}, Ti <: AbstractArray{<:Integrati
     sites::Sy
 end
 
+function Adapt.parent_type(::Type{<:Comrade.SiteArray{T, N, A}}) where {T, N, A}
+    return A
+end
+
+function Adapt.adapt(to, x::SiteArray)
+    data = Adapt.adapt(to, parent(x))
+    return SiteArray(data, x.times, x.frequencies, x.sites)
+end
+
+
 # function ChainRulesCore.rrule(::Type{SiteArray}, data::AbstractArray, args...)
 #     s = SiteArray(data, args...)
 #     pd = ProjectTo(data)
@@ -41,15 +51,15 @@ EnzymeRules.inactive(::(typeof(Base.size)), ::SiteArray) = nothing
 Base.parent(a::SiteArray) = getfield(a, :data)
 Base.size(a::SiteArray) = size(parent(a))
 Base.IndexStyle(::Type{<:SiteArray{T, N, A}}) where {T, N, A} = Base.IndexStyle(A)
-Base.@propagate_inbounds Base.getindex(a::SiteArray{T}, i::Integer) where {T} = getindex(parent(a), i)
-Base.@propagate_inbounds Base.getindex(a::SiteArray, I::Vararg{Integer, N}) where {N} = getindex(parent(a), I...)
+Base.@propagate_inbounds Base.getindex(a::SiteArray, i::Int) = rgetindex(parent(a), i)
+Base.@propagate_inbounds Base.getindex(a::SiteArray{T, N}, I::Vararg{Int, N}) where {T, N} = rgetindex(parent(a), I...)
 Base.@propagate_inbounds Base.setindex!(m::SiteArray, v, i::Integer) = setindex!(parent(m), v, i)
 Base.@propagate_inbounds Base.setindex!(m::SiteArray, v, i::Vararg{Integer, N}) where {N} = setindex!(parent(m), v, i...)
-Base.@propagate_inbounds function Base.getindex(m::SiteArray, I...)
-    return SiteArray(getindex(parent(m), I...), getindex(m.times, I...), getindex(m.frequencies, I...), getindex(m.sites, I...))
+Base.@propagate_inbounds function Base.getindex(m::SiteArray, I::AbstractArray...)
+    return SiteArray(rgetindex(parent(m), I...), rgetindex(m.times, I...), rgetindex(m.frequencies, I...), rgetindex(m.sites, I...))
 end
 
-Base.@propagate_inbounds function Base.view(A::SiteArray, I...)
+Base.@propagate_inbounds function Base.view(A::SiteArray, I::AbstractArray...)
     return SiteArray(view(A.data, I...), view(times(A), I...), view(frequencies(A), I...), view(sites(A), I...))
 end
 
diff --git a/src/instrument/utility.jl b/src/instrument/utility.jl
index 795bc2142..ad527745e 100644
--- a/src/instrument/utility.jl
+++ b/src/instrument/utility.jl
@@ -22,8 +22,8 @@ julia> site_tuple(sites, ScanSeg(); AA = FixedSeg(1.0))
 (AA = FixedSeg(1.0), AP = ScanSeg(), LM = ScanSeg(), PV = ScanSeg())
 julia> site_tuple(sites, ScanSeg(); AA = FixedSeg(1.0), PV = TrackSeg())
 (AA = FixedSeg(1.0), AP = ScanSeg(), LM = ScanSeg(), PV = TrackSeg())
-julia> site_tuple(sites, Normal(0.0, 0.1); reference=:AA, LM = Normal(0.0, 1.0))
-(AP = Normal(0.0, 0.1), LM = Normal(0.0, 1.0), PV = Normal(0.0, 0.1))
+julia> site_tuple(sites, VLBIGaussian(0.0, 0.1); reference=:AA, LM = VLBIGaussian(0.0, 1.0))
+(AP = VLBIGaussian(0.0, 0.1), LM = VLBIGaussian(0.0, 1.0), PV = VLBIGaussian(0.0, 0.1))
 ```
 """
 function site_tuple(sites::NTuple{N, Symbol}, default; kwargs...) where {N}
diff --git a/src/mrf_image.jl b/src/mrf_image.jl
index 40ff4c03c..6f772f8aa 100644
--- a/src/mrf_image.jl
+++ b/src/mrf_image.jl
@@ -141,7 +141,7 @@ function corr_image_prior(
     )
     rat = corr_length / step(grid.X)
     cmarkov = ConditionalMarkov(base, grid; order = order)
-    dρ = Dists.truncated(Dists.InverseGamma(1.0, -log(frac_below_beam) * rat); lower, upper)
+    dρ = VLBITruncated(VLBIInverseGamma(1.0, -log(frac_below_beam) * rat); lower, upper)
     cprior = HierarchicalPrior(cmarkov, dρ)
     return cprior
 end
diff --git a/src/observations/array.jl b/src/observations/array.jl
index 0656f1426..5dd55d98c 100644
--- a/src/observations/array.jl
+++ b/src/observations/array.jl
@@ -151,6 +151,14 @@ Base.@kwdef struct EHTArrayConfiguration{A <: AbstractBaselineDatum, F, T, S, D
     datatable::D
 end
 
+function Adapt.adapt_structure(to, c::EHTArrayConfiguration)
+    tarr = adapt(to, c.tarr)
+    scans = adapt(to, c.scans)
+    datatable = adapt(to, c.datatable)
+    return EHTArrayConfiguration(c.bandwidth, tarr, scans, c.mjd, c.ra, c.dec, c.source, c.timetype, datatable)
+end
+
+
 function Base.show(io::IO, config::EHTArrayConfiguration)
     println(io, "EHTArrayConfiguration:")
     println(io, "  source:      ", config.source)
diff --git a/src/observations/obstable.jl b/src/observations/obstable.jl
index 7fd889f41..5432d7f9c 100644
--- a/src/observations/obstable.jl
+++ b/src/observations/obstable.jl
@@ -182,3 +182,10 @@ function Base.view(obs::EHTObservationTable{F}, i::AbstractVector) where {F <: C
     s = @view noise(obs)[i, i]
     return EHTObservationTable{F}(m, s, conf)
 end
+
+function Adapt.adapt_structure(to, x::EHTObservationTable{T}) where {T}
+    m = adapt(to, x.measurement)
+    n = adapt(to, x.noise)
+    c = adapt(to, x.config)
+    return EHTObservationTable{T}(m, n, c)
+end
diff --git a/src/posterior/abstract.jl b/src/posterior/abstract.jl
index 37ecaea5d..6ad49b652 100644
--- a/src/posterior/abstract.jl
+++ b/src/posterior/abstract.jl
@@ -6,6 +6,11 @@ export vlbimodel, logdensityof, dimension, skymodel, instrumentmodel, dataproduc
     VLBIPosterior, logdensityof, loglikelihood, chi2, residuals
 
 
+function prepare_device(any, ex)
+    return any
+end
+
+
 """
     $(TYPEDEF)
 
@@ -58,7 +63,6 @@ EnzymeRules.inactive(::typeof(instrumentmodel), args...) = nothing
 
 function DensityInterface.logdensityof(post::AbstractVLBIPosterior, x)
     pr = logprior(post, x)
-    !isfinite(pr) && return -Inf
     l = loglikelihood(post, x)
     # isnan(l) && throw(ArgumentError("NaN in loglikelihood at $x"))
     return l + pr
diff --git a/src/posterior/vlbiposterior.jl b/src/posterior/vlbiposterior.jl
index f037ed698..391c3ca84 100644
--- a/src/posterior/vlbiposterior.jl
+++ b/src/posterior/vlbiposterior.jl
@@ -20,7 +20,7 @@ end
 """
     VLBIPosterior(skymodel::SkyModel, instumentmodel::InstrumentModel, 
                   dataproducts::EHTObservationTable...; 
-                  admode=set_runtime_activity(Reverse))
+                  admode=set_runtime_activity(set_strong_zero(Reverse)))
 
 Creates a VLBILikelihood using the `skymodel` its related metadata `skymeta`
 and the `instrumentmodel` and its metadata `instumentmeta`. The `model` is a 
@@ -51,12 +51,12 @@ function sky(θ, metadata)
     return m
 end
 
-skyprior = (r = Uniform(μas2rad(10.0), μas2rad(30.0)), a = Uniform(1.0, 10.0))
+skyprior = (r = VLBIUniform(μas2rad(10.0), μas2rad(30.0)), a = VLBIUniform(1.0, 10.0))
 g  = imagepixels(μas2rad(100.0), μas2rad(100.0), 256, 256)
 skym = SkyModel(sky, skyprior, g)
 
-G = SingleStokesGain(x->exp(x.lg + 1im*x.pg))
-intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
+G = SingleStokesGain(x->exp(complex(x.lg, x.pg)))
+intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
             pg = ArrayPrior(IIDSitePrior(ScanSeg(), DiagVonMises(0.0, inv(π^2))))
             )
 intmodel = InstrumentModel(G, intprior, array)
@@ -69,7 +69,7 @@ post = VLBIPosterior(skym, intmodel, dlcamp, dcphase)
         instrumentmodel::AbstractInstrumentModel,
         dataproducts::EHTObservationTable...;
         imgdata = nothing,
-        admode = EnzymeCore.set_runtime_activity(EnzymeCore.Reverse)
+        admode = EnzymeCore.set_runtime_activity(EnzymeCore.set_strong_zero(EnzymeCore.Reverse))
     )
 
 
@@ -93,7 +93,7 @@ end
 
 VLBIPosterior(
     skymodel::AbstractSkyModel, dataproducts::EHTObservationTable...;
-    admode = EnzymeCore.set_runtime_activity(EnzymeCore.Reverse), kwargs...
+    admode = EnzymeCore.set_runtime_activity(EnzymeCore.set_strong_zero(EnzymeCore.Reverse)), kwargs...
 ) =
     VLBIPosterior(skymodel, IdealInstrumentModel(), dataproducts...; admode, kwargs...)
 
diff --git a/src/skymodels/macro.jl b/src/skymodels/macro.jl
index 3742e3264..99a13ca97 100644
--- a/src/skymodels/macro.jl
+++ b/src/skymodels/macro.jl
@@ -182,9 +182,9 @@ this to plug in pre-built compound priors:
 ```julia
 @sky function imager(grid; ftot, mimg, cprior)
     c    ~ cprior                                       # e.g. HierarchicalPrior or corr_image_prior
-    σimg ~ truncated(Normal(0.0, 0.5); lower = 0.0)
-    fg   ~ Uniform(0.0, 1.0)
-    ρs   ~ ntuple(Returns(Uniform(0.01, 10.0)), 3)      # tuple-of-IID priors
+    σimg ~ VLBITruncated(VLBIGaussian(0.0, 0.5); lower = 0.0)
+    fg   ~ VLBIUniform(0.0, 1.0)
+    ρs   ~ ntuple(Returns(VLBIUniform(0.01, 10.0)), 3)  # tuple-of-IID priors
     rast = apply_fluctuations(CenteredLR(), mimg, σimg .* c.params)
     pimg = parent(rast)
     @. pimg = (ftot * (1 - fg)) * pimg
diff --git a/test/Core/bayes.jl b/test/Core/bayes.jl
index c19293cc6..af4656114 100644
--- a/test/Core/bayes.jl
+++ b/test/Core/bayes.jl
@@ -144,11 +144,11 @@ using FiniteDifferences
     @test gz ≈ gfd
 
     R = JonesR()
-    Gp = JonesG(x -> (exp(x.lg + 1im * x.gp), exp(x.lg + 1im * x.gp)))
+    Gp = JonesG(x -> (exp(complex(x.lg, x.gp)), exp(complex(x.lg, x.gp))))
     J = JonesSandwich(Gp, R)
     pr = (
-        lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
-        gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
+        lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0))),
+        gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
     )
     intm_coh = InstrumentModel(J, pr)
     skym = SkyModel(test_skymodel_polarized, test_prior(), imagepixels(μas2rad(150.0), μas2rad(150.0), 256, 256); metadata = (; lp = 0.1))
@@ -174,11 +174,11 @@ end
     post_lc = VLBIPosterior(skym, lcamp)
     post_vis = VLBIPosterior(skym, vis)
 
-    G = SingleStokesGain(x -> exp(x.lg + 1im .* x.gp))
+    G = SingleStokesGain(x -> exp(complex(x.lg, x.gp)))
     intm = InstrumentModel(
         G, (
-            lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
-            gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
+            lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0))),
+            gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
         )
     )
     post_gvis = VLBIPosterior(skym, vis)
@@ -242,8 +242,8 @@ end
     G = SingleStokesGain(x -> exp(x.lg + 1im .* x.gp))
     intm = InstrumentModel(
         G, (
-            lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
-            gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
+            lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0))),
+            gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
         )
     )
 
@@ -275,8 +275,8 @@ end
     G = SingleStokesGain(x -> exp(x.lg + 1im .* x.gp))
     intm = InstrumentModel(
         G, (
-            lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0))),
-            gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
+            lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0))),
+            gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0)), refant = SEFDReference(0.0), phase = true),
         )
     )
 
diff --git a/test/Core/macro.jl b/test/Core/macro.jl
index ef629341e..2b9e9067f 100644
--- a/test/Core/macro.jl
+++ b/test/Core/macro.jl
@@ -9,16 +9,16 @@ using VLBIImagePriors
     @testset "parity with hand-written SkyModel" begin
         # Macro form
         @sky function macro_gauss(grid; flo = μas2rad(1.0), fhi = μas2rad(40.0))
-            σ ~ Uniform(flo, fhi)
-            τ ~ Uniform(0.35, 0.65)
+            σ ~ VLBIUniform(flo, fhi)
+            τ ~ VLBIUniform(0.35, 0.65)
             return stretched(Gaussian(), σ * τ, σ)
         end
 
         # Hand-written equivalent
         manual_f(θ, meta) = stretched(Gaussian(), θ.σ * θ.τ, θ.σ)
         manual_prior = (
-            σ = Uniform(μas2rad(1.0), μas2rad(40.0)),
-            τ = Uniform(0.35, 0.65),
+            σ = VLBIUniform(μas2rad(1.0), μas2rad(40.0)),
+            τ = VLBIUniform(0.35, 0.65),
         )
 
         m_macro = macro_gauss(g)
@@ -41,7 +41,7 @@ using VLBIImagePriors
 
     @testset "metadata flow and grid in metadata" begin
         @sky function with_grid(grid; ftot = 1.0)
-            σ ~ Uniform(μas2rad(1.0), μas2rad(40.0))
+            σ ~ VLBIUniform(μas2rad(1.0), μas2rad(40.0))
             return ftot * stretched(Gaussian(), σ, σ)
         end
         m = with_grid(g; ftot = 2.5)
@@ -55,7 +55,7 @@ using VLBIImagePriors
 
     @testset "kwarg without default is required" begin
         @sky function needs_kw(grid; ftot)
-            σ ~ Uniform(μas2rad(1.0), μas2rad(40.0))
+            σ ~ VLBIUniform(μas2rad(1.0), μas2rad(40.0))
             return ftot * stretched(Gaussian(), σ, σ)
         end
         @test_throws UndefKeywordError needs_kw(g)
@@ -65,11 +65,11 @@ using VLBIImagePriors
 
     @testset "tuple-of-IID priors via ntuple" begin
         @sky function ntuple_prior(grid; n = 3, lo = 0.01, hi = 10.0)
-            ρs ~ ntuple(Returns(Uniform(lo, hi)), n)
+            ρs ~ ntuple(Returns(VLBIUniform(lo, hi)), n)
             return stretched(Gaussian(), ρs[1], ρs[2])
         end
         m = ntuple_prior(g)
-        @test m.prior.ρs isa NTuple{3, <:Uniform}
+        @test m.prior.ρs isa NTuple{3, <:VLBIImagePriors.AffineDistribution}
         x = rand(Comrade.HypercubeTransform.NamedDist(m.prior))
         @test x.ρs isa NTuple{3, Float64}
     end
@@ -80,7 +80,7 @@ using VLBIImagePriors
 
         @sky function hier_imager(grid; cprior, σscale = 0.1)
             c ~ cprior
-            σimg ~ Exponential(σscale)
+            σimg ~ VLBIExponential(σscale)
             rast = to_simplex(CenteredLR(), σimg .* c.params)
             return ContinuousImage(rast, grid, BSplinePulse{3}())
         end
@@ -93,9 +93,9 @@ using VLBIImagePriors
 
     @testset "VLBIPosterior compatibility" begin
         @sky function gauss2(grid)
-            σ ~ Uniform(μas2rad(1.0), μas2rad(40.0))
-            τ ~ Uniform(0.35, 0.65)
-            f1 ~ Uniform(0.8, 1.2)
+            σ ~ VLBIUniform(μas2rad(1.0), μas2rad(40.0))
+            τ ~ VLBIUniform(0.35, 0.65)
+            f1 ~ VLBIUniform(0.8, 1.2)
             return f1 * stretched(Gaussian(), σ * τ, σ)
         end
         m = gauss2(g)
@@ -114,20 +114,20 @@ using VLBIImagePriors
         # Tilde nested inside another expression is rejected
         @test_throws LoadError @eval @sky function nested_tilde(grid)
             if true
-                σ ~ Uniform(0.0, 1.0)
+                σ ~ VLBIUniform(0.0, 1.0)
             end
             return Gaussian()
         end
 
         # Name clash between sampled param and kwarg
         @test_throws LoadError @eval @sky function clash(grid; σ = 1.0)
-            σ ~ Uniform(0.0, 1.0)
+            σ ~ VLBIUniform(0.0, 1.0)
             return stretched(Gaussian(), σ, σ)
         end
 
         # Name clash with positional grid name
         @test_throws LoadError @eval @sky function gridclash(grid)
-            grid ~ Uniform(0.0, 1.0)
+            grid ~ VLBIUniform(0.0, 1.0)
             return Gaussian()
         end
 
@@ -138,14 +138,14 @@ using VLBIImagePriors
 
         # Multiple positional args rejected
         @test_throws LoadError @eval @sky function twopos(grid, extra)
-            σ ~ Uniform(0.0, 1.0)
+            σ ~ VLBIUniform(0.0, 1.0)
             return Gaussian()
         end
 
         # Duplicate tilde names
         @test_throws LoadError @eval @sky function dupes(grid)
-            σ ~ Uniform(0.0, 1.0)
-            σ ~ Uniform(0.0, 2.0)
+            σ ~ VLBIUniform(0.0, 1.0)
+            σ ~ VLBIUniform(0.0, 2.0)
             return Gaussian()
         end
     end
diff --git a/test/Core/models.jl b/test/Core/models.jl
index ef4e0f3a6..f256e78d0 100644
--- a/test/Core/models.jl
+++ b/test/Core/models.jl
@@ -199,8 +199,8 @@ end
 
         G = SingleStokesGain(x -> exp(x.lg + 1im * x.gp))
         intprior = (
-            lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-            gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, inv(π^2))); refant = SEFDReference(0.0)),
+            lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+            gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, inv(π^2))); refant = SEFDReference(0.0)),
         )
 
         intm = InstrumentModel(G, intprior)
@@ -248,8 +248,8 @@ end
 
         G = SingleStokesGain(x -> exp(x.lg + 1im * x.gp))
         intprior = (
-            lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-            gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, inv(π^2))); refant = SEFDReference(0.0)),
+            lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+            gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, inv(π^2))); refant = SEFDReference(0.0)),
         )
 
         intm = InstrumentModel(G, intprior)
@@ -339,20 +339,20 @@ end
 
 
         intprior = (
-            lgR = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-            gpR = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, inv(π^2))); phase = true, refant = SEFDReference(0.0)),
-            lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)), phase = false),
-            gprat = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)), refant = SingleReference(:AA, 0.0)),
-            dRx = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-            dRy = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-            dLx = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-            dLy = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
+            lgR = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+            gpR = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, inv(π^2))); phase = true, refant = SEFDReference(0.0)),
+            lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)), phase = false),
+            gprat = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)), refant = SingleReference(:AA, 0.0)),
+            dRx = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+            dRy = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+            dLx = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+            dLy = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
         )
 
 
         intm = InstrumentModel(J, intprior)
         intm2 = InstrumentModel(J2, intprior)
-        intjg = InstrumentModel(JG, (; lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)))))
+        intjg = InstrumentModel(JG, (; lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)))))
         show(IOBuffer(), MIME"text/plain"(), intm)
 
 
@@ -523,14 +523,14 @@ end
         F = JonesF()
 
         intprior = (
-            lgR = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-            gpR = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, inv(π^2))); phase = true, refant = SEFDReference(0.0)),
-            lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)), phase = false),
-            gprat = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-            dRx = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-            dRy = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-            dLx = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
-            dLy = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.2))),
+            lgR = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+            gpR = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, inv(π^2))); phase = true, refant = SEFDReference(0.0)),
+            lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)), phase = false),
+            gprat = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+            dRx = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+            dRy = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+            dLx = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
+            dLy = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.2))),
         )
 
 
@@ -650,13 +650,13 @@ end
 
         @testset "standard" begin
             intprior = (
-                lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-                gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, inv(π^2))); refant = SEFDReference(0.0)),
+                lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+                gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, inv(π^2))); refant = SEFDReference(0.0)),
             )
 
             intprior_mar = (
-                lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-                gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
+                lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+                gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
             )
 
 
@@ -678,8 +678,8 @@ end
         @testset "markov" begin
 
             intprior = (
-                lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1))),
-                gp = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
+                lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
+                gp = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, inv(π^2))); refant = SEFDReference(0.0), phase = true),
             )
 
             intm = InstrumentModel(G, intprior)
diff --git a/test/test_util.jl b/test/test_util.jl
index e8ac41b2e..2aaedc904 100644
--- a/test/test_util.jl
+++ b/test/test_util.jl
@@ -53,10 +53,10 @@ end
 
 function test_prior2()
     return (
-        f = Uniform(0.8, 1.2),
-        r = Uniform(μas2rad(5.0), μas2rad(30.0)),
-        τ = Uniform(0.0, 0.5),
-        α = Uniform(2.0, 8.0),
+        f = VLBIUniform(0.8, 1.2),
+        r = VLBIUniform(μas2rad(5.0), μas2rad(30.0)),
+        τ = VLBIUniform(0.0, 0.5),
+        α = VLBIUniform(2.0, 8.0),
     )
 end
 
@@ -70,7 +70,7 @@ end
 
 function testimg_prior(grid)
     cprior = corr_image_prior(grid, μas2rad(20.0))
-    prior = (c = cprior, σimg = Exponential(0.1))
+    prior = (c = cprior, σimg = VLBIExponential(0.1))
     return prior
 end
 
@@ -85,22 +85,22 @@ end
 
 function testimg_add_prior(grid)
     cprior = corr_image_prior(grid, μas2rad(20.0))
-    prior = (c = cprior, σimg = Exponential(0.1), fg = Uniform(0.0, 1.0))
+    prior = (c = cprior, σimg = VLBIExponential(0.1), fg = VLBIUniform(0.0, 1.0))
     return prior
 end
 
 
 function test_prior()
     return (
-        f1 = Uniform(0.8, 1.2),
-        σ1 = Uniform(μas2rad(1.0), μas2rad(40.0)),
-        τ1 = Uniform(0.35, 0.65),
-        ξ1 = Uniform(-π / 2, π / 2),
-        f2 = Uniform(0.3, 0.7),
-        σ2 = Uniform(μas2rad(1.0), μas2rad(40.0)),
-        τ2 = Uniform(0.35, 0.65),
-        ξ2 = Uniform(-π / 2, π / 2),
-        x = Uniform(-μas2rad(40.0), μas2rad(40.0)),
-        y = Uniform(-μas2rad(40.0), μas2rad(40.0)),
+        f1 = VLBIUniform(0.8, 1.2),
+        σ1 = VLBIUniform(μas2rad(1.0), μas2rad(40.0)),
+        τ1 = VLBIUniform(0.35, 0.65),
+        ξ1 = VLBIUniform(-π / 2, π / 2),
+        f2 = VLBIUniform(0.3, 0.7),
+        σ2 = VLBIUniform(μas2rad(1.0), μas2rad(40.0)),
+        τ2 = VLBIUniform(0.35, 0.65),
+        ξ2 = VLBIUniform(-π / 2, π / 2),
+        x = VLBIUniform(-μas2rad(40.0), μas2rad(40.0)),
+        y = VLBIUniform(-μas2rad(40.0), μas2rad(40.0)),
     )
 end