ErrorFunctions.cpp

#include <vector>
#include <map>
#include <utility>
#include <cmath>
#include <sstream>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int_distribution.hpp>



#include "ErrorFunctions.h"
#include "exceptions.h"
#include "message.h"
#include "../mpi_wrapper.h"

namespace lib4neuro {

    size_t ErrorFunction::get_dimension() {
        return this->dimension;
    }

    std::vector<NeuralNetwork*>& ErrorFunction::get_nets() {
        return nets;
    }

    DataSet* ErrorFunction::get_dataset() const {
        return this->ds;
    }

    void ErrorFunction::set_dataset(DataSet* ds) {
        this->ds = ds;
    }

    void ErrorFunction::divide_data_train_test(double percent_test) {
        size_t ds_size = this->ds->get_n_elements();

        /* Store the full data set */
        this->ds_full = this->ds;

        /* Choose random subset of the DataSet for training and the remaining part for validation */
        boost::random::mt19937                    gen;
        boost::random::uniform_int_distribution<> dist(0,
                                                       ds_size - 1);

        size_t test_set_size = ceil(ds_size * percent_test);

        std::vector<unsigned int> test_indices;
        test_indices.reserve(test_set_size);
        for (size_t i = 0; i < test_set_size; i++) {
            test_indices.emplace_back(dist(gen));
        }
        std::sort(test_indices.begin(),
                  test_indices.end(),
                  std::greater<unsigned int>());

        std::vector<std::pair<std::vector<double>, std::vector<double>>> test_data, train_data;

        /* Copy all the data to train_data */
        for (auto e : *this->ds_full->get_data()) {
            train_data.emplace_back(e);
        }

        /* Move the testing data from train_data to test_data */
        for (auto ind : test_indices) {
            test_data.emplace_back(train_data.at(ind));
        }
        for(auto ind : test_indices) {
            train_data.erase(train_data.begin() + ind);
        }

        /* Re-initialize data set for training */
        this->ds = new DataSet(&train_data,
                               this->ds_full->get_normalization_strategy());

        /* Initialize test data */
        this->ds_test = new DataSet(&test_data,
                                    this->ds_full->get_normalization_strategy());
    }

    size_t ErrorFunction::divide_data_worst_subset(
            std::vector<size_t> &subset_indices,
            std::vector<int> &active_subset,
            std::vector<float> &entry_errors,
			size_t expansion_len,
			float tolerance
    ) {
        if( this->ds_full == nullptr ){
            this->ds_full = this->ds;
        }
        size_t ds_size = this->ds_full->get_n_elements();

        if( entry_errors.size() != ds_size ){
            entry_errors.resize( ds_size );
        }

        if( active_subset.size() != ds_size ){
            active_subset.resize( ds_size );
            std::fill(active_subset.begin(), active_subset.end(), 0);
        }

        std::vector<double> error_vector( this->get_n_outputs());
        for( size_t i = 0; i < ds_size; ++i ) {
            entry_errors[ i ] = this->eval_single_item_by_idx( i, nullptr, error_vector );
        }

        std::vector<std::pair<std::vector<double>, std::vector<double>>> train_set;
		std::map<double, size_t> expanded_set;
		for( size_t i = 0; i < ds_size; ++i ){
            if( active_subset[ i ] > 0 ){
                continue;
            }

			// expanded_set[ entry_errors[ i ] ] = i;
			expanded_set.insert(std::pair<double, size_t>(entry_errors[ i ], i));

			if( expanded_set.size() > expansion_len ){
				expanded_set.erase( expanded_set.begin( ) );
			}
		}

		subset_indices.clear( );
		for( size_t i = 0; i < ds_size; ++i ){
			if( entry_errors[ i ] < tolerance ){
				/* entry is very well learned, we leave it aside, for now */
				active_subset[ i ] = 0;
			}

			if( active_subset[ i ] > 0 ){
				subset_indices.push_back( i );
			}
		}


        double max_error = -1.0;
        size_t max_error_entry_idx = 0;
		for( auto el: expanded_set ){
			max_error = el.first;
			max_error_entry_idx = el.second;

			subset_indices.push_back( el.second );
		}


        for( auto el: subset_indices ){
            train_set.emplace_back( this->ds_full->get_data( )->at( el ) );
            active_subset[el]++;

            if(active_subset[el] > 2){
                active_subset[el] = 2;
            }
        }

        if( this->ds != this->ds_full ){
            delete this->ds;
        }
        this->ds = new DataSet(&train_set,
                           this->ds_full->get_normalization_strategy());

        return train_set.size( );

    }

    void ErrorFunction::return_full_data_set_for_training() {
        if (this->ds_test || this->ds != this->ds_full) {
//            delete this->ds;
            this->ds = this->ds_full;
            this->ds_full = nullptr;
        }
    }
	
	void ErrorFunction::add_to_hessian_and_rhs_single(
		arma::Mat<double>& H,
		arma::Col<double>& rhs,
		size_t entry_idx
	){
		
        std::vector<double> partial_error(this->get_n_outputs());
		std::vector<std::vector<double>> jac_loc;
		
		this->nets[0]->get_jacobian( jac_loc, this->ds->get_data()->at( entry_idx ), partial_error );
		
		auto nr = jac_loc.size();
		auto nc = this->get_dimension();

		for (size_t ri = 0; ri < nr; ++ri) {
			auto alpha = partial_error[ri];
			for (size_t ci = 0; ci < nc; ++ci) {
				rhs.at(ci) += alpha * jac_loc[ri][ci];
			}
		}
		
		double rval, cval;
		for( auto d = 0; d < nr; ++d ){
			for( auto ci = 0; ci < nc; ++ci ){
				cval = jac_loc[d][ci];
				for( auto ri = 0; ri < nc; ++ri ){
					rval = jac_loc[d][ri];
					H.at(ri, ci) += rval * cval;
				}