#include "TimeSliceProcess.h" 

TimeSliceProcess::TimeSliceProcess(double time, double conditional_variance, double variance, double grid_spread, double grid_center, size_t number_of_gridpoints, double conditional_means_factor, const MFMBase* model)
{
	initialise(time, conditional_variance, variance, grid_spread, grid_center, number_of_gridpoints, conditional_means_factor, model);
}

void TimeSliceProcess::initialise(double time, double conditional_variance, double variance, double grid_spread, double grid_center, size_t number_of_gridpoints, double conditional_means_factor, const MFMBase* model)
{// TODO can now put all this logic in a constructor
	m_grid_center = grid_center;
	m_grid_sd_spread_factor = grid_spread;
	m_time = time;
	m_conditional_variance = conditional_variance;
	m_conditional_sigma = std::sqrt(conditional_variance);
	m_variance = variance;
	m_model = model;
	m_spot_fitted = false;

	populate_grid(number_of_gridpoints, grid_spread, grid_center);
	calculate_conditional_means(conditional_means_factor);

	m_J0.resize(number_of_gridpoints);
	m_spot_values.resize(number_of_gridpoints);
	m_spot_fit = std::move(std::unique_ptr<gelib::PiecewisePoly1d> (new gelib::PiecewisePoly1d()));
}

std::unique_ptr<TimeSliceProcess> TimeSliceProcess::construct(double time, double conditional_variance, double variance, double grid_spread, double grid_center, size_t number_of_gridpoints, double conditional_means_factor, const MFMBase* model)
{
	return std::unique_ptr<TimeSliceProcess>(new TimeSliceProcess(time, conditional_variance, variance, grid_spread, grid_center, number_of_gridpoints, conditional_means_factor, model));
}

TimeSliceProcess::TimeSliceProcess(TimeSliceProcess&& tsp, std::unique_ptr<OptionPricer> pricer)
	:
	TimeSliceProcess(std::move(tsp))
{
	m_pricer = std::move(pricer);
}

TimeSliceProcess::TimeSliceProcess(std::unique_ptr<TimeSliceProcess> tsp, std::unique_ptr<OptionPricer> pricer)
	:TimeSliceProcess(std::move(*tsp.get()), std::move(pricer))
{}

void TimeSliceProcess::populate_grid(size_t number_of_gridpoints, double grid_spread, double grid_center)
{
	(grid_center == 0.0) ?
		GridMaker1d::make_grid_zero_center(number_of_gridpoints, grid_spread, m_variance, m_grid)
		:
		GridMaker1d::make_grid(number_of_gridpoints, grid_spread, grid_center, m_variance, m_grid);
}

void TimeSliceProcess::calculate_conditional_means(double conditional_means_factor)
{
	m_conditional_grid_means.resize(total_grid_points());

	for (size_t i = 0; i < total_grid_points(); ++i)
		m_conditional_grid_means[i] = conditional_means_factor * m_grid[i];

}

double TimeSliceProcess::calculate_variance(double mean_reversion_parameter_times_2, double prev_time, double this_time, double prev_volatility_squared, double* prev_variance_intermediate)
{
	*prev_variance_intermediate = (*prev_variance_intermediate) +
		prev_volatility_squared *
		(1.0 / mean_reversion_parameter_times_2) *
		(std::exp(mean_reversion_parameter_times_2 * this_time) - std::exp(mean_reversion_parameter_times_2 * prev_time));

	double scale = std::exp(-mean_reversion_parameter_times_2 * this_time);

	// becomes total variance of x at this_time
	return (*prev_variance_intermediate) * scale;
}

double TimeSliceProcess::calculate_variance(double prev_time, double this_time, double prev_volatility_squared, double* prev_variance_intermediate)
{

	*prev_variance_intermediate = (*prev_variance_intermediate) + prev_volatility_squared * (this_time - prev_time);
	return *prev_variance_intermediate;

}

gelib::d_vector TimeSliceProcess::volatilities_to_variances(double mean_reversion_param, gelib::d_vector const& times, gelib::d_vector const& process_volatilities)
{
	gelib::d_vector process_variances(process_volatilities.size());
	return volatilities_to_variances(mean_reversion_param, times, process_volatilities, process_variances);
}

gelib::d_vector const& TimeSliceProcess::volatilities_to_variances(double mean_reversion_param, gelib::d_vector const& times, gelib::d_vector const& process_volatilities, gelib::d_vector& process_variances)
{
	gelib::d_vector process_volatilities_squared(process_volatilities.size());
	for (size_t j = 0; j < process_volatilities_squared.size(); ++j)
	{
		process_volatilities_squared[j] = process_volatilities[j] * process_volatilities[j];
	}// TODO back to transform
	//std::transform(process_volatilities.begin(), process_volatilities.end(), process_volatilities_squared.begin(), [](double sigma) {return sigma * sigma; });
	return squared_volatilities_to_variances(mean_reversion_param, times, process_volatilities_squared, process_variances);
}


gelib::d_vector const& TimeSliceProcess::squared_volatilities_to_variances(double mean_reversion_param, gelib::d_vector const& times, gelib::d_vector const& process_volatilities_squared, gelib::d_vector& process_variances)
{
	double prev_time, this_time, time_increment_volatility_squared, prev_variance_intermediate;
	std::function<double(void)> variance_function;

	process_variances.resize(process_volatilities_squared.size());

	if (mean_reversion_param != 0.0)
	{
		double mean_reversion_param_times_2 = 2 * mean_reversion_param;
		variance_function = [&]()
		{ return calculate_variance(mean_reversion_param_times_2, prev_time, this_time, time_increment_volatility_squared, &prev_variance_intermediate); };

	}
	else
	{
		variance_function = [&]()
		{ return calculate_variance(prev_time, this_time, time_increment_volatility_squared, &prev_variance_intermediate); };
	}

	prev_time = 0.0;
	prev_variance_intermediate = 0.0;

	for (size_t i = 0; i < times.size(); ++i)
	{
		this_time = times[i];
		time_increment_volatility_squared = process_volatilities_squared[i];

		process_variances[i] = variance_function();

		prev_time = this_time;
	}

	return process_variances;
}

gelib::d_vector const& TimeSliceProcess::variances_to_squared_volatilites(double mean_reversion_param, gelib::d_vector const& times, gelib::d_vector const& process_variances, gelib::d_vector& process_volatilities_squared)
{
	// TODO implement general case
	assert(mean_reversion_param == 0.0);

	process_volatilities_squared.resize(process_variances.size());

	process_volatilities_squared[0] = process_variances[0] / times[0];

	for (size_t j = 1; j < process_variances.size(); ++j)
	{
		process_volatilities_squared[j] = (process_variances[j] - process_variances[j - 1]) / (times[j] - times[j - 1]);
	}

	return process_volatilities_squared;
}

void TimeSliceProcess::allow_refit()
{
	m_spot_fitted = false;
}

double TimeSliceProcess::get_option_in_arrears_price(double k) const
{
	return m_pricer->digital_option_in_arrears_price(k);
}

double TimeSliceProcess::get_option_in_arrears_price_inverse(double price) const
{
	return m_pricer->digital_option_in_arrears_price_inverse(price);
}

void TimeSliceProcess::calculate_spot()
{
	// calculate spot values by inverting J0 using the digital option in arrears price formula for the process in question. 
	std::transform(m_J0.begin(), m_J0.end(), m_spot_values.begin(), [&](double j0i) -> double { return m_pricer->digital_option_in_arrears_price_inverse(j0i); });
	
	// we have new spot values, so should be allowed a new fit
	m_spot_fitted = false;
}

void TimeSliceProcess::fit_spot(const std::string& extrapolation)
{
	if (m_spot_fitted)
		return;

	char c = gelib::PiecewisePoly1d::parse_extrapolation_string(extrapolation);

	// important here that the fit function is given a constant pointer to the grid or will be modified
	m_spot_fit->fit(get_grid()->data(), &m_spot_values[0], m_grid.size(), m_model->get_order(), c);

	m_spot_fitted = true;
}

void TimeSliceProcess::interpolate_spot(gelib::d_vector const& x, gelib::d_vector& f_x) const
{
	if (!m_spot_fitted)
		throw std::runtime_error("Call the spot fitting function before interpolating");

	m_spot_fit->evaluate(x, f_x);
}

double TimeSliceProcess::interpolate_spot(double x) const
{
	if (!m_spot_fitted)
		throw std::runtime_error("Call the spot fitting function before interpolating");

	return m_spot_fit->evaluate(&x);
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// construct a seqeuence of contiguous random processes
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

std::vector<std::unique_ptr<TimeSliceProcess>> TimeSliceProcess::construct(double grid_spread, size_t number_of_gridpoints, const MFMBase* model, double mean_reversion_param, gelib::d_vector const& times, gelib::d_vector const& process_volatilities_squared, gelib::d_vector const& process_variances, gelib::d_vector const& grid_centers)
{
	std::vector<std::unique_ptr<TimeSliceProcess>> processes;

	size_t n = times.size();

	processes.resize(n + 1);

	double mean_reversion_param_times_2 = mean_reversion_param * 2;
	double one_over_mean_reversion_param_times_2 = 1.0;
	double prev_time, this_time, next_time, conditional_means_factor, increment_variance;


	std::function<double(double)> calc_piece_variance;

	if (mean_reversion_param != 0.0)
	{
		one_over_mean_reversion_param_times_2 = 1.0 / mean_reversion_param_times_2;

		calc_piece_variance = [&](double volatility_squared)
		{
			return volatility_squared * one_over_mean_reversion_param_times_2 * (1.0 - std::exp(-mean_reversion_param_times_2 * (next_time - this_time)));
		};
	}
	else
	{

		calc_piece_variance = [&](double volatility_squared)
		{
			return volatility_squared * (next_time - this_time);
		};
	}

	/////////////////////////////////////////////////
	// process 0
	/////////////////////////////////////////////////

	this_time = 0.0;
	next_time = times[0];

	// note that grids at T0 share the same center as grid at T1
	processes[0] = TimeSliceProcess::construct(this_time, calc_piece_variance(process_volatilities_squared[0]), 0.0, 1.0, 0.0, 1, 0.0, model);

	/////////////////////////////////////////////////
	// process 1 to n - 1
	/////////////////////////////////////////////////

	for (size_t i = 1; i < n; ++i)
	{

		// set the variables for the process at time slice i
		prev_time = this_time;
		this_time = next_time;
		next_time = times[i];
		conditional_means_factor = std::exp(-mean_reversion_param * (next_time - this_time));
		increment_variance = calc_piece_variance(process_volatilities_squared[i]);

		// construct process 
		processes[i] = TimeSliceProcess::construct(this_time, increment_variance, process_variances[i - 1], grid_spread, grid_centers[i - 1], number_of_gridpoints, conditional_means_factor, model);

	}

	prev_time = this_time;
	this_time = next_time;
	conditional_means_factor = 0.0;
	increment_variance = 0.0;
	processes[n] = TimeSliceProcess::construct(this_time, increment_variance, process_variances[n - 1], grid_spread, grid_centers[n - 1], number_of_gridpoints, conditional_means_factor, model);

	return std::move(processes);
}

size_t TimeSliceProcess::total_grid_points() const
{
	return m_grid.size();
}

double TimeSliceProcess::get_conditional_variance() const
{
	return m_conditional_variance;
}

double TimeSliceProcess::get_conditional_standard_deviation() const
{
	return m_conditional_sigma;
}

double TimeSliceProcess::get_process_variance() const
{
	return m_variance;
}

double TimeSliceProcess::time() const
{
	return m_time;
}

double TimeSliceProcess::get_grid_center() const
{
	return m_grid_center;
}

double TimeSliceProcess::get_grid_spread_factor() const
{
	return m_grid_sd_spread_factor;
}

void TimeSliceProcess::get_grid(gelib::d_vector& out) const
{
	out = m_grid;
}

void TimeSliceProcess::get_conditional_means(gelib::d_vector& out) const
{
	out = m_conditional_grid_means;
}

const gelib::d_vector* TimeSliceProcess::get_conditional_means() const
{
	return &m_conditional_grid_means;
}

void TimeSliceProcess::get_J0_grid(gelib::d_vector& out) const
{
	out = m_J0;
}

gelib::d_vector* TimeSliceProcess::get_J0_grid()
{
	return &m_J0;
}

gelib::d_vector* TimeSliceProcess::get_spot_grid()
{
	return &m_spot_values;
}

void TimeSliceProcess::get_spot(gelib::d_vector& out) const
{
	out = m_spot_values;
}

const gelib::d_vector* TimeSliceProcess::get_spot() const
{
	return &m_spot_values;
}

double TimeSliceProcess::get_spot(size_t i) const
{
	return m_spot_values[i];
}

const gelib::d_vector* TimeSliceProcess::get_grid() const
{
	return &m_grid;
}
