grale.util

Module meant for various utilities.

grale.util.calculateImagePredictions(imgList, lensModel, cosmology=None, reduceImages='average', reduceSources='average', maxPermSize=7, localTraceFunction='fsolve', localTraceFunctionOptions=None, useFSolve=None, useAverageBeta=None)

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 instance, or just a single model.

    Depending on the type of lensModel, a different kind of algorithm will be used. If it’s based on a pixelated version of the lens plane, a global approach can be used to calculated the predicted image positions: based on the back-projected pixel locations and the estimated source position, a number of image positions can be recovered and will be matched to the input image positions.

    If only a lens model is provided, a more local approach is used: for a source position that’s estimated first, an optimization procedure is used, specified in localTraceFunction. For each observed image position, this function is called, which can then adjust the starting position (observed one) so that, when back-projected, it coincides with the source position.

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

  • reduceSources: the point images that are either input or result from the previous reduceImages step, are projected onto the source plane and a source position is estimated from these back-projected positions. The value of this parameter determines how this is done:

    • "average": the back-projected positions are averaged.

    • "magweighted": a weighted average of the back-projected points is used, where the weight corresponds to the magnification calculated at the image position.

    • "noreduction": no reduction to a single source position occurs. Instead, each of the back-projected image positions is used as a possible source position.

    • a list is specified: this is considered to be a list of source positions, so this list must have the same length as the imgList.

    • a function is specified: to reduce the back-projected positions this function is called with three arguments: (imgPlaneObject, srcIdx, imgInfo). The first is an image plane-like object, for which traceTheta and getBetaAndDerivatives can be called if needed. The second argument is the index of the source, of the input list. The third argument is a list of dictionaries with some information about the image points. The "beta" entry contains the back-projected position, the "theta" entry is the observed position. The function should return a list of length one containing the final source position.

  • maxPermSize: when the more global approach is used to calculate the predicted image positions, these predicted positions need to be matched to the observed ones. 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.

  • localTraceFunction: when the local approach is used to calculate the predicted image positions starting from the observed ones, this specifies the specific method. Currently this can be:

    • “fsolve”: uses SciPi’s fsolve for the optimization, the localTraceFunctionOptions value should be None.

    • “pgraleretrace”: corresponds to “pgraleretrace,ExpandedMultiStepNewton”. These are CPU variants for the OpenCL code that’s used in the parametric inversion code. Use getDefaultsForRetraceType to see the default values for type names like "NoTrace" or "ExpandedMultiStepNewton". Such values can be changed using the localTraceFunctionOptions parameter.

    • “pgraleretrace,NoTrace”

    • “pgraleretrace,SingleStepNewton”

    • “pgraleretrace,MultiStepNewton”

    • “pgraleretrace,ExpandedMultiStepNewton”

    • a function to be called: in this case, the function is first called with arguments (allPoints, localTraceFunctionOptions). The first argument is a list of same length as imgList, and contains tuples of which the first entry is a list of dictionaries with "beta" and "theta" keys, for the back-projected and observed positions respectively. The second position in the tuple contains the redshift. The localTraceFunctionOptions argument will contain the value below. This function should return a dictionary, say newOptions.

      Next calls to this function, to actually calculate a predicted source position will have arguments (imgPlaneObject, beta, thetaStart, newOptions). The imgPlaneObject argument has members traceTheta and getBetaAndDerivatives, beta is the source position, thetaStart is the position to be adjusted so that this corresponds to beta when back-projected, and newOptions are the adjusted options from a first call to the function.

  • localTraceFunctionOptions: initial options that are passed to the previous function.

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 the "noreduction" setting was used 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.

grale.util.adjustShearMeasurements(pixelFrameCoords, pixelFrameGamma1, pixelFrameGamma2, centeredRaDecCoords, mirror=False, tol=None)

It is possible that shear measurements are specified in one coordinate frame, but you need to know them in a rotated, possibly mirrored frame. This function can perform the required calculation. In the names of the parameters we’re assuming that shear measurements need to be transformed from values that are given in a coordinate frame that’s aligned with the pixels of the CCD camera, and need to be transformed to a coordinate system that’s aligned with RA/Dec coordinates (but also recentered). Of course, these names just suggest coordinate systems, the transformation could be the other way around as well, for example.

Arguments:

  • pixelFrameCoords: The coordinates of the shear measurements in the coordinate frame that’s aligned with the CCD pixels.

  • pixelFrameGamma1: The first shear component as measured in the coordinate frame that’s aligned with the pixels.

  • pixelFrameGamma2: Same for the second shear component.

  • centeredRaDecCoords: The same points as specified in pixelFrameCoords, but now having their coordinates in the frame that’s being transformed to, for example one aligned with the RA/Dec axes.

  • mirror: In case the two coordinate systems don’t differ by simply a rotation, this flag can be set to indicate this.

  • tol: tolerance parameter that is passed to SciPy’s minimize function (“Nelder-Mead”) is used.

grale.util.readLenstoolBayesDat(fileName, columnsToRemove=['Nsample', 'ln(Lhood)', 'Chi2'])

Read the file with name fileName that contains the samples that a Lenstool inversion produced - typically this is called bayes.dat - and returns the result as a pandas DataFrame. Some columns are removed by default, but this can be changed by setting the columnsToRemove argument.

grale.util.readParametricSamples(parametersFileName, samplesFileName, renameColumnsFunction=<function _defaultRenameColumnsFunction>, renameColumnsFunctionParams={'angdist': (4.84813681109536e-06, 'arcsec'), 'angle': (1.0, 'degree'), 'density': (1.0, 'kgm2'), 'mass': (1.98855e+30, 'msun'), 'veldisp': (1000.0, 'kms')})

A parametric inversion with the eaType set to MCMC will produce a file containing the samples themselves, and a file describing the units these samples are represented in internally. This function reads these files, specified by parametersFileName and samplesFileName, and rescales the samples into common units like radians or kg. By default, the values will be rescaled further into more readable units, like arcsec or M_sun. The columns will be renamed as well to reflect this.

You can specify a custom function to do this second rescaling/renaming by setting the renameColumnsFunction argument. If None, this step will be disabled, otherwise it should be a function that accepts two arguments: the name of a column and the parameters set in renameColumnsFunctionParams. It should return either None if nothing should be changed about this column’s values, or a tuple (newName, scale) that will cause the values in the column to be divided by scale and will store these scaled values in a new column named newName.

grale.util.periodicWrapAround(x, a, b)

Function that does something modulus-like for real numbers, makes sure that x is adjusted to lie between a and b. Idea was that using goodman-weare no real periodic wrap-around could be used for angles (I think this would no longer satisty the detailed balance), and perhaps boundaries could be avoided entirely. This function could then be used to postprocess samples to lie within the boundaries.

Some testing indicated that having no boundaries could cause the angles to grow to very large values, causing the sampling to be not as good anymore.

grale.util.checkEllipseOverlap(A1, c1, A2, c2)

Checks if two N dimensional ellipsoids, where the A matrices describe the ellipse parameters and ‘c’ matrices define the center, overlap. If not, ‘None’ is returned, if so, a point in the overlapping region is returned.

The routine works by recentering the coordinates on the first ellipsoid, re-orienting the axes to align with the first ellipsiod’s axes, and then rescaling the axes so that the first ellipsoid is transformed into a unit sphere. An optimization routine then looks for the point inside the second ellipsoid that’s closest to the origin. If the distance to that point is less than unity, the ellipses overlap and the point is transformed to the original coordinates. The routine is repeated from the second ellipsoid’s point of view, and the average of the two points is what’s returned.

grale.util.getDataMuSigma(data, useQuartileEstimator=True)

Estimate mean and standard deviation for the data specified. If useQuartileEstimator is set to False, the np.mean and np.std functions are used, otherwise these values are estimated from the quartiles (using np.percentile).

grale.util.getDataRangeForHistogram(data, numSigma=5, useQuartileEstimator=True)

This is a helper function in creating a histogram: for the data specified, it determines a range that is to be used for the histogram. The idea is to determine a plot range that can exclude some outliers that may make the plot less clear.

If numSigma is a regular number, getDataMuSigma() is called with these data and the useQuartileEstimator parameter to get mean and standard deviation. The range is then numSigma times the standard deviation, on each side of the calculated mean

If numSigma is set to infinity, then minimum and maximum values of the data are used for the range.

The function returns the range, the mean, and the standard deviation.

grale.util.getAllDataRangesForHistogram(df, colNames, rangeNumSigma=5, useQuartileEstimator=True)

For each column name in colNames, the getDataRangeForHistogram() function is called on the corresponding column from the df DataFrame. The rangeNumSigma and userQuartileEstimator arguments are passed on to that function. A dictionary with the column names as keys and their corresponding ranges as values, is returned.

grale.util.mergeDataRanges(r1, r2)

Merge two data ranges r1 and r2. If one of them is None, the other range is returned. Otherwise a numpy array consisting of minimum and maximum values is returned.

grale.util.mergeAllDataRanges(r1s, r2s)

The inputs are expected to be dictionaries with the same keys, and data ranges for each entry. For each of the keys, mergeDataRanges() is called on the values, and a dictionary with the results is returned.