#include "HullWhiteCalibrator.h"

void HullWhiteCalibrator::initialise(const std::vector<size_t>& grid_sizes, const MFMBase* model)
{
	TimeSliceCalibrator::initialise(grid_sizes, model);
	m_arrow_debreu.resize(grid_sizes[0]);
	std::fill(m_arrow_debreu.begin(), m_arrow_debreu.end(), 1.0);
	m_calibration_factors.resize(grid_sizes[0]);
	m_analytic_xlibor = false; 
}

std::unique_ptr<HullWhiteCalibrator> HullWhiteCalibrator::construct()
{
	std::unique_ptr<HullWhiteCalibrator> ad = std::unique_ptr<HullWhiteCalibrator>(
		new HullWhiteCalibrator());
	return std::move(ad);
}

void HullWhiteCalibrator::calculate_next_spot_values(std::vector<std::unique_ptr<TimeSliceProcess>>& next_processes)
{
	if (!m_analytic_xlibor)
		next_processes[0]->calculate_spot();

	if (!m_analytic_zlibor)
		next_processes[2]->calculate_spot();

	// this will calculate analytic libors if the flags are set to do so, otherwise we'll have the numerical ones
	UsesHullWhiteJ0::calculate_next_spot_values(next_processes);

	// fx process
	next_processes[1]->calculate_spot();

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

	// fit foreign libors to a piecewise polynomialwith linear extrapolation
	next_processes[2]->fit_spot("linear");

}

void HullWhiteCalibrator::calibrate_x(const InterestRateTimeSliceProcess* currency)
{
	std::transform(currency->get_spot()->begin(), currency->get_spot()->end(), m_calibration_factors.begin(),
		[&](double L) -> double {return (1.0 / (1.0 + currency->get_accrual_factor() * L)); });

	// multiply each calibration factor by corresponding arrow debreu
	std::transform(m_calibration_factors.begin(), m_calibration_factors.end(), m_arrow_debreu.begin(), m_calibration_factors.begin(),
		std::multiplies<double>());
}

void HullWhiteCalibrator::calibrate_x_T0(const InterestRateTimeSliceProcess* currency)
{
	assert(m_calibration_factors.size() == 1);
	m_calibration_factors[0] = currency->get_t1_discount_factor();
}

void HullWhiteCalibrator::calibrate_at_T0(std::vector<std::unique_ptr<TimeSliceProcess>>& processes_at_T0)
{
	calibrate_x_T0(static_cast<const InterestRateTimeSliceProcess*>(processes_at_T0[0].get()));
}

void HullWhiteCalibrator::calibrate(std::vector<std::unique_ptr<TimeSliceProcess>>& processes)
{
	calibrate_x(static_cast<const InterestRateTimeSliceProcess*>(processes[0].get()));
}

void HullWhiteCalibrator::calibrate_next(
	TimeSliceCalibrator& next_calibrator,
	std::vector<std::unique_ptr<TimeSliceProcess>>& current_processes, std::vector<std::unique_ptr<TimeSliceProcess>>& next_processes)
{

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// J0x
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	gelib::d_vector* J0x = next_processes[0]->get_J0_grid();
	const gelib::d_vector* next_xgrid = next_processes[0]->get_grid();
	const gelib::d_vector* means_grid = current_processes[0]->get_conditional_means();

	// Arrow debreu for the next grid
	gelib::d_matrix b_in(next_processes[0]->total_grid_points(), current_processes[0]->total_grid_points());
	gelib::d_vector cf;

	// every row is identical
	for (size_t j = 0; j < b_in.rows(); ++j)
	{
		cf = m_calibration_factors;
		b_in.insert_vector(cf, (int)j, -1);
		static_cast<HullWhiteCalibrator&>(next_calibrator).m_arrow_debreu[j] = 0.0;
	}

	gelib::basis_expectation_sum_1d(static_cast<HullWhiteCalibrator&>(next_calibrator).m_arrow_debreu, b_in, *next_xgrid, *means_grid,
		current_processes[0]->get_conditional_standard_deviation(), MFMBase::minus_inf, MFMBase::plus_inf, m_model->get_order(), "constant", MFMBase::DEBUG_ON);

	// J0 for the next grid
	double recip_sigma = 1.0 / current_processes[0]->get_conditional_standard_deviation();


	for (size_t j = 0; j < next_processes[0]->total_grid_points(); ++j)
	{
		(*J0x)[j] = 0.0;
		for (size_t i = 0; i < m_calibration_factors.size(); ++i)
			(*J0x)[j] += m_calibration_factors[i] * Algorithm::normal_CDF(((*means_grid)[i] - (*next_xgrid)[j]) * recip_sigma);
	}

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// J0y
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	gelib::d_vector* J0y = next_processes[1]->get_J0_grid();
	const gelib::d_vector* next_ygrid = next_processes[1]->get_grid();

	for (size_t j = 0; j < J0y->size(); ++j)
	{
		(*J0y)[j] = m_gauss_J0[0]->e_beta_expectation((*next_ygrid)[j]) / static_cast<FXTimeSliceProcess*>(next_processes[1].get())->get_spot_fx();
	}

	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// J0z
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	gelib::d_vector* J0z = next_processes[2]->get_J0_grid();
	const gelib::d_vector* next_zgrid = next_processes[2]->get_grid();

	for (size_t j = 0; j < J0z->size(); ++j)
	{
		(*J0z)[j] = m_gauss_J0[1]->e_beta_expectation((*next_zgrid)[j]);
	}

}

std::vector<double>& HullWhiteCalibrator::get()
{
	return m_arrow_debreu;
}
