ccgsl  2.7.0
C++wrappersforGnuScientificLibrary
ccgsl Documentation

The ccgsl package is designed to provide simple C++ wrappers for the GNU Scientific Library.It grew out of some practical classes designed to make GSL a little easier to use without changing its basic features. It is not meant to be an object oriented version of GSL. Rather it tries to keep most of the efficiency of GSL while letting a C++ compiler simplify the creation and destruction of GSL structs. It also allows better interaction with the C++ standard template library and optionally changes errors into C++ exceptions.

This documentation should be read in conjuction with the GSL manual to work out what functions do. Most functions here are just wrappers for GSL functions designed to work without C++ users using malloc() and free() functions or pointers of any kind. So, typically, instead of passing an array as double*, for example, you can pass a gsl_vector or often even a std::vector<double> or std::array<double,N>.

Much of the package is created using Ruby scripts to rewrite the GSL header files.

Here are some of the features of ccgsl.

  • The only installed files are header files. This means that no extra libraries need to be installed to get ccgsl to work. Just link the gsl libraries as usual.
  • Structs that have _alloc() and _free() functions are replaced with classes that behave like shared pointers or Java classes. They are lightweight shared handles that should be created as automatic variables. When one of these classes goes out of scope, it deletes the struct it points to automatically unless it is shared by another class.
  • The ccgsl classes can be used safely and efficiently inside STL containers. So, for example, std::list<gsl::vector> and std::vector<gsl::rng> work properly.
  • The vector and block classes are STL containers. So you can use the STL algorithms with them.
  • Most GSL functions are reimplemented in ccgsl namespaces and have nearly the same form as the corresponding GSL function. So, for example, gsl_sf_gamma() becomes gsl::sf::gamma() with the same parameters. Replace pointers to GSL structs with references to ccgsl classes.
  • You can mix GSL functions freely with ccgsl. The get() function of a ccgsl class gives a GSL pointer. So, for example gsl::vector::get() returns a gsl_vector* object.
  • ccgsl can turn on exceptions instead of GSL errors. Thus you can use the simpler form of many functions within a try block instead of testing the result of a function call with an extra parameter. See the gsl::exception class for details.
  • GSL uses pointers to structs (usually with void* parameters) to pass functions to other functions. C++ users more usually use function objects: that is, objects (for the parameters) with methods for the functions. So ccgsl derives subclasses of the GSL structs and allow you to construct these with functions or function objects. These are normally passed by reference in ccgsl.
  • In C++ it is usually safer to use std::vector<double> std::array<double,N> or gsl::vector than double[]. So ccgsl uses template functions where GSL uses pointers to arrays. Usually the template parameter can be replaced with a pointer to an array, because ccgsl provides an overloaded function. Most of these overloaded functions are not given explicitly in the documentation to avoid clutter.

Example

The following example shows an example of some of the features of ccgsl. It is designed to create a function fitSigma_u() that fits data to a distribution. The distribution is the sum of a normal and the negative of a half-normal distribution. So it is a good place to apply GSL.

Think of the following code fragments as joined together in a file called fitSigma_u.cc. Then to compile with gcc on a Linux system, one would use

g++ -o fitSigma_u fitSigma_u.cc -lgsl -lgslcblas -lm

First, we need the headers to tell us what parts of ccgsl we will use. We also include vector so that we can use std::vector<double>, which we can mix easily with ccgsl.

The next block of code defines the density function of a simplified version of the distribution we want. Note here that we use functions from cdf.hpp and randist.hpp to access normal density and distrution functions in GSL.

class SFAhpdf {
public:
SFAhpdf(double const lambda=0) : lambda {lambda} {}
double function(double x) const {
double const a {gsl::cdf::ugaussian_P(-lambda*x)};
double const b {gsl::ran::ugaussian_pdf(x)};
return 2*a*b;
}
void setLambda( double const lambda ){
this->lambda = lambda;
}
private:
double lambda;
};
double ugaussian_P(double const x)
C++ version of gsl_cdf_ugaussian_P().
Definition: cdf.hpp:35
double ugaussian_pdf(double const x)
C++ version of gsl_ran_ugaussian_pdf().
Definition: randist.hpp:551
double b(int order, double qq)
C++ version of gsl_sf_mathieu_b().
Definition: sf_mathieu.hpp:298
double a(int order, double qq)
C++ version of gsl_sf_mathieu_a().
Definition: sf_mathieu.hpp:272

The third block of code defines a distribution function. Like SFApdf we define it as a class so that we can supply a parameter lambda.

The main things to note here are the following.

  1. The distribution function is obtained by numerical integration using gsl::integration::qagil from integration.hpp.
  2. gsl::integration::qagil requires a gsl::function::scl. This is the ccgsl equivalent of a gsl_function (in fact it is a derived class) and can be constructed directly from an SFApdf object using pdf {sfahpdf,&SFAhpdf::function} in the constructor.
  3. gsl::integration::qagil also requires a workspace, which we create in the constructor. Note that there is no need to free the workspace once we are finished.
class SFAhdist {
public:
SFAhdist(double const lambda=0, double const epsabs=1e-6, double const epsrel=1e-6,
size_t N=250)
: lambda {lambda}, sfahpdf {lambda}, pdf {sfahpdf,&SFAhpdf::function},
N {N}, workspace{N}, epsabs{epsabs}, epsrel{epsrel} {}
double f(double const x) const {
double result {0};
double abserr {0};
gsl::integration::qagil(pdf,x,epsabs,epsrel,N,workspace,result,abserr);
if(abserr > epsabs)
std::cout << "warning: abserr = " << abserr << "." << std::endl;
return result;
}
double df(double const x) const {
return sfahpdf.function(x);
}
std::pair<double,double> fdf(double const x) const {
return std::make_pair(f(x),df(x));
}
void setLambda( double const lambda ){
this->lambda = lambda;
sfahpdf.setLambda(lambda);
}
private:
double lambda;
mutable SFAhpdf sfahpdf;
mutable gsl::function_scl pdf;
size_t const N;
mutable gsl::integration::workspace workspace;
double const epsabs;
double const epsrel;
};
Class that extends gsl_function so that it can be constructed from arbitrary function objects.
Workspace for integration.
Definition: integration.hpp:37
int qagil(function_scl &f, double const b, double const epsabs, double const epsrel, size_t const limit, workspace &workspace, double &result, double &abserr)
C++ version of gsl_integration_qagil().
int df(double const h, gsl_multifit_nlinear_fdtype const fdtype, vector const &x, vector const &wts, gsl::multifit::nlinear::function_fdf &fdf, vector const &f, matrix &J, vector &work)
C++ version of gsl_multifit_nlinear_df().
gsl_multilarge_nlinear_fdf fdf
Typedef for shorthand.
gsl_sf_result result
Typedef for gsl_sf_result.
Definition: sf_result.hpp:30

The fourth block of code contains no ccgsl functions. It changes the distribution by using a different parameterisation.

class F : public SFAhdist {
public:
F(double const sigma_u=0, double const epsabs=1e-6, double const epsrel=1e-6, size_t N=250)
: SFAhdist {sigma_u,epsabs,epsrel,N}, r{std::sqrt(1+sigma_u*sigma_u)} {}
void setSigma_u( double const sigma_u ){
setLambda(sigma_u);
r = std::sqrt(1+sigma_u*sigma_u);
}
double operator()(double const z, double const x_i) const {
return f(z/(x_i*r));
}
private:
// Keep as sqrt(1+sigma_u^2)
double r;
};
complex sqrt(complex const &z)
C++ version of gsl_complex_sqrt().
Definition: complex.hpp:933
double F(double phi, double k, mode_t mode)
C++ version of gsl_sf_ellint_F().
Definition: sf_ellint.hpp:138

We fit the distribution F in the previous code block by minimising an Anderson–Darling test statistic. The only ccgsl code of note here is q[i] = F::operator()(e[i],x[i]);. gsl::vector is the equivalent of gsl_vector but does not require and alloc or free functions and we can access members using [] as an alternative to gsl::vector::get().

We pass the gsl::vector arguments by reference. But they are sufficiently small objects that we can also pass by value. In either case the elements are only available by reference. If you want a copy of a gsl::vector you have to make it explicitly, for example by using gsl::vector::clone().

class AndersonDarling : public F {
public:
AndersonDarling(double const sigma_u=0, double const epsabs=1e-6, double const epsrel=1e-6,
size_t N=250)
: F {sigma_u,epsabs,epsrel,N} {}
// Notice we still have setSigma().
double operator()(gsl::vector const& e, gsl::vector const& x) const {
size_t const N {e.size()};
if(x.size() != N)
throw std::length_error("AndersonDarling::operator(): e and x must have same size.");
q.resize(N); // We’ll store temporary results in q
for(size_t i {0}; i < N; ++i){
q[i] = F::operator()(e[i],x[i]);
}
// We need q in ascending order.
std::sort(q.begin(),q.end());
// Start computing the A–D test statistic
double result {0};
for(size_t i {0}; i < N; ++i){
result += (2*i+1)*std::log(q[i]);
result += (2*(N-i)+1)*std::log(1-q[i]);
}
result *= -1.0/N;
return result - N; // N is not really needed for optimisation but is correct.
}
private:
mutable std::vector<double> q;
};
This class handles vector objects as shared handles.
Definition: vector.hpp:74
size_type size() const
The size (number of elements) of the vector.
Definition: vector.hpp:1169
complex log(complex const &a)
C++ version of gsl_complex_log().
Definition: complex.hpp:970

The next block of code defines anothe function object. This contains no ccgsl functions, but it is constructed so that we can create a gsl::function_scl object from an object of it.

class SFATestFunction {
public:
SFATestFunction(gsl::vector const& e, gsl::vector const& x)
: e {e}, x{x} {
size_t const N {e.size()};
if(x.size() != N)
throw std::length_error("SFATestFunction::SFATestFunction(): "
"e and x must have same size.");
}
double function(double const sigma_u){
AD.setSigma_u(sigma_u);
return AD(e,x);
}
private:
gsl::vector const& e;
gsl::vector const& x;
// We will need an AD test
AndersonDarling AD;
};

Now we come to the main function. It finds the best fit of the distribution with some weights provided by the gsl::vector x.

Here are some things to note.

  1. We use gsl::exception::enable() to turn GSL errors into exceptions, which we can catch in the code.
  2. We use gsl::fminimizer to call gsl_fminimizer.
  3. gsl::fminimizer needs a function. We create it using function_scl testFunction {sfaTestFunction,&SFATestFunction::function};.
  4. We call member functions of gsl::fminimizer, though you could use, for example, gsl::fminimizer::iterate(fminimizer) instead of fminimizer.iterate() if you like your code to look more like gsl_fminimizer_iterate(fminimizer).
double
fitSigma_u(gsl::vector const& e, gsl::vector const&x, double const sigma_u_upper){
double constexpr epsabs {1e-4}; // we don’t need extreme accuracy
try {
if(sigma_u_upper <= 0)
throw std::domain_error("fitSigma_u(): sigma_u_upper must be positive.");
// Use GSL minimisation
using namespace gsl;
// Create a minimiser
min::fminimizer fminimizer {min::fminimizer::quad_golden()};
// Create a function
SFATestFunction sfaTestFunction {e,x};
function_scl testFunction {sfaTestFunction,&SFATestFunction::function};
// Get initial guess
double const guess {std::min(1.0,sigma_u_upper)};
// Set the minimiser
double x_lower {0};
double x_upper {sigma_u_upper};
fminimizer.set(testFunction,guess,0,sigma_u_upper);
// Now start iterating
while(x_upper-x_lower > epsabs){
auto result = fminimizer.iterate();
if(GSL_EBADFUNC == result)
throw exception("fitSigma_u(): bad function", __FILE__, __LINE__, result);
if(GSL_FAILURE == result)
throw exception("fitSigma_u(): Could not minimise", __FILE__, __LINE__, result);
x_lower = fminimizer.x_lower();
x_upper = fminimizer.x_upper();
}
return fminimizer.x_minimum();
} catch(gsl::exception& e){
std::cerr << e.get_reason() << std::endl;
throw;
} catch(...) {
std::cerr << "Unknown error in C++ code." << std::endl;
throw;
}
}
This class is used to handle gsl exceptions so that gsl can use these rather than the GSL error handl...
Definition: exception.hpp:387
char const * get_reason() const
Get the message explaining the reason for the error/exception.
Definition: exception.hpp:416
static void enable()
Set the handler function to handle exceptions.
Definition: exception.hpp:452
static type const * quad_golden()
Static type.
Definition: min.hpp:488
int min(movstat::end_t const endtype, vector const &x, vector &y, workspace &w)
C++ version of gsl_movstat_min().
Definition: movstat.hpp:581
The gsl package creates an interface to the GNU Scientific Library for C++.
Definition: blas.hpp:34

The last block of code shows a test of the code. The only additional ccgsl feature of note here is that we can allocate ccgsl vector objects using brace-enclosed lists.

int main(){
gsl::vector x {2.617764,3.969121,1.157941,3.64723,1.873826,1.91901,1.183213,3.441218,
2.050355,1.469476,2.534231,3.520795,3.51455,1.186722,2.789286,3.081629,
2.241102,2.164384,3.84475,2.038049,3.417877,3.53866,2.472255,3.553603,
3.262197,2.115331,2.647081,3.544415,2.35675,3.766609,2.676864,1.589442,
1.272145,2.932182,2.880038,2.721431,3.652606,1.213367,3.819442,1.696094};
gsl::vector e {1.781801,0.2975302,0.506332,-3.589495,0.6598619,-2.724952,0.5975643,
-1.268469,0.02219264,-0.4164864,0.7725084,1.368086,-4.855619,0.05086869,
-3.530058,2.622537,-0.7611229,-2.796356,1.881148,0.4252468,-2.94855,
-7.740664,-12.56288,-5.385373,-11.5698,-1.322548,-5.1494,-9.592563,
-1.927037,-0.1411757,-0.1148876,-0.9549714,-1.549748,-1.546938,-7.231412,
-0.07847876,-3.261605,-0.003874675,-4.874783,-5.254548};
// test
SFATestFunction sfaTestFunction {e,x};
// estimate fit
double const sigma_u {fitSigma_u(e,x,10)};
std::cout << "Estimated sigma_u = " << sigma_u << std::endl;
exit(0);
}