WIPNNetwork
// 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 Method for a generalized Neural Network.
library("WIPNNetwork")

// reference:
//      https://blog.primen.dk/understanding-backpropagation/
//      https://ml-cheatsheet.readthedocs.io/en/latest/backpropagation.html
//      https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/
//      https://www.anotsorandomwalk.com/backpropagation-example-with-numbers-step-by-step/
//      https://machinelearningmastery.com/implement-backpropagation-algorithm-scratch-python/
//      https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/
//      https://cedar.buffalo.edu/~srihari/CSE574/Chap5/Chap5.3-BackProp.pdf

import RicardoSantos/FunctionNNLayer/2 as layer
import RicardoSantos/MLLossFunctions/1 as loss
import RicardoSantos/MLActivationFunctions/1 as activation


// @function helper method to generate random weights of (A * B) size.
// @param previous_size int, number of nodes on the left.
// @param next_size int, number of nodes on the right.
// @returns float array.
generate_random_weights (int previous_size, int next_size) => //{
    _weights = matrix.new<float>(next_size, previous_size)
    for _i = 0 to next_size-1
        for _j = 0 to previous_size-1
            matrix.set(_weights, _i, _j, math.random() * 2.0 - 1.0)
    _weights
//}

mprint (M, format='0.000') => //{
    _rows = matrix.rows(M)
    _cols = matrix.columns(M)
    _s = ''
    for _i = 0 to _rows-1
        _s += '[ '
        for _j = 0 to _cols-1
            _s += ((_j == _rows-1) ? '-> ': '') + str.tostring(matrix.get(M, _i, _j), format) + ' '
        _s += ']\n'
    _s
//}

to_column_matrix (float[] a) => //{
    _m = matrix.new<float>()
    matrix.add_col(_m, 0, a)
    _m
//}

// @function helper method to propagate the signals forward.
// @param inputs float array, network inputs.
// @param weights float array, network weights.
// @param layer_sizes int array, layer sizes.
// @param weights_max_width int, the maximum row size of the weights matrix.
// @returns float array.
propagate_forward (float[] inputs, matrix<float> weights, int[] layer_sizes, string[] layer_functions, float[] layer_bias) => //{
	// TODO: might be good to remove output layer from loop? (inputs -> (loop hidden to hidden) -> outputs)
	//     : missing alpha, scale parameterfor activation functions that might require it..
    string _debug = '\n Forward Propagation:'

    int _size_i = array.size(inputs)
    int _size_ls = array.size(layer_sizes)
    int _size_lb = array.size(layer_bias)
    // use default parameters if empty array is provided:
    float[] _layer_bias = _size_lb > 0 ? layer_bias : array.new_float(_size_ls, 1.0)
    // * process (input -> 1st hidden layer) or (input -> output layer):
    //  TODO: bias, alpha, scale for each layer.
	float[] _layer_outputs = layer.function(
	         inputs                 = inputs, 
	         weights                = matrix.submatrix(weights, 0, array.get(layer_sizes, 0), 0, _size_i),
	         activation_function    = array.get(layer_functions, 0), 
	         bias = array.get(_layer_bias, 0), alpha = na, scale = na
	         )
    _debug += str.format('\n Current Layer: {0}, Current Layer size: {1}, Previous Layer size: {2}, Base index: {3}, f(): {4}, b: {5}', 0, array.get(layer_sizes, 0), _size_i, 0, array.get(layer_functions, 0), array.get(_layer_bias, 0))
    // if theres more than 1 (hidden layer + output)
	// * process hidden and outputs:
	if _size_ls > 1
        float[] _layer_output = _layer_outputs
    	for _layer = 1 to _size_ls - 1
    	    // find the base level of current layer weights:
    	    //      base = (sum of nodes up to previous layer) * (max weights in a node)
    	    //      (weights in layer) = base + (previous layer nodes) * (current layer nodes)
    	    int _previous_layer_size = array.get(layer_sizes, _layer - 1)
    	    int _this_layer_size = array.get(layer_sizes, _layer)
    	    int _w_base_idx = array.sum(array.slice(layer_sizes, 0, _layer))
    	    matrix<float> _layer_weights = matrix.submatrix(weights, _w_base_idx, _w_base_idx + _this_layer_size, 0, _previous_layer_size) //resize_weights_matrix(array.slice(weights, _w_base_idx, _w_base_idx + weights_max_width * _this_layer_size), weights_max_width, _previous_layer_size)
            // process hidden and output layer:
            //  TODO: bias, alpha, scale for each layer.
    	    _layer_output := layer.function(
    	             inputs              = _layer_output, 
    	             weights             = _layer_weights, 
    	             activation_function = array.get(layer_functions, _layer), 
    	             bias=array.get(_layer_bias, _layer), alpha=na, scale=na
    	             )
    	    array.concat(_layer_outputs, _layer_output)
    	    _debug += str.format('\n Current Layer: {0}, Current Layer size: {1}, Previous Layer size: {2}, Base index: {3}, f(): {4}, b: {5}', _layer, _this_layer_size, _previous_layer_size, _w_base_idx, array.get(layer_functions, _layer), array.get(_layer_bias, _layer))
    [_layer_outputs, _debug]
//}


// @function Generalized Neural Network Method.
// @param x TODO: add parameter x description here
// @returns TODO: add what function returns
export network (
     float[] inputs, float[] targets, matrix<float> weights, 
     int[] layer_sizes, float[] layer_biases, string[] layer_functions, 
     string loss_function='mse', float learning_rate=0.01
     ) => //{
    // TODO: Layer activation function extra parameters.
    
    // int _weights_max_width = array.max(layer_sizes) // highest number of nodes on each layer as base for weight rectangle matrix slicing.
    int _total_nodes = array.sum(layer_sizes)
    
    int _size_i = array.size(inputs)
    int _size_t = array.size(targets)
    int _size_ls = array.size(layer_sizes)
    int _size_lb = array.size(layer_biases)
    int _size_lf = array.size(layer_functions)

    int _Wrows = matrix.rows(weights)
    int _Wcols = matrix.columns(weights)
    
    // TODO: rework error messages.
    switch
        (_size_i < 1)   => runtime.error('FunctionNNetwork -> network(): "inputs" has wrong size.')
        (_size_t < 1)   => runtime.error('FunctionNNetwork -> network(): "targets" has wrong size.')
        // (_size_w < 1)   => runtime.error('FunctionNNetwork -> network(): "weights" has wrong size.')
        (_size_ls < 1)  => runtime.error('FunctionNNetwork -> network(): "layer_sizes" has wrong size.')
        (_size_lb < 1)  => runtime.error('FunctionNNetwork -> network(): "layer_biases" has wrong size.')
        (_size_lf < 1)  => runtime.error('FunctionNNetwork -> network(): "layer_functions" has wrong size.')
        (_size_t != array.get(layer_sizes, _size_ls - 1))  => runtime.error('FunctionNNetwork -> network(): "targets" size does not match output layer size in "layer_sizes"')
        (na(loss_function)) => runtime.error('FunctionNNetwork -> network(): "loss_function" is invalid.') //temp just for use sake
        (na(learning_rate)) => runtime.error('FunctionNNetwork -> network(): "learning_rate" is invalid.') //temp just for use sake
	//
	// propagate forward:
	[_layer_outputs, _debug] = propagate_forward(inputs, weights, layer_sizes, layer_functions, layer_biases)
// 	_layer_outputs = array.new<float>(10), _debug=''
    //
    // propagate backward:
    // error:
    // TODO: function selection:
    float _output_error = loss.mse(targets, array.slice(_layer_outputs, _total_nodes - (array.get(layer_sizes, _size_ls - 1)), _total_nodes))
    // reference for delta rule: https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/
    // node_delta_error = (nodeoutput - nodetarget) * derivative(nodeoutput)
    // update_weight = weight - learning_rate * node_error_delta
    //
    // reference: https://cedar.buffalo.edu/~srihari/CSE574/Chap5/Chap5.3-BackProp.pdf
    // Weight = weight + (learning rate * node output error * node weight input)
    //##########################################################################
    // Work in progress...
    //##########################################################################
    
    // * calculate delta from output layer:
    float _output_delta = _output_error * activation.derivative(array.get(layer_functions, _size_ls-1), _output_error)
    _debug += str.format('\n Output Error: {0}, Delta: {1}', _output_error, _output_delta)
    float[] _node_errors = array.new_float(0)//_total_nodes) // #######################################
    float[] _error_deltas = array.new_float(_total_nodes)
    int _last_delta_idx = -1
    
    // * process the error delta for the output layer:
    float[] _layer_deltas = array.new_float(_size_ls, 0.0)
    array.set(_layer_deltas, _size_ls-1, _output_delta)
    
    _debug += '\n\n Back Propagation:'
    for _layer_idx = (_size_ls - 1) to 0
        int _this_layer_size = array.get(layer_sizes, _layer_idx)
        int _previous_layer_size = _layer_idx == 0 ? _size_i : array.get(layer_sizes, _layer_idx - 1)
        int _number_of_nodes_up_to_this_layer = array.sum(array.slice(layer_sizes, 0, _layer_idx))
        // matrix<float> _layer_weights = matrix.submatrix(weights, _number_of_nodes_up_to_this_layer, (_number_of_nodes_up_to_this_layer + _this_layer_size), 0, _Wcols) //array.slice(weights, _weights_max_width * _number_of_nodes_up_to_this_layer, _weights_max_width * (_number_of_nodes_up_to_this_layer + _this_layer_size))
        string _layer_function = array.get(layer_functions, _layer_idx)
        //**********************************************************************>
        // dot product of layer weights and layer_delta
        float _current_layer_delta = array.get(_layer_deltas, _layer_idx)
        float _layer_error = 0.0
        for _j = 0 to _Wcols - 1
            _layer_error += matrix.get(weights, _number_of_nodes_up_to_this_layer, _j) * _current_layer_delta
        //**********************************************************************<
        _debug += str.format('\n Layer: {0}, Current size: {1}, Previous Size: {2}, e: {3}, weights:\n', _layer_idx, _this_layer_size, _previous_layer_size, _layer_error)
        // * calculate node output derivative:
        for _node_idx = 0 to (_this_layer_size - 1)
            float _node_output = array.get(_layer_outputs, _number_of_nodes_up_to_this_layer + _node_idx)
            float _node_derivative = activation.derivative(_layer_function, _node_output)
            if _layer_idx > 0
                array.set(_layer_deltas, _layer_idx - 1, _layer_error * _node_derivative)//*********************
            float _node_target = 1.0 //##########################
            float[] _node_inputs = _layer_idx == 0 ? inputs : array.slice(_layer_outputs, _number_of_nodes_up_to_this_layer - _previous_layer_size, _number_of_nodes_up_to_this_layer)
            //calculate layer error
            //update weight
            for _weight_idx = 0 to (_previous_layer_size - 1)
                float _weight = matrix.get(weights, _number_of_nodes_up_to_this_layer + _node_idx, _weight_idx)
                switch (_layer_idx)
                    (_size_ls - 1) => _node_target := array.get(targets, _node_idx)
                    => _node_target := array.get(_node_inputs, _weight_idx)
                float _error_delta = (_node_output - _node_target) * _node_derivative
                _weight += learning_rate * _error_delta
                _debug += str.format(' {0} ', _weight)
                matrix.set(weights, _number_of_nodes_up_to_this_layer + _node_idx, _weight_idx, _weight) //###########################
    
    [_layer_outputs, _layer_deltas, _debug]
//{ usage:
// 7, 5 = i3, h5 h3 o2
var float[] inputs = array.from(0.0, 1, 2)
var float[] expected_outputs = array.from(1.0, 0.0)
var matrix<float> weights = generate_random_weights(5, 10)
// weights = array.from(
//      0.000, 0.001, 0.002, float(na), float(na),
//      0.003, 0.004, 0.005, float(na), float(na), 
//      0.006, 0.007, 0.008, float(na), float(na),
//      0.009, 0.010, 0.011, float(na), float(na),
//      0.012, 0.013, 0.014, float(na), float(na),
//      0.100, 0.101, 0.102, 0.103, 0.104,
//      0.105, 0.106, 0.107, 0.108, 0.109,
//      0.110, 0.111, 0.112, 0.113, 0.114,
//      0.200, 0.201, 0.202, float(na), float(na),
//      0.203, 0.204, 0.205, float(na), float(na), 
//      0.206, 0.207, 0.208, float(na), float(na)
//      )
//                                              H           H           O
var int[]    layer_sizes     = array.from(        5,         3,         2)
var float[]  layer_biases    = array.from(      1.0,       1.0,       1.0)
var string[] layer_functions = array.from('sigmoid', 'sigmoid', 'sigmoid')
[o, d, _debug] = network(inputs, expected_outputs, weights, layer_sizes, layer_biases, layer_functions, 'mse')
//
if true//islastconfirmedhistory
    var label la = label.new(bar_index, 0.0, '')
    label.set_x(la, bar_index), label.set_text(la, str.format('output:{0}\ndelta:{1}\nweights:{2}, \n{3}', str.tostring(o), str.tostring(d), matrix.elements_count(weights), _debug))
    
    // print NN to table:
    var table T = table.new(position.bottom_right, 7, 3, bgcolor=#000000, frame_color=color.black, frame_width=2, border_width=2, border_color=color.gray)
    // run a simple NN forward cycle once:
    _wh1 = matrix.submatrix(weights, 0, 5, 0, 3)// next, previous, layers
    _wh2 = matrix.submatrix(weights, 5, 8, 0, 5)
    _wo = matrix.submatrix(weights, 8, 10, 0, 3)
    _h1 = array.slice(o, 0, 5)//function(i, _wh1) // 1st hidden layer with 4 nodes, from inputs
    _h2 = array.slice(o, 5, 8)//function(h1, _wh2) // 2nd level hidden layer with 5 nodes , from first hidden
    _o = array.slice(o, 8, 10)//function(h2, _wo) // output layer with 3 nodes 
    // label.new(bar_index, 0.0, str.tostring(matrix.mult(matrix.transpose(_wh1), i)))
    //
    table.cell(T, 0, 0, 'Input:', text_color=color.silver, text_halign=text.align_right)
    table.cell(T, 1, 0, 'Hidden 1:', text_color=color.silver, text_halign=text.align_right), 
    table.cell(T, 3, 0, 'Hidden 2:', text_color=color.silver, text_halign=text.align_right)
    table.cell(T, 5, 0, 'output:', text_color=color.silver, text_halign=text.align_right)
    table.cell(T, 0, 1, '', text_color=color.silver, text_halign=text.align_right)
    table.cell(T, 1, 1, 'weights:', text_color=color.silver, text_halign=text.align_right), 
    table.cell(T, 2, 1, 'node:', text_color=color.silver, text_halign=text.align_right), 
    table.cell(T, 3, 1, 'weights:', text_color=color.silver, text_halign=text.align_right)
    table.cell(T, 4, 1, 'node:', text_color=color.silver, text_halign=text.align_right), 
    table.cell(T, 5, 1, 'weights:', text_color=color.silver, text_halign=text.align_right)
    table.cell(T, 6, 1, 'node:', text_color=color.silver, text_halign=text.align_right), 
    table.merge_cells(T, 1, 0, 2, 0)
    table.merge_cells(T, 3, 0, 4, 0)
    table.merge_cells(T, 5, 0, 6, 0)
    table.cell(T, 0, 2, str.tostring(to_column_matrix(inputs), '0.000'), text_color=color.silver, text_halign=text.align_left)
    table.cell(T, 1, 2, str.tostring(_wh1, '0.000'), text_color=color.silver, text_halign=text.align_left)
    table.cell(T, 2, 2, str.tostring(to_column_matrix(_h1), '0.000'), text_color=color.silver, text_halign=text.align_left)
    table.cell(T, 3, 2, str.tostring(_wh2, '0.000'), text_color=color.silver, text_halign=text.align_left)
    table.cell(T, 4, 2, str.tostring(to_column_matrix(_h2), '0.000'), text_color=color.silver, text_halign=text.align_left)
    table.cell(T, 5, 2, str.tostring(_wo, '0.000'), text_color=color.silver, text_halign=text.align_left)
    table.cell(T, 6, 2, str.tostring(to_column_matrix(_o), '0.000'), text_color=color.silver, text_halign=text.align_left)
//{ remarks:
//}}}


