GitHub

Modeling with non-analytic Fourier transforms

using VLBISkyModels


using CairoMakie
  Activating project at `~/work/VLBISkyModels.jl/VLBISkyModels.jl/examples`

While most of the models implemented in VLBISkyModels have an analytic Fourier transform this is not required. In this notebook we will describe how a user can do Bayesian model fitting with a non-analytic model.

The ExtendedRing model is an example of a non-analytic model. The image structure is given by

\[I(r) = \frac{\beta^\alpha}{2\pi \Gamma(\alpha)} r^{-\alpha-2}e^{-\beta/r}\]

This can be created as follows

m = ExtendedRing(8.0)
ExtendedRing{Float64}(8.0)

The argument is \alpha in the above equation. beta is given by $(1+\alpha)$.

This model does not have a simple analytic Fourier transform, e.g.

VLBISkyModels.visanalytic(ExtendedRing)
ComradeBase.NotAnalytic()

Therefore, to find the Fourier transform of the image we need to revert to numerical methods. For this notebook we will use the fast Fourier transform or FFT. Specifically we will use FFTW. To compute a numerical Fourier transform we first need to specify the image domain

gim = imagepixels(10.0, 10.0, 256, 256)
RectiGrid(
executor: Serial()
Dimensions: 
(↓ X Sampled{Float64} LinRange{Float64}(-4.98046875, 4.98046875, 256) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.98046875, 4.98046875, 256) ForwardOrdered Regular Points)
)

Second we need to specify the list of points in the uv domain we are interested in. Since VLBI tends sparsely sample the UV plan we provide a specific type for this type called UnstructuredDomain that can be used to specify the UV points,

u = randn(1000) / 2
v = randn(1000) / 2
guv = UnstructuredDomain((U=u, V=v))
UnstructredDomain(
executor: Serial()
Dimensions: 
1000-element StructArray(::Vector{Float64}, ::Vector{Float64}) with eltype @NamedTuple{U::Float64, V::Float64}:
 (U = -0.5971716336757743, V = -0.6933076742949483)
 (U = -0.3860173366059538, V = -0.7032930330645235)
 (U = -0.0909533407647639, V = 0.060411579279338934)
 (U = 0.07826362521753258, V = -0.1654473740348924)
 (U = 0.07246288964690037, V = 0.2628028098957007)
 (U = 0.5944038377533504, V = -0.36686367248017443)
 (U = 0.781273696373959, V = 0.07866555689089193)
 (U = -0.6324281115244417, V = 0.11926581668102483)
 (U = -0.49438708895956124, V = 0.5971168020325188)
 (U = 0.19323382364928224, V = -0.5099122171095087)
 ⋮
 (U = 0.018790236438149064, V = -0.6515282666333648)
 (U = 0.3807717462795669, V = -0.5213912470130875)
 (U = -0.1172289817307216, V = -0.3812093646088345)
 (U = -0.8475339335144163, V = 1.073354239254877)
 (U = 0.8443191705644753, V = 0.35476822221018156)
 (U = -0.4378612330761853, V = 0.00300710062448588)
 (U = 0.28626116944855773, V = 1.195471611735109)
 (U = -0.22202237131945426, V = 0.5812799229302226)
 (U = -0.09458077021824099, V = 0.03368616690694055)
)

We can now create a dual domain that contains both the image and the UV points and the specific Fourier transform algorithm we want to use. The options for algorithms are:

  1. NFFTAlg: Uses the Non-Uniform Fast Fourier Transform. Fast and accurate, this is the recommended algorithm for most cases.
  2. DFTAlg: Uses the Discrete Time Non-Uniform Fourier transform. Slow but exact, so only use for small image grids.
  3. FFTAlg: Moderately fast and moderately accurate. Typically only used internally for testing.

For this example we will use NFFTAlg since it is the recommended algorithm for most cases.

gfour = FourierDualDomain(gim, guv, NFFTAlg())
FourierDualDomain(
Algorithm: NFFTAlg{Float64, AbstractNFFTs.PrecomputeFlags, UInt32}(1, 1.0e-9, AbstractNFFTs.TENSOR, 0x00000000)
Image Domain: RectiGrid(
executor: Serial()
Dimensions: 
(↓ X Sampled{Float64} LinRange{Float64}(-4.98046875, 4.98046875, 256) ForwardOrdered Regular Points,
→ Y Sampled{Float64} LinRange{Float64}(-4.98046875, 4.98046875, 256) ForwardOrdered Regular Points)
)
Visibility Domain: UnstructredDomain(
executor: Serial()
Dimensions: 
1000-element StructArray(::Vector{Float64}, ::Vector{Float64}) with eltype @NamedTuple{U::Float64, V::Float64}:
 (U = -0.5971716336757743, V = -0.6933076742949483)
 (U = -0.3860173366059538, V = -0.7032930330645235)
 (U = -0.0909533407647639, V = 0.060411579279338934)
 (U = 0.07826362521753258, V = -0.1654473740348924)
 (U = 0.07246288964690037, V = 0.2628028098957007)
 (U = 0.5944038377533504, V = -0.36686367248017443)
 (U = 0.781273696373959, V = 0.07866555689089193)
 (U = -0.6324281115244417, V = 0.11926581668102483)
 (U = -0.49438708895956124, V = 0.5971168020325188)
 (U = 0.19323382364928224, V = -0.5099122171095087)
 ⋮
 (U = 0.018790236438149064, V = -0.6515282666333648)
 (U = 0.3807717462795669, V = -0.5213912470130875)
 (U = -0.1172289817307216, V = -0.3812093646088345)
 (U = -0.8475339335144163, V = 1.073354239254877)
 (U = 0.8443191705644753, V = 0.35476822221018156)
 (U = -0.4378612330761853, V = 0.00300710062448588)
 (U = 0.28626116944855773, V = 1.195471611735109)
 (U = -0.22202237131945426, V = 0.5812799229302226)
 (U = -0.09458077021824099, V = 0.03368616690694055)
)

Given this FourierDualDomain everything now works as before with analytic models. For example we can compute the intensitymap of the model as

img = intensitymap(m, gfour)
imageviz(img)
Example block output

and the visibility map using

vis = visibilitymap(m, gfour)
fig, ax = scatter(hypot.(vis.U, vis.V), real.(vis); label="Real")
scatter!(ax, hypot.(vis.U, vis.V), imag.(vis); label="Imag")
axislegend(ax)
ax.xlabel = "uv-dist"
ax.ylabel = "Flux"
fig
Example block output

Additionally, you can also modify the models and add them in complete generality. For example

mmod = modify(m, Shift(2.0, 2.0)) + Gaussian()
mimg = intensitymap(mmod, gfour)
mvis = visibilitymap(mmod, gfour)
1000-element UnstructuredMap{ComplexF64, Vector{ComplexF64}, UnstructuredDomain{StructArrays.StructVector{@NamedTuple{U::Float64, V::Float64}, @NamedTuple{U::Vector{Float64}, V::Vector{Float64}}, Int64}, Serial, ComradeBase.NoHeader}}:
 -0.00020851741769414234 + 0.0001210878231282749im
   -0.019649562104313047 + 0.04089256991426133im
      1.5294852902970941 - 0.2983931450152608im
      0.7532759389186023 - 0.4633338551930699im
     0.14049340321813517 - 0.1663080942763054im
     0.09892660255925315 - 0.028672666475414138im
    0.010163026461097917 + 0.052579770213922614im
    -0.13218358276986972 + 0.022316308572378773im
   -0.016401349102311585 - 0.057065750579294555im
     0.11739612767047793 - 0.1269268484193734im
                         ⋮
      0.0133966122373225 + 0.12929875875648184im
    0.026141068439281625 + 0.13025453872319034im
    -0.05504334382895032 - 0.001232682621723308im
   -0.016837252814948436 + 0.005281756759685787im
 -0.00041053785018367123 + 0.00027893756497882955im
    -0.07288745278016265 - 0.10186839130785598im
    0.026653201744456834 - 0.006186107047573407im
    0.029072357186724105 + 0.14204391914034425im
       1.415181972630129 - 0.5712859773337271im

Plotting everything gives

fig = Figure(; size=(800, 400))
ax1 = Axis(fig[1, 1]; xreversed=true, xlabel="RA (radians)", ylabel="Dec (radians)",
           aspect=1)
ax2 = Axis(fig[1, 2]; xlabel="uv-dist", ylabel="Amplitude")
image!(ax1, mimg; colormap=:afmhot)
scatter!(ax2, hypot.(mvis.U, mvis.V), abs.(mvis); label="Real")
fig
Example block output

This page was generated using Literate.jl.