Example EasyLens Worksheet

Welcome to EasyLens! This notebook will take you through all the steps in modeling a strong lens system using in multi-band survey data.

1. Set up your workbench

First we need to make sure we have all the software we need, and also get the data that we will be working on. The first step is to pip install all the packages in requirements.txt, like this:

pip install -r requirements.txt

Next, you'll need to do

python setup.py install

from the top level EasyLens directory to get all the EasyLens code onto your PYTHONPATH.

Once this is all done, you can run this notebook, starting with the following import cell.



In [ ]:

    
# External modules - try "pip install <module>" if you get an error.
import astropy.io.fits as pyfits
import astropy.wcs as pywcs
import pickle
import numpy as np
import os
import easylens

# It'll make the notebook clearer if we get a few tools out and give them easy names:
from easylens.Data.lens_system import LensSystem
from easylens.Data.show_lens import ShowLens
from easylens.Data.exposure import Exposure
import easylens.util as util

# We'll be doing some plotting:
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.axes_grid1 import AxesGrid, make_axes_locatable

2. Which lens are we going to model?

In this example, we'll use the system whose cutout images have been checked into the data directory. In future, we'll provide code to download any set of images you want to model. In either case we need to specify the lens name, and the bandpasses we'll be using.



In [ ]:

    
LENS = 'DESJ033056.9318-522813.0188'
FILTERS = ['g','r','i','z']



In [ ]:

    
# When the data is local, we just need to specify the path to the data folders. 
DATA_DIR = '../data/'+LENS



In [ ]:

    
# Now we can set up our file names, stored in a dictionary to enable reference by filter:
SCIFILE = {}
WHTFILE = {}
PSFFILE = {}
for filter in FILTERS:
    SCIFILE[filter] = DATA_DIR+'/'+LENS+'_'+filter+'_sci.fits'
    WHTFILE[filter] = DATA_DIR+'/'+LENS+'_'+filter+'_wht.fits'
    PSFFILE[filter] = DATA_DIR+'/'+LENS+'_'+filter+'_psf.fits'

3. Make a lens system object

Here's where we start making our model. We'll need to give it some images, and tell it where on the sky it is.



In [ ]:

    
# Get the sky position from one of the science images:
ra_center, dec_center = util.get_ra_dec_center(SCIFILE['g'])

# Initialize the model:
lensSystem = LensSystem(name="object03", ra=ra_center, dec=dec_center)

The lens system now needs to be given the images we have of it. Each cutout image needs to be accompanied by a 'weight' (wht) and PSF (psf) image; the Exposure class stores all three of these images together.



In [ ]:

    
# Let's keep the images in a dictionary, so we can refer to them by filter.
image = {}

for filter in FILTERS:

    # Set up an image object for this filter:
    image[filter] = Exposure(lensSystem.ra, lensSystem.dec)    
    image[filter].load_all(SCIFILE[filter],PSFFILE[filter],WHTFILE[filter])

    # Add this image to the model:
    lensSystem.add_image_data(image[filter], filter+'_band')

show images



In [ ]:

    
showLens = ShowLens(lensSystem)
print "Cutout images"
f, axes = showLens.show_images()
plt.show()

print "PSF images"
f, axes = showLens.show_psf()
plt.show()

print "Looking for the mask center locations"
%matplotlib notebook
f = showLens.show_single(2) #use number of band as input value: 1 for g_band, 2 for r_band etc.

assign masks for different objects

Define the center of the mask according the pixel you find in the picture above. Feel free to use the zoom functions etc.



In [ ]:

    
%matplotlib inline
xpix_1 = 50
ypix_1 = 49

xpix_2 = 77
ypix_2 = 45

ra_pos_1, dec_pos_1 = lensSystem.get_angle_coord('r_band',50,49)
ra_pos_2, dec_pos_2 = lensSystem.get_angle_coord('r_band',74,45)

print('Source 1')
print('RA/DEC=',ra_pos_1,'/',dec_pos_1)
print('Source 2')
print('RA/DEC=',ra_pos_2,'/',dec_pos_2)

kwargs_mask = {"type": "circular", "ra_c": 2, "dec_c": 1, "width": 10, "radius": 12}
kwargs_mask_source1 = {"type": "circular", "ra_c": 0, "dec_c": 0, "width": 10, "radius": 4}
kwargs_mask_source2 = {"type": "circular", "ra_c": ra_pos_2, "dec_c": dec_pos_2, "width": 10, "radius": 2}

Plot pictures with newly defined masks



In [ ]:

    
f, axes = showLens.show_images(ra_pos_2, dec_pos_2, kwargs_mask)
plt.show()

f, axes = showLens.show_images(ra_pos_2, dec_pos_2, kwargs_mask_source1)
plt.show()

f, axes = showLens.show_images(ra_pos_2, dec_pos_2, kwargs_mask_source2)
plt.show()

f, axes = showLens.show_pixel_hist()
plt.show()



In [ ]:

    
from easylens.easylens import EasyLens
easyLens = EasyLens(lensSystem, frame_bool = {"g_band": True, "r_band": True, "i_band": True,
            "Y_band": False, "z_band": True}, subgrid_res=2, lens_type='SPEMD')

w_SED1 = lensSystem.get_sed_estimate(ra_pos_1, dec_pos_1)
w_SED2 = lensSystem.get_sed_estimate(ra_pos_2, dec_pos_2)
print w_SED1, w_SED2

Now choose one particular band, for which we will show some reference pictures.



In [ ]:

    
refb = 'g_band'

SED modeling of source 1 (primary the lensing galaxy)



In [ ]:

    
from easylens.easylens import Source
#ra_pos_1 = 0
#dec_pos_1 = 0
beta_1 = 1.5
n_max_1 = 20

#kwargs_mask_test = {"type": "square", "ra_c": 2, "dec_c": 1, "width": 20, "radius": 12}

source1 = Source(name="source1", ra_pos=ra_pos_1, dec_pos=dec_pos_1, beta=beta_1, n_max=n_max_1, w_SED=w_SED1, lens_bool=False)
easyLens.add_source(source1, over_write=True)
easyLens.del_source("source2")
easyLens.update_mask(kwargs_mask_source1)
numPix = easyLens.get_pixels_unmasked()
d = easyLens.get_data_vector()
C_D_inv = easyLens.get_C_D_inv_vector()

for i in range(2):
    A = easyLens.get_response()
    param_array, model_array = easyLens.get_inverted(A, C_D_inv, d)

    data_list = easyLens.get_data_list()
    model_list = easyLens.get_model_images(model_array)
    residual_list = easyLens.get_residuals(model_array)
    de_convolved_list = easyLens.get_deconvolved(param_array)
    source_list_ref = easyLens.get_sources_original(param_array, refb)    
    
    A_sed = easyLens.get_response_sed(param_array)
    param_sed_array, model_sed_array = easyLens.get_inverted(A_sed, C_D_inv, d)
    w_SED1 = {'g_band': param_sed_array[0], 'r_band': param_sed_array[1], 'i_band': param_sed_array[2], 'z_band': param_sed_array[3]}
    source1 = Source(name="source1", ra_pos=ra_pos_1, dec_pos=dec_pos_1, beta=beta_1, n_max=n_max_1, w_SED=w_SED1, lens_bool=False)
    easyLens.add_source(source1, over_write=True)



In [ ]:

    
frame_list = easyLens.frame_list
print i
chi2 = np.sum((model_array-d)**2*C_D_inv)/numPix
print chi2
# residuals
print("residual images")
f, axes = showLens.show_list(residual_list, frame_list)
plt.show()
# sources separate
print("sources")
f, axes = showLens.show_list(source_list_ref, [refb])
plt.show()

SED modeling of source 2 (primary the source galaxy)



In [ ]:

    
#ra_pos_2 = ra_pos_2         
#dec_pos_2 = dec_pos_2
beta_2 = 1.
n_max_2 = 10


source2 = Source(name="source1", ra_pos=ra_pos_2, dec_pos=dec_pos_2, beta=beta_2, n_max=n_max_2, w_SED=w_SED2, lens_bool=False)
easyLens.add_source(source2, over_write=True)
easyLens.del_source("source2")
easyLens.update_mask(kwargs_mask_source2)
numPix = easyLens.get_pixels_unmasked()
d = easyLens.get_data_vector()
C_D_inv = easyLens.get_C_D_inv_vector()

for i in range(3):
    A = easyLens.get_response()
    param_array, model_array = easyLens.get_inverted(A, C_D_inv, d)

    data_list = easyLens.get_data_list()
    model_list = easyLens.get_model_images(model_array)
    residual_list = easyLens.get_residuals(model_array)
    de_convolved_list = easyLens.get_deconvolved(param_array)
    source_list_ref = easyLens.get_sources_original(param_array, refb)    
    
    A_sed = easyLens.get_response_sed(param_array)
    param_sed_array, model_sed_array = easyLens.get_inverted(A_sed, C_D_inv, d)
    w_SED2 = {'g_band': param_sed_array[0], 'r_band': param_sed_array[1], 'i_band': param_sed_array[2], 'z_band': param_sed_array[3]}
    source2 = Source(name="source2", ra_pos=ra_pos_2, dec_pos=dec_pos_2, beta=beta_2, n_max=n_max_2, w_SED=w_SED2, lens_bool=False)
    easyLens.add_source(source2, over_write=True)



In [ ]:

    
frame_list = easyLens.frame_list
print i
chi2 = np.sum((model_array-d)**2*C_D_inv)/numPix
print chi2
# residuals
print("residual images")
f, axes = showLens.show_list(residual_list, frame_list)
plt.show()

# sources separate
print("sources")
f, axes = showLens.show_list(source_list_ref, [refb])
plt.show()

add sources together



In [ ]:

    
beta=1.3
n_max = 20
source1 = Source(name="source1", ra_pos=0, dec_pos=0, beta=beta, n_max=n_max, w_SED=w_SED1, lens_bool=False)
source2 = Source(name="source2", ra_pos=ra_pos_2, dec_pos=dec_pos_2, beta=beta, n_max=n_max, w_SED=w_SED2, lens_bool=False)
easyLens.add_source(source1, over_write=True)
easyLens.add_source(source2, over_write=True)
easyLens.update_mask(kwargs_mask)

joint fit (but unlensed)



In [ ]:

    
d = easyLens.get_data_vector()
C_D_inv = easyLens.get_C_D_inv_vector()
A = easyLens.get_response()
param_array, model_array = easyLens.get_inverted(A, C_D_inv, d)



In [ ]:

    
data_list = easyLens.get_data_list()
model_list = easyLens.get_model_images(model_array)
residual_list = easyLens.get_residuals(model_array)
de_convolved_list = easyLens.get_deconvolved(param_array)
source_list_ref = easyLens.get_sources_original(param_array, refb)
frame_list = easyLens.frame_list

numPix = easyLens.get_pixels_unmasked()
chi2 = np.sum((model_array-d)**2*C_D_inv)/numPix
print chi2

# original images
print("original images")
f, axes = showLens.show_list(data_list, frame_list)
plt.show()


# modeled images
print("modeled images")
f, axes = showLens.show_list(model_list, frame_list)
plt.show()

# residuals
print("residual images")
f, axes = showLens.show_list(residual_list, frame_list)
plt.show()

# de-convolved images
print("de-convolved images")
f, axes = showLens.show_list(de_convolved_list, frame_list)
plt.show()

# sources separate
print("sources")
f, axes = showLens.show_list(source_list_ref, [refb, refb])
plt.show()

joint and lensed fit

try to fit with the same $\beta$ and $n_{max}$ as above.



In [ ]:

    
beta = 1.3
n_max = 20
source1 = Source(name="source1", ra_pos=0, dec_pos=0, beta=beta, n_max=n_max, w_SED=w_SED1, lens_bool=False)
source2 = Source(name="source2", ra_pos=ra_pos_2, dec_pos=dec_pos_2, beta=0.6, n_max=5, w_SED=w_SED2, lens_bool=True)
easyLens.add_source(source1, over_write=True)
easyLens.add_source(source2, over_write=True)
easyLens.update_mask(kwargs_mask)
update = False
#{'e_2': 0.12185159487243634, 'e_1': 0.066152388530625786, 'center_x': -0.64365565246052192, 'center_y': 0.023685031811570757, 
#'phi_E': 5.5195718134302467, 'gamma': 1.8214123493075136}
if update is True:
    kwargs_lens = {"phi_E": np.mean(samples[:,0]), "center_x": np.mean(samples[:,1]), "center_y": np.mean(samples[:,2]), 
                   "gamma":np.mean(samples[:,3]), "e1":np.mean(samples[:,4]), "e2":np.mean(samples[:,5])}
    print kwargs_lens
#kwargs_lens = {"phi_E_sis": 4.7, "center_x_sis": -0.55, "center_y_sis": -0.3}
else:
    kwargs_lens = {"phi_E": 5.3, "center_x": 0.7, "center_y": 1.7,'gamma':1.7,'e1':0.2, 'e2':0.5}
easyLens.add_lens(kwargs_lens)



In [ ]:

    
d = easyLens.get_data_vector()
C_D_inv = easyLens.get_C_D_inv_vector()
A = easyLens.get_response()
param_array, model_array = easyLens.get_inverted(A, C_D_inv, d)

save_image = True
data_list = easyLens.get_data_list()
model_list = easyLens.get_model_images(model_array)
residual_list = easyLens.get_residuals(model_array)
de_convolved_list = easyLens.get_deconvolved(param_array)
source_list_original_ref = easyLens.get_sources_original(param_array, refb)
source_list_lensed_ref = easyLens.get_sources_lensed(param_array, refb)
source_list_image_ref = easyLens.get_sources_image(param_array, refb)
frame_list = easyLens.frame_list

numPix = easyLens.get_pixels_unmasked()
chi2 = np.sum((model_array-d)**2*C_D_inv)/numPix
print chi2
%matplotlib inline
# original images
print("original images")
f, axes = showLens.show_list(data_list, frame_list)
if save_image is True:
    plt.savefig('original_images.pdf', format='pdf')
plt.show()


# modeled images
print("modeled images")
f, axes = showLens.show_list(model_list, frame_list)
if save_image is True:
    plt.savefig('modeled_images.pdf', format='pdf')
plt.show()

# residuals
print("residual images")
f, axes = showLens.show_list(residual_list, frame_list)
if save_image is True:
    plt.savefig('residual_images.pdf', format='pdf')
plt.show()

# de-convolved images
print("de-convolved images")
f, axes = showLens.show_list(de_convolved_list, frame_list)
if save_image is True:
    plt.savefig('de_convolved_images.pdf', format='pdf')
plt.show()

# sources separate
print("sources original")
f, axes = showLens.show_list(source_list_original_ref, [refb, refb])
if save_image is True:
    plt.savefig('original_objects.pdf', format='pdf')
plt.show()


print("sources lensed")
f, axes = showLens.show_list(source_list_lensed_ref, [refb, refb])
if save_image is True:
    plt.savefig('lensed_objects.pdf', format='pdf')
plt.show()

print("sources lensed and convolved")
f, axes = showLens.show_list(source_list_image_ref, [refb, refb])
if save_image is True:
    plt.savefig('lensed_convolved_objects.pdf', format='pdf')
plt.show()



In [ ]:

    
from easylens.Fitting.mcmc import MCMC_sampler
kwargs_fixed = {}
kwargs_lens = {"phi_E": 5.3, "center_x": 0.7, "center_y": 1.7, "gamma":1.7, 'e1':0.2, "e2":0.5}
#kwargs_lens = {"phi_E": 4.7, "center_x": -0.55, "center_y": -0.3, "gamma":2., 'e1':0.1, "e2":0.2, 'e1_ext':0.01, 'e2_ext':0.01}
sampler = MCMC_sampler(easyLens, kwargs_fixed)
walkerRatio = 12
n_run = 100
n_burn = 100
mean_start = [kwargs_lens["phi_E"],  kwargs_lens["center_x"], kwargs_lens["center_y"], kwargs_lens["gamma"], kwargs_lens["e1"], kwargs_lens["e2"]]
#mean_start = [kwargs_lens["phi_E"],  kwargs_lens["center_x"], kwargs_lens["center_y"], kwargs_lens["gamma"], kwargs_lens["e1"], kwargs_lens["e2"], kwargs_lens['e1_ext'], kwargs_lens['e2_ext']]
#{"phi_E": 5.3, "center_x": 0.7, "center_y": 1.7,'gamma':1.7,'e1':0.2, 'e2':0.5}
sigma_start = [0.1, 0.1, .1,.1,.1,.1]
lowerLimit = [5., -10., -10.,1.5,-1.,-1.]
upperLimit = [6., 10., 10.,2.5,1.,1.]
samples, prob = sampler.mcmc_CH(walkerRatio, n_run, n_burn, mean_start, sigma_start, lowerLimit, upperLimit, threadCount=1, init_pos=None, mpi_monch=False)



In [ ]:

    
import corner

corner.corner(samples)
plt.savefig('spemd-12-100-100.eps')
plt.show()



In [ ]: