Come analyze HEASARC, IRSA, and MAST data in the cloud! The Fornax Initiative is now welcoming all interested beta users.
Next: Third-Party Libraries In Local Up: Adding Models to XSPEC Previous: Loading a new model Contents
Writing a new model function
A model function is a subroutine that calculates the model spectrum
given an input array of energy bins and an array of parameter
values. The input array of energy bins gives the boundaries of the
energy bins and hence has one more entry than the output flux
arrays. The energy bins are assumed to be contiguous and ascending, and will be
determined by the response matrix in use. The subroutine should thus
make no assumptions about the energy range and bin sizes. The output
flux array for an additive model should be in terms of photons/cm
/s
(not photons/cm/s/keV) i.e. it is the model spectrum integrated
over the energy bin. The output array for a multiplicative model is
the multiplicative factor for that bin. Convolution models are
operators on the output from additive or multiplicative models. Model
subroutines can be written in Fortran, either in single or double
precision, in C++ using either C++-style arguments or C-style
arguments, and in C.
The model.dat entry
In addition to the subroutine, XSPEC requires a text file describing the model and its parameters. The standard models are specified in the model.dat file so we usually refer to this text file by that name. A sample model.dat entry has the following form:
modelentry 4 0. 1.e20 funcName add 0 0 lowT keV 0.1 0.0808 0.0808 79.9 79.9 0.001 highT keV 4. 0.0808 0.0808 79.9 79.9 0.001 Abundanc " " 1. 0. 0. 5. 5. 0.01 *redshift " " 0.0
The first line for each model gives the model name, the number of parameters, the low and high energies for which the model is valid, the name of the subroutine to be called and the type of model (add, mul, mix, or con, or acn). The final two arguments are flags: the first should be set to 1 if model variances are calculated by funcName and the second should be set to 1 if the model should be forced to perform a calculation for each spectrum. This final flag is necessary because if multiple spectra have the same energy bins, the default behavior is to perform the model calculation for just one spectrum and copy the results for each of the others. However, if a model depends on information about the spectrum in addition to its energy ranges, it must be forced to perform a calculation for each spectrum.
The fifth entry in the first line, which includes the name of the subroutine to be called, must be identified in the model.dat file according to the language in which it is written. Table C.1 below shows the required format of the function names in model.dat.
The remaining lines in the text file specify each parameter in the model. For regular model parameters the first two fields are the parameter name followed by an optional units label. If there is no units label, then there must be a quoted blank (“ ”) placeholder. The remaining 6 numerical entries are the default parameter value, hard min, soft min, soft max, hard max, and fit delta, which are described in the newpar command section.
There are three special types of parameter which can be used. If the name of the parameter is prefixed with a “*” the parameter is a “scale” parameter and cannot be made variable or linked to any kind of parameter other than another scale parameter. Since the parameter value can never vary only the initial value need be given. If the name of the parameter is prefixed with a “$” the parameter is a “switch” parameter which is not used directly as part of the calculation, but switches the model component function's mode of operation (i.e. calculate or interpolate). Switch parameters only have 2 fields: the parameter name and an integer value.
Finally, if a P is added at the end of the line for a parameter then the parameter is defined to be periodic. During a fit, a periodic parameter will not be pegged if it tries to exceed its hard limits. Instead it will be assigned a value within its limits: f(max + delta) = f(min + delta), f(min-delta) = f(max-delta). The soft min and max settings are irrelevant for period parameters and will be ignored.
Another way to visualize the required fields in the model.dat file:
modelentry 4 0. 1.e20 C_funcName add 0 0
| | | | | | | | |
<model name> | | | | | | | |
| | | | | | | |
<num params> | | | | | | |
| | | | | | |
<valid low energy> | | | | | |
| | | | | |
<valid high energy> | | | | |
| | | | |
<function identifier> | | | |
[C_|c_|F_|U_|` '] | | | |
| | | |
<function name> | | |
| | |
<type> | |
[add|mul|con|mix|acn|amx] | |
| |
<flag> |
|
<flag>
lowT keV 0.1 0.0808 0.0808 79.9 79.9 0.001
| | | | | | | |
<param> | | | | | | |
| | | | | | |
<unit> | | | | | |
| | | | | |
<default> | | | | |
| | | | |
<hard min> | | | |
| | | |
<soft min> | | |
| | |
<soft max> | |
<hard max> |
<fit delta>
The model subroutine function
When writing the code for the new model, Xspec expects specific arguments. Table C.1 below lists the function arguments required for the different language options.
| Call Type and Specification Function Name Format in model.dat |
Arguments and Type | Meaning |
Single precision Fortran |
real*4 ear(0:ne) | Energy array |
| integer ne | Size of flux array | |
| real*4 param(*) | Parameter values (Dimension must be specified inside the function) | |
| integer ifl | The spectrum number of model component being calculated | |
| real*4 photar(ne) | Output flux array | |
| real*4 photer(ne) | Output flux error array (optional) | |
Double precision Fortran |
real*8 ear(0:ne) | Energy array |
| integer ne | Size of flux array | |
| real*8 param(*) | Parameter values (Dimension must be specified inside the function) | |
| integer ifl | The spectrum number of model component being calculated | |
| real*8 photar(ne) | Output flux array | |
| real*8 photer(ne) | Output flux error array (optional) | |
C and C++ in C style |
const Real* energy | Energy array (size Nflux+1) |
| int Nflux | Size of flux array | |
| const Real* parameter | Parameter values | |
| int spectrum | The spectrum number of model component being calculated | |
| Real* flux | Output flux array | |
| Real* fluxError | Output flux error array (optional) | |
| const char* init | Initialization string (see below) | |
C++ in C++ style |
const RealArray& energy | Energy array |
| const RealArray& parameter | Parameter values | |
| int spectrum | The spectrum number of model component being calculated | |
| RealArray& flux | Output flux array | |
| RealArray& fluxError | Output flux error array (optional) | |
| const string& init | Initialization string (see below) | |
The second entry in the first column is the way the function name should be included in the model.dat entry. Note that the prefix is only to be included in model.dat - the actual function name in the source code is the base funcName.
For example, a model component in double precision Fortran is specified by:
modelentry 5 0. 1.e20 F_funcName add 0 0XSPEC sees the F_ and picks out the right function definition, calling the Fortran function funcName, which expects double precision arguments. The C-style call can clearly be compiled and implemented by either a C or a C++ compiler. We recommend using the C++ call if the model is written in C++ as it will reduce overhead in copying C arrays in and out the XSPEC internal data structures. To prevent unresolved symbol linkage errors, we also recommend prefacing C++ local model function definitions with the extern "C" directive.
Example C/C++ function definitions:
/* C style */
extern "C"
void funcName(const Real* energy, int Nflux, const Real* parameter,
int spectrum, Real* flux, Real* fluxVariance,
const char* init)
{
/* Model code: Do not allocate memory for flux and fluxVariance arrays.
XSPEC's C-function wrapper will allocate arrays prior to calling the
function (and will free them afterwards). */
}
// C++ style
extern "C"
void funcName(const RealArray& energy, const RealArray& parameter,
int spectrum, RealArray& flux, RealArray& fluxVariance,
const string& init)
{
// Model code: Should resize flux RealArray to energy.size()-1.
// Do the same for fluxVariance array if calculating errors, otherwise
// leave it at size 0.
}
Note on type definitions for (C and C++): XSPEC provides a typedef for Real, in the xsTypes.h header file. The distributed code has
typedef double Real;i.e. all calculations are performed in double precision. This is used for C models and C++ models with C-style arguments.
The type RealArray is a dynamic (resizeable) array of type Real. XSPEC uses the std::valarray template class to implement RealArray. The internal details of XSPEC require that the RealArray typedef supports vectorized assignments and mathematical operations, and indirect addressing (see C++ documentation for details). However, we do not recommend using specific features of the std::valarray class, such as array slicing, in case the typedef is changed in future.
The input energies are set by the response matrices of the detectors in use. The parameter ifl / spectrum (for Fortran / C or C++) is an integer which specifies to which response (and therefore which spectrum) these energies correspond. It exists to allow multi-dimensional models where the function might also depend on e.g. pulse-phase in a variable source. The output flux array should not be assumed to have any particular values on input. It is assumed to contain previously calculated values only by convolution/pileup models, which have the nature of operators. The output flux error array allows the function to return model variances.
The C and C++ call types allow one extra argument, which is a character string that can be appended to the top line of the model component description. This string is read on initialization and available to the model during execution. An example of its use might be the name of a file with specific data used in the model calculation: this allows different models to be implemented the same way except for different input data by specifying different names and input strings.
Model functions in Python
Model functions may also be written in Python. In model.dat the function field then takes the form Py_<module>.<func> (for example Py_mymodels.lpow), naming a function func in a Python module <module>.py placed beside the compiled library; initpackage and lmod build and load such a package in the usual way. The function takes the same three-to-five arguments as the other call types — energies, parameters, flux, and optionally a flux-error array and the spectrum number — which XSPEC passes as Python tuples (energies and parameters) and a pre-sized list (flux) to be filled in place.
A Python model can opt into a faster vectorized interface by decorating the function with @xspec.vectorized: XSPEC then passes zero-copy numpy arrays — read-only energies and parameters, and a writable flux array to be filled in place with flux[:] = ... — which avoids the per-element Python-object overhead of the tuple/list interface. The decorator is available both in PyXspec and, via import xspec, inside a local model running in the standalone XSPEC. The same function may equally be registered from PyXspec with AllModels.addPyMod. See the “Local Models in Python” section of the PyXspec manual for the full interface, including the flux-error convention.
Function utility routines
There are a number of internal XSPEC routines which may be useful for local models written in either C++ or Fortran. These are defined in FunctionUtility.h and xsFortran.cxx, respectively, in the directory Xspec/src/XSFunctions/Utilities. Documentation is available at: https://heasarc.gsfc.nasa.gov/docs/software/xspec/internal/XspecInternalFunctionsGuide.html.
These functions allow the model to find out the Solar abundances in use, get any information defined using the xset command, find the cosmological parameters currently set, get the contents of the XFLT#### keywords in the spectrum files, and load information into an internal database which can be accessed using tclout.
There are also C++ classes and associated methods for equilibrium and non-equilibrium collisional plasmas, photo-electric opacity from neutral and ionized materials, and Compton scattering. Descriptions of these can be found in XSFunctions_guide.pdf in Xspec/src/help.
Multi-dimensional models
Typically, an XSPEC model depends only on the energy (and the model parameters) but more complicated models are sometimes required. The model might depend on energy and, for a variable source, the time at which the spectrum was obtained or, for an extended source, its spatial position. The extra information about the spectrum must be stored in the XFLT#### keywords in the spectrum files. For instance, if the model has a dependence on time then for one spectrum XFLT0001 could be set to “start: 1345.54” and XFLT0002 to “end: 1356.43” and another have XFLT0001 set to “start: 1356.43” and XFLT0002 to “end: 1368.21”. An example C++ routine might then include the code
#include <XSFunctions/Utilities/FunctionUtility.h>
...
extern "C"
void mymodel(const RealArray& energyArray, const RealArray& params,
int spectrum, RealArray& flux, RealArray& fluxVariance,
const string& initString)
...
Real tstart = FunctionUtility::getXFLT(spectrum,"start");
Real tend = FunctionUtility::getXFLT(spectrum,"end");
...
Supplying analytic gradients
A model without a registered gradient is differentiated by finite
differences: to build the gradient needed by fit (with the
Levenberg–Marquardt or migrad methods) and by hmc,
XSPEC perturbs each parameter in turn and re-folds the model, at a cost
of order 
forward evaluations per gradient for 
thawed parameters
(most built-in models register analytic gradients and avoid this;
see Appendix H). A local model can instead supply
closed-form derivatives, which makes the fast analytic-gradient path
available for fits that include it and lets the model be used with
hmc. The internals of that path are described in
Appendix H.1; this section covers only what a model
author has to provide. The derivative routines use the RealArray
interface, so they are written in C++ (the Real and
RealArray types come from xsTypes.h, as for a C++ model
function); they can, however, be attached to a model whose forward function
is written in any language.
A component may provide one or both of two routines:
- Forward-mode Jacobian (XSCCGrad).
- Returns, on the same
energy grid as the forward function, the partial derivative of the output
flux with respect to each model parameter. This is what the fast
Levenberg–Marquardt and migrad paths use.
- Reverse-mode product (XSCCGradVJP).
- Given the adjoint
of the fit statistic with respect to the output
flux, returns
in a single backward pass
without forming the full Jacobian. It is optional — a component with
only the forward routine still works, as the pipeline forms the product
itself — but it is the efficient choice for components with many
parameters or an expensive Jacobian, and it is the form hmc
prefers.
For an additive or multiplicative component the two signatures are
// Forward-mode Jacobian: dFlux_dParam[k][i] = d flux_i / d param[k]
void mymodelGradient(const RealArray& energyArray,
const RealArray& parameterValues,
int spectrumNumber,
std::vector<RealArray>& dFlux_dParam,
const std::vector<bool>& parThawed,
const std::string& initString);
// Reverse-mode VJP: dStat_dParam[k] = sum_i dStat_dFlux[i] * dflux_i/dparam[k]
void mymodelVJP(const RealArray& energyArray,
const RealArray& parameterValues,
int spectrumNumber,
const RealArray& dStat_dFlux,
RealArray& dStat_dParam,
const std::vector<bool>& parThawed,
const std::string& initString);
The forward routine resizes dFlux_dParam to one row per
parameter, each row an array the same length as the forward flux array
(energyArray.size()-1), and sets entry [k][i] to
. Differentiate only the spectral shape:
for an additive component the normalization derivative is supplied by the
pipeline. The parThawed mask (length nPar) flags the
parameters the fit will actually use — when parThawed[k] is
false the routine may leave column k at zero and skip its
work — and an empty mask means “treat every parameter as thawed”.
Convolution components instead use the XSCCConvGrad and XSCCConvGradVJP signatures, which additionally take the upstream flux (and, in the forward case, the upstream Jacobian) so the component can chain its own linear operator onto the derivatives flowing in from the components it convolves. All four typedefs are defined, with detailed comments, in Xspec/src/XSFunctions/Utilities/funcType.h.
Registering the routines. Add a grad= marker to the end of the model's first line in model.dat:
mymodel 4 0. 1.e20 C_mymodel add 0 0 grad=gv
The flag is g to register a forward Jacobian only, or gv to register both the forward Jacobian and the reverse-mode VJP. By default the routines are looked up as <model>Gradient and <model>VJP, formed from the model name (the first field of the line) — so the entry above expects mymodelGradient and mymodelVJP. If your routine names differ from that default, give them explicitly, as the built-in powerlaw does:
powerlaw 1 0. 1.e20 C_powerLaw add 0 grad=gv:powerLawGradient,powerLawVJP
(the VJP name may be omitted, in which case a trailing Gradient in the first name is replaced by VJP). The marker is accepted only on add, mul and con rows; it is not supported on mixing (mix, amx), auto-convolution (acn), or Python / user (Py_, U_) models. When initpackage finds any grad= marker in a package's model.dat it generates the gradient function-map and its registration code automatically; just rebuild the package and load it as usual.
When the analytic path is used. The analytic and finite-difference paths are combined per parameter, so a model need not be fully gradiented to benefit. For each thawed parameter whose owning component has a registered gradient, fit (Levenberg–Marquardt or migrad) builds that Jacobian column analytically; the remaining parameters — the shape parameters of an ungradiented additive or multiplicative component (its normalization stays analytic), and the master of a parameter linked across data groups — are finite-differenced one column at a time and merged in. Supplying a gradient for your component therefore speeds up its parameters even when other components in the model lack one. The whole fit reverts to finite differences only for a few structural cases: an ungradiented convolution component (a convolution must be gradiented to stay on the analytic path; acn/pileup never is), a mixing model (mix, amx), a parameter linked to more than one master (a multi-source link such as = p1 + p2), a free response/gain parameter, or the xset DISABLE_ANALYTIC_GRAD yes override.
hmc uses the same hybrid: a fully gradiented model takes the fastest reverse-mode (VJP) path, a partly gradiented one finite-differences the deferred parameters on the statistic, and only the hard structural rejects (an ungradiented convolution, or a mix/amx model) make it refuse and ask you to use chain instead. A new gradient should still be checked against the finite-difference result; xset DISABLE_ANALYTIC_GRAD yes forces the FD path for that comparison (see Appendix H.1).
Next: Third-Party Libraries In Local Up: Adding Models to XSPEC Previous: Loading a new model Contents