This is a wrapper class for easy use of GSL splines It is (OpenMP) thread safe as long as splines are not created by a subthread and later tried to be used by more than that thread! More...

#include <Spline.h>

Public Member Functions
	GSLSpline ()=default

	GSLSpline (std::string name)
	Construct giving the spline a name (useful for error/warning messages) More...

	GSLSpline (double x, double y, int nx, std::string splinename="NoName", const gsl_interp_type *interpoltype=SPLINE_FIDUCIAL_INTERPOL_TYPE)
	Construct a spline from pointers to x and y. Both must have nx elements. More...

	GSLSpline (const DVector &x, const DVector &y, std::string splinename="NoName", const gsl_interp_type *interpoltype=SPLINE_FIDUCIAL_INTERPOL_TYPE)
	Construct a spline from vectors x and y. Both must have the same size. More...

	GSLSpline (const GSLSpline &rhs)

GSLSpline &	operator= (const GSLSpline &rhs)

	~GSLSpline ()

	operator bool () const
	Is the spline created or not? More...

void	create (const double x, const double y, int nx, std::string splinename="NoName", const gsl_interp_type *interpoltype=SPLINE_FIDUCIAL_INTERPOL_TYPE)
	Create a spline from pointers to x and y. Both must have nx elements. More...

void	create (const DVector &x, const DVector &y, std::string splinename="NoName", const gsl_interp_type *interpoltype=SPLINE_FIDUCIAL_INTERPOL_TYPE)
	Create a spline from vectors x and y. Both must have the same size. More...

double	operator() (double x) const
	Overload of the () operator for easy evaluation of the spline. More...

double	eval (double x) const
	Get the value of the spline (if out of bounds we use the closest value) More...

double	eval_deriv (double x, int deriv) const

double	deriv_x (double x) const
	Get the value of the first derivative of the spline. More...

double	deriv_xx (double x) const
	Get the value of the second derivative of the spline. More...

std::pair< double, double >	get_xrange () const
	Get the range the spline was created on. More...

std::string	get_name () const
	Get the name of the spline. More...

void	set_out_of_bounds_warning (bool v)
	Turn on/off warnings if we try to evaluate out of bounds (we use closest value in that case) More...

DVector	get_x_data ()
	Get the raw GSL data for x used to create the spline. More...

DVector	get_y_data ()
	Get the raw GSL data for y=f(x) used to create the spline. More...

void	free ()
	Free up memory associated with the spline. More...

Detailed Description

This is a wrapper class for easy use of GSL splines It is (OpenMP) thread safe as long as splines are not created by a subthread and later tried to be used by more than that thread!

The fiducial splines are cubic splines for 1D and bicubic for 2D assuming natural (y''=0) boundary conditions at the end-points

The fiducial choice allows to get the function value and up to second order derivatives. NB: if you change the fiducial choice to a lower order spline, e.g to gsl_interp_linear, then the second derivative deriv_xx etc. will not work

If evaluating out of bounds it will return the closest value. To get a warning when this happens call set_out_of_bounds_warning(true) before using the spline or set the define SPLINE_FIDUCIAL_SPLINE_WARNING to be true

Errors are handled via the throw_error function

Compile time defines:

USE_MPI : Use MPI (only effect is in on how errors are handled)

USE_OMP : Use OpenMP

SPLINE_FIDUCIAL_INTERPOL_TYPE : Cubic spline is the fiducial choice

SPLINE_FIDUCIAL_INTERPOL_TYPE_2D : Bicubic is the fiducial choice

SPLINE_FIDUCIAL_SPLINE_WARNING : Show warnings if we evaluate out of bounds

Definition at line 80 of file Spline.h.

Constructor & Destructor Documentation

◆ GSLSpline() [1/5]

FML::INTERPOLATION::SPLINE::GSLSpline::GSLSpline ( )

default

◆ GSLSpline() [2/5]

FML::INTERPOLATION::SPLINE::GSLSpline::GSLSpline ( std::string name )

Construct giving the spline a name (useful for error/warning messages)

Definition at line 37 of file Spline.cpp.

37: GSLSpline() { this->name = name; }

FML::INTERPOLATION::SPLINE::GSLSpline::GSLSpline

GSLSpline()=default

◆ GSLSpline() [3/5]

FML::INTERPOLATION::SPLINE::GSLSpline::GSLSpline	(	double *	x,
		double *	y,
		int	nx,
		std::string	splinename = `"NoName"`,
		const gsl_interp_type *	interpoltype = `SPLINE_FIDUCIAL_INTERPOL_TYPE`
	)

Construct a spline from pointers to x and y. Both must have nx elements.

Definition at line 39 of file Spline.cpp.

                                                                       {
                create(x, y, nx, splinename, interpoltype);
            }

◆ GSLSpline() [4/5]

FML::INTERPOLATION::SPLINE::GSLSpline::GSLSpline	(	const DVector &	x,
		const DVector &	y,
		std::string	splinename = `"NoName"`,
		const gsl_interp_type *	interpoltype = `SPLINE_FIDUCIAL_INTERPOL_TYPE`
	)

Construct a spline from vectors x and y. Both must have the same size.

Definition at line 47 of file Spline.cpp.

                                                                       {
                // Check that dimensions match
                if(x.size() != y.size())
                  throw std::runtime_error("GSLSpline::create The input vectors x, y don't have incompatible lengths!");
                create(x, y, splinename, interpoltype);
            }

◆ GSLSpline() [5/5]

FML::INTERPOLATION::SPLINE::GSLSpline::GSLSpline ( const GSLSpline & rhs )

Definition at line 75 of file Spline.cpp.

                                                      {
                // We just create the spline from scratch instead of copying all the data
                if (rhs.spline != nullptr) {
                    create(rhs.spline->x, rhs.spline->y, rhs.size_x, rhs.name, rhs.interpoltype_used);
                } else {
                    free();
                }
            }

◆ ~GSLSpline()

FML::INTERPOLATION::SPLINE::GSLSpline::~GSLSpline ( )

Definition at line 284 of file Spline.cpp.

284{ free(); }

Member Function Documentation

◆ create() [1/2]

void FML::INTERPOLATION::SPLINE::GSLSpline::create	(	const double *	x,
		const double *	y,
		int	nx,
		std::string	splinename = `"NoName"`,
		const gsl_interp_type *	interpoltype = `SPLINE_FIDUCIAL_INTERPOL_TYPE`
	)

Create a spline from pointers to x and y. Both must have nx elements.

Definition at line 89 of file Spline.cpp.

                                                                         {
                // Clean up if we already have a spline allocated
                if (spline != nullptr) {
                    free();
                }
 
                // Check if array is increasing
                int sign_one = x[1] > x[0] ? 1.0 : -1.0;
                for (int i = 1; i < nx; i++) {
                    if ((x[i] - x[i - 1]) * sign_one <= 0) {
                        std::string errormessage =
                            "[GSLSpline::create] x-array for spline " + name + " is not monotone\n";
                        throw_error(errormessage);
                    }
                }
 
                // If array is reversed then reverse it since GSL
                // requires a monotonely increasing array
                DVector xx, yy;
                if (sign_one < 0) {
                    xx = DVector(nx);
                    yy = DVector(nx);
                    for (int i = 0; i < nx; i++) {
                        xx[i] = x[nx - 1 - i];
                        yy[i] = y[nx - 1 - i];
                    }
                    x = xx.data();
                    y = yy.data();
                }
 
                // Set class variables
                xmin = x[0];
                xmax = x[nx - 1];
                dx_min = (x[1] - x[0]) / 2.0;
                dx_max = (x[nx - 1] - x[nx - 2]) / 2.0;
                size_x = nx;
                name = splinename;
                interpoltype_used = interpoltype;
 
                // Make the spline
                spline = gsl_spline_alloc(interpoltype, nx);
                gsl_spline_init(spline, x, y, nx);
 
                // Make accelerators (one per thread if OpenMP)
                // If nthreads = 1 we are likely trying to create a spline 
                // inside a OMP region. In that case allocate as many as we have
                // threads availiable on the system
#ifdef USE_OMP
                int nthreads = 1;
#pragma omp parallel
                {
                    int id = omp_get_thread_num();
                    if (id == 0)
                        nthreads = omp_get_num_threads();
                }
                if(nthreads == 1)
                  nthreads = nmax_threads;
 
                // If we make a spline inside a nested openmp then
                // the stuff above does not work
                // nthreads = std::max(16, nthreads);
                xaccs = std::vector<gsl_interp_accel *>(nthreads);
                for (auto & xa : xaccs) {
                    xa = gsl_interp_accel_alloc();
                }
#else
                xacc = gsl_interp_accel_alloc();
#endif
            }

◆ create() [2/2]

void FML::INTERPOLATION::SPLINE::GSLSpline::create	(	const DVector &	x,
		const DVector &	y,
		std::string	splinename = `"NoName"`,
		const gsl_interp_type *	interpoltype = `SPLINE_FIDUCIAL_INTERPOL_TYPE`
	)

Create a spline from vectors x and y. Both must have the same size.

Definition at line 163 of file Spline.cpp.

                                                                         {
                if (x.size() != y.size()) {
                    std::string errormessage =
                        "[GSLSpline::create] x and y array must have the same number of elements for spline " + name +
                        "\n";
                    throw_error(errormessage);
                }
                create(x.data(), y.data(), int(x.size()), splinename, interpoltype);
            }

◆ deriv_x()

double FML::INTERPOLATION::SPLINE::GSLSpline::deriv_x ( double x ) const

Get the value of the first derivative of the spline.

Definition at line 275 of file Spline.cpp.

275{ return eval_deriv(x, 1); }

FML::INTERPOLATION::SPLINE::GSLSpline::eval_deriv

double eval_deriv(double x, int deriv) const

Definition: Spline.cpp:203

◆ deriv_xx()

double FML::INTERPOLATION::SPLINE::GSLSpline::deriv_xx ( double x ) const

Get the value of the second derivative of the spline.

Definition at line 276 of file Spline.cpp.

276{ return eval_deriv(x, 2); }

◆ eval()

double FML::INTERPOLATION::SPLINE::GSLSpline::eval ( double x ) const

Get the value of the spline (if out of bounds we use the closest value)

Definition at line 183 of file Spline.cpp.

                                                 {
                if (spline == nullptr) {
                    std::string errormessage = "[GSLSpline::eval] Spline " + name + " has not been created!\n";
                    throw_error(errormessage);
                }
 
                // If out of bounds show a warning and set x to boundary value
                out_of_bounds_check(x);
                x = std::max(x, xmin);
                x = std::min(x, xmax);
 
                // Return f, f' or f'' depending on value of deriv
#ifdef USE_OMP
                gsl_interp_accel * xacc_thread = xaccs[omp_get_thread_num()];
#else
                gsl_interp_accel * xacc_thread = xacc;
#endif
                return gsl_spline_eval(spline, x, xacc_thread);
            }

◆ eval_deriv()

double FML::INTERPOLATION::SPLINE::GSLSpline::eval_deriv	(	double	x,
		int	deriv
	)		const

Definition at line 203 of file Spline.cpp.

                                                                  {
                if (spline == nullptr) {
                    std::string errormessage = "[GSLSpline::eval_deriv] Spline " + name + " has not been created!\n";
                    throw_error(errormessage);
                }
                if (deriv < 0 || deriv > 2) {
                    std::string errormessage = "[GSLSpline::eval_deriv] Got deriv = " + std::to_string(deriv) +
                                               " for spline " + name + " Expected 0 = f, 1 = f' or 2 = f''\n";
                    throw_error(errormessage);
                }
 
                // If out of bounds show a warning and set x to boundary value
                out_of_bounds_check(x);
                x = std::max(x, xmin);
                x = std::min(x, xmax);
 
                // Return f, f' or f'' depending on value of deriv
#ifdef USE_OMP
                gsl_interp_accel * xacc_thread = xaccs[omp_get_thread_num()];
#else
                gsl_interp_accel * xacc_thread = xacc;
#endif
 
                double dydx = 0.0;
                if (deriv == 0) {
                    dydx = gsl_spline_eval(spline, x, xacc_thread);
                } else if (deriv == 1) {
                    dydx = gsl_spline_eval_deriv(spline, x, xacc_thread);
                } else if (deriv == 2) {
                    dydx = gsl_spline_eval_deriv2(spline, x, xacc_thread);
                }
                return dydx;
            }

◆ free()

void FML::INTERPOLATION::SPLINE::GSLSpline::free ( )

Free up memory associated with the spline.

Definition at line 240 of file Spline.cpp.

                                 {
                if (spline != nullptr) {
                    // Reset class variables
                    xmin = xmax = 0.0;
                    size_x = 0;
                    dx_min = dx_max = 0.0;
 
                    // Free the spline
                    gsl_spline_free(spline);
                    spline = nullptr;
 
                    // Free accelerators
#ifdef USE_OMP
                    for (auto & xa : xaccs) {
                        gsl_interp_accel_free(xa);
                    }
                    std::vector<gsl_interp_accel *>().swap(xaccs);
#else
                    gsl_interp_accel_free(xacc);
                    xacc = nullptr;
#endif
                }
            }

◆ get_name()

std::string FML::INTERPOLATION::SPLINE::GSLSpline::get_name ( ) const

Get the name of the spline.

Definition at line 278 of file Spline.cpp.

278{ return name; }

◆ get_x_data()

DVector FML::INTERPOLATION::SPLINE::GSLSpline::get_x_data ( )

inline

Get the raw GSL data for x used to create the spline.

Definition at line 165 of file Spline.h.

165{ return DVector(spline->x, spline->x + spline->size); }

◆ get_xrange()

std::pair< double, double > FML::INTERPOLATION::SPLINE::GSLSpline::get_xrange ( ) const

Get the range the spline was created on.

Definition at line 277 of file Spline.cpp.

277{ return {xmin, xmax}; }

◆ get_y_data()

DVector FML::INTERPOLATION::SPLINE::GSLSpline::get_y_data ( )

inline

Get the raw GSL data for y=f(x) used to create the spline.

Definition at line 167 of file Spline.h.

167{ return DVector(spline->y, spline->y + spline->size); }

◆ operator bool()

FML::INTERPOLATION::SPLINE::GSLSpline::operator bool ( ) const

explicit

Is the spline created or not?

Definition at line 84 of file Spline.cpp.

84{ return (spline != nullptr); }

◆ operator()()

double FML::INTERPOLATION::SPLINE::GSLSpline::operator() ( double x ) const

Overload of the () operator for easy evaluation of the spline.

Definition at line 181 of file Spline.cpp.

181{ return eval(x); }

FML::INTERPOLATION::SPLINE::GSLSpline::eval

double eval(double x) const

Get the value of the spline (if out of bounds we use the closest value)

Definition: Spline.cpp:183

◆ operator=()

GSLSpline & FML::INTERPOLATION::SPLINE::GSLSpline::operator= ( const GSLSpline & rhs )

Definition at line 60 of file Spline.cpp.

                                                                  {
                if (this != &rhs) {
                    // We just create the spline from scratch instead of copying all the data
                    if (rhs.spline != nullptr) {
                        create(rhs.spline->x, rhs.spline->y, rhs.size_x, rhs.name, rhs.interpoltype_used);
                    } else {
                        free();
                    }
                }
                return *this;
            }

◆ set_out_of_bounds_warning()

void FML::INTERPOLATION::SPLINE::GSLSpline::set_out_of_bounds_warning ( bool v )

Turn on/off warnings if we try to evaluate out of bounds (we use closest value in that case)

Definition at line 279 of file Spline.cpp.

279{ out_of_bounds_warning = v; }

The documentation for this class was generated from the following files:

/Users/hansw/local/FML/FML/Spline/Spline.h
/Users/hansw/local/FML/FML/Spline/Spline.cpp

Public Member Functions