API¶
Warning
This is still in pre-alpha stage, not tested yet on full scale! Use at your own risk!
A collection of tools widely used in Weak Gravitational Lensing data analyses
Convergence maps¶
-
class
lenstools.
ConvergenceMap
(data, angle, masked=False, **kwargs)[source]¶ A class that handles 2D convergence maps and allows to compute their topological descriptors (power spectrum, peak counts, minkowski functionals)
>>> from lenstools import ConvergenceMap >>> from lenstools.defaults import load_fits_default_convergence
>>> test_map = ConvergenceMap.load("map.fit",format=load_fits_default_convergence) >>> imshow(test_map.data)
-
boundary
¶ Computes the boundaries of the masked regions, defined as the regions in which the convergence is still well defined but the first and second derivatives are not
Returns: array of bool. of the same shape as the map, with True values along the boundaries
-
cross
(other, statistic='power_spectrum', **kwargs)¶ Measures a cross statistic between maps
Parameters: - other (ConvergenceMap instance) – The other convergence map
- statistic (string or callable) – the cross statistic to measure (default is ‘power_spectrum’)
- kwargs (dict.) – the keyword arguments are passed to the statistic (when callable)
Returns: tuple – (l – array,Pl – array) = (multipole moments, cross power spectrum at multipole moments) if the statistic is the power spectrum, otherwise whatever statistic() returns on call
Raises: AssertionError if the other map has not the same shape as the input one
>>> test_map = ConvergenceMap.load("map.fit",format=load_fits_default_convergence) >>> other_map = ConvergenceMap.load("map2.fit",format=load_fits_default_convergence)
>>> l_edges = np.arange(200.0,5000.0,200.0) >>> l,Pl = test_map.cross(other_map,l_edges=l_edges)
-
cutRegion
(extent)¶ Cut a subset of the map, with the same resolution (warning! tested on square cuts only!)
Parameters: extent (array with units) – region in the map to select (xmin,xmax,ymin,ymax), must have units Returns: new ConvergenceMap instance that encloses the selected region
-
getValues
(x, y)¶ Extract the map values at the requested (x,y) positions; this is implemented using the numpy fast indexing routines, so the formats of x and y must follow the numpy advanced indexing rules. Periodic boundary conditions are enforced
Parameters: - x (numpy array or quantity) – x coordinates at which to extract the map values (if unitless these are interpreted as radians)
- y (numpy array or quantity) – y coordinates at which to extract the map values (if unitless these are interpreted as radians)
Returns: numpy array with the map values at the specified positions, with the same shape as x and y
Raises: IndexError if the formats of x and y are not the proper ones
-
gradient
(x=None, y=None, save=True)¶ Computes the gradient of the map and sets the gradient_x,gradient_y attributes accordingly
Parameters: - x (array with units) – optional, x positions at which to evaluate the gradient
- y (array with units) – optional, y positions at which to evaluate the gradient
- save (bool.) – if True saves the gradient as attrubutes
Returns: tuple – (gradient_x,gradient_y)
>>> test_map = ConvergenceMap.load("map.fit") >>> gx,gy = test_map.gradient()
-
hessian
(x=None, y=None, save=True)¶ Computes the hessian of the map and sets the hessian_xx,hessian_yy,hessian_xy attributes accordingly
Parameters: - x (array with units) – optional, x positions at which to evaluate the hessian
- y (array with units) – optional, y positions at which to evaluate the hessian
- save (bool.) – if True saves the gradient as attrubutes
Returns: tuple – (hessian_xx,hessian_yy,hessian_xy)
>>> test_map = ConvergenceMap.load("map.fit") >>> hxx,hyy,hxy = test_map.hessian()
-
info
¶ Displays some of the information stored in the map (mainly resolution)
-
load
(filename, format=None, **kwargs)¶ This class method allows to read the map from a data file, in various formats
Parameters: - filename (str.) – name of the file in which the map is saved
- format (str. or callable) – the format of the file in which the map is saved (can be a callable too); if None, it’s detected automatically from the filename
- kwargs (dict.) – the keyword arguments are passed to the format (if callable)
Returns: Spin0 instance with the loaded map
-
locatePeaks
(thresholds, norm=False)¶ Locate the peaks in the map
Parameters: - thresholds (array) – thresholds extremes that define the binning of the peak histogram
- norm (bool.) – normalization; if set to a True, interprets the thresholds array as units of sigma (the map standard deviation)
Returns: tuple – (peak height – array, peak locations – array with units)
Raises: AssertionError if thresholds array is not provided
>>> test_map = ConvergenceMap.load("map.fit") >>> thresholds = np.arange(map.data.min(),map.data.max(),0.05) >>> height,positions = test_map.locatePeaks(thresholds)
-
mask
(mask_profile, inplace=False)¶ Applies a mask to the convergence map: all masked pixels are given a nan value because they cannot be used in the calculations
Parameters: - mask_profile (array. or ConvergenceMap instance) – profile of the mask, must be an array of 1 byte intergers that are either 0 (if the pixel is masked) or 1 (if the pixel is not masked). Must be of the same shape as the original map
- inplace – if True the masking is performed in place and the original map is lost, otherwise a new instance of ConvergenceMap with the masked map is returned
Returns: the masked convergence map if inplace is False, otherwise a float corresponding to the masked fraction of the map
-
maskBoundaries
()¶ Computes the mask boundaries defined in the following way: a boundary is a region where the convergence value is defined, but the gradients are not defined.
Returns: float. : perimeter/area ratio of the mask boundaries
-
maskedFraction
¶ Computes the masked fraction of the map
Returns: float, the masked fraction of the map
-
mean
()¶ Computes the mean value of the pixels, taking into account eventual masking
Returns: float.
-
minkowskiFunctionals
(thresholds, norm=False)¶ Measures the three Minkowski functionals (area,perimeter and genus characteristic) of the specified map excursion sets
Parameters: - thresholds (array) – thresholds that define the excursion sets to consider
- norm (bool.) – normalization; if set to a True, interprets the thresholds array as units of sigma (the map standard deviation)
Returns: tuple – (nu – array, V0 – array, V1 – array, V2 – array) nu are the bins midpoints and V are the Minkowski functionals
Raises: AssertionError if thresholds array is not provided
>>> test_map = ConvergenceMap.load("map.fit") >>> thresholds = np.arange(-2.0,2.0,0.2) >>> nu,V0,V1,V2 = test_map.minkowskiFunctionals(thresholds,norm=True)
-
moments
(connected=False, dimensionless=False)¶ Measures the first nine moments of the convergence map (two quadratic, three cubic and four quartic)
Parameters: - connected (bool.) – if set to True returns only the connected part of the moments
- dimensionless (bool.) – if set to True returns the dimensionless moments, normalized by the appropriate powers of the variance
Returns: array – (sigma0,sigma1,S0,S1,S2,K0,K1,K2,K3)
>>> test_map = ConvergenceMap.load("map.fit") >>> var0,var1,sk0,sk1,sk2,kur0,kur1,kur2,kur3 = test_map.moments() >>> sk0,sk1,sk2 = test_map.moments(dimensionless=True)[2:5] >>> kur0,kur1,kur2,kur3 = test_map.moments(connected=True,dimensionless=True)[5:]
-
pdf
(thresholds, norm=False)¶ Computes the one point probability distribution function of the convergence map
Parameters: - thresholds (array) – thresholds extremes that define the binning of the pdf
- norm (bool.) – normalization; if set to a True, interprets the thresholds array as units of sigma (the map standard deviation)
Returns: tuple – (threshold midpoints – array, pdf normalized at the midpoints – array)
Raises: AssertionError if thresholds array is not provided
>>> test_map = ConvergenceMap.load("map.fit") >>> thresholds = np.arange(map.data.min(),map.data.max(),0.05) >>> nu,p = test_map.pdf(thresholds)
-
peakCount
(thresholds, norm=False)¶ Counts the peaks in the map
Parameters: - thresholds (array) – thresholds extremes that define the binning of the peak histogram
- norm (bool.) – normalization; if set to a True, interprets the thresholds array as units of sigma (the map standard deviation)
Returns: tuple – (threshold midpoints – array, differential peak counts at the midpoints – array)
Raises: AssertionError if thresholds array is not provided
>>> test_map = ConvergenceMap.load("map.fit") >>> thresholds = np.arange(map.data.min(),map.data.max(),0.05) >>> nu,peaks = test_map.peakCount(thresholds)
-
powerSpectrum
(l_edges)¶ Measures the power spectrum of the convergence map at the multipole moments specified in the input
Parameters: l_edges (array) – Multipole bin edges Returns: tuple – (l – array,Pl – array) = (multipole moments, power spectrum at multipole moments) Raises: AssertionError if l_edges are not provided >>> test_map = ConvergenceMap.load("map.fit") >>> l_edges = np.arange(200.0,5000.0,200.0) >>> l,Pl = test_map.powerSpectrum(l_edges)
-
save
(filename, format=None, double_precision=False)¶ Saves the map to an external file, of which the format can be specified (only fits implemented so far)
Parameters: - filename (str.) – name of the file on which to save the plane
- format (str.) – format of the file, only FITS implemented so far; if None, it’s detected automatically from the filename
- double_precision (bool.) – if True saves the Plane in double precision
-
savefig
(filename)¶ Saves the map visualization to an external file
Parameters: filename (str.) – name of the file on which to save the map
-
setAngularUnits
(unit)¶ Convert the angular units of the map to the desired unit
Parameters: unit (astropy units) – astropy unit instance to which to perform the conversion
-
smooth
(scale_angle, kind='gaussian', inplace=False, **kwargs)¶ Performs a smoothing operation on the convergence map
Parameters: - scale_angle (float.) – size of the smoothing kernel (must have units)
- kind (str.) – type of smoothing to be performed (only implemented gaussian so far)
- inplace (bool.) – if set to True performs the smoothing in place overwriting the old convergence map
- kwargs (dict.) – the keyword arguments are passed to the filter function
Returns: ConvergenceMap instance (or None if inplace is True)
-
std
()¶ Computes the standard deviation of the pixels, taking into account eventual masking
Returns: float.
-
visualize
(fig=None, ax=None, colorbar=False, **kwargs)¶ Visualize the convergence map; the kwargs are passed to imshow
-
Shear maps¶
-
class
lenstools.
ShearMap
(data, angle)[source]¶ A class that handles 2D shear maps and allows to perform a set of operations on them
>>> from lenstools import ShearMap
>>> test = ShearMap.load("shear.fit",format=lenstools.defaults.load_fits_default_shear) >>> test.side_angle 1.95 >>> test.data #The actual map values
-
convergence
()[source]¶ Reconstructs the convergence from the E component of the shear
Returns: new ConvergenceMap instance
-
decompose
(l_edges, keep_fourier=False)¶ Decomposes the shear map into its E and B modes components and returns the respective power spectral densities at the specified multipole moments
Parameters: - l_edges (array) – Multipole bin edges
- keep_fourier (bool.) – If set to True, holds the Fourier transforms of the E and B mode maps into the E and B attributes of the ShearMap instance
Returns: returns: tuple – (l – array,P_EE,P_BB,P_EB – arrays) = (multipole moments, EE,BB power spectra and EB cross power) >>> test_map = ShearMap.load("shear.fit",format=load_fits_default_shear) >>> l_edges = np.arange(300.0,5000.0,200.0) >>> l,EE,BB,EB = test_map.decompose(l_edges)
-
fromEBmodes
(fourier_E, fourier_B, angle=<Quantity 3.14 deg>)¶ This class method allows to build a shear map specifying its E and B mode components
Parameters: - fourier_E (numpy 2D array, must be of type np.complex128 and must have a shape that is appropriate for a real fourier transform, i.e. (N,N/2 + 1); N should be a power of 2) – E mode of the shear map in fourier space
- fourier_B (numpy 2D array, must be of type np.complex128 and must have a shape that is appropriate for a real fourier transform, i.e. (N,N/2 + 1); N should be a power of 2) – B mode of the shear map in fourier space
- angle (float.) – Side angle of the real space map in degrees
Returns: the corresponding ShearMap instance
Raises: AssertionErrors for inappropriate inputs
-
getValues
(x, y)¶ Extract the map values at the requested (x,y) positions; this is implemented using the numpy fast indexing routines, so the formats of x and y must follow the numpy advanced indexing rules. Periodic boundary conditions are enforced
Parameters: - x (numpy array or quantity) – x coordinates at which to extract the map values (if unitless these are interpreted as radians)
- y (numpy array or quantity) – y coordinates at which to extract the map values (if unitless these are interpreted as radians)
Returns: numpy array with the map values at the specified positions, with shape (N,shape x) where N is the number of components of the map field
Raises: IndexError if the formats of x and y are not the proper ones
-
gradient
(x=None, y=None)¶ Computes the gradient of the components of the spin1 field at each point
Parameters: - x (array with units) – optional, x positions at which to evaluate the gradient
- y (array with units) – optional, y positions at which to evaluate the gradient
Returns: the gradient of the spin1 field in array form, of shape (4,:,:) where the four components are, respectively, 1x,1y,2x,2y; the units for the finite difference are pixels
-
info
¶ Displays some of the information stored in the map (mainly resolution)
-
load
(filename, format=None, **kwargs)¶ This class method allows to read the map from a data file, in various formats
Parameters: - filename (str.) – name of the file in which the map is saved
- format (str. or callable) – the format of the file in which the map is saved (can be a callable too); if None, it’s detected automatically from the filename
- kwargs (dict.) – the keyword arguments are passed to the format (if callable)
Returns: Spin1 instance with the loaded map
-
save
(filename, format=None, double_precision=False)¶ Saves the map to an external file, of which the format can be specified (only fits implemented so far)
Parameters: - filename (str.) – name of the file on which to save the plane
- format (str.) – format of the file, only FITS implemented so far; if None, it’s detected automatically from the filename
- double_precision (bool.) – if True saves the Plane in double precision
-
savefig
(filename)¶ Saves the map visualization to an external file
Parameters: filename (str.) – name of the file on which to save the map
-
setAngularUnits
(unit)¶ Convert the angular units of the map to the desired unit
Parameters: unit (astropy units) – astropy unit instance to which to perform the conversion
-
sticks
(fig=None, ax=None, pixel_step=10, multiplier=1.0)¶ Draw the ellipticity map using the shear components
Parameters: - ax (matplotlib ax object) – ax on which to draw the ellipticity field
- pixel_step (int.) – One arrow will be drawn every pixel_step pixels to avoid arrow overplotting
- multiplier (float.) – Multiplies the stick length by a factor
Returns: ax – the matplotlib ax object on which the stick field was drawn
>>> import matplotlib.pyplot as plt >>> test = ShearMap.load("shear.fit",loader=load_fits_default_shear) >>> fig,ax = plt.subplots() >>> test.sticks(ax,pixel_step=50)
-
visualize
(fig=None, ax=None, component_labels=('$\\gamma_1$', '$\\gamma_2$'), colorbar=False, **kwargs)¶ Visualize the shear map; the kwargs are passed to imshow
-
visualizeComponents
(fig, ax, components='EE, BB, EB', region=(200, 9000, -9000, 9000))¶ Plots the full 2D E and B mode power spectrum (useful to test statistical isotropicity)
Parameters: - fig (matplotlibfigure object) – figure on which to draw the ellipticity field
- ax (matplotlib ax object or array of ax objects, can be None in which case new ax are created) – ax on which to draw the ellipticity field
- components – string that contains the components to plot; the format is a sequence of {EE,BB,EB} separated by commas
- region (tuple (lx_min,lx_max,ly_min,ly_max)) – selects the multipole region to visualize
-
Statistics¶
-
class
lenstools.
Ensemble
(file_list=[], data=None, num_realizations=0, metric='chi2')[source]¶ A class that handles statistical operations on weak lensing maps; an ensemble is a collection of different statistical realization of the same random variable. This class has an attribute ‘data’ that is a numpy array which first axis corresponds to the realization number.
>>> from lenstools.statistics import Ensemble
-
bootstrap
(callback, bootstrap_size=10, resample=10, seed=None)[source]¶ Computes a custom statistic on the Ensemble using the bootstrap method
Parameters: - callback (callable) – statistic to compute on the ensemble; takes the resampled Ensemble data as an input
- bootstrap_size (int.) – size of the resampled ensembles used in the bootstraping; must be less than or equal to the number of realizations in the Ensemble
- resample (int.) – number of times the Ensemble is resampled
- seed (int.) – if not None, this is the random seed of the random resamples
Returns: the bootstraped Ensemble statistic
-
compare
(rhs, **kwargs)[source]¶ Computes a comparison score between two Ensembles: computes a chi2-style difference between two different ensembles to assert how different they are
-
covariance
(bootstrap=False, **kwargs)[source]¶ Computes the ensemble covariance matrix
Parameters: - bootstrap (bool.) – if True the covariance matrix is computed with a bootstrap estimate
- kwargs (dict.) – the keyword arguments are passed to the bootstrap method
Returns: ndarray with the covariance matrix, has shape (self.data[0],self.data[0])
-
cut
(min=None, max=None, feature_label=None, inplace=True)[source]¶ Allows to manually cut the ensemble along the second axis, if you want to select a feature subset; you better know what you are doing if you are using this function
Parameters: - min (int or float) – left extreme of the cut, included; if a list of indices is passed, the cut is performed on those indices, on the second axis and the remaining parameters are ignored
- max (int or float) – right extreme of the cut, included
- feature_label (array) – if not None, serves as a reference for min and max, in which the ensemble is cut according to the position of the min and max elements in the feature_label array
- inplace (bool.) – if True, the Ensemble cut is made in place, otherwise a new ensemble is returned
Returns: the min and max cut indices if feature_label is None, otherwise it returns the array of the cut feature; if min was a list of indices, these are returned back. Returns the new Ensemble if inplace is False
-
differentiate
(step=None, order=1)[source]¶ Compute a new ensemble, in which the feature is differentiated with respect to its label (e.g. the differentiation of the first minkowski functional is the PDF)
Parameters: step (float or array) – measure unit on the feature label axis Returns: new Ensemble instance with the differentiated feature
-
classmethod
fromdata
(npy_data)[source]¶ Builds the ensemble from data in numpy array format
Parameters: npy_data (array) – numpy array with the data, the first dimension must be the number of realizations
-
classmethod
fromfilelist
(file_list)[source]¶ Builds the ensemble from a file list: each file corresponds to a different realization
Parameters: file_list (list of str.) – List of files on which to define the ensemble
-
group
(group_size, kind='sparse', inplace=True)[source]¶ Sometimes it happens that different realizations in the ensemble need to be grouped together, for example when they belong to different subfields of the same observation field. With this function you can group different realizations together by taking the mean, and reduce the total number of realizations in the ensemble
Parameters: - group_size (int) – how many realizations to put in a group, must divide exactly the total number of realizations in the ensemble
- kind (str. or scipy.sparse) – specifies how to do the grouping; if set to “sparse” the groups are formed by taking one realizations every num_realizations/group_size (for example ([1,1001,...,9001],[2,1002,...,9002]) if num_realizations=10000 and group_size=10). If set to “contiguous” then the realizations are grouped as ([1,2,...,10],[11,12,...,20]). Otherwise you can set kind to your own sparse matrix scheme
-
load
(callback_loader=None, pool=None, from_old=False, **kwargs)[source]¶ Loads the ensemble into memory, can spread the calculations on multiple processors using a MPI pool
Parameters: - callback_loader (function) – This function gets executed on each of the files in the list and populates the ensemble. If None provided, it performs a numpy.load on the specified files. Must return a numpy array with the loaded data
- pool (MPI pool object) – MPI pool for multiprocessing (imported from emcee https://github.com/dfm/emcee)
- from_old (bool.) – If True, the loaded data are interpreted as an old, already existing ensemble, which means that only one file (in which the old ensemble is saved) is loaded, the first dimension of the data is 1 and hence it is discarded
- kwargs (dict.) – Any additional keyword arguments to be passed to callback_loader
>>> from lenstools import Ensemble >>> from lenstools.statistics import default_callback_loader
>>> map_list = ["conv1.fit","conv2.fit","conv3.fit"] >>> l_edges = np.arange(200.0,50000.0,200.0)
>>> conv_ensemble = Ensemble.fromfilelist(map_list) >>> conv_ensemble.load(callback_loader=default_callback_loader,pool=pool,l_edges=l_edges)
-
mean
()[source]¶ Computes the ensemble average over realizations
Returns: ndarray with the averages, has the same shape as self.data[0]
-
principalComponents
()[source]¶ Computes the principal components of the Ensemble
Returns: pcaHandler instance
-
classmethod
read
(filename, callback_loader=None, **kwargs)[source]¶ Reads a numpy file into an Ensemble
Parameters: - filename (str.) – name of the file to read
- callback_loader (function) – This function gets executed on each of the files in the list and populates the ensemble. If None provided, it performs a numpy.load on the specified file. Must return a numpy array with the loaded data
- kwargs (dict.) – Any additional keyword arguments to be passed to callback_loader
Returns: Ensemble instance read from the file
-
classmethod
readall
(filelist, callback_loader=None, **kwargs)[source]¶ Reads a list of files into an Ensemble
Parameters: - filelist (list.) – list of files to read
- callback_loader (function) – This function gets executed on each of the files in the list and populates the ensemble. If None provided, it performs a numpy.load on the specified file. Must return a numpy array with the loaded data
- kwargs (dict.) – Any additional keyword arguments to be passed to callback_loader
Returns: Ensemble instance read from the file
-
save
(filename, format=None, **kwargs)[source]¶ Save ensemble data in an external file (in arbitrary format)
Parameters: - filename (str.) – file name of the external file
- kwargs (dict.) – the keyword arguments are passed to the saver (or to format if callable)
Format: format in which to save the ensemble; if None the format is auto detected from the filename
-
scale
(weights)[source]¶ Set manually the units of each row of the ensemble by multiplying it by a row of weights
Parameters: weights (array) – row of weights used to scale the ensemble, must have the same length as the number of rows
-
selfChi2
()[source]¶ Computes the Ensemble distribution of chi squared values, defined with respect to the same Ensemble mean and covariance
Returns: array with the self chi squared for each realization
-
shuffle
(seed=None)[source]¶ Changes the order of the realizations in the Ensemble
Parameters: seed (int.) – random seed for the random shuffling
-
split
(index)[source]¶ Inverse of the * operator: this method uses an Indexer instance to break down a multiple descriptor ensemble in many, smaller, single descriptor ensembles
Parameters: index (Indexer instance) – index of descriptors with which to perform the split Returns: list of Ensemble instances, one for each element in index Raises: AssertionError if shape of the ensemble data is not suitable
-
subset
(realizations)[source]¶ Returns a new ensemble that contains only the selected realizations
Parameters: realizations (int. or array of int.) – realizations to keep in the subset Ensemble, must be compatible with numpy array indexing syntax Returns: new Ensemble instance that contains only the selected realizations
-
transform
(transformation, inplace=False, **kwargs)[source]¶ Allows a general transformation on the Ensemble by calling a function on its data
Parameters: - transformation (callable) – callback function called on the ensemble data
- inplace (bool.) – if True the transformation is performed in place, otherwise a new Ensemble is created
- kwargs (dict.) – the keyword arguments are passed to the transformation callable
Returns: new Ensemble instance if the transformation is nor performed in place
-
Constraining cosmology¶
-
class
lenstools.constraints.
Analysis
(parameter_set=None, training_set=None)[source]¶ The base class of this module; the idea in weak lensing analysis is that one has a set of simulated data, that serves as training model, and then uses that set to fit the observations for the best model parameters.
Parameters: - parameter_set (array) – the points in parameter space that the simulated data cover; the first axis refers to the model, the second to the number of parameters
- training_set (array) – the measured feature in each model; the first axis refers to the model, the others to the feature indices (or bins)
- observed_set (array) – the measured feature in the data, should be a one dimensional array
-
add_feature_label
(feature_label)[source]¶ Add a feature label to the current analysis, i.e. a set of multipole moments if the feature is the power spectrum, or a set of thresholds if the feature is a PDF, etc...
Parameters: feature_label (array.) – the feature label to add, must have the same shape as the training set
-
add_model
(parameters, feature)[source]¶ Add a model to the training set of the current analysis
Parameters: - parameters (array) – parameter set of the new model
- feature (array) – measured feature of the new model
-
find
(parameters, rtol=1e-05)[source]¶ Finds the index of the training model that has the specified combination of parameters
Parameters: - parameters (array.) – the parameters of the model to find
- rtol (float.) – tolerance of the search (must be less than 1)
Returns: array of int. with the indices of the corresponding models
-
classmethod
load
(filename)[source]¶ Unpickle a previously pickled instance: be sure the file you are unpickling comes from a trusted source, this operation is potentially dangerous for your computer!
Parameters: filename (str. or file object) – Name of the file from which you want to load the instance, or file object
-
principalComponents
()[source]¶ Computes the principal components of the training_set
Returns: pcaHandler instance
-
remove_model
(model_list)[source]¶ Remove one or more models from the Analysis instance
Parameters: model_list (int. or list of int.) – list of the indices of the models to remove (0 indicates the first model)
-
reparametrize
(formatter, *args, **kwargs)[source]¶ Reparametrize the parameter set of the analysis by calling the formatter handle on the current parameter set (can be used to enlarge/shrink/relabel the parameter set)
Parameters: formatter (callable) – formatter function called on the current parameter_set (args and kwargs are passed to it)
-
save
(filename)[source]¶ Save the current Analysis instance as a pickled binary file for future reuse as an emulator; useful after you trained the emulator with a simulated feature set and you want to reuse it in the future
Parameters: filename (str. or file object) – Name of the file to which you want to save the emulator, or file object
-
transform
(transformation, inplace=False, **kwargs)[source]¶ Allows a general transformation on the training_set of the analysis by calling an arbitrary transformation function
Parameters: - transformation (callable) – callback function called on the training_set
- inplace (bool.) – if True the transformation is performed in place, otherwise a new Analysis instance is created
- kwargs (dict.) – the keyword arguments are passed to the transformation callable
-
class
lenstools.constraints.
FisherAnalysis
(parameter_set=None, training_set=None)[source]¶ -
add_feature_label
(feature_label)¶ Add a feature label to the current analysis, i.e. a set of multipole moments if the feature is the power spectrum, or a set of thresholds if the feature is a PDF, etc...
Parameters: feature_label (array.) – the feature label to add, must have the same shape as the training set
-
check
()[source]¶ Asserts that the parameters are varied one at a time, and that a parameter is not varied more than once
Raises: AssertionError
-
chi2
(observed_feature, features_covariance)[source]¶ Computes the chi2 between an observed feature and the fiducial feature, using the provided covariance
Parameters: - observed_feature (array) – observed feature to fit, its last dimension must have the same shape as self.training_set[0]
- features_covariance (2 dimensional array (or 1 dimensional if diagonal)) – covariance matrix of the simulated features, must be provided for a correct fit!
Returns: chi2 of the comparison
Return type: float.
-
classify
(observed_feature, features_covariance, labels=[0, 1], confusion=False)[source]¶ Performs a Fisher classification of the observed feature, choosing the most probable label based on the value of the chi2
Parameters: - observed_feature (array) – observed feature to fit, the last dimenstion must have the same shape as self.training_set[0]
- features_covariance (2 dimensional array (or 1 dimensional if assumed diagonal)) – covariance matrix of the simulated features, must be provided for a correct classification!
- labels (iterable) – labels of the classification, must be the indices of the available classes (from 0 to training_set.shape[0])
- confusion (bool.) – if True, an array with the label percentage occurrences is returned; if False an array of labels is returned
Returns: array with the labels resulting from the classification
Return type: int.
-
compute_derivatives
()[source]¶ Computes the feature derivatives with respect to the parameter sets using one step finite differences; the derivatives are computed with respect to the fiducial parameter set
Returns: array of shape (p,N), where N is the feature dimension and p is the number of varied parameters
-
find
(parameters, rtol=1e-05)¶ Finds the index of the training model that has the specified combination of parameters
Parameters: - parameters (array.) – the parameters of the model to find
- rtol (float.) – tolerance of the search (must be less than 1)
Returns: array of int. with the indices of the corresponding models
-
fisher_matrix
(simulated_features_covariance, observed_features_covariance=None)[source]¶ Computes the Fisher matrix of the associated features, that in the end allows to compute the paramter confidence contours (around the fiducial value)
Parameters: - simulated_features_covariance (2 dimensional array (or 1 dimensional if assumed diagonal)) – covariance matrix of the simulated features, must be provided for a correct fit!
- observed_features_covariance (2 dimensional array (or 1 dimensional if assumed diagonal)) – covariance matrix of the simulated features, if different from the simulated one; if None the simulated feature covariance is used
Returns: 2 dimensional array with the Fisher matrix of the analysis
-
fit
(observed_feature, features_covariance)[source]¶ Maximizes the gaussian likelihood on which the Fisher matrix formalism is based, and returns the best fit for the parameters given the observed feature
Parameters: - observed_feature (array) – observed feature to fit, must have the same shape as self.training_set[0]
- features_covariance (2 dimensional array (or 1 dimensional if assumed diagonal)) – covariance matrix of the simulated features, must be provided for a correct fit!
Returns: array with the best fitted parameter values
-
load
(filename)¶ Unpickle a previously pickled instance: be sure the file you are unpickling comes from a trusted source, this operation is potentially dangerous for your computer!
Parameters: filename (str. or file object) – Name of the file from which you want to load the instance, or file object
-
observables2parameters
(features_covariance=None)[source]¶ Computes the conversion matrix M that allows to match a feature vector V to its best fit parameters P, in the sense P = P[fiducial] + MV
Parameters: features_covariance (2 dimensional array (or 1 dimensional if diagonal)) – covariance matrix of the simulated features, must be provided! Returns: the (p,N) conversion matrix Return type: array
-
principalComponents
()¶ Computes the principal components of the training_set
Returns: pcaHandler instance
-
remove_model
(model_list)¶ Remove one or more models from the Analysis instance
Parameters: model_list (int. or list of int.) – list of the indices of the models to remove (0 indicates the first model)
-
save
(filename)¶ Save the current Analysis instance as a pickled binary file for future reuse as an emulator; useful after you trained the emulator with a simulated feature set and you want to reuse it in the future
Parameters: filename (str. or file object) – Name of the file to which you want to save the emulator, or file object
-
set_fiducial
(n)[source]¶ Sets the fiducial model (with respect to which to compute the derivatives), default is 0 (i.e. self.parameter_set[0])
Parameters: n (int.) – the parameter set you want to use as fiducial
-
transform
(transformation, inplace=False, **kwargs)¶ Allows a general transformation on the training_set of the analysis by calling an arbitrary transformation function
Parameters: - transformation (callable) – callback function called on the training_set
- inplace (bool.) – if True the transformation is performed in place, otherwise a new Analysis instance is created
- kwargs (dict.) – the keyword arguments are passed to the transformation callable
-
variations
¶ Checks the parameter variations with respect to the fiducial cosmology
Returns: iterable with the positions of the variations
-
varied
¶ Returns the indices of the parameters that are varied
Returns: list with the indices of the varied parameters
-
-
class
lenstools.constraints.
LikelihoodAnalysis
(parameter_set=None, training_set=None)[source]¶ -
add_feature_label
(feature_label)¶ Add a feature label to the current analysis, i.e. a set of multipole moments if the feature is the power spectrum, or a set of thresholds if the feature is a PDF, etc...
Parameters: feature_label (array.) – the feature label to add, must have the same shape as the training set
-
add_model
(parameters, feature)¶ Add a model to the training set of the current analysis
Parameters: - parameters (array) – parameter set of the new model
- feature (array) – measured feature of the new model
-
chi2
(parameters, observed_feature, features_covariance, split_chunks=None, pool=None)[source]¶ Computes the chi2 part of the parameter likelihood with the usual sandwich product with the covariance matrix; the model features are computed with the interpolators
Parameters: - parameters ((N,p) array where N is the number of points and p the number of parameters) – new points in parameter space on which to compute the chi2 statistic
- observed_feature (array) – observed feature on which to condition the parameter likelihood
- features_covariance (array) – covariance matrix of the features, must be supplied
- split_chunks (int.) – if set to an integer bigger than 0, splits the calculation of the chi2 into subsequent chunks, each that takes care of an equal number of points. Each chunk could be taken care of by a different processor
Returns: array with the chi2 values, with the same shape of the parameters input
-
chi2Contributions
(parameters, observed_feature, features_covariance)[source]¶ Computes the individual contributions of each feature bin to the chi2; the model features are computed with the interpolators. The full chi2 is the sum of the individual contributions
Parameters: - parameters ((N,p) array where N is the number of points and p the number of parameters) – new points in parameter space on which to compute the chi2 statistic
- observed_feature (array) – observed feature on which to condition the parameter likelihood
- features_covariance (array) – covariance matrix of the features, must be supplied
Returns: numpy 2D array with the contributions to the chi2 (off diagonal elements are the contributions of the cross correlation between bins)
-
find
(parameters, rtol=1e-05)¶ Finds the index of the training model that has the specified combination of parameters
Parameters: - parameters (array.) – the parameters of the model to find
- rtol (float.) – tolerance of the search (must be less than 1)
Returns: array of int. with the indices of the corresponding models
-
likelihood
(chi2_value, **kwargs)[source]¶ Computes the likelihood value with the selected likelihood function, given the pre-computed chi2 value
Parameters: - chi2_value (array) – chi squared values
- kwargs – keyword arguments to be passed to your likelihood function
-
load
(filename)¶ Unpickle a previously pickled instance: be sure the file you are unpickling comes from a trusted source, this operation is potentially dangerous for your computer!
Parameters: filename (str. or file object) – Name of the file from which you want to load the instance, or file object
-
predict
(parameters)[source]¶ Predicts the feature at a new point in parameter space using the bin interpolators, trained with the simulated features
Parameters: parameters (array) – new points in parameter space on which to compute the chi2 statistic; it’a (N,p) array where N is the number of points and p the number of parameters, or array of size p if there is only one point
-
principalComponents
()¶ Computes the principal components of the training_set
Returns: pcaHandler instance
-
remove_model
(model_list)¶ Remove one or more models from the Analysis instance
Parameters: model_list (int. or list of int.) – list of the indices of the models to remove (0 indicates the first model)
-
save
(filename)¶ Save the current Analysis instance as a pickled binary file for future reuse as an emulator; useful after you trained the emulator with a simulated feature set and you want to reuse it in the future
Parameters: filename (str. or file object) – Name of the file to which you want to save the emulator, or file object
-
set_likelihood
(function=None)[source]¶ Sets the likelihood function to a custom function input by the user: the default is the usual exp(-0.5*chi^2)
-
train
(use_parameters='all', **kwargs)[source]¶ Builds the interpolators for each of the feature bins using a radial basis function approach
Parameters: - use_parameters (list. or “all”) – which parameters actually vary in the supplied parameter set (it doesn’t make sense to interpolate over the constant ones)
- kwargs – keyword arguments to be passed to the interpolator constructor
-
Confidence contour plotting¶
-
class
lenstools.contours.
ContourPlot
(fig=None, ax=None)[source]¶ A class handler for contour plots
-
expectationValue
(function, **kwargs)[source]¶ Computes the expectation value of a function of the parameters over the current parameter likelihood
-
getLikelihood
(likelihood_filename, parameter_axes={'w': 1, 'Omega_m': 0, 'sigma8': 2}, parameter_labels={'w': '$w$', 'Omega_m': '$\\Omega_m$', 'sigma8': '$\\sigma_8$'})[source]¶ Load the likelihood function from a numpy file
-
getLikelihoodValues
(levels, precision=0.001)[source]¶ Find the likelihood values that correspond to the selected p_values
-
getMaximum
(which='full')[source]¶ Find the point in parameter space on which the likelihood is maximum
-
marginal
(parameter_name='w', levels=None)[source]¶ Marginalize the likelihood over all parameters but one
-
plotContours
(colors=['red', 'green', 'blue'], display_percentages=True, display_maximum=True, fill=False, **kwargs)[source]¶ Display the confidence likelihood contours
-
plotMarginal
(parameter, levels=[0.684], colors=['red', 'blue', 'green'], alpha=0.5, fill=False)[source]¶ Plot the likelihood function marginalized over all parameters except one
-
point
(coordinate_x, coordinate_y, color='green', marker='o')[source]¶ Draws a point in parameter space at the specified physical coordinates
-
setUnits
(parameter, parameter_min, parameter_max, parameter_unit)[source]¶ Set manually the physical units for each of the likelihood axes
-
slice
(parameter_name='w', parameter_value=-1.0)[source]¶ Slice the likelihood cube by fixing one of the parameters
-
Noise¶
-
class
lenstools.
GaussianNoiseGenerator
(shape, side_angle, label)[source]¶ A class that handles generation of Gaussian simulated noise maps
-
classmethod
forMap
(conv_map)[source]¶ This class method generates a Gaussian noise generator intended to be used on a convergence map: i.e. the outputs of its methods can be added to the convergence map in question to simulate the presence of noise
Parameters: conv_map (ConvergenceMap instance) – The blueprint of the convergence map you want to generate the noise for Raises: AssertionError if conv_map is not a ConvergenceMap instance
-
fromConvPower
(power_func, seed=0, **kwargs)[source]¶ This method uses a supplied power spectrum to generate correlated noise maps in real space via FFTs
Parameters: - power_func (function with the above specifications, or numpy array (l,Pl) of shape (2,n)) – function that given a numpy array of l’s returns a numpy array with the according Pl’s (this is the input power spectrum); alternatively you can pass an array (l,Pl) and the power spectrum will be calculated with scipy’s interpolation routines
- seed (int.) – seed of the random generator
- kwargs – keyword arguments to be passed to power_func, or to the interpolate.interp1d routine
Returns: ConvergenceMap instance of the same exact shape as the one used as blueprint
-
getShapeNoise
(z=1.0, ngal=<Quantity 15.0 1 / arcmin2>, seed=0)[source]¶ This method generates a white, gaussian shape noise map for the given redshift of the map
Parameters: - z (float.) – single redshift of the backround sources on the map
- ngal (float.) – assumed angular number density of galaxies (must have units of angle^-2)
- seed (int.) – seed of the random generator
Returns: ConvergenceMap instance of the same exact shape as the one used as blueprint
-
classmethod
Indexing¶
-
class
lenstools.index.
Indexer
(descriptor_list)[source]¶ This class is useful for indexing an array of statistical descriptors that are hstacked together; the Indexer instance keeps track of the memory regions where the different descriptors are stored
Parameters: descriptor_list (list.) – list of Descriptor subinstances, such as PowerSpectrum; each of these sub instances must be a subclass of Descriptor with a ‘first’ and ‘last’ getter methods implemented
-
class
lenstools.index.
PowerSpectrum
(l_edges)[source]¶ Power spectrum indexing class
Parameters: l_edges (array) – bin edges of multipole moments
-
class
lenstools.index.
Peaks
(thresholds, norm=False)[source]¶ Peaks histogram indexing class
Parameters: thresholds (array) – peak histogram bin edges
-
class
lenstools.index.
PDF
(thresholds, norm=False)[source]¶ Probability density function indexing class, identical to Peaks
Simulation suites¶
-
class
lenstools.simulations.
IGS1
(H0=72.0, Om0=0.26, w0=-1.0, sigma8=0.798, ns=0.96, root_path=None, name=None)[source]¶ Class handler of the IGS1 simulations set, inherits the cosmological parameters from the astropy.cosmology.FlatwCDM class; the default parameter values are the fiducial ones
-
getNames
(realizations, z=1.0, kind='convergence', big_fiducial_set=False)[source]¶ Get the full name of the IGS1 maps, once a redshift, realization identificator and kind are specified
Parameters: - z (float.) – redshift plane of the maps, must be one of [1.0,1.5,2.0]
- realizations (list. or int.) – list of realizations to get the names of, the elements must be in [1,1000]
- kind (str.) – decide if retrieve convergence or shear maps, must be one of [convergence,shear1,shear2]
- big_fiducial_set (bool.) – set to True if you want to get the names of the bigger fiducial simulation based on 45 N-body simulations
-
load
(realization, **kwargs)[source]¶ Reads in a specific realization of the convergence field (in FITS format) and returns a ConvergenceMap instance with the loaded map
Parameters: - realization (int.) – the specific realization to read
- kwargs – the keyword arguments are passed to the getNames method
Returns: ConvergenceMap instance with the loaded map
-
-
class
lenstools.simulations.
CFHTemu1
(H0=72.0, Om0=0.26, w0=-1.0, sigma8=0.798, ns=0.96, root_path=None, name=None)[source]¶ Class handler of the weak lensing CFHTemu1 simulations set, inherits from IGS1; this simulation suite contains 91 different cosmological models based on 1 N-body simulation each. Each model has 1000 realizations for each of the 13 CFHT subfields
-
classmethod
getModels
(root_path='/default')[source]¶ This class method uses a dictionary file to read in all the cosmological model parameters and instantiate the corresponding CFHTemu1 objects for each one of them
Parameters: root_path (str.) – path of your CFHT emu1 simulations copy Returns: list of CFHTemu1 instances
-
getNames
(realizations, subfield=1, smoothing=0.5)[source]¶ Get the full name of the CFHT emu1 maps, once a subfield and smoothing scale are specified
Parameters: - subfield (int.) – the specific CFHT subfield you want to retrieve, must be between 1 and 13
- smoothing (float.) – smoothing scale of the maps you wish to retrieve, in arcmin
- realizations (list. or int.) – list of realizations to get the names of, the elements must be in [1,1000]
-
load
(realization, **kwargs)¶ Reads in a specific realization of the convergence field (in FITS format) and returns a ConvergenceMap instance with the loaded map
Parameters: - realization (int.) – the specific realization to read
- kwargs – the keyword arguments are passed to the getNames method
Returns: ConvergenceMap instance with the loaded map
-
squeeze
(with_ns=False)¶ Returns the cosmological parameters of the model in numpy array form
Parameters: with_ns (bool.) – if True returns also ns as the last parameter Returns: numpy array (Om0,w0,si8,ns–optionally)
-
classmethod
-
class
lenstools.simulations.
CFHTcov
(H0=72.0, Om0=0.26, w0=-1.0, sigma8=0.798, ns=0.96, root_path=None, name=None)[source]¶ Class handler of the weak lensing CFHTcov simulations set, inherits from CFHTemu1; this simulation suite contains 1000 realizations for each of the 13 CFHT subfields, based on 50 independent N-body simulations of a fiducial LambdaCDM universe. Useful to measure accurately descriptor covariance matrices
-
classmethod
getModels
(root_path='/default')[source]¶ On call, this class method returns a CFHTcov instance initialized with the cosmological parameters of the only available model in the suite
Parameters: root_path (str.) – path of your CFHTcov simulations copy Returns: CFHTcov instance initialized with the fiducial cosmological parameters
-
getNames
(realizations, subfield=1, smoothing=0.5)[source]¶ Get the full name of the CFHTcov maps, once a subfield and smoothing scale are specified
Parameters: - subfield (int.) – the specific CFHT subfield you want to retrieve, must be between 1 and 13
- smoothing (float.) – smoothing scale of the maps you wish to retrieve, in arcmin
- realizations (list. or int.) – list of realizations to get the names of, the elements must be in [1,1000]
-
load
(realization, **kwargs)¶ Reads in a specific realization of the convergence field (in FITS format) and returns a ConvergenceMap instance with the loaded map
Parameters: - realization (int.) – the specific realization to read
- kwargs – the keyword arguments are passed to the getNames method
Returns: ConvergenceMap instance with the loaded map
-
squeeze
(with_ns=False)¶ Returns the cosmological parameters of the model in numpy array form
Parameters: with_ns (bool.) – if True returns also ns as the last parameter Returns: numpy array (Om0,w0,si8,ns–optionally)
-
classmethod
Simulation design¶
-
class
lenstools.simulations.
Design
[source]¶ A class that proves useful in designing simulation sets: the main functionality provided is the uniform sampling of an arbirtarily high dimensional parameter space. The points in parameter space are chosen to be as spread as possible by minimizing a cost function, but enforcing a latin hypercube structure, in which each parameter value appears once
-
add_parameter
(parameter_name, min, max, label)[source]¶ Add a dimension to the design by specifying a parameter name, a range and a parameter label (can be in tex format)
Parameters: - parameter_name (str.) – the name of the parameter
- min (float.) – the lower range of the sample interval
- max (float.) – the higher range of the sample interval
- label (str.) – the parameter label you want displayed on a plot, can be in tex format
-
cost
(p, Lambda)[source]¶ Computes the cost function of the current configuration given the metric parameters (p,Lambda)
Parameters: - Lambda (float.) – metric parameter of the cost function; if set to 1.0 the cost function corresponds is the Coulomb potential energy
- p (float.) – metric parameter of the cost function; if set to 2.0 the distances between points are the Euclidean ones
Returns: the value of the cost function
-
diagonalCost
(Lambda)[source]¶ Computes the cost function of a diagonal configuration with a specified number of points and a metric parameter lambda; the cost function is calculated on the scaled parameter values, which are always between 0 and 1
Parameters: Lambda (float.) – metric parameter of the cost function; if set to 1.0 the cost function corresponds is the Coulomb potential energy Returns: the value of the cost function for a diagonal configuration
-
classmethod
load
(filename, labels)[source]¶ Load a pre-computed design from a file, only numpy format supported so far
Parameters: - filename (str. or array) – name of the file from which to load the design, or numpy array with the points
- labels (list.) – labels of the cosmological parameters included in the design
Returns: new Design instance
-
put_points
(npoints)[source]¶ Lay down a number of points on the empty Design: the points are initially layed down on the diagonal of the hypercube
Parameters: npoints (int.) – the number of points to lay down
-
sample
(p, Lambda, seed=0, maxIterations=10000)[source]¶ Evenly samples the parameter space by minimizing the cost function computed with the metric parameters (p,Lambda)
Parameters: - Lambda (float.) – metric parameter of the cost function; if set to 1.0 the cost function corresponds is the Coulomb potential energy
- p (float.) – metric parameter of the cost function; if set to 2.0 the distances between points are the Euclidean ones
- seed (int.) – random seed with which the sampler random generator will be initialized
- maxIterations (int.) – maximum number of iterations that the sampler can perform before stopping
Returns: the relative change of the cost function the last time it varied during the sampling
-
savefig
(filename)[source]¶ Save the visualization to an external file
Parameters: filename (str.) – name of the file on which to save the plot
-
set_title
(title)[source]¶ Give a title to the visualized design
Parameters: title (str.) – title of the figure
-
visualize
(fig=None, ax=None, parameters=None, **kwargs)[source]¶ Visualize the design configuration using matplotlib
Parameters: parameters (list.) – the parameters to visualize, you can specify two or three of them, by their names. If None, all parameters are visualized
-
write
(filename=None, max_rows=None, format='ascii.latex', column_format='{0:.3f}', **kwargs)[source]¶ Outputs the points that make up the design in a nicely formatted table
Parameters: - filename (str. or file descriptor) – name of the file to which the table will be saved; if None the contents will be printed
- max_rows (int.) – maximum number of rows in the table, if smaller than the number of points the different chunks are hstacked (useful if there are too many rows for display)
- format (str.) – passed to the Table.write astropy method
- column_format (str.) – format specifier for the numerical values in the Table
- kwargs (dict.) – the keyword arguments are passed to astropy.Table.write method
Returns: the Table instance with the design parameters
-
Nicaea bindings¶
-
class
lenstools.simulations.
NicaeaSettings
[source]¶ Class handler of the code settings (non linear modeling, tomography, transfer function, etc...)
-
knobs
¶ Lists available settings to tune
-
-
class
lenstools.simulations.
Nicaea
(H0=72.0, Om0=0.26, Ode0=0.74, Ob0=0.046, w0=-1.0, wa=0.0, sigma8=0.798, ns=0.96, name=None)[source]¶ Main class handler for the python bindings of the NICAEA cosmological code, written by M. Kilbinger & collaborators
-
convergencePowerSpectrum
(ell, z=2.0, distribution=None, distribution_parameters=None, settings=None, **kwargs)[source]¶ Computes the convergence power spectrum for the given cosmological parameters and redshift distribution using NICAEA
Parameters: - ell (array.) – multipole moments at which to compute the power spectrum
- z (float., array or None) – redshift bins for the sources; if a single float is passed, single redshift is assumed
- distribution (None,callable or list) – redshift distribution of the sources (normalization not necessary); if None a single redshift is expected; if callable, z must be an array, and a single redshift bin is considered, with the galaxy distribution specified by the call of distribution(z); if a list is passed, each element must be a NICAEA type
- distribution_parameters (str. or list.) – relevant only when distribution is a list or callable. When distribution is callable, distribution_parameters has to be one between “one” and “all” to decide if one or multiple redshift bins have to be considered. If it is a list, each element in it should be the tuple of parameters expected by the correspondent NICAEA distribution type
- settings (NicaeaSettings instance) – NICAEA code settings
- kwargs (dict.) – the keyword arguments are passed to the distribution, if callable
Returns: ( NlxNz array ) computed power spectrum at the selected multipoles (when computing the cross components these are returned in row major C ordering)
-
classmethod
fromCosmology
(cosmo)[source]¶ Builds a Nicaea instance from one of astropy.cosmology objects, from which it inherits all the cosmological parameter values
Parameters: cosmo (astropy FLRW) – one of astropy cosmology instances Returns: Nicaea instance with the cosmological parameters inherited from cosmo
-
shearTwoPoint
(theta, z=2.0, distribution=None, distribution_parameters=None, settings=None, kind='+', **kwargs)[source]¶ Computes the shear two point function for the given cosmological parameters and redshift distribution using NICAEA
Parameters: - theta (array. with units) – angles at which to compute the two point function
- z (float., array or None) – redshift bins for the sources; if a single float is passed, single redshift is assumed
- distribution (None,callable or list) – redshift distribution of the sources (normalization not necessary); if None a single redshift is expected; if callable, z must be an array, and a single redshift bin is considered, with the galaxy distribution specified by the call of distribution(z); if a list is passed, each element must be a NICAEA type
- distribution_parameters (str. or list.) – relevant only when distribution is a list or callable. When distribution is callable, distribution_parameters has to be one between “one” and “all” to decide if one or multiple redshift bins have to be considered. If it is a list, each element in it should be the tuple of parameters expected by the correspondent NICAEA distribution type
- settings (NicaeaSettings instance) – NICAEA code settings
- kind (str.) – must be “+” or “-“
- kwargs (dict.) – the keyword arguments are passed to the distribution, if callable
Returns: ( NtxNz array ) computed two point function at the selected angles (when computing the cross components these are returned in row major C ordering)
-
Gadget2 snapshot handling¶
-
class
lenstools.simulations.
Gadget2Snapshot
(fp=None, pool=None, length_unit=<Quantity 1.0 kpc>, mass_unit=<Quantity 10000000000.0 solMass>, velocity_unit=<Quantity 1.0 km / s>)[source]¶ A class that handles Gadget2 snapshots, mainly I/O from the binary format and spatial information statistics
-
cutPlaneAdaptive
(normal=2, center=<Quantity 7.0 Mpc>, left_corner=None, plane_resolution=<Quantity 0.1 Mpc>, neighbors=64, neighborDistances=None, kind='density', projectAll=False)[source]¶ Cuts a density (or gravitational potential) plane out of the snapshot by computing the particle number density using an adaptive smoothing scheme; the plane coordinates are cartesian comoving
Parameters: - normal (int. (0,1,2)) – direction of the normal to the plane (0 is x, 1 is y and 2 is z)
- center (float. with units) – location of the plane along the normal direction
- plane_resolution (float. with units (or int.)) – plane resolution (perpendicular to the normal)
- left_corner (tuple of quantities or None) – specify the position of the lower left corner of the box; if None, the minimum of the (x,y,z) of the contained particles is assumed
- neighbors (int.) – number of nearest neighbors to use in the adaptive smoothing procedure
- neighborDistances (array with units) – precomputed distances of each particle to its N-th nearest neighbor; if None these are computed
- kind (str. (“density” or “potential”)) – decide if computing a density or gravitational potential plane (this is computed solving the poisson equation)
- projectAll (bool.) – if True, all the snapshot is projected on a single slab perpendicular to the normal, ignoring the position of the center
Returns: tuple(numpy 2D array with the computed particle number density (or lensing potential),bin resolution along the axes,number of particles on the plane)
-
cutPlaneAngular
(normal=2, thickness=<Quantity 0.5 Mpc>, center=<Quantity 7.0 Mpc>, left_corner=None, plane_lower_corner=<Quantity [ 0., 0.] deg>, plane_size=<Quantity 0.15 deg>, plane_resolution=<Quantity 1.0 arcmin>, thickness_resolution=<Quantity 0.1 Mpc>, smooth=None, tomography=False, kind='density', space='real')[source]¶ Same as cutPlaneGaussianGrid(), except that this method will return a lens plane as seen from an observer at z=0; the spatial transverse units are converted in angular units as seen from the observer
Parameters: - normal (int. (0,1,2)) – direction of the normal to the plane (0 is x, 1 is y and 2 is z)
- thickness (float. with units) – thickness of the plane
- center (float. with units) – location of the plane along the normal direction; it is assumed that the center of the plane is seen from an observer with a redshift of self.header[“redshift”]
- left_corner (tuple of quantities or None) – specify the position of the lower left corner of the box; if None, the minimum of the (x,y,z) of the contained particles is assumed
- plane_lower_corner (float with units.) – lower left corner of the plane, as seen from the observer (0,0) corresponds to the lower left corner of the snapshot
- plane_size (float with units) – angular size of the lens plane (angles start from 0 in the lower left corner)
- plane_resolution (float. with units (or int.)) – plane angular resolution (perpendicular to the normal)
- thickness_resolution (float. with units (or int.)) – plane resolution (along the normal)
- smooth (int. or None) – if not None, performs a smoothing of the angular density (or potential) with a gaussian kernel of scale “smooth x the pixel resolution”
- tomography (bool.) – if True returns the lens plane angular density for each slab, otherwise a projected density (or lensing potential) is computed
- kind (str. (“density” or “potential”)) – decide if computing an angular density or lensing potential plane (this is computed solving the poisson equation)
- space (str.) – if “real” return the lens plane in real space, if “fourier” the Fourier transform is not inverted
Returns: tuple(numpy 2D or 3D array with the (unsmoothed) particle angular number density,bin angular resolution, total number of particles on the plane); the constant spatial part of the density field is subtracted (we keep the fluctuation only)
-
cutPlaneGaussianGrid
(normal=2, thickness=<Quantity 0.5 Mpc>, center=<Quantity 7.0 Mpc>, plane_resolution=<Quantity 0.1 Mpc>, left_corner=None, thickness_resolution=<Quantity 0.1 Mpc>, smooth=None, tomography=False, kind='density')[source]¶ Cuts a density (or gravitational potential) plane out of the snapshot by computing the particle number density on a slab and performing Gaussian smoothing; the plane coordinates are cartesian comoving
Parameters: - normal (int. (0,1,2)) – direction of the normal to the plane (0 is x, 1 is y and 2 is z)
- thickness (float. with units) – thickness of the plane
- center (float. with units) – location of the plane along the normal direction
- plane_resolution (float. with units (or int.)) – plane resolution (perpendicular to the normal)
- left_corner (tuple of quantities or None) – specify the position of the lower left corner of the box; if None, the minimum of the (x,y,z) of the contained particles is assumed
- thickness_resolution (float. with units (or int.)) – plane resolution (along the normal)
- smooth (int. or None) – if not None, performs a smoothing of the density (or potential) with a gaussian kernel of scale “smooth x the pixel resolution”
- kind (str. (“density” or “potential”)) – decide if computing a density or gravitational potential plane (this is computed solving the poisson equation)
Returns: tuple(numpy 3D array with the (unsmoothed) particle number density,bin resolution along the axes, number of particles on the plane)
-
getID
(first=None, last=None, save=True)[source]¶ Reads in the particles IDs, 4 byte ints, (read in of a subset is allowed): when first and last are specified, the numpy array convention is followed (i.e. getID(first=a,last=b)=getID()[a:b])
Parameters: - first (int. or None) – first particle in the file to be read, if None 0 is assumed
- last (int. or None) – last particle in the file to be read, if None the total number of particles is assumed
- save (bool.) – if True saves the particles IDs as attribute
Returns: numpy array with the particle IDs
-
getPositions
(first=None, last=None, save=True)[source]¶ Reads in the particles positions (read in of a subset is allowed): when first and last are specified, the numpy array convention is followed (i.e. getPositions(first=a,last=b)=getPositions()[a:b])
Parameters: - first (int. or None) – first particle in the file to be read, if None 0 is assumed
- last (int. or None) – last particle in the file to be read, if None the total number of particles is assumed
- save (bool.) – if True saves the particles positions as attribute
Returns: numpy array with the particle positions
-
getVelocities
(first=None, last=None, save=True)[source]¶ Reads in the particles velocities (read in of a subset is allowed): when first and last are specified, the numpy array convention is followed (i.e. getVelocities(first=a,last=b)=getVelocities()[a:b])
Parameters: - first (int. or None) – first particle in the file to be read, if None 0 is assumed
- last (int. or None) – last particle in the file to be read, if None the total number of particles is assumed
- save (bool.) – if True saves the particles velocities as attrubute
Returns: numpy array with the particle velocities
-
gridID
()[source]¶ Compute an ID for the particles in incresing order according to their position on a Nside x Nside x Nside grid; the id is computed as x + y*Nside + z*Nside**2
Returns: the gridded IDs Return type: array of float
-
header
¶ Displays the snapshot header information
Returns: the snapshot header information in dictionary form Return type: dict.
-
lensMaxSize
()[source]¶ Computes the maximum observed size of a lens plane cut out of the current snapshot
-
neighborDistances
(neighbors=64)[source]¶ Find the N-th nearest neighbors to each particle
Parameters: neighbors (int.) – neighbor order Returns: array with units
-
numberDensity
(resolution=<Quantity 0.5 Mpc>, smooth=None, left_corner=None, save=False)[source]¶ Uses a C backend gridding function to compute the particle number density for the current snapshot: the density is evaluated using a nearest neighbor search
Parameters: - resolution (float with units or int.) – resolution below which particles are grouped together; if an int is passed, this is the size of the grid
- smooth (int. or None) – if not None, performs a smoothing of the density (or potential) with a gaussian kernel of scale “smooth x the pixel resolution”
- left_corner (tuple of quantities or None) – specify the position of the lower left corner of the box; if None, the minimum of the (x,y,z) of the contained particles is assumed
- save (bool.) – if True saves the density histogram and resolution as instance attributes
Returns: tuple(numpy 3D array with the (unsmoothed) particle number density on a grid,bin resolution along the axes)
-
classmethod
open
(filename, pool=None)[source]¶ Opens a gadget snapshot at filename
Parameters: - filename (str. or file.) – file name of the gadget snapshot
- pool (MPIWhirlPool instance) – use to distribute the calculations on different processors; if not None, each processor takes care of one of the snapshot parts, appending as ”.n” to the filename
-
pos2R
(filename, variable_name='pos')[source]¶ Saves the positions of the particles in a R readable format, for facilitating visualization with RGL
Parameters: - filename (str.) – name of the file on which to save the particles positions
- variable_name (str.) – name of the variable that contains the (x,y,z) positions in the R environment
-
powerSpectrum
(k_edges, resolution=None)[source]¶ Computes the power spectrum of the relative density fluctuations in the snapshot at the wavenumbers specified by k_edges; a discrete particle number density is computed before hand to prepare the FFT grid
Parameters: - k_edges (array.) – wavenumbers at which to compute the density power spectrum (must have units)
- resolution (float with units, int. or None) – optional, fix the grid resolution to some value; to be passed to the numberDensity method. If none this is computed automatically from the k_edges
Returns: tuple(k_values(bin centers),power spectrum at the specified k_values)
-
savefig
(filename)[source]¶ Save the snapshot visualization to an external file
Parameters: filename (str.) – file name to which the figure will be saved
-
setHeaderInfo
(Om0=0.26, Ode0=0.74, w0=-1.0, wa=0.0, h=0.72, redshift=100.0, box_size=<Quantity 20.833333333333336 Mpc>, flag_cooling=0, flag_sfr=0, flag_feedback=0, flag_stellarage=0, flag_metals=0, flag_entropy_instead_u=0, masses=<Quantity [ 0.00000000e+00, 1.03000000e+10, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00] solMass>, num_particles_file_of_type=None, npartTotalHighWord=array([0, 0, 0, 0, 0, 0], dtype=uint32))[source]¶ Sets the header info in the snapshot to write
-
setPositions
(positions)[source]¶ Sets the positions in the current snapshot (with the intent of writing them to a properly formatted snapshot file)
Parameters: positions ((N,3) array with units) – positions of the particles, must have units
-
setVelocities
(velocities)[source]¶ Sets the velocities in the current snapshot (with the intent of writing them to a properly formatted snapshot file)
Parameters: velocities ((N,3) array with units) – velocities of the particles, must have units
-
visualize
(fig=None, ax=None, scale=False, first=None, last=None, **kwargs)[source]¶ Visualize the particles in the Gadget snapshot using the matplotlib 3D plotting engine, the kwargs are passed to the matplotlib scatter method
Parameters: scale (bool.) – if True, multiply all the (comoving) positions by the scale factor
-
write
(filename, files=1)[source]¶ Writes particles information (positions, velocities, etc...) to a properly formatter Gadget snapshot
Parameters: - filename (str.) – name of the file to which to write the snapshot
- files (int.) – number of files on which to split the writing of the snapshot (useful if the number of particles is large); if > 1 the extension ”.n” is appended to the filename
-
writeParameterFile
(filename, settings=<lenstools.simulations.gadget2.Gadget2Settings object>)[source]¶ Writes a Gadget2 parameter file to evolve the current snapshot using Gadget2
Parameters: - filename (str.) – name of the file to which to write the parameters
- settings (Gadget2Settings) – tunable settings of Gadget2 (see Gadget2 manual)
-
Ray tracing simulations¶
-
class
lenstools.simulations.
Plane
(data, angle, redshift=2.0, cosmology=None, comoving_distance=None, unit=Unit("rad2"), num_particles=None, masked=False)[source]¶ -
classmethod
load
(filename, format=None, init_cosmology=True)[source]¶ Loads the Plane from an external file, of which the format can be specified (only fits implemented so far)
Parameters: - filename (str.) – name of the file from which to load the plane
- format (str.) – format of the file, only FITS implemented so far; if None, it’s detected automatically from the filename
- init_cosmology (bool.) – if True, instantiates the cosmology attribute of the PotentialPlane
Returns: PotentialPlane instance that wraps the data contained in the file
-
randomRoll
(seed=None, lmesh=None)[source]¶ Randomly shifts the plane along its axes, enforcing periodic boundary conditions
Parameters: - seed (int.) – random seed with which to initialize the generator
- lmesh (array) – the FFT frequency meshgrid (lx,ly) necessary for the calculations in fourier space; if None, a new one is computed from scratch (must have the appropriate dimensions)
-
save
(filename, format=None, double_precision=False)[source]¶ Saves the Plane to an external file, of which the format can be specified (only fits implemented so far)
Parameters: - filename (str.) – name of the file on which to save the plane
- format (str.) – format of the file, only FITS implemented so far; if None, it’s detected automatically from the filename
- double_precision (bool.) – if True saves the Plane in double precision
-
classmethod
-
class
lenstools.simulations.
DensityPlane
(data, angle, redshift=2.0, cosmology=None, comoving_distance=None, unit=Unit("rad2"), num_particles=None, masked=False)[source]¶ Class handler of a lens density plane, inherits from the parent Plane class; additionally it defines redshift and comoving distance attributes that are needed for ray tracing operations
-
potential
(lmesh=None)[source]¶ Computes the lensing potential from the density plane solving the Poisson equation via FFTs
Parameters: lmesh (array) – the FFT frequency meshgrid (lx,ly) necessary for the calculations in fourier space; if None, a new one is computed from scratch (must have the appropriate dimensions) Returns: PotentialPlane instance with the computed lensing potential
-
-
class
lenstools.simulations.
PotentialPlane
(data, angle, redshift=2.0, cosmology=None, comoving_distance=None, unit=Unit("rad2"), num_particles=None, masked=False)[source]¶ Class handler of a lens potential plane, inherits from the parent Plane class; additionally it defines redshift and comoving distance attributes that are needed for ray tracing operations
-
deflectionAngles
(x=None, y=None, lmesh=None)[source]¶ Computes the deflection angles for the given lensing potential by taking the gradient of the potential map; it is also possible to proceed with FFTs
Parameters: - x (array with units) – optional; if not None, compute the deflection angles only for rays hitting the lens at the particular x positions (mainly for speedup in case there are less light rays than the plane resolution allows; must proceed in real space to allow speedup)
- y (array with units) – optional; if not None, compute the deflection angles only for rays hitting the lens at the particular y positions (mainly for speedup in case there are less light rays than the plane resolution allows; must proceed in real space to allow speedup)
- lmesh (array) – the FFT frequency meshgrid (lx,ly) necessary for the calculations in fourier space; if None, a new one is computed from scratch (must have the appropriate dimensions)
Returns: DeflectionPlane instance, or array with deflections of rays hitting the lens at (x,y)
-
density
(x=None, y=None)[source]¶ Computes the projected density fluctuation by taking the laplacian of the potential; useful to check if the potential is reasonable
Parameters: - x (array with units) – optional; if not None, compute the density only for rays hitting the lens at the particular x positions (mainly for speedup in case there are less light rays than the plane resolution allows; must proceed in real space to allow speedup)
- y (array with units) – optional; if not None, compute the density only for rays hitting the lens at the particular y positions (mainly for speedup in case there are less light rays than the plane resolution allows; must proceed in real space to allow speedup)
Returns: DensityPlane instance with the density fluctuation data (if x and y are None), or numpy array with the same shape as x and y
-
shearMatrix
(x=None, y=None, lmesh=None)[source]¶ Computes the shear matrix for the given lensing potential; it is also possible to proceed with FFTs
Parameters: - x (array with units) – optional; if not None, compute the shear matrix only for rays hitting the lens at the particular x positions (mainly for speedup in case there are less light rays than the plane resolution allows; must proceed in real space to allow speedup)
- y (array with units) – optional; if not None, compute the shear matrix only for rays hitting the lens at the particular y positions (mainly for speedup in case there are less light rays than the plane resolution allows; must proceed in real space to allow speedup)
- lmesh (array) – the FFT frequency meshgrid (lx,ly) necessary for the calculations in fourier space; if None, a new one is computed from scratch (must have the appropriate dimensions)
Returns: ShearTensorPlane instance, or array with deflections of rays hitting the lens at (x,y)
-
-
class
lenstools.simulations.raytracing.
DeflectionPlane
(data, angle, redshift=2.0, comoving_distance=None, cosmology=None, unit=Unit("rad"))[source]¶ Class handler of a lens deflection plane, inherits from the parent Spin1 class and holds the values of the deflection angles of the light rays that cross a potential plane
-
convergence
(precision='first')[source]¶ Computes the convergence from the deflection angles by taking the appropriate components of the jacobian
Parameters: precision (str.) – if “first” computes the convergence at first order in the lensing potential (only one implemented so far) Returns: ConvergenceMap instance with the computed convergence values
-
jacobian
()[source]¶ Computes the jacobian of the deflection angles, useful to compute shear and convergence; units are handled properly
Returns: the jacobian of the deflection field in array form, of shape (4,:,:) where the four components are, respectively, xx,xy,yx,yy
-
omega
(precision='first')[source]¶ Computes the omega field (i.e. the real space B mode) from the deflection angles by taking the appropriate components of the jacobian
Parameters: precision (str.) – if “first” computes omega at first order in the lensing potential (only one implemented so far) Returns: Spin0 instance with the computed omega values
-
shear
(precision='first')[source]¶ Computes the shear from the deflection angles by taking the appropriate components of the jacobian
Parameters: precision (str.) – if “first” computes the shear at first order in the lensing potential (only one implemented so far) Returns: ShearMap instance with the computed convergence values
-
-
class
lenstools.simulations.raytracing.
ShearTensorPlane
(data, angle, redshift=2.0, comoving_distance=None, cosmology=None, unit=Unit(dimensionless))[source]¶ Class handler of a plane of shear matrices, inherits from the parent Spin2 class and holds the 3 values of the symmetric shear matrix (2 diagonal + 1 off diagonal), for each pixel
-
class
lenstools.simulations.
RayTracer
(lens_mesh_size=None)[source]¶ Class handler of ray tracing operations: it mainly computes the path corrections of light rays that travel through a set of gravitational lenses
-
addLens
(lens_specification)[source]¶ Adds a gravitational lens to the ray tracer, either by putting in a lens plane, or by specifying the name of a file which contains the lens specifications
Parameters: lens_specification – specifications of the lens to add, either in tuple(filename,distance,redshift) or as a PotentialPlane instance
-
convergenceDirect
(initial_positions, z=2.0, save_intermediate=False)[source]¶ Computes the convergence directly integrating the lensing density along the unperturbed line of sight
Parameters: - initial_positions (numpy array or quantity) – initial angular positions of the light ray bucket, according to the observer; if unitless, the positions are assumed to be in radians. initial_positions[0] is x, initial_positions[1] is y
- z (float.) – redshift of the sources
- save_intermediate (bool.) – save the intermediate values of the convergence as successive lenses are crossed
Returns: convergence values at each of the initial positions
-
displayRays
(initial_positions, z=2.0, projection='2d', fig=None, ax=None)[source]¶ Uses matplotlib to display a visual of the lens system and of the deflection that the light rays which traverse it undergo
param initial_positions: initial angular positions of the light ray bucket, according to the observer; if unitless, the positions are assumed to be in radians. initial_positions[0] is x, initial_positions[1] is y :type initial_positions: numpy array or quantity
Parameters: - z (float. or array) – redshift of the sources; if an array is passes, a redshift must be specified for each ray, i.e. z.shape==initial_positions.shape[1:]
- projection (str.) – can be “2d” for the projections of the ray positions along x and y or “3d” for the full 3d view
- fig (matplotlib figure) – figure object that owns the plot
-
randomRoll
(seed=None)[source]¶ Randomly rolls all the lenses in the system along both axes
Parameters: seed (int.) – random seed with which to initialize the generator
-
reorderLenses
()[source]¶ Reorders the lenses in the ray tracer according to their comoving distance from the observer
-
shoot
(initial_positions, z=2.0, initial_deflection=None, precision='first', kind='positions', save_intermediate=False, compute_all_deflections=False, callback=None, **kwargs)[source]¶ Shots a bucket of light rays from the observer to the sources at redshift z (backward ray tracing), through the system of gravitational lenses, and computes the deflection statistics
Parameters: - initial_positions (numpy array or quantity) – initial angular positions of the light ray bucket, according to the observer; if unitless, the positions are assumed to be in radians. initial_positions[0] is x, initial_positions[1] is y
- z (float. or array) – redshift of the sources; if an array is passed, a redshift must be specified for each ray, i.e. z.shape==initial_positions.shape[1:]
- initial_deflection (numpy array or quantity) – if not None, this is the initial deflection light rays undergo with respect to the line of sight (equivalent to specifying the first derivative IC on the lensing ODE); must have the same shape as initial_positions
- precision (str.) – precision at which to compute weak lensing quantities, must be “first” for first order in the lensing potential, or “second” for added precision (not functional yet)
- kind (str.) – what deflection statistics to compute; “positions” will calculate the ray deflections after they crossed the last lens, “jacobian” will compute the lensing jacobian matrix after the last lens, “shear” and “convergence” will compute the omonimous weak lensing statistics
- save_intermediate (bool.) – save the intermediate positions of the rays too
- compute_all_deflections (bool.) – if True, computes the gradients of the lensing potential at every pixel on the lens(might be overkill if Nrays<<Npixels); must be True if the computation is done with FFTs
- callback (callable) – if not None, this callback function is called on the current ray positions array at each step in the ray tracing; the current raytracing instance and the step number are passed as additional arguments, hence callback must match this signature
- kwargs (dict.) – the keyword arguments are passed to the callback if not None
Returns: angular positions (or jacobians) of the light rays after the last lens crossing
-
shootForward
(source_positions, z=2.0, save_intermediate=False, grid_resolution=512, interpolation='nearest')[source]¶ Shoots a bucket of light rays from the source at redshift z to the observer at redshift 0 (forward ray tracing) and computes the according deflections using backward ray tracing plus a suitable interpolation scheme (KD Tree based)
Parameters: - source_positions (numpy array or quantity) – angular positions of the unlensed sources
- z (float.) – redshift of the sources
- save_intermediate (bool.) – if True computes and saves the apparent image distortions after each lens is crossed (can be computationally expensive)
- grid_resolution (int.) – the number of points on a side of the interpolation grid (must be choosen big enough according to the number of sources to resolve)
- interpolation (str.) – only “nearest” implemented so far
Returns: apparent positions of the sources as seen from the observer
-
Observation sets¶
-
class
lenstools.observations.
CFHTLens
(root_path=None)[source]¶ Class handler of the CFHTLens reduced data set, already split in 13 3x3 deg^2 subfields
Defaults¶
-
lenstools.defaults.
load_fits_default_convergence
(filename)[source]¶ This is the default convergence fits file loader, it assumes that the two components of the shear are stored in two different image FITS files, which have an ANGLE keyword in the header
Parameters: gamma_file – Name of the FITS file that contains the shear map Returns: tuple – (angle,ndarray – kappa; kappa is the convergence map) Raises: IOError if the FITS files cannot be opened or do not exist
-
lenstools.defaults.
load_fits_default_shear
(filename)[source]¶ This is the default shear fits file loader, it assumes that the two components of the shear are stored in a single image FITS file, which have an ANGLE keyword in the header
Parameters: gamma_file – Name of the FITS file that contains the shear map Returns: tuple – (angle,ndarray – gamma; gamma[0] is the gamma1 map, gamma[1] is the gamma2 map); the maps must follow matrix ordering, i.e. the first axis (0) is y and the second axis (1) is x. This matters for the E/B mode decomposition Raises: IOError if the FITS files cannot be opened or do not exist
-
lenstools.defaults.
default_callback_loader
(filename, l_edges)[source]¶ Default ensemble loader: reads a FITS data file containing a convergence map and measures its power spectrum
Parameters: args (Dictionary) – A dictionary that contains all the relevant parameters as keys. Must have a “map_id” key Returns: ndarray of the measured statistics Raises: AssertionError if the input dictionary doesn’t have the required keywords
Limber integration¶
-
class
lenstools.
LimberIntegrator
(cosmoModel=FlatLambdaCDM(name="WMAP9", H0=69.3 km / (Mpc s), Om0=0.286, Tcmb0=2.725 K, Neff=3.04, m_nu=[ 0. 0. 0.] eV, Ob0=0.0463))[source]¶ A 3D power spectrum integrator that will compute the convergence power spectrum using the Limber approximation.
Parameters: cosmoModel (astropy.cosmology) – One of astropy.cosmology objects (WMAP9 cosmology is set by default) -
convergencePowerSpectrum
(l)[source]¶ Computes the convergence power spectrum with the Limber integral of the 3D matter power spectrum; this still assumes a single source redshift at z0 = max(z)
Parameters: l (array) – multipole moments at which to compute the convergence power spectrum Returns: array – the convergence power spectrum at the l values specified
-
load3DPowerSpectrum
(loader, *args, **kwargs)[source]¶ Loads in the matter power spectrum from (pre-computed) external files; args and kwargs are passed to the loader
Parameters: loader (function) – must return, in order, k, z, P(k,z)
-
setUnits
(kappa_units=None, power_units=None)[source]¶ Set the physical units for wavenumber and matter power spectrum, default for length is Mpc
-
writeCAMBSettings
(l, z, powFileRoot='matterpower', transfer_high_precision=False, transfer_k_per_logint=0, transfer_interp_matterpower=True)[source]¶ Outputs a StringIO object that will contain the redshift settings of the CAMB parameter file that will needed in order for CAMB to produce the linear or non linear matter power spectra that will then be integrated by the computeConvergence() method
Parameters: - l (array) – multipole moments at which to compute the convergence power spectrum
- z (array) – redshift bins at which the matter power spectrum is calculated (assumed to be a monotonic array with more than 1 element)
- powFileRoot (str.) – root of the filename that you want to give to the CAMB power spectrum outputs
- transfer_high_precision (bool.) – read CAMB documentation (this sets the precision of the calculated transfer function)
- transfer_k_per_logint (int.) – read CAMB documentation (this sets the k wavenumber binning)
- transfer_interp_matterpower (bool.) – read CAMB documentation (this sets how the matter power is interpolated between different k’s)
Returns: str – the portion of the CAMB parameter file relevant to the 3D matter power spectrum
-