Sherpa User Models in C

Install user model tar package

Download the usermodel.tar.gz file and unpack it in a clean working directory.

mkdir my_working_dir
cd my_working_dir
wget "http://hea-www.harvard.edu/~jcm/usermodel.tar.gz"
tar xvzf usermodel.tar.gz

Set up environment variables

The C user model enviroment makes use of five environment variables:

• PYTHONPATH lets Sherpa find your new user function library at run time.
• CC and CXX give the path of the C and C++ compilers (gcc and g++).
• CFLAGS lets you define the location of include files needed by your user functions.
• LDFLAGS and CXXFLAGS are needed if you are compiling a 32-bit program on an x86_64 machine.

unix% setenv PYTHONPATH ${PWD}/lib/python2.7/site-packages
unix% setenv CC /usr/bin/gcc
unix% setenv CXX /usr/bin/g++
unix% setenv CFLAGS "-m32 -I/home/jcm/src/linux/inc -I/home/jcm/src/linux/inc32"
unix% setenv LDFLAGS "-m32"
unix% setenv CXXFLAGS "-m32"

In this case, we assume that our user function library will be found in ./lib/python2.7/site-packages. When running CIAO tools and applications, wrapper scripts append the location of the CIAO python packages to PYTHONPATH on the fly; setting up CIAO does not set the PYTHONPATH in your normal environment. If you've set PYTHONPATH yourself, you'll need to do 'setenv PYTHONPATH ${PYTHONPATH}:${PWD}/lib/python2.7/site-packages' instead here.

The -m32 is only needed for 32-bit compilation on x86_64.

We will be using the external library /home/jcm/src/linux/lib/libjcmnutils.a as well as the system libm.a and the python2.7 library in our CIAO installation. The CFLAGS and LDFLAGS settings reflect this.

The actual user functions

In this example, we imagine we want to fit one of several simple parameterized functions to some data points; no instrument model is involved. We create one more .c files containing the functions. In this case, I create a single source file multi.c:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "jcmutils.h"


double flatfunc( double* p, double x )
{
 double r;
 r = p[0];
 printf( "CALLED FUNC1 R=%f\n", r );
    return (r-6378.0);
}

double conicfunc( double* p, double x )
{
 double r;
 r = p[0]*(1-p[1]*p[1])/(1 + p[1]*utn_dcosd(x));
 printf( "CALLED FUNC2 R=%f\n", r );
    return (r-6378.0);
}

This file contains two trivial functions: flatfunc and conicfunc, each of which takes a vector of parameters p, a parameter value x, and returns a function value. An external function utn_dcosd is referenced; its prototype is in the external include file jcmutils.h (which the compiler will find in the path defined by CFLAGS) and its definition is in the library libjcmnutils.a which we will reference below in a call to the script setup.py.

Wrapper functions

Sherpa's actual user function interface expects you to return a whole vector of function values, one for each x-axis point in the model. So, we write a wrapper file to call the individual functions in a loop. (You may find it more efficient just to evaluate the vector of values directly, and omit the need for a separate multi.c file—some find the separated version easier to understand and maintain).

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
double flatfunc( double* p, double x );
double conicfunc( double* p, double x );
#include "userproto.h"

int flat( double* p, double* x, int xsize, double* y, int ysize )
{
 long i;
 for ( i = 0; i < xsize; i++ )
 {
  y[i] = flatfunc( p, x[i] );
 }
 return EXIT_SUCCESS;
}

int conic( double* p, double* x, int xsize, double* y, int ysize )
{
 long i;
 for ( i = 0; i < xsize; i++ )
 {
  y[i] = conicfunc( p, x[i] );
 }
 return EXIT_SUCCESS;
}

// Integral versions - dummy
int flat_int( double* p, double* x1, int xsize1, double* x2, int xsize2, double* y, int ysize )
{ 
 return EXIT_FAILURE;
}

int conic_int( double* p, double* x1, int xsize1, double* x2, int xsize2, double* y, int ysize )
{ 
 return EXIT_FAILURE;
}

Here we have TWO wrappers for each user function, one for functions which take single x values as input and return the function value at that x value, and one for functions which return integrals of the function values across bins (x1,x2). We are using the non-integral versions, but the other ones still need to be there as dummies.

The double-valued functions, flat and conic, operate on an array of values x of size xsize and return an array of values y of size ysize, which is always equal to xsize.

To keep the compiler happy, we also include function prototypes for both the wrappers and the functions they call.

Create the pyproto.h and userproto.h files

For each wrapper function foo with n parameters, we add a line USERMODELFCT1D( foo, n ), to the pyproto.h file. This example uses

USERMODELFCT1D( flat, 4 ),
USERMODELFCT1D( conic, 4 ),

In fact, flat only has 1 parameter, but I find it useful to define extra unused parameters so that all my functions have the same number of parameters—that makes it easier to switch which model I'm using in the final Sherpa script.

We also need to make the userproto.h file containing a C prototype for each function.

  int flat(const double* p, const double* x, const int xsize,
              double* y, const int ysize);
  int conic(const double* p, const double* x, const int xsize,
              double* y, const int ysize);
  int flat_int(const double* p, const double* xlo, const int xlosize,
                  const double* xhi, const int xhisize,
                  double* y, const int ysize);
  int conic_int(const double* p, const double* xlo, const int xlosize,
                  const double* xhi, const int xhisize,
                  double* y, const int ysize);