FRApy: Fitting Resolved Arcs with python¶
FRApy is python code that fits gravitational arcs in image plane with analytical modes, such as metallicity gradients and velocity fields, taking into account the lensing distortions. This is done in a Bayesian framework, by maximising a log-logarithm likelihood function using an MCMC sampler (emcee).
Intro¶
Gravitational arcs¶
Gravitational arcs are galaxies in the backgroud of a cluster of galaxies or a massive galaxy. These massive objects act as a magnifying lens, increasing the size of the backgroud galaxies in the sky, which allows us to resolve them at smaller spatial scales. However, the magnification is not uniform, and the lensed galaxies appear ‘distorted’ in the sky, typically in an arc-like shape (hence the name graviational ‘arcs’).
The spatial distortion makes the analysis of these objects more difficult, especially when several images of the same object (called multiple images) are available. FRApy deals with this issue through forward-modelling: we start with a model in source plane (the ‘undistorted’ galaxy), lens it to image plane using a lensing model and compare it with the data.
What you’ll need¶
- Deflection maps. These are 2D fits images produced from a lensing model that map how much deflection (in arcseconds along both directions RA DEC, 1 map each) a little photon landing in position RA,DEC in the sky has suffered due to the gravitational lensing effect. This is the core of FRApy’s modelling and allows one to generate lensed images that can be compared with the observations. One does this by calculating where each pixel in image plane originated from in the source plane. A number of Frontier Fields models provide these maps you can check them here here (example given for Abell 370). It is also possible to create these maps using LENSTOOL using the dpl command.
2. An analytical model that you hope describes your data, such as a gradient for metallicity or an arctangent model velocity fields. FRApy comes with a number of these models, that we used in our own publications, but we also show how you can create and fit your own model.
Install¶
>> pip install -i https://test.pypi.org/simple/ FRApy
The project is also available GitHub.
It requires the following packages:
- numpy==1.15.4
- matplotlib==3.0.2
- astropy==3.1
- reproject==0.4
- emcee==2.2.1
- pickle==4.0
- corner==2.0.1
Quick start¶
FRApy is a python module. A minimum working example would look something like this:
# Import FRApy
from frapy import Observation,Metallicity_Gradient,Output
from frapy import fit_model,make_input_parameters
# Load Observations
obs = Observation(z=0.611,
data_path='Demo_data/AS1063_map_metallicity.fits',
unc_path='Demo_data/AS1063_map_metallicity_unc.fits',
seeing = 1.03/0.2)
# Choose a Model, in this case a linear metallicity gradient
model = Metallicity_Gradient(zlens=0.322,
dfx_path='Demo_data/AS1063_dplx.fits',
dfy_path='Demo_data/AS1063_dply.fits')
model.create_projection_maps(obs)
# Fit the data
input_par = make_input_parameters(name = ('cx', 'cy', 'q', 'pa', 'z_grad', 'z_0'),
value = ( 29, 23, 0.7, 20, -0.02, 9.0),
minimum = ( 28, 22, 0.4, -20, -0.1, 8.5),
maximum = ( 33, 27, 0.9, 90, 0.0, 9.5))
out = fit_model(obs,model,input_par,'output_file',nsteps=2000,nwalkers=24)
# Inspect the fit
results = Output('output_file')
results.best_parameters()
Mode in depth demos¶
We have included demo data used in Patricio et al. in prep. and two notebooks with examples on how to fit your data.
Authors and Citations¶
This code was developed by:
- Vera Patricio (vera.patricio@dark-cosmology.dk)
- Johan Richard, lensing specialist
If you use FRApy in your science, please add the following citation:
and
Patrício et. al, in prep.
and don’t forget astropy and emcee!
Observations¶
-
class
frapy.observations.
Observation
(data_path, z, unc_path=None, seeing=0)[source]¶ Data to be fitted.
This class handles the data that is going to be fitted with the models. It includes the image (data) to be fit, the associated uncertainty (optionally ), the redshift of the object and the seeing.
Parameters: - data_path (str) – Path to the fits file containing the data.
- unc_path (str) – Path to the fits file containing the associated uncertainty (optional).
- z (float) – Redshift.
- seeing (float) – Seeing (in pixels). To be used in the model convolution.
- data (float array) – Data to be fitted. It is read from the data_path
- unc (float array) – Data uncertainted. It is read from the unc_path.
Models¶
The BaseModel handles all the lensing part, producing a distance map that is used by all the other models.
-
class
frapy.models.
BaseModel
(zlens, dfx_path, dfy_path, df_ang=0, cx=0, cy=0, q=1, pa=0)[source]¶ Global lensing model to be used in all other Model classes .
This class prepares the deflection maps to be used with a particular object (i.e. at a particular redshift) and observations (i.e. aligns the maps with the data).
The main output is a distance map, in kiloparsecs, and an azimuthal map that serve as base for all the models being fit (metallicity gradient, velocity…)
Parameters: - z_lens (float) – Redshift of the gravitational lens.
- dfx_path (str) – Path to the fits file with the x deflection.
- dfy_path (str) – Path to the fits file with the y deflection.
- cx (int) – x position of the centre (in pixels)
- cy (int) – y position of the centre (in pixels)
- q (float) – axis ratio (b/a)
- pa (float) – Position angle (0 North, +90 East )
- df_ang (float) – Angle between x axis and North (measured anti-clockwise) in the deflection maps.
- project_x (float array) – Lensing model (deflection in x direction) to be used to a particular object. Created with the ‘create_deflection_maps_for_object’ method.
- project_y (float array) – Lensing model (deflection in y direction) to be used to a particular object. Created with the ‘create_deflection_maps_for_object’ method.
- data (array) – An array with a realisation of a model made from the current parameter values.
- conv_data (array) – An array with a realisation of a model made from the current parameter values, convolved by the seeing of observations.
-
convolve_with_seeing
(seeing)[source]¶ Convolves a model with a Gaussian with width (sigma) ‘seeing’.
-
create_projection_maps
(Observation, correct_z=True)[source]¶ Takes the more global deflection maps produced by a graviatational lensing fitting code, and converts these maps to ‘projection’ maps, that maps where a pixel in source plane should be ‘projected’ in image plane, for this particular Observation. The project_x and project_y attributes are created with this function.
-
make_azimuthal_map
()[source]¶ Produces an azimuthal map, in kpc, centrered in ‘cx’,’cy’ and assuming a ratio of ‘q’ between the minor and major axis, with the major axis in the ‘pa’ direction.
-
class
frapy.models.
Metallicity_Gradient
(zlens, dfx_path, dfy_path, df_ang=0, cx=0, cy=0, q=1, pa=0, z_grad=-1, z_0=0)[source]¶ Linear metallicity gradient.
This model inherits the distance maps attributes (cx,cy,q and pa), from which the metallicity at each point is calculated assuming a gradient and a central metallicity value:
Z(r) = Delta Z * r + Z_0
with r the radius in kpc, Delta Z the gradient in dex/kpc, Z_0 the central metallicity.
Parameters:
-
class
frapy.models.
Metallicity_Gradient_Constant_Centre
(zlens, dfx_path, dfy_path, cx=0, cy=0, q=1, pa=0, z_grad=-1, z_0=0, r_flat=0.5, z_grad_inner=-1)[source]¶ Linear metallicity gradient with a flatenning of the centre at r_flat.
This model inherits the distance maps attributes (cx,cy,q and pa), from which the metallicity at each point is calculated assuming a gradient and a central metallicity value:
Z(r) = Delta Z * r + Z_0
with r the radius in kpc, Delta Z the gradient in dex/kpc, Z_0 the central metallicity.
For r < r_flat, the metallicity is constant.
Parameters: - cx (int) – x position of the centre (in pixels)
- cy (int) – y position of the centre (in pixels)
- q (float) – axis ratio (a/b)
- pa (float) – Position angle (0 North, +90 East )
- z_grad (float) – Gradient in dex/kpc.
- z_0 (float) – Central metallicity value (value at cx,cy)
- r_flat (float) – Radius that delimits the central zone where the metallicity is flat
-
class
frapy.models.
Velocity_Arctangent
(zlens, dfx_path, dfy_path, cx=0, cy=0, q=1, pa=0, v_t=100, r_t=10)[source]¶ Exponential velocity model.
This model inherits the distance and azimuthal maps, from which an arctangent model of the velocity at each point is calculated assuming the following formulae:
V(r) = v_t frac{2}{pi} arctan (frac{2r}{r_t})
with r the radius in kpc, v_t the terminal velocity and r_t the transition radius.
Parameters:
Fitting¶
-
frapy.fit_model.
fit_model
(obs, model, parameters, outputname, nsteps=1000, nwalkers=24, mask=None, binning_map=None)[source]¶ Routine that fits the observations using a given model and the emcee sampler.
We make use of the emcee sampler (http://dfm.io/emcee/current/) to fit the free parameters of the model to the observations. We are maximising the following log-probabiluty function:
$ln(probability) = ln(priors) + ln(likelihood)$
with the log likelohood function as:
$ln(likelihood) = -frac{1}{2} ( frac{(data-model)^2}{uncertainty^2} + ln(2 pi uncertainty^2))$
Both the model and the observations should be instances of the Observations and BaseModel classes from frapy.
The parameters input is a dictionary in the shape:
parameters = {parameter_name1:{‘value’:X, ‘min’:Y, ‘max’:Z}, parameter_name2:{‘value’:A, ‘min’:B, ‘max’:C},…}
where the parameter_name variables should correspond to the parameters in the model being used; value is the starting value of each parameter; and min and max the minimum and maximum values allowed. We assume uniform priors to all parameters (i.e. all values between min and max have the same prior probability). Parameters not present in this dictionary will not be varied and will be kept to the value of the input model.
It is possible to mask part of the data out by using a mask. This should be a 2D array, of the same shape as the data, with only zeros (masked values) and ones (valid values). The maximisation will be made using only the valid values.
If the data was binned when deriving the quantity being fit, i.e. if pixels were grouped and analysed as a single pixel and that value taken as the value of all the pixels grouped, it is possible to include this information using a binning_map. This should be a 2D array in which pixels of the same bin are given the same value. Pixels in the image that were not analysed (not included in the binning) should be given negative values. These are not included in the minimisation.
Parameters: - obs (Observation) – An instance of the Observation class
- model (Metallicity_Gradient,Velocity) – A frapy model (based in the BaseModel class)
- parameters (dictionary) – A dictionary containing the parameters of the model to be varied and their limits. Parameters not in this dictionary will not be varied.
- outputname (str) – Name of the output pickle file.
- nsteps (int) – number of steps of the emcee walkers. Default: 1000
- nwalkers (int) – Number of emcee walkers. Default: 24
- mask (array int) – Array of the same shape as the data containing only zeros (masked values) or ones (valid values). Optional.
- binning_map (array int) – Array of the same shape as the data containing encoding the pixels that were groupped togther. Optional.
Returns: - Returns a dictionary with – sampler chain sampler lnprobability parameter names in the correct order input parameters the mask used the binning map used the observations the model
- This is also saved as a pickles file.
Explore the output¶
-
class
frapy.check_fit.
Output
(outfile)[source]¶ Allows the output of fit_model to be inspected
Reads the pickle output and allows to plot:
. the walkers positions at each iteration to check for convergence . a corner plot of the results . the 50th, 16th and 84th percentiles (mean and +/- 1 sigma)Parameters: outfile (str) – The name of the pickle file being inspected (without the ‘.pickle’ extension) -
best_parameters
(start=0)[source]¶ Calculates the 16th, 50th and 84th percentiles for each parameter. Only uses iterations after ‘start’
-
check_convergence
()[source]¶ Plots the walkers positions at each iteration for each parameter as well as the value of the log-likelihood probability for each iteration.
-
goodness_of_fit
(best_parameters)[source]¶ Given a dictionary with parameter names and values, calculates the chi2, reduced chi2 (chi2/dof), the log-likelihood probability and the Bayesian Information Criteria (BIC) for the model with those parameters values.
Parameters: best_parameters (dictionary) – Dictionary in the shape {parameter_name1:{value:X,min:Z,max:Z},parameter_name2:{value:X,min:Z,max:Z}} (from the check_fit.best_parameters function, for example). Returns: chi2/dof – Reduced chi2 Return type: float
-
make_cornerplot
(start=0)[source]¶ Makes a corner plot of the results. Only uses iterations after ‘start’
-
plot_solution
(best_parameters)[source]¶ Given a dictionary with parameter names and values, plots the model.
Parameters: best_parameters (dictionary) – Dictionary in the shape {parameter_name1:{value:X,min:Z,max:Z},parameter_name2:{value:X,min:Z,max:Z}} (from the check_fit.best_parameters function, for example). Returns: - model (array float) – the model with the best parameters
- residuals (array float) – the residuals (data - model)
-
Miscelaneous¶
This module containts miscelaneous functions used in the fitting.
-
frapy.utils.
make_input_parameters
(name, value, minimum, maximum)[source]¶ Outputs a parameter dictionary to be used in the fit_model function.
This dictionary has the following form:
{parameter_name1:{‘value’:X, ‘min’:Y, ‘max’:Z}, parameter_name2:{‘value’:A, ‘min’:B, ‘max’:C}, … }Parameters: - name (array str) – An array of strings containing the names of the model’s parameters to be fitted.
- value (array float) – The initial values of these parameters (ACTUALLY NOT USED)
- minimum (array float) – The minimum value allowed for each parameter
- maximum (array float) – The maximum value allowed for each parameter
Returns: parameter
Return type: dictionary