-- Available @ https://fxcodebase.com/code/viewtopic.php?f=17&t=74958

--+------------------------------------------------------------------------------------------------+
--|                                                            Copyright © 2024, Gehtsoft USA LLC  | 
--|                                                                         http://fxcodebase.com  |
--+------------------------------------------------------------------------------------------------+
--|                                                                   Developed by : Mario Jemic   |                    
--|                                                                       mario.jemic@gmail.com    |
--|                                               https://appliedmachinelearning.systems/contact/  | 
--+------------------------------------------------------------------------------------------------+

--+------------------------------------------------------------------------------------------------+
--|                                           Our work would not be possible without your support. |
--+------------------------------------------------------------------------------------------------+
--|                                                               Paypal:  https://goo.gl/9Rj74e   |
--|                                                             Patreon :  http://tiny.cc/1ybwxz   |  
--|                                                      Buy Me a Coffee:  http://tiny.cc/bj7vxz   | 
--+------------------------------------------------------------------------------------------------+



-- Indicator profile initialization routine
-- Defines indicator profile properties and indicator parameters
-- TODO: Add minimal and maximal value of numeric parameters and default color of the streams
function Init()
    indicator:name("Polynomial differentiator");
    indicator:description("");
    indicator:requiredSource(core.Tick);
    indicator:type(core.Oscillator);

  
 	indicator.parameters:addGroup("Calculation");	
    indicator.parameters:addInteger("Period", "RSI Period", "", 14, 1, 2000);   
    indicator.parameters:addInteger("Degree", "Polynomial Order", "", 2, 1, 2000); 
    indicator.parameters:addInteger("Window", "Frame Length (must be greater than Order)", "", 18, 2, 2000);	 
    indicator.parameters:addGroup("Line Style");	
    indicator.parameters:addInteger("width", "Line width", "", 1, 1, 5);
    indicator.parameters:addInteger("style", "Line style", "", core.LINE_SOLID);
    indicator.parameters:setFlag("style", core.FLAG_LINE_STYLE);
	
	indicator.parameters:addColor("color1", "Up Trend Bar Color", "", core.rgb(0, 255, 0)); 
	indicator.parameters:addColor("color2", "Down Bar in Up Trend Color", "", core.rgb(0, 200, 0));
    indicator.parameters:addColor("color3", "Up Bar in Down Trend Color", "", core.rgb(255, 0, 0));
    indicator.parameters:addColor("color4", "Down Trend Bar Color", "", core.rgb(200, 0, 0));
   	
end

-- Indicator instance initialization routine
-- Processes indicator parameters and creates output streams
-- TODO: Refine the first period calculation for each of the output streams.
-- TODO: Calculate all constants, create instances all subsequent indicators and load all required libraries
-- Parameters block
 
	
local first;
local source = nil;
local Period, Degree,  Window, SignalLength, IsSignal; 
local Indicator;
	
-- Routine
 function Prepare(nameOnly)   
 
    
	Period=instance.parameters.Period;
	Degree=instance.parameters.Degree;
	Window=instance.parameters.Window; 
	source = instance.source
 
    local name = profile:id() .. "(" ..  instance.source:name().. "," ..  Period   .. "," .. Degree .. "," .. Window   .. ")";
    instance:name(name); 


    if   (nameOnly) then
        return;
    end


	
	RSI= core.indicators:create("RSI", source, Period);
	first=RSI.DATA:first() ; 
	 

  
    PolynomialDifferentiator = instance:addStream("PolynomialDifferentiator", core.Bar, name, "PolynomialDifferentiator", instance.parameters.color1, first );
    PolynomialDifferentiator:setPrecision(math.max(2, instance.source:getPrecision())); 
    PolynomialDifferentiator:addLevel(0);	
	
	 
 
end


function Update(period, mode)

	RSI:update(mode); 

	if period <= first
	or  not source:hasData(period) 
	then
	return;
	end
	  
    PolynomialDifferentiator[period] = Differentiator(RSI.DATA, Window, Degree, period) 
	if PolynomialDifferentiator[period] > 0 then
	    if PolynomialDifferentiator[period]>PolynomialDifferentiator[period-1] then
		PolynomialDifferentiator:setColor(period, instance.parameters.color1);	
        else		
		PolynomialDifferentiator:setColor(period, instance.parameters.color2);
		end
	else				
        if PolynomialDifferentiator[period]>PolynomialDifferentiator[period-1] then	
		PolynomialDifferentiator:setColor(period, instance.parameters.color3);
		else
		PolynomialDifferentiator:setColor(period, instance.parameters.color4);
		end
    end	
	
	 
	
end

-- Subroutines
function matrix_get(_A,_i,_j,_nrows)
    return _A[_i+_nrows*_j]
end

function matrix_set(_A,_value,_i,_j,_nrows)
    _A[_i+_nrows*_j] = _value
end
 

function transpose(_A,_nrows,_ncolumns)
    local _AT = {}
	for i= 1, _nrows*_ncolumns  do
	_AT[i]=0;
	end	
	
    for i = 0, _nrows-1 do
        for j = 0, _ncolumns-1 do
            matrix_set(_AT, matrix_get(_A,i,j,_nrows),j,i,_ncolumns)
        end
    end
    return _AT
end


function multiply(_A,_B,_nrowsA,_ncolumnsA,_ncolumnsB)
    local _C = {}
	for i= 1, _nrowsA*_ncolumnsB  do
	_C[i]=0;
	end		
	
    local _nrowsB = _ncolumnsA
    local elementC = 0.0
    for i = 0, _nrowsA-1 do
        for j = 0, _ncolumnsB-1 do
            elementC = 0
            for k = 0, _ncolumnsA-1 do
                elementC = elementC + matrix_get(_A,i,k,_nrowsA)*matrix_get(_B,k,j,_nrowsB)
            end
            matrix_set(_C,elementC,i,j,_nrowsA)
        end
    end
    return _C
end


 
-- Square norm of vector _X with size _n
function vnorm(_X, _n)
    local _norm = 0.0
    for i = 0, _n-1 do
        _norm = _norm + math.pow( _X[i], 2)
    end
    return math.sqrt(_norm)
end

function qr_diag(_A, _nrows, _ncolumns)
    local _Q = {}
	for i= 1, _nrows*_ncolumns do
	_Q[i]=0;
	end
	
    local _R = {}
	for i= 1, _nrows*_ncolumns do
	_R[i]=0;
	end
	
    local _a = {}
	for i= 1, _nrows  do
	_a[i]=0;
	end	
    local _q = {}
	for i= 1, _nrows  do
	_q[i]=0;
	end
    local _r = 0.0
    local _aux = 0.0

    for i = 0, _nrows-1 do
        _a[i] = matrix_get(_A, i, 0, _nrows)
    end
    _r = vnorm(_a, _nrows)

    matrix_set(_R, _r, 0, 0, _ncolumns)
    for i = 0, _nrows-1 do
        matrix_set(_Q, _a[i]/_r, i, 0, _nrows)
    end

    if _ncolumns ~= 1 then
        for k = 1, _ncolumns-1 do
            for i = 0, _nrows-1 do
                _a[i] = matrix_get(_A, i, k, _nrows)
            end
            for j = 0, k-1 do
                _r = 0
                for i = 0, _nrows-1 do
                    _r = _r + matrix_get(_Q, i, j, _nrows)*_a[i]
                end
                matrix_set(_R, _r, j, k, _ncolumns)
                for i = 0, _nrows-1 do
                    _aux = _a[i] - _r*matrix_get(_Q, i, j, _nrows)
                    _a[i] = _aux
                end
            end
            _r = vnorm(_a, _nrows)
            matrix_set(_R, _r, k, k, _ncolumns)
            for i = 0, _nrows-1 do
                matrix_set(_Q, _a[i]/_r, i, k, _nrows)
            end
        end
    end

    return _Q, _R
 
end
 

function pinv(_A, _nrows, _ncolumns)
    -- Pseudoinverse of matrix _A calculated using QR decomposition
    -- Input: 
    -- _A: array: implied as a (_nrows x _ncolumns) matrix _A = [[column_0],[column_1],...,[column_(_ncolumns-1) 
    -- Output: 
    -- _Ainv: array implied as a (_ncolumns x _nrows) matrix _A = [[row_0],[row_1],...,[row_(_nrows-1) 
    -- ----
    -- First find the QR factorization of A: A = QR,
    -- where R is upper triangular matrix.
    -- Then _Ainv = R^-1*Q^T.
    -- ----

    _Q, _R = qr_diag(_A, _nrows, _ncolumns)
    _QT = transpose(_Q, _nrows, _ncolumns)

    -- Calculate Rinv:
    _Rinv = {}
	for i= 1, _ncolumns*_ncolumns do
	_Rinv[i]=0;
	end	
	
    _r = 0.0
    matrix_set(_Rinv, 1/matrix_get(_R, 0, 0, _ncolumns), 0, 0, _ncolumns)
    if _ncolumns ~= 1 then
        for j = 1, _ncolumns-1 do
            for i = 0, j-1 do
                _r = 0.0
                for k = i, j-1 do
                    _r = _r + matrix_get(_Rinv, i, k, _ncolumns) * matrix_get(_R, k, j, _ncolumns)
                end
                matrix_set(_Rinv, _r, i, j, _ncolumns)
            end
            for k = 0, j-1 do
                matrix_set(_Rinv, -matrix_get(_Rinv, k, j, _ncolumns) / matrix_get(_R, j, j, _ncolumns), k, j, _ncolumns)
            end
            matrix_set(_Rinv, 1/matrix_get(_R, j, j, _ncolumns), j, j, _ncolumns)
        end
    end

    _Ainv = multiply(_Rinv, _QT, _ncolumns, _ncolumns, _nrows)
    return _Ainv
end




function Differentiator(_src, _window, _degree, period) 	

    -- Polynomial differentiator
    -- input:
    -- _src: input series
    -- _window: integer: width of the moving lookback window
    -- _degree: integer: degree of fitting polynomial
    -- output:
    -- _diff: series: time derivative

    -- Vandermonde matrix with implied dimensions (window x degree+1)
    -- Linear form: J = [ [z]^0, [z]^1, ... [z]^degree], with z = [ (1-window)/2 to (window-1)/2 ] 
    _J = {}
	for i= 1, _window*(_degree+1) do
	_J[i]=0;
	end		
	
    for i = 0, _window-1 do
        
        for j = 0, _degree do
             matrix_set(_J,math.pow(i,j),i,j,_window)
			--            
        end
    end

    -- Vector of raw datapoints:
    _Y_raw = {}
    for j = 0, _window-1 do
        _Y_raw[j] = _src[period-_window+1+j]
    end

    -- Calculate polynomial coefficients which minimize the loss function
    -- Note: Lua does not have a built-in function for pseudo-inverse or matrix multiplication.
    -- You would need to use a library or implement these functions yourself.
    _C = pinv(_J, _window, _degree+1)
    _a_coef = multiply(_C, _Y_raw, _degree+1, _window, 1)

    -- For smoothing, approximate the last point (i.e. z=window-1) by a0
    _diff = 0.0
    for i = 1, _degree do
        _diff = _diff + i * _a_coef[i] * math.pow(_window-1, i-1)
    end

    return _diff
end

--+------------------------------------------------------------------------------------------------+
--|                                                                    We appreciate your support. | 
--+------------------------------------------------------------------------------------------------+
--|                                                               Paypal:  https://goo.gl/9Rj74e   |
--|                                                             Patreon :  http://tiny.cc/1ybwxz   |  
--|                                                      Buy Me a Coffee:  http://tiny.cc/bj7vxz   | 
--+------------------------------------------------------------------------------------------------+
--|  Cryptocurrency  |  Network                    |  Address                                      |
--+------------------------------------------------+-----------------------------------------------+
--|  USDT            |  ERC20 (ETH Ethereum)       |  0xe53aab6bc468a963a02d1319660ee60cf80fc8e7   | 
--|  USDT            |  TRC20 (Tron)               |  TTBXsfuPm2rk36AkdemY7muNXGjyziC86g           |
--|  USDT            |  BEP20 (BSC BNB Smart Chain)|  0xe53aab6bc468a963a02d1319660ee60cf80fc8e7   |
--|  USDT            |  Matic Polygon              |  0xe53aab6bc468a963a02d1319660ee60cf80fc8e7   |
--|  USDT            |  SOL Solana                 |  3nh5rpUKopcYLNU4zGCdUFAkM3iRQq8VVUmuzVG6VDf2 |
--|  USDT            |  ARBITRUM Arbitrum One      |  0xe53aab6bc468a963a02d1319660ee60cf80fc8e7   |
--+------------------------------------------------+-----------------------------------------------+