Dear fellows,
I have a question concerning the gsl (spline) interpolation, especially
about the important 'initialization' process.
As I understand it, there are two possibilities: first, the
'gsl_interp_init' command and the higher-level-interface
'gsl_spline_init'.
According to the documentation, there is an 'interpolation object'
initialized, but - please forgive my rather naive question - where is this
object? Can it be accessed? In C++ I learned that I end up at least with
pointers to a (newly created) object.
My question is based on the attempt to parallelize a code, especially to
give each thread a copy of an interpolated data set. Thus my assumption
is, that I have to write a copy constructor for my interpolation class
(see appended) file, where each time the 'gsl_spline_ini' command is
executed?
Thanks for any help and have a pleasant weekend,
-- Klaus Huthmacher.
/* interpolation.cpp */
#include "interpolation.h"
/* default constructor */
interpolation::interpolation():
_x(), _y(),
_x_min(0.0), _x_max(0.0),
_acc(nullptr,gsl_interp_accel_free),
_spline(nullptr,gsl_spline_free)
{}
/* constructor */
interpolation::interpolation(std::vector<double> x, std::vector<double> y):
_x( x ), _y( y ),
_x_min( _x.front() ), _x_max( _x.back() ),
_acc(gsl_interp_accel_alloc(),gsl_interp_accel_free),
_spline(gsl_spline_alloc( gsl_interp_cspline, _x.size()),gsl_spline_free)
{
gsl_spline_init( _spline.get(), _x.data(), _y.data(), _x.size() );
}
/* set */
auto interpolation::set(std::vector<double> x, std::vector<double> y) -> void
{
_x = x;
_y = y;
_x_min = _x.front();
_x_max = _x.back();
_acc = std::unique_ptr<gsl_interp_accel,void(*)(gsl_interp_accel*)>
(gsl_interp_accel_alloc(),gsl_interp_accel_free);
_spline = std::unique_ptr<gsl_spline,void(*)(gsl_spline*)>
(gsl_spline_alloc( gsl_interp_cspline, _x.size()),gsl_spline_free);
gsl_spline_init( _spline.get(), _x.data(), _y.data(), _x.size() );
}
/* operator */
auto interpolation::operator()(double x) -> double
{
if( x < _x.front() )
return _y.front();
if( x > _x.back() )
return _y.back();
return gsl_spline_eval( _spline.get(), x , _acc.get() );
}/* interpolation.h */
#ifndef INTERPOLATION_H_
#define INTERPOLATION_H_
#include <iostream>
#include <gsl/gsl_spline.h>
#include <memory>
#include <vector>
class interpolation
{
public:
interpolation();
interpolation(std::vector<double>, std::vector<double>);
auto set(std::vector<double>, std::vector<double>) -> void;
auto operator()(double x) -> double;
auto min()-> double { return _x_min; }
auto max()-> double { return _x_max; }
private:
std::vector<double> _x, _y;
double _x_min, _x_max;
std::unique_ptr<gsl_interp_accel,void(*)(gsl_interp_accel*)> _acc;
std::unique_ptr<gsl_spline,void(*)(gsl_spline*)> _spline;
};
#endif
/* End of File */
