NeuralNetwork.cpp

/**
 * DESCRIPTION OF THE FILE
 *
 * @author Michal Kravčenko
 * @date 13.6.18 - 
 */
#include <armadillo>
#include <iostream>
#include <4neuro.h>

#include "message.h"
#include "NeuralNetwork.h"
#include "NeuralNetworkSerialization.h"
#include "exceptions.h"


namespace lib4neuro {
    NeuralNetwork::NeuralNetwork() {
        this->gen = boost::random::mt19937(std::time(0));
        this->neurons = new ::std::vector<Neuron *>(0);
        this->neuron_biases = new ::std::vector<double>(0);
        this->neuron_potentials = new ::std::vector<double>(0);
        this->neuron_bias_indices = new ::std::vector<int>(0);

        this->connection_weights = new ::std::vector<double>(0);
        this->connection_list = new ::std::vector<ConnectionFunctionGeneral *>(0);
        this->inward_adjacency = new ::std::vector<std::vector<std::pair<size_t, size_t>> *>(0);
        this->outward_adjacency = new ::std::vector<std::vector<std::pair<size_t, size_t>> *>(0);

        this->neuron_layers_feedforward = new ::std::vector<std::vector<size_t> *>(0);
        this->neuron_layers_feedbackward = new ::std::vector<std::vector<size_t> *>(0);

        this->input_neuron_indices = new ::std::vector<size_t>(0);
        this->output_neuron_indices = new ::std::vector<size_t>(0);

        this->delete_weights = true;
        this->delete_biases = true;
        this->layers_analyzed = false;
    }

    NeuralNetwork::NeuralNetwork(std::string filepath) {
        ::std::ifstream ifs(filepath);
        if(ifs.is_open()) {
            try {
                boost::archive::text_iarchive ia(ifs);
                ia >> *this;
            }catch(boost::archive::archive_exception& e) {
                THROW_RUNTIME_ERROR("Serialized archive error: '" + e.what() + "'! Please, check if your file is really "
                                                                               "the serialized DataSet.");
            }
            ifs.close();
        } else {
            THROW_RUNTIME_ERROR("File '" + filepath + "' couldn't be open!");
        }

        this->gen = boost::random::mt19937(std::time(0));
    }

    NeuralNetwork::~NeuralNetwork() {

        if (this->neurons) {
            for (auto n: *(this->neurons)) {
                delete n;
                n = nullptr;
            }
            delete this->neurons;
            this->neurons = nullptr;
        }

        if (this->neuron_potentials) {
            delete this->neuron_potentials;
            this->neuron_potentials = nullptr;
        }

        if (this->neuron_bias_indices) {
            delete this->neuron_bias_indices;
            this->neuron_bias_indices = nullptr;
        }

        if (this->output_neuron_indices) {
            delete this->output_neuron_indices;
            this->output_neuron_indices = nullptr;
        }

        if (this->input_neuron_indices) {
            delete this->input_neuron_indices;
            this->input_neuron_indices = nullptr;
        }

        if (this->connection_weights && this->delete_weights) {
            delete this->connection_weights;
            this->connection_weights = nullptr;
        }

        if (this->neuron_biases && this->delete_biases) {
            delete this->neuron_biases;
            this->neuron_biases = nullptr;
        }

        if (this->connection_list) {

            if (this->delete_weights) {
                for (auto c: *this->connection_list) {
                    delete c;
                    c = nullptr;
                }
            }
        }
        delete this->connection_list;
        this->connection_list = nullptr;

        if (this->inward_adjacency) {
            for (auto e: *this->inward_adjacency) {
                if (e) {
                    delete e;
                    e = nullptr;
                }
            }
            delete this->inward_adjacency;
            this->inward_adjacency = nullptr;
        }

        if (this->outward_adjacency) {
            for (
                auto e: *this->outward_adjacency) {
                if (e) {
                    delete e;
                    e = nullptr;
                }
            }
            delete this->
                    outward_adjacency;
            this->
                    outward_adjacency = nullptr;
        }

        if (this->neuron_layers_feedforward) {
            for (
                auto e: *this->neuron_layers_feedforward) {
                delete e;
                e = nullptr;
            }
            delete this->neuron_layers_feedforward;
            this->neuron_layers_feedforward = nullptr;
        }

        if (this->neuron_layers_feedbackward) {
            for (
                auto e: *this->neuron_layers_feedbackward) {
                delete e;
                e = nullptr;
            }
            delete this->neuron_layers_feedbackward;
            this->neuron_layers_feedbackward = nullptr;
        }

    }

    NeuralNetwork *NeuralNetwork::get_subnet(::std::vector<size_t> &input_neuron_indices,
                                             ::std::vector<size_t> &output_neuron_indices) {

        THROW_NOT_IMPLEMENTED_ERROR();

        NeuralNetwork *output_net = nullptr;
// TODO rework due to the changed structure of the class
//    Neuron * active_neuron, * target_neuron;
//
//    size_t n = this->neurons->size();
//    bool *part_of_subnetwork = new bool[n];
//    ::std::fill(part_of_subnetwork, part_of_subnetwork + n, false);
//
//    bool *is_reachable_from_source = new bool[n];
//    bool *is_reachable_from_destination = new bool[n];
//    ::std::fill(is_reachable_from_source, is_reachable_from_source + n, false);
//    ::std::fill(is_reachable_from_destination, is_reachable_from_destination + n, false);
//
//    bool *visited_neurons = new bool[n];
//    ::std::fill(visited_neurons, visited_neurons + n, false);
//
//    size_t active_set_size[2];
//    active_set_size[0] = 0;
//    active_set_size[1] = 0;
//    size_t * active_neuron_set = new size_t[2 * n];
//    size_t idx1 = 0, idx2 = 1;
//
//    /* MAPPING BETWEEN NEURONS AND THEIR INDICES */
//    size_t idx = 0, idx_target;
//    for(Neuron *v: *this->neurons){
//        v->set_idx( idx );
//        idx++;
//    }
//
//    /* INITIAL STATE FOR THE FORWARD PASS */
//    for(size_t i: input_neuron_indices ){
//
//        if( i < 0 || i >= n){
//            //invalid index
//            continue;
//        }
//        active_neuron_set[idx1 * n + active_set_size[idx1]] = i;
//        active_set_size[idx1]++;
//
//        visited_neurons[i] = true;
//    }
//
//    /* FORWARD PASS */
//    while(active_set_size[idx1] > 0){
//
//        //we iterate through the active neurons and propagate the signal
//        for(int i = 0; i < active_set_size[idx1]; ++i){
//            idx = active_neuron_set[i];
//
//            is_reachable_from_source[ idx ] = true;
//
//            active_neuron = this->neurons->at( idx );
//
//            for(Connection* connection: *(active_neuron->get_connections_out())){
//
//                target_neuron = connection->get_neuron_out( );
//                idx_target = target_neuron->get_idx( );
//
//                if( visited_neurons[idx_target] ){
//                    //this neuron was already analyzed
//                    continue;
//                }
//
//                visited_neurons[idx_target] = true;
//                active_neuron_set[active_set_size[idx2] + n * idx2] = idx_target;
//                active_set_size[idx2]++;
//            }
//        }
//        idx1 = idx2;
//        idx2 = (idx1 + 1) % 2;
//        active_set_size[idx2] = 0;
//    }
//
//
//    /* INITIAL STATE FOR THE FORWARD PASS */
//    active_set_size[0] = active_set_size[1] = 0;
//    ::std::fill(visited_neurons, visited_neurons + n, false);
//
//    for(size_t i: output_neuron_indices ){
//
//        if( i < 0 || i >= n){
//            //invalid index
//            continue;
//        }
//        active_neuron_set[idx1 * n + active_set_size[idx1]] = i;
//        active_set_size[idx1]++;
//
//        visited_neurons[i] = true;
//    }
//
//    /* BACKWARD PASS */
//    size_t n_new_neurons = 0;
//    while(active_set_size[idx1] > 0){
//
//        //we iterate through the active neurons and propagate the signal
//        for(int i = 0; i < active_set_size[idx1]; ++i){
//            idx = active_neuron_set[i];
//
//            is_reachable_from_destination[ idx ] = true;
//
//            active_neuron = this->neurons->at( idx );
//
//            if(is_reachable_from_source[ idx ]){
//                n_new_neurons++;
//            }
//
//            for(Connection* connection: *(active_neuron->get_connections_in())){
//
//                target_neuron = connection->get_neuron_in( );
//                idx_target = target_neuron->get_idx( );
//
//                if( visited_neurons[idx_target] ){
//                    //this neuron was already analyzed
//                    continue;
//                }
//
//                visited_neurons[idx_target] = true;
//                active_neuron_set[active_set_size[idx2] + n * idx2] = idx_target;
//                active_set_size[idx2]++;
//            }
//        }
//        idx1 = idx2;
//        idx2 = (idx1 + 1) % 2;
//        active_set_size[idx2] = 0;
//    }
//
//    /* FOR CONSISTENCY REASONS */
//    for(size_t in: input_neuron_indices){
//        if( !is_reachable_from_destination[in] ){
//            n_new_neurons++;
//        }
//        is_reachable_from_destination[in] = true;
//    }
//    /* FOR CONSISTENCY REASONS */
//    for(size_t in: output_neuron_indices){
//        if( !is_reachable_from_source[in] ){
//            n_new_neurons++;
//        }
//        is_reachable_from_source[in] = true;
//    }
//
//    /* WE FORM THE SET OF NEURONS IN THE OUTPUT NETWORK  */
//    if(n_new_neurons > 0){
////        printf("Number of new neurons: %d\n", n_new_neurons);
//        output_net = new NeuralNetwork();
//        output_net->set_weight_array( this->connection_weights );
//
//        ::std::vector<size_t > local_inputs(0), local_outputs(0);
//        local_inputs.reserve(input_neuron_indices.size());
//        local_outputs.reserve(output_neuron_indices.size());
//
//        ::std::vector<Neuron*> local_n_arr(0);
//        local_n_arr.reserve( n_new_neurons );
//
//        ::std::vector<Neuron*> local_local_n_arr(0);
//        local_local_n_arr.reserve( n_new_neurons );
//
//        int * neuron_local_mapping = new int[ n ];
//        ::std::fill(neuron_local_mapping, neuron_local_mapping + n, -1);
//        idx = 0;
//        for(size_t i = 0; i < n; ++i){
//            if(is_reachable_from_source[i] && is_reachable_from_destination[i]){
//                neuron_local_mapping[i] = (int)idx;
//                idx++;
//
//                Neuron *new_neuron = this->neurons->at(i)->get_copy( );
//
//                output_net->add_neuron( new_neuron );
//                local_local_n_arr.push_back( new_neuron );
//                local_n_arr.push_back( this->neurons->at(i) );
//            }
//        }
//        for(size_t in: input_neuron_indices){
//            local_inputs.push_back(neuron_local_mapping[in]);
//        }
//        for(size_t in: output_neuron_indices){
//            local_outputs.push_back(neuron_local_mapping[in]);
//        }
//
////        printf("%d\n", local_n_arr.size());
////        printf("inputs: %d, outputs: %d\n", local_inputs.size(), local_outputs.size());
//        int local_idx_1, local_idx_2;
//        for(Neuron* source_neuron: local_n_arr){
//            //we also add the relevant edges
//            local_idx_1 = neuron_local_mapping[source_neuron->get_idx()];
//
//            for(Connection* connection: *(source_neuron->get_connections_out( ))){
//                target_neuron = connection->get_neuron_out();
//
//                local_idx_2 = neuron_local_mapping[target_neuron->get_idx()];
//                if(local_idx_2 >= 0){
//                    //this edge is part of the subnetwork
//                    Connection* new_connection = connection->get_copy( local_local_n_arr[local_idx_1], local_local_n_arr[local_idx_2] );
//
//                    local_local_n_arr[local_idx_1]->add_connection_out(new_connection);
//                    local_local_n_arr[local_idx_2]->add_connection_in(new_connection);
//
////                    printf("adding a connection between neurons %d, %d\n", local_idx_1, local_idx_2);
//                }
//
//            }
//
//        }
//        output_net->specify_input_neurons(local_inputs);
//        output_net->specify_output_neurons(local_outputs);
//
//
//        delete [] neuron_local_mapping;
//    }
//
//    delete [] is_reachable_from_source;
//    delete [] is_reachable_from_destination;
//    delete [] part_of_subnetwork;
//    delete [] visited_neurons;
//    delete [] active_neuron_set;
//
//
        return output_net;
    }

    size_t NeuralNetwork::add_neuron(Neuron *n, BIAS_TYPE bt, size_t bias_idx) {

        if (bt == BIAS_TYPE::NO_BIAS) {
            this->neuron_bias_indices->push_back(-1);
        } else if (bt == BIAS_TYPE::NEXT_BIAS) {
            this->neuron_bias_indices->push_back((int) this->neuron_biases->size());
            this->neuron_biases->resize(this->neuron_biases->size() + 1);
        } else if (bt == BIAS_TYPE::EXISTING_BIAS) {
            if (bias_idx >= this->neuron_biases->size()) {
                ::std::cerr << "The supplied bias index is too large!\n" << ::std::endl;
            }
            this->neuron_bias_indices->push_back((int) bias_idx);
        }

        this->outward_adjacency->push_back(new ::std::vector<std::pair<size_t, size_t>>(0));
        this->inward_adjacency->push_back(new ::std::vector<std::pair<size_t, size_t>>(0));

        this->neurons->push_back(n);

        this->layers_analyzed = false;
        return this->neurons->size() - 1;
    }

    size_t
    NeuralNetwork::add_connection_simple(size_t n1_idx, size_t n2_idx, SIMPLE_CONNECTION_TYPE sct,
                                         size_t weight_idx) {

        ConnectionFunctionIdentity *con_weight_u1u2;
        if (sct == SIMPLE_CONNECTION_TYPE::UNITARY_WEIGHT) {
            con_weight_u1u2 = new ConnectionFunctionIdentity();
        } else {
            if (sct == SIMPLE_CONNECTION_TYPE::NEXT_WEIGHT) {
                weight_idx = this->connection_weights->size();
                this->connection_weights->resize(weight_idx + 1);
            } else if (sct == SIMPLE_CONNECTION_TYPE::EXISTING_WEIGHT) {
                if (weight_idx >= this->connection_weights->size()) {
                    ::std::cerr << "The supplied connection weight index is too large!\n" << ::std::endl;
                }
            }

            con_weight_u1u2 = new ConnectionFunctionIdentity(weight_idx);
        }

        size_t conn_idx = this->add_new_connection_to_list(con_weight_u1u2);

        this->add_outward_connection(n1_idx, n2_idx, conn_idx);
        this->add_inward_connection(n2_idx, n1_idx, conn_idx);

        this->layers_analyzed = false;

        return this->connection_list->size() - 1;
    }

    void NeuralNetwork::add_existing_connection(size_t n1_idx, size_t n2_idx, size_t connection_idx,
                                                NeuralNetwork &parent_network) {

        size_t conn_idx = this->add_new_connection_to_list(parent_network.connection_list->at(connection_idx));

        this->add_outward_connection(n1_idx, n2_idx, conn_idx);
        this->add_inward_connection(n2_idx, n1_idx, conn_idx);

        this->layers_analyzed = false;
    }

    void NeuralNetwork::copy_parameter_space(std::vector<double> *parameters) {
        if (parameters != nullptr) {
            for (unsigned int i = 0; i < this->connection_weights->size(); ++i) {
                (*this->connection_weights).at(i) = (*parameters).at(i);
            }

            for (unsigned int i = 0; i < this->neuron_biases->size(); ++i) {
                (*this->neuron_biases).at(i) = (*parameters).at(i + this->connection_weights->size());
            }
        }
    }

    void NeuralNetwork::set_parameter_space_pointers(NeuralNetwork &parent_network) {

        if (this->connection_weights) {
            delete connection_weights;
        }

        if (this->neuron_biases) {
            delete this->neuron_biases;
        }

        this->connection_weights = parent_network.connection_weights;
        this->neuron_biases = parent_network.neuron_biases;

        this->delete_biases = false;
        this->delete_weights = false;
    }

    void NeuralNetwork::add_to_jacobian(arma::Mat<double> &J, size_t row_idx,
                                        std::pair<std::vector<double>, std::vector<double>> &data) {

        std::vector<double> fv(this->get_n_outputs());
        std::vector<double> gradient(this->get_n_biases() + this->get_n_weights());

        this->eval_single( data.first, fv );

        std::vector<double> error_partial(this->get_n_outputs());
        std::fill(error_partial.begin(), error_partial.end(), 0.0);

        for( size_t i = 0; i < this->get_n_outputs(); ++i){
            error_partial[i] = 1;
            std::fill(gradient.begin(), gradient.end(), 0.0);
            this->add_to_gradient_single(data.first, error_partial, 1.0, gradient);
            error_partial[i] = 0;

            for( size_t ci = 0; ci < gradient.size(); ++ci){
                J.at(row_idx + i, ci) += gradient[ ci ];
            }
        }
    }

    void NeuralNetwork::eval_single(::std::vector<double>& input,
                                    ::std::vector<double>& output,
                                    ::std::vector<double>* custom_weights_and_biases) {

        if ((this->input_neuron_indices->size() * this->output_neuron_indices->size()) <= 0) {
            THROW_INVALID_ARGUMENT_ERROR("Input and output neurons have not been specified!");
        }

        if (this->input_neuron_indices->size() != input.size()) {
            THROW_INVALID_ARGUMENT_ERROR("Data input size != Network input size");
        }

        if (this->output_neuron_indices->size() != output.size()) {
            THROW_INVALID_ARGUMENT_ERROR("Data output size != Network output size");
        }

        double potential, bias;
        int bias_idx;

        this->copy_parameter_space(custom_weights_and_biases);

        this->analyze_layer_structure();

        /* reset of the output and the neuron potentials */
        ::std::fill(output.begin(), output.end(), 0.0);
        ::std::fill(this->neuron_potentials->begin(), this->neuron_potentials->end(), 0.0);

        /* set the potentials of the input neurons */
        for (size_t i = 0; i < this->input_neuron_indices->size(); ++i) {
            this->neuron_potentials->at(this->input_neuron_indices->at(i)) = input[i];
        }

        /* we iterate through all the feed-forward layers and transfer the signals */
        for (auto layer: *this->neuron_layers_feedforward) {
            /* we iterate through all neurons in this layer and propagate the signal to the neighboring neurons */

            for (auto si: *layer) {
                bias = 0.0;
                bias_idx = this->neuron_bias_indices->at(si);
                if (bias_idx >= 0) {
                    bias = this->neuron_biases->at(bias_idx);
                }
                potential = this->neurons->at(si)->activate(this->neuron_potentials->at(si), bias);

                for (auto c: *this->outward_adjacency->at(si)) {
                    size_t ti = c.first;
                    size_t ci = c.second;

                    this->neuron_potentials->at(ti) +=
                            this->connection_list->at(ci)->eval(*this->connection_weights) * potential;
                }
            }
        }

        unsigned int i = 0;
        for (auto oi: *this->output_neuron_indices) {
            bias = 0.0;
            bias_idx = this->neuron_bias_indices->at(oi);
            if (bias_idx >= 0) {
                bias = this->neuron_biases->at(bias_idx);
            }
            output[i] = this->neurons->at(oi)->activate(this->neuron_potentials->at(oi), bias);
            ++i;
        }
    }

    void NeuralNetwork::eval_single_debug(::std::vector<double> &input, ::std::vector<double> &output,
                                    ::std::vector<double> *custom_weights_and_biases) {
        if ((this->input_neuron_indices->size() * this->output_neuron_indices->size()) <= 0) {
            THROW_INVALID_ARGUMENT_ERROR("Input and output neurons have not been specified!");
        }

        if (this->input_neuron_indices->size() != input.size()) {
            THROW_INVALID_ARGUMENT_ERROR("Data input size != Network input size");
        }

        if (this->output_neuron_indices->size() != output.size()) {
            THROW_INVALID_ARGUMENT_ERROR("Data output size != Network output size");
        }

        double potential, bias;
        int bias_idx;

        this->copy_parameter_space(custom_weights_and_biases);

        this->analyze_layer_structure();

        /* reset of the output and the neuron potentials */
        ::std::fill(output.begin(), output.end(), 0.0);
        ::std::fill(this->neuron_potentials->begin(), this->neuron_potentials->end(), 0.0);

        /* set the potentials of the input neurons */
        for (size_t i = 0; i < this->input_neuron_indices->size(); ++i) {
            this->neuron_potentials->at(this->input_neuron_indices->at(i)) = input[i];
            std::cout << this->neuron_potentials->at(this->input_neuron_indices->at(i)) << ", ";
        }
        std::cout << std::endl;



        /* we iterate through all the feed-forward layers and transfer the signals */
        for (auto layer: *this->neuron_layers_feedforward) {
            /* we iterate through all neurons in this layer and propagate the signal to the neighboring neurons */

            for (auto si: *layer) {
                bias = 0.0;
                bias_idx = this->neuron_bias_indices->at(si);
                if (bias_idx >= 0) {
                    bias = this->neuron_biases->at(bias_idx);
                }
                potential = this->neurons->at(si)->activate(this->neuron_potentials->at(si), bias);
                std::cout << "  applying bias: " << bias << " to neuron potential: " << this->neuron_potentials->at(si) << " -> " << potential << std::endl;

                for (auto c: *this->outward_adjacency->at(si)) {
                    size_t ti = c.first;
                    size_t ci = c.second;

                    this->neuron_potentials->at(ti) +=
                            this->connection_list->at(ci)->eval(*this->connection_weights) * potential;

                    std::cout << "  adding input to neuron " << ti << " += " << this->connection_list->at(ci)->eval(*this->connection_weights) << "*" << potential << std::endl;
                }
            }
        }

        unsigned int i = 0;
        for (auto oi: *this->output_neuron_indices) {
            bias = 0.0;
            bias_idx = this->neuron_bias_indices->at(oi);
            if (bias_idx >= 0) {
                bias = this->neuron_biases->at(bias_idx);
            }
            output[i] = this->neurons->at(oi)->activate(this->neuron_potentials->at(oi), bias);
            std::cout << "setting the output[" << i << "] = " << output[i] << "(bias = " << bias << ")" << std::endl;
            ++i;
        }
    }

    void NeuralNetwork::add_to_gradient_single(std::vector<double> &input, ::std::vector<double> &error_derivative,
                                               double error_scaling, ::std::vector<double> &gradient) {


        ::std::vector<double> scaling_backprog(this->get_n_neurons());
        ::std::fill(scaling_backprog.begin(), scaling_backprog.end(), 0.0);

        size_t bias_shift = this->get_n_weights();
        size_t neuron_idx;
        int bias_idx;
        double neuron_potential, neuron_potential_t, neuron_bias, connection_weight;

        NeuronDifferentiable *active_neuron;

        /* initial error propagation */
        ::std::vector<size_t> *current_layer = this->neuron_layers_feedforward->at(
                this->neuron_layers_feedforward->size() - 1);
        //TODO might not work in the future as the output neurons could be permuted
        for (size_t i = 0; i < current_layer->size(); ++i) {
            neuron_idx = current_layer->at(i);
            scaling_backprog[neuron_idx] = error_derivative[i] * error_scaling;
        }

        /* we iterate through all the layers in reverse order and calculate partial derivatives scaled correspondingly */
        for (size_t j = this->neuron_layers_feedforward->size(); j > 0; --j) {

            current_layer = this->neuron_layers_feedforward->at(j - 1);

            for (size_t i = 0; i < current_layer->size(); ++i) {

                neuron_idx = current_layer->at(i);
                active_neuron = dynamic_cast<NeuronDifferentiable *> (this->neurons->at(neuron_idx));

                if (active_neuron) {
                    bias_idx = this->neuron_bias_indices->at(neuron_idx);
                    neuron_potential = this->neuron_potentials->at(neuron_idx);

                    if (bias_idx >= 0) {
                        neuron_bias = this->neuron_biases->at(bias_idx);
                        gradient[bias_shift + bias_idx] += scaling_backprog[neuron_idx] *
                                                           active_neuron->activation_function_eval_derivative_bias(
                                                                   neuron_potential, neuron_bias);
                        scaling_backprog[neuron_idx] *= active_neuron->activation_function_eval_derivative(
                                neuron_potential,
                                neuron_bias);
                    }

                    /* connections to lower level neurons */
                    for (auto c: *this->inward_adjacency->at(neuron_idx)) {
                        size_t ti = c.first;
                        size_t ci = c.second;

                        neuron_potential_t = this->neurons->at(ti)->get_last_activation_value( );
                        connection_weight = this->connection_list->at(ci)->eval(*this->connection_weights);

                        this->connection_list->at(ci)->eval_partial_derivative(*this->get_parameter_ptr_weights(),
                                                                               gradient,
                                                                               neuron_potential_t *
                                                                               scaling_backprog[neuron_idx]);

                        scaling_backprog[ti] += scaling_backprog[neuron_idx] * connection_weight;
                    }
                } else {
                    THROW_INVALID_ARGUMENT_ERROR(
                            "Neuron used in backpropagation does not contain differentiable activation function!\n");
                }
            }
        }
    }

    void NeuralNetwork::add_to_gradient_single_debug(std::vector<double> &input, ::std::vector<double> &error_derivative,
                                               double error_scaling, ::std::vector<double> &gradient) {


        ::std::vector<double> scaling_backprog(this->get_n_neurons());
        ::std::fill(scaling_backprog.begin(), scaling_backprog.end(), 0.0);

        size_t bias_shift = this->get_n_weights();
        size_t neuron_idx;
        int bias_idx;
        double neuron_potential, neuron_activation_t, neuron_bias, connection_weight;

        NeuronDifferentiable *active_neuron;

        /* initial error propagation */
        ::std::vector<size_t> *current_layer = this->neuron_layers_feedforward->at(
                this->neuron_layers_feedforward->size() - 1);
        //TODO might not work in the future as the output neurons could be permuted
        std::cout << "Error scaling on the output layer: ";
        for (size_t i = 0; i < current_layer->size(); ++i) {
            neuron_idx = current_layer->at(i);
            scaling_backprog[neuron_idx] = error_derivative[i] * error_scaling;

            std::cout << scaling_backprog[neuron_idx] << " [neuron " << neuron_idx << "], ";
        }
        std::cout << std::endl;

        /* we iterate through all the layers in reverse order and calculate partial derivatives scaled correspondingly */
        for (size_t j = this->neuron_layers_feedforward->size(); j > 0; --j) {

            current_layer = this->neuron_layers_feedforward->at(j - 1);

            for (size_t i = 0; i < current_layer->size(); ++i) {

                neuron_idx = current_layer->at(i);
                active_neuron = dynamic_cast<NeuronDifferentiable *> (this->neurons->at(neuron_idx));

                if (active_neuron) {
                    std::cout << "  [backpropagation] active neuron: " << neuron_idx << std::endl;

                    bias_idx = this->neuron_bias_indices->at(neuron_idx);
                    neuron_potential = this->neuron_potentials->at(neuron_idx);

                    if (bias_idx >= 0) {
                        neuron_bias = this->neuron_biases->at(bias_idx);
                        gradient[bias_shift + bias_idx] += scaling_backprog[neuron_idx] *
                                                           active_neuron->activation_function_eval_derivative_bias(
                                                                   neuron_potential, neuron_bias);
                        scaling_backprog[neuron_idx] *= active_neuron->activation_function_eval_derivative(
                                neuron_potential,
                                neuron_bias);
                    }

                    std::cout << "      [backpropagation] scaling coefficient: " << scaling_backprog[neuron_idx] << std::endl;

                    /* connections to lower level neurons */
                    for (auto c: *this->inward_adjacency->at(neuron_idx)) {
                        size_t ti = c.first;
                        size_t ci = c.second;

                        neuron_activation_t = this->neurons->at(ti)->get_last_activation_value( );
                        connection_weight = this->connection_list->at(ci)->eval(*this->connection_weights);

                        std::cout << "      [backpropagation] value ("<<ti<< "): " << neuron_activation_t << ", scaling: " << scaling_backprog[neuron_idx] << std::endl;

                        this->connection_list->at(ci)->eval_partial_derivative(*this->get_parameter_ptr_weights(),
                                                                               gradient,
                                                                               neuron_activation_t *
                                                                               scaling_backprog[neuron_idx]);

                        scaling_backprog[ti] += scaling_backprog[neuron_idx] * connection_weight;
                    }
                } else {
                    THROW_INVALID_ARGUMENT_ERROR(
                            "Neuron used in backpropagation does not contain differentiable activation function!\n");
                }
            }
        }
    }

    void NeuralNetwork::randomize_weights() {

        // Init weight guess ("optimal" for logistic activation functions)
        double r = 1.0 / (this->neuron_biases->size() + this->connection_weights->size());

        boost::random::uniform_real_distribution<> dist(-r, r);

        for (size_t i = 0; i < this->connection_weights->size(); i++) {
            this->connection_weights->at(i) = dist(gen);
//            std::cout << "weight[" << i <<"]" << this->connection_weights->at(i) << std::endl;
        }
    }

    void NeuralNetwork::randomize_biases() {



        double r = 1.0 / (this->neuron_biases->size() + this->connection_weights->size());
        // Init weight guess ("optimal" for logistic activation functions)
        boost::random::uniform_real_distribution<> dist(-r, r);
        for (size_t i = 0; i < this->neuron_biases->size(); i++) {
            this->neuron_biases->at(i) = dist(gen);
//            std::cout << "bias[" << i <<"]" << this->neuron_biases->at(i) << std::endl;
        }
    }

    void NeuralNetwork::randomize_parameters() {
        this->randomize_biases();
        this->randomize_weights();
    }

    void NeuralNetwork::scale_biases(double alpha) {
        for(size_t i = 0; i < this->get_n_biases(); ++i){
            this->neuron_biases->at( i ) *= alpha;
        }
    }

    void NeuralNetwork::scale_weights(double alpha) {
        for(size_t i = 0; i < this->get_n_weights(); ++i){
            this->connection_weights->at( i ) *= alpha;
        }
    }

    void NeuralNetwork::scale_parameters(double alpha) {
        this->scale_biases( alpha );
        this->scale_weights( alpha );
    }

    size_t NeuralNetwork::get_n_inputs() {
        return this->input_neuron_indices->size();
    }

    size_t NeuralNetwork::get_n_outputs() {
        return this->output_neuron_indices->size();
    }

    size_t NeuralNetwork::get_n_weights() {
        return this->connection_weights->size();
    }

    size_t NeuralNetwork::get_n_biases() {
        return this->neuron_biases->size();
    }

    int NeuralNetwork::get_neuron_bias_index(size_t neuron_idx) {
        return this->neuron_bias_indices->at(neuron_idx);
    }

    size_t NeuralNetwork::get_n_neurons() {
        return this->neurons->size();
    }

    void NeuralNetwork::specify_input_neurons(std::vector<size_t> &input_neurons_indices) {
        if (!this->input_neuron_indices) {
            this->input_neuron_indices = new ::std::vector<size_t>(input_neurons_indices);
        } else {
            delete this->input_neuron_indices;
            this->input_neuron_indices = new ::std::vector<size_t>(input_neurons_indices);
        }
    }

    void NeuralNetwork::specify_output_neurons(std::vector<size_t> &output_neurons_indices) {
        if (!this->output_neuron_indices) {
            this->output_neuron_indices = new ::std::vector<size_t>(output_neurons_indices);
        } else {
            delete this->output_neuron_indices;
            this->output_neuron_indices = new ::std::vector<size_t>(output_neurons_indices);
        }
    }

    void NeuralNetwork::write_weights() {
        std::cout << "Connection weights: ";
        if (this->connection_weights) {
            for (size_t i = 0; i < this->connection_weights->size() - 1; ++i) {
                std::cout << this->connection_weights->at(i) << ", ";
            }
            std::cout << this->connection_weights->at(this->connection_weights->size() - 1) << std::endl;
        }
    }

    void NeuralNetwork::write_weights(std::string file_path) {
        std::ofstream ofs(file_path);

        if(!ofs.is_open()) {
            THROW_RUNTIME_ERROR("File " + file_path + " can not be opened!");
        }

        ofs << "Connection weights: ";

        if (this->connection_weights) {
            for (size_t i = 0; i < this->connection_weights->size() - 1; ++i) {
                ofs << this->connection_weights->at(i) << ", ";
            }
            ofs << this->connection_weights->at(this->connection_weights->size() - 1) << std::endl;
        }
    }

    void NeuralNetwork::write_weights(std::ofstream* file_path) {
        *file_path << "Connection weights: ";

        if (this->connection_weights) {
            for (size_t i = 0; i < this->connection_weights->size() - 1; ++i) {
                *file_path << this->connection_weights->at(i) << ", ";
            }
            *file_path << this->connection_weights->at(this->connection_weights->size() - 1) << std::endl;
        }
    }

    void NeuralNetwork::write_biases() {
        std::cout << "Network biases: ";

        if(this->neuron_biases) {
            for(unsigned int i = 0; i < this->neuron_biases->size() - 1; i++) {
                std::cout << this->neuron_biases->at(i) << ", ";
            }
            std::cout << this->neuron_biases->at(this->neuron_biases->size() - 1) << std::endl;
        }
    }

    void NeuralNetwork::write_biases(std::string file_path) {
        std::ofstream ofs(file_path);

        if(!ofs.is_open()) {
            THROW_RUNTIME_ERROR("File " + file_path + " can not be opened!");
        }

        ofs << "Network biases: ";

        if(this->neuron_biases) {
            for(unsigned int i = 0; i < this->neuron_biases->size() - 1; i++) {
                ofs << this->neuron_biases->at(i) << ", ";
            }
            ofs << this->neuron_biases->at(this->neuron_biases->size() - 1) << std::endl;
        }
    }

    void NeuralNetwork::write_biases(std::ofstream* file_path) {
        *file_path << "Network biases: ";

        if(this->neuron_biases) {
            for(unsigned int i = 0; i < this->neuron_biases->size() - 1; i++) {
                *file_path << this->neuron_biases->at(i) << ", ";
            }
            *file_path << this->neuron_biases->at(this->neuron_biases->size() - 1) << std::endl;
        }
    }

    void NeuralNetwork::write_stats() {
        ::std::cout << std::flush
                    << "Number of neurons: " << this->neurons->size() << ::std::endl
                    << "Number of connections: " << this->connection_list->size() << ::std::endl
                    << "Number of active weights: " << this->connection_weights->size() << ::std::endl
                    << "Number of active biases: " << this->neuron_biases->size() << ::std::endl;

        if(this->normalization_strategy) {
            ::std::cout << std::flush
                        << "Normalization strategy maximum value: "
                        << this->normalization_strategy->get_max_value() << std::endl
                        << "Normalization strategy minimum value: "
                        << this->normalization_strategy->get_min_value()
                        << std::endl;
        }
    }

    void NeuralNetwork::write_stats(std::string file_path) {
        std::ofstream ofs(file_path);

        if(!ofs.is_open()) {
            THROW_RUNTIME_ERROR("File " + file_path + " can not be opened!");
        }

        ofs << "Number of neurons: " << this->neurons->size() << ::std::endl
            << "Number of connections: " << this->connection_list->size() << ::std::endl
            << "Number of active weights: " << this->connection_weights->size() << ::std::endl
            << "Number of active biases: " << this->neuron_biases->size() << ::std::endl;

        if(this->normalization_strategy) {
            ofs << "Normalization strategy maximum value: "
                << this->normalization_strategy->get_max_value() << std::endl
                << "Normalization strategy minimum value: "
                << this->normalization_strategy->get_min_value()
                << std::endl;
        }

        ofs.close();
    }

    void NeuralNetwork::write_stats(std::ofstream* file_path) {
        *file_path << "Number of neurons: " << this->neurons->size() << ::std::endl
                   << "Number of connections: " << this->connection_list->size() << ::std::endl
                   << "Number of active weights: " << this->connection_weights->size() << ::std::endl
                   << "Number of active biases: " << this->neuron_biases->size() << ::std::endl;