#include <iostream>
#include <memory>
#include <vector>

#include "gsl/gsl_integration.h"
#include "gsl/gsl_math.h"
#include "gsl/gsl_spline.h"

class data_set
{
	public:
		data_set(std::vector<double>, std::vector<double>);
		std::vector<double> _x, _y;
		std::unique_ptr<gsl_interp_accel,void(*)(gsl_interp_accel*)> _acc;
		std::unique_ptr<gsl_spline,void(*)(gsl_spline*)> _spline;
};

data_set::data_set(std::vector<double> x, std::vector<double> y):
_x(x), _y(y),
_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());
}

class gsl_params
{
	public:
		gsl_params(std::unique_ptr<gsl_interp_accel,void(*)(gsl_interp_accel*)>,
			std::unique_ptr<gsl_spline,void(*)(gsl_spline*)>);
		std::unique_ptr<gsl_interp_accel,void(*)(gsl_interp_accel*)> _acc;
		std::unique_ptr<gsl_spline,void(*)(gsl_spline*)> _spline;
};

gsl_params::gsl_params(
		std::unique_ptr<gsl_interp_accel,void(*)(gsl_interp_accel*)> acc,
		std::unique_ptr<gsl_spline,void(*)(gsl_spline*)> spline
		):
		_acc(std::move(acc)), _spline(std::move(spline))
		{}

auto gsl_kernel(double x, void* p) -> double
{
	std::unique_ptr<gsl_params> params(static_cast<gsl_params*>(p));
	return gsl_spline_eval(params->_spline.get(), x, params->_acc.get());
}

auto gsl_integral(gsl_params* my_params, double lower, double upper) -> double
{
	unsigned int size = 1000;
	double result(0.0), error(0.0);

	std::unique_ptr<gsl_integration_workspace,void(*)(gsl_integration_workspace*)>
		w(gsl_integration_workspace_alloc(size),gsl_integration_workspace_free);

	gsl_function f{ gsl_kernel, my_params };

	gsl_integration_qags(&f, lower, upper, 0.0, 1.0e-7, size, w.get(), &result, &error);

	return result;
}

int main()
{
	std::vector<double> x, y;
	
	for(unsigned int i = 0; i <= 100; ++i)
	{
		double x_value = static_cast<double>(i * 0.1);
		x.push_back(x_value);
		y.push_back(x_value * x_value);
	}

	data_set data(std::move(x), std::move(y));

	gsl_params params(std::move(data._acc), std::move(data._spline));

	std::cout << gsl_integral(&params, 0.0, 3.0) << std::endl;

	return 0;
}