#include "CalibratorTwoCurrency.h"

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//===================================================================================================================
//
// Constructor 
//
//===================================================================================================================
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

std::unique_ptr<TimeSliceCalibrator> TwoCurrencyArrowDebreuCalibrator::construct()
{
	return std::move(std::unique_ptr<TwoCurrencyArrowDebreuCalibrator>(
		new TwoCurrencyArrowDebreuCalibrator()));
}

void TwoCurrencyArrowDebreuCalibrator::initialise(const std::vector<size_t>& grid_sizes, const MFMBase* model)
{
	// initialise the single currency driving model
	SingleCurrencyCalibrator::initialise(grid_sizes, model);

	// Note that grid_sizes array here has x, y, z ordering 
	m_domestic_factors.resize(grid_sizes[0]);

	// arrow debreu cube
	m_arrow_debreu.resize(grid_sizes[0], grid_sizes[1], grid_sizes[2]);
	std::fill(m_arrow_debreu.begin(), m_arrow_debreu.end(), 1.0);

	// accumulated arrow debreu cube
	m_B0_C.resize(grid_sizes[0], grid_sizes[1], grid_sizes[2]);
	std::fill(m_B0_C.begin(), m_B0_C.end(), 1.0);

	m_larrow_grid_k.resize(grid_sizes.size());

	for (size_t i = 0; i < m_larrow_grid_k.size(); ++i)
		m_larrow_grid_k[i].resize(grid_sizes[i]);

	// Note - assumes 2 -> 3 -> 1 integration ordering and grids and means to forward equal dimensions
	m_fx_drift_grid.resize(m_larrow_grid_k[0].size(), m_larrow_grid_k[1].size(), m_larrow_grid_k[2].size());
	m_mu3.resize(m_larrow_grid_k[0].size(), m_larrow_grid_k[1].size(), m_larrow_grid_k[2].size());
	m_mu2.resize(m_larrow_grid_k[0].size(), 1, m_larrow_grid_k[2].size());
	m_mu1.resize(m_larrow_grid_k[0].size(), 1, 1);

	// Working arrays for J0y calculation
	m_fixedy_mu3_vector.resize(grid_sizes[1]);
	m_fixedy_arrow_debreu_factors.resize(grid_sizes[1]);

	// Working arrays for J0z calculation
	m_calN_grid.resize(grid_sizes[1]);
	m_calN_plusmu.resize(grid_sizes[1]);
	m_calN_function.resize(grid_sizes[1]);
	m_Hhat_function.resize(grid_sizes[1]);
	m_Hhat_interpolation.resize(grid_sizes[1]);

}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//===================================================================================================================
//
// All necessary work at the current time step
//
//===================================================================================================================
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void TwoCurrencyArrowDebreuCalibrator::calibrate_at_T0(std::vector<std::unique_ptr<TimeSliceProcess>>& processes_at_T0, SliceProcessCorrelations const& T0_to_T1_increment)
{

	SingleCurrencyCalibrator::calibrate_at_T0(processes_at_T0);

	assert(m_B0_C.size() == 1);

	m_domestic_factors[0] = static_cast<const InterestRateTimeSliceProcess*>(processes_at_T0[0].get())->get_t1_discount_factor();

	// B_0_C
	m_B0_C[0] = m_domestic_factors[0];

	// means
	m_mu1[0] = processes_at_T0[0]->get_grid_center(); // = 0
	m_mu2[0] = processes_at_T0[2]->get_grid_center(); // = 0

	// note we need access to the next fx process to calculate m_mu3[0] for T0 calibrator, so this is assigned in calibrate_next_at_T0 and not here
}

void TwoCurrencyArrowDebreuCalibrator::calibrate(InterestRateTimeSliceProcess* dom_process, FXTimeSliceProcess* fx_process, InterestRateTimeSliceProcess* fgn_process, SliceProcessCorrelations3d const& next_increment)
{
	std::transform(dom_process->get_spot()->begin(), dom_process->get_spot()->end(), m_domestic_factors.begin(),
		[&](double L) -> double {return (1.0 / (1.0 + dom_process->get_accrual_factor() * L)); });

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// set up the B_0_C cube
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	const auto mult_func = [](double& a, double& b)
	{
		double prod = a * b;
		std::swap(prod, b);
	};

	//multiply every row by 1 / (1 + alpha L(x))
	for (size_t i = 0; i < m_B0_C.cols(); ++i)
		for (size_t j = 0; j < m_B0_C.layers(); ++j)
			m_B0_C.vector_op(mult_func, m_domestic_factors, -1, (int)i, (int)j);

	// below line assumes arrow debreu and B0 have the same dimensions
	assert(gelib::d_cube::matching_dimensions(m_B0_C, m_arrow_debreu));

	// multiply B0 by arrow debreu prices (we will have them from calibration at previous slice)
	std::transform(m_B0_C.begin(), m_B0_C.end(), m_arrow_debreu.begin(), m_B0_C.begin(), std::multiplies<double>());

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// x means and z means
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	gelib::d_vector tmp_mu(fgn_process->total_grid_points());

	for (size_t i = 0; i < m_mu2.rows(); ++i)
	{
		std::copy(m_larrow_grid_k[2].begin(), m_larrow_grid_k[2].end(), tmp_mu.begin());

		// load tmp to mu2 cube layer at row i, note that mu2 just has one column
		m_mu2.insert_vector(tmp_mu, (int)i, 0, -1);
	}

	// mu1
	std::copy(m_larrow_grid_k[0].begin(), m_larrow_grid_k[0].end(), m_mu1.begin());

}

void TwoCurrencyArrowDebreuCalibrator::calibrate(std::vector<std::unique_ptr<TimeSliceProcess>>& processes, SliceProcessCorrelations const& next_increment)
{
	// calulate domestic currency factors
	SingleCurrencyCalibrator::calibrate(processes);

	calibrate(
		static_cast<InterestRateTimeSliceProcess*>(processes[0].get()),
		static_cast<FXTimeSliceProcess*>(processes[1].get()),
		static_cast<InterestRateTimeSliceProcess*>(processes[2].get()),
		static_cast<SliceProcessCorrelations3d const&>(next_increment));
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//===================================================================================================================
//
// All necessary work to calibrate the next time step
//
//===================================================================================================================
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void TwoCurrencyArrowDebreuCalibrator::calibrate_next(TimeSliceCalibrator& next_calibrator, std::vector<std::unique_ptr<TimeSliceProcess>>& current_processes, std::vector<std::unique_ptr<TimeSliceProcess>>& next_processes, SliceProcessCorrelations const& next_increment)
{
	SingleCurrencyCalibrator::calibrate_next(next_calibrator, current_processes, next_processes);

	// cast to the model-specific classes we need and run the actual next calibrate/J0 driver functions
	calibrate_next_T0plus(
		static_cast<TwoCurrencyArrowDebreuCalibrator&>(next_calibrator),
		static_cast<const InterestRateTimeSliceProcess*>(current_processes[0].get()),
		static_cast<FXTimeSliceProcess*>(current_processes[1].get()),
		static_cast<const InterestRateTimeSliceProcess*>(current_processes[2].get()),
		static_cast<InterestRateTimeSliceProcess*>(next_processes[0].get()),
		static_cast<FXTimeSliceProcess*>(next_processes[1].get()),
		static_cast<InterestRateTimeSliceProcess*>(next_processes[2].get()),
		static_cast<SliceProcessCorrelations3d const&>(next_increment));
}


void TwoCurrencyArrowDebreuCalibrator::calibrate_next_at_T0(TimeSliceCalibrator& next_calibrator, std::vector<std::unique_ptr<TimeSliceProcess>>& current_processes, std::vector<std::unique_ptr<TimeSliceProcess>>& next_processes, SliceProcessCorrelations const& next_increment)
{

	SingleCurrencyCalibrator::calibrate_next_at_T0(next_calibrator, current_processes, next_processes);

	FXTimeSliceProcess* fx_process = static_cast<FXTimeSliceProcess*>(current_processes[1].get());
	FXTimeSliceProcess* next_fx_process = static_cast<FXTimeSliceProcess*>(next_processes[1].get());

	fx_process->set_up_drift(static_cast<SliceProcessCorrelations3d const&>(next_increment), "linear", *next_fx_process->get_grid(), next_fx_process->get_drift_fitting_grid());

	m_mu3[0] = fx_process->mean_drift_at_T0();
	m_fx_drift_grid[0] = m_mu3[0];

	// do everything else as part of general calibrate next
	calibrate_next(
		static_cast<TwoCurrencyArrowDebreuCalibrator&>(next_calibrator),
		static_cast<const InterestRateTimeSliceProcess*>(current_processes[0].get()),
		static_cast<FXTimeSliceProcess*>(current_processes[1].get()),
		static_cast<const InterestRateTimeSliceProcess*>(current_processes[2].get()),
		static_cast<InterestRateTimeSliceProcess*>(next_processes[0].get()),
		static_cast<FXTimeSliceProcess*>(next_processes[1].get()),
		static_cast<InterestRateTimeSliceProcess*>(next_processes[2].get()),
		static_cast<SliceProcessCorrelations3d const&>(next_increment));

}

void TwoCurrencyArrowDebreuCalibrator::calibrate_next_T0plus(TwoCurrencyArrowDebreuCalibrator& next_calibrator, const InterestRateTimeSliceProcess* dom_process, FXTimeSliceProcess* fx_process, const InterestRateTimeSliceProcess* fgn_process, InterestRateTimeSliceProcess* next_dom_process, FXTimeSliceProcess* next_fx_process, InterestRateTimeSliceProcess* next_fgn_process, SliceProcessCorrelations3d const& next_increment)
{

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// y-means grid and fx drift
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	// set up the function f_{k-1} in the fx process at slice k-1 using grids stored fx process at slice k
	fx_process->set_up_drift(next_increment, "linear", *next_fx_process->get_grid(), next_fx_process->get_drift_fitting_grid());

	// calculate mean drift values into this vector
	gelib::d_vector tmp_drift(fx_process->total_grid_points());

	const auto y_to_mu3_vec = [&](size_t x_index, size_t z_index, gelib::d_vector const& y_pts, gelib::d_vector& mu3_pts)
	{
		// calculate f_{k-1}, with the result doing in to tmp_drift
		fx_process->mean_drift(next_increment, m_larrow_grid_k[0][x_index], m_larrow_grid_k[2][z_index], m_larrow_grid_k[1], tmp_drift);

		double z = m_larrow_grid_k[2][z_index] - next_increment.alpha() * m_larrow_grid_k[0][x_index];
		double lhs = next_increment.beta() * m_larrow_grid_k[0][x_index] + next_increment.gamma() * z;

		// calculate mu3(gi) at each gridpoint
		std::transform(tmp_drift.begin(), tmp_drift.end(), mu3_pts.begin(), [&lhs](double drift) { return lhs + drift; });

		// store the drift values at each point, will be useful later on in calibration
		m_fx_drift_grid.insert_vector(tmp_drift, (int)x_index, -1, (int)z_index);

	};

	// set tmp to size of the fx grid
	gelib::d_vector tmp_mu(fx_process->total_grid_points());

	for (size_t i = 0; i < m_mu3.rows(); ++i)
	{
		for (size_t j = 0; j < m_mu3.layers(); ++j)
		{
			// load mu3 values to tmp
			y_to_mu3_vec(i, j, m_larrow_grid_k[1], tmp_mu);

			// load tmp to mu3 cube at i'th row and jth layer
			m_mu3.insert_vector(tmp_mu, (int)i, -1, (int)j);
		}
	}


	// do everything else as part of general calibrate next
	calibrate_next(next_calibrator, dom_process, fx_process, fgn_process, next_dom_process, next_fx_process, next_fgn_process, next_increment);

}

void TwoCurrencyArrowDebreuCalibrator::calibrate_next(
	TwoCurrencyArrowDebreuCalibrator& next_calibrator,
	const InterestRateTimeSliceProcess* dom_process, FXTimeSliceProcess* fx_process, const InterestRateTimeSliceProcess* fgn_process,
	InterestRateTimeSliceProcess* next_dom_process, FXTimeSliceProcess* next_fx_process, InterestRateTimeSliceProcess* next_fgn_process,
	SliceProcessCorrelations3d const& next_increment)
{

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// Make s arrow grids for the next slice
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	// forward x grid
	GridMaker1d::make_grid_zero_center(
		next_dom_process->total_grid_points(), next_dom_process->get_grid_spread_factor(), next_increment.get_Mk_Ck_Mkt()[0],
		next_calibrator.m_larrow_grid_k[0]);

	// forward y grid - note spread not increased by 10% TODO
	GridMaker1d::make_grid(next_fx_process->total_grid_points(), next_fx_process->get_grid_spread_factor(),
		next_fx_process->get_grid_center(), next_fx_process->get_uncorrelated_grid_variance(), next_calibrator.m_larrow_grid_k[1]);

	// forward z grid
	GridMaker1d::make_grid_zero_center(next_fgn_process->total_grid_points(), next_fgn_process->get_grid_spread_factor(), next_increment.get_Mk_Ck_Mkt()[8],
		next_calibrator.m_larrow_grid_k[2]);

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// Do the integration
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	try
	{
		// check if we have a global calculator, if not throw to the default method
		if (!m_arrow_debreu_calculator_ref)
			throw 0;

		// this will throw if something went wrong with the transform
		gelib::BasisPiecewisePoly3d::gaussian_expectation_sum_3d(
			*m_arrow_debreu_calculator_ref, next_calibrator.m_arrow_debreu, next_calibrator.m_larrow_grid_k,
			m_B0_C, m_mu3, m_mu2, m_mu1,
			std::sqrt(next_increment.get_Mk_Ck_Ckminus1_Mkt()[0]),
			std::sqrt(next_increment.get_Mk_Ck_Ckminus1_Mkt()[4]),
			std::sqrt(next_increment.get_Mk_Ck_Ckminus1_Mkt()[8]),
			MFMBase::minus_inf, MFMBase::plus_inf,
			MFMBase::minus_inf, MFMBase::plus_inf,
			MFMBase::minus_inf, MFMBase::plus_inf, m_model->get_order(), "constant");
	}
	// Each calculation is independent at each time slice
	catch (...)
	{
		gelib::basis_expectation_sum_3d_triple(
			next_calibrator.m_arrow_debreu, next_calibrator.m_larrow_grid_k, m_B0_C,
			m_mu3, m_mu2, m_mu1,
			std::sqrt(next_increment.get_Mk_Ck_Ckminus1_Mkt()[0]),
			std::sqrt(next_increment.get_Mk_Ck_Ckminus1_Mkt()[4]),
			std::sqrt(next_increment.get_Mk_Ck_Ckminus1_Mkt()[8]),
			MFMBase::minus_inf, MFMBase::plus_inf,
			MFMBase::minus_inf, MFMBase::plus_inf,
			MFMBase::minus_inf, MFMBase::plus_inf, m_model->get_order(),
			"constant", MFMBase::DEBUG_ON, 2, 3, 1);
	}
	// TODO extrapolation only constant for J0

	// Calculate next J0y and J0z, already have J0x from the SingleCurrencyCalibrator 
	calculate_next_J0y(next_calibrator, dom_process, fx_process, fgn_process, next_dom_process, next_fx_process, next_fgn_process, next_increment);

	// note that the J0z are dependent on the FX values on the next grid
	next_fx_process->calculate_spot();

	calculate_next_J0z(next_calibrator, dom_process, fx_process, fgn_process, next_dom_process, next_fx_process, next_fgn_process, next_increment);

}

void TwoCurrencyArrowDebreuCalibrator::calculate_next_J0y(TwoCurrencyArrowDebreuCalibrator& next_calibrator,
	const InterestRateTimeSliceProcess* dom_process, const FXTimeSliceProcess* fx_process, const InterestRateTimeSliceProcess* fgn_process,
	InterestRateTimeSliceProcess* next_dom_process, FXTimeSliceProcess* next_fx_process, InterestRateTimeSliceProcess* next_fgn_process,
	SliceProcessCorrelations3d const& next_increment)
{

	// 1 / sqrt(Cyy^k - Cyy^k-1)
	double recip_sigma = 1.0 / std::sqrt(next_increment.get_increment_covariance()[4]);

	const gelib::d_vector* next_grid = next_fx_process->get_grid();
	gelib::d_vector* next_J0_grid = next_fx_process->get_J0_grid();

	// iterate over J0y(y*)
	for (size_t k = 0; k < next_grid->size(); ++k)
	{
		(*next_J0_grid)[k] = 0.0;

		// iterate over z
		for (size_t j = 0; j < m_B0_C.layers(); ++j)
		{
			// iterate over x
			for (size_t i = 0; i < m_B0_C.rows(); ++i)
			{

				m_fx_drift_grid.extract_vector(m_fixedy_mu3_vector, (int)i, -1, (int)j);

				// y vector of arrow debreu with x, z fixed
				// note that B_0_C is already arrow debreu multiplied by 1 / (1 + alpha L(x))
				m_B0_C.extract_vector(m_fixedy_arrow_debreu_factors, (int)i, -1, (int)j);

				for (size_t l = 0; l < m_fixedy_mu3_vector.size(); ++l)
					(*next_J0_grid)[k] += m_fixedy_arrow_debreu_factors[l] * Algorithm::normal_CDF((m_fixedy_mu3_vector[l] - (*next_grid)[k]) * recip_sigma);


			}
		}
	}

}

void TwoCurrencyArrowDebreuCalibrator::calculate_next_J0z(TwoCurrencyArrowDebreuCalibrator& next_calibrator,
	const InterestRateTimeSliceProcess* dom_process, const FXTimeSliceProcess* fx_process, const InterestRateTimeSliceProcess* fgn_process,
	InterestRateTimeSliceProcess* next_dom_process, FXTimeSliceProcess* next_fx_process, InterestRateTimeSliceProcess* next_fgn_process,
	SliceProcessCorrelations3d const& next_increment)
{

	const gelib::d_vector& next_z_grid = *next_fgn_process->get_grid();
	const gelib::d_vector& next_y_grid = *next_fx_process->get_grid();
	gelib::d_vector& next_J0z_grid = *next_fgn_process->get_J0_grid();

	// covariance matrix in increment
	const gelib::d_matrix Ck_Ckminus1 = next_increment.get_increment_covariance();

	// zero-centered grid with variance C^yy_{k} - C^yy_{k - 1}
	GridMaker1d::make_grid_zero_center(
		next_fx_process->total_grid_points(), next_fx_process->get_grid_spread_factor(), Ck_Ckminus1[4], next_calibrator.m_calN_grid);

	double sd = std::sqrt(Ck_Ckminus1[4]);

	//  (C^yz_{k} - C^yz_{k - 1}) / (C^yy_{k} - C^yy_{k - 1})
	double omega_k = Ck_Ckminus1[7] / Ck_Ckminus1[4];

	//  (C^zz_{k} - C^zz_{k - 1}) - ((C^yz_{k} - C^yz_{k - 1})^2 / (C^yy_{k} - C^yy_{k - 1}))
	double sigma_k = std::sqrt(Ck_Ckminus1[8] - Ck_Ckminus1[7] * omega_k);

	double recip_sigma_k = 1.0 / sigma_k;

	gelib::PiecewisePoly1d calN_function_poly;
	gelib::PiecewisePoly1d Hhat_function_poly;

	char extrapolation = gelib::PiecewisePoly1d::parse_extrapolation_string("linear");

	// make sure spot fx grid is fitted
	next_fx_process->fit_spot("linear");

	// uncorrelated z value
	double z;

	// iterate over z*
	for (size_t i = 0; i < next_z_grid.size(); ++i)
	{
		next_J0z_grid[i] = 0.0;

		// iterate over z in arrow debreu
		for (size_t l = 0; l < m_B0_C.layers(); ++l)
		{

			// iterate over x in arrow debreu
			for (size_t j = 0; j < m_B0_C.rows(); ++j)
			{
				// z* - z(gi)
				z = next_z_grid[i] - (m_larrow_grid_k[2][l] - next_increment.alpha() * m_larrow_grid_k[0][j]);

				//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
				// fit Hhat on a y-grid for this z* - z(gi)
				//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
				for (size_t k = 0; k < next_y_grid.size(); ++k)
				{
					
					// get a grid of mu + calN
					std::transform(next_calibrator.m_calN_grid.begin(), next_calibrator.m_calN_grid.end(), next_calibrator.m_calN_plusmu.begin(), [&](double N) {return N + next_y_grid[k]; });

					// interpolate FX values on (mu + calN) grid
					next_fx_process->interpolate_spot(next_calibrator.m_calN_plusmu, next_calibrator.m_calN_function);

					for (size_t t = 0; t < next_calibrator.m_calN_function.size(); ++t)
						next_calibrator.m_calN_function[t] *= Algorithm::normal_CDF((next_calibrator.m_calN_grid[t] * omega_k - z) * recip_sigma_k);

					// fit the function for fixed z*, mu
					calN_function_poly.fit(static_cast<const double*>(&next_calibrator.m_calN_grid[0]), &next_calibrator.m_calN_function[0], next_calibrator.m_calN_grid.size(), m_model->get_order(), extrapolation);

					// This is Hhat(z*, mu) for a single mu with fixed z*
					next_calibrator.m_Hhat_function[k] = calN_function_poly.gaussian_expectation(0.0, sd, MFMBase::minus_inf, MFMBase::plus_inf);

					//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
					// Uncomment below to fit Hhat to it's closed form solution in the Hull White case 
					//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

					/*
					next_calibrator.m_Hhat_function[k] = Gaussian2d::e_beta_expectation(
						1,
						z / sigma_k,
						0.0,
						next_y_grid[k],
						1 + omega_k * omega_k * Ck_Ckminus1[4] / (sigma_k * sigma_k),
						Ck_Ckminus1[4],
						omega_k * Ck_Ckminus1[4] / sigma_k);
						*/

				}

				// fit Hhat on the y-grid
				Hhat_function_poly.fit(static_cast<const double*>(&next_y_grid[0]), &next_calibrator.m_Hhat_function[0], next_y_grid.size(), m_model->get_order(), extrapolation);


				m_fx_drift_grid.extract_vector(m_fixedy_mu3_vector, (int)j, -1, (int)l);

					// y vector of arrow debreu with z, x fixed
					// note that B_0_C is already arrow debreu multiplied by 1 / (1 + alpha L(x))
				m_B0_C.extract_vector(m_fixedy_arrow_debreu_factors, (int)j, -1, (int)l);

				// interpolate y means on Hhat function
				Hhat_function_poly.evaluate(m_fixedy_mu3_vector, m_Hhat_interpolation);

				//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
				// Uncomment below to use closed form Hull White H(z*, mu(gi)) directly instead of interpolating
				//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

				/*
				std::transform(m_fixedy_mu3_vector.begin(), m_fixedy_mu3_vector.end(), m_Hhat_interpolation.begin(),
					[&](double x) -> double {return Gaussian2d::e_beta_expectation(
						1.0,
						z / sigma_k,
						0.0,
						x,
						1.0 + omega_k * omega_k * Ck_Ckminus1[4] / (sigma_k * sigma_k),
						Ck_Ckminus1[4],
						omega_k * Ck_Ckminus1[4] / sigma_k); });
				*/

				next_J0z_grid[i] += std::inner_product(m_fixedy_arrow_debreu_factors.begin(), m_fixedy_arrow_debreu_factors.end(), m_Hhat_interpolation.begin(), 0.0);

			}
			// end of iterate over x
		}
		// end of iterate over z

		// J0(z*) is set by this point. For safety, ensure we didn't exceed upper bound (foreign discount factor) on J0z with numerical error
		next_J0z_grid[i] = std::min(next_J0z_grid[i], next_fgn_process->get_discount_factor() - 1e-10);
	}
	// end of iterate over z*
}

double TwoCurrencyArrowDebreuCalibrator::roll_back(std::vector<double> const& spread) const
{
	return std::inner_product(m_B0_C.begin(), m_B0_C.end(), spread.begin(), 0.0);
}

void TwoCurrencyArrowDebreuCalibrator::calculate_next_spot_values(std::vector<std::unique_ptr<TimeSliceProcess>>& next_processes)
{
	// note that FX spot rates already calculated at this point

	// x libors
	SingleCurrencyCalibrator::calculate_next_spot_values(next_processes);

	// z libors
	next_processes[2]->calculate_spot();

	// fit x libors to a piecewise polynomial with linear extrapolation
	next_processes[0]->fit_spot("linear");

	// fit z libors to a piecewise polynomial with linear extrapolation
	next_processes[2]->fit_spot("linear");
}

void TwoCurrencyArrowDebreuCalibrator::set_arrow_debreu_calculator(gelib::BasisPiecewisePoly3d* arrow_debreu)
{
	m_arrow_debreu_calculator_ref = arrow_debreu;
}

gelib::d_cube const& TwoCurrencyArrowDebreuCalibrator::get_arrow_debreu() const
{
	return m_arrow_debreu;
}

std::vector<gelib::d_vector>& TwoCurrencyArrowDebreuCalibrator::get_s_l_arrow()
{
	return m_larrow_grid_k;
}

gelib::d_cube const& TwoCurrencyArrowDebreuCalibrator::get_B_0_i() const
{
	return m_B0_C;
}

gelib::d_cube const& TwoCurrencyArrowDebreuCalibrator::get_mu1() const
{
	return m_mu1;
}

gelib::d_cube const& TwoCurrencyArrowDebreuCalibrator::get_mu2() const
{
	return m_mu2;
}

gelib::d_cube const& TwoCurrencyArrowDebreuCalibrator::get_mu3() const
{
	return m_mu3;
}

std::vector<double>& TwoCurrencyArrowDebreuCalibrator::get()
{
	return m_B0_C;
}