EHTJulia / StationaryRandomFields.jl

A Julia package to create noise signals following a specific power spectrum
MIT License
2 stars 0 forks source link

Implementation guidelines for Comrade compatibility #5

Open ptiede opened 1 month ago

ptiede commented 1 month ago

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 commented 1 month ago

@rohandahale