net_test_pde_1.cpp

/**
 * Example solving the time dependent water flow 1D diffusion PDE:
 *
 * y_xx - y_t = 0, for (x, t) in [0, 1] x [0, 1]
 * y(0, t) = sin(t)
 * y(x, 0) = e^(-sqrt(0.5)x) * sin(-sqrt(0.5)x)
 *
 * -------------------------------------------
 * Analytical solution:
 * NN representation: sum over [a_i * (1 + e^(bi - x * w_ix - t * w_it))^(-1)]
 * -------------------------------------------
 * Optimal NN setting with biases (4 inner neurons)
 * Path   1. wx =      0.51954589, wt =     -0.48780445, b =      0.35656955, a =      1.69279158
 * Path   2. wx =     -1.24173503, wt =      1.13351300, b =      0.32528567, a =      1.69148458
 * Path   3. wx =      0.64754127, wt =      0.95758760, b =     -0.95852707, a =      2.77877453
 * Path   4. wx =      1.65439557, wt =     -0.31784248, b =     -1.81237586, a =     -3.96157108
 * @author Michal Kravčenko
 * @date 9.8.18
 */

#include <random>
#include <iostream>
#include <fstream>

#include <4neuro_public.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);

    solver.randomize_parameters();
    solver.solve(gd);
}

void export_solution(size_t n_test_points,
                     double te,
                     double ts,
                     l4n::DESolver& solver,
                     l4n::MultiIndex& alpha_00,
                     l4n::MultiIndex& alpha_01,
                     l4n::MultiIndex& alpha_20,
                     const std::string prefix) {
    l4n::NeuralNetwork* solution    = solver.get_solution(alpha_00);
    l4n::NeuralNetwork* solution_t  = solver.get_solution(alpha_01);
    l4n::NeuralNetwork* solution_xx = solver.get_solution(alpha_20);

    size_t i, j;
    double x, t;
    /* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
    /* first boundary condition & its error */

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

    printf("Exporting file '%s' : %7.3f%%\r",
           final_fn.c_str(),
           0.0);
    std::cout.flush();

    std::vector<double> input(2), output(1), output_t(1), output_xx(1);
    std::ofstream       ofs(final_fn,
                            std::ofstream::out);
    double              frac        = (te - ts) / (n_test_points - 1);
    for (i = 0; i < n_test_points; ++i) {
        x      = i * frac + ts;
        for (j = 0; j < n_test_points; ++j) {
            t     = j * frac + ts;
            input = {x, t};

            solution->eval_single(input,
                                  output);

            ofs << x << " " << t << " " << output[0] << std::endl;
            printf("Exporting file '%s' : %7.3f%%\r",
                   final_fn.c_str(),
                   (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
            std::cout.flush();
        }
    }
    printf("Exporting file '%s' : %7.3f%%\n",
           final_fn.c_str(),
           100.0);
    std::cout.flush();
    ofs.close();

    /* governing equation error */
    sprintf(buff,
            "%sdata_2d_pde1_first_equation_error.txt",
            prefix.c_str());
    final_fn = std::string(buff);

    ofs = std::ofstream(final_fn,
                        std::ofstream::out);
    printf("Exporting file '%s' : %7.3f%%\r",
           final_fn.c_str(),
           0.0);
    for (i = 0; i < n_test_points; ++i) {
        x      = i * frac + ts;
        for (j = 0; j < n_test_points; ++j) {
            t     = j * frac + ts;
            input = {x, t};

            solution_t->eval_single(input,
                                    output_t);
            solution_xx->eval_single(input,
                                     output_xx);

            ofs << x << " " << t << " " << std::fabs(output_xx[0] - output_t[0]) << std::endl;
            printf("Exporting file 'data_2d_pde1_first_equation_error.txt' : %7.3f%%\r",
                   (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
            std::cout.flush();
        }
    }
    printf("Exporting file '%s' : %7.3f%%\n",
           final_fn.c_str(),
           100.0);
    std::cout.flush();
    ofs.close();

    /* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
    /* first boundary condition & its error */
    sprintf(buff,
            "%sdata_1d_pde1_yt.txt",
            prefix.c_str());
    std::string final_fn_t(buff);

    sprintf(buff,
            "%sdata_1d_pde1_yx.txt",
            prefix.c_str());
    std::string final_fn_x(buff);

    ofs = std::ofstream(final_fn_t,
                        std::ofstream::out);
    std::ofstream ofs2(final_fn_x,
                       std::ofstream::out);
    printf("Exporting files '%s' and '%s' : %7.3f%%\r",
           final_fn_t.c_str(),
           final_fn_x.c_str(),
           0.0);
    for (i = 0; i < n_test_points; ++i) {
        x = frac * i + ts;
        t = frac * i + ts;

        double yt = std::sin(t);
        double yx = std::pow(l4n::E,
                             -0.707106781 * x) * std::sin(-0.707106781 * x);

        input = {0, t};
        solution->eval_single(input,
                              output);
        double evalt = output[0];

        input = {x, 0};
        solution->eval_single(input,
                              output);
        double evalx = output[0];

        ofs << i + 1 << " " << t << " " << yt << " " << evalt << " " << std::fabs(evalt - yt) << std::endl;
        ofs2 << i + 1 << " " << x << " " << yx << " " << evalx << " " << std::fabs(evalx - yx) << std::endl;

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

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

void test_pde(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) {

    /* do not change below */
    size_t        n_inputs    = 2;
    size_t        n_equations = 3;
    l4n::DESolver solver_01(n_equations,
                            n_inputs,
                            n_inner_neurons);

    /* SETUP OF THE EQUATIONS */
    l4n::MultiIndex alpha_00(n_inputs);
    l4n::MultiIndex alpha_01(n_inputs);
    l4n::MultiIndex alpha_20(n_inputs);

    alpha_00.set_partial_derivative(0,
                                    0);
    alpha_01.set_partial_derivative(1,
                                    1);
    alpha_20.set_partial_derivative(0,
                                    2);

    /* the governing differential equation */
    solver_01.add_to_differential_equation(0,
                                           alpha_20,
                                           "1.0");
    solver_01.add_to_differential_equation(0,
                                           alpha_01,
                                           "-1.0");

    /* dirichlet boundary condition */
    solver_01.add_to_differential_equation(1,
                                           alpha_00,
                                           "1.0");
    solver_01.add_to_differential_equation(2,
                                           alpha_00,
                                           "1.0");

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

    double frac, x, t;

    /* TRAIN DATA FOR THE GOVERNING DE */
    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_zero;

    /* GOVERNING EQUATION RHS */
    frac = (de - ds) / (train_size - 1);
    for (unsigned int i = 0; i < train_size; ++i) {
        for (unsigned int j = 0; j < train_size; ++j) {
            inp = {frac * j, frac * i};
            out = {0.0};
            data_vec_zero.emplace_back(std::make_pair(inp,
                                                      out));
        }
    }
    l4n::DataSet      ds_00(&data_vec_zero);

    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_t;
    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_x;
    /* ISOTROPIC TRAIN SET */
    frac = (de - ds) / (train_size - 1);
    for (unsigned int i = 0; i < train_size; ++i) {
        inp = {0.0, frac * i};
        out = {std::sin(inp[1])};
        data_vec_t.emplace_back(std::make_pair(inp,
                                               out));

        inp = {frac * i, 0.0};
        out = {std::pow(l4n::E,
                        -0.707106781 * inp[0]) * std::sin(-0.707106781 * inp[0])};
        data_vec_x.emplace_back(std::make_pair(inp,
                                               out));

    }
    l4n::DataSet ds_t(&data_vec_t);
    l4n::DataSet ds_x(&data_vec_x);
    std::cout << "Train data setup finished" << std::endl;



    /* 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_t);
    solver_01.set_error_function(2,
                                 l4n::ErrorFunctionType::ErrorFuncMSE,
                                 ds_x);
    std::cout << "Error function defined" << std::endl;

    /* Solving the equation */
    optimize_via_particle_swarm(solver_01,
                                alpha_00,
                                max_iters,
                                n_particles);
    export_solution(n_test_points,
                    te,
                    ts,
                    solver_01,
                    alpha_00,
                    alpha_01,
                    alpha_20,
                    "particle_");

    optimize_via_gradient_descent(solver_01,
                                  accuracy);
    export_solution(n_test_points,
                    te,
                    ts,
                    solver_01,
                    alpha_00,
                    alpha_01,
                    alpha_20,
                    "gradient_");
}

int main() {
    std::cout << "Running lib4neuro Partial Differential Equation example   1" << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;
    std::cout
        << "          Governing equation: y_xx - y_t = 0,                                   for (x, t) in [0, 1] x [0, 1]"
        << std::endl;
    std::cout << "Dirichlet boundary condition:    y(0, t) = sin(t),                              for t in [0, 1]"
              << std::endl;
    std::cout << "Dirichlet boundary condition:    y(x, 0) = exp(-sqrt(0.5)x) * sin(-sqrt(0.5)x), for x in [0, 1]"
              << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;
    std::cout
        << "Expressing solution as y(x, t) = sum over [a_i / (1 + exp(bi - wxi*x - wti*t))], 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      = 5;
    double       accuracy        = 1e-1;
    double       ds              = 0.0;
    double       de              = 1.0;

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

    size_t particle_swarm_max_iters = 10;
    size_t n_particles              = 5;
    test_pde(accuracy,
             n_inner_neurons,
             train_size,
             ds,
             de,
             test_size,
             ts,
             te,
             particle_swarm_max_iters,
             n_particles);

    return 0;
}