Pipeline

The usual way to perform a strong lensing analysis using Gravity is through the use of a series of functions. The functions involved had some special features:

When executed, they will generate Julia code which, when compiled, will run at full speed, thus making the analysis as fast as possible.
They take out from the end user most of the burden associated with the use of a complex gravitational lensing library such as Gravity.
Still, they retain a lot of flexibility, allowing an advanced user to ingest Julia expressions in the analysis, such as arbitrary scaling relations or non-standard parameters.

The pipeline is built upon two main functions:

gravitymodel is used to load a configuration (usually, but not necessarily, written in a YAML file). This function is also responsible for the validation of the configuration and for the definition of all the main building blocks needed for the analysis (such as the log-likelihood and many other associated functions).
runmodel performs the actual analysis following the specifications provided in the setup file. The output, for all Bayesian algorithms, will be generally provided following the MCMCChains specification.

Additionally, Gravity.jl also provide some ancillary functions, such as imageplainchain, which can be used to improve the convergence of many algorithms (see below).

One-pass pipeline

The simplest combination of commands involves the use of the two functions briefly described above:

julia> using Gravity

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

julia> chain = runmodel(model)

During the analysis performed by the last command, some information will be displayed (in a format that depends on the algorithm used for the parameter inference). The result will be saved in the chain variable: for Bayesian Monte Carlo algorithms, chain will be an MCMC-Chain. As a result, one can obtain various relevant statistical quantities on the involved variables in simple ways:

julia> mean(chain)

julia> quantile(chain)

Moreover, it is possible to obtain simple visualizations of the chains using the techniques described in the MCMCChains documentation.

Two-pass pipeline

Many inference algorithms can benefit from some approximate estimates of the inference parameters. In such cases, it might be wise to perform a two-step analysis

In the first step, the inference is carried out in the source plane (using sampling.scheme: sourceplane in the configuration), thus reducing significantly the number of free parameters and avoiding the inversion of the lens equation.
In the second step, the parameters values deduced in the first step are used as starting points for an analysis in the image plane (sampling.scheme: imageplane).

This two-pass pipeline can be easily adopted by specifying sampling.scheme: sourceplane in the configuration file and by adding an initialization keyword in the sampling section. This can be obtained for example from the following series of commands:

julia> using Gravity

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

julia> chain_sp = runmodel(model_sp);

julia> model_ip = gravitymodel("model.yml"; scheme="simplified-imageplane", initialization=chain_sp);

julia> chain_ip = runmodel(model_ip);

Since this series of commands is quite common, the model_ip definition can also be written more concisely as

julia> model_ip = simplifiedimageplane_model(chain_sp);

Note that in this case the configuration will be automatically obtained from the chain and modified appropriately.

Interface

Below we describe the main functions used in the pipeline. For a description of the internals, please see the dedicate page.

Gravity.gravitymodel — Function

gravitymodel()
model = gravitymodel("filename.yaml"; kwargs...)
model = gravitymodel(dictionary; kwargs...)
model = gravitymodel(chain; kwargs...)

Load a full gravitational lensing model stored in a filename (entered as A plain string), in a multiline string, or in a dictionary. The filename or the multiline string must be entered using the YAML markup language.

Alternatively, one can pass an (ordered) dictionary with the full model, as in the GravitySetup constructors. This can be useful if one needs to modify programmatically a model.

As a third option, one can create a new model from a previous MCMC chain.

In all cases, the function accepts keywords parameters which can be used to modify in-place the setup. The names reflect the recognized keywords in the sampling or sampling.likelihood-options sections of the YAML file: scheme, algorithm, method, bayesianfactor, matches, duplicates, missedpredictionspenalty, grid, gridrefinemeent_criticalcurves, gridrefinement_images, gridrefinement_lenses, rounds, chains, variationalrounds, variationalchains, explorer, points, iterations, samples, timelimit, threads, copies, randomseed, and initialization. A few of these keywords can be set to nothing or to missing to indicate that the corresponding setup parameter is unspecified: this is convenient for keywords such as grid or randomseed, for example.

The special parameters keyword can be set to additional parameter definitions to be included in the configuration file. Must be a dictionary with strings as keys: for example parameters=Dict("lens_redshift" => 0.78).

The special debug keyword can be used to show the code produced by this function. It accepts a symbol or a list of symbols, and will show the corresponding defitions: for example, one could pass debug=:loglikelihood to display the definition of the loglikelihood function. Use debug=:all to show all definitions (see the extended help). Since the output is often very long, you 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.

If called without any parameter, the function generates an example of YAML configuration. The displayed value of this configuration uses different colors for clarity. Use repl(gravitymodel()) to obtain a simple string representation of it.

Example

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

julia> result = runmodel(model)

Extended help

The function generates the code to evaluate the prior distribution of the parameters using [setup_code_prior] and the log-posterior using [setup_code_likelihood]. The generated code will be included in a module and will include, among other things, the following function definitions:

PriorDistribution(): a multivariate, continuous distribution, representing the prior; it has methods such rand, logpdf, and quantile defined. For convenience the constant prior is defined as the result of this call.
ReferenceDistribution(): a multivariate, continuous distribution, representing the reference distribution, i.e. the initial guesses of the parameters; it has methods such rand, logpdf, and quantile defined. For convenience the constant reference is defined as the result of this call.
likelihood_options(): returns the constant parameters to pass to the several functions (see below); for convenience the constant options is defined as the result of this call.
loglikelihood(x, options=likelihood_options()): computes the log of the likelihood for a given vector of parameters x.
logposterior(x, options=likelihood_options()): computes the log of the posterior for a given vector of parameters x.
make_cosmo(x, options=likelihood_options())): returns the cosmological model implied by the parameter vector x.
make_lenssys(x, options=likelihood_options()): returns the LensSystem associated the parameters x.
make_sources(x, options=likelihood_options(), sys=make_lensys(x, options)): returns the sources associated to the parameters x.
make_images(x, options=likelihood_options(), sys=make_lensys(x, options)): returns the images associated to the parameters x.
make_predictions(x, options=likelihood_options(), sys=make_lensys(x, options)): returns the predicted images associated to the parameters x.

Help for the individual functions will be available.

source

Gravity.runmodel — Function

result = runmodel(model [, initial_state]; show=true)

Perform model inference for the given model.

The parameter model must have been obtained using the function gravitymodel. The result produced depends on the specific algorithm and method used in the lensing model inference.

The optional parameter initial_state can be used to provide an initial state for the inference: this is used internally by extendchain.

The keyword show controls the display of informative messages.

This functions saves in result.info a number of relevant metadata, including

model: the input model
cputime: the total CPU time used (summed over all threads)
start_time and stop_time: the start- and stop-time of the execution
sysinfo: a collection of various system information data

Other metadata are specific to the algorithm used.

Info

runmodel, together with gravitymodel, are the two main user interface routines provided by Gravity. For MCMC algorithms, it creates an MCMCChains.Chains object which can then be plotted using various packages (such as StatsPlots.jl or PairPlots.jl).

source

Gravity.savereport — Function

savereport(path, chain [, model] [, comment]; kw...)

Generate a report of the MCMC execution and save it to the specified path.

The report includes best-fit results, model plots, MCMC diagnostics, and the configuration used for the Gravity model. The report can be saved in Markdown and/or HTML formats. The comment argument can be used to add a custom comment at the beginning of the report (it must be a valid Markdown string).

If Makie is available, the report also includes plots of the best-fit model. Moreover, if PairPlots is available, a corner plot of the MCMC samples is included. The plots are saved in the requested formats (i.e., SVG or PNG, as requested by the figure_formats keyword argument) in the assets subdirectory. Note, however, that only one of these formats will be embedded in the HTML report (SVG preferred over PNG).

By default, this function will also generate an Aladin-lite plot of the best-fit model, useful to visualize the lensing configuration. The Aladin-lite plot, if requested (:aladin in figure_formats), is always embedded in the HTML report, with no external dependency. Moreover, if :aladin is included in the formats keyword argument, the Aladin-lite plot is also saved as a standalone HTML file in the main report directory.

Arguments

path: the directory where the report will be saved
chain: the MCMC chain object containing the sampled parameters and metadata
model: (optional) the Gravity model used for the MCMC sampling (useful to update the report if some parameters not affecting the chain have been modified); if not provided, the model will be extracted from chain.info.model
comment: (optional) a custom comment (Markdown string) to include at the beginning of the report

Keyword Arguments

open_in_browser: if true, open the HTML report in the default web browser after generation (default: true)
auto_rename: if true, automatically rename parameters for better readability (default: true)
digits: number of significant digits to display (default: 4)
formats: tuple specifying the formats to save the report; the accepted values are :md, :html, :latex, and :aladin (default: (:md, :html, :latex, :aladin))
figure_formats: tuple specifying the formats for figures; the accepted values are:svg, :png, :aladin (default: (:svg, :png, :aladin))
savechain: if true, save the MCMC chain as a CSV file in the report directory (default: false)

source

Gravity.ptreport — Function

ptreport(chain, index=1)

Display a report for a MCMC chain.

Currently this function is only implemented for the pt algorithm. It shows a report that is similar to the one produced during the execution of the algorithm.

The index parameter is used to select the chain in the case of a multi-chain run.

source

Gravity.simplifiedimageplane_model — Function

newmodel = simplifiedimageplane_model(chain; kwargs...)

Create a new gravity model from a previous MCMC chain.

The new model will be initialized with the parameters from the last iteration and will have scheme: simplified-imageplane.

The optional keyword arguments are passed to gravitymodel and can be used to set additional sampling parameters.

source

Gravity.imageplane_model — Function

newmodel = imageplane_model(chain; kwargs...)

Create a new gravity model from a previous MCMC chain.

The new model will be initialized with the parameters from the last iteration and will have scheme: imageplane.

The optional keyword arguments are passed to gravitymodel and can be used to set additional sampling parameters.

source

Gravity.imageplanechain — Function

newchain = imageplanechain(sourceplanechain; xonly=true)

Extend an MCMC obtained with scheme: sourceplane to include source parameters

This function computes the source configuration for each combination of parameters of the chain. It returns a chain that includes all original parameters plus the parameters assocuated to each source.

If xonly=true, the resulting chain will only include the x coordinate of each source. This is useful if one performs an image-plane analysis without scheme: extended-imageplane as a likelihood option (the default). Otherwise, the chain will include all the parameters of each source.

See also: imageplanechain

source

Gravity.extendchain — Function

newchain = extendchain(oldchain)

Rerun an MCMC chain with a larger number of iterations.

This function extends an oldchain with more iterations; it uses the same algorithm used in the original run. The new chain will double the number of points of the original chain.

source

Gravity.loadmodel — Function

dict = loadmodel(filename; gmodule=Module(:GravityModel))

Load a model from a YAML file and and returns an OrderedDict.

This function effectly calls YAML.load_file with the appropriate dicttype parameter (an OrderedDict), and then expands loop: sections appropriately.

The optional keyword can be used to specify the module where the model will be created.

source

Gravity.GravitySetup — Type

GravitySetup

A full setup for gravitational lensing problem.

Members

path: the path to the source YAML file, if any, or nothing
ctime: the creation time of the YAML file, corresponding to the modification time of the source file
mtime: the modification time of the YAML file
rawsource: the original source of the setup, as a string
kwargs: the keyword arguments passed to gravitymodel to create the setup; together with rawsource, this is used to recreate the setup when loading a chain
source: an OrderedDict with the full setup, with implicit parameters made explicit
preamble: an OrderedDict with a list of variable definitions to use in the preamble
parameters: an OrderedDict with the free parameters as keys and their dimensionality as values
source_preamble: an OrderedDict with a list of source-related variable definitions to use in the preamble
source_parameters: an OrderedDict with the source-related free parameters as keys and their dimensionality as values; source paramers explicitly fit (thus appearing also in parameters) are included with negative values
samples: an OrderedDict with free and constant parameters as keys and a vector of samples (generated using the prior) as values
data: an OrderedDict with the list of data (FITS files), preloaded for speed

!!! warning The source_preamble and source_parameters field are initialized only for sourceplane and simplified-sourceplane schemas. Their main use is to keep track of source-related free parameters (such as the source positions) which needs to be ignored for sourceplane inference. This task is important to be able to generate a corresponding imageplane setup. Note that an exception is given by the source redshift: this quantity is relevant for sourceplane inference too, and therefore is saved in the preamble and parameters field in all cases.

Constructors

GravitySetup(something [; rawsource=nothig, gmodule=Main, filename=nothing,
             kw...])

Loads a filename, a string in the YAML markup language, or directly a dictionary, validates it using setup_validation, and parses it with setup_parse to produce the GravitySetup structure.

The optional rawsource parameter is just kept identical in the returned GravitySetup structure, and is used to recreate the setup when loading a chain. The optional gmodule parameter is the module where the model will be created. The optional filename parameter is the name of the file from which the setup was loaded, if any.

All other parameters are used to modify in-place the setup. The names reflect the recognized keywords in the sampling or sampling.likelihood-options sections of the YAML file: parameters, lens_lenses, scheme, algorithm, method, bayesianfactor, matches, duplicates, missedpredictionspenalty, grid, gridrefinement_criticalcurves, gridrefinement_images, gridrefinement_lenses, rounds, chains, variationalrounds, variationalchains, points, iterations, timelimit, grid, threads, copies, randomseed, and initialization. Additionally, this function recognizes the lens_lenses and use_threads special parameters.

The initialization keyword is treated in a special way. If it is an MCMCChains, it its value will be saved in a variable in the generated model with the name mcmc_init, and this variable will be used as initialization for the sampling.

source