#ifndef GELIB_IMPL_BASISFUNCTION3d_H_
#define GELIB_IMPL_BASISFUNCTION3d_H_

#include "../n_dimensional/IndexSet.h"
#include "../1_dimensional/BasisFunction.h"
#include "../1_dimensional/Spline1d.h"
#include "Gaussian3d.h"
#include "d_cube.h"

namespace gelib_impl
{
	class BasisSpline3d
	{
		///////////////////////////////////////////////////////////////////////////////
		// The working arrays for 3d expectation
		///////////////////////////////////////////////////////////////////////////////
		std::vector<d_cube> m_input_data;	
		std::vector<d_cube> m_output_data;

		///////////////////////////////////////////////////////////////////////////////
		// Wrapper for the indices of each spline to output to some form 
		// of n-dimensional storage.
		///////////////////////////////////////////////////////////////////////////////
		IndexSet<3> m_indices;

		///////////////////////////////////////////////////////////////////////////////
		// The basis functions in each dimension indexed by m_indices
		///////////////////////////////////////////////////////////////////////////////
		std::vector<std::unique_ptr<BasisSpline1d>> m_basis_fits;

		//////////////////////////////////////////////////////////////////
		// True if respective basis fits need to be re-fit
		//////////////////////////////////////////////////////////////////
		std::vector<bool> m_needs_fitting;

		//////////////////////////////////////////////////////////////////
		// Original grids in each dimension
		//////////////////////////////////////////////////////////////////
		std::vector<d_vector> m_grids;

		//////////////////////////////////////////////////////////////////
		// helper function for gaussian_expectation_sum_3d
		//////////////////////////////////////////////////////////////////
		static void prepare_expectation_sum_3d(
			std::vector<d_vector>& grid_data,
			std::vector<const d_cube*>& input_data,
			std::array<size_t, 3> ordering,
			d_cube& out_1, d_cube& out_2, d_cube& out_3);

		void construct_indices(
			std::vector<d_vector>& grid_data,
			std::vector<const IndexedData<3>*>& input_data, std::vector<IndexedData<3>*>& output_data,
			std::array<size_t, 3> ordering);

		void update_grids(std::vector<d_vector> const& grid_data);


	public:

		////////////////////////////////////////////////////////////////////////////////////
		// modify the transform parameters for new grid arrays while keeping the same fits
		////////////////////////////////////////////////////////////////////////////////////
		void update_transform();

		////////////////////////////////////////////////////////////////////////////////////
		// Constructors for single integration
		////////////////////////////////////////////////////////////////////////////////////
		BasisSpline3d(size_t varied_dimension, d_vector& grid, std::vector<const d_cube*>& input_data, d_cube& output_data);

		BasisSpline3d(size_t varied_dimension, d_vector& grid, std::vector<const d_cube*>& input_data, d_cube& output_data, std::array<size_t, 3> ordering);

		////////////////////////////////////////////////////////////////////////////////////
		// Constructor for triple integration
		////////////////////////////////////////////////////////////////////////////////////
		BasisSpline3d(std::vector<d_vector>& grid_data, std::vector<const d_cube*> const& input_data,
			d_cube const & output_data, std::array<size_t, 3> ordering);

		BasisSpline3d(std::vector<d_vector>& grid_data, const d_cube& B_0, const d_cube& means1, const d_cube& means2, const d_cube& means3,
			d_cube const& output_data, std::array<size_t, 3> ordering);

		size_t num_indices(size_t which);

		std::array<int, 3>& get_varied_indices(size_t which, size_t i);

		std::unique_ptr<BasisSpline1d>& get_basis_grid(size_t which);

		void gaussian_expectation(double sigma, double l, double h, size_t which_fit, size_t dimension, d_cube& E_X);

		void fit_input(size_t dimension, size_t order, char boundaries, double domain_tol, size_t which_fit);

		static void fit_and_expectation(
			size_t which_fit, size_t dimension, double sigma, double l, double h, d_cube& output,
			BasisSpline3d& b3d, size_t order, char extrapolation, double domain_tol);

		//////////////////////////////////////////////////////////////////
		// Original grids in each dimension
		//////////////////////////////////////////////////////////////////
		static void gaussian_expectation_sum_3d(
			d_cube& output,
			std::vector<d_vector>& grid_data,
			std::vector<const d_cube*>& input_data,
			double sigma1, double sigma2, double sigma3,
			double l1, double h1, double l2, double h2, double l3, double h3,
			size_t order, const std::string& extrapolation);

		static void gaussian_expectation_sum_3d(
			d_cube& output,
			std::vector<d_vector>& grid_data,
			std::vector<const d_cube*>& input_data,
			double sigma1, double sigma2, double sigma3,
			double l1, double h1, double l2, double h2, double l3, double h3,
			size_t order, const std::string& extrapolation, const std::array<size_t, 3>& ordering);

		//////////////////////////////////////////////////////////////////////////////////////////
		// If the grids are the same size as previous can preserve the basis fits for efficiency
		/////////////////////////////////////////////////////////////////////////////////////////
		static void gaussian_expectation_sum_3d(
			BasisSpline3d& b3d, d_cube& output,
			std::vector<d_vector>& grid_data,
			std::vector<const d_cube*>& input_data,
			double sigma1, double sigma2, double sigma3,
			double l1, double h1, double l2, double h2, double l3, double h3,
			size_t order, const std::string& extrapolation);

		static void gaussian_expectation_sum_1d(
			d_cube& output,
			size_t dimension,
			d_vector& grid_data,
			std::vector<const d_cube*>& input_data,
			double sigma1,
			double l1, double h1,
			size_t order, char extrapolation, double domain_tol, const std::array<size_t, 3>& ordering);

		static void gaussian_expectation_sum_1d(
			d_cube& output,
			size_t dimension,
			d_vector& grid_data,
			std::vector<const d_cube*>& input_data,
			double sigma1,
			double l1, double h1,
			size_t order, char extrapolation, double domain_tol);

		static void gaussian_expectation_sum_3d(
			BasisSpline3d& b3d, d_cube& output, std::vector<d_vector>& grid_data,
			const d_cube& B_0, const d_cube& means1, const d_cube& means2, const d_cube& means3,
			double sigma1, double sigma2, double sigma3,
			double l1, double h1, double l2, double h2, double l3, double h3, 
			size_t order, const std::string& extrapolation);

};
}
#endif