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=splat(*),] 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(splat(*), 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


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

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.

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=Optimization.NoAD(), 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.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