#include<cmath>
#include<vector>
#include<map>
#include"Date.h"
#include"TermStructure.h"
#include"TermStructureCubicSpline.h"
#include"Tridiagonal.h"

using namespace std;

TermStructureCubicSpline::TermStructureCubicSpline() {};

TermStructureCubicSpline::TermStructureCubicSpline(vector<Date> dates,vector<double> rates,double daycnt):TermStructure(dates,rates,daycnt)
{	
	map<double,double> F=TermStructure::getFunction();

	vector<double> xValues;
	vector<double> yValues;
	vector<double> xDifferences;
	vector<double> yDifferences;
	vector<double> subDiag;
	vector<double> diag;
	vector<double> superDiag;
	vector<double> slopes;

	map<double,double>::const_iterator node=F.begin();

	while(node!=F.end()){
		xValues.push_back(node->first);
		yValues.push_back(node->second);
		node++;
	}

	for(size_t k=0;k<xValues.size();k++){
		xDifferences.push_back(xValues[k+1]-xValues[k]);
		yDifferences.push_back(yValues[k+1]-yValues[k]);
	};

	for(size_t k=0;k<=xDifferences.size();k++){
		double slope=yDifferences[k]/xDifferences[k];
		slopes.push_back(slope);
		};

	for(size_t k=0;k<=xDifferences.size()-2;k++){

		diag.push_back((xDifferences[k+1]+xDifferences[k])/3);
		subDiag.push_back(xDifferences[k]/6);
		superDiag.push_back(xDifferences[k+1]/6);
	};

	A.setSubDiagonal(subDiag);
	A.setDiagonal(diag);
	A.setSuperDiagonal(superDiag);

	vector<double> rhs;

	for(size_t k=1;k<slopes.size();k++)
	{
		rhs.push_back(slopes[k+1]-slopes[k]);
	}
	
	doubleDerivatives.push_back(0.0);
	doubleDerivatives.insert(doubleDerivatives.end(),A.solver(rhs).begin(),A.solver(rhs).end());
	doubleDerivatives.push_back(0.0);
}


double TermStructureCubicSpline::buildCurve(const Date &dt){

	map<double,double> F=TermStructure::getFunction();

	map<double,pair<double,double> > Spline;
	map<double,double>::const_iterator t=F.begin();

	for(size_t k=0;k<=doubleDerivatives.size();k++)
	{
		Spline.insert(make_pair(t->first,make_pair(t->second,doubleDerivatives[k])));
		t++;
	};
	
	int start=getStart().daysSinceStart();
	double timeToDate= (dt.daysSinceStart()-start)/getDayCount();
	

	map<double,pair<double,double> >::const_iterator node=Spline.begin();

	while(node!=Spline.end() && (node->first)<timeToDate)
	{ node++;}



	double A = (node->first-timeToDate)/(node->first-(node--)->first);
	double B=1-A;
	double C=(pow(A,3)-A)*pow((node->first-(node--)->first),2)/6;
	double D=(pow(B,3)-B)*pow((node->first-(node--)->first),2)/6;


	double y=A*(node->second).first+B*((node++)->second).first+C*(node->second).second+D*((node++)->second).second;

	return y;

};