#include <fstream.h>
#include <math.h>
#include "matrix.cpp"

void read_datafile(char *filename, int n, long double *x, long double *y);
void write_results(char *filename, int n, long double *x, long double *y,
                   long double *dydx, long double *d2ydx2);
void fit_poly(int n, int degree, long double *x, long double *y, long double *coeff);
void create_output(int outpnts, int degree, long double *dydx, long double *d2ydx2,
			long double *coeff, long double *xnu, long double *ynu, long double xo, long double xf);
void fit_poly2(int n, int degree, long double *x, long double *y, long double *coeff);
void main(void)
{
   char output_file[50];
   char input_file[50];
   int n, outpnts, degree;
   long double *coeff,*x,*y,*dydx,*d2ydx2,*xnu,*ynu;

   cout << "Enter input filename: ";
   cin >> input_file;
   cout << "Enter output filename: ";
   cin >> output_file;
   cout << "Enter the number of data points to read: ";
   cin >> n;
   cout << "Enter the number of points to output: ";
   cin >> outpnts;
	cout << "Enter degree of polynomial: ";
	cin >> degree;

	x = new long double[n];
   y = new long double[n];
   dydx = new long double[outpnts];
   d2ydx2 = new long double[outpnts];
   xnu = new long double[outpnts];
   ynu = new long double[outpnts];
	coeff = new long double[degree+1];

   read_datafile(input_file,n,x,y);
   fit_poly2(n,degree,x,y,coeff);
   for (int i=0; i<=degree; i++)
   	cout << "Coefficient " << i << ": " << coeff[i] << endl;
   create_output(outpnts, degree, dydx, d2ydx2,coeff, xnu, ynu, x[0], x[n-1]);
   write_results(output_file, outpnts, xnu, ynu, dydx, d2ydx2);

   delete []x;
   delete []y;
   delete []dydx;
   delete []d2ydx2;
   delete []xnu;
   delete []ynu;
   delete []coeff;
	cout << "Pausing. Enter integer to exit: "; int wait; cin >> wait;
}

/*
   n = number data points
   degree = degree of polynomial to fit
   x = array of x values
   y = array of y values
   coeff = array for coefficients of fit polynomial (degree + 1 coefficients
				   are returned)
   This function fits a polynomial to the given data with a least squares
   method, using Gauss-Jordan elimination to solve the linear equations.
   */

void fit_poly2(int n, int degree, long double *x, long double *y, long double *coeff)
{
   // Consider sums a row matrix of degree+1 elements containing yk * xk^i
   long double *sums;
   sums = new long double[degree+1];

   for (int i=0; i<=degree; i++) // Sum over degree of term
   	{
      sums[i] = 0.0;
      for (int k=0; k<n; k++) sums[i] += y[k]*powl(x[k],i); // Sum over data
      }

   long double sum;
   matrix delta(degree+1, degree+1);
   for (int i=0; i<=degree; i++)       // Sum over degree in each dimension
   	{
      for (int j=0; j<=degree; j++)
      	{
         sum = 0;
         for (int k=0;k<n; k++) // sum over all data points
         	{
            sum += powl(x[k],i)*powl(x[k],j);
            }
         delta.set(i,j,sum);
         }
      }

   // The coefficients satisfy the matrix equation coeff = sums * delta^-1

   delta.inverse(); // takes inverse of matrix using Gauss-Jordan Elimination
   delta.times_row(sums, coeff); // returns sums * delta^-1 in coeff array
   delete [] sums;

}

/*
   n = number data points
   degree = degree of polynomial to fit
   x = array of x values
   y = array of y values
   coeff = array for coefficients of fit polynomial (degree + 1 coefficients
				   are returned)
   This function fits a polynomial to the given data using Cramer's rule to
   solve the linear equations.  The determinates are solved using a recursive
   cofactor method.
   */
void fit_poly(int n, int degree, long double *x, long double *y, long double *coeff)
	{
   // The coefficients are given by det(co-matrix) / det(delta matrix)
   // The equations to be solved are represented as matrices:
   //    sums = coeff * delta

   long double sum;
   matrix delta(degree+1, degree+1);
   for (int i=0; i<=degree; i++)
   	{
      for (int j=0; j<=degree; j++)
      	{
         sum = 0;
         for (int k=0;k<n; k++) // sum over all data points
         	{
            sum += powl(x[k],i)*powl(x[k],j);
            }
         delta.set(i,j,sum);
         }
      }

   // Create the co-matrices equal to delta originally
   matrix **co_mats;
   co_mats = new matrix*[degree+1];
   for (int i=0; i<=degree; i++) // copy delta into coeff. matrices
   	{
      co_mats[i] = new matrix(degree+1, degree+1);
      co_mats[i]->equals(&delta);
      }

   // Create the sums

   long double *sums;
   sums = new long double[degree+1];

   for (int i=0; i<=degree; i++) // create sums
   	{
      sums[i] = 0.0;
      for (int k=0; k<n; k++) sums[i] += y[k]*powl(x[k],i);
      }

   for (int i=0; i<=degree; i++) // INSERT sums in proper cols of co-matrices
   	{
      for (int j=0; j<=degree; j++) // setting a col
	   	{
         co_mats[i]->set(j,i,sums[j]);
         }
      }

   long double det_delta; // Compute determinant of delata a single time
   det_delta = delta.det();

   // Compute coeficients... display a message since the algorithm can increase
   // in computational time to unreasonable lengths very quickly.
   for (int i=0; i<=degree; i++)
   	{
      coeff[i] = co_mats[i]->det() / det_delta;
      cout << "Finsihed coeff[" << i << "]" << endl;
      delete co_mats[i];
      }

   delete [] sums;
   delete [] co_mats;

   }

/*    Opens the text file named "filename" and reads n number of x,y data points
    where the input is tab delimited. If the end of the file is reached it
    will not read the full amount of data and display an error message.   */
void read_datafile(char *filename, int n, long double *x, long double *y)
	{
   fstream input_file;
   long double xpoint,ypoint;
   input_file.open(filename, ios::in);

   for(int i=0; ( i<n ) && ( input_file >> xpoint)  ; i++)
   	{

      input_file >> ypoint;
      x[i] = xpoint;
      y[i] = ypoint;
      }
   input_file.close();

	}

/*    Opens text file named "filename" and outputs n datapoints in the format:
    x tab y tab dydx tab d2ydx2 endline.   */
void write_results(char *filename, int n, long double *x, long double *y, long double *dydx,
                   long double *d2ydx2)
	{
   fstream output_file;
   output_file.open(filename, ios::out);
   if(!output_file) {cerr << "Error opening output file."; return;}
   for(int i=0; i<n; i++)
   	{
      output_file << x[i] << '\t' << y[i] << '\t' << dydx[i] << '\t'
                  << d2ydx2[i] << endl;
      }
   output_file.close();

	}


   /*   outpnts = number of points to evaluate on the polynomial

        degree = degree of polynomial
    outpnts number of points are computed for xnu, ynu, dydx, d2ydx2.
    coeff = array of (degree+1) coefficients for the polynomial
    The xnu values are evenly spaced from xo to xf, with points on the
       boundaries.   */

void create_output(int outpnts, int degree, long double *dydx,
	long double *d2ydx2, long double *coeff, long double *xnu, long double *ynu,
   long double xo, long double xf)
  	{
  	long double f(int degree, long double *coeff, long double x);
	long double dfdy(int degree, long double *coeff, long double x);
   long double d2fdy2(int degree, long double *coeff, long double x);

   long double dx;
  	dx = (xf - xo) / (long double)(outpnts-1);

   for (int i=0;i<outpnts; i++)
     	{
      xnu[i] = xo + (long double)i*dx;
      ynu[i] = f(degree, coeff, xnu[i]);
      dydx[i] = dfdy(degree, coeff, xnu[i]);
      d2ydx2[i] = d2fdy2(degree, coeff, xnu[i]);
      }

   }

   /* These 3 functions evaulate the polynomial of "degree" (or its derivatives)
      with the (degree+1) coefficients stored in coeff at the value x.
   */
long double f(int degree, long double *coeff, long double x)
	{
   long double val = 0.0;
   for (int i=0;i<=degree;i++)
   	val += coeff[i]*powl(x,i);
   return val;
   }

long double dfdy(int degree, long double *coeff, long double x)
	{
   long double val = 0.0;
   for (int i=0;i<=degree;i++)
   	{
	   if (i != 0)
   		val += (long double)i*coeff[i]*powl(x,i-1);
      }
   return val;
   }

long double d2fdy2(int degree, long double *coeff, long double x)
	{
   long double val = 0.0;
   for (int i=0;i<=degree;i++)
   	{
	   if (i != 0 && i-1 != 0)
   		val += (long double)(i*(i-1))*coeff[i]*powl(x,i-2);
      }
   return val;
   }