#include "UsesHullWhiteJ0.h"

void UsesHullWhiteJ0::calculate_next_spot_values(std::vector<std::unique_ptr<TimeSliceProcess>>& next_processes)
{
	// if this flag is true, replace the x libors with their analytic solutions
	if (m_analytic_xlibor)
	{
		gelib::d_vector R(next_processes[0]->total_grid_points());
		gelib::d_vector* next_libors = next_processes[0]->get_spot_grid();
		const gelib::d_vector* next_xgrid = next_processes[0]->get_grid();

		for (size_t i = 0; i < R.size(); ++i)
			R[i] = m_mu_R + std::sqrt(m_C_RR / next_processes[0]->get_process_variance()) * (*next_xgrid)[i];

		for (size_t i = 0; i < next_libors->size(); ++i)
			(*next_libors)[i] = (std::exp(R[i]) - 1.0) / static_cast<InterestRateTimeSliceProcess*>(next_processes[0].get())->get_accrual_factor();

		// refit
		next_processes[0]->allow_refit();
		next_processes[0]->fit_spot("linear");
	}


	// if this flag is true, replace the z libors with their analytic solutions
	if (m_analytic_zlibor)
	{
		gelib::d_vector Rhat(next_processes[2]->total_grid_points());
		gelib::d_vector* next_zlibors = next_processes[2]->get_spot_grid();
		const gelib::d_vector* next_zgrid = next_processes[2]->get_grid();

		for (size_t i = 0; i < Rhat.size(); ++i)
			Rhat[i] = m_mu_Rhat + std::sqrt(m_C_RhatRhat / next_processes[2]->get_process_variance()) * (*next_zgrid)[i];

		for (size_t i = 0; i < next_zlibors->size(); ++i)
			(*next_zlibors)[i] = (std::exp(Rhat[i]) - 1.0) / static_cast<InterestRateTimeSliceProcess*>(next_processes[2].get())->get_accrual_factor();

		// refit
		next_processes[2]->allow_refit();
		next_processes[2]->fit_spot("linear");
	}
}

void UsesHullWhiteJ0::prepare_analytic_xlibor(double mu_R, double C_RR)
{
	m_mu_R = mu_R;
	m_C_RR = C_RR;
	m_analytic_xlibor = true;
}

void UsesHullWhiteJ0::prepare_analytic_zlibor(double mu_Rhat, double C_RhatRhat)
{
	m_mu_Rhat = mu_Rhat;
	m_C_RhatRhat = C_RhatRhat;
	m_analytic_zlibor = true;
}

void UsesHullWhiteJ0::pre_calibrate(double mu_x_1, double mu_y_1, double Cxx_1, double Cyy_1, double Cxy_1, double mu_x_2, double mu_y_2, double Cxx_2, double Cyy_2, double Cxy_2)
{
	m_gauss_J0[0] = Gaussian2d::construct(mu_x_1, mu_y_1, Cxx_1, Cyy_1, Cxy_1);
	m_gauss_J0[1] = Gaussian2d::construct(mu_x_2, mu_y_2, Cxx_2, Cyy_2, Cxy_2);

	// J0y
	m_gauss_J0[0]->fix_beta(-1.0);
	// J0z
	m_gauss_J0[1]->fix_beta(1.0);
}
