Tutorial : Fitting a line profile using deep learning
In this tutorial you will, hopefully, learn how fit Lyman-alpha line profiles using deep learning with zELDA.
Getting started
Let’s start by loading zELDA creating a mock line profile that we will fit later. For more details on how to create a mock line profile go to Mock line profiles
import Lya_zelda_II as Lya
your_grids_location = '/This/Folder/Contains/The/Grids/'
Lya.funcs.Data_location = your_grids_location
Geometry = 'Thin_Shell_Cont'
LyaRT_Grid = Lya.load_Grid_Line( Geometry )
# Defining the model parameters:
z_t = 0.5 # redshift of the source
V_t = 50.0 # Outflow expansion velocity [km/s]
log_N_t = 20. # Logarithmic of the neutral hydrogen column density [cm**-2]
t_t = 0.01 # Dust optical depth
log_EW_t = 1.5 # Logarithmic the intrinsic equivalent width [A]
W_t = 0.5 # Intrinsic width of the line [A]
F_t = 1. # Total flux of the line
# Defining the quality of the line profile:
PNR_t = 15.0 # Signal to noise ratio of the maximum of the line.
FWHM_t = 0.2 # Full width half maximum diluting the line. Mimics finite resolution. [A]
PIX_t = 0.1 # Wavelength binning of the line. [A]
w_Arr , f_Arr , s_Arr = Lya.Generate_a_real_line( z_t , V_t, log_N_t, t_t, F_t, log_EW_t, W_t , PNR_t, FWHM_t, PIX_t, LyaRT_Grid, Geometry )
where /This/Folder/Contains/The/Grids/ is the place where you store the LyaRT data grids, as shown in the installation section. And… It’s done! w_Arr is a numpy array that contains the wavelength where the line profile is evaluated. Meanwhile, f_Arr is the actual line profile. s_Arr is the uncertainty of the flux density. Remeber that if you want to use the line profile grid with lower RAM memory occupation you must pass MODE=’LIGHT’ to Lya.load_Grid_Line.
Let’s have a look to how the line looks:
w_Arr , f_Arr , s_Arr = Lya.Generate_a_real_line( z_t , V_t, log_N_t, t_t, F_t, log_EW_t, W_t , PNR_t, FWHM_t, PIX_t, LyaRT_Grid, Geometry )
w_pix_Arr , f_pix_Arr = Lya.plot_a_rebinned_line( w_Arr , f_Arr , PIX_t )
import pylab as plt
plt.plot( w_pix_Arr , f_pix_Arr )
plt.xlabel('wavelength[A]' , size=15 )
plt.ylabel('Flux density [a.u.]' , size=15 )
plt.xlim(1815,1835)
plt.show()
Now that we have our mock line profile. Let’s load the neural network. As we have produce a line profile for an outflow (V_t>0) we are going to load the deep neural network for outflows
machine_data = Lya.Load_NN_model( 'Outflow' )
In case you want to do the analysis for inflows just call Lya.Load_NN_model( ‘Inflow’ ). machine_data is a python dictionary that contains all the necessary data for the deep neural network approach. Let’s pick up from it two variables:
machine = machine_data['Machine' ]
w_rest_Arr = machine_data[ 'w_rest' ]
machine is an object from skitlearn with the trained neural network and w_rest_Arr is the rest frame wavelength where the line profiles used for the training were evaluated. w_rest_Arr is important to check that the neural networks is working in the same wavelength array that the line profiles will be evaluated. In principle you don’t have to do anything with w_rest_Arr, but we need to pass it to other functions.
Using the DNN in the un-perturbed line profile
Let’s start by simple evaluating the DNN using the mock line profile without perturbing it:
Sol , z_sol = Lya.NN_measure( w_Arr , f_Arr , s_Arr , FWHM_t , PIX_t , machine , w_rest_Arr , N_iter=None )
Done! . Sol is a matrix that contains the prediction by the DNN and z_sol is the predicted redshift. You can print the predicted properties doing:
print( 'The measured redshift is' , z_sol )
print( 'The measured logarithm of the expasion velocity is' , Sol[0,1] )
print( 'The measured logarithm of the HI column density is' , Sol[0,2] )
print( 'The measured logarithm of the dust optical depth is' , Sol[0,3] )
print( 'The measured logarithm of the intrinsic equivalent width is' , Sol[0,4] )
print( 'The measured logarithm of the intrinsic width is' , Sol[0,5] )
print( 'The measured shift of the true Lya wavelgnth from the maximum of the line is' , Sol[0,0] )
This should give something like
The measured redshift is 0.49994403239322693
The measured logarithm of the expasion velocity is 1.5821419036064905
The measured logarithm of the HI column density is 20.149247231711733
The measured logarithm of the dust optical depth is -3.310662004999448
The measured logarithm of the intrinsic equivalent width is 1.458352960574508
The measured logarithm of the intrinsic width is -0.804093047888869
The measured shift of the true Lya wavelgnth from the maximum of the line is -1.2773994188976223
Let’s see how this new spectrum compares with the target:
PNR = 100000. # let's put infinite signal to noise in the model line
V_sol = 10**Sol[0,1] # Expansion velocity km/s
logN_sol = Sol[0,2] # log of HI column density cm**-2
t_sol = 10**Sol[0,3] # dust optical depth
logE_sol = Sol[0,4] # log intrinsic EW [A]
W_sol = 10**Sol[0,5] # intrinsic width [A]
# creates the line
w_One_Arr , f_One_Arr , _ = Lya.Generate_a_real_line( z_sol , V_sol, logN_sol, t_sol, F_t, logE_sol, W_sol, PNR, FWHM_t, PIX_t, LyaRT_Grid, Geometry )
# plot the target and the predicted line
w_pix_One_Arr , f_pix_One_Arr = Lya.plot_a_rebinned_line( w_One_Arr , f_One_Arr , PIX_t )
plt.plot( w_pix_Arr , f_pix_Arr , label='Target' )
plt.plot( w_pix_One_Arr , f_pix_One_Arr , label='1 iter' )
plt.legend(loc=0)
plt.xlabel('wavelength[A]' , size=15 )
plt.ylabel('Flux density [a.u.]' , size=15 )
plt.xlim(1815,1835)
plt.show()
You should get something like:
Using the DNN with Monte Carlo perturbations
Normally, it is better to do more than one iteration, as it leads to better results. These iterations basically perturb the flux density f_Arr by adding gaussian noise with the amplitude of s_Arr in each wavelength bin. Then, this new perturbed spectrum is send to the DNN. For each of these iterations the output of the DNN is stored. For example for 1000 iterations :
>>> Sol , z_sol , log_V_Arr , log_N_Arr , log_t_Arr , z_Arr , log_E_Arr , log_W_Arr = Lya.NN_measure( w_Arr , f_Arr , s_Arr , FWHM_t , PIX_t , machine , w_rest_Arr , N_iter=1000 )
The arrays log_V_Arr, log_N_Arr, log_t_Arr, z_Arr, log_E_Arr and log_W_Arr contain the output of the DNN for the iterations for the logarithms of the expansion velocity, the logarithm of the neutral hydrogen column density, the logarithm of the dust optical depth, the redshift, the logarithm of the intrinsic equivalent width and the logarithm of the intrinsic width of the line. From these arrays we can compute the result from the DNN analysis by taking the 50th percentile. The +-1 sigma uncertainty can be computed as the 16th and 84th percentile.
import numpy as np
# Redshitft
z_50 = np.percentile( z_Arr , 50 )
z_16 = np.percentile( z_Arr , 16 )
z_84 = np.percentile( z_Arr , 84 )
# Expansion velocity
V_50 = 10 ** np.percentile( log_V_Arr , 50 )
V_16 = 10 ** np.percentile( log_V_Arr , 16 )
V_84 = 10 ** np.percentile( log_V_Arr , 84 )
# Logarithmic of HI column density
log_N_50 = np.percentile( log_N_Arr , 50 )
log_N_16 = np.percentile( log_N_Arr , 16 )
log_N_84 = np.percentile( log_N_Arr , 84 )
# Dust optical depth
t_50 = 10 ** np.percentile( log_t_Arr , 50 )
t_16 = 10 ** np.percentile( log_t_Arr , 16 )
t_84 = 10 ** np.percentile( log_t_Arr , 84 )
# Logarithmic of intrinsic equivalent width
log_E_50 = np.percentile( log_E_Arr , 50 )
log_E_16 = np.percentile( log_E_Arr , 16 )
log_E_84 = np.percentile( log_E_Arr , 84 )
# Intrinsic width
W_50 = 10 ** np.percentile( log_W_Arr , 50 )
W_16 = 10 ** np.percentile( log_W_Arr , 16 )
W_84 = 10 ** np.percentile( log_W_Arr , 84 )
let’s see how the line profiles look:
# Compute the 100 iterations line profile
w_50th_Arr , f_50th_Arr , _ = Lya.Generate_a_real_line( z_50 , V_50, log_N_50, t_50, F_t, log_E_50, W_50, PNR, FWHM_t, PIX_t, LyaRT_Grid, Geometry )
# Get cooler profiles
w_pix_50th_Arr , f_pix_50th_Arr = Lya.plot_a_rebinned_line( w_50th_Arr , f_50th_Arr , PIX_t )
# Plot
plt.plot( w_pix_Arr , f_pix_Arr , label='Target' )
plt.plot( w_pix_One_Arr , f_pix_One_Arr , label='1 iter' )
plt.plot( w_pix_50th_Arr , f_pix_50th_Arr , label='1000 iter')
plt.legend(loc=0)
plt.xlabel('wavelength[A]' , size=15 )
plt.ylabel('Flux density [a.u.]' , size=15 )
plt.xlim(1815,1835)
plt.show()
finally, let’s compare the parameters that we got with the input:
print( 'The true redshift is' , z_t , 'and the predicted is' , z_50 , '(-' , z_50-z_16 , ', +' , z_84-z_50 , ')' )
print( 'The true expansion velocity is' , V_t , 'and the predicted is' , V_50 , '(-' , V_50-V_16 , ', +' , V_84-V_50 , ')' )
print( 'The true dust optical depth is' , t_t , 'and the predicted is' , t_50 , '(-' , t_50-t_16 , ', +' , t_84-t_50 , ')' )
print( 'The true intrinsic width is' , W_t , 'and the predicted is' , W_50 , '(-' , W_50-W_16 , ', +' , W_84-W_50 , ')' )
print( 'The true log of HI column density is' , log_N_t , 'and the predicted is' , log_N_50 , '(-' , log_N_50-log_N_16 , ', +' , log_N_84-log_N_50 , ')' )
print( 'The true log of equivalent width is' , log_EW_t , 'and the predicted is' , log_E_50 , '(-' , log_E_50-log_E_16 , ', +' , log_E_84-log_E_50 , ')' )
This should give something like:
The true redshift is 0.5 and the predicted is 0.49999833428137275 (- 0.00017321665235831007 , + 0.0003615214512187048 )
The true expansion velocity is 50.0 and the predicted is 47.070589157142614 (- 16.100374040796254 , + 48.27234502291723 )
The true dust optical depth is 0.01 and the predicted is 0.00379679848371737 (- 0.003483235501588427 , + 0.049396128990436335 )
The true intrinsic width is 0.5 and the predicted is 0.280484205908298 (- 0.12228181625600373 , + 0.2150273326940031 )
The true log of HI column density is 20.0 and the predicted is 20.019139948537997 (- 0.5728866241916535 , + 0.207985045834004 )
The true log of equivalent width is 1.5 and the predicted is 1.5595962407058306 (- 0.09992888862396399 , + 0.16009784914990055 )
The particular values that you print will be slightly different when you run it, but more or less it should go in the same direction.
That was fun, hah? Now you know how to use the deep neural network scheme in zELDA.