#include "MFMSingleCurrency.h"

MFMSingleCurrency::MFMSingleCurrency(
	double mean_reversion_param, size_t gridpoints_per_slice,
	double sd_spread, double discount_factor_t1,
	size_t poly_order, gelib::d_vector const& times, gelib::d_vector const& forward_libors,
	gelib::d_vector const& caplet_volatilities, gelib::d_vector const& x_volatilities, gelib::d_vector const& accrual_factors)
	:
	m_t1_discount_factor(discount_factor_t1),
	m_mean_reversion_param(mean_reversion_param),

	// Call MFMGeneral constructor specialised with a single process type
	MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>(
	poly_order, gridpoints_per_slice, times,

	// provide a constructor for the calibrator used by this model 
	&SingleCurrencyCalibrator::construct,

	// Construct an array of interest rate processes to drive the model 
	std::move(
		InterestRateTimeSliceProcess::construct(
			sd_spread, gridpoints_per_slice, this, mean_reversion_param, times,
			square(x_volatilities),
			TimeSliceProcess::volatilities_to_variances(mean_reversion_param, times, x_volatilities),
			forward_libors, caplet_volatilities, accrual_factors, discount_factor_t1)
	))
{}

gelib::d_vector MFMSingleCurrency::square(gelib::d_vector const& x)
{
	gelib::d_vector out;
	out.resize(x.size());
	std::transform(x.begin(), x.end(), out.begin(), [](double xi) {return xi * xi; });
	return out;
}

double MFMSingleCurrency::value_caplet(double strike_price, size_t which_slice)
{
	std::string extrapolation = "linear";

	if (!m_bls_calibrated)
		throw std::runtime_error("Model must be calibrated before pricing caplets");

	const TimeSlice<SingleCurrencyCalibrator>* prev_slice = m_time_slices[which_slice - 1];
	const TimeSlice<SingleCurrencyCalibrator>* slice = m_time_slices[which_slice];

	gelib::d_vector y1(slice->total_grid_points(0));
	gelib::d_vector y2(slice->total_grid_points(0));

	gelib::d_vector libors;
	get_process<InterestRateTimeSliceProcess>(0, which_slice)->get_spot(libors);

	for (size_t j = 0; j < slice->total_grid_points(0); ++j)
		y1[j] = (libors[j] - strike_price) / (1 + get_process<InterestRateTimeSliceProcess>(0, which_slice)->get_accrual_factor() * libors[j]);

	gelib::d_vector V_prev_slice(prev_slice->total_grid_points(0));

	gelib::gaussian_expectation_max_1d(V_prev_slice, slice->get_grid(0), y1, y2, prev_slice->get_conditional_means(0), prev_slice->get_conditional_sigma(0), minus_inf, plus_inf, m_poly_order, extrapolation, DEBUG_ON);

	return static_cast<SingleCurrencyCalibrator&>(prev_slice->get_calibrator()).weight_vector(V_prev_slice);
}

//////////////////////////////////////////////////////////////////////////////////////
// Extract data from the model using these functions
//////////////////////////////////////////////////////////////////////////////////////

void MFMSingleCurrency::get_discount_factors(gelib::d_vector& out) const
{
	out.resize(m_number_of_slices);
	for (size_t i = 0; i < m_number_of_slices; ++i)
		out[i] = MFMGeneral::get_process<InterestRateTimeSliceProcess>(0, i)->get_discount_factor();
}

void MFMSingleCurrency::get_variances(gelib::d_vector& out) const
{
	MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_variances(out, 0);
}

double MFMSingleCurrency::digital_caplet_price(size_t which_slice, double k) const
{
	return MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_option_in_arrears_price(which_slice, 0, k);
}

double MFMSingleCurrency::digital_caplet_price_inverse(size_t which_slice, double dcp) const
{
	return MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_option_in_arrears_price_inverse(which_slice, 0, dcp);
}

void MFMSingleCurrency::get_grids(gelib::d_matrix& out) const
{
	MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_grids(out, 0);
}

void MFMSingleCurrency::get_conditional_means(gelib::d_matrix& out) const
{
	MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_conditional_means(out, 0);
}

void MFMSingleCurrency::get_ad_grids(gelib::d_matrix& out) const
{
	MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_ad_grids(out);
}

void MFMSingleCurrency::get_J0_grids(gelib::d_matrix& out) const
{
	MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_J0_grids(out, 0);
}

void MFMSingleCurrency::get_libor_grids(gelib::d_matrix& out) const
{
	MFMGeneral<SingleCurrencyCalibrator, InterestRateTimeSliceProcess>::get_spot_grids(out, 0);
}

//////////////////////////////////////////////////////////////////////////////////////
// End of getter functions
//////////////////////////////////////////////////////////////////////////////////////

const std::string MFMSingleCurrency::constructor_check(
	double mean_reversion_param, size_t gridpoints_per_slice, double sd_spread, double discount_factor_t1,
	size_t poly_order, gelib::d_vector const& times, gelib::d_vector const& forward_libors, gelib::d_vector const& caplet_volatilities,
	gelib::d_vector const& x_volatilities, gelib::d_vector const& accrual_factors)
{
	const size_t max_poly_order = gelib_impl::Polynomial1d::MAXIMUM_DEGREE;

	std::array<size_t, 5> sizes = { forward_libors.size(), caplet_volatilities.size(), times.size(), x_volatilities.size(), accrual_factors.size() };

	for (size_t i = 0; i < accrual_factors.size(); i++)
	{
		if (accrual_factors[i] < 0.0)
			return "Accrual factor " + std::to_string(i) + " is negative.";
	}
	for (size_t i = 0; i < times.size() - 1; i++)
	{
		if (times[i + 1] - times[i] < 0.0)
			return "At i = " + std::to_string(i) + " T_{i+1} - T_{i} is negative.";
	}
	if (times.back() > 1000)
		return "Times should not exceed 1000";
	if (sd_spread <= 0)
		return "Spread should be positive.";
	if (!((0.0 <= discount_factor_t1) && (discount_factor_t1 <= 1.0)))
		return "Discount factor should be between 0 and 1";
	if (poly_order > max_poly_order)
		return "Polynomial order is too large or a negative parameter was supplied.";
	if (times.size() < 2)
		return "There should be at least 2 times";
	if (!std::equal(sizes.begin() + 1, sizes.end(), sizes.begin()))
		return "The vectors for accrual factors, times, forward libors, caplet volatilities and x volatilities should all be equal lengths.";
	if (gridpoints_per_slice > MFMBase::MAXIMUM_GRID_POINTS)
		return "Maximum grid points exceeded, was a negative parameter supplied?";

	return "";

}

const std::string MFMSingleCurrency::caplet_valuation_check(double strike_price, size_t which_slice, size_t number_of_slices)
{
	if (strike_price <= 0)
		return "Strike price must be positive.";
	if ((which_slice == 0) || (which_slice > number_of_slices))
		return "Time slice " + std::to_string(which_slice) + " can not be used to value caplet.";
	return "";
}