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");
    l4n::GradientDescent gd( accuracy, 1000 , 50);

    solver.randomize_parameters( );
    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, out;

    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 = {frac * i};
        out = {0.0};
        data_vec_g.emplace_back(std::make_pair(inp, out));

        test_points[i] = inp[0];
    }

    /* 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));
//
//        test_points[i] = inp[0];
//    }
   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-6;
    double ds = 0.0;
    double de = 4.0;

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

    size_t particle_swarm_max_iters = 1000;
    size_t n_particles = 100;

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


    return 0;
}