GSL C++ Wrappers

The GNU Scientific Library is wicked. It has well-written C code and pretty good documentation for almost every numerical function I've needed so far. However, using raw C functions amongst oject-oriented C++ code can be difficult, ugly or repetitive. I write OO code because I think in terms of objects when I'm programming. So I wrote a few classes to help my translate, and I've used all of them many times in my GMF simulation.

Maximization


aon_maximize.h
aon_maximize.cpp

Usage:

First, define the mathematical function that you wish to maximize (if there are multiple arguments, the one to be maximized goes in 'Arg', the others are passed using a ParamStruct. 

	double myFunction(double Arg, void* Params)
	{
		// get the parameters from the void-pointer
		myParamStruct myParams = *(myParamStruct *) Params;
		double y = myParams.param1;
		double z = myParams.param2;
		double x = Arg;
	
		// do your function, whatever it is
		return x*y*y/exp(-z);
	}

Then in your code, create a ParamStruct if necessary and setup the maximization. 

	myParamStruct myParams;
	myParams.param1 = 0.5; myParams.param2 = 77.3;
      
	Maximization myMax;
	if (myMax.setFunction (&myFunction, &myParams) )
		log_error();

Now, you can easily get the maximum of the function: the arguments to FindMaximum are Guess, LowerLimit and UpperLimit. 

	cout << "maximum: " << myMax.FindMaximum(25, 0, 50) << endl;

Interpolation


aon_interpolation.h
aon_interpolation.cpp

Usage:

Define a derived class and over-ride the abstract BuildSpline() function to produce the function you wish to interpolate. Pass any extra arguments as a list<double>. 

   MyInterpClass : public Interpolation
   {
   public:
      MyInterpClass (const int NumPoints, 
                     const gsl_interp_type *InterpType);
      virtual int BuildSpline (const double& LowerBound, 
                               const double& UpperBound,
                               list& ArgList);
   };

The constructor is dead simple, choose the right interpretation type and the number of points you want. ( I use gsl_interp_cspline). 

   MyInterpClass::MyInterpClass(const int NumPoints,
                                const gsl_interp_type *InterpType)
   :  Interpolation (NumPoints, InterpType)
   {
      itsSpline = gsl_spline_alloc (InterpType, NumPoints);
      itsAccel = gsl_interp_accel_alloc ();
   }

The complicated part is building the spline. Here you have to understand how GSL interpolations work - the object orientation is in the way that you access the interpolation, rather than how you build it.

Create the over-riding BuildSpline function, in which we calculate the values for the spline, then use the gsl function to set it up. 

   int MyInterpClass::BuildSpline(const double& LowerBound,
                                  const double& UpperBound,
                                  list& ArgList)
   {
      // Get any arguments you may need from ArgList
      double z = ArgList.pop_front();

      // create the arrays that will hold our argument/value pairs
      double *xa = new double[itsNumPoints];
      double *ya = new double[itsNumPoints];
	
      // Fill the arrays using the function y(x, [z1, z2...]) that 
      // you wish to interpolate
      double Arg = LowerBound;
      for (int i = 0; i < itsNumPoints; i++)
      {
         xa[i] = Arg;
         ya[i] = some_expensive_operation (x, z);
         Arg += (UpperBound - LowerBound)/ (itsNumPoints + 0.0);
      }
      // call the GSL function that sets up the interpolation
      Status = gsl_spline_init (itsSpline, xa, ya, itsNumPoints);	
 
      delete []xa;
      delete []ya;
      if (Status)
         return EXIT_FAILURE;

      // set the status of the interpolation so that we can use it.
      itsInitializedFlag = 1;
   
      return 0;
   }

OK, so we spent a long time writing the interpolation - now it is really easy to use:

Create the interpolation 

   // create interpolation.
   MyInterpClass MyInterp(1000, gsl_interp_cspline);

   // call BuildSpline
   list MyArgList;
   MyArgList.push_back(z);
   if ( MyInterp.BuildSpline(LowerBound, UpperBound, MyArgList) )
   	log_error();

And to access the interpolation: 

   if ( x < LowerBound )
   {
      return SomeApproxOfFunction(x);
   }
   if ( x > UpperBound )
   {
      return AnotherApproxOfFunction(x);
   }
   else
   {
      return MyInterpolation.Evaluate(x);
   }

Bessel Functions


gmf_bessel.h
gmf_bessel.cpp

Usage:

This one is really easy - although it does use my debug functions. 
   double K = BesselK(Nu, x);