grale.util

Module meant for various utilities, but currently only RMS calculation.

grale.util.calculateImagePredictions(imgList, lensModel, cosmology=None, reduceImages='average', useAverageBeta=True, maxPermSize=7, useFSolve=True)

For a set of images and a lens model, the predicted images for each observed image are calculated. The results can subsequently be used in the calculateRMS() function to obtain the RMS.

Parameters:

imgList: a list of dictionaries with entries
- "imgdata": ImagesData instance that describes the observed images.
- "z": the redshift for these observed images.
lensModel: this can be a MultiLensPlane instance, a MultiImagePlane instance, a list of (lens model, redshift) tuples, a LensPlane <grale.images.LensPlane>` instance, or just a single model.

In case only one or more lens models are specified, the predicted image positions are estimated by calculating the derivatives of \(\vec{\beta}(\vec{\theta})\) and using them to compensate for small differences in source plane positions \(\vec{\beta}\); or if useFSolve is True, several iterations will be used to approximate the true solution using scipy’s fsolve routine.

On the other hand, if e.g. a MultiLensPlane instance is used, then the traceBetaApproximately function is called to obtain estimates of image plane positions corresponding to a source plane position. This is more computationally demanding, but should yield a more correct result.
cosmology: must only be specified somehow (e.g. “default” if a default cosmology was set if a (list of) lens models was specified. If e.g. a MultiLensPlane was used, the internally stored cosmological model will be used instead.
reduceimages: each image in an images data set will be converted to a single point. By default, this is done by averaging the image point positions. If something else needs to be done, you can specify a function here, which will be called with two arguments: the ImagesData instance and the image index corresponding to the particular image.
useAverageBeta: This can be either a boolean or a list. If True, then the back-projected image points are averaged to estimate a single source position. If False, then each back-projected image is used as an estimate of the source position. If this is a list, the number of entries should match imgList, and each entry is used as the source position for that imgList entry.
maxPermSize: when a source plane position is used to estimate the corresponding image plane positions, these predicted positions are grouped with the observed positions. If the derivatives-based method is used, then this is automatically possible, but if the more accurate tracing is used, we need to find out which predicted point corresponds to which observed point. As a first attempt, the procedure tries to use the points closest to the observed ones, but this may not be possible and different permutations of the positions will then be explored. Since this can become quite slow, an error is generated if more elements than this would need to be permutated.
useFSolve: if only a lens model is used, this cause SciPy’s fsolve to be used to zoom in on the real image plane vectors. Otherwise, a one-step approach is used based on the deflection angle derivatives.

grale.util.calculateRMS(predictions, angularUnit, avgImage=False)

Using the output of calculateImagePredictions(), the RMS value is calculated. Using the notation from Appendix A of Williams et al (2018), in general there are \(i=1..I\) sources, with \(j=1..J_i\) images each. When using each back-projected observed image position \(\vec{\theta}_{i,j}\) as a separate estimate of the source position, there will be \(k=1..J_i\) image predictions \(\vec{\theta}_{i,j,k}\) for each observed position.

When treating each predicted-observed difference of all sources equally, the RMS will be given by

\[{\rm RMS}_{\rm full}^2 = \frac{1}{K}\sum_{i=1}^I \sum_{j=1}^{J_i} \sum_{k=1}^{J_i} \left|\vec{\theta}_{i,j,k} - \vec{\theta}_{i,j}\right|^2\]

where \(K = \sum_{i=1}^{I} J_i^2\) is the total number of terms in this summation. Alternatively, we could take the average of the right-most sum, thereby averaging the squared differences for each observed image point. Sources with more images will still have more terms in the rest of the summation, but only one per image (and not \(J_i\) for each image). This results in

\[{\rm RMS}_{\rm equalimages}^2 = \frac{1}{J}\sum_{i=1}^I \sum_{j=1}^{J_i} \frac{1}{J_i} \sum_{k=1}^{J_i} \left|\vec{\theta}_{i,j,k} - \vec{\theta}_{i,j}\right|^2\]

Here, \(J = \sum_{i=1}^I J_i\). We could also average all \(J_i^2\) terms per source first, and then average over the remaining \(I\) terms, one for each source:

\[{\rm RMS}_{\rm equalsources}^2 = \frac{1}{I}\sum_{i=1}^I \frac{1}{J_i^2} \sum_{j=1}^{J_i} \sum_{k=1}^{J_i} \left|\vec{\theta}_{i,j,k} - \vec{\theta}_{i,j}\right|^2\]

This function returns a dictionary with these three RMS values (non-squared), expressed in the specified angularUnit.

If avgImage is True, then all image predictions for each observed position are averaged out first. This will only make sense in case useAverageBeta was False in the call to calculateImagePredictions() since otherwise there’s only one predicted position.

Depending on the settings of such flags, and depending on the use of the derivatives-based prediction or the more accurate tracing, some of these results can yield the same value.

grale.util.createMonopoleBasisFunctions(avoidSources, Dd, subDiv, size, center=[0, 0], widthFactor=3.0, rangeFactor=4.0, centralDensity=1.0, overlapNeeded=None, randomOffset=True, cellCenterCallback=None, cellCenterCallbackState=None)

This creates a set of the monopole basis functions from 2008MNRAS.389..415L with which mass can be redistributed in between the images in a lensing system. The use from this article is illustrated in the example_monopoledegen example.

It is also used in a different way in 2012MNRAS.425.1772L, where the deflection field at some image points is modified. This kind of use is demonstrated in example example_gen_msd.

Arguments:

avoidSources: a list of ImagesData instances, describing the images to avoid. None of the generated monopole basis functions will overlap with these images.
Dd: the angular diameter distance to the lens plane, to be used in these basis functions.
subDiv, size, center: the square region of size size, centered on center, will be subdivided into subDiv cells along each axis.
widthFactor: if a monopole can be added for a certain grid cell, the width of the positive part will be this factor times the cell width.
rangeFactor: for each cell, the minimal distance to the images will be determined. Only if this is larger than this factor times the cell width, will a monopole be added.
centralDensity: the monopole basis function will have this central density.
overlapNeeded: in the case that we’re trying to modify properties at some points, we need to make sure that a basis function actually covers these points, otherwise it will not have any effect. This parameter is a set of ImagesData instances, with which overlap is required.
randomOffset: if set to True (the default), the specified center will not be used exactly, but some randomness will be added.
cellCenterCallback: if set, for each cell of the resulting grid this callback function will be called. It will be passed four parameters: the center of the cell, the size of the cell, a boolean indicating if a monopole basis function could be added for this position, and the cellCenterCallbackState argument that’s specified next.
cellCenterCallbackState: when the previous callback function is set, this value will be used as it’s last argument.

grale.util.createThetaGrid(bottomLeft, topRight, NX, NY): This creates a grid of xy-position, from bottomLeft to topRight, where the x-direction is evenly divided into NX points, and the y-direction into NY points. The result is a numpy array of shape (NY, NX, 2)

grale.util.createThetaGridAndImagesMask(bottomLeft, topRight, NX, NY, maskRegions, enlargements=2, enlargeDiagonally=False, circleToPolygonPoints=10000)

The bottomLeft, topRight, NX and NY parameters are passed to createThetaGrid. Then, for this grid a mask is created based on the other arguments.

The maskRegions argument is a list of region descriptions that are combined into one final mask for the grid, a boolean numpy array of shape (NY, NX). The entries of the list can be the following:

an ImagesData instance: the grid positions covered by the images in this instance are set to True in the mask.

a 2D boolean array, also of shape (NY, NX): the positions where this mask is True are also set in the final mask.

a dictionary, with at least a key called type, which can take on the following values:

hull: in this case, a key imgdata should be present as well. It should contain either a single ImagesData instance, or a list of these instances. All points will be collected, and the convex hull of this region will be used.

circle: for this type, center and radius should be present as well. Internally this will be converted into a polygon, of which the number of points can be controlled by the circleToPolygonPoints argument.

rect: describes a rectangle, and bottomleft and topright properties should also be present.

square: for this type, center and size should also be present.

polygon: this describes a polygon, there should be a coord property that represents a list of points that make up the polygon.

point: for this type, a coord property describes a single point.

For each of these types, an invert key may be present as well. If set to True, then the region will be inverted.

The enlargements value determines how many times the mask should grow. If a value is True, then the values to the left, right, below and above will be set to True as well. If enlargeDiagonally is set to True, the four diagonal values will also be set to True. This procedure is repeated the number of times specified by enlargements.

The return value is a tuple of the resulting grid, and the final mask.

grale.util.findOptimizedSourcePositions(imgList, lensModel, cosmology=None, reduceImages='average', optRoutine=<function nelderMeadSourcePositionOptimizer>, optParams={})

For a given lens model, which typically will not back-project the images of a source onto the exact same position, this routine looks for a ‘best’ source position.

This function returns a list of optimized source positions, one for each entry in imgList.

Arguments:

imgList: a list of ImagesData instances describing the images in the strong lensing system. Each image will be reduced to a point image position, see also the reduceImages argument.
lensModel: the lens model, which can be represented as a number of things:
- a GravitationalLens instance. In this case, the cosmological model cosmology must be set as well.
- a :class: LensPlane <grale.images.LensPlane> instance. In this case as well, the cosmology parameter must be set.
- a MultiLensPlane or MultiImagePlane instance. The cosmology parameter must be None, the interally set cosmology will be used.
- a list of (GravitationalLens, redshift) tuples, where the cosmology parameter must be set as well.
cosmology: depending on the way the lens model is specified (see lensModel argument), this should either be None or set to the cosmological model to use.
reduceImages: each image will need to be reduced to a point image position. If set to "average", then the average position of the image points will be used for this. This can also be set to a function, and if so it will be called with a ImagesData instance as parameter and an index of the image. The function should then return a single point that corresponds to the specified image.
optRoutine: this is the optimization routine used to calculate the optimal source position for a set of image positions. The default is the nelderMeadSourcePositionOptimizer(). The routine is called with the following arguments: the index into a list of all images (this is useful for the parallel version parallelFindOptimizedSourcePositions()), the point image positions for a single source point (a dictionary containing theta, beta, z and depending on the input model also betaderivs and invbetaderivs), the function that’s created based on the input model to trace an image plane point to its source plane position, and the expanded dictionary optParams from below.
optParams: this can be used to pass additional parameters to the optRoutine.

grale.util.parallelFindOptimizedSourcePositions(imgList, lensModel, cosmology=None, reduceImages='average', optRoutine=<function nelderMeadSourcePositionOptimizer>, optParams={}, numThreads=0)

This is a version of findOptimizedSourcePositions() that tries to speed up the procedure using multiple threads. It handles several sources in parallel.

Arguments:

imgList, lensModel, cosmology, reduceImages, optRoutine and optParams have the same meaning as in findOptimizedSourcePositions()
numThreads: the number of threads to use to run in parallel. If set to zero or None, this will be determined automatically.

grale.util.nelderMeadSourcePositionOptimizer(imgIdx, imgPos, traceFunction, feedbackObject=None, betaAvgToleranceArcsec=1e-05, betaAvgPenalty=10000, nmRounds=3, nmMaxFev=50)

This is the default optimizer used by findOptimizedSourcePositions() and parallelFindOptimizedSourcePositions(). It uses the Nelder-Mead method as implemented in scipy.

Arguments:

imgIdx, imgPos, traceFunction: as described in the optRoutine argument of findOptimizedSourcePositions().
feedbackObject: can be used to specify a particular feedback mechanism
betaAvgToleranceArcsec: depending on the lens model, it is possible that for a specific source plane position beta, there simply isn’t an image plane position near the one considered that maps well to this. To filter out such beta positions, there is this threshold that specifies that if the back-projected (adjusted) image plane position still differs more than betaAvgToleranceArcsec (so in arc seconds) from beta, then a penalty should be added.
betaAvgPenalty: in the case that a penalty (described above) needs to be added, so the source plane difference exceeds betaAvgToleranceArcsec, this difference times betaAvgPenalty is added to the function evaluation result.
nmRounds: the Nelder-Mead routine can be run a number of times, and the best result will be kept. This argument specifies that number.
nmMaxFev: each Nelder-Mead run can have at most this many function evaluations.