FunctionMatrixSolve
// 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 Matrix Equation solution for Ax = B, finds the value of x.
library(title='FunctionMatrixSolve')

//reference:
//      https://introcs.cs.princeton.edu/java/95linear/Matrix.java.html


// @function helper to create and return a random M-by-N matrix with values between 0 and 1
random (int M, int N) => //{
    matrix<float> _A = matrix.new<float>(M, N)
    for _i = 0 to M-1
        for _j = 0 to N-1
            matrix.set(_A, _i, _j, math.random())
    _A
//}


// @function Helper to swap rows i and j
swap_rows (matrix<float> M, int i, int j) => //{
    int _cols = matrix.columns(M)
    for _k = 0 to _cols - 1
        _t = matrix.get(M, i, _k)
        matrix.set(M, i, _k, matrix.get(M, j, _k))
        matrix.set(M, j, _k, _t)
//}

// @function Solves Matrix Equation for Ax = B, finds value for x.
// @param A matrix<float>, Square matrix with data values.
// @param B matrix<float>, One column matrix with data values.
// @returns matrix<float> with X, x = A^-1 b, assuming A is square and has full rank
// https://introcs.cs.princeton.edu/java/95linear/Matrix.java.html
export solve (matrix<float> A, matrix<float> B) => //{
    int _Arows = matrix.rows(A)
    int _Acols = matrix.columns(A)
    int _Brows = matrix.rows(B)
    int _Bcols = matrix.columns(B)
    //
    if (_Arows != _Acols or _Brows != _Acols or _Bcols != 1)
        runtime.error("Illegal matrix dimensions.")
    //
    // create copies of the data
    matrix<float> _A = matrix.copy(A)
    matrix<float> _B = matrix.copy(B)

    // Gaussian elimination with partial pivoting
    for _i = 0 to _Acols-1
        // find pivot row and swap
        int _max = _i
        int _j = _i + 1
        while _j < _Acols
            float _Aji = matrix.get(_A, _j, _i)
            float _Ami = matrix.get(_A, _max, _i)
            if (math.abs(_Aji) > math.abs(_Ami))
                _max := _j
            _j += 1
        swap_rows(_A, _i, _max)
        swap_rows(_B, _i, _max)
        // singular
        float _Aii = matrix.get(_A, _i, _i)
        if (_Aii == 0.0) 
            runtime.error("Matrix is singular.")
        // pivot within b
        _j := _i + 1
        while _j < _Acols
            float _Bj0 = matrix.get(_B, _j, 0) - matrix.get(_B, _i, 0) * matrix.get(_A, _j, _i) / matrix.get(_A, _i, _i)
            matrix.set(_B, _j, 0, _Bj0)
            _j += 1
        // pivot within A
        _j := _i + 1
        while _j < _Acols
            float _m = matrix.get(_A, _j, _i) / matrix.get(_A, _i, _i)
            for _k = _i + 1 to _Acols-1
                float _Ajk = matrix.get(_A, _j, _k)
                matrix.set(_A, _j, _k, _Ajk - matrix.get(_A, _i, _k) * _m)
            matrix.set(_A, _j, _i, 0.0)
            _j += 1
    // back substitution
    matrix<float> _X = matrix.new<float>(_Acols, 1)
    for _j = _Acols-1 to 0
        float _t = 0.0
        _k = _j + 1
        while _k < _Acols
            _t += matrix.get(_A, _j, _k) * matrix.get(_X, _k, 0)
            _k += 1
        matrix.set(_X, _j, 0, (matrix.get(_B, _j, 0) - _t) / matrix.get(_A, _j, _j))
    _X
//}

// A = random(5, 5)
// B = random(5, 1)

A = matrix.new<float>(0, 5)
matrix.add_row(A, 0, array.from( 1.0,  0.0,  0.0,  0.0, 0.0))
matrix.add_row(A, 1, array.from(-1.0,  1.0,  0.0,  0.0, 0.0))
matrix.add_row(A, 2, array.from(-1.0, -1.0,  1.0,  0.0, 0.0))
matrix.add_row(A, 3, array.from(-1.0, -1.0, -1.0,  1.0, 0.0))
matrix.add_row(A, 4, array.from(-1.0, -1.0, -1.0, -1.0, 1.0))
B = matrix.new<float>(5, 0)
matrix.add_col(B, 0, array.from(0.5, 0.0, 0.0, 0.0, 0.5))
// // should be [0.5, 0.5, 1, 2, 4.5] // https://www.mathstools.com/section/main/system_equations_solver#.Yj9zcDzLe2I

// A = matrix.new<float>(0, 5)
// matrix.add_row(A, 0, array.from(1.4350, 0.6996, 1.5083, 1.1921, 1.0485 ))
// matrix.add_row(A, 1, array.from(0.6996, 1.2069, 1.6564, 1.3865, 0.7041 ))
// matrix.add_row(A, 2, array.from(1.5083, 1.6564, 3.2516, 2.5462, 1.9850 ))
// matrix.add_row(A, 3, array.from(1.1921, 1.3865, 2.5462, 2.5376, 1.4350 ))
// matrix.add_row(A, 4, array.from(1.0485, 0.7041, 1.9850, 1.4350, 1.5252 ))
// B = matrix.new<float>(5, 0)
// matrix.add_col(B, 0, array.from(0.9636, 0.2268, 0.2594, 0.3278, 0.1377))
// should be [1.2279, 0.1246, -0.3779, 0.0943, -0.4083 ] // 
// https://www.online-java.com/ 
//      with files: https://introcs.cs.princeton.edu/java/stdlib/StdOut.java.html
//                  https://introcs.cs.princeton.edu/java/95linear/Matrix.java.html

X = solve(A, B) 

table T = table.new(position.bottom_right, 5, 2)
table.cell(T, 0, 0, 'A', bgcolor=#ffffff)
table.cell(T, 0, 1, str.tostring(A, '#.0000'), bgcolor=#ffffff)
// table.cell(T, 1, 1, '*', bgcolor=#ffffff)
table.cell(T, 4, 0, 'X', bgcolor=#ffffff)
table.cell(T, 4, 1, str.tostring(X, '#.0000'), bgcolor=#ffffff)
// table.cell(T, 3, 1, '=', bgcolor=#ffffff)
table.cell(T, 2, 0, 'B', bgcolor=#ffffff)
table.cell(T, 2, 1, str.tostring(B, '#.0000'), bgcolor=#ffffff)

