Last modified: 13 Jul 2010

URL: https://cxc.cfa.harvard.edu/sherpa/threads/cusermodels/

Sherpa User Models in C

Sherpa Threads (CIAO 4.12 Sherpa v1)


Overview

Synopsis:

This thread explains how to write your own Sherpa parameterized models in C and include them in a Sherpa fitting session.

Last Update: 13 Jul 2010 - updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.


Contents


Setup

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.


Writing your functions

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);

Compiling the user library

We now make the user library by

python setup.py cfiles=multi.c,multijcm.c libs=jcmnutils,m libdirs=/home/jcm/src/linux/lib 
install  --prefix=`pwd`

Here we pass the names of our files using the cfiles argument; setup.py is the standard file provided with the distribution. Note the backwards-quotes around pwd, which is a Unix trick to evaluate the pwd command on the fly—if you prefer you can just explicitly do

python setup.py cfiles=multi.c,multijcm.c install  --prefix=/home/jcm/test/usermodel

or whatever your working directory is.

If successful, this command will create a 'lib/python2.7/site-packages' subdirectory in the directory pointed to by 'prefix=', and a file 'userfuncs.so' inside it. The command also creates a 'build' subdirectory in your current directory; this subdirectory may be deleted. (You should delete this directory if it already exists, before running the command). To check that the userfuncs library has been created, type

unix> ls lib/*/*/*.so
lib/python2.7/site-packages/userfuncs.so

Using your user model

Now we are ready to fit your data in Sherpa using the new models! Here is an example Sherpa session. You can either type the following lines interactively, or put them in a file user.py and type 'sherpa user.py'.

sherpa> import userfuncs
sherpa> load_data('jcm.data')
sherpa> plot_data()

sherpa> load_user_model(userfuncs.flat, "myf")
sherpa> add_user_pars("myf", ["a","p2","p3","p4"], [1,0,0,0], [0,0,0,0])
sherpa> set_model(myf)

sherpa> fit()
sherpa> plot_fit()

sherpa> load_user_model(userfuncs.conic, "myconic")
sherpa> add_user_pars("myconic", ["a","e","p3","p4"], [1,0,0,0], [0,0,0,0])
sherpa> set_model(myconic)

sherpa> fit()
sherpa> plot_fit()

In this session, the command 'import userfuncs' looks for the userfuncs.so file created by setup.py. This makes the models userfunc.flat and userfuncs.conic, which you can use with the load_user_model function to create Sherpa functions with the names you prefer—in this case, myf and myconic. We then give names and initial values to the parameters with add_user_pars, and then we are ready—we can use set_model, fit, and plot_fit as for a normal Sherpa session. The output of the second fit in this case is Figure 1.

Figure 1: Fit plot

[Sherpa plot of conic function fit]
[Print media version: Sherpa plot of conic function fit]

Figure 1: Fit plot

The fit of the conic user model to some test data, jcm.data


Summary

In this thread we have shown you how you can fit your own C-based function to 1-dimensional data.


History

15 Dec 2008 original version
13 Jul 2010 updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.