net_test_ode_1.cpp

/**
 * Example solving the following ODE:
 *
 * g(t) = (d^2/d^t)y(t) + 4 (d/dt)y(t) + 4y(t) = 0, for t in [0, 4]
 * y(0) = 1
 * (d/dt)y(0) = 1
 *
 * -------------------------------------------
 * Analytical solution: e^(-2x) * (3x + 1)
 *
 * @author Michal Kravčenko
 * @date 17.7.18 -
 */
#include <random>
#include <vector>
#include <utility>
#include <iostream>

#include "../include/4neuro.h"
#include "Solvers/DESolver.h"

void test_analytical_gradient_y(unsigned int n_inner_neurons) {
    unsigned int n_inputs = 1;
}

void test_analytical_y(unsigned int n_inner_neurons){
    unsigned int n_inputs = 1;
    unsigned int n_outputs = 1;
    unsigned int n_equations = 1;

    /* SETUP OF THE TRAINING DATA */
    std::vector<double> inp, out;

    double d1_s = 0.0, d1_e = 1.0, frac, alpha, x;

    /* TRAIN DATA FOR THE GOVERNING DE */
    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_g;
    unsigned int train_size = 150;

    /* ISOTROPIC TRAIN SET */
    frac = (d1_e - d1_s) / (train_size - 1);
//    for(unsigned int i = 0; i < train_size; ++i){
//        x = frac * i;
//        inp = {x};
//        out = {std::pow(E, -2 * x) * (3 * x + 1)};
//        data_vec_g.emplace_back(std::make_pair(inp, out));
//
//        std::cout << i + 1 << " " << x << " " << out[0] << std::endl;
//    }

    /* CHEBYSCHEV TRAIN SET */
    alpha = PI / (train_size - 1);
    frac = 0.5 * (d1_e - d1_s);
    for(unsigned int i = 0; i < train_size; ++i){
        x = (std::cos(PI - alpha * i) + 1.0) * frac + d1_s;
        inp = {x};
        out = {std::pow(E, -2 * x) * (3 * x + 1)};
        data_vec_g.emplace_back(std::make_pair(inp, out));

        std::cout << i + 1 << " " << x << " " << out[0] << std::endl;
    }
    DataSet ds_00(&data_vec_g);

//    return;

    DESolver solver_01( n_equations, n_inputs, n_inner_neurons, n_outputs );

    NeuralNetwork * solution = solver_01.get_solution();



    /* PARTICLE SWARM TRAINING METHOD SETUP */
    unsigned int max_iters = 150;

    //must encapsulate each of the partial error functions
    double *domain_bounds = new double[ 2 * n_inner_neurons * (n_inputs + n_outputs) ];
    for(unsigned int i = 0; i < n_inner_neurons * (n_inputs + n_outputs); ++i){
        domain_bounds[2 * i] = -800.0;
        domain_bounds[2 * i + 1] = 800.0;
    }

    double c1 = 0.5, c2 = 0.8, w = 0.8;

    unsigned int n_particles = 100;

    double gamma = 0.5, epsilon = 0.001, delta = 0.9;

    MSE mse(solution, &ds_00);
    ParticleSwarm swarm_01(&mse, domain_bounds, c1, c2, w, n_particles, max_iters);
    swarm_01.optimize( gamma, epsilon, delta );

    unsigned int n_test_points = 100;
    std::vector<double> input(1), output(1);
    double error = 0.0;
    for(unsigned int i = 0; i < n_test_points; ++i){
        input[0] = i * ((d1_e - d1_s) / (n_test_points - 1)) + d1_s;
        solution->eval_single(input, output);
        std::cout << i + 1 << " " << std::abs(output[0] - std::pow(E, -2 * input[0]) * (3 * input[0] + 1))<< " " << output[0] << std::endl;
    }

    delete [] domain_bounds;
}

void test_odr(unsigned int n_inner_neurons){
    unsigned int n_inputs = 1;
    unsigned int n_outputs = 1;
    unsigned int n_equations = 3;

    DESolver solver_01( n_equations, n_inputs, n_inner_neurons, n_outputs );

    /* SETUP OF THE EQUATIONS */
    MultiIndex alpha_0( n_inputs );
    MultiIndex alpha_1( n_inputs );
    MultiIndex alpha_2( n_inputs );
    alpha_2.set_partial_derivative(0, 2);
    alpha_1.set_partial_derivative(0, 1);

    /* the governing differential equation */
    solver_01.add_to_differential_equation( 0, alpha_2, 1.0 );
    solver_01.add_to_differential_equation( 0, alpha_1, 4.0 );
    solver_01.add_to_differential_equation( 0, alpha_0, 4.0 );

    /* dirichlet boundary condition */
    solver_01.add_to_differential_equation( 1, alpha_0, 1.0 );

    /* neumann boundary condition */
    solver_01.add_to_differential_equation( 2, alpha_1, 1.0 );

//    double weights[2] = {0.5, 1.0};
//    std::vector<double> inp = { 1.0 };
//    solver_01.eval_equation( 0, weights, inp );


    /* SETUP OF THE TRAINING DATA */
    std::vector<double> inp, out;

    double d1_s = 0.0, d1_e = 4.0, frac, alpha;

    /* TRAIN DATA FOR THE GOVERNING DE */
    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_g;
    unsigned int train_size = 100;

    /* ISOTROPIC TRAIN SET */
    frac = (d1_e - d1_s) / (train_size - 1);
    for(unsigned int i = 0; i < train_size; ++i){
        inp = {frac * i};
        out = {0.0};
        data_vec_g.emplace_back(std::make_pair(inp, out));
    }

    /* CHEBYSCHEV TRAIN SET */
//    alpha = PI / (train_size - 1);
//    frac = 0.5 * (d1_e - d1_s);
//    for(unsigned int i = 0; i < train_size; ++i){
//        inp = {(std::cos(alpha * i) + 1.0) * frac + d1_s};
//        out = {0.0};
//        data_vec_g.emplace_back(std::make_pair(inp, out));
//    }
    DataSet ds_00(&data_vec_g);

    /* TRAIN DATA FOR DIRICHLET BC */
    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_y;
    inp = {0.0};
    out = {1.0};
    data_vec_y.emplace_back(std::make_pair(inp, out));
    DataSet ds_01(&data_vec_y);

    /* TRAIN DATA FOR NEUMANN BC */
    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_dy;
    inp = {0.0};
    out = {1.0};
    data_vec_dy.emplace_back(std::make_pair(inp, out));
    DataSet ds_02(&data_vec_dy);

    /* Placing the conditions into the solver */
    solver_01.set_error_function( 0, ErrorFunctionType::ErrorFuncMSE, &ds_00 );
    solver_01.set_error_function( 1, ErrorFunctionType::ErrorFuncMSE, &ds_01 );
    solver_01.set_error_function( 2, ErrorFunctionType::ErrorFuncMSE, &ds_02 );





    /* PARTICLE SWARM TRAINING METHOD SETUP */
    unsigned int max_iters = 100;

    //must encapsulate each of the partial error functions
    double *domain_bounds = new double[ 2 * n_inner_neurons * (n_inputs + n_outputs) ];
    for(unsigned int i = 0; i < n_inner_neurons * (n_inputs + n_outputs); ++i){
        domain_bounds[2 * i] = -800.0;
        domain_bounds[2 * i + 1] = 800.0;
    }

    double c1 = 0.5, c2 = 1.5, w = 0.8;

    unsigned int n_particles = 100;

    double gamma = 0.5, epsilon = 0.02, delta = 0.9;

    solver_01.solve_via_particle_swarm( domain_bounds, c1, c2, w, n_particles, max_iters, gamma, epsilon, delta );

    NeuralNetwork *solution = solver_01.get_solution();


    unsigned int n_test_points = 100;
    std::vector<double> input(1), output(1);
    for(unsigned int i = 0; i < n_test_points; ++i){

        input[0] = i * ((d1_e - d1_s) / (n_test_points - 1)) + d1_s;

        solution->eval_single(input, output);

        std::cout << i + 1 << " " << output[0] << std::endl;
    }


    delete [] domain_bounds;
}

int main() {

    unsigned int n_inner_neurons = 20;

//    test_odr(n_inner_neurons);

    test_analytical_y(n_inner_neurons);
    test_analytical_gradient_y(n_inner_neurons);


    return 0;
}