Lens inversion

The inversion of the non-linear lens equation can be a non-trivial task and, often, finding all possible solutions can be difficult.

This problem is tackled in Gravity using two general strategies: a simple inversion that starts from an initial position close to the solution (for example, an observed image), and a more robust inversion that tries to find all solutions using information of the lens equation at the corners of a Gravity.MaskedGrid.

In both cases, iterative algorithms are used to solve the non-linear lensing equation. The implemented algorithms are

Newton's method, based on the use of the local derivative of the lensing equation. It is generally quite fast and robust.
Broyden's methood, a quasi-Newton method based on rank-one updates of the Jacobian matrix. It does not use therefore derivatives of the lensing equation.
Klement's method, a modification of Broyden's method, also not using derivatives.

Gravity.ConvergenceError — Type

An error raised by non-linear equation solvers in case of missing convergence.

source

Gravity.lensunmap — Function

lensunmap(lenssys, ς, β [, θ₀ [, lensinversion])
lensunmap(lenssys, ς, β [, θ₀]; lensinversion=default_lens_inversion)

Perform a lens inversion and find counter-images of the source position β.

The inversion uses a cascade of methods, as specified by the lens_inversion keyword: if one method fails (because it reaches the maximum alloowed number of iterations), the next metthod is picked. The starting point used is θ₀ (by default taken to be equal to β).

Method-specific keywords parameters can be provided using the lens_inversion parameter. That must be an NamedTuple with the format method = flags, where the method can be one of :newton, :broyden, or :klement, and the flags are the method-specific flags to use. The special method :success can be used to avoid the raise of a ConvergenceError: if specified, the solution returned is the last point of the last algorithm.

The recognized parameters for each method are the following:

trustregion: the squared size of the trust region, i.e. the maximum ratio between size of the position correction (|δx|) and the error on the source plane (|δy|), squared: $|\delta x|^2 \le \tt{trustregion} \cdot |\delta y|^2$. If set to nothing, no trust region is used.
ftol: the tolerance on the squared distance between the mapped point and β. If set to nothing, this tolerance is ignored.
xtol: the tolerance on the squared distance between two iterations. If set to nothing, this tolerance is ignored.
maxiters: the maximum number of iterations.

source

lensunmap(lenssys, ς, β, [θs,] grid [, mapped_corners]; keywords)

Find the counter-images of a source position β in a gridded lens system.

The solutions are found by iterating all cells in grid and by finding cells that, when mapped through the lens equation, contain the point β. In case the mapped grid corners are already available, they can be passed in mapped_corners; if not, they are (effectively) computed as mapped_corners = lensmap.(sys, ς, grid.corners).

The function returns an array of SPoint values, each corresponding to one solution of the lensing equation (see also lensunmapdistortion).

Any additional keyword present is directly passed to lensunmapdistortion.

source

lensunmap(lenssys, ς, source [; kw])
lensunmap(lenssys, ς, source, (image|θ) [; kw])
lensunmap(lenssys, ς, source, images, grid [, mapped_corners] [; kw])

Compute a counter-image corresponding to a given source.

This method inverts the lens equation for a LensSystem lenssys, with associated SourceCoeff information in ς. source can be any Gravity.SingleSource, such as a Gravity.PointSource or a Gravity.LuminousSource, or any Gravity.ImageFamily.

If provided, image must be an image measurement corresponding to the given source: so, for example, if source is a Gravity.PointSource, image should be a Gravity.PointImageMeasurement. For SingleImageMeasurement's, one can alternatively pass the image position, i.e. θ = point(image).

In all cases the returned value will be of the same type of source.

Alternatively, image (or θ) can be an vectors, representing the approximate positions of various counter-images. The returned value in this case will be an array of images. Note that no effort is taken to make sure that two returned images are not identical.

As a third possibility, one can pass a AbstractGrid grid and optionally its mapped corners: in this case the approximate positions to seed the inversion are computed from the mapped grid cells. The returned value is an array of sources.

In any case, all optional keywords kw are passed to the inversion routine (lensunmap for points).

source

lensunmap(lenssys, source_plane [; kw])
lensunmap(lenssys, source_plane, θs [; kw])
lensunmap(lenssys, source_plane, image_plane [; kw])
lensunmap(lenssys, source_plane, grid [, mapped_corners] [; kw])

Map a source-plane into an image-plane.

This method calls the appropriate lens-unmapping routines for all images in the plane and collects the results into an image-plane. Various versions of this routine are avaiable: they differ in the way the seeds for the inversion of the lens equation are obtained.

If no starting position are available (first line above), the inversion just starts using the source positions as seeds. The output will be an image-plane, that is a SourcePlane having ImageFamily'ies consisting of single images.
Optionally, one can pass an array of SPoint's θs: this will be used as starting position for each family. A simple array will result in ImageFamily'ies consisting of single images. If, instead, θs is an array of arrays of SPoint's, these will be interpreted as different seeds for different images of the various ImageFamily'ies present, and the output will be correspondingly rich.
The same effect can be obtained by passing a full image_plane: if the positions of the various images there will be used as seeds for the inversion.
Finally, one can pass a AbstractGrid grid, with optional mapped_corners to perform a grid-inversion of the lens equation.

Depending on the value of Gravity.use_threads, the lensing operations for the various sources will be performed using threads.

source

Gravity.lensunmapdistortion — Function

lensunmapdistortion(lenssys, ς, β [, θ₀ [, lens_inversion]])
lensunmapdistortion(lenssys, ς, β [, θ₀]; lens_inversion=default_lens_inversion)

Perform a lens inversion and find counter-images of the source position β.

The result is a PointDiff structure providing both the image position and the Jacobian of the ray-tracing equation at the found position (but see below).

The recognized parameters for each method are the following:

trustregion: the squared size of the trust region, i.e. the maximum ratio between size of the position correction (|δx|) and the error on the source plane (|δy|), squared: $|\delta x|^2 \le \tt{trustregion} \cdot |\delta y|^2$. If set to nothing, no trust region is used.
ftol: the tolerance on the squared distance between the mapped point and β. If set to nothing, this tolerance is ignored.
xtol: the tolerance on the squared distance between two iterations. If set to nothing, this tolerance is ignored.
maxiters: the maximum number of iterations.

source

lensunmapdistortion(lenssys, ς, β [, θs] [, grid [, mapped_corners]]; keywords)

Find the counter-images of the source position(s) β in a lens system.

The solutions are found by iterating all cells in grid and by finding cells that, when mapped through the lens equation, contain the point β. In case the mapped grid corners are already available, they can be passed in mapped_corners; if not, they are computed as mapped_corners = lensmap(lenssys, ς, grid.corners).

Duplicate solutions will be removed: by default, the tolerance (squared distance) for duplicate solutions will be equal to 16 * xtol.

The function returns an array of SPointDiff values, each corresponding to one solution of the lensing equation.

The keyword parameters fix the stopping criteria: the function stops searching for a solution when the mapped point is closer than sqrt(ftol) to β, or when the correction between two iterations is smaller than sqrt(xtol).

Method-specific keywords parameters can be provided using the lens_inversion keyword. That must be a NamedTuple with the format :method = flags, where the method can be one of :newton, :broyden, or :klement, and the flags are the method-specific flags to use. The special method :success can be used to avoid the raise of a ConvergenceError: if specified, the solution returned is the last point of the last algorithm.

Keywords

lens_inversion: method-specific keywords parameters (see above)
xtol: the tolerance on the squared distance between two solutions. Default is taken the smaller xtol set for among the fitting methods.
failsafe: this keyword is used to define some critical cases. If it is :always (the default), the algorithm will always perform a first search of an image position using the source position β as starting point. If it is :maybe, this action will be carried out only in case no image is found within the grid. These settings will generally avoid the case where no predicted image is found, which will in turn generally lead to complications in the likelihood estimates.
matches: if :all (the default), all matches founds are returned; if :best, the closest solution to each point of θs is returned.

source