Comrade.jl

Appearance

Comrade API ​

Index ​

Comrade.Comrade Module
julia
Comrade

Composable Modeling of Radio Emission

source

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.EHTArrayBaselineDatum, F, T, S, D<:(AbstractArray{A<:Comrade.EHTArrayBaselineDatum})} <: Comrade.AbstractArrayConfiguration{A<:Comrade.EHTArrayBaselineDatum}

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

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.AbstractSkyModel Type
julia
AbstractSkyModel

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

Any subtype must implement the following methods

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

The following methods have default implementations:

  • idealvisibilities(m::AbstractSkyModel, x): Computes the ideal 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.

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.idealvisibilities Function
julia
idealvisibilities(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.CalTable Type
julia
struct CalTable{T, 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 stamp. 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.sites Method
julia
sites(g::CalTable)

Return the sites in the calibration table

source

Comrade.IIDSitePrior Type
julia
IIDSitePrior(seg::Segmentation, dist)

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.

Example

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

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

source

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

Construct a prior for an entire array of sites.

  • The default_dist is the default distribution for all sites. Currently only IIDSitePrior is supported.

  • 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.

Example

julia
p = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0, 0.1)); LM = IIDSitePrior(ScanSeg(), Normal(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.Segmentation

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.Segmentation

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

source

Comrade.TrackSeg Type
julia
struct TrackSeg <: Comrade.Segmentation

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

source

Comrade.timestamps Function
julia
timestamps(seg::Segmentation, 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.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.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(x.lg1 + 1im.*x.gp1)
    g2 = g1*exp(x.lgratio + 1im.*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.

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())

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.

Example

A Stokes I example is

julia
julia> G = SingleStokesGain(x->exp(x.lg + 1im*x.pg))
julia> intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(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.(x.lgr + 1im*x.gpr)
        gL = gr*exp.(x.lgrat + 1im*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(), Normal(0.0, 0.1)),
                    gpr = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))),
                    lgrat = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(0.0, 0.1)),
                    gprat = ArrayPrior(IIDSitePrior(ScanSeg(), DiagonalVonMises(0.0, inv(π^2))),
                    dRre = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.1)),
                    dRim = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.1)),
                    dLre = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(0.0, 0.1)),
                    dLim = ArrayPrior(IIDSitePrior(TrackSeg(), Normal(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

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

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.forward_model 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.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=nothing)

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. We recommend using Enzyme.set_runtime_activity(Enzyme.Reverse) 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 = Uniform(μas2rad(10.0), μas2rad(30.0)), a = Uniform(1.0, 10.0))
g  = imagepixels(μas2rad(100.0), μas2rad(100.0), 256, 256)
skym = SkyModel(sky, skyprior, g)

G = SingleStokesGain(x->exp(x.lg + 1im*x.pg))
intprior = (lg = ArrayPrior(IIDSitePrior(ScanSeg(), Normal(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

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 y from parameter space to the transformed space (e.g. unit hypercube if using ascube).

For the inverse transform see transform

source

HypercubeTransform.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

HypercubeTransform.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 = ascube(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

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)

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.

source

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

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.

Keyword Arguments

  • table: A string specifying the table you wish to read in. There are two options: :samples which corresponds to the actual MCMC chain, and stats which corresponds to additional information about the sampler, 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.

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, 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)

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.

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, Normal(0.0, 0.1); reference=:AA, LM = Normal(0.0, 1.0))
(AP = Normal(0.0, 0.1), LM = Normal(0.0, 1.0), PV = Normal(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.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

eht-imaging interface (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