Chandra X-Ray Observatory
	(CXC)
Skip to the navigation links
Last modified: 16 Nov 2016

URL: http://cxc.harvard.edu/sherpa/threads/template_model/

Sherpa Template Models

Sherpa Threads (CIAO 4.9 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: 16 Nov 2016 - reviewed for CIAO 4.9; linked tarball with template models, and updated examples to use the tarball directory struction. Fixed typos.


Contents


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: sherpa.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 Log frequency in Hz, and the second is the luminosity in log of the rest-frame, integrated luminosity—Log(νLν)—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 Mbh/Msun), accretion rate (Eddington), and inclination angle (cos θ) 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 

Loading the Source Data

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
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.9 Sherpa version 1 Friday, December 2, 2016

sherpa> load_ascii('source.dat', ncols=3, dstype=Data1D)

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]")

Loading the Template Models

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
sherpa> load_template_interpolator("KNN",KNNInterpolator)

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

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]
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.
           http://adsabs.harvard.edu/abs/1986ApJ...303..336G



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]
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> load_template_model("kerr_templ", "templ_index.dat", KNNInterpolator)
sherpa> set_model(kerr_templ)
sherpa> set_method("neldermead")

sherpa> fit()
Dataset               = 1
Method                = neldermead
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
Method                = neldermead
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> load_template_model("kerr_templ","templ_index.dat")
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> load_template_model("kerr_templ","templ_index.dat",template_interpolator_name=None)
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.


Last modified: 16 Nov 2016
Smithsonian Institute Smithsonian Institute

The Chandra X-Ray Center (CXC) is operated for NASA by the Smithsonian Astrophysical Observatory. 60 Garden Street, Cambridge, MA 02138 USA.   Email:   cxchelp@head.cfa.harvard.edu Smithsonian Institution, Copyright © 1998-2017. All rights reserved.