Comrade.jl

Appearance

Comrade API ​

Index ​

Model Definitions ​

Models ​

For the description of the model API see VLBISkyModels.

Data Interface ​

Data Tables ​

Comrade.AbstractVLBITable Type
julia
abstract type AbstractVLBITable{F}

An abstract VLBI table that is used to store data for a VLBI observation. To implement your own table you just need to specify the VLBISkyModels.rebuild function.

source
Comrade.datatable Function
julia
datatable(obs::AbstractVLBITable)

Construct a table from the observation obs. The table is usually a StructArray of fields

source
julia
datatable(obs::AbstractObservationTable)

Returns a tabular representation of the data. Note that for closures this ignores the covariance between quantities, which is otherwise included in the full EHTObservationTable.

source
Comrade.AbstractArrayConfiguration Type
julia
abstract type AbstractArrayConfiguration{F<:Comrade.AbstractBaselineDatum} <: Comrade.AbstractVLBITable{F<:Comrade.AbstractBaselineDatum}

This defined the abstract type for an array configuration. Namely, baseline times, SEFD's, bandwidth, observation frequencies, etc.

source
Comrade.EHTArrayBaselineDatum Type
julia
struct EHTArrayBaselineDatum{T, P, V} <: Comrade.AbstractBaselineDatum

A single datum of an ArrayConfiguration

source
Comrade.EHTArrayConfiguration Type
julia
struct EHTArrayConfiguration{A<:Comrade.AbstractBaselineDatum, F, T, S, D<:(AbstractArray{A<:Comrade.AbstractBaselineDatum})} <: Comrade.AbstractArrayConfiguration{A<:Comrade.AbstractBaselineDatum}

Table that specified pertinent details about the EHT array during an observation. These are typically items that are known before the observation is made.

Fields

  • bandwidth: Observing bandwith (Hz)

  • tarr: Telescope array file

  • scans: Scan times

  • mjd: modified julia date of the observation

  • ra: RA of the observation in J2000 (deg)

  • dec: DEC of the observation in J2000 (deg)

  • source: Common source name

  • timetype: Time zone used.

  • datatable: A struct array of EHTArrayBaselineDatum

source
Comrade.ClosureConfig Type
julia
struct ClosureConfig{F, A<:Comrade.AbstractArrayConfiguration{F}, D, V, E} <: Comrade.AbstractArrayConfiguration{F}

Array config file for closure quantities. This stores the design matrix designmat that transforms from visibilties to closure products.

Fields

  • ac: Array configuration for visibilities

  • designmat: Closure design matrix

  • vis: visibilities to closure design matrix

  • noise: visibility noises to closure design matrix

source
Comrade.sites Method
julia
sites(d::AbstractArrayConfiguration)

Get all the sites in a observation. The result is a vector of symbols.

source
Comrade.domain Method
julia
domain(ac; executor, header)

Get the u, v, time, freq domain of the array.

source
Comrade.beamsize Method
julia
beamsize(ac::AbstractArrayConfiguration)

Calculate the approximate beam size of the array ac as the inverse of the longest baseline distance.

source
Comrade.logclosure_amplitudes Function
julia
logclosure_amplitudes(vis::AbstractArray, d::DesignMatrix)

Compute the log-closure amplitudes for a set of visibilities with a design matrix d.

Notes

This uses a closure design matrix for the computation.

source
Comrade.closure_phases Function
julia
closure_phases(vis::AbstractArray, d::DesignMatrix)

Compute the closure phases for a set of visibilities and design matrix d

Notes

This uses a closure design matrix for the computation.

source
Comrade.AbstractObservationTable Type
julia
abstract type AbstractObservationTable{F<:Comrade.AbstractVisibilityDatum} <: Comrade.AbstractVLBITable{F<:Comrade.AbstractVisibilityDatum}

The abstract obervation table. This contains actual data plus the array configuration.

source
Comrade.EHTObservationTable Type
julia
struct EHTObservationTable{T<:Comrade.AbstractVisibilityDatum, S<:AbstractArray, E<:AbstractArray, A<:Comrade.AbstractArrayConfiguration} <: Comrade.AbstractObservationTable{T<:Comrade.AbstractVisibilityDatum}

The main data product type in Comrade this stores the data which can be a StructArray of any AbstractInterferometryDatum type. Note that the underlying structure is not part of the public API. Users should typically construct tables from the extract_table function.

Fields

  • measurement: Obervation measurement

  • noise: Observation thermal noise

  • config: Array config holds ancillary information about array

source
Comrade.domain Method
julia
domain(obs::AbstractObservationTable; executor=Serial(), header=ComradeBase.NoHeader()

Returns the u, v, time, frequency domain of the observation.

source
Comrade.arrayconfig Method
julia
arrayconfig(obs::AbstractObservationTable)
arrayconfig(obs::AbstractObservationTable, p::Symbol)

Returns the array configuration for a given observation. If p is provided then only the property p is returned.

source
Comrade.beamsize Method
julia
beamsize(obs::AbstractObservationTable)

Calculate the approximate beam size of the observation obs as the inverse of the longest baseline distance.

source
Comrade.sites Method
julia
sites(d::AbstractObservationTable)

Get all the sites in a observation. The result is a vector of symbols.

source
Comrade.TimeTable Type
julia
struct TimeTable{O<:Comrade.AbstractVLBITable, T, S}

Wraps EHTObservationTable in a table that separates the observation into scans. This implements the table interface. You can access scans by directly indexing into the table. This will create a view into the table not copying the data.

Example

julia
julia> st = timetable(obs)
julia> st[begin] # grab first scan
julia> st[end]   # grab last scan
source
Comrade.Scan Type
julia
struct Scan{T, I, S}

Composite type that holds information for a single scan of the telescope.

Fields

  • time: Scan time

  • index: Scan indices which are (scan index, data start index, data end index)

  • scan: Scan data usually a StructArray of a <:AbstractVisibilityDatum

source
Comrade.timetable Function
julia
timetable(obs::AbstractArrayConfiguration)

Reorganizes the observation into a table of scans, where scan are defined by unique timestamps. To access the data you can use scalar indexing

Example

julia
st = timetable(obs)
# Grab the first scan
scan1 = st[1]

# Acess the detections in the scan
scan1[1]

# grab e.g. the baselines
scan1[:baseline]
source
Comrade.flag Function
julia
flag(condition, data::AbstractObservationTable)

Flags the data in vis that satisfy the condition condition. For the opposite effect call filter(condition, data)

condition is a function that accepts a tuple (baseline, measurement, noise) and returns a boolean. baseline are the fields that define the measurement like in datatable(vis).baseline.

Example

To flag all data with a uvdistance less than 0.1e9 use:

julia
cond(x) = uvdist(x) < 0.1e9
flag_data(vis, cond)
source
Base.filter Method
julia
filter(condition, data::AbstractObservationTable)

Filters all the data in data that satisfy the condition condition, i.e. if condition returns true for a given datum, it is kept in the returned object.

This is the converse of flag.

source
Comrade.select_baseline Function
julia
select_baseline(data::AbstractObservationTable, bl::Tuple{Symbol, Symbol})

Selects the baseline bl from the data object. The baseline is a tuple of symbols, but the order of the symbols does not matter.

source
Comrade.add_fractional_noise Function
julia
add_fractional_noise(dvis::AbstractObservationTable, ferr)

Adds fractional noise to the data in dvis with the fractional error ferr, i.e. the returned data is identical except the noise is replaced with

√(σ² + ferr² * |V|²)

where σ is dvis[:noise] and V is dvis[:measurement].

Warning

When the datum type <: ClosureProducts, we do not update the noise covariance matrix. This means that cross-correlation may be slightly incorrect.

source
Comrade.load_array_txt Function
julia
load_array_txt(path)

Read an ehtim-style antenna text file and return a dictionary keyed by site symbol. Each value is a NamedTuple with the columns present in the file: (SEFD1, SEFD2, fr_parallactic, fr_elevation, fr_offset) plus the d-term fields (DR, DL) if present. fr_offset is converted from degrees (the ehtim file convention) to radians.

The file must begin with a commented header naming the columns, e.g.

julia
# site  X  Y  Z  SEFDR  SEFDL  DR_RE  DR_IM  DL_RE  DL_IM  FR_PAR  FR_ELEV  FR_OFF

Columns are matched by name (case-insensitive), so the order of fields and the relative position of the D-term and FR blocks are flexible. X, Y, Z are ignored here — antenna positions come from the uvfits AIPS-AN table. The result can be passed to any extract_* as array_overrides=load_array_txt("array.txt").

source
Comrade.reset_mounts! Function
julia
reset_mounts!(dvis::AbstractObservationTable, overrides)
reset_mounts!(dvis::AbstractObservationTable, file::String)

Update the mount-related columns of dvis's array table in place from a per-site overrides dict. Each entry is a Symbol => NamedTuple whose keys may be any of :SEFD1, :SEFD2, :fr_parallactic, :fr_elevation, :fr_offset. Missing keys keep the existing value.

The dict is most easily produced by load_array_txt. This is the recommended way to apply an ehtim-style antenna text file to a freshly extracted observation:

julia
dvis = extract_table(uvd, Coherencies(; time_average = VLBI.GapBasedScans()))
reset_mounts!(dvis, load_array_txt("array.txt"))
source

Datums ​

Comrade.AbstractVisibilityDatum Type
julia
AbstractVisibilityDatum

An abstract type for all VLBI data types. See Comrade.EHTVisibilityDatum for an example.

source
Comrade.EHTCoherencyDatum Type
julia
struct EHTCoherencyDatum{S, B<:Comrade.AbstractBaselineDatum, M<:(StaticArraysCore.SArray{Tuple{2, 2}, Complex{S}, 2}), E<:(StaticArraysCore.SArray{Tuple{2, 2}, S, 2})} <: Comrade.AbstractVisibilityDatum{S}

A Datum for a single coherency matrix

Fields

  • measurement: coherency matrix, with entries in Jy

  • noise: visibility uncertainty matrix, with entries in Jy

  • baseline: baseline information

source
Comrade.EHTVisibilityDatum Type
julia
struct EHTVisibilityDatum{Pol, S<:Number, B<:Comrade.AbstractBaselineDatum} <: Comrade.AbstractSinglePolDatum{Pol, S<:Number}

A struct holding the information for a single measured complex visibility.

FIELDS

  • measurement: Complex Vis. measurement (Jy)

  • noise: noise of the complex vis (Jy)

  • baseline: baseline information

source
Comrade.EHTVisibilityAmplitudeDatum Type
julia
struct EHTVisibilityAmplitudeDatum{P, S<:Number, B<:Comrade.AbstractBaselineDatum} <: Comrade.AbstractSinglePolDatum{P, S<:Number}

A struct holding the information for a single measured visibility amplitude.

FIELDS

  • measurement: amplitude (Jy)

  • noise: noise of the visibility amplitude (Jy)

  • baseline: baseline information

source
Comrade.EHTLogClosureAmplitudeDatum Type
julia
struct EHTLogClosureAmplitudeDatum{P, S<:Number, B<:Comrade.AbstractBaselineDatum} <: Comrade.ClosureProducts{P, S<:Number}

A Datum for a single log closure amplitude.

  • measurement: log-closure amplitude

  • noise: log-closure amplitude noise in the high-snr limit

  • baseline: baselines for the closure phase

source
Comrade.EHTClosurePhaseDatum Type
julia
struct EHTClosurePhaseDatum{P, S<:Number, B<:Comrade.AbstractBaselineDatum} <: Comrade.ClosureProducts{P, S<:Number}

A Datum for a single closure phase.

Fields

  • measurement: closure phase (rad)

  • noise: noise of the closure phase assuming the high-snr limit

  • baseline: baselines for the closure phase

source
Comrade.triangle Function
julia
triangle(b::EHTClosurePhaseDatum)

Returns the sites used in the closure phase triangle.

source
Comrade.baseline Function
julia
baseline(scan::Scan)

Return the baselines for each datum in a scan

source
Comrade.quadrangle Function
julia
quadrangle(b::EHTClosurePhaseDatum)

Returns the sites used in the closure amplitude quadrangle.

source

Data Products ​

Comrade.extract_table Function
julia
extract_table(obs, dataproducts::VLBIDataProducts)

Extract an Comrade.EHTObservationTable table of data products dataproducts. To pass additional keyword for the data products you can pass them as keyword arguments to the data product type. For a list of potential data products see subtypes(Comrade.VLBIDataProducts).

Example

julia
julia> dlcamp, dcphase = extract_table(obs, LogClosureAmplitudes(;snrcut=3.0), ClosurePhases(;snrcut=3.0, cut_trivial=true))
julia> dcoh = extract_table(obs, Coherencies())
julia> dvis = extract_table(obs, VisibilityAmplitudes())
source
Comrade.Visibilities Type
julia
Visibilities(;kwargs...)

Type to specify to extract the complex visibilities table in the extract_table function. Optional keywords are passed through extract_table to specify additional option.

Special keywords for eht-imaging with Pyehtim.jl

Any keyword arguments are ignored for now. Use eht-imaging directly to modify the data.

source
Comrade.VisibilityAmplitudes Type
julia
Visibilities(;kwargs...)

Type to specify to extract the log closure amplitudes table in the extract_table function. Optional keywords are passed through extract_table to specify additional option.

Special keywords for eht-imaging with Pyehtim.jl

For a list of potential keyword arguments see eht-imaging and add_amp command for obsdata.

source
Comrade.ClosurePhases Type
julia
ClosuresPhases(;kwargs...)

Type to specify to extract the closure phase table in the extract_table function. Optional keywords are passed through extract_table to specify additional option.

Special keywords for eht-imaging with Pyehtim.jl

For a list of potential keyword arguments see eht-imaging and add_cphase command for obsdata. In addition note we have changed the following:

  • count: How the closures are formed, the available options are "min-correct", "min", "max"

Warning

The count keyword argument is treated specially in Comrade. The default option is "min-correct" and should almost always be used. This option construct a minimal set of closure phases that is valid even when the array isn't fully connected. For testing and legacy reasons we ehtim other count options are also included. However, the current ehtim count="min" option is broken and does construct proper minimal sets of closure quantities if the array isn't fully connected.

source
Comrade.LogClosureAmplitudes Type
julia
LogClosureAmplitudes(;kwargs...)

Type to specify to extract the log closure amplitudes table in the extract_table function. Optional keywords are passed through extract_table to specify additional option.

Special keywords for eht-imaging with Pyehtim.jl

For a list of potential keyword arguments see eht-imaging and add_cphase command for obsdata. In addition note we have changed the following:

  • count: How the closures are formed, the available options are "min-correct", "min", "max"

Returns an EHTObservation with log-closure amp. datums

Warning

The count keyword argument is treated specially in Comrade. The default option is "min-correct" and should almost always be used. This option construct a minimal set of closure phases that is valid even when the array isn't fully connected. For testing and legacy reasons we ehtim other count options are also included. However, the current ehtim count="min" option is broken and does construct proper minimal sets of closure quantities if the array isn't fully connected.

source
Comrade.Coherencies Type
julia
Coherencies(;kwargs...)

Type to specify to extract the coherency matrices table in the extract_table function. Optional keywords are passed through extract_table to specify additional option.

Special keywords for eht-imaging with Pyehtim.jl

Any keyword arguments are ignored for now. Use eht-imaging directly to modify the data.

source

VLBI Modeling ​

Sky Models ​

Comrade.@sky Macro
julia
@sky function name(grid; kwargs...)
    param₁ ~ dist₁
    param₂ ~ dist₂
    # ... model body using params and kwargs ...
    return AbstractModel
end

Define a SkyModel constructor in a single block. Each name ~ dist line contributes an entry to the prior NamedTuple; everything else becomes the model body. The macro emits three definitions in the calling module:

  • _<name>_sky(θ, metadata) — the model function passed to SkyModel.

  • _<name>_prior(; kwargs..., grid) — builds the prior NamedTuple.

  • <name>(grid; kwargs...) — the user-facing constructor, returning a SkyModel.

The first positional argument is the image grid. Keyword arguments are stored in the model metadata (along with the grid) and are in scope for both the prior RHS expressions and the model body. The metadata destructure inside the model body only unpacks names that the body actually references.

Hierarchical priors

The RHS of ~ is any Julia expression evaluated in the prior-builder's scope. Use this to plug in pre-built compound priors:

julia
@sky function imager(grid; ftot, mimg, cprior)
    c    ~ cprior                                       # e.g. HierarchicalPrior or corr_image_prior
    σimg ~ VLBITruncated(VLBIGaussian(0.0, 0.5); lower = 0.0)
    fg   ~ VLBIUniform(0.0, 1.0)
    ρs   ~ ntuple(Returns(VLBIUniform(0.01, 10.0)), 3)  # tuple-of-IID priors
    rast = apply_fluctuations(CenteredLR(), mimg, σimg .* c.params)
    pimg = parent(rast)
    @. pimg = (ftot * (1 - fg)) * pimg
    return ContinuousImage(rast, BSplinePulse{3}())
end

Conditional priors where one tilde RHS references another tilde's drawn value are not supported directly. Express such cases by wrapping the chain in a HierarchicalPrior on a single tilde RHS — HierarchicalPrior nests, so HierarchicalPrior(base, HierarchicalPrior(...)) covers arbitrary-depth hyperprior chains.

source
Comrade.AbstractSkyModel Type
julia
AbstractSkyModel

The abstract type for Comrade Sky Models. For a concrete implementation see SkyModel.

Any subtype must implement the following methods

  • ObservedSkyModel(m::AbstractSkyModel, array::AbstractArrayConfiguration): Constructs an observed sky model given the sky model m and the array configuration array. This method is used to compute the visibilities and image of the sky model.

The following methods have default implementations:

  • idealmaps(datatype, m::AbstractSkyModel, x): Computes the ideal intensities and visibilities of the sky model m given the model parameters x.

  • skymodel(m::AbstractSkyModel, x): Returns the sky model image given the model parameters x.

  • domain(m::AbstractSkyModel): Returns the domain of the sky model m.

  • set_array(m::AbstractSkyModel, array::AbstractArrayConfiguration): Sets the array configuration for the sky model m and returns the observed sky model and prior.

  • set_prior(m::AbstractSkyModel, array::AbstractArrayConfiguration): Sets the prior for the sky model m given the array configuration array. This is used to set the priors for the model parameters.

source
Comrade.SkyModel Type
julia
SkyModel(f, prior, grid::AbstractRectiGrid; algorithm = NFFTAlg(), metadata=nothing)

Construct a sky model using the function map f with parameter priors prior, where the image model is defined on the domain grid. If the underlying model produced by f is non-analytic, then algorithm is used to numerically Fourier transform the model image. The metadata option contains additional information needed by the model f.

Arguments

  • f(x, p) : A function must be two arguments, where x are the model parameters and p is the metadata. typically x and p are named tuples.

  • prior : A named tuple of priors for the model parameters defined in x. Each model parameter used in x must have a defined element in the prior.

  • grid : The domain on which the model is defined. This defines the field of view and resolution of the model. Note that if f produces a analytic model then this field of view isn't used directly in the computation of the visibilities.

Optional Arguments

  • algorithm : The Fourier transform algorithm used to compute the visibilities. The default is NFFTAlg() which uses a non-uniform fast Fourier transform. Other options can be found by using subtypes(VLBISkyModels.FourierTransform)

  • metadata : Additional information needed by the model f. These are the addtional arguments passed to the model function f.

source
Comrade.FixedSkyModel Type
julia
FixedSkyModel(m::AbstractModel, grid::AbstractRectiGrid; algorithm = NFFTAlg())

Construct a sky model that has no free parameters. This is useful for models where the image structure is known apriori but the instrument model is unknown.

Arguments

  • m : The fixed sky model.

  • grid : The domain on which the model is defined. This defines the field of view and resolution of the model. Note that if f produces a analytic model then this field of view isn't used directly in the computation of the visibilities.

Optional Arguments

  • algorithm : The Fourier transform algorithm used to compute the visibilities. The default is NFFTAlg() which uses a non-uniform fast Fourier transform. Other options can be found by using subtypes(VLBISkyModels.FourierTransform)
source
Comrade.MultiSkyModel Type
julia
MultiSkyModel(skymodels::NamedTuple)

Create a sky model that is a collection of multiple components. This is useful when you, e.g., want to decompose the sky into multiple grids, such as for wide-field imaging where the is a core component and a very far away component.

Arguments

  • skymodels : A named tuple of sky models, where each model is a SkyModel.

Example

julia
julia> s1 = SkyModel(...)
julia> s2 = SkyModel(...)
julia> stot = MultiSkyModel((core = m1, far = m2))
julia> skymodel(stot, x)
(core = ..., far = ...)
source
Comrade.idealmaps Function
julia
idealmaps(D, m::AbstractSkyModel, x)

Computes the ideal non-corrupted visibilities of the sky model m given the model parameters x.

source
Comrade.skymodel Method
julia
skymodel(post::AbstractVLBIPosterior, θ)

Returns the sky model or image of a posterior using the parameter valuesθ

source

Instrument Models ​

Comrade.@instrument Macro
julia
@instrument function name(; refbasis = CirBasis(), kwargs...)
    G = @jones begin ... end
    D = @jones begin ... end
    # ... build and compose the Jones terms ...
    return AbstractJonesMatrix
end

Define an InstrumentModel constructor in a single block, mirroring @sky. Each Jones term is built with a @jones block, which carries its own ~ priors and param_map; the priors of all blocks are collected into the instrument prior. The body then composes the built Jones terms (e.g. with JonesSandwich) and returns the result.

The instrument model takes no positional arguments. Keyword arguments are tunable constants in scope for the ~ prior expressions and the @jones bodies. refbasis is treated specially: it is forwarded only to InstrumentModel and defaults to CirBasis() if not listed. Priors are declared only inside @jones blocks — a top-level ~ is an error. ~ lines are optional overall, so a parameter-free response model is allowed.

The JonesSandwich combination function is lifted to a named function automatically (named functions, unlike closures, serialize reliably).

Example:

julia
@instrument function fullpol(; refbasis = CirBasis(), gain_std = 0.1, dterm_std = 0.2)
    G = @jones begin
        lgR   ~ ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, gain_std)))
        gpR   ~ ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2)));
                           refant = SEFDReference(0.0), phase = true)
        lgrat ~ ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, gain_std)))
        gprat ~ ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)))
        gR = exp(complex(lgR, gpR))
        gL = gR * exp(complex(lgrat, gprat))
        return JonesG((gR, gL))
    end
    D = @jones begin
        dRx ~ ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, dterm_std)))
        dRy ~ ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, dterm_std)))
        dLx ~ ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, dterm_std)))
        dLy ~ ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, dterm_std)))
        return JonesD((complex(dRx, dRy), complex(dLx, dLy)))
    end
    R = @jones begin
        return JonesR(; add_fr = true)
    end
    return JonesSandwich(G, D, R) do g, d, r
        adjoint(r) * g * d * r
    end
end

intm = fullpol(; gain_std = 0.05)

@instrument selfcal = true function ... builds the model with selfcal = true (see InstrumentModel). Its Jones terms must then omit the station phase, which the model sets to the self-calibration solution plus the sampled residual selfcal_residual.

A Stokes-I example is a single block:

julia
@instrument function gainsonly(; refbasis = CirBasis())
    return @jones begin
        lg ~ ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)))
        gp ~ ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))))
        return SingleStokesGain(exp(complex(lg, gp)))
    end
end
source
Comrade.@jones Macro
julia
@jones begin
    param₁ ~ ArrayPrior(...)
    # ... statements building the matrix entries from the params ...
    return JonesG(entries)
end

Build one Jones term inside an @instrument block. Each name ~ ArrayPrior(...) line is a prior for this term; the rest of the block is how the we paraemeterized the Jones matrix.

A block with no ~ lines is parameter-free and is emitted verbatim (e.g. JonesR):

julia
R = @jones begin
    return JonesR(; add_fr = true)
end

@jones is only meaningful inside @instrument (its ~ priors have nowhere else to go); using it standalone with ~ lines is an error.

source
Comrade.CalTable Type
julia
struct CalTable{T, F, G<:(AbstractVecOrMat)}

A Tabes of calibration quantities. The columns of the table are the telescope sites codes. The rows are the calibration quantities at a specific time and frequency. This user should not use this struct directly. Instead that should call caltable.

source
Comrade.caltable Method
julia
caltable(s::SiteArray)

Creates a calibration table from a site array

source
Comrade.siteparams Function
julia
siteparams(x)

Extract the instrument parameter values from a single instrument-parameter sample. For hierarchical samples (params = SiteArray, hyperparams = NamedTuple) (produced by GaussMarkovSitePrior with fitted hyperparameters) this returns x.params; for plain SiteArray samples it is the identity.

source
Comrade.sites Method
julia
sites(g::CalTable)

Return the sites in the calibration table

source
Comrade.IIDSitePrior Type
julia
IIDSitePrior(seg::Segmentation, dist; init = nothing)

Create a site prior that is independent and identically distributed (IID) across all times and frequencies. The seg argument is a segmentation object that defines how fine the time segmentation is. The dist argument is the distribution of the site prior.

init = FixedInit(v) fixes each site's first time stamp (per frequency channel) to v instead of sampling it, as the same option does for GaussMarkovSitePrior. Use it when the parameter is added to a separate per-site offset term: otherwise the offset and a common shift of the site's values are degenerate, and for a phase the 2π periodicity of the likelihood turns that degeneracy into extra modes. A stamp that the referencing scheme (refant) also fixes must get the same value.

Example

julia
A = IIDSitePrior(ScanSeg(), VLBIGaussian(0, 1))

creates a site prior that is constant across scans and each scan has a unit Normal prior.

source
Comrade.GaussMarkovSitePrior Type
julia
GaussMarkovSitePrior(seg::TimeSegmentation, process::AbstractGaussMarkovProcess; init = StationaryInit(), param = NonCentered())

A site prior correlated in time, following the continuous-time Gauss-Markov process (e.g. OrnsteinUhlenbeck) at the times given by the segmentation seg. The log-density factorizes over each site's time chain using the process's exact discretization, so it costs O(n) in the number of times and handles irregular time sampling exactly. With multiple frequency channels, each (site, channel) pair is an independent chain.

Process fields that are Distributions.Distributions are fitted hyperparameters. When any are present the prior sample is a NamedTuple (params = SiteArray, hyperparams = NamedTuple), as for HierarchicalPrior; otherwise the sample is a plain SiteArray. Each site gets its own hyperparameters: the process passed to ArrayPrior is a per-site template, and each site's fitted fields are nested under the site's name (e.g. x.instrument.lg.hyperparams.AA.σ). The usual ArrayPrior per-site overrides replace the template, or the whole prior, for that site. caltable/plotcaltable accept the hierarchical sample directly; for elementwise postprocessing use the params field, e.g. exp.(x.instrument.lg.params).

Parameterization

The param keyword picks the chain's latent coordinates: NonCentered, Centered, or, for wrapped processes only, AngleEmbedded. centered = true / centered = false is shorthand for Centered() / NonCentered(); passing both param and centered throws. The default is NonCentered for every process that supports it, which is all of them except WrappedOrnsteinUhlenbeck.

Centered coordinates give the sampler the chain's prior correlation. For a random walk the condition number of that correlation grows with the number of time stamps, and a diagonal mass matrix cannot absorb it. A common shift of one station's whole phase chain, which trades against source position and structure, is then a diagonal direction through every coordinate of the chain; in non-centered coordinates it is a single coordinate (the chain's level), and the increments are prior-standardized. Centered can mix better when every point is strongly constrained by the data and the prior is weak.

For a real-line process (OrnsteinUhlenbeck, BrownianMotion), NonCentered also removes the funnel between fitted hyperparameters and gains and makes the transport to the Std spaces exact, so both asflat and ascube work. Centered supports asflat only.

For a wrapped (circular) process every parameterization supports asflat only. The chain density is 2π-periodic, so a real coordinate read directly as an angle gives an improper flat target, with identical modes spaced 2π apart in every coordinate. Centered and NonCentered therefore multiply each free phase's flat density by the sheet weight N(δ; 0, Q) / WN(δ; Q), where δ is its step from the previous chain value and Q that step's variance. The weights sum to one over the 2π sheets, so the circular model is unchanged while the flat target becomes proper; a 2π step costs exp(−(2π)²/2Q).

Under Centered the flat target is the unwrapped Gaussian random walk in the raw angles. Under NonCentered each free phase's coordinate is its increment divided by the step's standard deviation, the flat target is N(0, I) in those coordinates, and the phase is the wrapped running sum. In both, the chain's first free point keeps a raw-angle coordinate. For a non-centered wrapped chain, logprior of a constrained sample contains only the chain-opening terms, and the rest of the density is in the flat transform's log-Jacobian; only their sum, which samplers and optimizers evaluate, is the chain density. On the circle NonCentered requires a transition mean that carries the previous point's sheet: WrappedBrownian (Φ ≡ 1) or the 2π-equivariant sine drift of VonMisesProcess, whose recurrence is resolved sequentially. The shortest-arc drift of WrappedOrnsteinUhlenbeck does not, so that process is Centered only.

AngleEmbedded embeds each free phase as two latent reals through the angle transform of DiagonalVonMises. It has no sheet structure (a phase wrap is a continuous winding of the latent point), at the cost of twice the phase dimension, a per-pair geometry that a diagonal mass matrix cannot capture, and the chain correlation of the centered form.

Phases

For an absolute phase use WrappedBrownian with init = UniformInit(): the prior is 2π-periodic, needs no separate circular offset term, and its coherence-time posterior is not biased by phase wraps. For a phase held near a level, such as an R/L gain ratio phase, use WrappedOrnsteinUhlenbeck, whose WN(μ, σ²) marginal keeps the chain from wandering. A real-line process (e.g. OrnsteinUhlenbeck with init = FixedInit(0.0) on top of a circular offset) models the phase as an unwrapped real and is appropriate only when excursions stay well below 2π; larger excursions produce 2π-shifted posterior modes from the likelihood's periodicity in exp(iθ). The phase = true option of ArrayPrior is not supported with this prior and throws.

Initial distribution

init (see AbstractInitialPrior) sets the distribution of each chain's first time stamp. StationaryInit (the default) uses the process's stationary marginal. GaussianInit uses an explicit N(μ0, σ0²), e.g. a wide one so the starting value is nearly free. FixedInit conditions the chain to an exact value at its first time, the continuous-time analogue of phase = true. For an unwrapped phase chain on top of a separate circular offset, FixedInit(0.0) removes the level freedom that would otherwise trade against the offset and, through the likelihood's 2π periodicity, produce 2π-shifted posterior modes. If the reference scheme also fixes a chain's first point, the two values must agree, or construction throws.

Reference stations

On the real line a reference-fixed value is a value of the process, and the free points are drawn from p(free | fixed) exactly. Scattered fixed indices, such as those from SEFDReference, are supported.

On the circle a reference-fixed value is a gauge choice. A station's phase is defined only up to the per-time global phase that the reference scheme removes, so being the reference says nothing about the station's own gain phase, and conditioning on it would pull the station's whole chain toward the reference value. A wrapped chain is therefore built over the time stamps where the station is free, with transitions taken over the combined gap across skipped stamps: a point after a reference-fixed stamp is an offset from the last free point. The fixed values still enter the likelihood but carry no prior term. For WrappedBrownian this chain is the exact marginal of the full one.

Example

julia
lg ~ ArrayPrior(GaussMarkovSitePrior(IntegSeg(), OrnsteinUhlenbeck(σ = VLBIExponential(0.1), τ = VLBIExponential(2.0))))
gp ~ ArrayPrior(GaussMarkovSitePrior(ScanSeg(), OrnsteinUhlenbeck(σ = 3.0, τ = 1.0)); refant = SEFDReference(0.0))
## phase-fluctuation term meant to be combined with a circular per-track offset
dgp ~ ArrayPrior(GaussMarkovSitePrior(ScanSeg(), OrnsteinUhlenbeck(σ = 1.0, τ = 1.0); init = FixedInit(0.0)); refant = SEFDReference(0.0))

Note

The example uses the VLBI* distributions (VLBIExponential, VLBIInverseGamma, ...) for fitted hyperparameters. These accept a scalar traced argument and so work on both the CPU and GPU/Reactant backends; the plain Distributions.jl types work only on the CPU.

source
Comrade.AbstractChainParameterization Type
julia
AbstractChainParameterization

Abstract type for the latent coordinates of a GaussMarkovSitePrior chain: Centered, NonCentered, or AngleEmbedded. See the site prior's docstring for which to use when.

source
Comrade.Centered Type
julia
Centered()

Chain coordinates equal to the free points' values: the gain itself or, for a wrapped process, the raw angle with a sheet weight. The flat target keeps the chain's prior correlation. Only asflat is supported. See GaussMarkovSitePrior.

source
Comrade.NonCentered Type
julia
NonCentered()

Whitened chain coordinates: each free point's coordinate is an iid standard normal variate, so the flat target has none of the chain's prior correlation. On the real line the coordinate is mapped to the value through the chain's exact conditional moments, including the reference-station bridge and the initial marginal. On the circle it is the step onto the previous phase's transition mean divided by the step's standard deviation, and the chain's first free point keeps a raw-angle coordinate. Requires a transition mean that is affine in the previous point or 2π-equivariant (like the sine drift of VonMisesProcess), which excludes WrappedOrnsteinUhlenbeck. This is the default. See GaussMarkovSitePrior.

source
Comrade.AngleEmbedded Type
julia
AngleEmbedded()

Chain coordinates that embed each free phase as two latent reals through the angle transform of DiagonalVonMises, doubling the dimension. These coordinates have no 2π sheet structure. Valid only for a wrapped process. See GaussMarkovSitePrior.

source
Comrade.AbstractGaussMarkovProcess Type
julia
AbstractGaussMarkovProcess

Abstract type for continuous-time Gauss-Markov processes in state-space form, used to build time-correlated instrument priors (see GaussMarkovSitePrior) whose log-density factorizes over the time chain and so costs O(n) in the number of times.

Interface

A concrete process p must implement

  • isstationary(p) -> Bool: whether the process has a stationary distribution.

  • stationary_moments(p) -> (μ, P): the stationary mean and variance (stationary processes only).

  • transition_moments(p, Δt) -> (Φ, Q): the exact discretization over a gap Δt (in hours), so that x(t + Δt) | x(t) ~ N(μ + Φ*(x(t) - μ), Q).

The distribution of each chain's first time stamp is set by the site prior's AbstractInitialPrior, not by the process.

Any field of the process that is a Distributions.Distribution is a fitted hyperparameter with that distribution as its prior (see hyperprior); it is replaced by a sampled value via materialize before stationary_moments/transition_moments are called. Numeric fields are fixed. The generic materialize uses ConstructionBase.setproperties; a concrete type can customize reconstruction via ConstructionBase.constructorof.

Note

Only scalar-state processes are supported.

source
Comrade.OrnsteinUhlenbeck Type
julia
OrnsteinUhlenbeck(; σ, τ, μ = 0.0)

A stationary Ornstein-Uhlenbeck process, i.e. the continuous-time AR(1) process. This is the mean-reverting Gaussian process with marginal x(t) ~ N(μ, σ²) and exponential correlation cov(x(t), x(t+Δ)) = σ² exp(-|Δ|/τ), whose exact discretization over a gap Δt is

julia
x(t + Δt) | x(t) ~ N(μ + φ*(x(t) - μ), σ²*(1 - φ²)),   φ = exp(-Δt/τ).

Arguments

  • σ: the marginal (stationary) standard deviation. Pass a Number to fix it or a Distributions.Distribution to fit it as a hyperparameter with that prior.

  • τ: the correlation timescale in hours (the times of a VLBI observation are in UTC hours). A Number fixes it; a Distributions.Distribution fits it.

  • μ: the stationary mean (default 0.0). Not fittable.

Note

Fitted hyperparameters are evaluated pointwise, so on GPU/Reactant backends their priors must accept a scalar traced argument. Use the VLBI* distributions (e.g. VLBIExponential, VLBIInverseGamma) rather than the plain Distributions.jl types for a model you intend to run through Reactant.

Example

julia
# fixed hyperparameters
OrnsteinUhlenbeck(σ = 0.1, τ = 2.0)
# fitted marginal std and correlation time (VLBI* variants work on CPU and Reactant)
OrnsteinUhlenbeck(σ = VLBIExponential(0.1), τ = VLBIInverseGamma(3.0, 6.0))
source
Comrade.BrownianMotion Type
julia
BrownianMotion(; D)

A driftless Brownian motion (Wiener process) with diffusion coefficient D: increments over a gap Δt (in hours) are independent N(0, D*Δt), so √D is the typical drift per √hour. This is the non-reverting limit of OrnsteinUhlenbeck (τ → ∞ at fixed D = 2σ²/τ), for quantities that wander without a preferred level.

The process is nonstationary (the marginal variance grows linearly in time), so a chain needs an explicit initial prior, GaussianInit or FixedInit; StationaryInit throws.

Arguments

  • D: the diffusion coefficient, in units of the parameter squared per hour. Pass a Number to fix it or a Distributions.Distribution to fit it as a hyperparameter with that prior.
source
Comrade.WrappedBrownian Type
julia
WrappedBrownian(; τ)

Brownian motion on the circle: angular increments over a gap Δt (hours) are wrapped normal WN(0, 2Δt/τ), the wrap of BrownianMotion onto (−π, π]. Use it for gain phases: the prior is 2π-periodic, so 2π-shifted posterior modes are equivalent (the likelihood sees only exp(iθ)), and the coherence-time posterior is not biased by phase wraps.

In the visibility domain the process decorrelates exponentially, E[e^{i(θ(t+Δ) − θ(t))}] = e^{−Δ/τ}, so τ is the phase coherence time in hours (the diffusion coefficient is D = 2/τ rad² per hour). The circular stationary law is uniform, so with UniformInit the chain also absorbs the per-track phase offset and no separate circular offset term is needed.

The transition mean is the previous phase, so this process supports param = NonCentered() (the default): each free phase's coordinate is its increment divided by the step's standard deviation, and the phase is the wrapped running sum. param = Centered() uses the raw angles. Both use sheet weights to keep the flat target proper without changing the circular model; AngleEmbedded avoids them at twice the phase dimension. Only asflat is supported. See GaussMarkovSitePrior.

A reference-fixed time stamp is a gauge choice, so it is dropped from the chain and the free points step over the combined gap (wrapped normal transitions compose: their variances add). A reference station's phases at other times are offsets from its last free phase, not pulled toward the reference value.

Arguments

  • τ: the phase coherence time in hours (the times of a VLBI observation are in UTC hours): the gap over which the visibility-domain coherence E[e^{iΔθ}] falls by 1/e. Pass a Number to fix it or a Distributions.Distribution to fit it as a hyperparameter with that prior.

Example

julia
gp ~ ArrayPrior(
    GaussMarkovSitePrior(IntegSeg(), WrappedBrownian(τ = VLBIInverseGamma(3.0, 6.0)); init = UniformInit());
    refant = SEFDReference(0.0)
)
source
Comrade.WrappedOrnsteinUhlenbeck Type
julia
WrappedOrnsteinUhlenbeck(; σ, τ, μ = 0.0)

The mean-reverting process on the circle: the wrap of OrnsteinUhlenbeck onto (−π, π], with the drift taken along the shortest arc so that the law is 2π-periodic in both of its arguments and is a circular Markov kernel,

julia
θ(t + Δt) | θ(t) ~ WN(μ + φ*wrap(θ(t) - μ), σ²*(1 - φ²)),   φ = exp(-Δt/τ),

where wrap maps to (−π, π]. Its circular marginal is the wrapped normal WN(μ, σ²), so, unlike WrappedBrownian, the phase stays near a level, and the visibility-domain coherence

julia
E[e^{i(θ(t+Δ) − θ(t))}] = exp(−σ²(1 − e^{−Δ/τ}))

levels off at exp(−σ²) rather than decaying to zero. Taking σ, τ → ∞ at fixed 2σ²/τ = 2/τ_coh recovers WrappedBrownian(τ = τ_coh).

Use it for a gain ratio (R/L) phase, an instrumental offset that is stable up to slow drift: one chain replaces a WrappedBrownian chain with FixedInit plus a separate circular offset, and σ sets how far the ratio phase may stray. The shortest-arc drift drops the previous phase's sheet, so this process supports only param = Centered() and asflat; VonMisesProcess is a non-centered alternative. As for every wrapped process, reference-fixed time stamps are dropped from the chain. See GaussMarkovSitePrior.

Arguments

  • σ: the circular marginal spread in radians, i.e. the WN(μ, σ²) scale. Pass a Number to fix it or a Distributions.Distribution to fit it as a hyperparameter.

  • τ: the reversion timescale in hours (the times of a VLBI observation are in UTC hours), the same units as the τ of OrnsteinUhlenbeck and WrappedBrownian. A Number fixes it; a Distributions.Distribution fits it.

  • μ: the circular mean (default 0.0). Not fittable.

Valid for σ well below π

Three properties hold only up to the WN(μ, σ²) mass beyond μ ± π: stationarity of the marginal, exact composition of the reference-antenna conditioning of GaussMarkovSitePrior over skipped time stamps, and continuity of the kernel at the antipode of μ. That mass is 1.7e-3 at σ = 1 rad and e^{-123} at σ = 0.2 rad. Do not use this process as an uninformative circular prior (σ ≳ 2); use WrappedBrownian with UniformInit instead.

Example

julia
## an R/L ratio phase: pinned near zero, drifting on a 12 hour timescale
gprat ~ ArrayPrior(
    GaussMarkovSitePrior(
        ScanSeg(), WrappedOrnsteinUhlenbeck(σ = VLBIExponential(0.2), τ = VLBIInverseGamma(1.0, 12.0))
    );
    refant = SingleReference(:AA, 0.0)
)
source
Comrade.VonMisesProcess Type
julia
VonMisesProcess(; σ, τ, μ = 0.0)

The mean-reverting circular process with a smooth drift. It has the same σ, τ, μ, and transition variance σ²(1 − φ²), φ = exp(−Δt/τ), as WrappedOrnsteinUhlenbeck, with the shortest-arc drift replaced by the sine drift

julia
θ(t + Δt) | θ(t) ~ WN(θ(t) − (1 − φ)·sin(θ(t) − μ), σ²*(1 − φ²)),

a discretization of the von Mises (circular Langevin) diffusion, whose stationary law is vM(μ, 1/σ²) ≈ WN(μ, σ²) for small σ. For deviations well inside (−π, π) the two drifts agree (sin ε ≈ ε), so this process can replace WrappedOrnsteinUhlenbeck in the σ ≪ π regime, with the same warning: do not use it as an uninformative circular prior; use WrappedBrownian there.

The sine mean is smooth and 2π-equivariant in the previous point, whereas the shortest-arc drift jumps by 2πφ at the antipode of μ, so this process supports param = NonCentered() (the default), whose coordinates decouple the chain values from fitted hyperparameters and from each other. The drift also vanishes at the antipode, an unstable equilibrium where the stationary density is smallest; for σ ≪ π the chain almost never reaches it.

The drift is nonlinear, so a non-centered chain is resolved as a sequential recurrence rather than the batched affine scan; under Reactant this is a serialized loop. It suits scan-level chains, e.g. a gain ratio phase about a separate circular offset with init = FixedInit(0.0), or carrying its own level with the default StationaryInit. For an integration-level absolute phase use WrappedBrownian.

Arguments

  • σ: the circular marginal spread in radians. Pass a Number to fix it or a Distributions.Distribution to fit it as a hyperparameter.

  • τ: the reversion timescale in hours (the times of a VLBI observation are in UTC hours). A Number fixes it; a Distributions.Distribution fits it.

  • μ: the circular mean (default 0.0). Not fittable.

Example

julia
## an R/L ratio phase: pinned near zero, drifting on a 12 hour timescale, whitened
gprat ~ ArrayPrior(
    GaussMarkovSitePrior(
        ScanSeg(), VonMisesProcess(σ = VLBIExponential(0.2), τ = VLBIInverseGamma(1.0, 12.0))
    );
    refant = SingleReference(:AA, 0.0)
)
source
Comrade.stationary_moments Function
julia
stationary_moments(p::AbstractGaussMarkovProcess) -> (μ, P)

Return the stationary mean μ and variance P of the process p.

source
Comrade.transition_moments Function
julia
transition_moments(p::AbstractGaussMarkovProcess, Δt) -> (Φ, Q)

Return the transition coefficient Φ and innovation variance Q for the exact discretization of p over a time gap Δt in hours, so that x(t + Δt) | x(t) ~ N(μ + Φ*(x(t) - μ), Q).

source
Comrade.isstationary Function
julia
isstationary(p::AbstractGaussMarkovProcess) -> Bool

Whether p has a stationary distribution, i.e. implements stationary_moments. StationaryInit, the default initial prior of GaussMarkovSitePrior, requires it. Every process must implement this; there is no fallback.

source
Comrade.AbstractInitialPrior Type
julia
AbstractInitialPrior

Abstract type for the initial distribution p(x(t₁)) of a GaussMarkovSitePrior chain, i.e. how each (site, frequency) chain's first time stamp is treated:

  • StationaryInit: the process's stationary marginal (the default; requires a stationary process).

  • GaussianInit: an explicit normal marginal N(μ0, σ0²).

  • FixedInit: the first point is held at an exact value by conditioning.

  • UniformInit: uniform on the circle, for nonstationary wrapped processes.

On the real line, initial priors must be Gaussian (a delta is the zero-variance limit): the exact chain density, the reference-antenna conditioning, and the whitened transports use closed-form Gaussian conditioning.

Interface

  • initial_moments(init, p) -> (m₀, P₀): mean and variance of the first time stamp.
source
Comrade.StationaryInit Type
julia
StationaryInit()

Start each chain in the process's stationary marginal from stationary_moments: N(μ, P) for a real-line process and the wrapped normal WN(μ, P) for a circular one (e.g. WrappedOrnsteinUhlenbeck). Every time stamp then has the same marginal. The default of GaussMarkovSitePrior. Requires isstationary(process).

source
Comrade.FixedInit Type
julia
FixedInit(value)

Hold each chain at value at its first time stamp by exact conditioning, so the prior describes the quantity's evolution from a known start. Composes with reference-fixed (refant) indices; if the reference scheme fixes a chain's first point to a different value, posterior construction throws.

FixedInit(0.0) on a phase-fluctuation chain on top of a separate circular offset term removes the level redundancy between the two, which the likelihood's 2π periodicity would otherwise turn into posterior modes at 2π-shifted levels.

source
Comrade.GaussianInit Type
julia
GaussianInit(μ0, σ0)

Start each chain with x(t₁) ~ N(μ0, σ0²), independent of the process's stationary marginal (if any). A large σ0 leaves the first point nearly free. σ0 must be positive, or the constructor throws; for an exact starting value use FixedInit.

source
Comrade.UniformInit Type
julia
UniformInit()

Start each chain uniformly on the circle (−π, π]. Valid only for wrapped processes that are not mean-reverting, such as WrappedBrownian, whose transitions leave the uniform law invariant; the chain then absorbs the per-track phase offset. A mean-reverting circular process (WrappedOrnsteinUhlenbeck) uses StationaryInit instead.

source
Comrade.initial_moments Function
julia
initial_moments(init::AbstractInitialPrior, p::AbstractGaussMarkovProcess) -> (m₀, P₀)

Mean and variance of a chain's first time stamp under init.

source
Comrade.marginal_moments Function
julia
marginal_moments(init, p, Δt) -> (m, P)

Moments of the chain marginal a gap Δt after the chain start: the initial moments propagated through the transition law, m = μ + Φ(m₀ − μ), P = Φ²P₀ + Q. At Δt == 0 this is initial_moments exactly; for StationaryInit the marginal is the stationary one at every gap.

source
Comrade.hyperprior Function
julia
hyperprior(p::AbstractGaussMarkovProcess) -> NamedTuple

Return the priors of the fitted hyperparameters of p, i.e. the fields of p that are Distributions.Distributions. Returns an empty NamedTuple when every field is a number, in which case the process has no fitted hyperparameters.

source
Comrade.materialize Function
julia
materialize(p::AbstractGaussMarkovProcess, hp::NamedTuple)

Substitute the fitted hyperparameters of p with the values in hp, returning a fully numeric process. Fields of p that are Distributions.Distributions are replaced by hp.<fieldname>; numeric fields are kept. When p has no fitted hyperparameters this is the identity.

source
Comrade.ArrayPrior Type
julia
ArrayPrior(default_dist; refant=NoReference(), phase=false, gauge=:none, kwargs...)

Construct a prior for an entire array of sites.

  • The default_dist is the default distribution for all sites. This may be an

IIDSitePrior or a time-correlated GaussMarkovSitePrior.

  • Different priors for specified sites can be set using kwargs.

  • The refant set the reference antennae to be used and is typically only done for priors that

correspond to gain phases.

  • The phase argument is a boolean that specifies if

the prior is for a phase or not. The phase argument is experimental and we recommend setting it to false currently.

  • The gauge argument is :none (the default) or :phase. gauge = :phase marks the

parameter as one term of the station phase. Baseline phases depend only on differences of station phases, so the sum of the gauge = :phase terms at a site is determined only up to one constant per (time, frequency). When several such terms are added (for example a track-level offset plus a per-integration residual), references placed on each term can leave the sum undetermined. The InstrumentModel then finds the remaining flat directions with Comrade.gauge_pins and, depending on its gaugefix setting, throws or pins them.

Example

julia
p = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0, 0.1)); LM = IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 1.0)) refant=SEFDReference())

means that every site has a normal prior with mean 0 and 0.1 std. dev. except LM which is mean zero and unit std. dev. Finally the refant is using the SEFDReference scheme.

source
Comrade.Segmentation Type
julia
abstract type Segmentation

The data segmentation scheme to use. This is used specify how often we want various instrument hyperparameters to change. A user subtyping this expression must implement the following functions:

  • timestamps: Specifies the time region for each segmentation scheme given an array
source
Comrade.IntegSeg Type
julia
struct IntegSeg <: Comrade.TimeSegmentation

Data segmentation such that the quantity is constant over the time stamps in the data. If the data is scan-averaged before then IntegSeg will be identical to ScanSeg.

source
Comrade.ScanSeg Type
julia
struct ScanSeg <: Comrade.TimeSegmentation

Data segmentation such that the quantity is constant over a scan.

source
Comrade.TrackSeg Type
julia
struct TrackSeg <: Comrade.TimeSegmentation

Data segmentation such that the quantity is constant over a track, i.e., the observation "night".

source
Comrade.timestamps Function
julia
timestamps(seg::TimeSegmentation, array::AbstractArrayConfiguration)

Return the time stamps or really a vector of integration time regions for a given segmentation scheme seg and array configuration array.

source
Comrade.SingleReference Type
julia
SingleReference(site::Symbol, val)

Selects a single reference site for all scans. The value of the site is set to val.

source
Comrade.MultiReference Type
julia
MultiReference(sites, val)

Selects several reference sites, fixing every entry of each site in sites to val. sites is any iterable of site names.

Use it for a parameter with more than one gauge freedom. The usual case is a track-segmented phase offset when the schedule splits the array into groups that share no station, for example when a subset of the array observes for a stretch without the reference site. Each such group needs its own reference.

Notes

This pins the sites it is given and does not check that the result leaves the model identified. Too few pins leave a flat direction in the posterior; too many over-constrain a scan, since a scan carries only one global phase.

source
Comrade.SEFDReference Type
julia
SEFDReference(val::Number, sefd_index = 1)

Selects the reference site based on the SEFD of each telescope, where the smallest SEFD is preferentially selected. The reference gain is set to val and the user can select to use the n lowest SEFD site by passing sefd_index = n.

Notes

This is done on a per-scan basis so if a site is missing from a scan the next highest SEFD site will be used.

source
Comrade.gauge_pins Function
julia
gauge_pins(terms, array::AbstractArrayConfiguration; prefer = default_gauge_preference)
gauge_pins(terms, stamps, sites_at, sefd; prefer = default_gauge_preference)

Return the entries that must be pinned, as a Vector{Tuple{Symbol, Int}} of (term name, entry index), so that a station phase written as the sum of terms is fully determined by the data. An empty result means the parameterization is already identified.

terms is a NamedTuple mapping each term's name to (; smap, fixed), where smap is the term's SiteLookup and fixed holds the entry indices that its referencing scheme and init pins hold constant. The first form reads the schedule off an array configuration. The second takes stamps, a vector of (time, frequency) pairs sorted by time; sites_at, the sites present at each stamp; and sefd, a site-keyed NamedTuple or dictionary of system equivalent flux densities. The sites of one stamp must be connected by baselines: a (time, frequency) whose baselines split into disconnected groups is listed once per group.

Each datum (site s at stamp τ) gives one equation relating the stamp's phase constant to the free entries covering (s, τ). An entry covers a datum when its site matches and its integration time and frequency channel contain the stamp. The solver:

  1. Removes each equation with three or more unknowns that contains an unknown appearing in no other equation, and marks that unknown as determined (a stamp constant if there is one, otherwise the least preferred entry). This repeats until no equation can be removed.

  2. Treats the remaining equations as graph edges: an equation with one unknown determines it, and one with two unknowns links them. A remaining equation with three or more unknowns throws an ArgumentError.

  3. Pins one entry in each connected group of unknowns that no equation determines, chosen by prefer; see default_gauge_preference.

A group containing only stamp constants needs no pin.

Notes

InstrumentModel applies this analysis, through its gaugefix setting, to the terms declared with gauge = :phase.

source
Comrade.default_gauge_preference Function
julia
default_gauge_preference(entry)

Rank the entries that gauge_pins may pin: the entry with the smallest key is the one it pins, and the largest key is the one it treats as determined when a datum can only determine one of its entries.

entry is a NamedTuple (; term, index, site, dt, sefd, occupancy): the term's name, the entry's position in the term's SiteLookup, its site, the width of its integration time, the site's SEFD (Inf when the array table has no such site), and the number of stamps at which the site appears.

The key is (-dt, sefd, -occupancy): pin the entry with the longest integration time first, then the most sensitive site, then the site with the most data.

source
Comrade.SingleStokesGain Type
julia
SingleStokesGain(param_map)

Construct a gain term that is applicable to a single measured visibility. This is useful for pure stokes I modeling. The param_map is a function that maps the set of parameters to the gain term.

The arguments of param_map is typically a named tuple, where some of the elements are the parameter values needed to compute the gain. The return value of the param_map should be a single number or complex gain.

Example

julia
G = SingleStokesGain(x->exp(x.lg + xp.gp))
source
Comrade.JonesG Type
julia
JonesG(param_map)

Describes a gain Jones matrix with layout

g1 0 0 g2

where g1 and g2 are the gains for first and second feed of the telescope.

The param_map is a function that maps the set of parameters to the gain term. The arguments of param_map is typically a named tuple, where some of the elements are the parameter values needed to compute the gain. The return value of the param_map should be a two element tuple where the first element is the complex gain g1 and the second element is the complex gain g2.

Example

julia
G = JonesG() do
    g1 = exp(complex(x.lg1, x.gp1))
    g2 = g1*exp(complex(x.lgratio, x.gpratio))
    return g1, g2
end
source
Comrade.JonesD Type
julia
JonesD(param_map)

Describes a leakage Jones matrix with layout

1 d1 d2 1

where d1 and d2 are the d-terms for first and second feed of the telescope.

The param_map is a function that maps the set of parameters to the gain term. The arguments of param_map is typically a named tuple, where some of the elements are the parameter values needed to compute the gain. The return value of the param_map should be a two element tuple where the first element is the complex d-term d1 and the second element is the complex d-term d2.

Example

julia
D = JonesD() do
    d1 = complex(x.d1real, x.d1imag)
    d2 = complex(x.d2real, x.d2imag)
    return d1, d2
end
source
Comrade.JonesR Type
julia
JonesR(;add_fr=true)

The response Jones matrix. This is the reponse the telescope has to the incoming electric field, if the telescope was ideal. If add_fr=true then feed rotation are included.

This Jones matrix has no parameters so it doesn't accept a param_map. The add_fr argument is a boolean that specifies if feed rotation should be included.

source
Comrade.JonesF Type
julia
JonesF(;add_fr=true)

The feed rotation Jones matrix. This matrix describes the orientation of the feeds of the telescope.

This Jones matrix has no parameters so it doesn't accept a param_map. The add_fr argument is a boolean that specifies if feed rotation should be included.

source
Comrade.GenericJones Type
julia
GenericJones(param_map)

Construct a generic dense jones matrix with four parameterized elements.

The param_map is a function that maps the set of parameters to the gain term. The arguments of param_map is typically a named tuple, where some of the elements are the parameter values needed to compute the gain. The return value of the param_map should be a four element tuple where the elements are the entries of the jones matrix in column major order.

Example

julia
J = GenericJones() do
    return x.j11, x.j21, x.j12, x.j22
end
source
Comrade.JonesSandwich Type
julia
JonesSandwich([decomp_function=*,] matrices::AbstractJonesMatrix...)

Constructs a Jones matrix that is the results combining multiple Jones matrices together. The specific composition is determined by the decomp_function. For example if the decomp function is * then the matrices are multiplied together, if it is + then they are added.

Examples

julia
G = JonesG(x->(x.gR, x.gL)) # Gain matrix
D = JonesD(x->(x.dR, x.dL)) # leakage matrix
F = JonesF()                # Feed rotation matrix

J = JonesSandwich(*, G, D, F) # Construct the full Jones matrix as G*D*F

# Or if you want to include FR calibration
J = JonesSandwich(G, D, F) do g, d, f
    return adjoint(f)*g*d*f
end
source
Comrade.AbstractInstrumentModel Type
julia
abstract type AbstractInstrumentModel

The abstract instrument model. For a concrete implementation see IdealInstrumentModel and InstrumentModel.

Any subtype must implement the following methods

  • set_array(m::AbstractInstrumentModel, array::AbstractArrayConfiguration): Sets the array configuration for the instrument model m and returns the observed instrument model and prior.

  • apply_instrument(vis, m::AbstractInstrumentModel, x): Applies the instrument model m to the visibilities vis given the model parameters x.

A model that depends on the data also implements set_array(m, array, data), where data are the data products the posterior fits; the default ignores them.

source
Comrade.IdealInstrumentModel Type
julia
IdealInstrument(array::AbstractArrayConfiguration)

Constructs an ideal instrument that has no corruptions or feed rotation.

source
Comrade.InstrumentModel Type
julia
InstrumentModel(jones, prior; refbasis = CirBasis(), gaugefix = :error, selfcal = false)

Builds an instrument model using the jones matrix jones, with priors prior. The reference basis is refbasis and is used to define what the ideal basis is. Namely, the basis that you have the ideal visibilties to be represented in. For classical VLBI refbasis = CirBasis is a good default option, sinc the majority of the array uses circular feeds. For linear feed arrays like VGOS a user should switch to LinBasis, although failure to do so will not cause any errors, and is just a less efficient representation of the visibilities.

Arguments

  • jones : The jones matrix that represents the instrument. This is a function that takes in the parameters of the instrument and returns a jones matrix. See SingleStokesGain for a Stokes I example and JonesG or JonesD for polarized examples.

  • prior: A named tuple of ArrayPrior that specify what the priors are for each component used to construct the jones matrix using the function jones

Optional Arguments

  • refbasis: The reference basis used for the computation. The default is CirBasis() which are circular feeds.

  • gaugefix: What to do when the priors declared with gauge = :phase (see ArrayPrior) add up to a station phase the data cannot determine. :error (the default) throws, listing the entries that would fix it; :pin adds those entries to each affected prior's referencing scheme and logs one line per pin. The check uses Comrade.gauge_pins and runs when the model is given an array configuration.

  • selfcal: whether each station's common phase is written as φ = φ̂ + selfcal_residual, where φ̂ is the self-calibration solution. jones must then omit the station phase. See the section below.

Self-calibrated station phase

With selfcal = true, the model adds the instrument parameter selfcal_residual: one angle per station and data time, with prior WrappedNormal(0, 2π), whose density differs from the uniform density on the circle by a relative 2exp(-2π²) ≈ 5e-9 at most. In each group of stations that baselines link at one time and frequency, the station with the largest total data weight Σ 1/σ² is the reference, and its residual is fixed to 0. The priors of all other terms are those in prior.

φ̂ is the phase of the leading eigenvector of each group's Hermitian station matrix

julia
A_ij = Σ_p V_ij,p conj(M_ij,p) / σ²_ij,p      (i ≠ j)
A_ii = Σ_j Σ_p |M_ij,p|² / σ²_ij,p

where V and σ are the data and their noise, M the model visibilities without the station phase, and p runs over the correlation products. The off-diagonal part is the phase-dependent term of the Gaussian log-likelihood; the diagonal makes A positive semidefinite when the model matches the data. The eigenvector is approximated by A^32 e_r, 32 power iterations from the unit vector e_r of the reference station, so φ̂ is zero there.

φ̂ depends on the sky and the other instrument parameters but not on the residual, so the change of variables from φ to the residual has unit Jacobian. Because the prior is uniform on the circle, the prior of φ does not depend on φ̂, and the posterior is that of a uniform station phase. The residual is the station phase in coordinates that depend on the data, so compare or reuse draws through φ (selfcal_phase). Because φ̂ depends on the data, the model can only be used in a VLBIPosterior with one table of complex visibilities or coherencies.

Example

A Stokes I example is

julia
julia> G = SingleStokesGain(x->exp(complex(x.lg, x.pg)))
julia> intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
            pg = ArrayPrior(IIDSitePrior(ScanSeg(), DiagVonMises(0.0, inv(π^2))))
            )

julia> intm = InstrumentModel(G, intprior)

A standard polarized example is

julia
julia> G = JonesG() do
        gR = exp.(complex(x.lgr, x.gpr))
        gL = gr*exp.(complex(x.lgrat, x.gprat))
        return gR, gL
    end
julia> D = JonesD() do
        dR = complex.(x.dRre, x.dRim)
        dL = complex.(x.dLre, x.dLim)
        return gR, gL
    end
julia> R = JonesR()
julia> J = JonesSandwich(G, D, R)
julia> intprior = (lgr = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)),
                    gpr = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))),
                    lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1)),
                    gprat = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))),
                    dRre = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1)),
                    dRim = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1)),
                    dLre = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1)),
                    dLim = ArrayPrior(IIDSitePrior(TrackSeg(), VLBIGaussian(0.0, 0.1))
                    )
julia> intm = InstrumentModel(J, intprior)

which construct the gain matrix from R and ratios, and D is the small leakage matrix. JonesR is the response matrix that controls how the site responds to the ideal visibility in the reference basis.

source
Comrade.SiteArray Type
julia
SiteArray(data, times, frequencies, sites)

A SiteArray is an array of data that has a specified ordering of times, frequencies, and sites. Each data point is assigned a unique time, frequency, and site code. This allows for easy selection of data points based on these criteria and forms the base array for instrument modeling.

To select a subset of the data based on a specifid site, time and frequency you can use

julia
sarr[S=:ALMA, Ti=1:10, Fr=1:10]

which will grab the first 10 time and frequency points for the ALMA site.

Otherwise indexing into the array will return an element whose time, frequency, and site are the element of the times, frequencies, and sites arrays respectively.

source
Comrade.SiteLookup Type
julia
SiteLookup(s::SiteArray)

Construct a site lookup dictionary for a site array.

source
Comrade.forward_jones Function
julia
forward_jones(J::AbstractJonesMatrix, xs::NamedTuple{N})

Construct the forward model for the jones matrix model J with the parameters xs.

The xs is a named tuple where the keys are the parameter names and the values are SiteArrays with the parameter values. The return value is a SiteArray whose dimension is the largest of the elements of `xs, and whose elements are the jones matrices for the specific parameters.

source
Comrade.selfcal_phase Function
julia
selfcal_phase(post::VLBIPosterior, θ)

The station phase φ = φ̂ + residual at θ of a posterior whose instrument model self-calibrates the station phase (see InstrumentModel), as a SiteArray.

source

Posterior Constructions ​

Comrade.AbstractVLBIPosterior Type
julia
abstract type AbstractVLBIPosterior

An abstract VLBI posterior. See VLBIPosterior for a concrete implementation. This implements the DensityInterface and LogDensityProblem interfaces.

Default methods include:

  • logprior(d::AbstractVLBIPosterior, θ): Computes the log-prior of the posterior.

  • loglikelihood(d::AbstractVLBIPosterior, θ): Computes the log-likelihood of the posterior.

  • dimension(d::AbstractVLBIPosterior): Returns the dimension of the posterior.

  • skymodel(d::AbstractVLBIPosterior, θ): Returns the sky model of the posterior.

  • prior_sample(rng::AbstractRandom, d::AbstractVLBIPosterior, dims...): Samples from the prior of the posterior.

  • forward_model(d::AbstractVLBIPosterior, θ): Computes the forward model visibilities of the posterior.

  • dataproducts(d::AbstractVLBIPosterior): Returns the data products you are fitting as a tuple.

source
Comrade.logprior Function
julia
logprior(d::AbstractVLBIPosterior, θ)

Computes the log-prior of the posterior d with parameters θ.

source
Comrade.loglikelihood Function
julia
loglikelihood(d::AbstractVLBIPosterior, θ)

Computes the log-likelihood of the posterior d with parameters θ.

source
Comrade.dataproducts Function
julia
dataproducts(d::AbstractVLBIPosterior)

Returns the data products you are fitting as a tuple. The order of the tuple corresponds to the order of the dataproducts argument in VLBIPosterior.

source
Comrade.skymodel Function
julia
skymodel(d::AbstractVLBIPosterior)

Returns the sky model of the posterior d.

source
julia
skymodel(post::AbstractVLBIPosterior, θ)

Returns the sky model or image of a posterior using the parameter valuesθ

source
Comrade.instrumentmodel Method
julia
instrumentmodel(d::AbstractVLBIPosterior)

Returns the instrument model of the posterior d.

source
Comrade.instrumentmodel Method
julia
instrumentmodel(post::AbstractVLBIPosterior, θ)

Returns the instrument model of a posterior using the parameter values θ. The output will be a SiteArray of the Jones matrices for each site, time, and frequency.

source
Comrade.prior_sample Function
julia
prior_sample([rng::AbstractRandom], post::AbstractVLBIPosterior, [dims=1])

Returns sample from the prior distribution of the posterior. If dims is specified then multiple independent draws are returned with shape dims.

source
Comrade.likelihood Function
julia
likelihood(d::ConditionedLikelihood, μ)

Returns the likelihood of the model, with parameters μ. That is, we return the distribution of the data given the model parameters μ. Samples from this distribution are simulated data.

source
Comrade.VLBIPosterior Type
julia
VLBIPosterior(skymodel::SkyModel, instumentmodel::InstrumentModel, 
              dataproducts::EHTObservationTable...; 
              admode=set_runtime_activity(set_strong_zero(Reverse)))

Creates a VLBILikelihood using the skymodel its related metadata skymeta and the instrumentmodel and its metadata instumentmeta. The model is a function that converts from parameters θ to a Comrade AbstractModel which can be used to compute visibilitymap and a set of metadata that is used by model to compute the model.

To enable automatic differentiation, the admode keyword argument can be set to any EnzymeCore.Mode type of if no AD is desired then nothing. The default is Enzyme.set_runtime_activity(Enzyme.Reverse) which should work for essentially every problem. Note that runtime activity does have a perfomance cost, and as Enzyme and Comrade matures we expect this to not need runtime activity.

Warning

The model itself must be a two argument function where the first argument is the set of model parameters and the second is a container that holds all the additional information needed to construct the model. An example of this is when the model needs some precomputed cache to define the model.

Example

julia
dlcamp, dcphase = extract_table(obs, LogClosureAmplitude(), ClosurePhases())
array = arrayconfiguration(dlcamp)

function sky(θ, metadata)
    (; r, a) = θ
    m = stretched(ExtendedRing(a), r, r)
    return m
end

skyprior = (r = VLBIUniform(μas2rad(10.0), μas2rad(30.0)), a = VLBIUniform(1.0, 10.0))
g  = imagepixels(μas2rad(100.0), μas2rad(100.0), 256, 256)
skym = SkyModel(sky, skyprior, g)

G = SingleStokesGain(x->exp(complex(x.lg, x.pg)))
intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), VLBIGaussian(0.0, 0.1))),
            pg = ArrayPrior(IIDSitePrior(ScanSeg(), DiagVonMises(0.0, inv(π^2))))
            )
intmodel = InstrumentModel(G, intprior, array)

post = VLBIPosterior(skym, intmodel, dlcamp, dcphase)
source
Comrade.simulate_observation Function
julia
simulate_observation([rng::Random.AbstractRNG], post::VLBIPosterior, θ; add_thermal_noise=true)

Create a simulated observation using the posterior and its data post using the parameter values θ. In Bayesian terminology this is a draw from the posterior predictive distribution.

If add_thermal_noise is true then baseline based thermal noise is added. Otherwise, we just return the model visibilities.

source
Comrade.residuals Function
julia
residuals(post::VLBIPosterior, θ)

Compute the residuals for each data product in post using the parameter values θ. The resturn objects are EHTObservationTables, where the measurements are the residuals.

source
Comrade.TransformedVLBIPosterior Type
julia
struct TransformedVLBIPosterior{P<:VLBIPosterior, T} <: Comrade.AbstractVLBIPosterior

A transformed version of a VLBIPosterior object. This is an internal type that an end user shouldn't have to directly construct. To construct a transformed posterior see the asflat, ascube.

source
Comrade.ImgNormalData Type
julia
ImgNormalData(reduction, measurement, noise)

Container for image-domain data assuming a normal likelihood with a measurement and noise.

source
TransformVariables.transform Method
julia
transform(posterior::TransformedVLBIPosterior, x)

Transforms the value x from the transformed space (e.g. unit hypercube if using ascube) to parameter space which is usually encoded as a NamedTuple.

For the inverse transform see inverse

source
TransformVariables.inverse Method
julia
inverse(posterior::TransformedVLBIPosterior, x)

Transforms the value x from parameter space to the transformed space (e.g. unit hypercube if using ascube).

For the forward transform see transform

source
VLBIImagePriors.ascube Method
julia
ascube(post::VLBIPosterior)

Construct a flattened version of the posterior where the parameters are transformed to live in (0, 1), i.e. the unit hypercube.

This returns a TransformedVLBIPosterior that obeys the DensityInterface and can be evaluated in the usual manner, i.e. logdensityof. Note that the transformed posterior automatically includes the terms log-jacobian terms of the transformation.

Example

julia
julia> tpost = ascube(post)
julia> x0 = prior_sample(tpost)
julia> logdensityof(tpost, x0)

Notes

This is the transform that should be used if using typical NestedSampling methods, i.e. ComradeNested. For the transformation to unconstrained space see asflat

source
VLBIImagePriors.asflat Method
julia
asflat(post::VLBIPosterior)

Construct a flattened version of the posterior where the parameters are transformed to live in (-∞, ∞).

This returns a TransformedVLBIPosterior that obeys the DensityInterface and can be evaluated in the usual manner, i.e. logdensityof. Note that the transformed posterior automatically includes the terms log-jacobian terms of the transformation.

Example

julia
julia> tpost = asflat(post)
julia> x0 = prior_sample(tpost)
julia> logdensityof(tpost, x0)

Notes

This is the transform that should be used if using typical MCMC methods, i.e. NUTS. For the transformation to the unit hypercube see ascube

source
ProbabilityTransports.transport_to Method
julia
transport_to(post::VLBIPosterior, space)

Build a TransformedVLBIPosterior whose parameters live in the latent space (TVFlat() for unconstrained ℝⁿ, StdUniform() for the unit hypercube). The asflat / ascube helpers call this with the respective spaces.

Also accepts an already-transformed posterior: it is returned unchanged when it already lives in space, and rebuilt from the underlying VLBIPosterior otherwise — an explicit space request is always honored.

The space argument is typed (rather than left as ::Any) so these methods stay strictly more specific than ProbabilityTransports' own transport_to(dist, ::AbstractStdDist) / transport_to(dist, ::TVFlat) and don't collide with them by ambiguity.

source

Inference ​

Comrade.comrade_opt Function
julia
comrade_opt(post::VLBIPosterior, opt, adtype=nothing, args...; initial_params=nothing, kwargs...)

Optimize the posterior post using the opt optimizer. The adtype specifies the automatic differentiation. The args/kwargs are forwarded to `the specific optimization package.

Warning

This function won't have any methods until the optimization package is loaded. We recommend loading the Optimization.jl package.

source
Comrade.MemoryStore Type
julia
MemoryStore

Stored the HMC samplers in memory or ram.

source
Comrade.DiskStore Type
julia
DiskStore(;name::String = "Results", stride::Int = 100, callback = default_disk_callback)

Type that specifies to save the samples results to disk.

Fields

  • name: Path of the directory where the results will be saved. If the path does not exist it will be automatically created.

  • stride: The output stride, i.e. every stride steps the MCMC output will be dumped to disk.

  • callback: A callback info -> ... run after every batch is written to disk. Shared by the AdvancedHMC and ReactantNUTS disk-sampling paths; see default_disk_callback for the fields of info. The default logs a one-line summary per batch.

source
Comrade.load_samples Function
julia
load_samples(out::DiskOutput, indices::Union{Base.Colon, UnitRange, StepRange}=Base.Colon(); table=:both)
load_samples(out::String, indices::Union{Base.Colon, UnitRange, StepRange}=Base.Colon(); table=:both)

The the results from a HMC run saved to disk. To read in the output the user can either pass the resulting out object, or the path to the directory that the results were saved, i.e. the path specified in DiskStore.

Arguments

  • out::Union{String, DiskOutput}: If out is a string is must point to the direct that the DiskStore pointed to. Otherwise it is what is directly returned from sample.

  • indices: The indices of the that you want to load into memory. The default is to load the entire table. Only the sample files holding these draws are read, each once.

Keyword Arguments

  • table: The table to read: :both (the default) gives the MCMC chain with its sampler statistics, :samples the chain alone, and :stats the sampler statistics alone, e.g., the log density of each sample and tree statistics.
source
Comrade.PosteriorSamples Type

PosteriorSamples(chain, stats; metadata=Dict(:sampler=>:unknown))

This is the default sampler output from Comrade MCMC extensions. The object contains the posterior samples, the sampler statistics, and metadata about the sampler used.

Indexing this array behaves like indexing the samples organized as a StructArray. By default all NamedTuples and Tuples are unwrapped.

The samples of several chains are a draws × chains PosteriorSamples, made with hcat of the chains' samples (reduce(hcat, chains) for a vector of them). Then s[i, k] is draw i of chain k, s[:, k] is chain k, each parameter (e.g. s.sky.σ) is a draws × chains array, and linear indexing, iteration and Comrade.rmap run over the draws of all chains. With MCMCDiagnosticTools loaded, ess, rhat and mcse of a PosteriorSamples give the diagnostic of every parameter, nested like the output of rmap.

To access the samples use postsamples and to access the sampler statistics use samplerstats, or the acesss sampler specific information use samplerinfo.

To recursively map a function over the samples the unexported Comrade.rmap.

source
Comrade.postsamples Function
julia
postsamples(s::PosteriorSamples)

Return the samples from the PosteriorSamples object s

source
Comrade.samplerstats Function
julia
samplerstats(s::PosteriorSamples)

Return the sampler statistics from the PosteriorSamples object s

source
Comrade.samplerinfo Function
julia
samplerinfo(s::PosteriorSamples)

Return the metadata from the PosteriorSamples object s.

source
Comrade.resample_equal Function
julia
resample_equal(post::PosteriorSamples, nsamples::Int)

Resample the posterior samples so you have nsamples of equal weight. In order for this method to be applicable a :weights field must be present in the sampler statistics and the weight must correspond to the probability weights of the samples.

source
Comrade.residual Function
julia
residual(post::AbstractVLBIPosterior, p)

Plots the normalized residuals for the posterior post given the parameters p.

source
Comrade.residual_data Function
julia
residual_data(vis::AbstractArray, data::EHTObservationTable)

Compute the residuals for the model visibilities vis and the data data. The residuals are not normalized and the returned object is an EHTObservationTable.

source
Comrade.chi2 Function
julia
chi2(post::AbstractVLBIPosterior, p; reduce=false)

Returns a tuple of the chi-squared values for each data product in the posterior post given the parameters p. Note that the chi-square is not reduced. If you want to reduce it by dividing by the number of data points then set reduce=true.

Note that reduce doesn't take into account the number of model parameters which for non-linear models is difficult to define globally.

source
Comrade.ReactantNUTS Type
julia
ReactantNUTS(; kwargs...)

Experimental NUTS sampler that runs through Reactant's ProbProg mcmc_logpdf. Like AdvancedHMC's NUTS, all sampler configuration lives in the struct.

Warmup runs Stan's windowed adaptation: fixed init_buffer = 75, term_buffer = 50, base window 25, doubling, with the same 15%/10% fallback when the buffers don't fit — one continuous Nesterov dual-averaging step-size process (γ=0.05, t₀=10, κ=0.75, target accept 0.8) plus a Welford diagonal mass matrix updated at window ends.

Warmup is executed in chunks, by default the same size used for sampling (saveto.stride for a DiskStore, otherwise the chunk_size passed to sample; override with sample's warmup_chunk). The ProbProg backend anchors the windowed schedule to the global warmup length (total_warmup/warmup_offset), so chunking is bit-identical to one fused warmup whatever the chunk length. When sampling to a DiskStore, the adaptation state is checkpointed every checkpoint_every warmup steps, so an interrupted warmup can be resumed with restart = true, and one draw per chunk is logged to <name>/warmup for inspection.

Keyword arguments

n_adapts = 1000 total warmup steps; the internal Stan windowed schedule is derived from this by the ProbProg backend. init_step_size = 0.01 initial leapfrog step size; seeds dual averaging (see note) max_tree_depth = 10 max_delta_energy = 1000.0 strong_zero = true turn 0_Inf / 0_NaN in the gradient into 0; REQUIRED for stiff image models or every proposal NaN-rejects and the chain freezes – keep true unless you know better.

Differences from AdvancedHMC's NUTS

Two pieces of AdvancedHMC's setup are not (yet) reachable through ProbProg:

  • No find_reasonable_step_size. AdvancedHMC seeds dual averaging from find_good_stepsize; ProbProg has no such heuristic and instead anchors dual averaging at init_step_size (its prox-center is log(10*init_step_size)). Pick init_step_size in the right ballpark for your problem.

  • Target acceptance is fixed at 0.8 in the backend (matching NUTS(0.8)'s default) and is not currently configurable.

source
Comrade.default_disk_callback Function
julia
default_disk_callback(info) -> NamedTuple

Default per-batch callback used by DiskStore. It logs one summary line for the batch just written to disk and returns a small NamedTuple (collected into the sampler's sample_history metadata).

info is a NamedTuple describing the batch. The fields guaranteed by every sampler backend (AdvancedHMC's sample and ReactantNUTS) — i.e. the ones this default callback uses — are:

  • round::Int : batch index (1-based)

  • nrounds::Int : total number of batches

  • num_samples::Int : draws in this batch

  • time::Real : wall-clock seconds spent drawing this batch

  • numerical_error::Vector{Bool} : per-draw divergence flags

  • params : last draw, as the constrained (; sky[, instrument]) NamedTuple

  • step_size::Real : current leapfrog step size

  • extras::NamedTuple : backend-specific data — present on every backend, but its CONTENTS are NOT guaranteed to match across backends

info.extras is where each backend exposes whatever it can cheaply provide (AdvancedHMC's samplerstats, ReactantNUTS's host-side MCMCState view, ...). A callback that reaches into extras is therefore backend-specific; the fields are documented on the corresponding sample method (AdvancedHMC's sample and ReactantNUTS's sample).

A custom callback can do anything with info, e.g. plot the current draw:

julia
DiskStore(; name = "Results", callback = info -> display(imageviz(intensitymap(skymodel(post, info.params.sky), g))))
source

Moves ​

Comrade.AbstractMove Type
julia
AbstractMove

Supertype of Metropolis–Hastings moves on the latent point of a posterior. A move implements propose, step_distribution and move_name, and may implement likelihood_invariant, step_label, component_label, context_data and traceable.

propose must not modify x. Proposing with the reverse step from the proposal must return to x with the opposite log-determinant. check_move tests a move.

source
Comrade.MoveSet Type
julia
MoveSet(post::VLBIPosterior, moves; space = nothing, θ0 = nothing, rounds = 1,
        initial_scale = 0.1, target_accept = 0.45, output = nothing,
        rng = Random.default_rng())

A set of AbstractMoves in the latent space space of post, with their step scales and acceptance statistics.

  • rounds: proposals of each move per call.

  • initial_scale: the starting step scale of moves with a NormalStep or ComponentwiseNormalStep; ignored for other moves, whose entries may be nothing.

  • target_accept: the acceptance rate the step scales adapt toward during warmup.

  • θ0: a parameter point. If given, construction checks there that every likelihood_invariant move leaves the likelihood unchanged, and that the components of every componentwise move are independent, with steps drawn from rng.

  • output: a file path. If given, move_summary is serialized there after every call.

rounds, initial_scale and target_accept take one number for all moves or a vector with one value per move. A componentwise move must be likelihood_invariant.

julia
ms(x, tpost, phase, rng; step = nothing, total = nothing, report = true) -> x′

Run the moves on the latent point x of tpost, which must use the same space, in the sampler phase phase (:warmup or :sampling). Call it between the steps of any MCMC sampler. Returns x itself if no proposal was accepted. A proposal with a NaN acceptance ratio is rejected. During warmup the step scales adapt, and with report the call logs the acceptance rates and step scales since the last log, labeled with step out of total. Throws if an accepted likelihood_invariant move changed the likelihood.

julia
ms(state, tpost, info::NamedTuple, rng) -> state

The between_chunks hook of ReactantNUTS: moves state.position, reading phase, step, total and report from info.

The manual page "Writing and using moves" describes how proposals are accepted and run on the device.

source
Comrade.NormalStep Type
julia
NormalStep()

A scalar step u ~ N(0, τ²), where τ is the move's step scale. The reverse step is −u.

source
Comrade.ComponentwiseNormalStep Type
julia
ComponentwiseNormalStep(n)

A vector step of n components, uᵢ ~ N(0, τᵢ²), each with its own step scale. The reverse step is −u. Each component is accepted or rejected on its own, so the move must implement component_coords, the components must change disjoint coordinates, and component i of the proposal must depend only on uᵢ and the coordinates of component i. propose returns one log-determinant per component.

source
Comrade.SymmetricDiscreteStep Type
julia
SymmetricDiscreteStep(draw, reverse)

A step drawn by draw(rng). reverse(step) is the reverse step, which draw must return with the same probability as step.

source
Comrade.CompensatedMove Type
julia
CompensatedMove(name, view, shift, block, compensate; index = 1,
                logdet = (vb, x, x′, ctx) -> 0.0, likelihood_invariant = false,
                context = (;), traceable = false)

A NormalStep move that adds the step to the index-th latent coordinate of the parameter at path shift, and sets the parameter at path block to compensate(vb, x, x′, ctx). Here vb is the value of the block at x, and x′ is x with only the shift applied. compensate keeps fixed something the likelihood depends on, such as a product of parameters.

  • logdet(vb, x, x′, ctx): log|det ∂x′/∂x| of the whole move, given the complete proposal x′.

  • likelihood_invariant: see likelihood_invariant.

  • context: a NamedTuple of arrays added to ctx.

  • traceable: see traceable. compensate and logdet must then read ctx.view and the context arrays from ctx.

Throws if shift or block is not in view, or if the shifted coordinate lies in block.

source
Comrade.ChainHyperMove Type
julia
ChainHyperMove

A move of one fitted hyperparameter field (e.g. σ or τ) of a Gauss–Markov chain term, with one ComponentwiseNormalStep component per site; build it with chain_hyper_moves. Component i changes the hyperparameter h of site i and recomputes the site's whitened chain coordinates z, so the chain values g, and with them the likelihood, stay the same. The log-determinant log p(g | h′) − log p(g | h) + (‖z′‖² − ‖z‖²)/2 makes the acceptance that of the hyperparameter given the chain values.

source
Comrade.PhaseSheetMove Type
julia
PhaseSheetMove(post::VLBIPosterior, terms; space = nothing)

A move that adds ±2π to one site's real-line Gauss–Markov phase chain, from a free point up to the next point the model fixes (or the end of the chain). terms names the instrument phase terms, e.g. (:gp,). The likelihood depends on a phase only through e^{iφ}, so it is unchanged; the chain prior changes. A step picks a term, site, starting point and sign uniformly, and the opposite sign reverses it. The shift is a constant in latent space, so the log-determinant is 0.

Throws if a term is not an instrument parameter, has no free point, or is a wrapped chain.

source
Comrade.chain_hyper_moves Function
julia
chain_hyper_moves(view::CoordinateView, path) -> Vector{ChainHyperMove}

One ChainHyperMove per fitted hyperparameter field of the Gauss–Markov chain term at path (e.g. (:instrument, :lg1)), named "chain_hyper[term.field]".

Throws if path is not a chain term with fitted hyperparameters, or if a site's chain is grouped or wrapped non-centered: the move needs the chain density of each site.

source
Comrade.flux_gain_move Function
julia
flux_gain_move(view::CoordinateView; flux, gains, power = 2)

A CompensatedMove that trades the scalar flux at path flux against the gain log-amplitudes at path gains (a Gauss–Markov chain or per-site constants). A step changes the flux from F to F′ and subtracts log(F′/F) / power from every gain. The likelihood is unchanged when each visibility depends on these only through F · exp(power · g), e.g. power = 2 when both stations of a baseline share the log-amplitude g. The log-determinant is 0, which holds when the gain shift is a constant in latent space: a whitened chain with fixed hyperparameters, or Gaussian per-site constants.

source
Comrade.CoordinateView Type
julia
CoordinateView(post::VLBIPosterior, space = nothing)

The parameters of post in a latent space: the flat space of asflat (space = nothing) or ProbabilityTransports.StdNormal(). A parameter is named by its path, such as (:sky, :σa) or (:instrument, :lg1), and is read and written with coords, value and latent. view.tbase is the transformed posterior.

source
Comrade.check_move Function
julia
check_move(move, post, θs; space = nothing, rng = Random.default_rng(), τ = 0.1,
           rtol = 1e-8, h = 1e-6, logdet_atol = 1e-5, nprior = 0, nsteps = 5,
           ncoords = 5, ks_alpha = 1e-3, ngeweke = 0, geweke_steps = 10) -> NamedTuple

Test move on post in the latent space space, at each parameter point in θs, with steps at scale τ. Throws an error naming the first failed check:

  • reversal: the reverse step from the proposal returns to the start, with the opposite log-determinant;

  • log-determinant: logdet matches a finite-difference Jacobian (step h);

  • invariance: if likelihood_invariant(move), the log-likelihood is unchanged;

  • components: for a ComponentwiseNormalStep move, the components are independent: disjoint coordinates, and component i of the step changes only component i.

Two optional tests compare samples with two-sample Kolmogorov–Smirnov tests at level ks_alpha, on up to ncoords of the coordinates the move changes:

  • nprior > 0: nprior prior draws each take nsteps steps of the move targeting the prior, and are compared with fresh prior draws.

  • ngeweke > 0: Geweke's (2004) joint-distribution test of the kernel MoveSet runs. Each of ngeweke replicates draws parameters from the prior and data from the model (simulate_observation), then alternates geweke_steps move steps with new data. The results are compared with fresh draws of parameters and data, also along the move's direction and on the log-likelihood. Arrays in the move context are not rebuilt for the new data.

Returns the largest reversal error, log-determinant error and log-likelihood change; the KS statistics ks and geweke with their critical values; and geweke_accept, the acceptance rate in the Geweke test. A move that rarely moves passes the Geweke test untested, so check that geweke_accept is not near zero.

source
Comrade.check_moveset Function
julia
check_moveset(ms::MoveSet; ngeweke, geweke_steps = 10, rng = Random.default_rng(),
              ncoords = 5, ks_alpha = 1e-3) -> NamedTuple

Geweke's joint-distribution test (see check_move) of one call of ms at its current step scales, with ngeweke replicates of geweke_steps steps. It compares up to ncoords of the coordinates each move changes, and the position along each move's direction. Returns geweke, geweke_critical and geweke_accept (one rate per move). Throws if a KS statistic exceeds its critical value.

source
Comrade.move_summary Function
julia
move_summary(ms::MoveSet) -> Vector{NamedTuple}

Statistics for each move: name; the current step scale τ (a vector for a componentwise move, nothing for a move without one); rounds; the proposed and accepted counts in warmup and sampling; and bystep, the counts per step_label (per component_label for a componentwise move) as label => [proposed, accepted].

source
Comrade.move_seconds Function
julia
move_seconds(ms::MoveSet) -> Float64

The total wall time of the calls of ms so far, in seconds.

source

The protocol a move implements, and the helpers for writing and checking one:

Comrade.propose Function
julia
propose(move, x, step, ctx)

The proposal of move from the latent point x: (x′, logdet) with logdet = log|det ∂x′/∂x| (one per component for a ComponentwiseNormalStep). ctx comes from move_context.

source
Comrade.step_distribution Function
julia
step_distribution(move)

The step distribution of move: a NormalStep, ComponentwiseNormalStep or SymmetricDiscreteStep.

source
Comrade.move_name Function
julia
move_name(move) -> String

The name of move in logs and statistics.

source
Comrade.likelihood_invariant Function
julia
likelihood_invariant(move) -> Bool

Whether move leaves the likelihood unchanged (default false). MoveSet accepts such a move from the change in the prior alone, without evaluating the likelihood, then checks that the likelihood did not change.

source
Comrade.context_data Function
julia
context_data(move, view) -> NamedTuple

Arrays that move reads from ctx in propose (default none). update_context! replaces their contents.

source
Comrade.traceable Function
julia
traceable(move) -> Bool

Whether propose for move compiles with Reactant (default false). On a device position, MoveSet runs such a proposal on the device, with ctx.view over the device posterior. A traceable move has a NormalStep and reads every array it needs from ctx rather than capturing it.

source
Comrade.move_context Function
julia
move_context(view::CoordinateView, moves) -> NamedTuple

The ctx argument of propose: (; view, ...) merged with the context_data of every move.

source
Comrade.update_context! Function
julia
update_context!(ms::MoveSet, key::Symbol, value) -> MoveSet

Copy value into the move-context array key of ms (same size). The moves read the new value from the next call on; a compiled device program is updated without recompiling.

source
Comrade.step_label Function
julia
step_label(move, step)

A key under which move_summary counts the proposals and acceptances of step, or nothing (default) to count only the move's totals.

source
Comrade.component_coords Function
julia
component_coords(move) -> Vector{Vector{Int}}

The latent coordinates each component of a ComponentwiseNormalStep move changes.

source
Comrade.component_label Function
julia
component_label(move, i)

The name of component i of a ComponentwiseNormalStep move in its statistics (default i).

source
Comrade.node Function
julia
node(view::CoordinateView, path) -> transform node

The transform node of the parameter at path.

source
Comrade.coords Function
julia
coords(view::CoordinateView, path) -> UnitRange

The latent coordinates of the parameter at path.

source
Comrade.value Function
julia
value(view::CoordinateView, x, path)

The value of the parameter at path at the latent point x.

source
Comrade.latent Function
julia
latent(view::CoordinateView, path, v) -> Vector

The latent coordinates of the parameter at path that give it the value v.

source

Misc ​

Comrade.site_tuple Function
julia
site_tuple(sites, default; reference=nothing kwargs...)
site_tuple(obs::AbstractObservationTable, default; reference=nothing, kwargs...)
site_tuple(obs::AbstractArrayConfiguration, default; reference=nothing, kwargs...)

Convienence function that will construct a NamedTuple of objects whose names are the sites in the observation obs or explicitly in the argument sites. The NamedTuple will be filled with default if no kwargs are defined otherwise each kwarg (key, value) pair denotes a sites and value pair.

Optionally the user can specify a reference sites that will be dropped from the tuple. This is useful for selecting a reference sites for gain phases

Examples

julia
julia> sites = (:AA, :AP, :LM, :PV)
julia> site_tuple(sites, ScanSeg())
(AA = ScanSeg(), AP = ScanSeg(), LM = ScanSeg(), PV = ScanSeg())
julia> site_tuple(sites, ScanSeg(); AA = FixedSeg(1.0))
(AA = FixedSeg(1.0), AP = ScanSeg(), LM = ScanSeg(), PV = ScanSeg())
julia> site_tuple(sites, ScanSeg(); AA = FixedSeg(1.0), PV = TrackSeg())
(AA = FixedSeg(1.0), AP = ScanSeg(), LM = ScanSeg(), PV = TrackSeg())
julia> site_tuple(sites, VLBIGaussian(0.0, 0.1); reference=:AA, LM = VLBIGaussian(0.0, 1.0))
(AP = VLBIGaussian(0.0, 0.1), LM = VLBIGaussian(0.0, 1.0), PV = VLBIGaussian(0.0, 0.1))
source
Comrade.dirty_image Function
julia
dirty_image(fov::Real, npix::Int, obs::EHTObservation{<:EHTVisibilityDatum}) where T

Computes the dirty image of the complex visibilities assuming a field of view of fov and number of pixels npix using the complex visibilities found in the observation obs.

The dirty image is the inverse Fourier transform of the measured visibilties assuming every other visibility is zero.

source
Comrade.dirty_beam Function
julia
dirty_beam(fov::Real, npix::Int, obs::EHTObservation{<:EHTVisibilityDatum})

Computes the dirty beam of the complex visibilities assuming a field of view of fov and number of pixels npix using baseline coverage found in obs.

The dirty beam is the inverse Fourier transform of the (u,v) coverage assuming every visibility is unity and everywhere else is zero.

source
Comrade.beamsize Function
julia
beamsize(ac::AbstractArrayConfiguration)

Calculate the approximate beam size of the array ac as the inverse of the longest baseline distance.

source
julia
beamsize(obs::AbstractObservationTable)

Calculate the approximate beam size of the observation obs as the inverse of the longest baseline distance.

source
Comrade.apply_fluctuations Function
julia
apply_fluctuations(f, mimg::IntensityMap, δ::AbstractArray)

Apply multiplicative fluctuations to an image mimg with fluctuations δ. The function f is applied to the fluctuations and then the the transfored δ are multiplicatively applied to the image.

source
Comrade.apply_fluctuations! Function
julia
apply_fluctuations!(t::UnitFluxMap{F}, out::IntensityMap, mimg::IntensityMap, δ::AbstractArray)

Apply multiplicative fluctuations to an image mimg with fluctuations δ using the transformation t, by broadcasting the function F over the array δ and then normalizing the resulting array to have unit flux.

source
Comrade.UnitFluxMap Type
julia
UnitFluxMap(f)

A transformation that broadcasts a function f over an array, and the normalizes the resulting array to have unit flux. This is used with the function apply_fluctuations while imaging.

source
Comrade.corr_image_prior Function
julia
corr_image_prior(grid::AbstractRectiGrid, corr_length::Real; base=GMRF, order=1, lower=1.0, upper=2*max(size(grid)...))
corr_image_prior(grid::AbstractRectiGrid, obs::EHTObservationTable; base=GMRF, order=1, lower=1.0, upper=2*max(size(grid)...))

Construct a correlated image prior, for the image with grid grid, and using the observation dvis. The correlation will be a Markov Random Field (MRF) of order order, with the base distribution base. For base you can choose any of the Markov random fields defined in VLBIImagePriors, the default is GMRF which is a Gaussian MRF.

As part the prior will be a hierarchical prior with the correlation length as a hyperparameter. By default the correlation parameter uses a first order inverse gamma distribution for its prior. The frac_below_beam parameter is the fraction of the correlation prior mass that is below the beam size of the observation dvis. The lower and upper parameters are the lower and upper bounds of the correlation length, we don't let the correlation length to be too small or large for numerical reasons.

Arguments

  • grid::AbstractRectiGrid: The grid of the image to be reconstructed.

  • corr_length: The correlation length of the MRF. If this is an EHTObservationTable then the corr_length will be the approximate beam size of the observation.

Keyword Arguments

  • base: The base distribution of the MRF. Options include GMRF, EMRF, and CMRF

  • order: The order of the MRF. Default is first order

  • frac_below_beam: The fraction of the correlation prior mass that is below the beam size of the observation dvis. the default is 0.01 which means only 1% of the log-image correlation length is below the beam size.

  • lower: The lower bound of the correlation length. Default is 1.0

  • upper: The upper bound of the correlation length. Default is 2*max(size(grid)...)

Warn

An order > 2 will be slow since we switch to a sparse matrix representation of the MRF.

source
Comrade.rmap Function
julia
rmap(f, x::PosteriorSamples)

Recursively map a function f over the elements of x. For instance to compute the mean of all fields you can do Comrade.rmap(mean, chain)

source

Internal (Not Public API) ​

Comrade.build_datum Function
julia
build_datum(data::AbstractObservationTable, i::Int)

Build the datum F for a given observation table data. This is an internal method that users shouldn't have to deal with directly unless they are implementing a new AbstractObservationTable.

source
Comrade.ObservedSkyModel Type
julia
ObservedSkyModel(sky::AbstractSkyModel, array::AbstractArrayConfiguration)

Constructs a sky model that has been theoretically observed by an array with configuration array. Note that this is not a public facing method and is used internally to construct the observed sky model for use in VLBIPosterior. Users should construct a SkyModel and pass that to a VLBIPosterior object instead.

source
Comrade.forward_model_map Function
julia
forward_model(d::AbstractVLBIPosterior, θ)

Computes the forward model visibilities of the posterior d with parameters θ. Note these are the complex visiblities or the coherency matrices, not the actual data products observed.

source
Comrade.fixed_indices Function
julia
fixed_indices(d)

The entry indices an observed instrument prior holds constant: the pins of its referencing scheme and its init pins.

source

Data Internals (Internal) ​

Comrade.extract_amp Function
julia
extract_amp(obs; kwargs...)

Extracts the visibility amplitudes from an obs. This is an internal method for dispatch. Only use this if interfacing Comrade with a new data type.

source
Comrade.extract_cphase Function
julia
extract_cphase(obs; kwargs...)

Extracts the closure phases from an obs. This is an internal method for dispatch. Only use this if interfacing Comrade with a new data type.

source
Comrade.extract_lcamp Function
julia
extract_lcamp(obs; kwargs...)

Extracts the log-closure amplitudes from an obs. This is an internal method for dispatch. Only use this if interfacing Comrade with a new data type.

source
Comrade.extract_vis Function
julia
extract_vis(obs; kwargs...)

Extracts the stokes I complex visibilities from an obs. This is an internal method for dispatch. Only use this if interfacing Comrade with a new data type.

source
Comrade.extract_coherency Function
julia
extract_coherency(obs; kwargs...)

Extracts the full coherency matrix from an observation. This is an internal method for dispatch. Only use this if interfacing Comrade with a new data type.

source
Comrade.closure_designmat Function
julia
closure_designmat(type::Symbol, closures, scanvis::Scan)

Dispatch to closurephase_designmat or closureamp_designmat based on type (:cphase or :lcamp).

source
Comrade.elevation_parallactic Function
julia
elevation_parallactic(antenna_xyz, ra_deg, dec_deg, jd)
    -> (elevation_rad, parallactic_rad)

Compute the source elevation and parallactic angle at a single antenna and time.

Arguments:

  • antenna_xyz: ECEF coordinates of the antenna in metres (3-tuple, vector, or SVector{3}).

  • ra_deg, dec_deg: J2000 right ascension and declination of the source in degrees.

  • jd: Julian date (UT) of the observation.

Returns the elevation angle (radians, above horizon) and parallactic angle (radians, atan2 form) at that antenna.

Uses AstroLib.ct2lst for local mean sidereal time. The parallactic angle uses the standard atan2 form

julia
p = atan2(sin(HA) cos(lat), sin(lat) cos(δ) - cos(lat) sin(δ) cos(HA))

which matches the convention used by both eht-imaging and Astropy.

source
Comrade.ecef_to_geodetic Function
julia
ecef_to_geodetic(x, y, z) -> (lat, lon, h)

Convert WGS84 Earth-centered Earth-fixed coordinates (x, y, z) in metres to geodetic latitude (rad), longitude (rad), and ellipsoidal height (m).

Uses the closed-form Bowring formula. Result is accurate to better than a millimetre for any point on or near the Earth's surface.

source
Comrade.build_closure_config Function
julia
build_closure_config(dvis::EHTObservationTable, closures::StructArray;
                     type::Symbol, count::String="min")

Build a ClosureConfig from a visibility observation table dvis and a closure-spec table closures whose rows have at least the fields (T, F, noise, baseline) (with baseline a NTuple{3} for closure phases or NTuple{4} for log-closure amplitudes).

Arguments

  • type: :cphase or :lcamp

  • count: "min" (compute a minimal closure set per scan) or "max" (use the supplied closures as-is).

source

Plotting ​

Warning A user must first load Makie or a Makie backend, e.g., CairoMakie to use this functionality

Comrade.plotfields Function
julia
plotfields(obsdata::EHTObservationTable, field1, field2; 
            site1=nothing, site2=nothing, axis_kwargs=(;), legend_kwargs=(;), scatter_kwargs=(;))

Plots two data fields against each other.

Arguments

  • obsdata : EHTObservationTable containing the data to plot (closure quantities not supported yet)

  • field1 and field2 : The fields to plot. field1 - x axis, field2 - y axis If field1 or field2 is a function it will apply to datatable(obsdata) to get the value If field1 or field2 is a symbol, it will look for a predefined function:

    • :U - baseline u coordinate

    • :V - baseline v coordinate

    • :Ti - time

    • :Fr - frequency

    • :measure - measurement

    • :noise - noise

    • :amp - visibility amplitude

    • :phase - visibility phase

    • :uvdist - projected baseline length

    • :snr - signal to noise ratio

    • :res - normalized residual visibilities (only if obsdata contains the residuals)

  • site1 and site2 : Keywords for the sites forming the baseline being plotted, e.g. :ALMA, :APEX.

  • axis_kwargs : Keyword arguments for each subplot's Axis.

  • legend_kwargs : Keyword arguments passed to the figure legend.

  • scatter_kwargs : Keyword arguments passed to scatter! in each subplot.

source
Comrade.plotfields! Function
julia
plotfields!(f::Makie.GridPosition, obsdata::Comrade.EHTObservationTable, field1, field2;
            kwargs...)

Plots the two data fields in the grid position f. For more information about the arguments see the plotfields function.

source
Comrade.axisfields Function
julia
axisfields(fig, obsdata::EHTObservationField, field1, field2;
            legend=true, conjugate=true, axis_kwargs=(;), legend_kwargs=(;), scatter_kwargs=(;))

Plots two data fields from obsdata against each other on fig, returns a Makie Axis which can be used to configure subplots.

! warning This function is deprecated, use plotfields! instead.

source
Comrade.plotcaltable Function
julia
plotcaltable(gt...; width=150, height=125, layout=nothing, markers=nothing, labels=nothing, 
                    axis_kwargs=(;), legend_kwargs=(;), figure_kwargs=(;), scatter_kwargs=(;))

Automatically generate a grid of subplots plotting the Comrade.CalTable information. Each subplot corresponds to a different station in the array.

Argments

  • gt : The CalTables to plot.

Keyword Arguments

  • width : Subplot width

  • height : Subplot height

  • layout : Subplot layout (optional).

  • axis_kwargs : Keyword arguments for each subplot's Axis.

  • legend_kwargs : Keyword arguments passed to the figure legend.

  • figure_kwargs : Keyword arguments passed to Figure().

  • scatter_kwargs : Keyword arguments passed to scatter! in each subplot.

source
Comrade.baselineplot Function
julia
baselineplot(data, bl, fieldx, fieldy; kwargs...)
baselineplot(data, fieldx, fieldy; kwargs...)

Plots the baseline bl with fieldx on the x axis and fieldy on the y axis.

If bl is a Colon or bl is not specified all baselines are plotted.

Arguments

  • data : The data to plot

  • field1 and field2 : The fields to plot. field1 - x axis, field2 - y axis If field1 or field2 is a function, baselineplot will broadcast them over datatable(obsdata) to get the x and y values as the outputs from field1.(datatable(obsdata)) and field2.(datatable(obsdata)). If field1 or field2 is a symbol, it will look for a predefined function:

    • :U - baseline u coordinate

    • :V - baseline v coordinate

    • :Ti - time

    • :Fr - frequency

    • :measurement - measurement

    • :measwnoise - measurement and noise

    • :noise - noise

    • :amp - visibility amplitude

    • :phase - visibility phase

    • :uvdist - projected baseline length

    • :snr - signal to noise ratio

    • :res - normalized residual visibilities (only if obsdata contains the residuals)

Examples

julia
julia> baselineplot(log_closure_amplitudes, Comrade.uvdist, Comrade.measwnoise, error=true)

Plot type

The plot type alias for the baselineplot function is BaselinePlot.

Attributes

alpha = 1.0 — The alpha value of the colormap or color attribute. Multiple alphas like in plot(alpha=0.2, color=(:red, 0.5)), will get multiplied.

clip_planes = @inherit clip_planes automatic — Clip planes offer a way to do clipping in 3D space. You can set a Vector of up to 8 Plane3f planes here, behind which plots will be clipped (i.e. become invisible). By default clip planes are inherited from the parent plot or scene. You can remove parent clip_planes by passing Plane3f[].

color = @inherit markercolor — The color of the marker

colorim = @inherit linecolor — The color to use for the imaginary component if fieldy returns a complex number. Note this does not cycle with the color attribute

colormap = @inherit colormap :viridis — Sets the colormap that is sampled for numeric colors. PlotUtils.cgrad(...), Makie.Reverse(any_colormap) can be used as well, or any symbol from ColorBrewer or PlotUtils. To see all available color gradients, you can call Makie.available_gradients().

colorrange = automatic — The values representing the start and end points of colormap.

colorscale = identity — The color transform function. Can be any function, but only works well together with Colorbar for identity, log, log2, log10, sqrt, logit, Makie.pseudolog10, Makie.Symlog10, Makie.AsinhScale, Makie.SinhScale, Makie.LogScale, Makie.LuptonAsinhScale, and Makie.PowerScale.

cycle = [:color] — No docs available.

depth_shift = 0.0 — Adjusts the depth value of a plot after all other transformations, i.e. in clip space, where -1 <= depth <= 1. This only applies to GLMakie and WGLMakie and can be used to adjust render order (like a tunable overdraw).

error = true — Deprecated will automatically plot errors if the field is a Measurements.Measurement

fxaa = true — Adjusts whether the plot is rendered with fxaa (fast approximate anti-aliasing, GLMakie only). Note that some plots implement a better native anti-aliasing solution (scatter, text, lines). For them fxaa = true generally lowers quality. Plots that show smoothly interpolated data (e.g. image, surface) may also degrade in quality as fxaa = true can cause blurring.

highclip = automatic — The color for any value above the colorrange.

inspectable = @inherit inspectable — Sets whether this plot should be seen by DataInspector. The default depends on the theme of the parent scene.

inspector_clear = automatic — Sets a callback function (inspector, plot) -> ... for cleaning up custom indicators in DataInspector.

inspector_hover = automatic — Sets a callback function (inspector, plot, index) -> ... which replaces the default show_data methods.

inspector_label = automatic — Sets a callback function (plot, index, position) -> string which replaces the default label generated by DataInspector.

linecap = @inherit linecap — Cap of the error bar

linecapwidth = 0 — Width of the error bar cap

linewidth = @inherit linewidth — Width of the error bar

lowclip = automatic — The color for any value below the colorrange.

marker = @inherit marker — The marker type to use

markersize = @inherit markersize — Sets the size of the markers

model = automatic — Sets a model matrix for the plot. This overrides adjustments made with translate!, rotate! and scale!.

nan_color = :transparent — The color for NaN values.

overdraw = false — Controls if the plot will draw over other plots. This specifically means ignoring depth checks in GL backends

rasterize = @inherit rasterize false — Controls whether CairoMakie rasterizes a plot (true, >0) or not (false, 0). If an integer is given, it sets the pixel density for rasterization.

space = :data — Sets the transformation space for box encompassing the plot. See Makie.spaces() for possible inputs.

ssao = false — Adjusts whether the plot is rendered with ssao (screen space ambient occlusion). Note that this only makes sense in 3D plots and is only applicable with fxaa = true.

strokecolor = @inherit markerstrokecolor — Sets the color around the markers

strokewidth = @inherit markerstrokewidth — Sets the width of the color around the markers

transformation = :automatic — Controls the inheritance or directly sets the transformations of a plot. Transformations include the transform function and model matrix as generated by translate!(...), scale!(...) and rotate!(...). They can be set directly by passing a Transformation() object or inherited from the parent plot or scene. Inheritance options include:

  • :automatic: Inherit transformations if the parent and child space is compatible

  • :inherit: Inherit transformations

  • :inherit_model: Inherit only model transformations

  • :inherit_transform_func: Inherit only the transform function

  • :nothing: Inherit neither, fully disconnecting the child's transformations from the parent

Another option is to pass arguments to the transform!() function which then get applied to the plot. For example transformation = (:xz, 1.0) which rotates the xy plane to the xz plane and translates by 1.0. For this inheritance defaults to :automatic but can also be set through e.g. (:nothing, (:xz, 1.0)).

transparency = false — Adjusts how the plot deals with transparency. In GLMakie transparency = true results in using Order Independent Transparency.

visible = true — Controls whether the plot gets rendered or not.

source
Comrade.baselineplot! Function

baselineplot! is the mutating variant of plotting function baselineplot. Check the docstring for baselineplot for further information.

source
ComradeMakieExt.BaselinePlot Type

BaselinePlot is the plot type associated with plotting function baselineplot. Check the docstring for baselineplot for further information.

source