Scripting Sculptor 01 - Modelling the example spectrum in a script

In this notebook we will introduce on how to use Sculptor’s SpecFit and SpecModel classes to generate a Sculptor fit within a python script.

The code in this example that generates the model fit is also available as a pure python script in the examples folder: example_fit_setup_script.py

We begin by importing the modules we need for this example. In addition to the SpecFit, SpecModel, and SpecOneD classes we need pkg_resources to easily access the SDSS quasar example spectrum, which is provided in sculptor/data/example_spectra.

[1]:

# Matplotlib statements to make plots nice
# %matplotlib notebook
from IPython.display import set_matplotlib_formats
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 300
plt.rcParams['savefig.dpi'] = 300
# set_matplotlib_formats('svg')

[2]:

import pkg_resources

from sculptor import specfit as scfit
from sculptor import speconed as scspec
from sculptor import specmodel as scmod

[INFO] Import "sculptor_extensions" package: my_extension
[INFO] Import "sculptor_extensions" package: qso
[INFO] SWIRE library found.
[INFO] FeII iron template of Vestergaard & Wilkes 2001 found. If you will be using these templates in your model fit and publication, please add the citation to the original work, ADS bibcode: 2001ApJS..134....1V
[INFO] FeII iron template of Tsuzuki et al. 2006 found. If you will be using these templates in your model fit and publication, please add the citation to the original work, ADS bibcode: 2006ApJ...650...57T
[INFO] FeII iron template of Boroson & Green 1992 found. If you will be using these templates in your model fit and publication, please add the citation to the original work, ADS bibcode: 1992ApJS...80..109B

These imports automatically intialize some global variables, which allow you to use the model masks and functions defined in Sculptor and its extensions. We will take a look at how to write Sculptor extensions in a later tutorial. Sculptor comes with two extensions “my_extension”, which provides an example on how to set up your own extension, and the “qso” extension, which we will use heavily in this example. The “qso” extension also comes with additional data, e.g. FeII iron templates. If you want to use them in your work Sculptor reminds you where to find the original papers, so you can cite them appropriately.

Initializing the SpecFit object

Let us begin by importing the SDSS quasar spectrum. The quasar is at a redshift of z=3.227, which we will use for initializing the SpecFit class.

[3]:

# Initialize a new SpecOneD object and read in the quasar spectrum
spec = scspec.SpecOneD()
filename = pkg_resources.resource_filename('sculptor', 'data/example_spectra/J030341.04-002321.8_0.fits')
spec.read_sdss_fits(filename)

redshift = 3.227

# Plot the quasar spectrum for confirmation that everything worked well.
spec.plot()

In the next step we initialize the SpecFit object with the SpecOneD spectrum and also supply the redshift information.

[4]:

# Initialize SpecFit object
fit = scfit.SpecFit(spec, redshift)

The SpecFit object currently does not include any models and if we use the .plot() method, it will only show us the quasar spectrum. However, we can check the rest-frame dispersion axis, if the redshift keyword was used correctly.

[5]:

fit.plot()

Adding our first SpecModel object - Fitting the continuum

In the next step we want to add a simple continuum model to our fit. However, we currently don’t really know which models and masks are available to us based on Sculptor and its loaded extensions. Let’s get a list of all model names and mask names. This information is available as a global variable in the sculptor.specmodel module:

[6]:

print('MODEL FUNCTIONS:')
for model_func in scmod.model_func_list:
    print (model_func)
print('\n')
print('MASKS:')
for mask in scmod.mask_presets:
    print(mask)

MODEL FUNCTIONS:
Constant (amp)
Power Law (amp, slope)
Gaussian (amp, cen, sigma, shift)
Lorentzian (amp, cen, gamma, shift)
My Model
Power Law (2500A)
Power Law (2500A) + BC
Power Law (2500A) + BC (fractional)
Line model Gaussian
SiIV (2G components)
CIV (2G components)
MgII (2G components)
HBeta (2G components)
HAlpha (2G components)
[OIII] doublet (2G)
[NII] doublet (2G)
[SII] doublet (2G)
CIII] complex (3G components)
SWIRE Ell2
SWIRE NGC6090
FeII template 1200-2200 (VW01, cont)
FeII template 1200-2200 (VW01, split)
FeII template 2200-3500 (VW01, cont)
FeII template 2200-3500 (VW01, split)
FeII template 2200-3500 (T06, cont)
FeII template 2200-3500 (T06, split)
FeII template 3700-7480 (BG92, cont)
FeII template 3700-5600 (BG92, split)


MASKS:
My mask
QSO Continuum+FeII
QSO Cont.W. VP06
QSO Fe+Cont.W. CIV Shen11
QSO Fe+Cont.W. MgII Shen11
QSO Fe+Cont.W. HBeta Shen11
QSO Fe+Cont.W. HAlpha Shen11

Using these models when writing Sculptor scripts usually requires to fully understand their parameters and functionality. Therefore, it is adviseable to study their source code, before attempting to write scripts.

In order to fit a model to the spectrum, we need to first add a SpecModel object to our fit (SpecFit object). Then we can access this SpecModel and add model functions to it. In order to fit the SpecModel

[7]:

# Add the continuum SpecModel
fit.add_specmodel()
# We access the initialized SpecModel object by using the first item in the SpecFit.specmodels list
contmodel = fit.specmodels[0]
# Let us rename this SpecModel 'Continuum'
contmodel.name = 'Continuum'

# Define the model function name
model_name = 'Power Law (2500A)'
# Define the model prefix
# It is important to keep track of the prefix to later access this model in the analysis
model_prefix = 'PL_'
# Add the model function to the SpecModel
contmodel.add_model(model_name, model_prefix)

We have now added the model function to the SpecModel ‘Continuum’. However, we cannot yet fit the model successfully as we have not defined the regions to which the continuum model should be fit. Let us add them manually based. The exact fit regions for the continuum model of this quasar have been determined beforehand.

[8]:

contmodel.add_wavelength_range_to_fit_mask(8300, 8620)
contmodel.add_wavelength_range_to_fit_mask(7105, 7205)
contmodel.add_wavelength_range_to_fit_mask(5400, 5450)
contmodel.add_wavelength_range_to_fit_mask(5685, 5750)
contmodel.add_wavelength_range_to_fit_mask(6145, 6212)

[INFO] Manual mask range 8300 8620
[INFO] Manual mask range 7105 7205
[INFO] Manual mask range 5400 5450
[INFO] Manual mask range 5685 5750
[INFO] Manual mask range 6145 6212

Now that the fit mask has been updated we can fit the continuum model and plot the resulting fit. We can either use the contmodel.plot() function to only show the SpecModel components or use the fit.plot() function to show all SpecModels fits.

[9]:

# Fit the continuum model
contmodel.fit()
# Plot the fitted continuum model and the spectrum
contmodel.plot()
# Plot the SpecFit fit with all SpecModels (We only have 1 at the moment)
fit.plot()

Adding and manipulating a SpecModel - Fitting the SiIV line

In our next step we will add an emission line model to fit the SiIV emission line at a rest-frame wavelength of ~1400A. To do this we start by adding another SpecModel to our fit and specify the wavelength regions (in observed-frame) we want to use for the fit.

[10]:

# Add the SiIV emission line model
fit.add_specmodel()
# Access the SpecModel object by choosing the second SpecModel object in the SpecFit
# specmodels list.
siiv_model = fit.specmodels[1]
# Rename the SpecModel
siiv_model.name = 'SiIV_line'

# Add wavelength regions to the SpecModel for the fit
siiv_model.add_wavelength_range_to_fit_mask(5790, 5870)
siiv_model.add_wavelength_range_to_fit_mask(5910, 6015)

[INFO] Manual mask range 5790 5870
[INFO] Manual mask range 5910 6015

However, we will use the pre-defined ‘SiIV (2G components)’ model function from the Sculptor qso extension, which has pre-defined model prefixes. Therefore, we pass a None as the model prefix here.

The model functions allow to pass additional keyword arguments that modify the redshift or amplitude of the emission line model. The redshift is automatically passed by the SpecModel class if we don’t specify a different value here. In this case we only want to pass an amplitude.

[11]:

model_name = 'SiIV (2G components)'
model_prefix = None
siiv_model.add_model(model_name, model_prefix, amplitude=20)

The ‘SiIV (2G components)’ model function added two emission line models with the prefixes ‘SiIV_A_’ and ‘SiIV_B_’ to our SpecModel object. We can check whether the models were successfully added to the SpecModel object by accessing the SpecModel ‘model_list’. The items in the model list are LMFIT Model objects.

[12]:

for model in siiv_model.model_list:
    print(model)

<lmfit.Model: Model(line_model_gaussian, prefix='SiIV_A_')>
<lmfit.Model: Model(line_model_gaussian, prefix='SiIV_B_')>

We see that two model functions of the type ‘line_model_gaussian’ have been added to the SiIV SpecModel. By default the model function holds the redshift parameter fixed during the fit. However, for our purposes we want the redshift to be a variable. To change this we need to access the parameters of the model function.

The parameters for each of the models in the ‘model_list’ are stored in the ‘params_list’. Each item in the ‘params_list’ is a LMFIT Parameters object, holding the individual parameters of the associated model from the ‘model_list’.

[13]:

for params in siiv_model.params_list:
    print(params)
    print('\n')

Parameters([('SiIV_A_z', <Parameter 'SiIV_A_z', value=3.227 (fixed), bounds=[3.0656499999999998:3.38835]>), ('SiIV_A_flux', <Parameter 'SiIV_A_flux', value=106446.7019431226, bounds=[0:106446701.94312261]>), ('SiIV_A_cen', <Parameter 'SiIV_A_cen', value=1399.8 (fixed), bounds=[-inf:inf]>), ('SiIV_A_fwhm_km_s', <Parameter 'SiIV_A_fwhm_km_s', value=2500, bounds=[300:20000]>)])


Parameters([('SiIV_B_z', <Parameter 'SiIV_B_z', value=3.227 (fixed), bounds=[3.0656499999999998:3.38835]>), ('SiIV_B_flux', <Parameter 'SiIV_B_flux', value=106446.7019431226, bounds=[0:106446701.94312261]>), ('SiIV_B_cen', <Parameter 'SiIV_B_cen', value=1399.8 (fixed), bounds=[-inf:inf]>), ('SiIV_B_fwhm_km_s', <Parameter 'SiIV_B_fwhm_km_s', value=2500, bounds=[300:20000]>)])

To access a single LMFIT Parameter object from the Parameters we need to use its name, which consists of the model prefix and the parameter name, which are shown above.

By accessing a specific parameters we can change its attributes, defined in the LMFIT documentation: * name (str) – Name of the Parameter. * value (float, optional) – Numerical Parameter value * vary (bool, optional) – Whether the Parameter is varied during a fit (default is True). * min (float, optional) – Lower bound for value (default is -numpy.inf, no lower bound). * max (float, optional) – Upper bound for value (default is numpy.inf, no upper bound). * expr (str, optional) – Mathematical expression used to constrain the value during the fit (default is None).

We will now change the vary attribute of the redshift parameters ‘SiIV_A_z’ of model function 0, and ‘SiIV_B_z’ of model function 1 to True.

[14]:

# Make the redshifts variable parameters
params = siiv_model.params_list[0]
params['SiIV_A_z'].vary = True
params = siiv_model.params_list[1]
params['SiIV_B_z'].vary = True

Then we fit the SiIV SpecModel.

[15]:

# Fit the SiIV SpecModel
siiv_model.fit()
# Display the fitted SpecModel
siiv_model.plot()

Fitting the CIV line

In the next step we add the CIV model, basically repeating the same procedure.

[16]:

# Add the CIV emission line model
fit.add_specmodel()
civ_model = fit.specmodels[2]
civ_model.name = 'CIV_line'

civ_model.add_wavelength_range_to_fit_mask(6240, 6700)

model_name = 'CIV (2G components)'
model_prefix = None
civ_model.add_model(model_name, model_prefix, amplitude=10)

# Make the redshift a variable parameter
params = civ_model.params_list[0]
params['CIV_A_z'].vary = True
params = civ_model.params_list[1]
params['CIV_B_z'].vary = True

civ_model.fit()

civ_model.plot()

[INFO] Manual mask range 6240 6700

Fitting the CIII] line

After we successfully added the CIV line to the fit, we will now add a model for the SiIII], AlIII, and CIII] lines. While the previous SiIV and CIV model functions were comprised of two individual model functions called ‘line_model_gaussian’, the CIII] complex model is one model function comprised of three Gaussians with set central wavelength values.

[17]:

# Add the CIII] complex emission line model
fit.add_specmodel()
ciii_model = fit.specmodels[3]
ciii_model.name = 'CIII]_complex'

ciii_model.add_wavelength_range_to_fit_mask(7800, 8400)

model_name = 'CIII] complex (3G components)'
model_prefix = None
ciii_model.add_model(model_name, model_prefix, amplitude=2)

params = ciii_model.params_list[0]
params['CIII_z'].vary = True

# Initial fit
ciii_model.fit()
# Second fit to make sure the model converged
ciii_model.fit()

# Plot the model
ciii_model.plot()

[INFO] Manual mask range 7800 8400

Adding a basic line model (Gaussian) - Fitting absorption lines

Just blueward of the SiIV line this quasar has two prominent absorption features, which we also want to model. For this purpose we use the basic ‘Line model Gaussian’ model function and provide the parameters values as keyword arguments ourselves.

[18]:

# Add absorption line models
fit.add_specmodel()
abs_model = fit.specmodels[4]
abs_model.name = 'Abs_lines'

abs_model.add_wavelength_range_to_fit_mask(5760, 5790)

model_name = 'Line model Gaussian'
model_prefix = 'Abs_A'
abs_model.add_model(model_name, model_prefix, amplitude=-15,
                    cenwave=5766, fwhm=200, redshift=0)
model_name = 'Line model Gaussian'
model_prefix = 'Abs_B'
abs_model.add_model(model_name, model_prefix, amplitude=-15,
                    cenwave=5776, fwhm=200, redshift=0)

abs_model.fit()

abs_model.plot()

[INFO] Manual mask range 5760 5790

The default .plot() function dispersion and flux ranges may not always be appropriate to check whether the model was actually fit appropriately. Using the matplotlib notebook functionality one can zoom into the region of the two absorption lines to check whether the SpecModel fit did a good job.

Visualizing the full quasar model fit

Let us now visualize the entire fit by calling the SpecFit plot function:

[19]:

fit.plot()

Accessing SpecModel fit results and saving them

The fitted quasar model looks reasonable. However, in order to analyze the model we need to understand how to access the fitted parameters. Once a SpecModel has been fit, we can access the fit results simply by calling the .fit_result attribute of the SpecModel. It returns the LMFIT fit report for the SpecModel, giving us insight into the model that was fit, the fit statistics, the variables, and the correlations.

[20]:

# Printing the fit result for the CIV model
civ_model.fit_result

[20]:

Model

(Model(line_model_gaussian, prefix='CIV_A_') + Model(line_model_gaussian, prefix='CIV_B_'))

Fit Statistics

fitting method	leastsq
# function evals	121
# data points	309
# variables	6
chi-square	431.940013
reduced chi-square	1.42554460
Akaike info crit.	115.498142
Bayesian info crit.	137.898189

Variables

name	value	standard error	relative error	initial value	min	max	vary
CIV_A_z	3.20972764	0.00149858	(0.05%)	3.227	3.06565000	3.38835000	True
CIV_A_flux	1820.82409	100.100791	(5.50%)	100	0.00000000	26611675.5	True
CIV_A_cen	1549.06000	0.00000000	(0.00%)	1549.06	-inf	inf	False
CIV_A_fwhm_km_s	12122.2834	611.964256	(5.05%)	2500	300.000000	20000.0000	True
CIV_B_z	3.20984418	7.2015e-04	(0.02%)	3.227	3.06565000	3.38835000	True
CIV_B_flux	1151.33782	114.160367	(9.92%)	100	0.00000000	26611675.5	True
CIV_B_cen	1549.06000	0.00000000	(0.00%)	1549.06	-inf	inf	False
CIV_B_fwhm_km_s	4791.00832	222.068758	(4.64%)	2500	300.000000	20000.0000	True

Correlations (unreported correlations are < 0.100)

CIV_A_flux	CIV_B_flux	-0.9571
CIV_A_fwhm_km_s	CIV_B_flux	0.9133
CIV_B_flux	CIV_B_fwhm_km_s	0.9105
CIV_A_flux	CIV_B_fwhm_km_s	-0.9067
CIV_A_flux	CIV_A_fwhm_km_s	-0.8203
CIV_A_fwhm_km_s	CIV_B_fwhm_km_s	0.7812
CIV_A_z	CIV_B_z	-0.5566
CIV_A_z	CIV_A_fwhm_km_s	0.2409
CIV_A_z	CIV_B_flux	0.2348
CIV_A_z	CIV_A_flux	-0.2255
CIV_A_z	CIV_B_fwhm_km_s	0.2012
CIV_A_fwhm_km_s	CIV_B_z	-0.1323
CIV_B_z	CIV_B_flux	-0.1230
CIV_B_z	CIV_B_fwhm_km_s	-0.1201
CIV_A_flux	CIV_B_z	0.1194

We can save the fit result of each individual SpecModel using the save_fit_report function. One has to supply the folder path where the fit report should be saved. As an example we can save the fit report for the CIV SpecModel to the current folder.

[21]:

# Save the CIV SpecModel fit report to the current folder
civ_model.save_fit_report('.')

# Check if the fit report was saved
! ls

TestSpectralBroadening.ipynb      specmodel_CIV_line_fit_report.txt
scripting_sculptor_1.ipynb        speconed_demonstration.ipynb
scripting_sculptor_2.ipynb        spectrum_preparation.ipynb
scripting_sculptor_3.ipynb

Let’s have a brief look at what was saved exactly:

[22]:

h = open("specmodel_CIV_line_fit_report.txt", "r")
for line in h:
    print (line)
h.close()

[[Fit Statistics]]

    # fitting method   = leastsq

    # function evals   = 121

    # data points      = 309

    # variables        = 6

    chi-square         = 431.940013

    reduced chi-square = 1.42554460

    Akaike info crit   = 115.498142

    Bayesian info crit = 137.898189

[[Variables]]

    CIV_A_z:          3.20972764 +/- 0.00149858 (0.05%) (init = 3.227)

    CIV_A_flux:       1820.82409 +/- 100.100791 (5.50%) (init = 100)

    CIV_A_cen:        1549.06 (fixed)

    CIV_A_fwhm_km_s:  12122.2834 +/- 611.964256 (5.05%) (init = 2500)

    CIV_B_z:          3.20984418 +/- 7.2015e-04 (0.02%) (init = 3.227)

    CIV_B_flux:       1151.33782 +/- 114.160367 (9.92%) (init = 100)

    CIV_B_cen:        1549.06 (fixed)

    CIV_B_fwhm_km_s:  4791.00832 +/- 222.068758 (4.64%) (init = 2500)

[[Correlations]] (unreported correlations are < 0.100)

    C(CIV_A_flux, CIV_B_flux)           = -0.957

    C(CIV_A_fwhm_km_s, CIV_B_flux)      =  0.913

    C(CIV_B_flux, CIV_B_fwhm_km_s)      =  0.910

    C(CIV_A_flux, CIV_B_fwhm_km_s)      = -0.907

    C(CIV_A_flux, CIV_A_fwhm_km_s)      = -0.820

    C(CIV_A_fwhm_km_s, CIV_B_fwhm_km_s) =  0.781

    C(CIV_A_z, CIV_B_z)                 = -0.557

    C(CIV_A_z, CIV_A_fwhm_km_s)         =  0.241

    C(CIV_A_z, CIV_B_flux)              =  0.235

    C(CIV_A_z, CIV_A_flux)              = -0.225

    C(CIV_A_z, CIV_B_fwhm_km_s)         =  0.201

    C(CIV_A_fwhm_km_s, CIV_B_z)         = -0.132

    C(CIV_B_z, CIV_B_flux)              = -0.123

    C(CIV_B_z, CIV_B_fwhm_km_s)         = -0.120

    C(CIV_A_flux, CIV_B_z)              =  0.119

Wonderful! The full fit report, which we displayed above, was saved in the .txt file ‘specmodel_CIV line_fit_report.txt’. Let’s remove the file before we proceed.

[23]:

! rm specmodel_CIV_line_fit_report.txt
! ls

TestSpectralBroadening.ipynb scripting_sculptor_3.ipynb
scripting_sculptor_1.ipynb   speconed_demonstration.ipynb
scripting_sculptor_2.ipynb   spectrum_preparation.ipynb

Full fits, fitting algorithms, and saving your results

Performing a global consecutive fit and saving the SpecModel fit results

In some cases it might be useful to fit individual SpecModels and save their results. However, in most cases we want to fit all SpecModels consecutively and save all results. For this purpose the SpecFit object has a function called fit. It can take a keyword argument save_results, which defaults to False. If we set it to True the fit results will automatically be saved to the current folder. By providing a different foldername keyword argument the user can choose the folder, where the fit results will be saved.

[24]:

# Fit all SpecModels consecutively and save their results to the current folder
fit.fit(save_results=True)

# Check the current folder for the .txt files
! ls

TestSpectralBroadening.ipynb      specmodel_2_FitAll_fit_report.txt
scripting_sculptor_1.ipynb        specmodel_3_FitAll_fit_report.txt
scripting_sculptor_2.ipynb        specmodel_4_FitAll_fit_report.txt
scripting_sculptor_3.ipynb        speconed_demonstration.ipynb
specmodel_0_FitAll_fit_report.txt spectrum_preparation.ipynb
specmodel_1_FitAll_fit_report.txt

Note that in this case, the SpecModels have not been saved with their given names, but rather with their list inidices in the specfit.specmodels list. However, if you used easy to understand model prefixes a simple look into the .txt files will let you recover the fit parameters easily.

To keep the notebook directory clean, let’s remove the files again:

[25]:

# Remove the SpecModel fit report files
! rm *.txt
# Check the directory again
! ls

TestSpectralBroadening.ipynb scripting_sculptor_3.ipynb
scripting_sculptor_1.ipynb   speconed_demonstration.ipynb
scripting_sculptor_2.ipynb   spectrum_preparation.ipynb

Choosing the fit algorithm available in LMFIT

The SpecFit object has a string attribute called fitting_method, which allows you to specify with which algorithm LMFIT will fit your model to the data. An overview over which algorithms are implemented in LMFIT can be found here. The list of names available using Sculptor is saved as a global variable dictionary fitting_methods in the SpecModel module, that translates the human readable names of the algorithms into the LMFIT method options. Note that not all of them are fully tested in the Sculptor framework and may lead to errors, if used inappropriately.

[26]:

for key in scmod.fitting_methods:
    print('Name: {} \n Method {}'.format(key, scmod.fitting_methods[key]))

Name: Levenberg-Marquardt
 Method leastsq
Name: Nelder-Mead
 Method nelder
Name: Maximum likelihood via Monte-Carlo Markov Chain
 Method emcee
Name: Least-Squares minimization
 Method least_squares
Name: Differential evolution
 Method differential_evolution
Name: Brute force method
 Method brute
Name: Basinhopping
 Method basinhopping
Name: Adaptive Memory Programming for Global Optimization
 Method ampgo
Name: L-BFGS-B
 Method lbfgsb
Name: Powell
 Method powell
Name: Conjugate-Gradient
 Method cg
Name: Cobyla
 Method cobyla
Name: BFGS
 Method bfgs
Name: Truncated Newton
 Method tnc
Name: Newton GLTR trust-region
 Method trust-krylov
Name: Trust-region for constrained obtimization
 Method trust-constr
Name: Sequential Linear Squares Programming
 Method slsqp
Name: Simplicial Homology Global Optimization
 Method shgo
Name: Dual Annealing Optimization
 Method dual_annealing

For more information on the different fitting algorithms, please refer to the LMFIT documentation.

The default method of Sculptor’s SpecFit class is always ‘Levenberg-Marquardt’. However, this can be easily changed:

[27]:

# Current fitting method
print('Default fitting method: ', fit.fitting_method)

# Change fitting method to Nelder-Mead
fit.fitting_method = 'Nelder-Mead'

# Check fitting method again
print('New fitting method: ', fit.fitting_method)

# Fit all SpecModels
fit.fit()
# Display result
fit.plot()

Default fitting method:  Levenberg-Marquardt
New fitting method:  Nelder-Mead

Some of the fitting methods need extra parameters. A prominent example for this is the ‘Maximum likelihood via Monte-Carlo Markov Chain’, which uses emcee to perform a maximum likelihood fit using MCMC. For this particular method default parameters are implemented. However, we will devote an entire notebook to show how you can use this fitting option to produce science grade results.

Saving the model fit

One of the goals of the Sculptor package is to enable easy reproducibility of model fits to astronomical spectra and their analysis. Therefore, you can save your entire fit (SpecFit object) to a folder with the SpecFit save function. It takes the folderpath+foldername as an attribute and will create the folder if it does not find it in the specified directory.

[28]:

# Quick look into the current directory
! ls

# Save the SpecFit object
fit.save('example_fit_notebook')

# Check if the folder was created and contents were saved
! ls

TestSpectralBroadening.ipynb scripting_sculptor_3.ipynb
scripting_sculptor_1.ipynb   speconed_demonstration.ipynb
scripting_sculptor_2.ipynb   spectrum_preparation.ipynb

/opt/anaconda3/envs/sculptor-env/lib/python3.9/site-packages/pandas/core/generic.py:2606: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed-integer,key->block0_values] [items->Index(['value'], dtype='object')]

  pytables.to_hdf(

[INFO] Saving SpecModel fit result
[INFO] Saving new model file: example_fit_notebook/0_PL__model.json
[INFO] Saving SpecModel fit result
[INFO] Saving new model file: example_fit_notebook/1_SiIV_A__model.json
[INFO] Saving new model file: example_fit_notebook/1_SiIV_B__model.json
[INFO] Saving SpecModel fit result
[INFO] Saving new model file: example_fit_notebook/2_CIV_A__model.json
[INFO] Saving new model file: example_fit_notebook/2_CIV_B__model.json
[INFO] Saving SpecModel fit result
[INFO] Saving new model file: example_fit_notebook/3_CIII__model.json
[INFO] Saving SpecModel fit result
[INFO] Saving new model file: example_fit_notebook/4_Abs_A_model.json
[INFO] Saving new model file: example_fit_notebook/4_Abs_B_model.json
TestSpectralBroadening.ipynb scripting_sculptor_3.ipynb
example_fit_notebook         speconed_demonstration.ipynb
scripting_sculptor_1.ipynb   spectrum_preparation.ipynb
scripting_sculptor_2.ipynb

[29]:

# Check the contents of the saved SpecFit folder
! ls ./example_fit_notebook/

0_PL__model.json          2_fitresult.json          specmodel_0_specdata.hdf5
0_fitresult.json          3_CIII__model.json        specmodel_1_specdata.hdf5
1_SiIV_A__model.json      3_fitresult.json          specmodel_2_specdata.hdf5
1_SiIV_B__model.json      4_Abs_A_model.json        specmodel_3_specdata.hdf5
1_fitresult.json          4_Abs_B_model.json        specmodel_4_specdata.hdf5
2_CIV_A__model.json       4_fitresult.json          spectrum.hdf5
2_CIV_B__model.json       fit.hdf5

At this point we do not want to go into detail into all the saved files. The fit.hdf5 saved the main attributes of the SpecFit and SpecModel objects. The individual models and their fit results were saved as .json files via the LMFIT functionality for saving models, parameters, and results. The input spectrum is saved as the spectrum.hdf5 file and the spectra of the SpecModel objects, along with the mask regions and model spectra are saved in the specmodel_X_specdata.hdf5 files.

Loading a SpecFit object from disk

We can now use the Sculptor GUI to now load the full model fit (SpecFit object) into the GUI from the folder we saved it to and then manipulate it to tune the fit.

Of course we can also instantiate a new SpecFit object and then load the previous result. This is quite simple:

[30]:

# Instantiate an empty SpecFit object
new_fit = scfit.SpecFit()

# Load the saved model fit using the folder name
new_fit.load('example_fit_notebook')

# Fit all SpecModels and then display the fit to see if everything worked
new_fit.fit()
new_fit.plot()

# We can also check if it saved the fitting method we changed earlier
print(new_fit.fitting_method)

Nelder-Mead

The full fit you see displayed above should resemble the one a few cells earlier, where we tested changing the fitting method. With this we conclude the demonstration on how to use scripts to construct, fit and save Sculptor model fits.

(We delete the save Sculptor fit to keep the notebook directory clean)

[31]:

# Delete the Sculptor model fit folder
! rm -r example_fit_notebook
# Check the directory
! ls

TestSpectralBroadening.ipynb scripting_sculptor_3.ipynb
scripting_sculptor_1.ipynb   speconed_demonstration.ipynb
scripting_sculptor_2.ipynb   spectrum_preparation.ipynb

[ ]: