Sherpa Template Models

Sherpa Threads (CIAO 4.11 Sherpa v1)

Overview

Synopsis:

The Sherpa template model is available for comparing a source data spectrum against a set of user-provided template models, in order to find the single template model which best matches the source data. A special grid-search fit optimization method, gridsearch, is used to fit a template model to a data set in Sherpa.

Starting with a directory of template model files conforming to a specific format, and a single index file listing the file contents of that directory, the models are loaded into the Sherpa session using the load_template_model function, an extension of load_table_model. After fitting the template model to the source data using the gridsearch fit optimization method, the parameters of the best matched template model are returned. Fit results may be examined in the usual way with the appropriate plotting commands.

The Sherpa template model supports linear, nearest-neighbor, and polynomial interpolation. Interpolation is used by the template model to match the data grid to the model grid—which must match before the fit statistics can be calculated for fitting.

Last Update: 14 Dec 2018 - reviewed for CIAO 4.11, no content change.

Getting Started: Setting up the Template Model Files

Model files

A collection of template models may be read into Sherpa from a directory full of template files using the load_template_model function, in order to compare a data set to all of the templates in that collection. The Sherpa model library accommodates ASCII template model files containing separated columns of x- and y-coordinates, such as wavelength, energy, or time, versus flux, counts, or intensity. In this thread, we consider a set of accretion disk models in the Kerr geometry (Siemiginowska et al., ApJ, 454, p.77, 1995), where there is a separate model file for each combination of the black hole mass, accretion rate, and inclination angle parameter values. All the data used in this thread can be also found in the tar file: template_model.tar.gz. Download and unpack this file, and go into the "template_model" directory before continuing.

The contents of the directory of Kerr model files is shown below:

% ls data/Kerr/

k10_001_01.dat   k6_01_05.dat   k7_01_1.dat    k8_03_025.dat   k9_01_075.dat
k10_001_025.dat  k6_01_075.dat  k7_03_01.dat   k8_03_05.dat    k9_01_1.dat
k10_001_05.dat   k6_01_1.dat    k7_03_025.dat  k8_03_075.dat   k9_03_01.dat
k10_001_075.dat  k6_03_01.dat   k7_03_05.dat   k8_03_1.dat     k9_03_025.dat
k10_001_1.dat    k6_03_025.dat  k7_03_075.dat  k8_08_01.dat    k9_03_05.dat
k10_01_01.dat    k6_03_05.dat   k7_03_1.dat    k8_08_025.dat   k9_03_075.dat
k10_01_025.dat   k6_03_075.dat  k7_08_01.dat   k8_08_05.dat    k9_03_1.dat
k10_01_05.dat    k6_03_1.dat    k7_08_025.dat  k8_08_075.dat   k9_08_01.dat
k10_01_075.dat   k6_08_01.dat   k7_08_05.dat   k8_08_1.dat     k9_08_025.dat
k10_01_1.dat     k6_08_025.dat  k7_08_075.dat  k8_1_01.dat     k9_08_05.dat
k10_03_01.dat    k6_08_05.dat   k7_08_1.dat    k8_1_025.dat    k9_08_075.dat
k10_03_025.dat   k6_08_075.dat  k7_1_01.dat    k8_1_05.dat     k9_08_1.dat
k10_03_05.dat    k6_08_1.dat    k7_1_025.dat   k8_1_075.dat    k9_1_01.dat
k10_03_075.dat   k6_1_01.dat    k7_1_05.dat    k8_1_1.dat      k9_1_025.dat
k10_03_1.dat     k6_1_025.dat   k7_1_075.dat   k9_001_01.dat   k9_1_05.dat
k10_08_01.dat    k6_1_05.dat    k7_1_1.dat     k9_001_025.dat  k9_1_075.dat
k10_08_025.dat   k6_1_075.dat   k8_01_01.dat   k9_001_05.dat   k9_1_1.dat
k10_08_05.dat    k6_1_1.dat     k8_01_025.dat  k9_001_075.dat  templ_index.dat
k10_08_075.dat   k7_01_01.dat   k8_01_05.dat   k9_001_1.dat
k10_08_1.dat     k7_01_025.dat  k8_01_075.dat  k9_01_01.dat
k6_01_01.dat     k7_01_05.dat   k8_01_1.dat    k9_01_025.dat
k6_01_025.dat    k7_01_075.dat  k8_03_01.dat   k9_01_05.dat



and the contents of a selected file, k10_001_05.dat is shown below.

% more data/Kerr/k10_001_05.dat

13.3830    43.6970
13.4330    43.8500
13.4830    43.9320
13.5330    44.0030
13.5830    44.0720
13.6330    44.1400
13.6830    44.2080
13.7330    44.2760
.
.
.
16.4830    33.9740
16.5330    32.2860
16.5830    30.3930
16.6330    28.2710
16.6830    25.8880
16.7330    23.2100
16.7830    20.1980
16.8330    16.8090
16.8830    12.9930


Each of the example ASCII-format template files contains columns of frequency and luminosity values defining the predicted spectra emitted by an accretion disk surrounding a supermassive black hole (see section 4.1 of the referenced paper), where local modified blackbody emission is assumed. The first column has units of log10 frequency in Hz, and the second is the luminosity in logarithm of the rest-frame, integrated luminosity—$$\log_{10}\left(\nu L_{\nu}\right)$$—in erg/s.

Index File

The template model index file should contain a table with one line per template data file, with three groups of columns in the following order:

• Model parameter columns, an arbitrary number of columns
• The MODELFLAG column which separates the parameter list from the filenames/model arrays, and marks lines which are to be used or not: MODELFLAG = 1—use the file; MODELFLAG = 0—do not use the file
• The FILENAME column which points to the data file for that instance, including the absolute path to the file.

Sherpa reads the index file in order to set up the template model with the parameters specified in the first line, and the arrays from the columns given by the data files.

Our example index file appears below, with model flags (all having value 1) and filenames listed in the two right-most columns, and the black hole mass ($$\log_{10} M_{bh}/M_{\odot}$$), accretion rate (Eddington), and inclination angle ($$\cos \theta$$) model parameters listed from the left. The combination of the three parameter values sets the spectral shape of the model.

% more templ_index.dat

# mass rate angle modelflag filename

6.0 0.1 0.1    1   data/Kerr/k6_01_01.dat
6.0 0.1 0.25   1   data/Kerr/k6_01_025.dat
6.0 0.1 0.5    1   data/Kerr/k6_01_05.dat
6.0 0.1 0.75   1   data/Kerr/k6_01_075.dat
6.0 0.1 1.0    1   data/Kerr/k6_01_1.dat

6.0 0.3 0.1    1   data/Kerr/k6_03_01.dat
6.0 0.3 0.25   1   data/Kerr/k6_03_025.dat
6.0 0.3 0.5    1   data/Kerr/k6_03_05.dat
6.0 0.3 0.75   1   data/Kerr/k6_03_075.dat
6.0 0.3 1.0    1   data/Kerr/k6_03_1.dat

.
.
.

10.0 0.3 0.1    1   data/Kerr/k10_03_01.dat
10.0 0.3 0.25   1   data/Kerr/k10_03_025.dat
10.0 0.3 0.5    1   data/Kerr/k10_03_05.dat
10.0 0.3 0.75   1   data/Kerr/k10_03_075.dat
10.0 0.3 1.0    1   data/Kerr/k10_03_1.dat

10.0 0.8 0.1    1   data/Kerr/k10_08_01.dat
10.0 0.8 0.25   1   data/Kerr/k10_08_025.dat
10.0 0.8 0.5    1   data/Kerr/k10_08_05.dat
10.0 0.8 0.75   1   data/Kerr/k10_08_075.dat
10.0  0.8 1.0    1   data/Kerr/k10_08_1.dat


After starting Sherpa, we read in our source data spectrum to be fitted by the template model in the usual way, using the appropriate load command. In this example, we use the load_ascii command to read the first three columns of ASCII data from file source.dat, namely the x (frequency), y (integrated luminosity), and y error columns.

% ciao

% sherpa
CIAO 4.11 Sherpa version 1 Wednesday, December 5, 2018

Python 3.5.4 (default, Oct 15 2018, 13:47:46)
IPython 6.5.0 -- An enhanced Interactive Python. Type '?' for help.

IPython profile: sherpa



The show_data and get_indep/get_dep commands may be used to check that the data was properly read:

sherpa> show_data()
Data Set: 1
Filter: 17.6983-7.4149 x
name      = source.dat
x         = Float64[46]
y         = Float64[46]
staterror = Float64[46]
syserror  = None

sherpa> get_indep()
(array([ 17.69831459,  15.22778313,  15.14764634,  15.04849851,
14.69471047,  14.57931705,  14.4512198 ,  14.2154123 ,
10.4876855 ,  10.34388301,  10.01346923,  10.01346923,
10.00024097,  10.00024097,   9.74586299,   9.74586299,
9.46062726,   9.46062726,   9.46062726,   9.46062726,
9.46062726,   9.46062726,   9.18956049,   9.18956049,
9.18956049,   8.92515939,   8.87679209,   8.56491923,
8.56491923,   8.56491923,   8.51861921,   8.49634282,
8.49634282,   8.24899768,   8.21758921,   8.09265048,
7.89428282,   7.89428282,   7.73445498,   7.73445498,
7.71243924,   7.6608522 ,   7.61552922,   7.5372157 ,
7.48181656,   7.41486977]),)

sherpa> get_dep()
array([ 44.83025145,  46.20275955,  46.03082425,  46.04305871,
45.63597532,  45.47709506,  45.46993128,  45.41821408,
43.73912372,  43.74154589,  44.06194405,  44.0665398 ,
43.96097247,  43.95510354,  43.98644329,  43.9809108 ,
44.02583231,  44.03575947,  44.05023729,  44.05023729,
44.08004227,  44.04785759,  44.07782081,  44.10165342,
44.10165342,  44.14189464,  44.16799744,  44.18745992,
44.22553677,  44.15139818,  44.18572621,  44.16069821,
44.16436316,  44.24415224,  44.19602617,  44.18918674,
44.05729915,  44.13020888,  43.94007455,  43.94007455,
44.07334351,  43.86696745,  43.87414603,  44.02331465,
44.05495793,  43.96676315])


We may also visualize the data before fitting with plot_data, and utilize several ChIPS commands to customize the data plot.

sherpa> plot_data()

sherpa> linear_scale()
sherpa> set_plot_xlabel(r"log \nu [Hz]")
sherpa> set_plot_ylabel(r"log \nu L_{\nu} [erg/s]")


We may now use the load_template_model function to read in our collection of templates and assign them to a Sherpa model instance, as demonstrated below.

sherpa> load_template_model("kerr_templ", "templ_index.dat")
sherpa> set_model(kerr_templ)


Here, we load the Kerr disk model templates listed in the index file templ_index.dat into the Sherpa model instance kerr_templ, and assign this model to our source data set with set_model.

The template interpolator may also be disabled

sherpa> load_template_model("kerr_templ", "templ_index.dat", template_interpolator_name=None)


but requires the fit method to be gridsearch or explicitly set, by importing the interpolator and setting the 'template_interpolator_name' argument:

sherpa> from sherpa.utils import neville, nearest_interp, linear_interp
sherpa> from sherpa.models import KNNInterpolator



We can examine the model parameters using show_model, where the parameter names contained in our index file are shown, along with the parameter values matching the first template in our list.

sherpa> show_model()
Model: 1
templatemodel.kerr_templ
Param               Type          Value          Min          Max      Units
-----               ----          -----          ---          ---      -----
kerr_templ.mass     thawed            6            6           10
kerr_templ.rate     thawed          0.1         0.01            1
kerr_templ.angle    thawed          0.1          0.1            1


Using the Grid-search Optimization Method

In order to fit a Sherpa template model to a data set, gridsearch must be chosen as the fit optimization method. We do this using set_method, after checking the current method setting with show_method:

sherpa> show_method()
Optimization Method: LevMar
name    = levmar
ftol    = 1.19209289551e-07
xtol    = 1.19209289551e-07
gtol    = 1.19209289551e-07
maxfev  = None
epsfcn  = 1.19209289551e-07
factor  = 100.0
verbose = 0

sherpa> set_method("gridsearch")
sherpa> set_method_opt('sequence', kerr_templ.parvals)

sherpa> show_method()
Optimization Method: GridSearch
name     = gridsearch
num      = 16
sequence = Float64[315]
numcores = 1
maxfev   = None
ftol     = 1.19209289551e-07
method   = None
verbose  = 0


The grid-search method evaluates the fit statistic for each point in the parameter-space grid provided by the user, which is specified as a list in the 'sequence' parameter (the list which the method should iterate through to evaluate the model function). In this example, we have set the sequence parameter to use the full list of mass, rate, and angle parameter values associated with the template model files in our collection. The list of parameter values may be accessed as follows:

sherpa> print(kerr_templ.parvals)

[[  6.     0.1    0.1 ]
[  6.     0.1    0.25]
[  6.     0.1    0.5 ]
[  6.     0.1    0.75]
[  6.     0.1    1.  ]
[  6.     0.3    0.1 ]
[  6.     0.3    0.25]
[  6.     0.3    0.5 ]
[  6.     0.3    0.75]
[  6.     0.3    1.  ]

.       .      .
.       .      .
.       .      .

[ 10.     0.3    0.1 ]
[ 10.     0.3    0.25]
[ 10.     0.3    0.5 ]
[ 10.     0.3    0.75]
[ 10.     0.3    1.  ]
[ 10.     0.8    0.1 ]
[ 10.     0.8    0.25]
[ 10.     0.8    0.5 ]
[ 10.     0.8    0.75]
[ 10.     0.8    1.  ]]



The best match to the source spectrum in this sequence is the grid point with the lowest value of the fit statistic; gridsearch returns the parameter values associated with this point.

Note

If the sequence parameter is kept at the default value of None, the gridsearch num parameter is used to generate a sequence of parameters to evaluate. A uniform grid of size nparnum elements is generated.

Fitting the Template Model to the Spectrum

Once the source data has been loaded, the template model properly set, and the fit optimization method switched to gridsearch and customized, we may fit the template model to the source spectrum, with the supplied errors included, in the usual way using the fit() command.

sherpa> show_stat()
Statistic: Chi2Gehrels
Chi Squared with Gehrels variance.

The variance is estimated from the number of counts in each bin,
but unlike Chi2DataVar, the Gaussian approximation is not
used. This makes it more-suitable for use with low-count data.

The standard deviation for each bin is calculated using the
approximation from [1]_:

sigma(i,S) = 1 + sqrt(N(i,s) + 0.75)

where the higher-order terms have been dropped. This is accurate
to approximately one percent. For data where the background has
not been subtracted then the error term is:

sigma(i) = sigma(i,S)

whereas with background subtraction,

sigma(i)^2 = sigma(i,S)^2 + [A(S)/A(B)]^2 sigma(i,B)^2

Notes
-----
The accuracy of the error term when the background has been
subtracted has not been determined. A preferable approach to
background subtraction is to model the background as well as the
source signal.

References
----------

.. [1] "Confidence limits for small numbers of events in
astrophysical data", Gehrels, N. 1986, ApJ, vol 303,
p. 336-346.

sherpa> notice(13.,17.)

sherpa> fit()
Dataset               = 1
Method                = gridsearch
Statistic             = chi2
Initial fit statistic = 426180
Final fit statistic   = 97.4647 at function evaluation 106
Data points           = 7
Degrees of freedom    = 4
Probability [Q-value] = 3.40753e-20
Reduced statistic     = 24.3662
Change in statistic   = 426083
kerr_templ.mass   9
kerr_templ.rate   0.3
kerr_templ.angle   1
WARNING: parameter value kerr_templ.angle is at its maximum boundary 1.0



The initial fit results indicate that, of all the template models in our collection, the one contained in template model file k9_03_1.dat represents the best match to the source data.

We can further examine the fit results using one of the plotting commands designed for this puropse, such as plot_fit_delchi or plot_fit_resid, to visualize the quality of the fit.

sherpa> plot_fit_resid()

sherpa> set_plot_xlabel(r"log \nu [Hz]")
sherpa> set_current_plot("plot1")
sherpa> set_plot_ylabel(r"log \nu L_{\nu} [erg/s]")
sherpa> set_current_plot("plot2")
sherpa> set_plot_ylabel("Residuals")

Note

The confidence and projection commands available for estimating errors on model parameters, such as confidence and region projection, are incompatible with the gridsearch optimization method used for fitting template models to source data.

Similarly, with the interpolator enabled, other fitting optimization methods are useable besides gridsearch, as demonstrated below.

sherpa> clear()
sherpa> set_model(kerr_templ)

sherpa> fit()
Dataset               = 1
Statistic             = chi2
Initial fit statistic = 426180
Final fit statistic   = 91.4691 at function evaluation 311
Data points           = 7
Degrees of freedom    = 4
Probability [Q-value] = 6.4175e-19
Reduced statistic     = 22.8673
Change in statistic   = 426089
kerr_templ.mass   8.9899
kerr_templ.rate   0.306754
kerr_templ.angle   0.995986


The fit results may also be saved to an ASCII file, in this case fit.out, which can be examined to see the parameter space that was explored by Sherpa.

sherpa> fit(outfile="fit.out")
Dataset               = 1
Statistic             = chi2
Initial fit statistic = 91.4691
Final fit statistic   = 91.4691 at function evaluation 225
Data points           = 7
Degrees of freedom    = 4
Probability [Q-value] = 6.4175e-19
Reduced statistic     = 22.8673
Change in statistic   = 1.70631e-08
kerr_templ.mass   8.98967
kerr_templ.rate   0.306241
kerr_templ.angle   0.995821

sherpa> !cat fit.out
# nfev statistic kerr_templ.mass kerr_templ.rate kerr_templ.angle
0.000000e+00 9.146914e+01 8.989896e+00 3.067536e-01 9.959857e-01
1.000000e+00 9.263496e+04 7.000000e+00 3.067536e-01 9.959857e-01
2.000000e+00 1.363064e+02 8.989896e+00 2.575000e-01 9.959857e-01
3.000000e+00 9.499398e+01 8.989896e+00 2.821268e-01 9.959857e-01
4.000000e+00 2.263918e+04 7.994948e+00 3.067536e-01 9.959857e-01
5.000000e+00 4.914621e+04 8.326597e+00 2.903357e-01 3.959857e-01
...
2.200000e+02 9.146914e+01 8.989844e+00 3.065936e-01 9.958963e-01
2.210000e+02 9.146914e+01 8.989859e+00 3.065486e-01 9.957764e-01
2.220000e+02 9.146914e+01 8.989950e+00 3.067572e-01 9.958588e-01
2.230000e+02 9.146914e+01 8.989925e+00 3.067376e-01 9.958894e-01
2.240000e+02 9.146914e+01 8.989673e+00 3.062410e-01 9.958209e-01


In the next example, we add a continuous model component to the template model, in this particular case, simply a constant. The default KNN-interpolation method is used and the fit is varied on the grid parameters, with the constant model component frozen.

sherpa> reset()
sherpa> notice(13.5,16.)
sherpa> set_method("gridsearch")
sherpa> set_method_opt("sequence",None)
sherpa> set_model(kerr_templ+const1d.c1)
sherpa> set_par(c1.c0, min=0.1, max=10)
sherpa> freeze(c1)

sherpa> fit()
Dataset               = 1
Method                = gridsearch
Statistic             = chi2
Initial fit statistic = 287265
Final fit statistic   = 83.1809 at function evaluation 4097
Data points           = 7
Degrees of freedom    = 4
Probability [Q-value] = 3.68823e-17
Reduced statistic     = 20.7952
Change in statistic   = 287182
kerr_templ.mass   8.93333
kerr_templ.rate   0.01
kerr_templ.angle   0.76
WARNING: parameter value kerr_templ.rate is at its minimum boundary 0.01

sherpa> show_model()
Model: 1
(template.kerr_templ + const1d.c1)
Param        Type          Value          Min          Max      Units
-----        ----          -----          ---          ---      -----
kerr_templ.mass thawed      8.93333            6           10
kerr_templ.rate thawed         0.01         0.01            1
kerr_templ.angle thawed         0.76          0.1            1
c1.c0        frozen            1          0.1           10


It should be noted that when using continuous model components, an interpolator is required in order to perform the grid-search, otherwise an error will be raised:

sherpa> reset()
sherpa> set_model(kerr_templ+const1d.c2)

sherpa> fit()

ModelErr: Interpolation of template parameters was disabled for this
model, but parameter values not in the template library have been
requested. Please use gridsearch method and make sure the sequence
option is consistent with the template library



History

 17 Jan 2012 The Sherpa template model is new in CIAO 4.4. 10 Dec 2013 updated for CIAO 4.6, new examples with continuous model component added, demonstrated setting interpolator. 17 Mar 2015 updated for CIAO 4.7, added missing commands dealing with loading in an interpolator. 09 Dec 2015 reviewed for CIAO 4.8, no content change. 16 Nov 2016 reviewed for CIAO 4.9; linked tarball with template models, and updated examples to use the tarball directory struction. Fixed typos. 29 May 2018 reviewed for CIAO 4.10, no content change. 14 Dec 2018 reviewed for CIAO 4.11, no content change.