//
// Created by quentin on 9/13/22.
//

#include "math.h"

//Polynomial Fit
#include <iostream>
#include <cmath>

using namespace Budget::Utilities;

PolynomialFunction Math::polynomialFit(int degree, int numberOfPoints, const double x[], const double y[]) {
	//Array that will store the values of sigma(xi),sigma(xi^2),sigma(xi^3)....sigma(xi^2numberOfPoints)
	double X[2 * degree + 1];
	for (int i = 0; i < 2 * degree + 1; i++) {
		X[i] = 0;
		for (int j = 0; j < numberOfPoints; j++) {
			//consecutive positions of the array will store numberOfPoints,sigma(xi),sigma(xi^2),sigma(xi^3)....sigma(xi^2n)
			X[i] = X[i] + pow(x[j], i);
		}
	}
	//B is the Normal matrix(augmented) that will store the equations, 'a' is for value of the final coefficients
	double B[degree + 1][degree + 2], a[degree + 1];
	for (int i = 0; i <= degree; i++) {
		for (int j = 0; j <= degree; j++) {
			//Build the Normal matrix by storing the corresponding coefficients at the right positions except the last column of the matrix
			B[i][j] = X[i + j];
		}
	}
	
	//Array to store the values of sigma(yi),sigma(xi*yi),sigma(xi^2*yi)...sigma(xi^degree*yi)
	double Y[degree + 1];
	for (int i = 0; i < degree + 1; i++) {
		Y[i] = 0;
		for (int j = 0; j < numberOfPoints; j++) {
			//consecutive positions will store sigma(yi),sigma(xi*yi),sigma(xi^2*yi)...sigma(xi^degree*yi)
			Y[i] = Y[i] + pow(x[j], i) * y[j];
		}
	}
	for (int i = 0; i <= degree; i++) {
		//load the values of Y as the last column of B(Normal Matrix but augmented)
		B[i][degree + 1] = Y[i];
	}
	
	//degree is made degree+1 because the Gaussian Elimination part below was for degree equations, but here degree is the degree of polynomial and for degree degree we get degree+1 equations
	degree = degree + 1;
	//From now Gaussian Elimination starts(can be ignored) to solve the set of linear equations (Pivotisation)
	for (int i = 0; i < degree; i++) {
		for (int k = i + 1; k < degree; k++) {
			if (B[i][i] < B[k][i]) {
				for (int j = 0; j <= degree; j++) {
					double temp = B[i][j];
					B[i][j] = B[k][j];
					B[k][j] = temp;
				}
			}
		}
	}
	
	//loop to perform the gauss elimination
	for (int i = 0; i < degree - 1; i++) {
		for (int k = i + 1; k < degree; k++) {
			double t = B[k][i] / B[i][i];
			for (int j = 0; j <= degree; j++) {
				//make the elements below the pivot elements equal to zero or elimnate the variables
				B[k][j] = B[k][j] - t * B[i][j];
			}
		}
	}
	
	//back-substitution
	for (int i = degree - 1; i >= 0; i--) {
		//x is an array whose values correspond to the values of x,y,z..
		//make the variable to be calculated equal to the rhs of the last equation
		a[i] = B[i][degree];
		for (int j = 0; j < degree; j++) {
			//then subtract all the lhs values except the coefficient of the variable whose value is being calculated
			if (j != i) {
				a[i] = a[i] - B[i][j] * a[j];
			}
		}
		//now finally divide the rhs by the coefficient of the variable to be calculated
		a[i] = a[i] / B[i][i];
	}
	std::cout << "\nThe values of the coefficients are as follows:\n";
	for (int i = 0; i < degree; i++)
		std::cout << "x^" << i << "=" << a[i] << std::endl;            // Print the values of x^0,x^1,x^2,x^3,....
	std::cout << "\nHence the fitted Polynomial is given by:\ny=";
	for (int i = 0; i < degree; i++)
		std::cout << " + (" << a[i] << ")" << "x^" << i;
	std::cout << "\n";
	
	std::vector<double> va(a, a + degree);
	return PolynomialFunction(va);
}