Comrade

Composable Modeling of Radio Emission

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.

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.

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

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

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.

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.

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.

AbstractVisibilityDatum

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

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

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.

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.

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

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.

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.

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.

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

A single datum of an ArrayConfiguration

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

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

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

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

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

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

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

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)

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

J = GenericJones() do
    return x.j11, x.j21, x.j12, x.j22
end

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

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

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

IdealInstrument(array::AbstractArrayConfiguration)

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

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

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.

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

D = JonesD() do
    d1 = complex(x.d1real, x.d1imag)
    d2 = complex(x.d2real, x.d2imag)
    return d1, d2
end

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.

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

G = JonesG() do
    g1 = exp(x.lg1 + 1im.*x.gp1)
    g2 = g1*exp(x.lgratio + 1im.*x.gpratio)
    return g1, g2
end

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.

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

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

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.

MemoryStore

Stored the HMC samplers in memory or ram.

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.

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.

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.

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

struct ScanSeg <: Comrade.Segmentation

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

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

SingleReference(site::Symbol, val)

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

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

G = SingleStokesGain(x->exp(x.lg + xp.gp))

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.

SiteArray(arr, sitelookup::SiteLookup)

Construct a site array with the entries arr and the site ordering implied by sitelookup.

SiteLookup(s::SiteArray)

Construct a site lookup dictionary for a site array.

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.

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> st = timetable(obs)
julia> st[begin] # grab first scan
julia> st[end]   # grab last scan

struct TrackSeg <: Comrade.Segmentation

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

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.

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

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)

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.

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.

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.

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.

baseline(scan::Scan)

Return the baselines for each datum in a scan

beamsize(ac::AbstractArrayConfiguration)

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

beamsize(obs::AbstractObservationTable)

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

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.

caltable(s::SiteArray)

Creates a calibration table from a site array

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.

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.

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.

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.

datatable(obs::AbstractVLBITable)

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

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.

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.

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.

domain(ac; executor, header)

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

domain(obs::AbstractObservationTable; executor=Serial(), header=ComradeBase.NoHeader()

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

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.

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.

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.

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.

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

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.

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.

idealvisibilities(m::AbstractSkyModel, x)

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

instrumentmodel(d::AbstractVLBIPosterior)

Returns the instrument model of the posterior d.

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.

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.

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.

loglikelihood(d::AbstractVLBIPosterior, θ)

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

logprior(d::AbstractVLBIPosterior, θ)

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

postsamples(s::PosteriorSamples)

Return the samples from the PosteriorSamples object s

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.

quadrangle(b::EHTClosurePhaseDatum)

Returns the sites used in the closure amplitude quadrangle.

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.

residual(post::AbstractVLBIPosterior, p)

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

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.

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.

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)

samplerinfo(s::PosteriorSamples)

Return the metadata from the PosteriorSamples object s.

samplerstats(s::PosteriorSamples)

Return the sampler statistics from the PosteriorSamples object s

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.

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

sites(d::AbstractArrayConfiguration)

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

sites(d::AbstractObservationTable)

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

sites(g::CalTable)

Return the sites in the calibration table

skymodel(post::AbstractVLBIPosterior, θ)

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

skymodel(d::AbstractVLBIPosterior)

Returns the sky model of the posterior d.

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.

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

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]

triangle(b::EHTClosurePhaseDatum)

Returns the sites used in the closure phase triangle.

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

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

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

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