VLBIImagePriors

AdditiveLR <: LogRatioTransform

Defines the additive log-ratio transform. The clr transformation moves from the simplex Sⁿ → $R^{n-1}$ and is given by

alr(x) = [log(x₁/xₙ) ... log(xₙ/xₙ)],

where g(x) = (∏xᵢ)ⁿ⁻¹ is the geometric mean. The inverse transformation is given by

alr⁻¹(x) = exp.(x)./(1 + sum(x[1:n-1])).

AngleTransform

A transformation that moves two vector x and y to an angle θ. Note that is x and y are normally distributed then the resulting distribution in θ is uniform on the circle.

CenterImage(X, Y)
CenterImage(img::IntensityMap)
CenterImage(grid::RectiGrid)

Constructs a projections operator that will take an arbritrary image and return a transformed version whose centroid is at the origin.

Arguments

• X, Y: the iterators for the image pixel locations
• img: A IntensityMap whose grid is the same as the image you want to center
• grid: The grid that the image is defined on.

Example

julia> using ComradeBase: centroid, imagepixels
julia> grid = imagepixels(μas2rad(100.0), μas2rad(100.0), 256, 256)
julia> K = CenterImage(grid)
julia> img = IntensityMap(rand(256, 256), grid)
julia> centroid(img)
(1.34534e-10, 5.23423e-11)
julia> cimg = K(img)
julia> centroid(cimg)
(1.231e-14, -2.43e-14)

Note

This center image works using a linear projection operator. This means that is does not necessarily preserve the image total flux and the positive definiteness of the image. In practise we find that the deviation from the original flux, and the amount of negative flux is small.

Explanation

The centroid is constructed by

X ⋅ I = 0
Y ⋅ I = 0

where I is the flattened image of length N, and X and Y are the image pixel locations. This can be simplified into the matrix equation

XY ⋅ I = 0

where XY is the 2×N matrix whose first for is given by X and second is Y. The space of centered images is then given by the null space of XY. Given this we can form a matrix C which is the kernel matrix of the nullspace whose columns are the orthogonal basis vectors of the null space of XY. Using this we can construct a centered image projection opertor by

K = CC'.

Our centered image is then given by I₀ = KI.

CenteredLR <: LogRatioTransform

Defines the centered log-ratio transform. The clr transformation moves from the simplex Sⁿ → Rⁿ and is given by

clr(x) = [log(x₁/g(x)) ... log(xₙ/g(x))]

where g(x) = (∏xᵢ)ⁿ⁻¹ is the geometric mean. The inverse transformation is given by the softmax function and is only defined on a subset of the domain otherwise it is not injective

clr⁻¹(x) = exp.(x)./sum(exp, x).

Notes

As mentioned above this transformation is bijective on the entire codomain of the function. However, unlike the additive log-ratio transform it does not treat any pixel as being special.

CenteredRegularizer(x, y, σ, p)

Regularizes a general image prior p such that the center of light is close the the origin of the imag. After regularization the log density of the prior is modified to

\[ \log p(I) \to \log p(I) - \frac{(x_C^2 + y_C^2)^2}{2\sigma\^2} N_x N_y\]

where N_x and N_y are the number of pixels in the x and y direction of the image, and $x_C, y_C$ are the center of light of the image I.

ConditionalMarkov(D, args...)

Creates a Conditional Markov measure, that behaves as a Julia functional. The functional returns a probability measure defined by the arguments passed to the functional.

Arguments

• D: The <: MarkovRandomField that defines the underlying measure
• args: Additional arguments used to construct the Markov random field cache. See MarkovRandomFieldGraph for more information.

Example

julia> grid = imagepixels(10.0, 10.0, 64, 64)
julia> ℓ = ConditionalMarkov(GaussMarkovRandomField, grid)
julia> d = ℓ(16) # This is now a distribution
julia> rand(d)

DiagonalVonMises(μ::Real, κ::Real)
DiagonalVonMises(μ::AbstractVector{<:Real}, κ::AbstractVector{<:Real})

Constructs a Von Mises distribution, with mean μ and concentraion parameter κ. If μ and κ are vectors then this constructs a independent multivariate Von Mises distribution.

Notes

This is a custom implementation since the version in Distributions.jl has certain properties that do not play well (having an support only between [-π+μ, π+μ]) with usual VLBI problems. Additionally this distribution has a special overloaded product_distribution method so that concatenating multiple DiagonalVonMises together preserves the type.

Alias for `ExpMarkovRandomField`

struct ExpMarkovRandomField{T<:Number, C} <: VLBIImagePriors.MarkovRandomField

A image prior based off of the zero mean unit variance Exponential Markov random field. The order of the Markov random field is specified

Fields

• ρ: The correlation length of the random field.

• graph: The Markov Random Field graph cache used to define the specific Markov random field class used.

Example

julia> ρ = 10.0
julia> d = ExpMarkovRandomField(ρ, (32, 32))
julia> cache = MarkovRandomFieldGraph(Float64, (32, 32)) # now instead construct the cache
julia> d2 = ExpMarkovRandomField(ρ, cache)
julia> scalematrix(d) ≈ scalematrix(d2)
true

ExpMarkovRandomField(ρ, img::AbstractArray; order::Integer=1)

Constructs a orderᵗʰ order Exponential Markov random field with dimensions size(img), correlation ρ and unit covariance.

The order parameter controls the smoothness of the field with higher orders being smoother. We recommend sticking with either order=1,2. For more information about the impact of the order see MarkovRandomFieldGraph.

ExpMarkovRandomField(ρ, cache::MarkovRandomFieldGraph)

Constructs a first order zero-mean and unit variance Exponential Markov random field using the precomputed cache cache.

ExpMarkovRandomField(ρ, dims; order=1)

Constructs a first order zero-mean unit variance Exponential Markov random field with dimensions dims, correlation ρ.

The order parameter controls the smoothness of the field with higher orders being smoother. We recommend sticking with either order=1,2. For more information about the impact of the order see MarkovRandomFieldGraph.

Alias for `GaussMarkovRandomField`

GaussMarkovRandomField(ρ, img::AbstractArray; order::Integer=1)

Constructs a orderᵗʰ order Gaussian Markov random field with dimensions size(img), correlation ρ and unit covariance.

The order parameter controls the smoothness of the field with higher orders being smoother. We recommend sticking with either order=1,2. Noting that order=1 is equivalent to the usual TSV and L₂ regularization from RML imaging. For more information about the impact of the order see MarkovRandomFieldGraph.

GaussMarkovRandomField(ρ, dims; order=1)
GaussMarkovRandomField(ρ, g::AbstractRectiGrid; order=1)

Constructs a orderᵗʰ order Gaussian Markov random field with dimensions size(img), correlation ρ and unit covariance.

The order parameter controls the smoothness of the field with higher orders being smoother. We recommend sticking with either order=1,2. Noting that order=1 is equivalent to the usual TSV and L₂ regularization from RML imaging. For more information about the impact of the order see MarkovRandomFieldGraph.

ImageDirichlet(α::AbstractMatrix)
ImageDirichlet(α::Real, ny, nx)

A Dirichlet distribution defined on a matrix. Samples from this produce matrices whose elements sum to unity. This is a useful image prior when you want to separately constrain the flux. The α parameter defines the usual Dirichlet concentration amount.

Notes

Much of this code was taken from Distributions.jl and it's Dirichlet distribution. However, some changes were made to make it faster. Additionally, we use define a custom rrule to speed up derivatives.

ImageSimplex(ny,nx)

This defines a transformation from ℝⁿ⁻¹ to the n probability simplex defined on an matrix with dimension ny×nx. This is a more natural transformation for rasterized images, which are most naturally represented as a matrix.

Notes

Much of this code was inspired by TransformVariables. However, we have specified custom rrules using Enzyme as a backend. This allowed the simplex transform to be used with Zygote and we achieved an order of magnitude speedup when computing the pullback of the simplex transform.

ImageSphericalUniform(nx, ny)

Construct a distribution where each image pixel is a 3-sphere uniform variable. This is useful for polarization where the stokes parameters are parameterized on the 3-sphere.

Currently we use a struct of vectors memory layout. That is the image is described by three matrices (X,Y,Z) grouped together as a tuple, where each matrix is one direction on the sphere, and we require norm((X,Y,Z)) == 1.

ImageUniform(a::Real, b::Real, nx, ny)

A uniform distribution in image pixels where a/b are the lower/upper bound for the interval. This then concatenates ny×nx uniform distributions together.

abstract type MarkovRandomField <: Distributions.Distribution{Distributions.Matrixvariate, Distributions.Continuous}

An abstract type for a MarkovRandomField. We assume that the distribution is of the form

p(x | ρ) = N(detQ(ρ)) f(xᵀQ(ρ)x),

where f is a function and N is the normalization of the distribution, and ρ is the correlation parameter.

To implement the informal interface e.g., MyRF <: MarkovRandomField, the user must implement

• lognorm(d::MyRF): Which computes the log of the normalization constant N
• unnormed_logpdf(d::MyRF, x::AbstractMatrix): Which computes f(xᵀQx)
• Distributions._rand!(rng::AbstractRNG, d::MyRF, x::AbstractMatrix): To enable sampling from the prior

Additionally, there are a number of auto-generated function that can be overwritten:

• graph(d::MyRF): Which returns the graph structure of the Markov Random Field. The default returns the property d.graph.
• corrparam(d::MyRF): Which returns the correlation parameter ρ of the Markov Random Field. The default returns the property d.ρ.
• Base.size(d::MyRF): Which returns the size of the distribution. This defaults to the size of the graph cache.
• scalematrix(d::MyRF): Which computes the scale matrix Q, of the random field. The default is to forward to the scalematrix(graph(d), corrparm(d)).
• (c::ConditionalMarkov{<:MyRF})(ρ): To map from a correlation to the distribution
• HypercubeTransform.asflat(d::MyRF): To map from parameter space to a flattened version. The default is TransformVariables.as(Matrix, size(d)...)
• Distributions.insupport(d::MyRF, x::AbstractMatrix) which checks if x is in the support of d. The default is to always return true.
• LinearAlgebra.logdet(d::MyRF) which computes the log determinant of Q. This defaults to logdet(graph(d), corrparam(d)).

Finally additional optional methods are:

• Distributions.mean(d::MyRF): Which computes the mean of the process
• Distributions.cov(d::MyRF): Which computes the covariance matrix of the process.
• Distributions.invcov(d::MyRF): Computes the precision matrix of the random field

For an example implementation see e.g., GaussMarkovRandomField

MarkovRandomFieldGraph(grid::AbstractRectiGrid; order=1)
MarkovRandomFieldGraph(img::AbstractMatrix; order=1)

Create a order Markov random field using the grid or image as its dimension.

The order keyword argument specifies the order of the Markov random process and is generally given by the precision matrix

Qₙ = τ(κI + G)ⁿ

where n = order, I is the identity matrix, G is specified by the first order stencil

.  -1  .
-1  4  -1
.   4  .

κ is the Markov process hyper-parameters. For n=1 κ is related to the correlation length ρ of the random field by

ρ = 1/κ

while for n>1 it is given by

ρ = √(8(n-1))/κ

Note that κ isn't set in the MarkovRandomFieldGraph, but rather once the noise process is set, i.e. one of the subtypes of MarkovRandomField.

Finally τ is chosen so that the marginal variance of the random field is unity. For n=1

τ = 1

for n=2

τ = 4πκ²

and for n>2 we have

τ = (N+1)4π κ²⁽ⁿ⁺¹⁾

Example

julia> m = MarkovRandomFieldGraph(imagepixels(10.0, 10.0, 64, 64))
julia> ρ = 10 # correlation length
julia> d = GaussMarkovRandomField(ρ, m) # build the Gaussian Markov random field

SphericalUnitVector{N}()

A transformation from a set of N+1 vectors to the N sphere. The set of N+1 vectors are inherently assumed to be N+1 a distributed according to a unit multivariate Normal distribution.

Notes

For more information about this transformation see the Stan manual. In the future this may be depricated when is merged.

struct TDistMarkovRandomField{T<:Number, C} <: VLBIImagePriors.MarkovRandomField

A image prior based off of the first-order Multivariate T distribution Markov random field.

Fields

• ρ: The correlation length of the random field.

• ν: The student T "degrees of freedom parameter which ≥ 1 for a proper prior

• graph: The Markov Random Field graph cache used to define the specific Markov random field class used.

Examples

julia> ρ, ν = 16.0, 1.0
julia> d = TDistMarkovRandomField(ρ, ν, (32, 32))
julia> cache = MarkovRandomFieldGraph(Float64, (32, 32)) # now instead construct the cache
julia> d2 = TDistMarkovRandomField(ρ, ν, cache)
julia> invcov(d) ≈ invcov(d2)
true

TDistMarkovRandomField(ρ, ν, img::AbstractArray; order=1)

Constructs a first order TDist Markov random field with zero median dimensions size(img), correlation ρ and degrees of freedom ν.

Note ν ≥ 1 to be a well-defined probability distribution.

The order parameter controls the smoothness of the field with higher orders being smoother. We recommend sticking with either order=1,2. For more information about the impact of the order see MarkovRandomFieldGraph.

TDistMarkovRandomField(ρ, ν, cache::MarkovRandomFieldGraph)

Constructs a first order TDist Markov random field with zero mean ,correlation ρ, degrees of freedom ν, and the precomputed MarkovRandomFieldGraph cache.

TDistMarkovRandomField(ρ, ν, dims)

Constructs a first order TDist Markov random field with zero mean ,correlation ρ, degrees of freedom ν, with dimension dims.

The order parameter controls the smoothness of the field with higher orders being smoother. We recommend sticking with either order=1,2. For more information about the impact of the order see MarkovRandomFieldGraph.

Alias for `TDistMarkovRandomField`

WrappedUniform(period)

Constructs a potentially multivariate uniform distribution that is wrapped a given period. That is

d = WrappedUniform(period)
logpdf(d, x) ≈ logpdf(d, x+period)

for any x.

If period is a vector this creates a multivariate independent wrapped uniform distribution.

alr!(x, y)

Compute the inverse alr transform. That is x lives in ℜⁿ and y, lives in Δⁿ

alrinv!(x, y)

Computes the additive logit transform inplace. This converts from ℜⁿ → Δⁿ where Δⁿ is the n-simplex

Notes

This function is mainly to transform the GaussMarkovRandomField to live on the simplex. In order to preserve the nice properties of the GRMF namely the log det we only use y[begin:end-1] elements and the last one is not included in the transform. This shouldn't be a problem since the additional parameter is just a dummy in that case and the Gaussian prior should ensure it is easy to sample from.

clr!(x, y)

Compute the inverse alr transform. That is x lives in ℜⁿ and y, lives in Δⁿ

clrinv!(x, y)

Computes the additive logit transform inplace. This converts from ℜⁿ → Δⁿ where Δⁿ is the n-simplex

Notes

This function is mainly to transform the GaussMarkovRandomField to live on the simplex.

corrparam(m::MarkovRandomField)

Returns the correlation parameter of the Markov Random field m. For details about the correlation parmeter see MarkovRandomFieldGraph.

graph(m::MarkovRandomField)

Returns the graph or graph cache of the Markov Random field m.

matern([T=Float64], dims::Dims{2})
matern([T=Float64], dims::Int...)

Creates an approximate Matern Gaussian process that approximates the Matern process on a regular grid which cyclic boundary conditions. This function returns a tuple of two objects

• A functor f of type StationaryMatern that iid-Normal noise to a draw from the Matern process. The functor call arguments are f(s, ρ, ν) where s is the random white Gaussian noise with dimension dims, ρ is the correlation length, and ν is Matern smoothness parameter
• The a set of prod(dims) standard Normal distributions that can serve as the noise generator for the process.

Example

julia> transform, dstd = matern((32, 32))
julia> draw_matern = transform(rand(dstd), 10.0, 2.0)
julia> draw_matern_aniso = transform(rand(dstd), (10.0, 5.0), π/4 2.0) # anisotropic Matern
julia> ones(32, 32) .+ 5.* draw_matern # change the mean and variance of the field

matern(img::AbstractMatrix)

Creates an approximate Matern Gaussian process with dimension size(img)

scalematrix(m::MarkovRandomField)

Return the scale matrix for the Markov Random field. For a Gaussian Markov random field this corresponds to the precision matrix of the Gaussian field.

For other random processes this is the metric of the inner product, i.e. Q in

xᵀQx

which is the distance from the origin to x using the metric Q.

to_real(t::LogRatioTransform, y)

Transform the value u that lives on the simplex to a value in the real vector space. See subtypes(LogRatioTransform) for a list of possible transformations.

The inverse of this transform is given by to_simplex(t, x).

Example

julia> y = randn(100)
julia> y .= y./sum(y)
julia> x = similar(y)
julia> to_real!(CenteredLR(), x, y)
julia> xar = similar(y, 99)
julia> to_real!(AdditiveLR(), xar, y)

to_real(t::LogRatioTransform, u)

Transform the value u that lives on the simplex to a value x in the real vector space. See subtypes(LogRatioTransform) for a list of possible transformations.

The inverse of this transform is given by to_simplex!(t, y, x).

Example

julia> y = randn(100)
julia> y .= y./sum(y)
julia> to_real(CenteredLR(), y)
julia> to_real(AdditiveLR(), y)

to_simplex!(t::LogRatioTransform, y, x)

Transform the vector x assumed to be a real valued array to the simplex using the log-ratio transform t and stores the value in y.

The inverse of this transform is given by to_real!(t, x, y) where y is a vector that sums to unity, i.e. it lives on the simplex.

Example

julia> x = randn(100)
julia> to_simplex(CenteredLR(), x)
julia> to_simplex(AdditiveLR(), x)

to_simplex(t::LogRatioTransform, x)

Transform the vector x assumed to be a real valued array to the simplex using the log-ratio transform t. See subtypes(LogRatioTransform) for a list of possible transformations.

The inverse of this transform is given by to_real(t, y) where y is a vector that sums to unity, i.e. it lives on the simplex.

Example

julia> x = randn(100)
julia> to_simplex(CenteredLR(), x)
julia> to_simplex(AdditiveLR(), x)