#pragma once

#include "BasisFunction3d.h"

namespace gelib
{
	using d_cube = gelib_impl::d_cube;
	using d_matrix = gelib_impl::d_matrix;
	using d_vector = gelib_impl::d_vector;

	/// <summary>
	/// Let \f$ b_{1}, b_{2}, b_{3} \f$ be the dimension lengths of \f$ B^{in} \f$. Fix two dimensions,
	/// say without loss of generalisation these are 2 and 3, leaving a vector varying in dimension 1.
	/// Let \f$ v_{i}, w_{i} \f$ be each of these vectors for \f$ B^{in} \f$ and \f$ M \f$ respectively,
	/// where \f$ 1 \leq i \leq b_{1} \f$.
	/// Now generate a matrix \f$ E \f$ of expectations \f$ \mathbb{E}[\mathcal{I}_{j}(X)] \f$ for \f$ X \sim N(\mu, \sigma^{2}) \f$
	/// where each \f$ \mu \f$ is an element of \f$ w_{i} \f$ and corresponds to a row of \f$ E \f$,
	/// and the j'th column of \f$ E \f$ corresponds to basis function \f$  \mathcal{I}_{j}(X)  \f$ 
	/// valued 1 at the j'th gridpoint of input grid \f$ x \f$ and 0 elsewhere.
	/// The output cube comprises of \f$ b_{2}b_{3} \f$ vectors, one for each i, given by the vector-matrix product \f$ v_{i}E \f$.
	/// Its final dimensions are \f$ n, b_{2}, b_{3} \f$ where \f$ n \f$ is the length of \f$ x \f$ 
	/// </summary>
	/// <param name="output">The output cube.</param>
	/// <param name="x">The grid of which to calculate expectations at each point.</param>
	/// <param name="B_in">\f$ B^{in} \f$ cube.</param>
	/// <param name="means">Means cube \f$ M \f$, identical dimensions to \f$ B^{in} \f$ unless they are length 1.</param>
	/// <param name="standard_deviation">The standard deviation of all random variables \f$ \sigma \f$.</param>
	/// <param name="l">Integration lower limit.</param>
	/// <param name="h">Integration upper limit.</param>
	/// <param name="dimension">The dimension to perform integration Either 1, 2, or 3.</param>
	/// <param name="order">Order of polynomials used in fitting.</param>
	/// <param name="extrapolation">Extrapolation type, one of "zero", "constant", or "linear".</param>
	/// <param name="check_input">if set to <c>true</c> the function will throw exceptions on invalid input.</param>
	void basis_expectation_sum_3d_single(
		d_cube& output,
		d_vector& x, d_cube const& B_in, d_cube const& means,
		double standard_deviation, double l, double h,
		size_t dimension, size_t order, const std::string& extrapolation, bool check_inputs);

	const std::string basis_function_expectation_sum_3d_single_check(
		d_vector const& x, d_cube const& B_in, d_cube const& means,
		double standard_deviation, double l, double h,
		size_t dimension, size_t order, const std::string& extrapolation);

	/// <summary>
	/// Three successive calls to basis_expectation_sum_3d_single, one over each of three dimensions x, y and z.
	/// Initially, \f$ B^{0} \f$ is used as the \f$ B^{in} \f$ parameter, while the successive two calls use the output 
	/// from the previous call as \f$ B^{in} \f$ (see basis_expectation_sum_3d_single). The \f$ M^{i} \f$ are used as 
	/// inputs in the order they are numbered.
	/// x, y and z dimensions correspond to rows, columns and layers in the input cubes respectively.
	/// Dimensional ordering of integration is provided by arguments.
	/// </summary>
	/// <param name="output">The output cube which will have the same dimensions as the grid data.</param>
	/// <param name="grid_data">Array of three grids to calculate expectations on, in the order x,y and z.</param>
	/// <param name="B_0_C">The cube \f$ B^{0} \f$.</param>
	/// <param name="means_1">First means cube \f$ M^{0} \f$, identical dimensions to \f$ B^{0} \f$.</param>
	/// <param name="means_2">Second means matrix \f$ M^{1} \f$, shares two dimension lengths with \f$ B^{0} \f$.</param>
	/// <param name="means_3">Third means vector \f$ M^{2} \f$, the same length of one dimension of \f$ B^{0} \f$.</param>
	/// <param name="sigma1">Standard deviation of x.</param>
	/// <param name="sigma2">Standard deviation of y.</param>
	/// <param name="sigma3">Standard deviation of z.</param>
	/// <param name="l1">The lower limit of integration for x.</param>
	/// <param name="h1">The upper limit of integration for x.</param>
	/// <param name="l2">The lower limit of integration for y.</param>
	/// <param name="h2">The upper limit of integration for y.</param>
	/// <param name="l3">The lower limit of integration for z.</param>
	/// <param name="h3">The upper limit of integration for z.</param>
	/// <param name="order">Order of polynomials used in fitting.</param>
	/// <param name="extrapolation">Extrapolation type, one of "zero", "constant", or "linear".</param>
	/// <param name="check_inputs">if set to <c>true</c> the function will throw exceptions on invalid input.</param>
	/// <param name="P_1">First dimension to integrate (1 = x, 2 = y, or 3 = z).</param>
	/// <param name="P_2">Second dimension in integration (1 = x, 2 = y, or 3 = z).</param>
	/// <param name="P_3">Third dimension in integration (1 = x, 2 = y, or 3 = z).</param>
	void basis_expectation_sum_3d_triple(d_cube& output, std::vector<d_vector>& grid_data, d_cube& B_0_C, d_cube& means_1, d_cube& means_2, d_cube& means_3, 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, bool check_inputs, size_t P_1, size_t P_2, size_t P_3);

	const std::string basis_function_expectation_sum_3d_triple_check(
		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);
}