Description

This document describes how to use the S-Lang programming language to create a user-defined Sherpa model. Such a model can then be fit to your data in exaclty the same manner as any of the pre-defined models in Sherpa.

Creating a S-lang based model in Sherpa is a two step process. First, a S-lang function with specific inputs must be compiled, and then a special Sherpa/S-lang function must be called to register the model with Sherpa. Once created, a S-lang model is entirely equivalent to any built-in model in terms of usability. A simple example is shown below, for a power law model:

define slang_pow() {

  variable p, norm, Emin, Emax, dE, Result;

  if ( _NARGS == 3 ) {
     (p, norm, Emin) = ();
     Result = norm*(Emin^(-p));

  } else if ( _NARGS == 4 ) {
     (p, norm, Emin, Emax) = ();
     dE = Emax - Emin;
     Result = norm*(Emin^(-p))*dE;
  }

  Result = typecast(Result, _typeof(Emin));
  return Result;
}

and the Sherpa/S-lang command to register the model after it has been compiled with the evalfile() command is:

  sherpa> () = register_model("slang_pow",["power","norm"], 1,
                              [1.0,1.e-2],  % default values
                              [-10,1.e-20], % Minimum values
                              [10,1.e5],    % Maximum values
                              [1,1]);       % Both thawed by default

S-lang User Model Function Parameters

The S-lang function takes a variable number of arguments, depending upon the number of model parameters, and returns a vector of the evaluated model on the S-lang stack. The input arguments are:

Integrated and Differential Models

One complication in writing a S-lang usermodel for Sherpa is that Sherpa models may be integrated (for use with binned data) or not (for unbinned or differential data). Sherpa, by default, creates all S-lang usermodels with integrate off, although this can be changed after the model is created using the <modelname> integrate on command. If integrate is off, the model is expected to caluculate a value at each input point X. If integrate is on, two vectors, Xmin and Xmax, will be passed to the function. In this case, the function should return the value of the model integrated over the each bin. The accuracy of the integration is up to the model author. Depending on the model and the input data, a simple evalution at one bin edge (or the center of the bin) and multiplication by the bin width (as shown in the example above) or an analytic integral may be needed. Note that the units output from the model are different in each case; if the model is not integrated, the units have a "per X-axis unit length", i.e. per keV or per Angstrom, while in the binned case they do not. Since the user can switch between wavelength and energy units at will in Sherpa (via analysis wavelength or analysis energy), the model writer must handle all these cases to have a complete model.

Dealing with binned data

A common mistake when using S-lang usermodels with PHA-type data is that PHA data are inherently binned, so the model have integrate set to "on". The code must integrated between the lower and upper edge of each detector channel to get the total photons in the channel. If this is not done, Sherpa will still calculate a result, with parameters that frequently will be close to the correct values except for a very wrong normalization.

Dealing with unbinned data

Sherpa can also use unbinned data, such as a list of X and Y points which has no inherent binning. Is is possible that the data may be binned, but the user wants to use the model with "integrate off." In which case the normalization will be incorrect, since the model will have different units than the data. In either case, Sherpa passes a single vector X to the function. In this case, the model should simply evaluate the function at each X point and return the resulting vector. Note that in the case of binned data with integrate off, the bin minimum (not the bin center) is passed to the function; this is how Sherpa treats all such cases.

Coding details

Since Sherpa allows any model to be set to integrate off, any S-lang usermodel must be able to deal with either binned or unbinned data. The simplest method is to check how many values have been passed to the routine; for a model with M parameters, there will be either M+1 (in the unbinned case) or M+2 (for binned data). If the M+1 values are passed, the routine should return the model evaluated at the requested points, and if M+2 values are given, the code should return the model integrated over each low to high bin.

Units

As noted above, the units of the independent axis (or axes) are vital. The routine get_axes(datasetNum) can be used to retrieve a structure containing the units for dataset datasetNum. When a model is evaluated in a fit, the currently active dataset number can be determined with the Sherpa routine get_dataset_number(). In the case of basic spectral models, the input units with be either keV (if analysis energy is set) or Angstroms (if analysis wave is set). The expected output in these cases is photons/cm**2/s in each bin (for binned data) or photons/cm**2/s/keV or photons/cm**2/s/A in the unbinned case. For more complex models, the values of the independent axes should be checked carefully.

S-lang types

One Sherpa requirement is that the datatype of the returned result be the same as the input vector. This can be enforced by using the following command before returning the result via

Result = typecast(Result, _typeof(Emin));

where Emin the input vector and Result the model's calculated values. This S-lang command simply forces the returned value to have the same type as the input.

Compiling the user model

The model can be compiled with the evalfile command, a S-Lang function which loads the given file into the S-Lang interpreter. To read in the S-Lang code stored in the file myusermodel.sl, the syntax is:

  sherpa> () = evalfile("./myusermodel.sl")

The initial '() =' is used to ensure that the stack is properly cleared after execution, avoiding unwanted messages printed to the screen.

Registering the new model with Sherpa

After the new model code is compiled, it must be registered with Sherpa using the register_model function. The syntax is:

  () = register_model("new model name",
                      ["parameter1","parameter2",...,"parmeterN",],
                      ModelDimension,
                      [p1_default,p2_default,...,pN_default], 
                      [p1_minimum,p2_minimum,...,pN_minimum], 
                      [p1_maximum,p2_maximum,...,pN_maximum], 
                      [p1_thawed, p2_thawed, ...,pN_thawed]);

for a total of seven arguments. The first argument is a S-lang string which is the compiled function's name and will become the name of the new model (for example some of Sherpa's built-in model names are ngauss, xsapec, and delta2d). The second argument is an N-element string vector giving the names of each of the N parameters, such as "pos", "ampl" or "kT". The third argument gives the dimensionality of the model; spectral models are all one-dimensional, but a spatial model (such as the built-in model gauss2d) could have two. The next three arguments are all N-element vectors giving the default, minimum, and maximum values of each model parameter. The last argument is a N-element integer vector which gives the default status (either thawed or frozen) for each model parameter; either 1 for thawed or 0 for frozen is allowed. These last four arguments are optional; if they are set to NULL values, standard values with be used.