Implementation guidelines for Comrade compatibility

[ptiede]

This is some pseudo-code that could make StationaryRandomFields.jl compatible with Comrade.

Comrade expects its priors to be a Distributions.jl object. Since we want a general multivariate distribution that means we will need to use the ArrayLikeVariate{N} stuff

import Distributions as Dist
struct PowerSpectrumDist{N, P<:AbstractPowerSpectrumModel{N},...} <: Dist.Distribution{Dist.ArrayLikeVariate{N}, Dist.Continuous}
      ps::P
     others?...
end

# To make compatible with Distributions we will need to implement

# Evaluating the logpdf
function Dists._logpdf(d::PowerSpectrumDist{N}, x::AbstractArray{<:Real, N}) where {N}
    # hint for this we likely will need an ifft to evaluate the results
...
end

# generating random numbers in-place
function Dists._rand!(rng::AbstractRNG, d::PowerSpectrumDist{N}, A::AbstractArray{<:Real, N}) where {N}
...
end

# The "size of the distribution" essentially the size of the input `x` in logpdf
Base.size(d::PowerSpectrumDist) = ...

# sampler (likely  optional)
Dists.sampler(d::PowerSpectrumDist) = ...


# To make this work with Comrade we also need to specify how to transform this. 
# Since we are dealing with Power spectrum models we probably want to transform
# to the standardized representation using an FFT let me know if you need help with this
function HypercubeTransform.asflat(d::PowerSpectrumDist)

end

For an idea for how this stuff would work you can check https://github.com/ptiede/VLBIImagePriors.jl/blob/main/src/matern.jl
which is a quick and dirty implementation of the Matern process, and isn't nearly as flexible as this package.

[ptiede]

@rohandahale

[kazuakiyama]

@ptiede @anilipour @rohandahale

I had a block of time this weekend, to take a look at this, and I made a substantial overhaul of the package. Please check the branch kazuakiyama_reinfrastructure

Now everything is compatible with Distributions.jl and stationary random fields handled in the library is completely generalized. Everything is now under StationaryRandomField type. Following @ptiede 's comment, now this package handles general fields not limited to Gaussian. Indeed, I implement a Laplace field which may be used by @anilipour for Bayeisan L1 RML.

I still plan to revise the package by myself as also there are still bugs and documentation is not updated. At this moment, I would like to get your feedback.

#---------------------------------------------
Here is the hierarchical structures.

Upstream End
AbstractStationaryRandomField{N} <: Distributions.Distribution{Distributions.ArrayLikeVariate{N}, Distributions.Continuous}

• The most upstream abstract type for N-dimensional Random Fields.
• Set up most of the distribution-agnostic methods for Distribution, and useful methods for Base
• See here

Simple Uncorrelated Random Field
AbstractUnivariateRandomField{N} <: AbstractStationaryRandomField{N}

• This specific random field type will sample uncorrelated field from a univariate distribution.
• This can be used to make a Bayesian version of some RML regularizers in VLBISkyRegularizers.jl (so relevant to @anilipour ). It implements
  • Gaussian Fields, which can be used for l2-norm regularisation
  • Laplace Fields, which can be used for l1-norm regularization
• see here

Uncorrelated Random Field on Fourier Domain
AbstractRandomFourierField{N} <: AbstractStationaryRandomField{N}

• This specific random field type will sample uncorrelated field on the Fourier domain. The field is supposed to be a real function of the signal domain.
• see here
• I implemented two versions.
  • univariate Fourier field that just samples from a univariate distribution, but ensures that the fields will be real on the signal domain.
  • powerspectrum Fourier field that has an extra-scaling with AbstractPowerSpectrumModel (which doesn't change from the original Anna's implementation).

Power Spectrum Random Field
AbstractPowerSpectrumRandomField{M} <: AbstractStationaryRandomField{N}

• This specific random field type will sample correlated field on the signal domain, using AbstractPowerSpectrumFourierField on a Fourier basis
• This is of course for @ptiede, @rohandahale for general imaging and ScatteringOptics.jl
• This has a forward and inverse function to transform signals between Fourier and Singal planes.

Finally @anilipour

• I didn't implement wavelet field but @anilipour may add a Wavelet Laplace Field for wavelet L1 to make a full Bayesian compressed sensing imaging in VLBISkyRegularizers.jl in a similar way with Power Spectrum Random Field.