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 "lib4neuro.h"

void optimize_via_particle_swarm( DESolver &solver, 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;

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

    solver.solve( swarm );

}

void optimize_via_gradient_descent( DESolver &solver, double accuracy ){
    printf("Solution via a gradient descent method!\n");
    GradientDescent gd( accuracy, 1000 );

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

void export_solution( size_t n_test_points, double te, double ts, DESolver &solver, MultiIndex &alpha_00, MultiIndex &alpha_01, MultiIndex &alpha_20, const std::string prefix ){
    NeuralNetwork *solution = solver.get_solution( alpha_00 );
    NeuralNetwork *solution_t = solver.get_solution( alpha_01 );
    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(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;
    DESolver solver_01( n_equations, n_inputs, n_inner_neurons );

    /* SETUP OF THE EQUATIONS */
    MultiIndex alpha_00( n_inputs );
    MultiIndex alpha_01( n_inputs );
    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<double> test_bounds_2d = {ds, de, ds, de};

    /* GOVERNING EQUATION RHS */
    auto f1 = [](std::vector<double>&input) -> std::vector<double> {
        std::vector<double> output(1);
        output[0] = 0.0;
        return output;
    };
    DataSet ds_00(test_bounds_2d, train_size, f1, 1);

    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(E, -0.707106781 * inp[0]) * std::sin( -0.707106781 * inp[0] )};
        data_vec_x.emplace_back(std::make_pair(inp, out));

    }
    DataSet ds_t(&data_vec_t);
    DataSet ds_x(&data_vec_x);




    /* 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_t );
    solver_01.set_error_function( 2, ErrorFunctionType::ErrorFuncMSE, &ds_x );

    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 = 4;
    unsigned int train_size = 50;
    double accuracy = 1e-3;
    double ds = 0.0;
    double de = 1.0;

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

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

    return 0;
}