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::ComradeBase.AbstractDims)

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.

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. This is helpful for Zygote autodiff.

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

A image prior based off of the first-order zero mean Gaussian Markov random field. This prior is similar to the combination of total squared variation TSV and L₂ norm, and is given by

(π²Σ)⁻¹ TSV(I) + (ρ^2π²Σ)L₂(I) + lognorm(ρ, Σ)

where ρ and Σ are the correlation length and variance of the random field and lognorm(ρ,Σ) is the log-normalization of the random field. This normalization is needed to jointly infer I and the hyperparameters ρ, Σ.

Fields

• ρ

  The correlation length of the random field.

• Σ

  The variance of the random field

• cache

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

Example

julia> ρ, Σ = 2.0, 1.0
julia> d = GaussMarkovRandomField(ρ, Σ, (32, 32))
julia> cache = MarkovRandomFieldCache(Float64, (32, 32)) # now instead construct the cache
julia> d2 = GaussMarkovRandomField(ρ, Σ, cache)
julia> invcov(d) ≈ invcov(d2)
true

GaussMarkovRandomField(ρ, Σ, img::AbstractArray)

Constructs a first order zero-mean Gaussian Markov random field with dimensions size(img), correlation ρ and diagonal covariance Σ.

GaussMarkovRandomField(ρ, Σ, cache::MarkovRandomFieldCache)

Constructs a first order Gaussian Markov random field using the precomputed cache cache. ρ and diagonal covariance Σ and the precomputed MarkovRandomFieldCache cache. If mean is not included then it is assume the mean is identically zero.

GaussMarkovRandomField(ρ, Σ, dims)

Constructs a first order zero-mean Gaussian Markov random field with dimensions dims, correlation ρ and diagonal covariance Σ.

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.

struct MarkovRandomFieldCache{A, TD, M}

A cache for a Markov random field.

Fields

• Λ

  Intrinsic Gaussian Random Field pseudo-precison matrix. This is similar to the TSV regularizer

• D

  Gaussian Random Field diagonal precision matrix. This is similar to the L2-regularizer

• λQ

  The eigenvalues of the Λ matrix which is needed to compute the log-normalization constant of the GMRF.

Examples

julia> mean = zeros(2,2) # define a zero mean
julia> cache = GRMFCache(mean)
julia> prior_map(x) = GaussMarkovRandomField(mean, x[1], x[2], cache)
julia> d = HierarchicalPrior(prior_map, product_distribution([Uniform(-5.0,5.0), Uniform(0.0, 10.0)]))
julia> x0 = rand(d)

MarkovRandomFieldCache(grid::AbstractDims)

Create a GMRF cache out of a Comrade model as the mean image.

Arguments

• m: Comrade model used for the mean image.
• grid: The grid of points you want to create the model image on.

Keyword arguments

• transform = identity: A transform to apply to the image when creating the mean image. See the examples

Example

julia> m = MarkovRandomFieldCache(imagepixels(10.0, 10.0, 64, 64))

MarkovRandomFieldCache(mean::AbstractMatrix)

Contructs the MarkovRandomFieldCache from the mean image mean. This is useful for hierarchical priors where you change the hyperparameters of the GaussMarkovRandomField, ρ and Σ.

NamedDist(d::NamedTuple{N})
NamedDist(;dists...)

A Distribution with names N. This is useful to construct a set of random variables with a set of names attached to them.

julia> d = NamedDist((a=Normal(), b = Uniform(), c = MvNormal(randn(2), rand(2))))
julia> x = rand(d)
(a = 0.13789342, b = 0.2347895, c = [2.023984392, -3.09023840923])
julia> logpdf(d, x)

Note that NamedDist values passed to NamedDist can also be abstract collections of distributions as well

julia> d = NamedDist(a = Normal(),
                     b = MvNormal(ones(2)),
                     c = (Uniform(), InverseGamma())
                     d = (a = Normal(), Beta)
                    )

How this is done internally is considered an implementation detail and is not part of the public interface.

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 variance of the random field

• ν

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

• cache

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

Examples

julia> ρ, Σ = 2.0, 1.0
julia> d = TDistMarkovRandomField(ρ, Σ, (32, 32))
julia> cache = MarkovRandomFieldCache(Float64, (32, 32)) # now instead construct the cache
julia> d2 = TDistMarkovRandomField(ρ, Σ, cache)
julia> invcov(d) ≈ invcov(d2)
true

TDistMarkovRandomField(ρ, Σ, img::AbstractArray)

Constructs a first order TDist Markov random field with mean image mean and correlation ρ and diagonal covariance Σ.

TDistMarkovRandomField(ρ, Σ, cache::MarkovRandomFieldCache)

Constructs a first order TDist Markov random field with zero mean ,correlation ρ, diagonal covariance Σ, and the precomputed MarkovRandomFieldCache cache.

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.

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 julia> y = randn(100) julia> y .= y./sum(y) julia> toreal(CenteredLR(), y) julia> toreal(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)