#ifndef GELIB_IMPL_BASISFUNCTION_H_
#define GELIB_IMPL_BASISFUNCTION_H_

#include <vector>
#include <memory>
#include "Spline1d.h"
#include "../2_dimensional/d_matrix.h"

namespace gelib_impl
{
	class BasisSpline1d
	{
		//////////////////////////////////////////////////////////////////
		// gradient of linear transform for internal grid points
		//////////////////////////////////////////////////////////////////
		double m_transform_a;

		//////////////////////////////////////////////////////////////////
		// intercept of linear transform for internal grid points
		//////////////////////////////////////////////////////////////////
		double m_transform_b;

		//////////////////////////////////////////////////////////////////
		// one spline per grid point which collectively represent a smooth
		// function
		//////////////////////////////////////////////////////////////////
		std::vector<std::unique_ptr<Spline1d>> m_splines;

		//////////////////////////////////////////////////////////////////
		// overall grid functions are defined on
		//////////////////////////////////////////////////////////////////
		double* m_internal_grid;

		//////////////////////////////////////////////////////////////////
		// function with 1 at a single grid point and 0 elsewhere
		//////////////////////////////////////////////////////////////////
		double* m_I;

		//////////////////////////////////////////////////////////////////
		// number of gridpoints of the function this represents
		//////////////////////////////////////////////////////////////////
		size_t m_num_gridpoints;

	public:

		////////////////////////////////////////////////////////////////////////////////////
		// modify the transform parameters for a new grid array while keeping the same fit
		////////////////////////////////////////////////////////////////////////////////////
		void update_transform(double* x, size_t data_spacing_x = 1);

		void reverse_transform(
			double* x) const;

		void reverse_transform_centered(
			double* x) const;

		//////////////////////////////////////////////////////////////////
		// gaussian expectation of this basis function, outputting to vector E
		// where I ~ N (mu, sigma).
		//////////////////////////////////////////////////////////////////
		void gaussian_expectation(double mu, double sigma, double l, double h, d_vector& E) const;

		double gaussian_expectation(size_t grid_j, d_vector const& mu, const double* sigma, const double* l, const double* h, d_vector const& in) const;

		void gaussian_expectation(d_vector const& mu, double sigma, double l, double h, d_matrix const& in, d_vector& out) const;

		//////////////////////////////////////////////////////////////////
		// gaussian expectation of this basis function, outputting to matrix E
		// where I ~ N (mu_i, sigma) so each row E_i corresponds to each mean in
		// vector mu
		//////////////////////////////////////////////////////////////////
		void gaussian_expectation(d_vector const& mu, double sigma, double l, double h, d_matrix& E) const;

		// TODO description
		void fit(double* x, size_t order, char boundaries = '0',
			double domain_tol = Structure::Domain_1d::ERROR_PRECISION(), size_t data_spacing_x = 1);

		static std::unique_ptr<BasisSpline1d> BasisFunction1d_from_fit(double* x, size_t num_data, size_t order, char boundaries = '0',
			double domain_tol = Structure::Domain_1d::ERROR_PRECISION(), size_t data_spacing_x = 1);

		~BasisSpline1d();

		BasisSpline1d(size_t grid_size);

	};
}

#endif /* GELIB_IMPL_BASISFUNCTION_H_ */


