/* This class is used to generate uniformely distributed points
 * for the forward rates vs maturity and volatility vs maturity
 * curves, given few points.
*/


#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <math.h>
#include "Curve.cpp"
#include "Point.cpp"

class Spline{
    public: Spline(){}
    /** Cubic Spline fit to data, based on Press et al, adapted by
    Cristian Constantin Bordeianu, Romania 
    input array x[n], y[n] represents tabulation function y(x)
    with x0 < x1 ... < x(n-1). n = # of tabulated points
    output yout for given xout (here xout via loop at end)
    yp1 and ypn: 1st derivatives at endpoints, evaluated internally
    y2[n]  is array of second derivatives
    (setting yp1 or  ypn >0.99e30 produces natural spline)**/
    // input data, here scaled for graph: 0 < x < 1, -100 y < 100

    //input x[] and y[] data, xout is the x value we want to evaluate
    // yout is the spline version of F(xout)
    public: double spline_one(vector<double> x1,
				              vector<double> y1,
							  double xout)
    {
        double* x = new double[x1.size()];
        double* y = new double[x1.size()];
        double yout=0.0; 
        double y2[9];
        int n,np; 
		for(int i=0;i<x1.size();i++){
			x[i]=0.0;
			y[i]=0.0;
		}
        n = x1.size(); 
        np = 15; 
        int klo,khi,k;
        double h,b,a;
        double p,qn,sig,un,yp1,ypn;
        double* u = new double[n];
        
        for (i=0;i<x1.size();i++){
            x[i] = x1.at(i);
            y[i] = y1.at(i);
        }
        if(x1.size()>3){
            yp1  =  (y[1]-y[0])/(x[1]-x[0]);  // 1st deriv at initial point
            yp1 = yp1 -(y[2]-y[1])/(x[2]-x[1]); 
            yp1 = yp1 +(y[2]-y[0])/(x[2]-x[0]);
            ypn = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]) // 1st deriv at end point
                -(y[n-2]-y[n-3])/(x[n-2]-x[n-3])
                +(y[n-1]-y[n-3])/(x[n-1]-x[n-3]);
            if (yp1 > 0.99e30) y2[0] = u[0] = 0.0; //natural spline test
            else {
                y2[0]  =  (-0.5);
                u[0] = (3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
            }
            for (int i=1; i<=n-2; i++)   // decomposition loop; y2, u are temps
            {
                sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
                p = sig*y2[i-1]+2.0;
                y2[i] = (sig-1.0)/p;
                u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
                u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
            }
            if (ypn > 0.99e30) qn = un = 0.0;   //test for natural
            else                       // else evaluate second derivative
            { 
                qn = 0.5;
                un = ((3.0)/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
            }
            y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2]+1.0);
            for (k = n-2;k>= 0;k--) {
                y2[k] = y2[k]*y2[k+1]+u[k];
            }    // back substitution
            // splint (initialization) ends
            // Parameters determined, Begin *spline* fit
            klo = 0;                //Bisection algor to find place in table
            khi = n-1;              // klo, khi bracket xout value 
            while (khi-klo > 1) {
                k = (khi+klo) >> 1;
                if (x[k] > xout) khi = k;
                else klo = k;     
            }
            h = x[khi]-x[klo];
            if (x[k] > xout) khi = k;
            else klo = k;
            h = x[khi]-x[klo];
            a = (x[khi]-xout)/h;
            b = (xout-x[klo])/h;
            yout = (a*y[klo]+b*y[khi]+((a*a*a-a)*y2[klo]
                   +(b*b*b-b)*y2[khi])*(h*h)/6.0);
        }
        return yout;
    }
    public: Curve spline(vector<double> x,vector<double> y,
				         int n,double maturity) {
        // TODO code application logic here
        Curve curve;
        double yield=0;
        double dt=maturity/(1.0*n);
        double t=0.0;
        for (int i=0;i<(n+1);i++){
            double yout=spline_one(x,y,t);
            Point* a = new Point(t,yout);
            curve.push_back(a);
            t=t+dt;
        }        
        return curve;
    }
    public: double* spline_array(vector<double> x,vector<double> y,int n_points,double maturity) {
        // TODO code application logic here
        double* curve=new double[n_points+1];
        double yield=0;
        double dt=maturity/(1.0*n_points);
        double t=0.0;
        for (int i=0;i<=n_points;i++){
            curve[i] = spline_one(x,y,t);
            t=t+dt;
        }        
        return curve;
    }
    public: double* spline_time_array(int n_points,double maturity) {
        // TODO code application logic here
        double* curve = new double[n_points+1];
        double dt=maturity/(1.0*n_points);
        double t=0.0;
        for (int i=0;i<=n_points;i++){
            curve[i] = t;
            t=t+dt;
        }        
        return curve;
    } 
    public: double* spline_prime_array(vector<double> x,vector<double> y,int n_points,double maturity) {
        // TODO code application logic here
        double* curve = new double[n_points+1];
        double yield=0;
        double dt=maturity/(1.0*n_points);
        double dt1=dt/10.0;
        double t=0.0;
        for (int i=0;i<=n_points;i++){
            double y2=spline_one(x,y,t+dt1);
            double y1=spline_one(x,y,t);
            double der = (y2-y1)/dt1;
            curve[i] = der;
            t=t+dt;
        }        
        return curve;
    } 
    public: Curve spline_prime(vector<double> x,vector<double> y,int n_points,double maturity) {
        // TODO code application logic here
        Curve curve;
        double yield=0;
        double dt=maturity/(1.0*n_points);
        double dt1=dt/10.0;
        double t=0.0;
        for (int i=0;i<n_points+1;i++){
            double y2=spline_one(x,y,t+dt1);
            double y1=spline_one(x,y,t);
            double der = (y2-y1)/dt1;
            Point* a = new Point(t,der);
            curve.push_back(a);
            t=t+dt;
        }        
        return curve;
    }            
};
/*void main() {
	vector<double> x;
	x.push_back(0.);x.push_back(0.12);x.push_back(0.25);x.push_back(0.37); 
    vector<double>y;
	y.push_back(10.6);y.push_back(16.0);y.push_back(45.0);y.push_back(83.5); 
    Spline s;
    double* yout = s.spline_array(x,y,5,1.0);
}
*/