FunctionPolynomialFit
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © RicardoSantos

//@version=5

// @description Performs Polynomial Regression fit to data.
// In statistics, polynomial regression is a form of regression analysis in which 
// the relationship between the independent variable x and the dependent variable 
// y is modelled as an nth degree polynomial in x. 
// reference: 
//      https://en.wikipedia.org/wiki/Polynomial_regression
//      https://www.bragitoff.com/2018/06/polynomial-fitting-c-program/
library(title='FunctionPolynomialFit', overlay=true)

// @function Perform Gauss-Elimination and returns the Upper triangular matrix and solution of equations.
// @param A float matrix, data samples.
// @param m int, defval=na, number of rows.
// @param n int, defval=na, number of columns.
// @returns float array with coefficients.
export gauss_elimination (matrix<float> A, int m=na, int n=na) => //{
    int _m = na(m) ? matrix.rows(A)    : m      , int _pm = _m - 1
    int _n = na(n) ? matrix.columns(A) : n      , int _pn = _n - 1
    float[] _X = array.new_float(_m + 1, 0.0)
    //
    for _i = 0 to _pm-1
        // partial pivoting:
        for _k = _i + 1 to _pm
            // if diagonal element(absolute value) is smaller than any of the terms below it
            if math.abs(matrix.get(A, _i, _i)) < math.abs(matrix.get(A, _k, _i))
                for _j = 0 to _pn
                    float _Aij = matrix.get(A, _i, _j)
                    matrix.set(A, _i, _j, matrix.get(A, _k, _j))
                    matrix.set(A, _k, _j, _Aij)
        // begin Gauss Elimination:
        for _k = _i + 1 to _pm
            float _term = matrix.get(A, _k, _i) / matrix.get(A, _i, _i)
            for _j = 0 to _pn
                matrix.set(A, _k, _j, matrix.get(A, _k, _j) - _term * matrix.get(A, _i, _j))
    // begin Back-Substitution:
    for _i = _pm to 0
        array.set(_X, _i, matrix.get(A, _i, _pn))
        for _j = _i + 1 to _pn
            array.set(_X, _i, array.get(_X, _i) - matrix.get(A, _i, _j) * array.get(_X, _j))
        array.set(_X, _i, array.get(_X, _i) / matrix.get(A, _i, _i))
    _X
//}

// @function Fits a polynomial of a degree to (x, y) points.
// @param X float array, data sample x point.
// @param Y float array, data sample y point.
// @param degree int, defval=2, degree of the polynomial.
// @returns float array with coefficients.
// note:
//      p(x) = p[0] * x**deg + ... + p[deg]
export polyfit (float[] X, float[] Y, int degree=2) => //{
    int _m = array.size(X)        , int _pm = _m - 1
    //
    float[] _X1 = array.new_float(1 + degree * 2, 0.0)
    for _i = 0 to degree * 2
        for _j = 0 to _pm
            array.set(_X1, _i, array.get(_X1, _i) + math.pow(array.get(X, _j), _i))

    //the normal augmented matrix
    matrix<float> _B = matrix.new<float>(degree + 1, degree + 2)
    // rhs
    float[] _Y1 = array.new_float(degree + 1, 0.0)
    for _i = 0 to degree
        for _j = 0 to _pm
            array.set(_Y1, _i, array.get(_Y1, _i) + math.pow(array.get(X, _j), _i) * array.get(Y, _j))
    for _i = 0 to degree
        for _j = 0 to degree
            matrix.set(_B, _i, _j, array.get(_X1, _i + _j))
    for _i = 0 to degree
        matrix.set(_B, _i, degree + 1, array.get(_Y1, _i))
    
    // _str = 'The polynomial fit is given by the equation:\n'
    // _str += str.tostring(_B) + '\n'
    _A = gauss_elimination(_B)
    // _str += str.tostring(_A) + '\n'
// usage:
// if barstate.isfirst
//     _x = array.from(1, 2, 3, 4, 5.0)
//     _y = array.from(1.0, 7.9, 26.99, 64, 125)
//     label.new(bar_index, 0.0, str.tostring(polyfit(_x, _y, 3)))
int length = input.int(10)
int degree = input.int(3)
var x = array.new_float(length, float(bar_index)), array.shift(x), array.push(x, float(bar_index))
var y = array.new_float(length, float(close)), array.shift(y), array.push(y, float(close))
coeffs = polyfit(x, y, degree)
float solve = 0.0
for _i = 0 to array.size(coeffs)-1
    solve += array.get(coeffs, _i) * math.pow(bar_index, _i)
plot(series=solve, color=color.teal)
//}

// @function interpolate the y position at the provided x.
// @param coeffs float array, coefficients of the polynomial.
// @param x float, position x to estimate y.
// @returns float.
export interpolate (float[] coeffs, float x) => //{
    float _y = 0.0
    for _i = 0 to array.size(coeffs) - 1
        _y += array.get(coeffs, _i) * math.pow(x, _i)
    _y
// usage:
if barstate.islast
    _py = interpolate(coeffs, bar_index[length])
    for _i = 1 to length
        _y = interpolate(coeffs, bar_index[length - _i])
        line.new(bar_index[length-_i]-1, _py, bar_index[length-_i], _y, color=color.yellow, width=2)
        _py := _y
    // forecast
    for _i = 1 to 5
        _y = interpolate(coeffs, bar_index + _i)
        line.new(bar_index + _i - 1, _py, bar_index + _i, _y, color=color.yellow, width=2)
        _py := _y
//}


