GitHub

Introspection

When one loads a model through, say, a command like

julia> model = gravitymodel("model.yml")

then, Gravity.jl creates and load a new module (named, by default GravityModel), that defines a number of functions that are needed for the analysis of the lensing system.

Note that the use of a module separates model-related constants and functions, and makes it possible to load at the same time more than a model. To obtain a summary of the model loaded, one can just use the help and type

help?> model

[use the question mark ? to enter the help mode in the REPL].

Most relevant generated functions and constants are described in this section. One can display the generated source-code of the module using the debug keyword in gravitymodel: that is,

julia> model = gravitymodel("model.yml"; debug=:all)

Since this command is generating a long output, one might want to load the package TerminalPager.jl and prepend the call to this function with the macro @stdout_to_pager; more simply, one can just enter | at the beginning of the REPL to activate the pager.

Alternatively, one can display the definition of a specific function by using its name (as a symbol) in place of :all. Hence, for example,

julia> model = gravitymodel("model.yml"; debug=:loglikelihood)

Basic functions and constants

The functions listed below are generally strictly needed for the Bayesian inference of the model parameters.

Main.model.PriorDistribution — Type
PriorDistribution()

A multivariate distribution representing the prior on the lensing model.

This distribution can be used to compute the log-prior, to generate random samples, and to compute the quantiles. The distribution has no parameters and thus can be instatiated just with PriorDistribution(). As a convenience, the constant prior is a singleton instance of PriorDistribution.

source
Main.model.ReferenceDistribution — Type
ReferenceDistribution()

A multivariate distribution representing the reference on the lensing model.

This distribution can be used to generate random samples for the initial proposal in the MCMC algorithm. The distribution has no parameters and thus can be instatiated just with ReferenceDistribution().

As a convenience, the constant reference is a singleton instance of ReferenceDistribution.

Note that, contrary to the prior distribution, the reference distribution can be chosen to reflect the expected distribution of the parameters in the model after looking at the data.

source
Main.model.likelihood_options — Function
likelihood_options()

Perform preliminary computations for all likelihood-related functions.

This function returns a LikelihoodOptions structure, which has to be passed to all likelihood-related functions such as loglikelihood or logposterior. Essentially, it performs calculations that can be factorized out the likelihood, i.e. that do not depend on the particular choice of the free parameters but only on the fixed constants.

For convenience, the constant options is defined as the result of a call to this function.

source
Main.model.loglikelihood — Function
loglikelihood(x [, options])

Compute the loglikelihood as a function of the x parameters.

In the standard (scalar) use, x is a vector with the free parameters, as specified in the setup file (the parameters' names and dimensionality can be checked using setup.parameters).

If x is a n × d matrix, the result is vector of length n, with the loglikelihood values computed in for each of the n sets of the d parameters.

The options can be pre-computed using likelihood_options.

source
Main.model.logposterior — Function
logposterior(x [, options])

Compute the logposterior as a function of the x parameters.

The logposterior is effectively obtained as the sum of the log of prior and of the log of the likelihood: so

@assert logposterior(x) ≈ logpdf(PriorDistribution(), x) + loglikelihood(x)

In the standard (scalar) use, x is a vector with the free parameters, as specified in the setup file (the parameters' names and dimensionality can be checked using setup.parameters).

If x is a n × d matrix, the result is vector of length n, with the logposterior values computed in for each of the n sets of the d parameters.

The options can be pre-computed using likelihood_options.

source
Main.model.make_cosmo — Function
make_cosmo(x [, options])

Return the cosmological model associated with the parameter vector x.

The returned model is an instance of Cosmology. The options can be pre-computed using likelihood_options. If the used configuration does not include free parameters in the cosmology section, the cosmological model returned is directly as options.cosmo (and, clearly, will be independent of x).

source
Main.model.make_grid — Function
make_grid(x [, options])

Return the grid created for a given parameter vector x.

The grid, if present, is used for the lens inversion. Therefore, it is created only if a grid: specification is present in the configuration file and if scheme: [simplfied-]imageplane. In all other cases this function returns nothing.

source
Main.model.make_images — Function
make_images(x [, options, sys])

Returns the images associated to the parameters x.

The function retunrs a vector of SourcePlane built on ImageFamily'ies. It just reports all image families specified in the configuration, organized as the code has parsed them.

Generally, the images do not depend on the parameters x, but the code still requires the use of the x argument to build the entire ImageFamily structure.

The options can be pre-computed using likelihood_options, while sys can be pre-computed using make_lenssys.

See also: make_sources, make_predictions, make_lenssys

source
Main.model.make_predictions — Function
make_predictions(x [, options, sys])

Compute predicted images associated to the parameters x.

The predictions is a vector of SourcePlane built on ImageFamily'ies, but it lists the various families predicted from the lensing model associated to the parameters x. The inversion of the lens equation for the various sources is carried out using the same technique employed by the likelihood: therefore, it will be grid-based if a sampling.grid option is present and sampling.scheme: imageplane.

The options can be pre-computed using likelihood_options, while sys can be pre-computed using make_lenssys.

Warning

This function ignores all passed keywords: in particuar, it does not take into account the matches keyword, and just returns all predicted images found, with no duplicates. Please use flat_predictions with matches=:best to select the best-matching prediction associated with each image.

See also: make_sources, make_images, make_lenssys

source
Main.model.make_sources — Function
make_sources(x [, options, sys]; kw...)

Compute sources associated to the parameters x.

The function returns a vector of SourcePlane's built on SingleSource's, containing the data of all sources. If scheme: imageplane, generally some of the parameters x will be used to describe the sources (for example, their positions on the lens plane). If, instead, scheme: sourceplane, the sources are created by lensing the images, and taking appropriate weigthed averages of all lensed images belonging to the same ImageFamily.

The options can be pre-computed using likelihood_options, while sys can be pre-computed using make_lenssys.

The optional keywords parameters can be used to modify the computation of the source.

See also: make_images, make_predictions, make_lenssys

source
Main.model.make_fiducial_sources — Function
make_fiducial_sources(x [, options, sys])

Compute fiducial sources associated to the observed imagess.

The function returns a vector of SourcePlane's built on SingleSource's, containing the data of all sources. This functionn is identical to make_sources for scheme: sourceplane, i.e. it returns, for each ImageFamily, an appropriate weigthed averages of all lensed images.

The options can be pre-computed using likelihood_options, while sys can be pre-computed using make_lenssys.

See also: make_sources, make_images, make_predictions, make_lenssys

source

Helper functions

Structures returned by commands such as Main.model.make_lenssys and Main.model.make_sources are hierarchical. This is mainly done for optimization reasons, but this might make it complicated to access the individual lenses, sources, or images. Therefore, the Gravity.jl provides some helper functions which "flatten" the hierarchical structure into simple dictionaries. This is generally obtained using the Gravity.flatten function.

Main.model.flat_lenses — Function
flat_lenses(x [, prop1, prop2...])

Return a (flat) dictionary with all predicted lenses in the model.

x is a vector with the free parameters, as specified in the setup file.

The optional arguments can be used to extract specific properties from the lenses: they must be symbols corresponding to the properties desired.

See also: flat_sources, flat_images, flat_predictions.

source
Main.model.flat_sources — Function
flat_sources(x [, options, sys])

Return a (flat) dictionary with all sources in the model.

x is a vector with the free parameters, as specified in the setup file.

The optional arguments can be used to extract specific properties from the sources: they must be symbols corresponding to the properties desired.

See also: flat_images, flat_predictions, flat_diffs, flat_lenses.

source
Main.model.flat_images — Function
flat_images(x [, prop1, prop2...])

Return a (flat) dictionary with all images in the model.

x is a vector with the free parameters, as specified in the setup file.

The optional arguments can be used to extract specific properties from the images: they must be symbols corresponding to the properties desired.

See also: flat_sources, flat_predictions, flat_diffs, flat_lenses.

source
Main.model.flat_predictions — Function
flat_predictions(x [, prop1, prop2...]; matches=:best, tolerant=false)

Return a (flat) dictionary with all predicted sources in the model.

x is a vector with the free parameters, as specified in the setup file.

The optional arguments can be used to extract specific properties from the predicted images: they must be symbols corresponding to the properties desired.

The two keywrods matches and tolerant modify the behaviour of this function. If matches=:all, the function returns all predicted images, regardless of the observed ones. In this case the names of the predicted images are modified to end in prediction_n, where n is a running integer. The default setting, matches=:best, associates to each image the closest prediction. In this case, if tolerant=true the function silently ignores images for which the lens inversion procedure failed.

See also: flat_sources, flat_images, flat_diffs, flat_lenses.

source
Main.model.flat_dataset_predictions — Function
flat_dataset_predictions(x)

Return a (flat) dictionary with all predicted datasets in the model.

x is a vector with the free parameters, as specified in the setup file.

The function returns a dictionary with keys equal to the various dataset names, and values with the predicted dataset, i.e. objects of type Gravity.Dataset. Each predicted dataset is identical to the input dataset (which can be obtained from the dictionary model.data), with the difference that its fluxes correspond to the predicted fluxes.

Tip

Datasets are rather complex objects, but one can easily extract their main properties with the functions Gravity.flux, Gravity.mask, and Gravity.ivar. Note also that one can save a dataset to a FITS file with a simple call to Gravity.fits_write.

source
Main.model.flat_unlenses — Function
flat_unlenses(x [, prop1, prop2...])

Return a (flat) dictionary with all predicted lenses in the model.

x is a vector with the free parameters, as specified in the setup file.

The optional arguments can be used to extract specific properties from the lenses: they must be symbols corresponding to the properties desired.

This function is identical to flat_lenses, but uses make_unlenssys to compute the lens system. This way, all input lens parameters are taken as if lens_lenses=false has been set.

See also: flat_sources, flat_images, flat_predictions.

source
Main.model.flat_mahalanobis — Function
flat_mahalanobis(x, prop1 [, prop2...])

Return a (flat) dictionary with Mahalanobis distances between observed and predicted images.

x is a vector with the free parameters, as specified in the setup file.

The function returs a dictionary with the Mahalanobis distances between the observed and the predicted images with respected to the requested properties. The Mahalanobis distance is computed as

\[ d = \sqrt{(\mathbf{x} - \mathbf{y})^T \mathbf{C}^{-1} (\mathbf{x} - \mathbf{y})} \; ,\]

where $C$ is the (co)variance of the image property, and $\mathbf{x}$ and $\mathbf{y}$ are the predicted and observed values.

See also: flat_sources, flat_images, flat_predictions, flat_lenses.

source
Main.model.rms — Function
rms(x, [prop=:x])

Compute the root-mean-square of the model predictions.

x is a vector with the free parameters, as specified in the setup file. The function computes the root-mean-square of the differences between the observed and the predicted images with respected to the requested property. By default, the property is :x, i.e. the positions of the images.

source
Main.model.flat_chisq — Function
flat_chisq(x)

Return a dictionary with the chi-square contribution of each family.

This function returns a dictionary (indexed with the source names) with all chi-square contributions from the associated images, saved in a named tuples. The contributions are separated by type:

For the first four contribution one can obtain a similar result (divided by each image, instead of collected for each image-family) with using Main.model.flat_mahalanobis: for example, the contribution from image positions can be computed as

julia> sum(abs2, values(model.flat_mahalanobis(x, :x)))

and similarly for the other terms (mag, Q, and t). The chi-square terms for Δt and ΔmagΔt, instead, cannot be computed this way, since they are associated with a covariance matrix: as a result, the corresponding chi-square term is not just a sum of (weighted) squares. The pixel term contains contributions from pixel scatters of individual extended images. Finally, the image part is collective computed for all extended sources: this is so because extended sources can have overlapping pixels, which prevent individual computation of chi-square contributions.

Warning

Since extended images might overlap, the pixel term should be taken with a lot of care. In particular, for images with overlapping masks this term will return wrong results. The image collective term is always correct, but cannot discriminate among different families.

source
Main.model.chisq — Function
chisq(x)

Return a tuple with the chi-square contributions from various terms.

The names returned depend on the source types, and can be

Warning

For non-overlapping sources, the pixel and the image terms should be identical. More generally, ignore the pixel term and use just the image one, which correctly takes into account the effects of overlapping extended sources.

source
Main.model.counterimages — Function
counterimages(pars; kw...)

Return a dictionary of all counterimages for each source in the model.

This function compute counterimages associated to each source, and returns a dictionary whose values are vectors of tuples. See flat_counterimages for details.

The returned object has keys associated to each source, with entries such as "sources2Point".

Keywords

  • grid: either :auto (default) or a Gravity.MaskedGrid. If :auto, the grid is automatically created using an educated guess based on the lens system properties.
  • criticalcurves: integer (default: 3). If grid=:auto, the maximum refinement level around critical curves.
  • images: integer (default: 3). If grid=:auto, the maximum refinement level around observed image positions.
  • lenses: integer (default: 3). If grid=:auto, the maximum refinement level around lens positions.
  • tolerance: float (default: -1.0). If positive, counter-images that are within this tolerance (in arcseconds) from predicted images are not included in the output.
Tip

Counterimages in general differ from predicted images, as they are always computed using a grid-based lens inversion. By using a positive value for the tolerance keyword, one check that a model does not predict un-observed images.

source
Main.model.flat_counterimages — Function
flat_counterimages(pars; tolerance=0,kw...)

Return a dictionary of all counterimages for each source in the model.

Counterimages are computed using the function counterimages and are returned as named tuples with the following fields:

  • x: the coordinates of the counterimage, as a 2-tuple;
  • Δmag: the image magnification, expressed as a difference in magnitudes with respect to the source (if positive, it means that the image is magnified);
  • t: the time-delay with respect to the first image, in days.

The optional keyword tolerance (default: 0) can be used to specify a tolerance (in arcseconds). If set to a positive value, counterimages that are within this tolerance from predicted images are not included in the output. This can be useful to check that a model does not predict un-observed images.

The returned object is a ordered dictionary with entries such as "sources2Pointcounterimages1".

!!! tip Counterimages are quite similar to predictions, with the important difference that a grid is always used to compute their positions. By default the grid is automatically chosen using an educated guess based on the lens system properties. Keywords control the grid definition and refinement and are directly passwd to counterimages and Gravity.counterimages.

source
Main.model.count_counterimages — Function
count_counterimages(pars; abstol=0.1, reltol=3.0, kw...)

Count mismatches in the observed and predicted images.

This function uses counterimages to compute all counterimages associated to the sources in the model, and compares them to the observed images. For each source, it returns a tuple with two integers: the first is the number of observed images that cannot be matched to any predicted image (within the specified tolerances), while the second is the number of predicted images that cannot be matched to any observed image.

The two optional keywords abstol and reltol specify the tolerances for the matching procedure: abstol is in units of the model (as specified for example by the angular_unit parameter in the configuration), while reltol is a dimensionless factor multiplying the covariance matrix of each observed image.

All other keywords are passed to counterimages to define the grid used to compute the counterimages.

source
Main.model.get_lensnames — Function
get_lensnames([x])

Return a flat string vector with all lens names.

In principle, in case of non-stable lenses, this is a function of the parameters x. However, if the lens redshift are fixed, the lens names are always the same, and this function can be called without arguments.

source
Main.model.get_sourcenames — Function
get_sourcenames([x])

Return a flat string vector with all source names.

In principle, in case of non-stable sources, this is a function of the parameters x. However, if the source redshift are fixed, the source names are always the same, and this function can be called without arguments.

source
Main.model.get_imagenames — Function
get_imagenames([x])

Return a flat string vector with all image names.

In principle, in case of non-stable sources, this is a function of the parameters x. However, if the source redshift are fixed, the image names are always the same, and this function can be called without arguments.

source
Main.model.checkparameters — Function
checkparameters(pars)

Check if the given parameters are within the prior support.

Returns an (ordered) dictionary with the names of the parameters that are outside the prior support as keys, and their values as values. If all parameters are within the prior support, an empty dictionary is returned.

source

Conversion functions and constants

The conversion between celestial and pixel coordinates is carried out by taking into account the WCS parameters specified in the configuration file. To help this process, the model has a number of associated constants and functions.

Main.model.wcsmap — Function
wcsmap(α, δ)

Map the celestial coordinates (α, δ) into pixel coordinates (x, y).

The function returns a tuple (x, y). It uses global_wcs to perform the transformation.

See also: wcsunmap.

source
Main.model.wcsunmap — Function
wcsunmap(x, y)

Map the celestial coordinates (x, y) into celestial coordinates (α, δ).

The function returns a tuple (x, y). It uses global_wcs to perform the transformation.

See also: wcsmap.

source

In particular, the function wcsmap can be used directly in the configuration file to convert celestial coordinates $(\alpha, \delta)$ into the global coordinate system used by Gravity.jl. A typical usage would be to set, for a lens or for an image:

x: wcsmap(ra, dec)

Other conversion functions are associated to the use of quadrupole moments.

Gravity.rqθ_to_Q — Function
rqθ_to_Q(r, q, θ)

Convert the ellipticity parameters into a quadrupole moment.

The ellipticity parameters are defined as follows:

  • q: the axis ratio of the ellipse, i.e. the ratio of the minor to the major axis;
  • θ: the position angle of the major axis, in radians;
  • r: the geometric mean of the semiaxes, i.e. the square root of their product.

The function returns a quadrupole moment Q as an SVector{3}.

source
Gravity.Q_to_rqθ — Function
Q_to_rqθ(Q)

Convert the quadrupole moment into ellipticity parameters.

The quadrupole moment Q is an SPointMap or an SVector{3} and the function returns the ellipticity parameters (r, q, θ) as a tuple.

The ellipticity parameters are defined as follows:

  • r: the geometric mean of the semiaxes, i.e. the square root of their product.
  • q: the axis ratio of the ellipse, i.e. the ratio of the minor to the major axis;
  • θ: the position angle of the major axis, in radians;
source

Note that these functions, instead, should not be used inside the configuration file, as this is unnecessary and can produce inefficient code. Rather, when specifying the measured ellipticity of an elliptical source, one can just use the notation

Q: Wishart(a, b, θ, ν)

where a and b are the measured semi-axes, θ is the position angle, and ν is estimated as $\nu \approx 4 a^2/\sigma_a^2 ≈ 4 b^2/\sigma_b^2$, i.e. is associated to the square inverse of the relative errors on the measurements of a and b.

Similarly, when defining the shape of a source, use

Q: IWishart(a, b, θ, ν)

Note that this will leave Q as free parameter. A posteriori, after the MCMC has been completed through runmodel, one can convert the sampled quadrupole moments Q into more common parameters using Gravity.Q_to_rqθ.