text:
// Copyright (C) 2019 T. Zachary Laine
//
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#include <boost/stl_interfaces/container_interface.hpp>

#include "ill_formed.hpp"

#include <gtest/gtest.h>

#include <array>
#include <deque>
#include <vector>

namespace bsi = boost::stl_interfaces::v2;

// Just like std::array, except for the 0-size specialization, and the fact
// that the base class makes brace-initialization wonky.
template<typename T, std::size_t N>
struct array : bsi::container_interface<array<T, N>>
{
    using value_type = T;
    using pointer = T *;
    using const_pointer = T const *;
    using reference = value_type &;
    using const_reference = value_type const &;
    using size_type = std::size_t;
    using difference_type = std::ptrdiff_t;
    using iterator = T *;
    using const_iterator = T const *;
    using reverse_iterator = boost::stl_interfaces::reverse_iterator<iterator>;
    using const_reverse_iterator =
        boost::stl_interfaces::reverse_iterator<const_iterator>;

    void fill(T const & x) { std::fill(begin(), end(), x); }

    iterator begin() noexcept { return elements_; }
    iterator end() noexcept { return elements_ + N; }

    size_type max_size() const noexcept { return N; }

    void swap(array & other)
    {
        using std::swap;
        T * element = elements_;
        for (auto & x : other) {
            swap(*element++, x);
        }
    }

    using base_type = bsi::container_interface<array<T, N>>;
    using base_type::begin;
    using base_type::end;

    T elements_[N];
};

using arr_type = array<int, 5>;

TEST(array, comparisons)
{
    arr_type sm;
    sm[0] = 1;
    sm[1] = 2;
    sm[2] = 3;
    sm[3] = 0;
    sm[4] = 0;
    arr_type md;
    md[0] = 1;
    md[1] = 2;
    md[2] = 3;
    md[3] = 4;
    md[4] = 0;
    arr_type lg;
    lg[0] = 1;
    lg[1] = 2;
    lg[2] = 3;
    lg[3] = 4;
    lg[4] = 5;

    EXPECT_TRUE(sm == sm);
    EXPECT_FALSE(sm == md);
    EXPECT_FALSE(sm == lg);

    EXPECT_FALSE(sm != sm);
    EXPECT_TRUE(sm != md);
    EXPECT_TRUE(sm != lg);

    EXPECT_FALSE(sm < sm);
    EXPECT_TRUE(sm < md);
    EXPECT_TRUE(sm < lg);

    EXPECT_TRUE(sm <= sm);
    EXPECT_TRUE(sm <= md);
    EXPECT_TRUE(sm <= lg);

    EXPECT_FALSE(sm > sm);
    EXPECT_FALSE(sm > md);
    EXPECT_FALSE(sm > lg);

    EXPECT_TRUE(sm >= sm);
    EXPECT_FALSE(sm >= md);
    EXPECT_FALSE(sm >= lg);

    EXPECT_FALSE(md == sm);
    EXPECT_TRUE(md == md);
    EXPECT_FALSE(md == lg);

    EXPECT_FALSE(md < sm);
    EXPECT_FALSE(md < md);
    EXPECT_TRUE(md < lg);

    EXPECT_FALSE(md <= sm);
    EXPECT_TRUE(md <= md);
    EXPECT_TRUE(md <= lg);

    EXPECT_TRUE(md > sm);
    EXPECT_FALSE(md > md);
    EXPECT_FALSE(md > lg);

    EXPECT_TRUE(md >= sm);
    EXPECT_TRUE(md >= md);
    EXPECT_FALSE(md >= lg);

    EXPECT_FALSE(lg == sm);
    EXPECT_FALSE(lg == md);
    EXPECT_TRUE(lg == lg);

    EXPECT_FALSE(lg < sm);
    EXPECT_FALSE(lg < md);
    EXPECT_FALSE(lg < lg);

    EXPECT_FALSE(lg <= sm);
    EXPECT_FALSE(lg <= md);
    EXPECT_TRUE(lg <= lg);

    EXPECT_TRUE(lg > sm);
    EXPECT_TRUE(lg > md);
    EXPECT_FALSE(lg > lg);

    EXPECT_TRUE(lg >= sm);
    EXPECT_TRUE(lg >= md);
    EXPECT_TRUE(lg >= lg);
}

TEST(array, swap)
{
    {
        arr_type v1;
        v1[0] = 3;
        v1[1] = 4;
        v1[2] = 0;
        v1[3] = 0;
        v1[4] = 0;
        arr_type v2;
        v2[0] = 4;
        v2[1] = 3;
        v2[2] = 0;
        v2[3] = 0;
        v2[4] = 0;

        arr_type const v1_copy = v1;
        arr_type const v2_copy = v2;

        static_assert(std::is_same<decltype(v1.swap(v2)), void>::value, "");
        static_assert(std::is_same<decltype(swap(v1, v2)), void>::value, "");

        v1.swap(v2);

        EXPECT_EQ(v1, v2_copy);
        EXPECT_EQ(v2, v1_copy);
    }

    {
        arr_type v1;
        v1[0] = 3;
        v1[1] = 4;
        v1[2] = 0;
        v1[3] = 0;
        v1[4] = 0;
        arr_type v2;
        v2[0] = 4;
        v2[1] = 3;
        v2[2] = 0;
        v2[3] = 0;
        v2[4] = 0;

        arr_type const v1_copy = v1;
        arr_type const v2_copy = v2;

        swap(v1, v2);

        EXPECT_EQ(v1, v2_copy);
        EXPECT_EQ(v2, v1_copy);
    }
}

template<typename Iter>
using writable_iter_t = decltype(
    *std::declval<Iter>() =
        std::declval<typename std::iterator_traits<Iter>::value_type>());

static_assert(
    !ill_formed<writable_iter_t, decltype(std::declval<arr_type &>().begin())>::
        value,
    "");
static_assert(
    !ill_formed<writable_iter_t, decltype(std::declval<arr_type &>().end())>::
        value,
    "");

static_assert(
    ill_formed<
        writable_iter_t,
        decltype(std::declval<arr_type const &>().begin())>::value,
    "");
static_assert(
    ill_formed<
        writable_iter_t,
        decltype(std::declval<arr_type const &>().end())>::value,
    "");
static_assert(
    ill_formed<writable_iter_t, decltype(std::declval<arr_type &>().cbegin())>::
        value,
    "");
static_assert(
    ill_formed<writable_iter_t, decltype(std::declval<arr_type &>().cend())>::
        value,
    "");

static_assert(
    ill_formed<
        writable_iter_t,
        decltype(std::declval<arr_type const &>().rbegin())>::value,
    "");
static_assert(
    ill_formed<
        writable_iter_t,
        decltype(std::declval<arr_type const &>().rend())>::value,
    "");
static_assert(
    ill_formed<
        writable_iter_t,
        decltype(std::declval<arr_type &>().crbegin())>::value,
    "");
static_assert(
    ill_formed<writable_iter_t, decltype(std::declval<arr_type &>().crend())>::
        value,
    "");

TEST(array, iterators)
{
    arr_type v0;
    v0[0] = 3;
    v0[1] = 2;
    v0[2] = 1;
    v0[3] = 0;
    v0[4] = 0;

    {
        arr_type v = v0;

        static_assert(
            std::is_same<decltype(v.begin()), arr_type::iterator>::value, "");
        static_assert(
            std::is_same<decltype(v.end()), arr_type::iterator>::value, "");
        static_assert(
            std::is_same<decltype(v.cbegin()), arr_type::const_iterator>::value,
            "");
        static_assert(
            std::is_same<decltype(v.cend()), arr_type::const_iterator>::value,
            "");
        static_assert(
            std::is_same<decltype(v.rbegin()), arr_type::reverse_iterator>::
                value,
            "");
        static_assert(
            std::is_same<decltype(v.rbegin()), arr_type::reverse_iterator>::
                value,
            "");
        static_assert(
            std::is_same<
                decltype(v.crbegin()),
                arr_type::const_reverse_iterator>::value,
            "");
        static_assert(
            std::is_same<
                decltype(v.crbegin()),
                arr_type::const_reverse_iterator>::value,
            "");

        std::array<int, 5> const a = {{3, 2, 1, 0, 0}};
        std::array<int, 5> const ra = {{0, 0, 1, 2, 3}};

        EXPECT_TRUE(std::equal(v.begin(), v.end(), a.begin(), a.end()));
        EXPECT_TRUE(std::equal(v.cbegin(), v.cend(), a.begin(), a.end()));

        EXPECT_TRUE(std::equal(v.rbegin(), v.rend(), ra.begin(), ra.end()));
        EXPECT_TRUE(std::equal(v.crbegin(), v.crend(), ra.begin(), ra.end()));

        arr_type v2;
        v2[0] = 8;
        v2[1] = 2;
        v2[2] = 1;
        v2[3] = 0;
        v2[4] = 9;

        *v.begin() = 8;
        *v.rbegin() = 9;
        EXPECT_EQ(v, v2);
    }

    {
        arr_type const v = v0;

        static_assert(
            std::is_same<decltype(v.begin()), arr_type::const_iterator>::value,
            "");
        static_assert(
            std::is_same<decltype(v.end()), arr_type::const_iterator>::value,
            "");
        static_assert(
            std::is_same<decltype(v.cbegin()), arr_type::const_iterator>::value,
            "");
        static_assert(
            std::is_same<decltype(v.cend()), arr_type::const_iterator>::value,
            "");
        static_assert(
            std::is_same<
                decltype(v.rbegin()),
                arr_type::const_reverse_iterator>::value,
            "");
        static_assert(
            std::is_same<
                decltype(v.rbegin()),
                arr_type::const_reverse_iterator>::value,
            "");
        static_assert(
            std::is_same<
                decltype(v.crbegin()),
                arr_type::const_reverse_iterator>::value,
            "");
        static_assert(
            std::is_same<
                decltype(v.crbegin()),
                arr_type::const_reverse_iterator>::value,
            "");

        std::array<int, 5> const a = {{3, 2, 1, 0, 0}};
        std::array<int, 5> const ra = {{0, 0, 1, 2, 3}};

        EXPECT_TRUE(std::equal(v.begin(), v.end(), a.begin(), a.end()));
        EXPECT_TRUE(std::equal(v.cbegin(), v.cend(), a.begin(), a.end()));

        EXPECT_TRUE(std::equal(v.rbegin(), v.rend(), ra.begin(), ra.end()));
        EXPECT_TRUE(std::equal(v.crbegin(), v.crend(), ra.begin(), ra.end()));
    }
}

template<
    typename Container,
    typename ValueType = typename Container::value_type>
using lvalue_push_front_t = decltype(
    std::declval<Container &>().push_front(std::declval<ValueType const &>()));
template<
    typename Container,
    typename ValueType = typename Container::value_type>
using rvalue_push_front_t = decltype(std::declval<Container &>().push_front(0));
template<typename Container>
using pop_front_t = decltype(std::declval<Container &>().pop_front());

static_assert(ill_formed<lvalue_push_front_t, arr_type>::value, "");
static_assert(ill_formed<rvalue_push_front_t, arr_type>::value, "");
static_assert(ill_formed<pop_front_t, arr_type>::value, "");

using std_deq_int = std::deque<int>;

static_assert(!ill_formed<lvalue_push_front_t, std_deq_int>::value, "");
static_assert(!ill_formed<rvalue_push_front_t, std_deq_int>::value, "");
static_assert(!ill_formed<pop_front_t, std_deq_int>::value, "");

template<
    typename Container,
    typename ValueType = typename Container::value_type>
using lvalue_push_back_t = decltype(
    std::declval<Container &>().push_back(std::declval<ValueType const &>()));
template<
    typename Container,
    typename ValueType = typename Container::value_type>
using rvalue_push_back_t = decltype(std::declval<Container &>().push_back(0));
template<typename Container>
using pop_back_t = decltype(std::declval<Container &>().pop_back());

static_assert(ill_formed<lvalue_push_back_t, arr_type>::value, "");
static_assert(ill_formed<rvalue_push_back_t, arr_type>::value, "");
static_assert(ill_formed<pop_back_t, arr_type>::value, "");

using std_vec_int = std::vector<int>;

static_assert(!ill_formed<lvalue_push_back_t, std_vec_int>::value, "");
static_assert(!ill_formed<rvalue_push_back_t, std_vec_int>::value, "");
static_assert(!ill_formed<pop_back_t, std_vec_int>::value, "");

TEST(array, front_back)
{
    {
        arr_type v;
        v[0] = 0;
        v[1] = 0;
        v[2] = 0;
        v[3] = 0;
        v[4] = 0;

        static_assert(std::is_same<decltype(v.front()), int &>::value, "");
        static_assert(std::is_same<decltype(v.back()), int &>::value, "");

        v.front() = 9;
        v.back() = 8;
        EXPECT_EQ(v[0], v.front());
        EXPECT_EQ(v[4], v.back());
    }

    {
        arr_type v0;
        v0[0] = 3;
        v0[1] = 0;
        v0[2] = 2;
        v0[3] = 0;
        v0[4] = 1;

        arr_type const v = v0;
        EXPECT_EQ(v.front(), 3);
        EXPECT_EQ(v.back(), 1);

        static_assert(
            std::is_same<decltype(v.front()), int const &>::value, "");
        static_assert(std::is_same<decltype(v.back()), int const &>::value, "");
    }
}

TEST(array, index_at)
{
    arr_type v0;
    v0[0] = 3;
    v0[1] = 2;
    v0[2] = 1;
    v0[3] = 0;
    v0[4] = 0;

    {
        arr_type v = v0;
        EXPECT_EQ(v[0], 3);
        EXPECT_EQ(v[1], 2);
        EXPECT_EQ(v[2], 1);
        EXPECT_NO_THROW(v.at(0));
        EXPECT_NO_THROW(v.at(1));
        EXPECT_NO_THROW(v.at(2));
        EXPECT_THROW(v.at(5), std::out_of_range);

        static_assert(std::is_same<decltype(v[0]), int &>::value, "");
        static_assert(std::is_same<decltype(v.at(0)), int &>::value, "");

        v[0] = 8;
        v.at(1) = 9;
        EXPECT_EQ(v[0], 8);
        EXPECT_EQ(v[1], 9);
    }

    {
        arr_type const v = v0;
        EXPECT_EQ(v[0], 3);
        EXPECT_EQ(v[1], 2);
        EXPECT_EQ(v[2], 1);
        EXPECT_NO_THROW(v.at(0));
        EXPECT_NO_THROW(v.at(1));
        EXPECT_NO_THROW(v.at(2));
        EXPECT_THROW(v.at(5), std::out_of_range);

        static_assert(std::is_same<decltype(v[0]), int const &>::value, "");
        static_assert(std::is_same<decltype(v.at(0)), int const &>::value, "");
    }
}

template<typename Container>
using resize_t = decltype(std::declval<Container &>().resize(0));

static_assert(ill_formed<resize_t, arr_type>::value, "");

static_assert(!ill_formed<resize_t, std_vec_int>::value, "");

template<
    typename Container,
    typename ValueType = typename Container::value_type>
using lvalue_insert_t = decltype(std::declval<Container &>().insert(
    std::declval<Container &>().begin(), std::declval<ValueType const &>()));
template<
    typename Container,
    typename ValueType = typename Container::value_type>
using rvalue_insert_t = decltype(std::declval<Container &>().insert(
    std::declval<Container &>().begin(), std::declval<ValueType &&>()));
template<
    typename Container,
    typename ValueType = typename Container::value_type>
using insert_n_t = decltype(std::declval<Container &>().insert(
    std::declval<Container &>().begin(), 2, std::declval<ValueType const &>()));
template<
    typename Container,
    typename ValueType = typename Container::value_type>
using insert_il_t = decltype(std::declval<Container &>().insert(
    std::declval<Container &>().begin(), std::initializer_list<int>{}));

static_assert(ill_formed<lvalue_insert_t, arr_type>::value, "");
static_assert(ill_formed<rvalue_insert_t, arr_type>::value, "");

// TODO: Broken.  Adding the proper constraint ICE's GCC 8 and 9.
#if 0
static_assert(ill_formed<insert_n_t, arr_type>::value, "");
#endif

static_assert(ill_formed<insert_il_t, arr_type>::value, "");

static_assert(!ill_formed<lvalue_insert_t, std_vec_int>::value, "");
static_assert(!ill_formed<rvalue_insert_t, std_vec_int>::value, "");
static_assert(!ill_formed<insert_n_t, std_vec_int>::value, "");
static_assert(!ill_formed<insert_il_t, std_vec_int>::value, "");

template<typename Container>
using erase_t = decltype(
    std::declval<Container &>().erase(std::declval<Container &>().begin()));

static_assert(ill_formed<erase_t, arr_type>::value, "");

static_assert(!ill_formed<erase_t, std_vec_int>::value, "");

template<typename Container>
using assign_t = decltype(std::declval<Container &>().assign(
    std::declval<int *>(), std::declval<int *>()));
template<typename Container>
using assign_n_t = decltype(std::declval<Container &>().assign(5, 5));
template<typename Container>
using assign_il_t =
    decltype(std::declval<Container &>().assign(std::initializer_list<int>{}));
template<typename Container>
using assignment_operator_il_t =
    decltype(std::declval<Container &>() = std::initializer_list<int>{});

static_assert(ill_formed<assign_t, arr_type>::value, "");
static_assert(ill_formed<assign_n_t, arr_type>::value, "");
static_assert(ill_formed<assign_il_t, arr_type>::value, "");
static_assert(ill_formed<assignment_operator_il_t, arr_type>::value, "");

static_assert(!ill_formed<assign_t, std_vec_int>::value, "");
static_assert(!ill_formed<assign_n_t, std_vec_int>::value, "");
static_assert(!ill_formed<assign_il_t, std_vec_int>::value, "");
static_assert(!ill_formed<assignment_operator_il_t, std_vec_int>::value, "");

template<typename Container>
using clear_t = decltype(std::declval<Container &>().clear());

static_assert(ill_formed<clear_t, arr_type>::value, "");

static_assert(!ill_formed<clear_t, std_vec_int>::value, "");


text:
#include <iostream>
#include <armadillo>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <chrono>
#include <ctime>
#include "helper_functions.h"
#include "solve_network.h"
#include "solve_dde.h"
#include "backprop.h"
#include "numerical_gradient.h"
#include "katana_get_params.hh"
#include "global_constants.h"

using namespace std;
using namespace arma;
using namespace globalconstants;

int main(int argc, char const *argv[])
{
	// for cpu time measuring
	clock_t start_overall = clock();
	double cumulative_solve_time = 0.0;
	double cumulative_backprop_time = 0.0;

	arma_rng::set_seed_random();

	// ### ### ### --- OPTIONS --- ### ### ###

	// file name of the text file where parameters and results will be saved
	string results_file_name = katana::getCmdOption(argv, argv + argc, "-filename", "results.txt");
	// print message to results file
	string print_msg = "Simulation of the deep learning delay system.";

	// task, possible options are "MNIST", "Fashion-MNIST", "CIFAR-10", "SVHN"
	string task = "Fashion-MNIST-denoising";
	// modify global_constants.h accordingly!

	// number of example images to save
	int save_examples = katana::getCmdOption(argv, argv + argc, "-save_examples", 0);

	// print weights (and diagonals) after each training epoch to text files in weights folder
	bool print_weights_to_file = false;

	// record data for time-signal video (or plot), only for MNIST or Fashion-MNIST
	bool record_data = katana::getCmdOption_bool(argv, argv + argc, "-record_data", false);
	int number_of_examples = katana::getCmdOption(argv, argv + argc, "-record_examples", 10);

	// If the following option is true, the program will not train the machine learning system.
	// Instead it will just choose diag indices n_prime according to the given paramters
	// and save them in the text file diag.txt. This will be done 6 times, i.e. the text file
	// will contain 6 lines of diag indices, which can be used for 6-fold cross validation
	bool make_diags = katana::getCmdOption_bool(argv, argv + argc, "-make_diags", false);
	// If the following option is true, the program will create initial weight files and abort.
	bool save_init_weights = katana::getCmdOption_bool(argv, argv + argc, "-save_init_weights", false);
	// If the following option is true, the program will load the initial weights from files instead of creating random weights.
	// Note that the corresponding diag file must be loaded as well.
	bool init_weights_from_file = katana::getCmdOption_bool(argv, argv + argc, "-init_weights_from_file", false);

	// option for system simulation
	// "dde_ibp":	semi-analytic heun method with modified trapezoidal rule employing integration-by-parts
	// "dde_heun":	semi-analytic heun method with standard trapezoidal rule
	// "network":	solves the equivalent network equations
	// "network_decoupled": solves network ignoring the additional direct linear connections of neighboring nodes and layers
	string system_simu = katana::getCmdOption(argv, argv + argc, "-system_simu", "dde_ibp");

	// option for gradient computation
	// backprop_standard: newly derived backpropagation algorithm for the deep learning delay system (was derived for equivalent network)
	// backprop_classic: classical backpropagation algorithm, ignores the additional direct linear connections of neighboring nodes and layers
	string grad_comp = katana::getCmdOption(argv, argv + argc, "-grad_comp", "backprop_standard");

	// gradient_check: compute numerical gradient (for dde_ibp) for comparision
	bool gradient_check = katana::getCmdOption_bool(argv, argv + argc, "-gradient_check", false);
	double epsilon_gradient_check = 1e-9;
	bool print_gradients_to_file = false;

	// ### ### ### --- PARAMETERS --- ### ### ###

	// M and P are defined in global_constants.h

	// task parameter
	double noise_sigma = 1.0;

	// ... for the system/network architectur:

	int N = katana::getCmdOption(argv, argv + argc, "-N", 100);  // number of nodes per hidden layer
	int L = katana::getCmdOption(argv, argv + argc, "-L", 2);  // number of hidden layers
	int D = katana::getCmdOption(argv, argv + argc, "-D", 50);  // number of delays (NOTE: must be the correct number also if diag_method == "from_file")

	double theta = katana::getCmdOption(argv, argv + argc, "-theta", 0.5);  // node separation
	double alpha = -1.0;  // factor in linear part of delay system

	string diag_method = katana::getCmdOption(argv, argv + argc, "-diag_method", "uniform");  // method to choose diagonals
	// "uniform": uniform distribution with the additional condition that there is at least one n'_d < 0 and one n'_d > 0
	// "equi_dist": diagonals with equal distances in between, if D is an odd number: main diagonal and same number of upper and lower,
	//              if D is an even number and diag_distance is even: centered, e.g. -15, -5, 5, 15,
	//              if D is even and diag_distance is odd: centered around -0.5, e.g. -8, -3, 2, 7.
	// "from_file": diag indices are taken from a given text file
	// For a banded matrix as connection matrix choose equi_dist and diag_distance = 1
	// For a full connection matrix choose D = 2 * N - 1 (and either uniform or equi_dist with diag_distance = 1)
	int diag_distance = katana::getCmdOption(argv, argv + argc, "-diag_distance", 0);  // distance between diagonals if diag_method == "equi_dist",
							// minimum distance between diagonals if diag_method is "uniform", otherwise ignored
							// Do not choose diag_distance too large. If (D - 1) * (2 * diag_distance - 1) >= 2 * (N - diag_margin) - 1 and method is uniform, the program will abort.
	string diag_file_path = katana::getCmdOption(argv, argv + argc, "-diag_file_path", "diag.txt"); // path to text file containing D interger numbers n_prime_d, ignored if method is not "from_file"

	// ... for the training:
	int number_of_epochs = katana::getCmdOption(argv, argv + argc, "-number_of_epochs", 100);
	double eta_0 = katana::getCmdOption(argv, argv + argc, "-eta0", 0.001);
	double eta_1 = katana::getCmdOption(argv, argv + argc, "-eta1", 1000.0);  // learning rate eta = min(eta_0, eta_1 / step)
	bool pixel_shift = katana::getCmdOption_bool(argv, argv + argc, "-pixel_shift", false);  // on-off switch for training input random 1-pixel shift
	int max_pixel_shift = katana::getCmdOption(argv, argv + argc, "-max_pixel_shift", 1);
	bool input_noise = katana::getCmdOption_bool(argv, argv + argc, "-training_noise", false);  // on-off switch for training input gaussian noise
	double training_noise_sigma = katana::getCmdOption(argv, argv + argc, "-sigma", 0.01);  // standard deviation of gaussian noise to disturb training input
	bool rotation = katana::getCmdOption_bool(argv, argv + argc, "-rotation", false);
	double max_rotation_degrees = katana::getCmdOption(argv, argv + argc, "-max_rotation_degrees", 15.0);
	bool horizontal_flip = katana::getCmdOption_bool(argv, argv + argc, "-flip", false);  // only for CIFAR-10
	// dropout
	double dropout_rate = katana::getCmdOption(argv, argv + argc, "-dropout", 0.0);

	// ... for weight initialization:
	double initial_input_weigt_radius = sqrt(6.0 / ((double)D/2.0 + (double)M + 1.0));
	double initial_hidden_weigt_radius =  sqrt(6.0 / ((double)D + 1.0));
	double initial_output_weigt_radius =  sqrt(6.0 / ((double)D/2.0 + (double)P + 1.0));

	// ... for numerics:

	int N_h = max(32 ,(int)(16 * theta));  // computation steps per virtual node for solving DDE

	// computational precision for sums with exp(alpha * theta * n) factor in summands
	// in the functions "get_deltas" and "get_gradient" in "backprop.cpp".
	// exp_precision = -35.0 means that terms, where exponential factor
	// is smaller than exp(-35), are ignored.
	// Since exp(-35) is approximately 6.3e-16, the gradient will still be computed with double precision.
	double exp_precision = -35.0;

	// ### ### ### --- ETC. --- ### ### ###

	// make diag text file and end program if make_diags option is true
	if (make_diags){
		ofstream diag_file;
		diag_file.open("diag.txt");
		for (int i = 0; i < 6; ++i){
			vector<int> diag_indices = get_diag_indices(N, D, diag_method, diag_distance, diag_file_path, 0);
			for (int d = 0; d < D - 1; ++d){
				diag_file << diag_indices[d] << " ";
			}
			diag_file << diag_indices[D - 1] << endl;
		}
		diag_file.close();
		cout << "The option make_diags was true." << endl;
		cout << "Program made (or overrode) diag.txt" << endl;
		return 0;
	}

	// make initial weight files and diag file and end program if make_diags save_init_weights option is true
	if (save_init_weights){
		vector<int> diag_indices;
		ofstream diag_file;
		diag_file.open("diag.txt");
		for (int i = 0; i < 6; ++i){
			diag_indices = get_diag_indices(N, D, diag_method, diag_distance, diag_file_path, 0);
			for (int d = 0; d < D - 1; ++d){
				diag_file << diag_indices[d] << " ";
			}
			diag_file << diag_indices[D - 1] << endl;
		}
		diag_file.close();
		cout << "The option save_init_weights was true." << endl;
		cout << "Program made (or overrode) diag.txt" << endl;
		mat input_weights(N, M + 1);
		mat output_weights(P, N + 1);
		cube hidden_weights(L - 1, N, N + 1);
		initialize_weights(input_weights, hidden_weights, output_weights, D, L, N, diag_indices,
						   initial_input_weigt_radius, initial_hidden_weigt_radius, initial_output_weigt_radius, true);
		cout << "Save initial weights and end program because the option save_init_weights is true." << endl;
		return 0;
	}

	// The following lines are to create a look up table "exp_table"
	// which contains the values exp(alpha * theta * n)
	// which are often needed by the functions "get_deltas" and "get_gradient" in "backprop.cpp".
	vector<double> exp_table;
	int n = 0;
	while (alpha * theta * double(n) >= - 35.0){
		double exp_value = exp(alpha * theta * double(n));
		exp_table.push_back(exp_value);
		++n;
	}

	// make results file and print information about parameters
	print_parameters(results_file_name, print_msg, task,
					system_simu, grad_comp, gradient_check,
					N, L, D, theta, alpha,
					diag_method, diag_distance, diag_file_path,
					number_of_epochs, eta_0, eta_1,
					pixel_shift, input_noise, training_noise_sigma,
					N_h, exp_precision);

	// determine data directory
	string data_dir = "data-Fashion-MNIST";

	// read image data from files to arrays:
	cube train_images(number_of_training_batches, training_batch_size, M);
	mat test_images(test_batch_size, M);
	int train_labels[number_of_training_batches][training_batch_size];
	int test_labels[test_batch_size];
	read_files(train_images, test_images, train_labels, test_labels, data_dir);

	// for video:
	vector<int> step_indices_for_video;
	if (record_data){
		fstream index_file;
		string index_string;
		index_file.open("video/step_indices.txt", ios::in);
		do {
			getline(index_file, index_string);
			if (index_string != ""){
				int i = stoi(index_string);
				step_indices_for_video.push_back(i);
			}
		} while (index_string != "");
		index_file.close();
	}

	// ### ### ### --- INITIALIZATION --- ### ### ###

	vector<int> training_batch_indices;
	cout << "Begin training." << endl;
	for (int i = 0; i < number_of_training_batches; ++i){
		training_batch_indices.push_back(i);
	}

	// initialize arrays which are used below to store the current system states
	mat activations(L, N);
	mat node_states(L, N);
	double output_activations[P];
	double outputs[P];
	vec g_primes(N);

	// initialize arrays to store deltas and gradient
	mat deltas(L, N);
	double output_deltas[P];
	mat input_weight_gradient(N, M + 1, fill::zeros);
	mat output_weight_gradient(P, N + 1, fill::zeros);
	cube weight_gradient(L - 1, N, N + 1, fill::zeros);

	mat deltas_0(L, N);
	double output_deltas_0[P];
	mat input_weight_gradient_0(N, M + 1, fill::zeros);
	mat output_weight_gradient_0(P, N + 1, fill::zeros);
	cube weight_gradient_0(L - 1, N, N + 1, fill::zeros);

	// initialize arrays to eventually store numerical gradient
	mat num_input_weight_gradient(N, M + 1, fill::zeros);
	mat num_output_weight_gradient(P, N + 1, fill::zeros);
	cube num_weight_gradient(L - 1, N, N + 1, fill::zeros);

	// The function "get_diag_indices" returns D intergers n'_d between - N + 1 and N - 1.
	// The delays are then tau_d = N - n'_d.
	vector<int> diag_indices;
	diag_indices = get_diag_indices(N, D, diag_method, diag_distance, diag_file_path, 0);

	// initialize weights.
	mat input_weights(N, M + 1, fill::zeros);
	mat output_weights(P, N + 1, fill::zeros);
	cube hidden_weights(L - 1, N, N + 1, fill::zeros);
	if (init_weights_from_file){
		load_initial_weights(input_weights, hidden_weights, output_weights);
	} else {
		initialize_weights(input_weights, hidden_weights, output_weights, D, L, N, diag_indices,
					   initial_input_weigt_radius, initial_hidden_weigt_radius, initial_output_weigt_radius, false);
	}

	// weights for test_runs.
	mat input_weights_scaled(N, M + 1, fill::zeros);
	mat output_weights_scaled(P, N + 1, fill::zeros);
	cube hidden_weights_scaled(L - 1, N, N + 1, fill::zeros);

	// dropout mask for weights:
	mat input_weights_mask(N, M + 1, fill::ones);
	cube hidden_weights_mask(L - 1, N, N + 1, fill::ones);
	mat output_weights_mask(P, N + 1, fill::ones);

	// ### ### ### --- STOCHASTIC GRADIENT DESCENT TRAINING --- ### ### ###

	// vectors to track training accuracy and validition accuracy (and eventually cosine similarity):
	vector<double> training_accuracy_vector;
	vector<double> accuracy_vector;
	vector<double> similarity_vector;

	// vector and variables to measure and save cpu time needed for each epoch
	vector<double> time_vector;
	clock_t start;
	clock_t ende;
	double epoch_time;

	// loop over training epochs:
	for (int epoch = 0; epoch < number_of_epochs; ++epoch){

		start = clock();

		// make vector with randomly shuffled indices between 0 and 59999 for each epoch of the stochastic gradient descent
		vector<int> index_vector;
		for (int index = 0; index < number_of_training_batches*training_batch_size; ++index){
			index_vector.push_back(index);
		}
		shuffle(begin(index_vector), std::end(index_vector), rng);

		// loop over training steps:
		int step_index = 0;
		for (int index : index_vector){
			++step_index;

			// choose random dropout nodes
			input_weights_mask = mat(N, M + 1, fill::ones);
			hidden_weights_mask = cube(L - 1, N, N + 1, fill::ones);
			output_weights_mask = mat(P, N + 1, fill::ones);
			// input layer
			for (int m = 0; m < M; ++m){
				double random_num = uniform(0.0, 1.0);
				if (random_num < dropout_rate){
					for (int n = 0; n < N; ++n){
						input_weights_mask(n, m) = 0.0;
					}
				}
			}
			// first hidden layer
			if (L > 1){
				for (int n = 0; n < N; ++n){
					double random_num = uniform(0.0, 1.0);
					if (random_num < dropout_rate){
						for (int m = 0; m < M + 1; ++m){
							input_weights_mask(n, m) = 0.0;
						}
						for (int j = 0; j < N; ++j){
							hidden_weights_mask(0, j, n) = 0.0;
						}
					}
				}
			}
			// hidden layers except first and last
			for (int l = 1; l < L - 1; ++l){
				for (int n = 0; n < N; ++n){
					double random_num = uniform(0.0, 1.0);
					if (random_num < dropout_rate){
						for (int i = 0; i < N + 1; ++i){
							hidden_weights_mask(l - 1, n, i) = 0.0;
						}
						for (int j = 0; j < N; ++j){
							hidden_weights_mask(l, j, n) = 0.0;
						}
					}
				}
			}
			// last hidden layer
			if (L > 1){
				for (int n = 0; n < N; ++n){
					double random_num = uniform(0.0, 1.0);
					if (random_num < dropout_rate){
						for (int i = 0; i < N + 1; ++i){
							hidden_weights_mask(L - 2, n, i) = 0.0;
						}
						for (int p = 0; p < P; ++p){
							output_weights_mask(p, n) = 0.0;
						}
					}
				}
			}

			// record data for video: step 0
			if (record_data && epoch == 0 && step_index == 1){
				int rec_step = 0;
				cout << "record data after step: " << rec_step << endl;
				//save weights
				string input_file_name = "video/weights_input_step_" + to_string(rec_step) + ".txt";
				string hidden_file_name = "video/weights_hidden_step_" + to_string(rec_step) + ".txt";
				string output_file_name = "video/weights_output_step_" + to_string(rec_step) + ".txt";
				input_weights.save(input_file_name, csv_ascii);
				hidden_weights.save(hidden_file_name, raw_ascii);
				output_weights.save(output_file_name, csv_ascii);
				//run system with examples from validation batch and save input, output and x_states
				for (int example_index = 0; example_index < number_of_examples; ++example_index){
					// select image as input
					vec input_data = test_images.row(example_index).t();
					int label = test_labels[index];
					// run system and record data for x_states
					solve_dde_ibp(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights, hidden_weights, output_weights, diag_indices, theta, alpha,
							  N, L, N_h, example_index+1, rec_step);
					// save input and output
					string input_file_name = "video/vector_input_step_" + to_string(rec_step) + "_example_" + to_string(example_index + 1) + ".txt";
					input_data.save(input_file_name, csv_ascii);
					string output_file_name = "video/vector_output_step_" + to_string(rec_step) + "_example_" + to_string(example_index + 1) + ".txt";
					ofstream output_file;
					output_file.open(output_file_name);
					for (double y_p : outputs){
						output_file << y_p << endl;
					}
					output_file.close();
				}
			}

			double eta = learning_rate(epoch, step_index, eta_0, eta_1);

			// select image as input
			div_t div_result = div(index, training_batch_size);
			int batch_index = training_batch_indices[div_result.quot];
			int image_index = div_result.rem;
			vec input_data = train_images.tube(batch_index, image_index);
			int label = train_labels[batch_index][image_index];

			// data augmentation or regression
			if (rotation && M == 3072){
				double rotation_degrees = uniform(-max_rotation_degrees, max_rotation_degrees);
				rotation32(input_data, rotation_degrees);
			}

			if (pixel_shift && M == 784){
				pixel_shift28(input_data, max_pixel_shift);
			}
			if (pixel_shift && M == 3072){
				pixel_shift32(input_data, max_pixel_shift);
			}
			if (input_noise){
				vec training_noise = training_noise_sigma * vec(M, fill::randn);
				input_data += training_noise;
			}
			if (horizontal_flip && M == 3072){
				double random_val = uniform(0,1);
				if (random_val > 0.5){
					horizontal_flip32(input_data);
				}
			}

			// get target
			double targets[P];
			for (int p = 0; p < P; ++p){
				targets[p] = input_data(p);
			}

			// add noise to input:
			vec noise = noise_sigma * vec(M, fill::randn);
			input_data += noise;
			for (int m = 0; m < M; ++m){
				if (input_data[m] < 0.0){
					input_data[m] = 0.0;
				}
				if (input_data[m] > 1.0){
					input_data[m] = 1.0;
				}
			}

			// solve the DDE (or network)
			clock_t start_solve = clock();
			if (system_simu == "dde_ibp"){
				solve_dde_ibp(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights%input_weights_mask, hidden_weights%hidden_weights_mask, output_weights%output_weights_mask, diag_indices, theta, alpha,
							  N, L, N_h);
			} else if (system_simu == "dde_heun"){
				solve_dde_heun(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights%input_weights_mask, hidden_weights%hidden_weights_mask, output_weights%output_weights_mask, diag_indices, theta, alpha,
							  N, L, N_h);
			} else if (system_simu == "network"){
				solve_network(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights%input_weights_mask, hidden_weights%hidden_weights_mask, output_weights%output_weights_mask, theta, alpha,
							  N, L);
			} else if (system_simu == "network_decoupled"){
				solve_network_decoupled(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights%input_weights_mask, hidden_weights%hidden_weights_mask, output_weights%output_weights_mask, N, L);
			} else {
				cout << system_simu << " is not a valid value for the system_simu option." << endl;
				abort();
			}
			clock_t end_solve = clock();
			cumulative_solve_time += (end_solve - start_solve) / (double)CLOCKS_PER_SEC;

			// compute deltas and gradient
			clock_t start_backprop = clock();
			mat f_prime_activations = f_prime_matrix(activations);
			if (grad_comp == "backprop_standard"){
				get_gradient_node_by_node(input_weight_gradient, weight_gradient, output_weight_gradient, input_data, node_states, f_prime_activations, g_primes, outputs, targets, hidden_weights%hidden_weights_mask, output_weights%output_weights_mask, diag_indices, N, L, theta, alpha);
			} else if (grad_comp == "backprop_classic"){
				get_gradient_classical_backprop(input_weight_gradient, weight_gradient, output_weight_gradient, input_data, node_states, f_prime_activations, g_primes, outputs, targets, hidden_weights%hidden_weights_mask, output_weights%output_weights_mask, diag_indices, N, L);
			} else {
				cout << grad_comp << " is not a valid value for the grad_comp option." << endl;
				abort();
			}
			clock_t end_backprop = clock();
			cumulative_backprop_time += (end_backprop - start_backprop) / (double)CLOCKS_PER_SEC;
			// eventually compute numerical gradient for comparision and save numerical gradient and backprop gradient in files
			if (gradient_check && step_index == 1){
				if (system_simu != "dde_ibp"){
					cout << "The option to compute the gradient numerically is only available if the method 'dde_ibp' is used to solve the delay system. " << endl;
					cout << "Set the option 'system_simu' to 'dde_ibp' or set the option 'num_gradient' to false." << endl;
					abort();
				}
				get_num_gradient(num_input_weight_gradient, num_weight_gradient, num_output_weight_gradient,
								 activations, node_states, output_activations, outputs, g_primes, label,
								 input_data, input_weights, hidden_weights, output_weights, diag_indices, theta, alpha, N, L, N_h, epsilon_gradient_check);
				if (print_gradients_to_file){
					write_gradients_to_file(input_weight_gradient, weight_gradient, output_weight_gradient, num_input_weight_gradient, num_weight_gradient, num_output_weight_gradient, epoch + 1, step_index, 0);
				}
				double cos_sim = cosine_similarity(input_weight_gradient, weight_gradient, output_weight_gradient, num_input_weight_gradient, num_weight_gradient, num_output_weight_gradient);
				similarity_vector.push_back(cos_sim);
			}

			input_weight_gradient = input_weight_gradient % input_weights_mask;
			weight_gradient = weight_gradient % hidden_weights_mask;
			output_weight_gradient = output_weight_gradient % output_weights_mask;

			// perform weight updates
			input_weights += - eta * input_weight_gradient;
			hidden_weights += - eta * weight_gradient;
			output_weights += - eta * output_weight_gradient;

			// record data for video: step 1 to end
			if (record_data && find(step_indices_for_video.begin(), step_indices_for_video.end(), epoch * 50000 + step_index) != step_indices_for_video.end()){
				int rec_step = epoch * 50000 + step_index;
				cout << "record data after step: " << rec_step << endl;
				//save weights
				string input_file_name = "video/weights_input_step_" + to_string(rec_step) + ".txt";
				string hidden_file_name = "video/weights_hidden_step_" + to_string(rec_step) + ".txt";
				string output_file_name = "video/weights_output_step_" + to_string(rec_step) + ".txt";
				input_weights.save(input_file_name, csv_ascii);
				hidden_weights.save(hidden_file_name, raw_ascii);
				output_weights.save(output_file_name, csv_ascii);
				//run system with examples from validation batch and save input, output and x_states
				for (int example_index = 0; example_index < number_of_examples; ++example_index){
					// select image as input
					vec input_data = test_images.row(example_index).t();
					int label = test_labels[index];
					// run system and record data for x_states
					solve_dde_ibp(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights, hidden_weights, output_weights, diag_indices, theta, alpha,
							  N, L, N_h, example_index+1, rec_step);
					// save input and output
					string input_file_name = "video/vector_input_step_" + to_string(rec_step) + "_example_" + to_string(example_index + 1) + ".txt";
					input_data.save(input_file_name, csv_ascii);
					string output_file_name = "video/vector_output_step_" + to_string(rec_step) + "_example_" + to_string(example_index + 1) + ".txt";
					ofstream output_file;
					output_file.open(output_file_name);
					for (double y_p : outputs){
						output_file << y_p << endl;
					}
					output_file.close();
				}
			}
		}

		//weight scaling
		input_weights_scaled = input_weights / (1.0 - dropout_rate);
		hidden_weights_scaled = hidden_weights / (1.0 - dropout_rate);
		output_weights_scaled = output_weights / (1.0 - dropout_rate);

		// loop to get accuracy on training set:
		double mse_sum = 0.0; ;
		for (int index = 0; index < number_of_training_batches *training_batch_size; ++index){
			//cout << "validation step (on training set)" << index + 1 << endl;

			// select image as input
			div_t div_result = div(index, training_batch_size);
			int batch_index = training_batch_indices[div_result.quot];
			int image_index = div_result.rem;
			vec input_data = train_images.tube(batch_index, image_index);
			int label = train_labels[batch_index][image_index];

			// get target
			double targets[P];
			for (int p = 0; p < P; ++p){
				targets[p] = input_data(p);
			}

			// add noise to input:
			vec noise = noise_sigma * vec(M, fill::randn);
			input_data += noise;
			for (int m = 0; m < M; ++m){
				if (input_data[m] < 0.0){
					input_data[m] = 0.0;
				}
				if (input_data[m] > 1.0){
					input_data[m] = 1.0;
				}
			}

			// solve the DDE (or network)
			clock_t start_solve = clock();
			if (system_simu == "dde_ibp"){
				solve_dde_ibp(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, diag_indices, theta, alpha,
							  N, L, N_h);
			} else if (system_simu == "dde_heun"){
				solve_dde_heun(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, diag_indices, theta, alpha,
							  N, L, N_h);
			} else if (system_simu == "network"){
				solve_network(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, theta, alpha,
							  N, L);
			} else if (system_simu == "network_decoupled"){
				solve_network_decoupled(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, N, L);
			} else {
				cout << system_simu << " is not a valid value for the system_simu parameter." << endl;
				abort();
			}
			clock_t end_solve = clock();
			cumulative_solve_time += (end_solve - start_solve) / (double)CLOCKS_PER_SEC;

			double mse = 0;
			for (int p = 0; p < P; ++p){
				mse += pow(outputs[p] - targets[p], 2.0);
			}
			mse = mse/(double)P;
			mse_sum += mse;
		}
		training_accuracy_vector.push_back(mse_sum / 50000.0);

		// loop for validation:
		mse_sum = 0.0; ;
		for (int index = 0; index < test_batch_size; ++index){

			vec input_data = test_images.row(index).t();
			int label = test_labels[index];

			// get target
			double targets[P];
			for (int p = 0; p < P; ++p){
				targets[p] = input_data(p);
			}

			// add noise to input:
			vec noise = noise_sigma * vec(M, fill::randn);
			input_data += noise;
			for (int m = 0; m < M; ++m){
				if (input_data[m] < 0.0){
					input_data[m] = 0.0;
				}
				if (input_data[m] > 1.0){
					input_data[m] = 1.0;
				}
			}

			// solve the DDE (or network)
			clock_t start_solve = clock();
			if (system_simu == "dde_ibp"){
				solve_dde_ibp(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, diag_indices, theta, alpha,
							  N, L, N_h);
			} else if (system_simu == "dde_heun"){
				solve_dde_heun(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, diag_indices, theta, alpha,
							  N, L, N_h);
			} else if (system_simu == "network"){
				solve_network(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, theta, alpha,
							  N, L);
			} else if (system_simu == "network_decoupled"){
				solve_network_decoupled(activations, node_states, output_activations, outputs, g_primes,
							  input_data, input_weights_scaled, hidden_weights_scaled, output_weights_scaled, N, L);
			} else {
				cout << system_simu << " is not a valid value for the system_simu parameter." << endl;
				abort();
			}
			clock_t end_solve = clock();
			cumulative_solve_time += (end_solve - start_solve) / (double)CLOCKS_PER_SEC;

			// save example images
			if (index < save_examples){
				string original_file_name = "example_results/epoch_" + to_string(epoch + 1) + "_example_" + to_string(index + 1) + "_original.txt";
				string input_file_name = "example_results/epoch_" + to_string(epoch + 1) + "_example_" + to_string(index + 1) + "_input.txt";
				string output_file_name = "example_results/epoch_" + to_string(epoch + 1) + "_example_" + to_string(index + 1) + "_output.txt";
				ofstream original_file;
				original_file.open(original_file_name);
				for (int p = 0; p < P; ++p){
					original_file << targets[p] << endl;
				}
				original_file.close();
				ofstream input_file;
				input_file.open(input_file_name);
				for (int p = 0; p < P; ++p){
					input_file << input_data(p) << endl;
				}
				input_file.close();
				ofstream output_file;
				output_file.open(output_file_name);
				for (int p = 0; p < P; ++p){
					output_file << outputs[p] << endl;
				}
				output_file.close();
			}

			double mse = 0;
			for (int p = 0; p < P; ++p){
				mse += pow(outputs[p] - targets[p], 2.0);
			}
			mse = mse/(double)P;
			mse_sum += mse;
		}
		accuracy_vector.push_back(mse_sum / 10000.0);

		cout << "epoch " << epoch + 1 << ": validation MSE = " << mse_sum / 10000.0 << endl;
		// eventually print weights to file at end of each epoch
		if (print_weights_to_file){
			print_weights(input_weights_scaled, hidden_weights_scaled, output_weights_scaled, diag_indices, 0, epoch);
		}

		ende = clock();
		epoch_time = ((double) (ende - start)) / (double)CLOCKS_PER_SEC;
		time_vector.push_back(epoch_time);
	}

	// print result to file:
	print_results(results_file_name, 0, diag_indices, training_accuracy_vector, accuracy_vector, similarity_vector);

	// for cpu time measuring
	clock_t end_overall = clock();
	double cpu_time_overall = (end_overall - start_overall) / (double)CLOCKS_PER_SEC;
	double cpu_time_residual = cpu_time_overall - cumulative_solve_time - cumulative_backprop_time;
	double cpu_time_solve_percentage = 100.0 * cumulative_solve_time / cpu_time_overall;
	double cpu_time_backprop_percentage = 100.0 * cumulative_backprop_time / cpu_time_overall;
	double cpu_time_residual_percentage = 100.0 * cpu_time_residual / cpu_time_overall;

	// get current time and date and print to results text file
	auto current_clock = chrono::system_clock::now();
	time_t current_time = chrono::system_clock::to_time_t(current_clock);
	ofstream results_file;
	results_file.open(results_file_name, ios_base::app);
	results_file << endl;
	results_file << endl;
	results_file << "total cpu time (in seconds): " << cpu_time_overall << endl;
	results_file << "cumulative cpu time for solving the DDE or network (in seconds): " << cumulative_solve_time << " (" << cpu_time_solve_percentage << "%)" << endl;
	results_file << "cumulative cpu time for backpropagation (in seconds): " << cumulative_backprop_time << " (" << cpu_time_backprop_percentage << "%)" << endl;
	results_file << "residual cpu time (in seconds): " << cpu_time_residual << " (" << cpu_time_residual_percentage << "%)" << endl;
	results_file << endl;
	results_file << "end of simulation: " << ctime(&current_time);
	results_file.close();

	return 0;
}

ext: cpp


text:
#include "layer_serializer.h"
#include "vi/nn/activation_function.h"
#include "vi/nn/layer.h"
#include <boost/property_tree/ptree.hpp>
#include <sstream>
#include <typeinfo>

namespace pt = boost::property_tree;

namespace vi {
namespace io {
layer_serializer::layer_serializer(const vi::nn::layer& layer) : layer_(layer) {}

void layer_serializer::serialize(boost::property_tree::ptree& layer_node) {
  const vi::nn::activation_function* a = layer_.activation().get();
  std::string activation_name("");
  if (typeid(*a) == typeid(vi::nn::sigmoid_activation)) {
    activation_name = "sigmoid";
  } else if (typeid(*a) == typeid(vi::nn::softmax_activation)) {
    activation_name = "softmax";
  } else if (typeid(*a) == typeid(vi::nn::hyperbolic_tangent)) {
    activation_name = "tanh";
  } else if (typeid(*a) == typeid(vi::nn::linear_activation)) {
    activation_name = "linear";
  } else {
    std::stringstream description;
    description << "invalid activation type: '" << typeid(*a).name() << "'";
    throw serializer::exception(description.str());
  }

  layer_node.put("activation_function", activation_name);
}
}
}

max_stars_repo_name: ville-k/vinn

max_issues_repo_name: ville-k/vinn

max_forks_repo_name: ville-k/vinn


text:
/*
 * Copyright 2014-present Facebook, Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#include <gtest/gtest.h>
#include <thrift/lib/cpp2/async/RequestChannel.h>
#include <thrift/lib/cpp2/async/FutureRequest.h>
#include <thrift/lib/cpp2/test/gen-cpp2/TestService.h>
#include <thrift/lib/cpp2/server/ThriftServer.h>
#include <thrift/lib/cpp2/async/HeaderClientChannel.h>
#include <thrift/lib/cpp2/test/gen-cpp2/DuplexService.h>
#include <thrift/lib/cpp2/test/gen-cpp2/DuplexClient.h>
#include <thrift/lib/cpp2/async/DuplexChannel.h>
#include <thrift/lib/cpp2/test/util/TestThriftServerFactory.h>

#include <thrift/lib/cpp/util/ScopedServerThread.h>
#include <folly/io/async/EventBase.h>
#include <thrift/lib/cpp/async/TAsyncSocket.h>

#include <thrift/lib/cpp2/async/GssSaslClient.h>

#include <folly/init/Init.h>

#include <boost/cast.hpp>
#include <boost/lexical_cast.hpp>
#include <memory>
#include <atomic>

using namespace apache::thrift;
using namespace apache::thrift::test::cpp2;
using namespace apache::thrift::util;
using namespace apache::thrift::async;
using namespace folly;
using std::shared_ptr;
using std::unique_ptr;
using std::make_shared;

class TestServiceInterface : public TestServiceSvIf {
public:
  void async_tm_sendResponse(
      unique_ptr<HandlerCallback<unique_ptr<std::string>>> callback,
      int64_t size) override {
    const auto& headers = callback->getConnectionContext()->
      getHeader()->getHeaders();
    if (headers.find("security_test")->second == "1") {
      EXPECT_NE(callback->getConnectionContext()->
                  getSaslServer()->getClientIdentity(), "");
    }
    callback.release()->resultInThread(folly::to<std::string>(size));
  }
};

void enableSecurity(HeaderClientChannel* channel,
                    const apache::thrift::SecurityMech mech,
                    bool failSecurity = false) {
  char hostname[256];
  EXPECT_EQ(gethostname(hostname, 255), 0);

  // Construct the service identity prefix from hostname.
  // Service identity should be sys.<hostname_scheme_prefix>@hostname.
  std::string hostnameStr(hostname);
  size_t end = hostnameStr.find_first_not_of("abcdefghijklmnopqrstuvwxyz");
  std::string serviceIdentityPrefix = std::string("sys.") +
    hostnameStr.substr(0, end) + std::string("@");

  std::string clientIdentity = std::string("host/") + hostname;
  std::string serviceIdentity = serviceIdentityPrefix + hostname;
  if (failSecurity) {
    serviceIdentity = "bogus";
  }

  channel->setSecurityPolicy(THRIFT_SECURITY_REQUIRED);

  auto saslClient = std::make_unique<GssSaslClient>(channel->getEventBase());
  saslClient->setClientIdentity(clientIdentity);
  saslClient->setServiceIdentity(serviceIdentity);
  saslClient->setSaslThreadManager(make_shared<SaslThreadManager>(
      make_shared<SecurityLogger>()));
  saslClient->setCredentialsCacheManager(
      make_shared<krb5::Krb5CredentialsCacheManager>());
  saslClient->setSecurityMech(mech);
  channel->setSaslClient(std::move(saslClient));
  channel->setSaslTimeout(5000);
}

HeaderClientChannel::Ptr getClientChannel(EventBase* eb,
                                          const folly::SocketAddress& address,
                                          bool failSecurity = false) {
  auto socket = TAsyncSocket::newSocket(eb, address);
  auto channel = HeaderClientChannel::newChannel(socket);

  enableSecurity(channel.get(), apache::thrift::SecurityMech::KRB5_GSS,
                 failSecurity);

  return channel;
}

class Countdown {
public:
  Countdown(int count, std::function<void()> f)
    : count_(count), f_(f)
  {}
  void down() {
    if (--count_ == 0) {
      f_();
    }
  }
private:
  int count_;
  std::function<void()> f_;
};

void runTest(std::function<void(HeaderClientChannel* channel)> setup) {
  apache::thrift::TestThriftServerFactory<TestServiceInterface> factory;
  factory.useStubSaslServer(false);
  ScopedServerThread sst(factory.create());

  EventBase base;
  auto channel = getClientChannel(&base, *sst.getAddress());
  setup(channel.get());
  auto sp = channel->getSecurityPolicy();
  TestServiceAsyncClient client(std::move(channel));
  Countdown c(3, [&base](){base.terminateLoopSoon();});

  RpcOptions rpcOptions;
  rpcOptions.setWriteHeader("security_test",
                            sp == THRIFT_SECURITY_REQUIRED ? "1" : "0");
  client.sendResponse(rpcOptions, std::make_unique<FunctionReplyCallback>(
                      [&c,&sp](ClientReceiveState&& state) {
    EXPECT_FALSE(state.isException());
    if (sp == THRIFT_SECURITY_REQUIRED) {
      EXPECT_TRUE(state.isSecurityActive());
    } else {
      EXPECT_FALSE(state.isSecurityActive());
    }
    std::string res;
    try {
      TestServiceAsyncClient::recv_sendResponse(res, state);
    } catch(const std::exception&) {
      ADD_FAILURE();
    }
    EXPECT_EQ(res, "10");
    c.down();
  }), 10);

  // fail on time out
  base.tryRunAfterDelay([] { ADD_FAILURE(); }, 5000);

  base.tryRunAfterDelay([&client,&c,&sp] {
    RpcOptions rpcOptions1;
    rpcOptions1.setWriteHeader("security_test",
                              sp == THRIFT_SECURITY_REQUIRED ? "1" : "0");
    client.sendResponse(rpcOptions1, std::make_unique<FunctionReplyCallback>(
          [&c,&sp](ClientReceiveState&& state) {
      EXPECT_FALSE(state.isException());
      if (sp == THRIFT_SECURITY_REQUIRED) {
        EXPECT_TRUE(state.isSecurityActive());
      } else {
        EXPECT_FALSE(state.isSecurityActive());
      }
      std::string res;
      try {
        TestServiceAsyncClient::recv_sendResponse(res, state);
      } catch(const std::exception&) {
        ADD_FAILURE();
      }
      EXPECT_EQ(res, "10");
      c.down();
    }), 10);

    RpcOptions rpcOptions2;
    rpcOptions2.setWriteHeader("security_test",
                              sp == THRIFT_SECURITY_REQUIRED ? "1" : "0");
    client.sendResponse(rpcOptions2, std::make_unique<FunctionReplyCallback>(
          [&c,&sp](ClientReceiveState&& state) {
      EXPECT_FALSE(state.isException());
      if (sp == THRIFT_SECURITY_REQUIRED) {
        EXPECT_TRUE(state.isSecurityActive());
      } else {
        EXPECT_FALSE(state.isSecurityActive());
      }
      std::string res;
      try {
        TestServiceAsyncClient::recv_sendResponse(res, state);
      } catch(const std::exception&) {
        ADD_FAILURE();
      }
      EXPECT_EQ(res, "10");
      c.down();
    }), 10);
  }, 1);

  base.loopForever();
}

TEST(Security, Basic) {
  runTest([](HeaderClientChannel*) {});
}

TEST(Security, CompressionZlib) {
  runTest([](HeaderClientChannel* channel) {
    channel->setTransform(transport::THeader::ZLIB_TRANSFORM);
  });
}

TEST(Security, CompressionSnappy) {
  runTest([](HeaderClientChannel* channel) {
    channel->setTransform(transport::THeader::SNAPPY_TRANSFORM);
  });
}

TEST(Security, CompressionZstd) {
  runTest([](HeaderClientChannel* channel) {
    channel->setTransform(transport::THeader::ZSTD_TRANSFORM);
  });
}

TEST(Security, ProtocolBinary) {
  runTest([](HeaderClientChannel* channel) {
    channel->setProtocolId(protocol::T_BINARY_PROTOCOL);
  });
}

TEST(Security, ProtocolCompact) {
  runTest([](HeaderClientChannel* channel) {
    channel->setProtocolId(protocol::T_COMPACT_PROTOCOL);
  });
}

TEST(Security, SASL) {
  runTest([](HeaderClientChannel* channel) {
    channel->getSaslClient()->setSecurityMech(
      apache::thrift::SecurityMech::KRB5_SASL);
  });
}

TEST(Security, GSS_NO_MUTUAL) {
  runTest([](HeaderClientChannel* channel) {
    channel->getSaslClient()->setSecurityMech(
      apache::thrift::SecurityMech::KRB5_GSS_NO_MUTUAL);
  });
}

TEST(Security, GSS) {
  runTest([](HeaderClientChannel* channel) {
    channel->getSaslClient()->setSecurityMech(
      apache::thrift::SecurityMech::KRB5_GSS);
  });
}

TEST(Security, Fallback) {
  runTest([](HeaderClientChannel* channel) {
    channel->setSecurityPolicy(THRIFT_SECURITY_PERMITTED);
    channel->setSaslTimeout(1);
    channel->getSaslClientCallback()->setSendServerHook([]{
      /* sleep override */
      usleep(20000);
    });
  });
}

class DuplexClientInterface : public DuplexClientSvIf {
public:
  DuplexClientInterface(int32_t first, int32_t count, bool& success)
    : expectIndex_(first), lastIndex_(first + count), success_(success)
  {}

  void async_tm_update(unique_ptr<HandlerCallback<int32_t>> callback,
                       int32_t currentIndex) override {
    auto callbackp = callback.release();
    EXPECT_EQ(currentIndex, expectIndex_);
    expectIndex_++;
    EventBase *eb = callbackp->getEventBase();
    callbackp->resultInThread(currentIndex);
    if (expectIndex_ == lastIndex_) {
      success_ = true;
      eb->runInEventBaseThread([eb] { eb->terminateLoopSoon(); });
    }
  }
private:
  int32_t expectIndex_;
  int32_t lastIndex_;
  bool& success_;
};

class Updater {
public:
  Updater(shared_ptr<DuplexClientAsyncClient> client,
          EventBase* eb,
          int32_t startIndex,
          int32_t numUpdates,
          int32_t interval)
    : client_(client)
    , eb_(eb)
    , startIndex_(startIndex)
    , numUpdates_(numUpdates)
    , interval_(interval)
  {}

  void update() {
    int32_t si = startIndex_;
    client_->update([si](ClientReceiveState&& state) {
      EXPECT_FALSE(state.isException());
      EXPECT_TRUE(state.isSecurityActive());
      try {
        int32_t res = DuplexClientAsyncClient::recv_update(state);
        EXPECT_EQ(res, si);
      } catch (const std::exception&) {
        ADD_FAILURE();
      }
    }, startIndex_);
    startIndex_++;
    numUpdates_--;
    if (numUpdates_ > 0) {
      Updater updater(*this);
      eb_->tryRunAfterDelay([updater]() mutable {
        updater.update();
      }, interval_);
    }
  }
private:
  shared_ptr<DuplexClientAsyncClient> client_;
  EventBase* eb_;
  int32_t startIndex_;
  int32_t numUpdates_;
  int32_t interval_;
};

class DuplexServiceInterface : public DuplexServiceSvIf {
  void async_tm_registerForUpdates(unique_ptr<HandlerCallback<bool>> callback,
                                   int32_t startIndex,
                                   int32_t numUpdates,
                                   int32_t interval) override {

    EXPECT_NE(callback->getConnectionContext()->
                getSaslServer()->getClientIdentity(), "");

    auto callbackp = callback.release();
    auto ctx = callbackp->getConnectionContext()->getConnectionContext();
    CHECK(ctx != nullptr);
    auto client = ctx->getDuplexClient<DuplexClientAsyncClient>();
    auto eb = callbackp->getEventBase();
    CHECK(eb != nullptr);
    if (numUpdates > 0) {
      Updater updater(client, eb, startIndex, numUpdates, interval);
      eb->runInEventBaseThread([updater]() mutable {updater.update();});
    };
    callbackp->resultInThread(true);
  }

  void async_tm_regularMethod(unique_ptr<HandlerCallback<int32_t>> callback,
                              int32_t val) override {
    EXPECT_NE(callback->getConnectionContext()->
                getSaslServer()->getClientIdentity(), "");
    callback.release()->resultInThread(val * 2);
  }
};

void duplexTest(const apache::thrift::SecurityMech mech) {
  enum {START=1, COUNT=3, INTERVAL=1};
  apache::thrift::TestThriftServerFactory<DuplexServiceInterface> factory;
  factory.useStubSaslServer(false).duplex(true);
  ScopedServerThread duplexsst(factory.create());
  EventBase base;
  std::shared_ptr<TAsyncSocket> socket(
    TAsyncSocket::newSocket(&base, *duplexsst.getAddress()));

  auto duplexChannel =
      std::make_shared<DuplexChannel>(DuplexChannel::Who::CLIENT, socket);
  enableSecurity(duplexChannel->getClientChannel().get(), mech);
  DuplexServiceAsyncClient client(duplexChannel->getClientChannel());

  bool success = false;
  ThriftServer clients_server(duplexChannel->getServerChannel());
  clients_server.setInterface(std::make_shared<DuplexClientInterface>(
      START, COUNT, success));
  clients_server.serve();

  client.registerForUpdates([](ClientReceiveState&& state) {
    EXPECT_FALSE(state.isException());
    EXPECT_TRUE(state.isSecurityActive());
    try {
      bool res = DuplexServiceAsyncClient::recv_registerForUpdates(state);
      EXPECT_TRUE(res);
    } catch (const std::exception&) {
      ADD_FAILURE();
    }
  }, START, COUNT, INTERVAL);

  // fail on time out
  base.tryRunAfterDelay([] { ADD_FAILURE(); }, 5000);

  base.loopForever();

  EXPECT_TRUE(success);
}

TEST(Security, DuplexSASL) {
  duplexTest(apache::thrift::SecurityMech::KRB5_SASL);
}

TEST(Security, DuplexGSS) {
  duplexTest(apache::thrift::SecurityMech::KRB5_GSS);
}

TEST(Security, DuplexGSSNoMutual) {
 duplexTest(apache::thrift::SecurityMech::KRB5_GSS_NO_MUTUAL);
}

// Test if multiple requests are pending in a queue, for security to establish,
// then we flow RequestContext correctly with each request.
void runRequestContextTest(bool failSecurity) {
  apache::thrift::TestThriftServerFactory<TestServiceInterface> factory;
  factory.useStubSaslServer(false);
  ScopedServerThread sst(factory.create());
  EventBase base;
  auto channel = getClientChannel(&base, *sst.getAddress(), failSecurity);
  TestServiceAsyncClient client(std::move(channel));
  Countdown c(2, [&base](){base.terminateLoopSoon();});

  {
    // Send first request with a unique RequestContext. This would trigger
    // security. Rest of the request would queue behind it.
    folly::RequestContextScopeGuard rctx;
    folly::RequestContext::get()->setContextData("first", nullptr);
    client.sendResponse(
        [&c](ClientReceiveState&&) {
          EXPECT_TRUE(folly::RequestContext::get()->hasContextData("first"));
          c.down();
        },
        10);
  }

  {
    // Send another request with a unique RequestContext. This request would
    // queue behind the first one inside HeaderClientChannel.
    folly::RequestContextScopeGuard rctx;
    folly::RequestContext::get()->setContextData("second", nullptr);
    client.sendResponse(
        [&c](ClientReceiveState&&) {
          EXPECT_FALSE(folly::RequestContext::get()->hasContextData("first"));
          EXPECT_TRUE(folly::RequestContext::get()->hasContextData("second"));
          c.down();
        },
        10);
  }

  // Now start looping the eventbase to guarantee that all the above requests
  // would always queue.
  base.loopForever();
}

TEST(SecurityRequestContext, Success) {
  runRequestContextTest(false);
}

TEST(SecurityRequestContext, Fail) {
  runRequestContextTest(true);
}

int main(int argc, char** argv) {
  setenv("KRB5_CONFIG", "/etc/krb5-thrift.conf", 0);
  ::testing::InitGoogleTest(&argc, argv);
  folly::init(&argc, &argv);
  return RUN_ALL_TESTS();
}

max_stars_repo_name: lucyge/FBThrift

max_issues_repo_name: lucyge/FBThrift

max_forks_repo_name: lucyge/FBThrift


text:
// part of stage 1 of gazebo simulation
// Reading the marker poses from a text file(marker_position.txt) and saving it as an ArucoMessage in a bagfile (aruco.bag)
#include <vector>
#include <string>
#include <fstream>

#include <ros/ros.h>
#include <tf/transform_broadcaster.h>
#include <rosbag/bag.h>
#include <aruco_mapping.h>
#include <geometry_msgs/PoseWithCovarianceStamped.h>
#include <geometry_msgs/Pose.h>
#include <boost/lexical_cast.hpp>

int main(int argc, char** argv){

  ros::init(argc, argv, "saving_aruco_wrt_map_from_data_file");

  std::string bag_file_name = "/home/fauzan/bags/aruco_from_data.bag";
  std::string data_file_name = "marker_positions.txt";

  ros::NodeHandle node("~");

  node.getParam("bag_file_name", bag_file_name);
  node.getParam("data_file_name", data_file_name);

  std::ifstream data_file(data_file_name.c_str());

  aruco_mapping::ArucoMarker aruco_msg;
  std::vector<int> marker_ids_vector;

  int marker_id;
  double d_x, d_y, d_z;
  double r_r, r_p, r_y;
  double q_x, q_y, q_z, q_w;
  geometry_msgs::Pose marker_pose;

  /**
   * marker_positions.txt can have Euler or Quarternion angles.
   * Uncomment the required format
   */

  int count = 0;

  // Euler
  // while(data_file >> marker_id >> d_x >> d_y >> d_z >> r_r >> r_p >> r_y){

  // Quarterion
  while(data_file >> marker_id >> d_x >> d_y >> d_z >> q_x >> q_y >> q_z >> q_w){

    marker_pose.position.x = d_x;
    marker_pose.position.y = d_y;
    marker_pose.position.z = d_z;

    // Euler
    // marker_pose.orientation = tf::createQuaternionMsgFromRollPitchYaw(r_r, r_p, r_y);
    // std::cout << marker_id << "\t";
    // std::cout << marker_pose.orientation.x << "\t";
    // std::cout << marker_pose.orientation.y << "\t";
    // std::cout << marker_pose.orientation.z << "\t";
    // std::cout << marker_pose.orientation.w << std::endl;

    // Quarternion
    marker_pose.orientation.x = q_x;
    marker_pose.orientation.y = q_y;
    marker_pose.orientation.z = q_z;
    marker_pose.orientation.w = q_w;

    aruco_msg.global_marker_poses.push_back(marker_pose);
    marker_ids_vector.push_back(marker_id);

    count++;
  }

  std::cout << "Found " << count << " markers" << std::endl;

  aruco_msg.num_of_visible_markers = count;
  aruco_msg.marker_ids = marker_ids_vector;

  rosbag::Bag bag;
  bag.open(bag_file_name, rosbag::bagmode::Write);
  bag.write("map_aruco_pose", ros::Time::now(), aruco_msg);
  bag.close();

  std::cout << "Done" << std::endl;

  return 0;
}