#ifndef GELIB_IMPL_Spline1d_h__
#define GELIB_IMPL_Spline1d_h__

#include <vector>
#include "Polynomial1d.h"
#include "../generic_classes/Containers.h"
#include "../n_dimensional/IndexedObject.h"
#include "Function1d.h"

namespace gelib_impl
{

	enum class ExtrapolationType
	{
		invalid,
		zero,
		constant,
		linear,
	};

	class Spline1d : public Function1d
	{
	private:

		///////////////////////////////////////////////////////////////////////////////////////////////
		//
		// Ensure this is called before any calls to gaussian_expectation, as it precalculates required
		// constants for an expectation calculation
		//
		///////////////////////////////////////////////////////////////////////////////////////////////
		void prepare_expectation();

		/*
		 * gradient of linear transform for underlying grid points
		 */
		double m_exp_transform_a;


		/*
		* intercept of linear transform for underlying grid points
		*/
		double m_exp_transform_b;

		/*
		 * The polynomials which make up this spline, ordered by domain
		 */
		std::vector<std::unique_ptr<Polynomial1d>> m_polynomials;

		/*
		 * The point in (finite) domain space which each polynomial begins
		 */
		std::vector<double> m_spline_points;

		/*
		 * return i where m_spline_points[i - 1] <= x < m_spline_points[i]
		 */
		size_t locate(double x);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// transform y, defined on original grid, to internal grid. Return location on internal grid
		////////////////////////////////////////////////////////////////////////////////////////////////
		size_t transform_and_locate(double* y);

		/*
		 * return the pointer to the polynomial with domain containing x
		 */
		Polynomial1d* poly_at(double x);

		/*
		 * Tolerance for letting two neighbouring polynomials be in the same domain
		 */
		double m_gridtol;

		/*
		 * gaussian_expectation_max(mu, sigma, a, b) if the location (indices in polynomial array) of a and b are known
		 * @Assumption: the variables are transformed to internal grid
		 */
		double gaussian_expectation_driver(double mu, double sigma, double a, double b, int a_location, int b_location);

		static double gaussian_expectation_max_searchroot_driver(Spline1d& s1, Spline1d& s2, double mu, double sigma, double a, double b, double tol, int* root_count, double* root, bool* p1_minus_p2_h_positive);

		static double gaussian_expectation_max_searchroot(Spline1d& s1, Spline1d& s2, double mu, double sigma, double a, double b, double tol, double* root, bool* p1_minus_p2_h_positive);

		static double gaussian_expectation_max_foundroot(Spline1d& s1, Spline1d& s2, const double* mu_t, const double* sigma_t, const double* l_t, const double* h_t, double tol, double root, bool p1_minus_p2_h_positive);

		static double gaussian_expectation_max(Spline1d& s1, Spline1d& s2, double mu, double sigma, double a, double b, double tol, double* root, bool* p1_minus_p2_h_positive);

		static void transform_grid_vector(double* x, double x1, double xn, size_t n, size_t data_gap, double* transform_a, double* transform_b);

	public:

		double get_gridtol() const;

		////////////////////////////////////////////////////////////////////////////////////
		// Check if the grid x matches the internal grid of this spline,
		// throw an exception if not
		////////////////////////////////////////////////////////////////////////////////////
		void check_matches_internal_grid(double* x, size_t data_gap, size_t n);

		////////////////////////////////////////////////////////////////////////////////////
		// modify the transform parameters for a new grid array while keeping the same fit
		////////////////////////////////////////////////////////////////////////////////////
		void update_transform(double a_transform, double b_transform);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// linear transform arbitrary grid vector length n in to integer grid points centered at 0
		////////////////////////////////////////////////////////////////////////////////////////////////
		static void transform_grid_vector(
			double* x, size_t n, size_t data_gap, double* transform_a, double* transform_b);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// transform expectation parameters so that the expectation is for original,
		// un-transformed grid points
		////////////////////////////////////////////////////////////////////////////////////////////////
		void reverse_transform_expectation(
			double* mu, double* sigma, double* l, double* h);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// transform x to it's position on the internal grid points
		////////////////////////////////////////////////////////////////////////////////////////////////
		void reverse_transform(
			double* x);

		static void reverse_transform(
			double* x, double a_transform, double b_transform);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// set coefficients under threshold = 0.0
		////////////////////////////////////////////////////////////////////////////////////////////////
		void round_coefficients(double tol, double recip_tol);

		/*
		 * number of polynomials in the spline
		 */

		size_t get_no_polynomials() const;

		/*
		 * retrieve the i'th polynomial pointer in this spline
		 */
		Polynomial1d* get_poly_pointer(size_t i) const;

		/*
		 * retrieve the i'th domain point of this spline
		 */
		double get_domain_point(size_t i) const;

		Spline1d spline_op(const Spline1d& s, std::unique_ptr<Polynomial1d>(Polynomial1d::* op) (const Polynomial1d&)) const;

		/*
		 * Add (or subtract) each polynomial of two splines creating a new spline.
		 * Note that the domains of the polynomials should match perfectly
		 */
		Spline1d operator+(const Spline1d& s) const;
		Spline1d operator-(const Spline1d& s) const;

		/*
		 * =======================Input===================================
		 *
		 * x, y: equi-length arrays of doubles containing data pairs to fit.
		 *
		 * n: length of x and y arrays
		 *
		 * order: the order of every polynomial in the spline
		 * NOTE - Required that n => order + 1
		 *
		 * boundaries ('0', 'c', 'l'): the type of values to use for y values outside of the
		 * fit domain  of provided x values
		 * '0' (default) - y values outside domain are set equal to 0
		 * 'c' - y values outside domain are set equal to the y values at the lower and
		 * upper limits of the domain
		 * 'l' - A linear (extrapolation) fit is created at both sides from the two outermost domain points
		 * to give y values outside domain
		 *
		 * =======================Output===================================
		 *
		 * n - 1 polynomials stored in the m_polynomials vector, each domain between x values
		 * TODO update new parameters to documentation
		 */
		void fit(const double* x, const double* y, size_t n,
			double transform_a, double transform_b,
			size_t order = 5, char boundaries = '0',
			size_t data_spacing_x = 1, size_t data_spacing_y = 1);

		void fit(const double* x, const double* y, size_t n,
			double transform_a, double transform_b,
			size_t order = 5, char left_extrapolation = '0', char right_extrapolation = '0',
			size_t data_spacing_x = 1, size_t data_spacing_y = 1);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// this fit will modify the array x to the internal grid
		////////////////////////////////////////////////////////////////////////////////////////////////
		void fit(double* x, const double* y, size_t n, size_t order = 5, char boundaries = '0',
			size_t data_spacing_x = 1, size_t data_spacing_y = 1);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// this fit won't modify the vector x, instead storing it's data in a temporary array
		////////////////////////////////////////////////////////////////////////////////////////////////
		void fit(const double* x, const double* y, size_t n, size_t order = 5, char boundaries = '0',
			size_t data_spacing_x = 1, size_t data_spacing_y = 1);

		////////////////////////////////////////////////////////////////////////////////////////////////
		// number of data points to left and right of each spline point using a given polynomial order
		////////////////////////////////////////////////////////////////////////////////////////////////
		static size_t left_margin(size_t data_per_fit, size_t data_per_fit_div_2, size_t n);
		static size_t right_margin(size_t left_margin, size_t data_per_fit, size_t data_per_fit_div_2, size_t n);

		/*
		 * return the value of the spline at point y on the internal grid, not necessarily the grid which
		 * the spline was fitted on.
		 */
		double evaluate(double y) override;

		/*
		 * transform x, defined on the original grid points of the fit, to equivalent on internal grid
		 * return the value of the spline at the point
		 */
		double evaluate(double* x);

		/*
		 * values of the spline over vector of x values
		 * outputting to the vector f_x
		 */
		void evaluate(d_vector const& x, d_vector& f_x);

		/*
		 * values of the spline over vector of x values
		 * outputting vector of results to some indexed container
		 */
		template <const size_t N>
		void evaluate(d_vector& x, IndexedData<N>& f_x, double* data_start, size_t data_gap)
		{
			double xi;
			d_vector f_x_vec;
			f_x_vec.resize(x.size());
			for (size_t i = 0; i < x.size(); ++i)
			{
				xi = x[i];
				f_x_vec[i] = evaluate(&xi);
			}
			f_x.insert_vector(f_x_vec, data_start, data_gap);
		}

		/*
		 * =======================Input===================================
		 *
		 * Splines s1 and s2, which are fitted from identical data at gridpoints 'domain_points'
		 * doubles mu, sigma, a, b. The orders of polynomials in s1 and s2 should be equivalent
		 *
		 * =======================Output===================================
		 *
		 * IF (s1 - s2 is monotonic over all domain_points)
		 *	returns E[max(s1(X), s2(X))]
		 * where E is expectation over interval [a, b] for a < b
		 * and X~(mu, sigma^{2})
		 *
		 * ELSE
		 *	returns NAN-value
		 *
		 */
		static void gaussian_expectation_max(d_vector& output, Spline1d& s1, Spline1d& s2, d_vector& mu, double sigma, double l, double h, double tol);

		static void gaussian_expectation_max(d_vector& output, Spline1d& s1, Spline1d& s2, d_vector const& mu, const double* sigma, const double* l, const double* h, double tol);

		//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		// expectation variables already transformed
		//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		static double gaussian_expectation_max(Spline1d& s1, Spline1d& s2, const double* mu_t, const double* sigma_t, const double* l_t, const double* h_t, double tol, double* root, bool* p1_minus_p2_h_positive);

		/*
		 * Calculate the expectation of this spline within the interval [l,h],
		 * i.e: sum {E[p(X)]} for p in m_polynomials where the domain of p is
		 * within [l, h]
		 * where X ~ N(mu, sigma^{2})
		 */
		double gaussian_expectation(double mu, double sigma, double l, double h);

		//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		// expectation variables already transformed
		//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		double gaussian_expectation(const double* mu_t, const double* sigma_t, const double* l_t, const double* h_t);

		/*
		 * locate the indices of the polynomial(s) where a and b reside in the spline. If they are in the same
		 * polynomial within the prescribed tolerance, m_grid_tol, then they are forced to the same location.
		 */
		void locate_pair(double a, double b, int* a_location, int* b_location);

		static const char UNRECOGNISED_EXTRAPOLATION;
		static bool valid_extrapolation_string(const std::string* input);
		static char parse_extrapolation_string(const std::string& input);

		//////////////////////////////////////////////////////////////////////////////////////////////
		// CONSTRUCTORS                                               
		//////////////////////////////////////////////////////////////////////////////////////////////

		//////////////////////////////////////////////////////////////////////////////////////////////
		// create a spline fit and return a unique_ptr to it, transforming the grid points x to be 
		// zero-centered
		//////////////////////////////////////////////////////////////////////////////////////////////
		static std::unique_ptr<Spline1d> Spline1d_from_fit(double* x, const double* y, size_t no_data, size_t order = 5, char boundaries = '0',
			double domain_tol = Structure::Domain_1d::ERROR_PRECISION(), size_t data_spacing_x = 1, size_t data_spacing_y = 1);

		//////////////////////////////////////////////////////////////////////////////////////////////
		// create a spline fit and return a unique_ptr, assuming that the grid points x are 
		// zero-centered using transformation parameters provided
		//////////////////////////////////////////////////////////////////////////////////////////////
		static std::unique_ptr<Spline1d> Spline1d_from_fit(const double* x, const double* y, size_t no_data,
			double transform_a, double transform_b,
			size_t order = 5, char boundaries = '0',
			double domain_tol = Structure::Domain_1d::ERROR_PRECISION(), size_t data_spacing_x = 1, size_t data_spacing_y = 1);

		static std::unique_ptr<Spline1d> Spline1d_from_fit(const double* x, const double* y, size_t no_data,
			double transform_a, double transform_b,
			size_t order, char left_extrapolation = '0', char right_extrapolation = '0',
			double domain_tol = Structure::Domain_1d::ERROR_PRECISION(), size_t data_spacing_x = 1, size_t data_spacing_y = 1);

		Spline1d(double gridtol = Structure::Domain_1d::ERROR_PRECISION());

		Spline1d(const std::vector<std::unique_ptr<Polynomial1d>>& polynomials,
			const std::vector<double>& spline_points,
			double transform_a, double transform_b,
			double gridtol = Structure::Domain_1d::ERROR_PRECISION());

		Spline1d(const Spline1d& s);

		~Spline1d();

	};

}

#endif // Spline1d_h__