API

VIDA

Is a image feature extraction tool for use with EHT images of black holes. It assumes that the image is close to one of the templates we have implemented and then tries to extract that feature from the image using one of the probability divergences implemented.

abstract type AbstractDivergence

Defines an AbstractDivergence type. This is the basic cost function of VIDA.

A subtype of AbstractDivergence must be a struct with at least two fields

• img::IntensityMap which holds the image you wish to run feature extraction on
• mimg::IntensityMap which is an internal buffer that stores the model image.

The divergence is evaluated roughly as

    normalize_div(div, sum(divergence_point(div, image, model)))

Therefore a user must implement the following methods

• divergence_point: Which evaluates the divergence at a single pixel
• normalize_div: Which normalizes and modifies the

abstract type AbstractMovie

Creates an abstract movie class

struct Bhattacharyya <: VIDA.AbstractDivergence

Type for the Bhattacharyya divergence. It constructed from an IntensityMap i.e. data. Additionally to evaluate the divergence we use a functor approach where if θ is your

Details

This computes the Bhattacharyya divergence which is related to Hellinger distance between two distributions. In fact, they are both minimized at the same point. The Bhattacharyya divergence is defined as

\[Bh(f_\theta||\hat{I}) = -\log\int \sqrt{f_\theta(x,y)\hat{I}(x,y)}dxdy,\]

where $\hat{I}$ is defined as the image normalized to unit flux.

struct KullbackLeibler <: VIDA.AbstractDivergence

Type for the KL divergence. It constructed from an IntensityMap i.e. data. Additionally to evaluate the divergence we use a functor approach where if θ is your

Details

This computes the KL divergence which is related to Hellinger distance between two distributions. In fact, they are both minimized at the same point. The Bhattacharyya divergence is defined as

\[KL(f_\theta||\hat{I}) = -\log\int f_{\theta}(x,y)\log \left(\frac{f_{\theta}(x,y)}{\hat{I}(x,y)}\right)dxdy,\]

where $\hat{I}$ is defined as the image normalized to unit flux.

This struct is also a functor.

struct LeastSquares <: VIDA.AbstractDivergence

Type for the least squares divergence. It constructed from an IntensityMap i.e. data. Additionally to evaluate the divergence we use a functor approach where if θ is your

Details

This computes the squared 2 norm between your image and template, both of which are normalized to unit flux.

To construct this just pass it an image object

julia> ls = LeastSquares(img::IntensityMap)

Notes

We have a template internal matrix the prevents internal allocations during computation This is a internal feature that the user does not need to worry about.

NxCorr(img::IntensityMap)

Construct the normalized cross correlation (NXCORR) divergence with respect to the image img. To maximize the NXCorr we instead compute the -log(|NxCorr|) as the divergence

NxCorr is defined as: NXCORR(n, m) = (Nσₙσₘ) Σᵢ (nᵢ - μₙ)(mᵢ - μₘ) where n and m are the two images μ is the mean and σ is the pixelwise standard deviation

Renyi(img::IntensityMap, α)

Construct the Renyi divergence with parameter α. It constructed from an IntensityMap i.e. data.

Details

This computes the KL divergence which is related to Hellinger distance between two distributions. In fact, they are both minimized at the same point. The Bhattacharyya divergence is defined as

\[Ry(f_\theta||\hat{I}) = \frac{1}{α-1}\log\int \left(\frac{f_{\theta}(x,y)^\alpha}{\hat{I}(x,y)^{\alpha-1}}\right)dxdy,\]

where $\hat{I}$ is defined as the image normalized to unit flux.

This is a very flexible divergence that reduces to many of the other divergences implemented.

• α = 1 corresponds to the KL divergence
• α = 1/2 corresponds to the Bhattacharyya divergence up to a multiplicative factor of 2

Typically we find that α=1.5 works well, as it focusses on the bright regions of the images moreso than the Bh and KL divergence. For α>2 the measure tends to devolve in something akin the to sup norm and fails to match the image structure.

struct VIDAMovie{T, F<:(IntensityMap{T, 3, D, G, A} where {D, G<:(ComradeBase.AbstractRectiGrid{D}), A<:AbstractArray{T, 3}}), I} <: VIDA.AbstractMovie

Details

Holds a X,Y,T IntensityMap plus an interpolator that lets you make a continuous movie

VIDAMovie(mov::IntensityMap{T,3})

Create an VIDAMovie class for easy interpolation between frames.

Arguments

• mov: An IntensityMap with axes (:X, :Y, :T) that represent the frames of a movie. Note that the time dimension does not have to be equi-spaced.

Returns

An VIDAMovie object that behaves like IntensityMap but lets you interpolate between frames with get_image(vida_movie, time).

struct VIDAProblem{D<:VIDA.AbstractDivergence, F, N, B}

A composite type that holds various properties that define the optimization process to extract the optimal filter.

Fields

• div: Divergence you will fit

• f: Function that defines the parametric family of templates

• autodiff: Type of autodiff to use when optimizing if any

• lb: The lower bounds of the parameter ranges to search over

• ub: The upper bounds of parameter ranges to search over

VIDAProblem(div, f, lb, ub)

Defines a VIDAProblem for optimization.

Arguments

• div: The divergence you wish to fit. This is an instantiation of VIDA.AbstractDivergence
• f: The function that defines the parametric family of templates you will fit. The function must accept a named tuple as its first argument, whose names define the parameters.
• lb: A NamedTuple whose names match the argument of f and whose values define the lower bounds of the parameters range you want to search over
• lb: A NamedTuple whose names match the argument of f and whose values define the upper bounds of the parameters range you want to search over

Example

julia> f(x) = SlashedGaussianRing(x.σ, x.s)
julia> div = Renyi(img, 0.75)
julia> lb = (x = 0.1, s = 0.001)
julia> ub = (x = 10.0, s = 0.999)
julia> prob = VIDAProblem(div, f, lb, ub)

flux(mov, t)

Returns the flux of the mov at the times time in fractional hours

blur(img, fwhm)

Blurs the img with a gaussian kernel with fwhm in μas. If fwhm is a scalar then the kernel is assumed to be symmetric, otherwise you the first entry is the fwhm in the EW direction and second the NS direction.

Returns the blurred image.

blur(mov, fwhm)

Blurs the mov with a gaussian kernel with fwhm in μas. If fwhm is a scalar then the kernel is assumed to be symmetric, otherwise you the first entry is the fwhm in the EW direction and second the NS direction.

Returns the blurred movie.

clipimage(clip, image)
clipimage(clip, image, mode)

Clips the image im according to the value clip. There are two modes for image clipping: - :relative which zeros the pixels whose intensity are below clip relative to the max. - :absolute which zeros the pixels whose intensity is below clip in Jy/pixel

divergence(d, m)

Computes the divergence of d with respect to the template m. This measure how similar the divergence with respect to the image stored in d is compared to the template.

Arguments

• d: An divergence or a subtype of VIDA.AbstractDivergence
• m: A template which ComradeBase.AbstractModel or more commonly a VIDA.AbstractImageTemplate

divergence_point(div::AbstractDivergence, p, q)

Evaluate the divergence div at a single pixel with value p and q for image and template. The total divergence

get_frames(mov)

Gets all the frames of the movie object mov. This returns a array of IntensityMap objects.

get_image(mov, t; keeptime)

Gets the frame of the movie object mov at the time t. This returns an IntensityMap object at the requested time. The returned object is found by linear interpolation.

get_times(mov)

Returns the times that the movie object mov was created at. This does not have to be uniform in time.

join_frames(times, images)

Joins an array of IntensityMap at specified times to form an VIDAMovie object.

Inputs

• times: An array of times that the image was created at
• images: An array of IntensityMap objects

Outputs

VIDAMovie object

load_hdf5(filename; polarization, style)

Loads an hdf5 movie file where filename should be a HDF5 file.

Details

This reads in a hdf5 movie file and outputs and VIDAMovie object. If polarization=true we will read in all stokes parameters.

Notes

Currently this only works with movies created by ehtim.

load_im_h5(fname; polarization)

where fname should be a hdf5 image file generated using illinois hdf5 process

Details

The function returns an IntensityMap object that contains the relevant image and parameters extracted from the fits file. It also ensures that we are astronomers and that the image using sky-left coordinates.

load_image(fname; polarization)

where fname is a image file. This will call a method based on the file extension. Curerntly we have .fits, and .h5 implemented.

Details

This reads in a fits file that is more robust to the various imaging algorithms in the EHT, i.e. is works with clean, smili, eht-imaging.

The function returns an IntensityMap object that contains the relevant image and parameters extracted from the fits file. It also ensures that we are astronomers and that the image using sky-left coordinates.

normalize_div(div, val)

Returns the proper divergence whose values have been normalized so that the value is always postive semi-definite and is only zero then the template and the image are equal.

nxcorr(img1, img2)

Returns the normalized cross correlation (NXCORR) between img1 and img2. NXCORR is defined as

NXCORR(n, m) = (Nσₙσₘ) Σᵢ (nᵢ - μₙ)(mᵢ - μₘ)

where n and m are the two images μ is the mean and σ is the pixelwise standard deviation

save_hdf5(filename, mov; style)

Saves and hdf5 file where filename is the write out location. Currently style only works with ehtim, namely we save HDF5 files that only work with ehtim.

triptic(img::SpatialIntensityMap, template)

Plots a triptic where the left panel is the img middle the template and the right a two cross-sections of the image and template

Runs the VIDA algorithm to find the optimal template defined by the problem prob. The optimization will be run using the Optimization.jl optimizer. You can optionally pass a random number generator used to specify the initial points or can explicitly pass the initial location of the optimization using init_params. By default vida first transforms the parameter space to the unit hypercube. If you do not wish to do this set unit_cube = false.

The remaining kwargs... are forwarded to the Optimization.solve function.

Arguments

• prob: Defines the problem you wish to solve.
• optimizer: Specifies the optimizer you want to use. Any optimizer from Optimization.jl works.

Optional Keyword arguments

• rng: Specifies the random number generator you want to use to select the initial points. Note that is not forwarded to the optimizers since not all can use a specific rng.
• init_params: Specify the initial point of the optimization routine. The default nothing will randomly draw a starting point from the parameter space
• unit_cube: If true the parameters are first transformed to the unit hypercube. If false they are transformed to ℝⁿ using TransformVariables
• kwargs...: Additional options to be passed to the solve function from Optimization.jl

regrid(mov::VIDAMovie, npix, xlim, ylim)

Inputs

• mov::VIDAMovie : Movie you want to regrid
• npix : Number of pixels in x and y direction
• xlim : Tuple with the limits of the image in the RA
• ylim : Tuple with the limits of the image in DEC