Post processing utilities

mean(chains[, params; kwargs...])

Calculate the mean of a chain.

std(chains[, params; kwargs...])

Calculate the standard deviation of a chain.

median(chains[, params; kwargs...])

Calculate the median of a chain.

quantile(chains[; q = [0.025, 0.25, 0.5, 0.75, 0.975], append_chains = true, kwargs...])

Compute the quantiles for each parameter in the chain.

Setting append_chains=false will return a vector of dataframes containing the quantiles for each chain.

savechain("chain.chn", chain [, format]; show=true, bare=false, type=nothing)

Save a chain to a file.

This function saves the model, the chain, and the source code of the model to a file. The format can be either :serialization (the default) or :jld2 (if the filename ends with .jld2).

The :serialization format uses the Serialization.jl standard library to serialize the chain and the model. This is the default format and it does not require any additional package. Note however that the Serialization format is not guaranteed to be compatible across different versions of Julia.

The :jld2 format requires the JLD2.jl package to be installed and imported. Files saved in this format are guaranteed to be compatible across different versions of Julia. Note however that saving and loading chains in this format can be significantly slower than using the :serialization format.

This function saves the model associated to the chain, so that when the chain is loaded again with loadchain the model is also restored.

The saved file can be loaded again with loadchain.

Keywords

• show::Bool=true: show a message when the chain is saved.
• bare::Bool=false: save only the chain, without the model and the source code. This is useful when the model is available somewhere else (for example, as an external YAML configuration file).
• type::Type: save the chain with the specified type. The use of type=Float32 or even type=Float16 can reduce the disk space, but it can also lead to loss of precision. The default is nothing, which saves the chain with the same type as the original chain.

See also: loadcheckpoint, loadchain.

Extended help

This function deals with the complication that the chain is generally associated to a gravity model, i.e. a Julia module. Specifically, if bare=false, it saves for each chain the following:

• All possible subchains, i.e. all chains defined in the mcmc_init field of the model (this is done recursively)
• The name of the model, which is the name of the module
• The raw source of the model, i.e. the original YAML content of the file passed to gravitymodel
• The filename passed to gravitymodel
• The ctime (creation time) and mtime (modification time) of the gravity model source file
• All keywords passed to gravitymodel
• The chain itself.

Essentially, all these data are ensure that we are able to rebuild the module associated with the chain when we load the chain again.

chain = loadchain("chain.chn" [, format]; show=true, bare=false, type=nothing)

Load a chain from a file.

This function loads a chain that was saved with savechain.

The optional format argument can be either :serialization or :jld2. If format is not specified, the format is inferred from the file magic number (the initial bytes of the file). If the file was saved with the :jld2 format, the JLD2.jl package must be installed and imported.

If originally present in the file and if bare=false (the default), the model associated to the chain is saved back in chain.info.model.

Keywords

• show::Bool=true: show a message when the chain is loaded.
• bare::Bool=false: loads only the chain, without the associated model and any other additional information. This is useful when the model is available somewhere else (for example, as an external YAML configuration file) or to speed up significantly the loading process (since the model is not loaded and therefore also not compiled).
• type::Type: converts the chain to the specified type. The default is nothing, which loads the chain with the same type as the original chain.

See also: loadcheckpoint, savechain.

chain = loadcheckpoint(model [, path])

Load a chain from a checkpoint and create a MCMCChains.Chains object.

To use this function one must have set checkpoint: true in the configuration file associated to the model. This allows one to load intermediate chains and check the progress of the MCMC run.

If path is not specified, the latest checkpoint saved in the results/ directory is loaded ("results/latest").

lp, x = bestlp(chains)

Return the best log-posterior and the corresponding parameter vector.

print_values([io::IO,] model, x; auto_rename=false)
print_values([io::IO,] gravitysetup, x; auto_rename=false)

Print parameter names and values for a GravitySetup and parameter vector x.

This function is generally used to associate the elements of a parameter vector x to the corresponding parameter names. It can be used in conjuction with bestlp to print the best-fit parameter values:

julia> print_values(model, Gravity.bestlp(chain)[2])

If the optional keyword auto_rename is set to true, the parameter names are simplified similarly to auto_rename_chain before printing.

flatten(x [, :prop1, :prop2, ...])

Return a flattened vector for hierarchical structures.

This function goes recursively through all elements of a hierarchical structure x and returns a flattened vector with all final elements. Optionally, one can provide a single property (as a symbol, like :x) or a list of properties to be extracted from each element. If a property is not found, the value missing is returned.

This function can be applied to a MCMCChains.Chains object, to a LensSystem, to a SourcePlane, and to an ImageFamily; it can also be applied to tuples or vectors of such objects.

chain_sources(chain[, :prop1, :prop2, ...])

Return a flattened vector of sources for each chain element.

This function goes through all vectors of parameters in the chain and returns a flattened vector with all sources. Optionally, one can provide a single property or a list of properties to be extracted from each source.

See also: chain_lenses, chain_predictions, chain_diffs, flatten.

chain_predictions(chain[, :prop1, :prop2, ...])

Return a flattened vector of predicted images for each chain element.

This function goes through all vectors of parameters in the chain and returns a flattened vector with all predicted images. Optionally, one can provide a single property or a list of properties to be extracted from each lens.

See also: chain_lenses, chain_sources, chain_diffs, flatten.

chain_diffs(chain, :prop1, [:prop2, ...])

Return a flattened vector of predicted images for each chain element.

This function goes through all vectors of parameters in the chain and returns a flattened vector with the difference between properties of observed and predicted images.

See also: chain_lenses, chain_sources, chain_predictions, flatten.

chain_lenses(chain[, :prop1, :prop2, ...])

Return a flattened vector of lenses for each chain element.

This function goes through all vectors of parameters in the chain and returns a flattened vector with all lenses. Optionally, one can provide a single property or a list of properties to be extracted from each lens.

See also: chain_sources, chain_predictions, flatten.

auto_rename_chain(chain)

Return chain with simplified names.

This function tries to replace the names of chain parameters returned by runmodel with significantly shorter ones. Specifically, it goes through all cosmology, source, and image parameter names, and tries to keep the most significant parts of a column name.

For example, a parameter name as lenses_dm1_LensTool-NIE_σ might be shortened to just σ. However, if there is also another parameter called (say) lenses_5_LensTool-NIS_σ, this function might choose to call the two parameters dm1_σ and lens-5_σ.

It is often recommenable to use this function before creating a pairplot from chain, to avoid creating confused results.

massdensity(lenssys, θ)
massdensity(lenssys)

Compute the dimensional mass density at the position θ.

In second version, return a closure that, when called with a point θ as an argument, will compute the mass density at that location.

This method sums over all lens planes in lenssys and computes the final mass density in units of M⊙ / pc². The computation is performed along the line of sight θ, ignoring any lensing effect.

See also: lensed_massdensity, mass_sqarcsec

lensed_massdensity(lenssys, θ)
lensed_massdensity(lenssys)

Compute the dimensional mass density at the position θ.

In second version, return a closure that, when called with a point θ as an argument, will compute the mass density at that location.

This method sums over all lens planes in lenssys and computes the final mass density in units of M⊙ / pc². The computation is performed along the line of sight θ, taking into account the lensing effects.

See also: massdensity, lensed_mass_sqarcsec.

mass_sqarcsec(sys, θ)
mass_sqarcsec(sys)

Compute the angular mass density at the position θ.

In second version, return a closure that, when called with a point θ as an argument, will compute the mass density at that location.

This method sums over all lens planes in lenssys and computes the final mass density in units of M⊙ / arcsec². The computation is performed along the line of sight θ, ignoring any lensing effect.

See also: lensed_mass_sqarcsec, massdensity.

lensed_mass_sqarcsec(lenssys, θ)

Compute the angular mass density at the position θ.

In second version, return a closure that, when called with a point θ as an argument, will compute the mass density at that location.

This method sums over all lens planes in lenssys and computes the final mass density in units of M⊙ / arcsec². The computation is performed along the line of sight θ, taking innto account the lensing effect.

See also: mass_sqarcsec, lensed_massdensity.

radialprofile(lenssys, x₀, r [; kw...])

Compute the average density radial profile at radius r.

This function returns the radial profile of the lens system lenssys at the radius r from the center x₀. The result is the average density in units of M⊙ / arcssec². Note that, as with any other angular parameter, r must be expressed in units given by Gravity.angular_unit.

The keywords (such as rtol or atol) are passed directly to the integration routine, QuadGK.quadgk.

In case one desires to have the radial profile at different radii, one can use the standard broadcasting syntax, e.g. radialprofile.(sys, Ref(x₀), [1, 2, 3]) (note the use of Ref to avoid broadcasting over the components of x₀).

See also: cumradialprofile, mass_sqarcsec.

cumradialprofile(lenssys, x₀, r [; kw...])

Compute the cumulative radial profile up to radius r.

This function returns the mass of the lens system lenssys within a disk of radius r centered at x₀. The result is the mass in units of M⊙. Note that, as with any other angular parameter, r must be expressed in units given by Gravity.angular_unit.

In case one desires to have the radial profile at different radii, one can use the standard broadcasting syntax, e.g. cumradialprofile.(sys, Ref(x₀), [1, 2, 3]) (note the use of Ref to avoid broadcasting over the components of x₀).

The keywords (such as rtol or atol) are passed directly to the integration routine, HCubature.hcubature.

See also: cumradialprofile, mass_sqarcsec.

fits_lensmodel(filename, model, pars, xrange, yange, [z=1.0]; planes)

Write results of a model, computed in a grid, in a FITS file.

This procedure can be used to create FITS files containing the results of a given lensing model, computed at a given redshift. The quantities computed and written to the file can be specified with the planes keyword, which must be a list of symbols. The valid planes are

• :massdensity: the total mass density Σ in units of M⊙ / pc² (see massdensity).
• :lensed_massdensity: the total mass density Σ in units of M⊙ / pc², taking into account lensing effects (see lensed_massdensity).
• :mass_sqarcsec: the total mass density Σ in units of M⊙ / arcsec² (see mass_sqarcsec).
• :lensed_mass_sqarcsec: the total mass density Σ in units of M⊙ / arcsec², taking into account lensing effects (see lensed_massdensity).
• :magnification: the inverse of the magnification, i.e. μ⁻¹ (see magnification)
• :timedelay: the time delay τ, in days (see timedelay)
• :lensmap: the lens map, stored as two planes, β₁ along x and β₂ along y (see lensmap)
• :distortion: the lens distortion, as a (possibly non-symmetrical) 2×2 matrix, stored in four planes J₁₁, J₁₂, J₂₁, and J₂₂ (see distortion)

All these planes are computed at the redshift z specified. Additionally, it is possible to compute lens quantities associated to the last lensing plane before z:

• :convergence: the plane dimensionless convergence κ (see convergence)
• :potential: the lensing potential ψ̂, see potential
• :deflection: the deflection angle in the angular unit of the model, which by default is 1", stored as two planes as α̂₁ along x and α̂₂ along y (see deflection)
• :jacobian: the jacobian of the deflection angle, as a symmetric 2×2 matrix, stored as three planes as A₁₁, A₁₂, and A₂₂ (see jacobian)

The default is planes = [:massdensity, :magnification, :timedelay, :deflection].

The plane names will be saved in the FITS header under the keywords PLANEn, where n is a positive integer: therefore, the default configuration will results in a series of keywords PLANE1 = "massdensity", ... PLANE4 = "deflection1", PLANE5 = "deflection2".

fits_linear_sources(filename, model, pars; filter=nothing)

Write the linear sources of a model to a FITS file.

This procedure can be used to create a FITS file containing the sum of all linear sources defined in the model, evaluated at the parameters pars.

The sources are summed together, taking into account their grids and offsets. Therefore, the source must have compatible grids: i.e., the same pixel sizes, and pixels that are aligned on a common grid. This can be easily obtained by using a snapped_dx_dy geometry when defining sources in the model. If some sources have incompatible grids, they will be skipped with a warning.

The filter keyword can be used to specify a filter for the sources to be included in the sum: only sources whose name matches the filter will be included. The filter can be a string, a regular expression, or a function that takes a source name as input and returns true if the source must be included, and false otherwise. If filter is nothing, all sources are included.

lensing_area(lenssys; factor=2, npts=10_000)

Compute the bounding box of a lensing system.

The computation takes into account the positions and einstein radii of all lenses. The area is finally enlarged by factor.

It returns a tuple with two ranges. The total number of points is kept close to npts.

fiducialposition(sys, z, x₁ [, x₂...])

Return the fiducial position of a source at redshift z observed at xᵢ.

The first argument, sys, is a LensSystem. Each of the xᵢ must be either a multivariate distribution, or a tuple with the coordinates of the image position. The function returns the fiducial position as a SPoint.

fiducialposition(sys, z, β, x₁ [, x₂...])

Return a parametrized fiducial position of a source.

The first argument, sys, is a LensSystem, while the second argument is the source redshift z. Each of the xᵢ must be either a distribution obtained, for example, with a syntax such as (1.0, 2.0) ± 0.1.

The β argument is a SPoint representing the deviation of the source from the mean position, in units of the source covariance. This function is used internally to convert a configuration file fiducialposition entry into a source position.

counterimages(sys, z, x₁ [, x₂...]; kw...)

Find all counterimages of a source at redshift z observed at xᵢ.

The first argument, sys, is a LensSystem. Each of the xᵢ must be either a multivariate distribution, or a tuple with the coordinates of the image position. The function returns a named tuple with a number of data for all counterimages of a source at redshift z observed at x₁, x₂... In particular, the returned named tuple contains the following field names:

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

If the keyword argument grid=:auto (the default), the procedure will find a suitable grid for the inversion of the lens equation using the function lensing_area. In this case the criticalcurves parameter controls the maximum level of refinement of the grid around the critical curves (default=4). The images and lenses similarly control the refinement level around the images and the lenses centers (default=3).

Alternatively, one can pass directly a Gravity.MaskedGrid object for grid.

If include_source=true, the output will include, as first element, the source position; the fields Δmag and t will be set to zero.

Example

julia> lp, pars = Gravity.bestlp(chain);

julia> Gravity.counterimages(model.make_lenssys(pars), 1.0, (3.14, 2.72), (4.78, -3.29))

julia> Gravity.counterimages(model.make_lenssys(pars), 1.0, (3.14, 2.72) ± 0.05, (4.78, -3.29) ± 0.08)