PyAutoLens

PyAutoLens is software for analysing strong gravitational lenses, an astrophysical phenomenon where a galaxy appears multiple times because its light is bent by the gravitational field of an intervening foreground lens galaxy.

Here is a schematic of a strong gravitational lens:

Schematic of Gravitational Lensing Credit: F. Courbin, S. G. Djorgovski, G. Meylan, et al., Caltech / EPFL / WMKO https://www.astro.caltech.edu/~george/qsolens/

This notebook gives an overview of PyAutoLens's features and API.

Lets first import autolens, its plotting module and the other libraries we'll need.

In [1]:
%matplotlib inline

import autolens as al
import autolens.plot as aplt

import matplotlib.pyplot as plt
from os import path

Lensing

First, we illustrate lensing calculations in PyAutoLens by creating an an image of a strong lens.

To describe the deflection of light, PyAutoLens uses Grid2D data structures, which are two-dimensional Cartesian grids of (y,x) coordinates.

Below, we make and plot a uniform Cartesian grid:

In [2]:
grid_2d = al.Grid2D.uniform(
    shape_native=(150, 150),
    pixel_scales=0.05,  # <- The pixel-scale describes the conversion from pixel units to arc-seconds.
)

grid_2d_plotter = aplt.Grid2DPlotter(grid=grid_2d)
grid_2d_plotter.figure_2d()

Our aim is to create an image of the source galaxy after its light has been deflected by the mass of the foreground lens galaxy. We therefore need to ray-trace the Grid2D's coordinates from the 'image-plane' to the 'source-plane'.

We therefore need analytic functions representing a galaxy's light and mass distributions. For this, PyAutoLens uses Profile objects, for example an ellipical sersic LightProfile object which represents a light distribution:

In [3]:
sersic_light_profile = al.lp.Sersic(
    centre=(0.0, 0.0),
    ell_comps=(0.2, 0.1),
    intensity=0.005,
    effective_radius=2.0,
    sersic_index=4.0,
)

By passing this profile the Grid2D, we can evaluate the light emitted at every (y,x) coordinate on the Grid2D and create an image of the LightProfile.

In [4]:
image_2d = sersic_light_profile.image_2d_from(grid=grid_2d)

plt.imshow(image_2d.native) # The use of 'native' is described at the start of the HowToLens tutorials.
Out[4]:
<matplotlib.image.AxesImage at 0x7f38fed1f520>

The PyAutoLens plot module provides methods for plotting objects and their properties, like the image of a LightProfile.

Note how, unlike the matplotlib method above, this figure is displayed with axis units of arc-seconds, a colorbar, labels, a title, etc. The PyAutoLens plot module takes care of all the heavy lifting that comes with making figures!

In [5]:
light_profile_plotter = aplt.LightProfilePlotter(
    light_profile=sersic_light_profile, grid=grid_2d
)
light_profile_plotter.figures_2d(image=True)

PyAutoLens uses MassProfile objects to represent a galaxy's mass distribution and perform ray-tracing calculations.

Below we create an elliptical isothermal MassProfile and compute its deflection angles on our Cartesian grid, where the deflection angles describe how the lens galaxy's mass bends the source's light:

In [6]:
isothermal_mass_profile = al.mp.Isothermal(
    centre=(0.0, 0.0), ell_comps=(0.1, 0.0), einstein_radius=1.6
)

deflections = isothermal_mass_profile.deflections_yx_2d_from(grid=grid_2d)

The deflection angles are easily plotted using the PyAutoLens plot module.

In [7]:
mass_profile_plotter = aplt.MassProfilePlotter(
    mass_profile=isothermal_mass_profile, grid=grid_2d
)
mass_profile_plotter.figures_2d(
    deflections_y=True, deflections_x=True
)

Many other lensing quantities are easily plotted with the MassProfilePltoter.

For readers unsure what the plotted quantities below mean, they are described in detail in chapter 1 of the HowToLens Jupyter notebook lectures:

In [8]:
mass_profile_plotter.figures_2d(
    convergence=True,
    magnification=True,
)

In PyAutoLens, a Galaxy object is a collection of LightProfile and MassProfile objects at a given redshift.

The code below creates two galaxies representing the lens and source galaxies shown in the strong lensing diagram above.

In [9]:
lens_galaxy = al.Galaxy(
    redshift=0.5, light=sersic_light_profile, mass=isothermal_mass_profile
)

source_light_profile = al.lp.Exponential(
    centre=(0.3, 0.2),
    ell_comps=(0.1, 0.0),
    intensity=0.1,
    effective_radius=0.5
)

source_galaxy = al.Galaxy(redshift=1.0, light=source_light_profile)

We can use a GalaxyPlotter to plot the properties of the lens and source galaxies.

In [10]:
lens_galaxy_plotter = aplt.GalaxyPlotter(galaxy=lens_galaxy, grid=grid_2d)
lens_galaxy_plotter.figures_2d(image=True, deflections_y=True, deflections_x=True)

source_galaxy_plotter = aplt.GalaxyPlotter(galaxy=source_galaxy, grid=grid_2d)
source_galaxy_plotter.figures_2d(image=True)

The geometry of the strong lens system depends on the cosmological distances between the Earth, the lens galaxy and the source galaxy. It therefore depends on the redshifts of the Galaxy objects.

By passing these Galaxy objects to a Tracer, PyAutoLens uses these galaxy redshifts and a cosmological model to create the appropriate strong lens system.

In [11]:
tracer = al.Tracer.from_galaxies(
    galaxies=[lens_galaxy, source_galaxy], cosmology=al.cosmo.Planck15()
)

image_2d = tracer.image_2d_from(grid=grid_2d)

tracer_plotter = aplt.TracerPlotter(tracer=tracer, grid=grid_2d)
tracer_plotter.figures_2d(image=True)

The TracerPlotter includes the MassProfile quantities we plotted previously, which can be plotted as a subplot that plots all these quantities simultaneously.

All of the objects introduced above are extensible. Galaxy objects can take many Profile's and Tracer's many Galaxy's.

If the galaxies are at different redshifts a strong lensing system with multiple lens planes will be created, performing complex multi-plane ray-tracing calculations.

To finish, lets create a Tracer with 3 galaxies at 3 different redshifts, forming a system with two distinct Einstein rings! The mass distribution of the first galaxy has separate components for its stellar mass and dark matter, where the stellar components use a LightAndMassProfile via the lmp module.

In [12]:
lens_galaxy_0 = al.Galaxy(
    redshift=0.5,
    bulge=al.lmp.Sersic(
        centre=(0.0, 0.0),
        ell_comps=(0.0, 0.05),
        intensity=0.5,
        effective_radius=0.3,
        sersic_index=3.5,
        mass_to_light_ratio=0.6,
    ),
    disk=al.lmp.Exponential(
        centre=(0.0, 0.0),
        ell_comps=(0.0, 0.1),
        intensity=1.0,
        effective_radius=2.0,
        mass_to_light_ratio=0.2,
    ),
    dark=al.mp.NFWSph(centre=(0.0, 0.0), kappa_s=0.08, scale_radius=30.0),
)

lens_galaxy_1 = al.Galaxy(
    redshift=1.0,
    bulge=al.lp.Exponential(
        centre=(0.00, 0.00),
        ell_comps=(0.05, 0.05),
        intensity=1.2,
        effective_radius=0.1,
    ),
    mass=al.mp.Isothermal(
        centre=(0.0, 0.0), ell_comps=(0.05, 0.05), einstein_radius=0.6
    ),
)

source_galaxy = al.Galaxy(
    redshift=2.0,
    bulge=al.lp.Sersic(
        centre=(0.0, 0.0),
        ell_comps=(0.0, 0.111111),
        intensity=0.7,
        effective_radius=0.1,
        sersic_index=1.5,
    ),
)

tracer = al.Tracer.from_galaxies(galaxies=[lens_galaxy_0, lens_galaxy_1, source_galaxy])

tracer_plotter = aplt.TracerPlotter(tracer=tracer, grid=grid_2d)
tracer_plotter.figures_2d(image=True)

Chapter 1 of the HowToLens Jupyter notebook tutorials gives new users a full run through of lensing calculations with PyAutoLens and introduces those not familiar with gravitational lensing to the core theory.

A more detailed overview of these calculations can also be found in notebooks/overview/overview_1_lensing.ipynb.

Lens Modeling

Lens modeling is the process of taking data of a strong lens (e.g. imaging data from the Hubble Space Telescope or interferometer data from ALMA) and fitting it with a lens model, to determine the LightProfile's and MassProfile's that best represent the observed strong lens.

Lens modeling with PyAutoLens uses the probabilistic programming language PyAutoFit, an open-source project that allows complex model fitting techniques to be straightforwardly integrated into scientific modeling software. Check it out if you are interested in developing your own software to perform advanced model-fitting!

We import PyAutoFit separately to PyAutoLens

In [13]:
import autofit as af

In this example, we consider Hubble Space Telescope imaging of the strong lens SLACS2303+1422. First, lets load this imaging dataset and plot it.

Note that the luminous emission of the foreground lens galaxy has been pre-subtracted in this image, to better highlight the source.

In [14]:
dataset_path = path.join("dataset", "slacs", "slacs2303+1422")

imaging = al.Imaging.from_fits(
    image_path=path.join(dataset_path, "image.fits"),
    psf_path=path.join(dataset_path, "psf.fits"),
    noise_map_path=path.join(dataset_path, "noise_map.fits"),
    pixel_scales=0.05,
)

imaging_plotter = aplt.ImagingPlotter(imaging=imaging)
imaging_plotter.figures_2d(image=True)

We next mask the dataset, to remove the exterior regions of the image that do not contain emission from the lensed source galaxy.

Note how when we plot the Imaging below, the figure now zooms into the masked region.

In [15]:
mask = al.Mask2D.circular(
    shape_native=imaging.shape_native,
    pixel_scales=imaging.pixel_scales,
    radius=3.0
)

imaging = imaging.apply_mask(mask=mask)

imaging_plotter = aplt.ImagingPlotter(imaging=imaging)
imaging_plotter.figures_2d(image=True)
2022-02-13 14:02:39,127 - autoarray.dataset.imaging - INFO - IMAGING - Data masked, contains a total of 11281 image-pixels

We now compose the lens model that we fit to the data using Model objects. This behave analogously to the Galaxy, LightProfile and MassProfile objects above, however their parameters are not specified and are instead determined by a fitting procedure.

We will fit our strong lens data with two galaxies:

  • A lens galaxy with an Isothermal MassProfile representing its mass.
  • A source galaxy with an Exponential LightProfile representing a disk.

The redshifts of the lens (z=0.5) and source(z=1.0) are fixed.

In [16]:
lens_galaxy_model = af.Model(
    al.Galaxy,
    redshift=0.5,
    mass=al.mp.Isothermal
)

source_galaxy_model = af.Model(al.Galaxy, redshift=1.0, disk=al.lp.Exponential)

We combine the lens and source model galaxies above into a Collection, which is the model we will fit. Note how we could easily extend this object to compose highly complex models containing many galaxies.

In [17]:
model = af.Collection(galaxies=af.Collection(lens=lens_galaxy_model, source=source_galaxy_model))

By printing the Model's we see that each parameters has a prior associated with it, which is used by the model-fitting procedure to fit the model.

In [18]:
print(lens_galaxy_model)
print()
print(source_galaxy_model)
Galaxy (centre_0, GaussianPrior, mean = 0.0, sigma = 0.1), (centre_1, GaussianPrior, mean = 0.0, sigma = 0.1), (ell_comps_0, GaussianPrior, mean = 0.0, sigma = 0.3), (ell_comps_1, GaussianPrior, mean = 0.0, sigma = 0.3), (einstein_radius, UniformPrior, lower_limit = 0.0, upper_limit = 8.0)

Galaxy (centre_0, GaussianPrior, mean = 0.0, sigma = 0.3), (centre_1, GaussianPrior, mean = 0.0, sigma = 0.3), (ell_comps_0, GaussianPrior, mean = 0.0, sigma = 0.5), (ell_comps_1, GaussianPrior, mean = 0.0, sigma = 0.5), (intensity, LogUniformPrior, lower_limit = 1e-06, upper_limit = 1000000.0), (effective_radius, UniformPrior, lower_limit = 0.0, upper_limit = 30.0)

The info attribute shows the model information in a more readable format:

In [ ]:
print(model.info)

We now choose the 'non-linear search', which is the fitting method used to determine the set of LightProfile and MassProfile parameters that best-fit our data.

In this example we use dynesty, a nested sampling algorithm that in our experience has proven very effective at lens modeling.

In [19]:
search = af.DynestyStatic(name="overview_example")
2022-02-13 14:02:39,604 - autofit.non_linear.abstract_search - INFO - Creating search

To perform the model-fit, we create an AnalysisImaging object which contains the log likelihood function that the non-linear search calls to fit the lens model to the data.

In [20]:
analysis = al.AnalysisImaging(dataset=imaging)

To perform the model-fit we pass the model and analysis to the search's fit method. This will output results (e.g., dynesty samples, model parameters, visualization) to hard-disk.

However, the lens modeling of this system takes 10-20 minutes. Therefore, to save time in this introduction, we have commented out the fit function below so you can skip through to the next section of the notebook. Feel free to uncomment the code and run the lens modeling yourself!

Once a model-fit is running, PyAutoLens outputs the results of the search to hard-disk on-the-fly. This includes lens model parameter estimates with errors non-linear samples and the visualization of the best-fit lens model inferred by the search so far.

In [21]:
# result = search.fit(model=model, analysis=analysis)

The animation below shows a slide-show of the lens modeling procedure. Many lens models are fitted to the data over and over, gradually improving the quality of the fit to the data and looking more and more like the observed image.

Lens Modeling Animation

Credit: Amy Etherington

The fit returns a Result object, which contains the best-fit Tracer and FitImaging as well as the full posterior information of the non-linear search, including all parameter samples, log likelihood values and tools to compute the errors on the lens model.

Chapter 2 of the HowToLens Jupyter notebook tutorials gives new users a full run through of lens modeling, including a primer on Bayesian statistics, how a non-linear search works and effective strategies for fitting data of a strong lens.

More examples can be found in the notebooks/overview/modeling.py and notebooks/imaging/modeling packages.

Simulating Lenses

PyAutoLens provides tool for simulating strong lens data-sets, which can be used to test lens modeling pipelines and train neural networks to recognise and analyse images of strong lenses.

Simulating strong lens images uses a SimulatorImaging object, which models the process that an instrument like the Hubble Space Telescope goes through observe a strong lens. This includes accounting for the exposure time to determine the signal-to-noise of the data, accounting for diffraction of the observed light by the telescope optics and accounting for the background sky in the exposure which acts as a source of noise.

In [22]:
psf_2d = al.Kernel2D.from_gaussian(shape_native=(11, 11), sigma=0.1, pixel_scales=0.05)

simulator = al.SimulatorImaging(
    exposure_time=300.0, background_sky_level=1.0, psf=psf_2d, add_poisson_noise=True
)

Once we have a simulator, we can use it to create an imaging dataset which consists of an image, noise-map and Point Spread Function (PSF) by passing it a tracer and grid.

This uses the tracer above to create the image of the strong lens and then add the effects that occur during data acquisition.

In [23]:
imaging = simulator.via_tracer_from(tracer=tracer, grid=grid_2d)

By plotting a subplot of the Imaging dataset, we can see this object includes the observed image of the strong lens (which has had noise and other instrumental effects added to it) as well as a noise-map and PSF:

In [24]:
imaging_plotter = aplt.ImagingPlotter(imaging=imaging)
imaging_plotter.subplot_imaging()