#pragma once

#include "Spline3d.h"

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//--------------------------------------------------------------------------------------------------------------------
// LIBRARY FUNCTIONS 
//--------------------------------------------------------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

namespace gelib
{

	using d_cube = gelib_impl::d_cube;
	using d_matrix = gelib_impl::d_matrix;
	using d_vector = gelib_impl::d_vector;

	/// <summary>
	/// Suppose we have a 3 dimensional random vector \f$ x_{1} = \mu(x_{0}) + \Delta \f$ where \f$ \Delta \sim N(0, C) \f$ For some covariance matrix C.
	/// \f$ \mu(x_{0}) \f$ is given for each dimension on a series of points in each dimension, by inputs a_hat_i. 
	/// This function computes \f$ E[\f$max\f$(f_{1}(x_{1}), f_{2}(x_{1}))|x_{0}] \f$ for each combination of points from a_hat_i,
	/// where \f$ f_{1}  \f$ and  \f$ f_{2} \f$ are arbitrary smooth functions defined on a series of three-dimensional grid points. 
	/// The ordering of integration is specified by the parameters \f$ P_{1}, P_{2} \f$ and \f$ P_{3} \f$, for instance, with 
	/// \f$ P_{1} = 3, P_{2} = 2, P_{3} = 1 \f$ the directions of integration will be over layers -> columns -> rows.
	/// </summary>
	/// <param name="A_out">Output cube.</param>
	/// <param name="b_1">Grid points (rows).</param>
	/// <param name="b_2">Grid points (columns).</param>
	/// <param name="b_3">Grid points (layers).</param>
	/// <param name="A_in_1">Function 1 represented on grid points.</param>
	/// <param name="A_in_2">Function 2 represented on grid points.</param>
	/// <param name="a_hat_1">Means for first variable.</param>
	/// <param name="a_hat_2">Means for Second variable.</param>
	/// <param name="a_hat_3">Means for Third variable.</param>
	/// <param name="C">The covariance matrix.</param>
	/// <param name="l1">The lower limit of integration (first integral).</param>
	/// <param name="h1">The upper limit of integration (first integral).</param>
	/// <param name="l2">The lower limit of integration (second integral).</param>
	/// <param name="h2">The upper limit of integration (second integral).</param>
	/// <param name="l3">The lower limit of integration (third integral).</param>
	/// <param name="h3">The upper limit of integration (third integral).</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="P_1">First dimension in integration.</param>
	/// <param name="P_2">Second dimension in integration.</param>
	/// <param name="P_3">Third dimension in integration.</param>
	/// <param name="check_input">if set to <c>true</c> the function will throw exceptions on invalid input.</param>
	void correlated_gaussian_expectation_max_3d(
		d_cube& A_out,
		d_vector const& b_1, d_vector const& b_2, d_vector const& b_3,
		d_cube& A_in_1, d_cube& A_in_2,
		d_vector& a_hat_1, d_vector& a_hat_2, d_vector& a_hat_3,
		d_matrix& C,
		double l1, double h1, double l2, double h2, double l3, double h3,
		size_t order, std::string& extrapolation,
		size_t P_1, size_t P_2, size_t P_3,
		bool check_input = true);

	const std::string correlated_gaussian_expectation_max_3d_check(
		d_cube& A_out,
		d_vector const& b_1, d_vector const& b_2, d_vector const& b_3,
		d_cube& A_in_1, d_cube& A_in_2,
		d_vector& a_hat_1, d_vector& a_hat_2, d_vector& a_hat_3,
		d_matrix& C,
		double l1, double h1, double l2, double h2, double l3, double h3,
		size_t order, std::string& extrapolation,
		size_t P_1, size_t P_2, size_t P_3
	);

	/// <summary>
	/// Suppose we have three dimensional random vector \f$ x_{1} = \mu(x_{0}) + \Delta \f$ where \f$ \Delta \sim N(0, C) \f$ For some covariance matrix C.
	/// \f$ \mu(x_{0}) \f$ is given for each dimension on a series of points in each dimension, by inputs a_hat_i. 
	/// This function computes \f$ E[f(x_{1})|x_{0}] \f$ for each possible combination of a_hat_i points,
	/// where \f$ f \f$ is an arbitrary smooth function defined on a series of three-dimensional grid points. 
	/// The ordering of integration is specified by the parameters \f$ P_{1}, P_{2} \f$ and \f$ P_{3} \f$, for instance, with 
	/// \f$ P_{1} = 3, P_{2} = 2, P_{3} = 1 \f$ the directions of integration will be over layers -> columns -> rows.
	/// </summary>
	/// <param name="A_out">Output cube.</param>
	/// <param name="b_1">Grid points (rows).</param>
	/// <param name="b_2">Grid points (columns).</param>
	/// <param name="b_3">Grid points (layers).</param>
	/// <param name="A_in">Function represented on grid point.</param>
	/// <param name="a_hat_1">Means for first variable.</param>
	/// <param name="a_hat_2">Means for Second variable.</param>
	/// <param name="a_hat_3">Means for Third variable.</param>
	/// <param name="C">The covariance matrix.</param>
	/// <param name="l1">The lower limit of integration (first integral).</param>
	/// <param name="h1">The upper limit of integration (first integral).</param>
	/// <param name="l2">The lower limit of integration (second integral).</param>
	/// <param name="h2">The upper limit of integration (second integral).</param>
	/// <param name="l3">The lower limit of integration (third integral).</param>
	/// <param name="h3">The upper limit of integration (third integral).</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="P_1">First dimension to integrate (1, 2, or 3).</param>
	/// <param name="P_2">Second dimension in integration (1, 2, or 3).</param>
	/// <param name="P_3">Third dimension in integration (1, 2, or 3).</param>
	/// <param name="check_input">if set to <c>true</c> the function will throw exceptions on invalid input.</param>
	void correlated_gaussian_expectation_3d(
		d_cube& A_out,
		d_vector const& b_1, d_vector const& b_2, d_vector const& b_3,
		d_cube& A_in,
		d_vector& a_hat_1, d_vector& a_hat_2, d_vector& a_hat_3,
		d_matrix& C,
		double l1, double h1, double l2, double h2, double l3, double h3,
		size_t order, std::string& extrapolation,
		size_t P_1, size_t P_2, size_t P_3,
		bool check_input = true);

	const std::string correlated_gaussian_expectation_3d_check(
		d_cube& A_out,
		d_vector const& b_1, d_vector const& b_2, d_vector const& b_3,
		d_cube& A_in,
		d_vector& a_hat_1, d_vector& a_hat_2, d_vector& a_hat_3,
		d_matrix& C,
		double l1, double h1, double l2, double h2, double l3, double h3,
		size_t order, std::string& extrapolation,
		size_t P_1, size_t P_2, size_t P_3);

	/// <summary>
	/// Let \f$ f(x, y, z) \f$ be a 3-dimensional function defined on 3 sets of gridpoints by a cube object
	/// with dimensions of \f$ n_{1} \times n_{2} \times n_{3} \f$.
	/// \f$ f \f$ can be varied in one dimension to produce a piecewise polynomial \f$ p_{ij} \f$ for each unique combination
	/// of points in the other two fixed dimensions \f$ i \f$ and \f$ j \f$.
	///
	/// A mean vector is given for each dimension \f$ \mu_{i} \f$, which have lengths \f$ m_{1},  m_{2} \f$ and \f$ m_{3} \f$,
	/// where \f$ m_{i} = n_{i} \f$ with the exception of the varied dimension, along with a length-3 weights vector \f$ w \f$.
	///
	/// Denote the varied dimension as \f$ v_{1} \f$ and fixed dimensions by \f$  v_{2},  v_{3}  \f$ where  \f$ v_{i} \in \{x=1, y=2, z=3\}\f$
	/// The weighted mean for each fixed i-j combination is 
	/// \f$ \mu_{ijk} = w_{v_{1}}\mu_{v_{1} k} + w_{v_{2}}\mu_{v_{2} i} + w_{v_{3}}\mu_{v_{3} j} \f$.
	/// The output is a cube with dimensions of \f$ m_{v_{1}} \times n_{v_{2}} \times n_{v_{3}} \f$, populated with vectors of expectations.
	/// That is, for each piecewise polynomial on the fixed gridpoints \f$ p_{ij} \f$, calculate \f$ m_{v_{1}} \f$ expectations,
	/// such that each output entry is \f$ \mathbb{E}_{k}[p_{ij}(K)] \f$ where \f$ K \sim N(\mu_{ijk}, \sigma) \f$
	/// </summary>
	/// <param name="output">The output cube.</param>
	/// <param name="x">The grid points to calculate expectations on. Must be equivalent to one dimension of function cube.</param>
	/// <param name="F_X">The points of the function defined at points on a 3-dimensional grid.</param>
	/// <param name="no_layers">The number of layers (matrices) in the function sube.</param>
	/// <param name="mu_1">Means for first argument of function (rows).</param>
	/// <param name="mu_2">Means for second argument of function (columns).</param>
	/// <param name="mu_3">Means for third argument of function (layers).</param>
	/// <param name="standard_deviation">The standard deviation.</param>
	/// <param name="l">The lower limit of integration.</param>
	/// <param name="h">The upper limit of integration.</param>
	/// <param name="weights">Length 3 weights vector</param>
	/// <param name="dimension">The dimension to calculate expectations on.</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 weighted_gaussian_expectation_3d(d_cube& output, d_vector const& x, d_cube const& F_X, size_t no_layers,
		d_vector const& mu_1, d_vector const& mu_2, d_vector const& mu_3,
		double standard_deviation, double l, double h, d_vector const& weights, size_t dimension, size_t order, const std::string& extrapolation, bool check_input = true);

	const std::string weighted_gaussian_expectation_3d_check(
		d_vector const& x, d_cube const& F_X,
		size_t no_layers,
		d_vector const& mu_1, d_vector const& mu_2, d_vector const& mu_3,
		double standard_deviation,
		double l, double h,
		d_vector const& weights,
		size_t dimension, size_t order, const std::string& extrapolation);


	/// <summary>
	/// Let \f$ f_{1}(x, y, z) \f$ and \f$ f_{2}(x, y, z) \f$ be 3-dimensional functions defined on 3 sets of gridpoints by a cube object
	/// with dimensions of \f$ n_{1} \times n_{2} \times n_{3} \f$.
	/// Each function can be varied in one dimension to produce a pair of piecewise polynomials \f$ p_{ij} ^ 1, p_{ij} ^ 2 \f$ for each unique combination
	/// of points in the other two fixed dimensions \f$ i \f$ and \f$ j \f$.
	///
	/// A mean vector is given for each dimension \f$ \mu_{i} \f$, which have lengths \f$ m_{1},  m_{2} \f$ and \f$ m_{3} \f$,
	/// where \f$ m_{i} = n_{i} \f$ with the exception of the varied dimension, along with a length-3 weights vector \f$ w \f$. 
	///
	/// Denote the varied dimension as \f$ v_{1} \f$ and fixed dimensions by \f$  v_{2},  v_{3}  \f$ where  \f$ v_{i} \in \{x=1, y=2, z=3\}\f$
	/// The weighted mean for each fixed i-j combination is 
	/// \f$ \mu_{ijk} = w_{v_{1}}\mu_{v_{1} k} + w_{v_{2}}\mu_{v_{2} i} + w_{v_{3}}\mu_{v_{3} j} \f$.
	/// The output is a cube with dimensions of \f$ m_{v_{1}} \times n_{v_{2}} \times n_{v_{3}} \f$, populated with vectors of expectations.
	/// That is, for each pair of piecewise polynomials on the fixed gridpoints \f$ p_{ij} ^ 1, p_{ij} ^ 2 \f$,
	/// calculate \f$ m_{v_{1}} \f$ expectations, such that each output entry is \f$ \mathbb{E}_{k}[ \f$max\f$ (p_{ij} ^ 1(K), p_{ij} ^ 2(K))] \f$ 
	/// where \f$ K \sim N(\mu_{ijk}, \sigma) \f$
	/// </summary>
	/// <param name="output">The output cube.</param>
	/// <param name="x">The grid points to calculate expectations on. Must be equivalent to one dimension of function cubes.</param>
	/// <param name="F_X1">The points of  \f$ f_{1} \f$  defined at points on a 3-dimensional grid.</param>
	/// <param name="F_X2">The points of  \f$ f_{2} \f$  defined at points on a 3-dimensional grid.</param>
	/// <param name="no_layers">The number of layers (matrices) in the function sube.</param>
	/// <param name="mu_1">Means for first argument of function (rows).</param>
	/// <param name="mu_2">Means for second argument of function (columns).</param>
	/// <param name="mu_3">Means for third argument of function (layers).</param>
	/// <param name="standard_deviation">The standard deviation.</param>
	/// <param name="l">The lower limit of integration.</param>
	/// <param name="h">The upper limit of integration.</param>
	/// <param name="weights">Length 3 weights vector</param>
	/// <param name="dimension">The dimension to calculate expectations on.</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 weighted_gaussian_expectation_max_3d(
		d_cube& output, d_vector const& x, d_cube const& F_X1, d_cube const& F_X2, size_t no_layers,
		d_vector const& mu_1, d_vector const& mu_2, d_vector const& mu_3,
		double standard_deviation, double l, double h, d_vector const& weights, size_t dimension, size_t order, const std::string& extrapolation, bool check_input = true);

	const std::string weighted_gaussian_expectation_max_3d_check(
		d_vector const& x, d_cube const& F_X1, d_cube const& F_X2,
		size_t no_layers,
		d_vector const& mu_1, d_vector const& mu_2, d_vector const& mu_3,
		double standard_deviation,
		double l, double h,
		d_vector const& weights,
		size_t dimension, size_t order, const std::string& extrapolation);

}
