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)
 * NN representation: sum over [a_i * (1 + e^(-x * w_i + b_i))^(-1)]
 * -------------------------------------------
 * Optimal NN setting with biases (2 inner neurons)
 * Path   1. w =     -1.66009975, b =     -0.40767447, a =      2.46457042
 * Path   2. w =     -4.38622765, b =      2.75707816, a =     -8.04752347
 * @author Michal Kravčenko
 * @date 17.7.18 -
 */

#include <random>
#include <iostream>
#include <chrono>
#include <4neuro.h>

void optimize_via_particle_swarm(l4n::DESolver& solver,
                                 l4n::MultiIndex& alpha,
                                 size_t max_iters,
                                 size_t n_particles) {

    printf("Solution via the particle swarm optimization!\n");
    std::vector<double> domain_bounds(
        2 * (solver.get_solution(alpha)->get_n_biases() + solver.get_solution(alpha)->get_n_weights()));

    for (size_t i = 0; i < domain_bounds.size() / 2; ++i) {
        domain_bounds[2 * i]     = -10;
        domain_bounds[2 * i + 1] = 10;
    }

    double c1 = 1.7;
    double c2 = 1.7;
    double w  = 0.700;

    /* if the maximal velocity from the previous step is less than 'gamma' times the current maximal velocity, then one
     * terminating criterion is met */
    double gamma = 0.5;

    /* if 'delta' times 'n' particles are in the centroid neighborhood given by the radius 'epsilon', then the second
     * terminating criterion is met ('n' is the total number of particles) */
    double epsilon = 0.02;
    double delta   = 0.7;

    l4n::ParticleSwarm swarm(
        &domain_bounds,
        c1,
        c2,
        w,
        gamma,
        epsilon,
        delta,
        n_particles,
        max_iters
    );

    solver.solve(swarm);

}

void optimize_via_gradient_descent(l4n::DESolver& solver,
                                   double accuracy) {
    printf("Solution via a gradient descent method!\n");
    solver.randomize_parameters();
// TODO does not work (poor design of netsum)
//    l4n::LevenbergMarquardt leven(10000, 0, 1e-6, 1e-6, 1e-6);
//    solver.solve(leven);


    l4n::GradientDescent gd(accuracy,
                            1000,
                            500000);
    solver.solve(gd);

}

void export_solution(size_t n_test_points,
                     double te,
                     double ts,
                     l4n::DESolver& solver,
                     l4n::MultiIndex& alpha_0,
                     l4n::MultiIndex& alpha_1,
                     l4n::MultiIndex& alpha_2,
                     const std::string prefix) {
    l4n::NeuralNetwork* solution    = solver.get_solution(alpha_0);
    l4n::NeuralNetwork* solution_d  = solver.get_solution(alpha_1);
    l4n::NeuralNetwork* solution_dd = solver.get_solution(alpha_2);

    /* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
    /* first boundary condition & its error */

    char buff[256];
    sprintf(buff,
            "%sdata_1d_ode1.txt",
            prefix.c_str());
    std::string final_fn(buff);

    std::ofstream ofs(final_fn,
                      std::ofstream::out);
    printf("Exporting files '%s': %7.3f%%\r",
           final_fn.c_str(),
           0.0);
    double frac = (te - ts) / (n_test_points - 1);

    std::vector<double> inp(1), out(1);

    for (size_t i = 0; i < n_test_points; ++i) {
        double x = frac * i + ts;
        inp[0] = x;

        solution->eval_single(inp,
                              out);
        double F = out[0];

        solution_d->eval_single(inp,
                                out);
        double DF = out[0];

        solution_dd->eval_single(inp,
                                 out);
        double DDF = out[0];

        ofs << i + 1 << " " << x << " " << std::pow(l4n::E,
                                                    -2 * x) * (3 * x + 1) << " " << F << " "
            << std::pow(l4n::E,
                        -2 * x) * (1 - 6 * x) << " " << DF << " " << 4 * std::pow(l4n::E,
                                                                                  -2 * x) * (3 * x - 2)
            << " " << DDF << std::endl;

        printf("Exporting files '%s': %7.3f%%\r",
               final_fn.c_str(),
               (100.0 * i) / (n_test_points - 1));
        std::cout.flush();
    }
    printf("Exporting files '%s': %7.3f%%\r",
           final_fn.c_str(),
           100.0);
    std::cout.flush();
    ofs.close();

    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;

}

void test_ode(double accuracy,
              size_t n_inner_neurons,
              size_t train_size,
              double ds,
              double de,
              size_t n_test_points,
              double ts,
              double te,
              size_t max_iters,
              size_t n_particles) {

    std::cout << "Finding a solution via the Particle Swarm Optimization and Gradient descent method!" << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;

    /* SOLVER SETUP */
    size_t        n_inputs    = 1;
    size_t        n_equations = 3;
    l4n::DESolver solver_01(n_equations,
                            n_inputs,
                            n_inner_neurons);

    /* SETUP OF THE EQUATIONS */
    l4n::MultiIndex alpha_0(n_inputs);
    l4n::MultiIndex alpha_1(n_inputs);
    l4n::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_0,
                                           "4.0");
    solver_01.add_to_differential_equation(0,
                                           alpha_1,
                                           "4.0");
    solver_01.add_to_differential_equation(0,
                                           alpha_2,
                                           "1.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");

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

    double d1_s = ds, d1_e = de, frac;

    /* TRAIN DATA FOR THE GOVERNING DE */
    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_g;
    std::vector<double>                                              test_points(train_size);


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

        test_points[i] = inp[0];
    }

    /* CHEBYSCHEV TRAIN SET */
    l4n::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));
    l4n::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));
    l4n::DataSet ds_02(&data_vec_dy);

    /* Placing the conditions into the solver */
    solver_01.set_error_function(0,
                                 l4n::ErrorFunctionType::ErrorFuncMSE,
                                 ds_00);

    solver_01.set_error_function(1,
                                 l4n::ErrorFunctionType::ErrorFuncMSE,
                                 ds_01);
    solver_01.set_error_function(2,
                                 l4n::ErrorFunctionType::ErrorFuncMSE,
                                 ds_02);


    /* TRAINING METHOD SETUP */
    /*  optimize_via_particle_swarm( solver_01, alpha_0, max_iters, n_particles );
      export_solution( n_test_points, te, ts, solver_01 , alpha_0, alpha_1, alpha_2, "particle_" );*/
    auto start = std::chrono::system_clock::now();

    optimize_via_gradient_descent(solver_01,
                                  accuracy);
    export_solution(n_test_points,
                    te,
                    ts,
                    solver_01,
                    alpha_0,
                    alpha_1,
                    alpha_2,
                    "gradient_");

    auto                          end             = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end - start;
    std::cout << "elapsed time: " << elapsed_seconds.count() << std::endl;
}

int main() {
    std::cout << "Running lib4neuro Ordinary Differential Equation example   1" << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;
    std::cout << "          Governing equation: y''(x) + 4y'(x) + 4y(x) = 0.0, for x in [0, 4]" << std::endl;
    std::cout << "Dirichlet boundary condition:                  y(0.0) = 1.0" << std::endl;
    std::cout << "  Neumann boundary condition:                 y'(0.0) = 1.0" << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;
    std::cout
        << "Expressing solution as y(x) = sum over [a_i / (1 + exp(bi - wxi*x ))], i in [1, n], where n is the number of hidden neurons"
        << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;

    unsigned int n_inner_neurons = 2;
    unsigned int train_size      = 10;
    double       accuracy        = 1e-1;
    double       ds              = 0.0;
    double       de              = 4.0;

    unsigned int test_size = 10;
    double       ts        = ds;
    double       te        = de + 2;

    size_t particle_swarm_max_iters = 10;
    size_t n_particles              = 2;

    test_ode(accuracy,
             n_inner_neurons,
             train_size,
             ds,
             de,
             test_size,
             ts,
             te,
             particle_swarm_max_iters,
             n_particles);


    return 0;
}