Skip to content
Snippets Groups Projects
net_test_pde_1.cpp 5 KiB
Newer Older
  • Learn to ignore specific revisions
  • Michal Kravcenko's avatar
    Michal Kravcenko committed
    /**
     * Example solving the water flow 1D diffusion PDE:
     *
     * (d^2/d^xx^x)y(x, t) - (d/dt)y(x, t) = 0, for t in [0, 1] x [0, 1]
     * y(0, t) = sin(t)
     * y(x, 0) = e^(-(0.5)^(0.5)x) * sin(-(0.5)^(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 (2 inner neurons)
     *
     * @author Michal Kravčenko
     * @date 9.8.18
     */
    
    #include <random>
    #include <iostream>
    
    
    Michal Kravcenko's avatar
    Michal Kravcenko committed
    #include "../../include/4neuro.h"
    #include "../Solvers/DESolver.h"
    
    Michal Kravcenko's avatar
    Michal Kravcenko committed
    
    
    void test_odr(size_t n_inner_neurons){
    
        /* solution properties */
        size_t train_size = 10;
        double d1_s = 0.0, d1_e = 1.0;
    
    
    Michal Kravcenko's avatar
    Michal Kravcenko committed
        /* swarm optimizer properties */
    
    Michal Kravcenko's avatar
    Michal Kravcenko committed
        unsigned int max_iters = 100;
        unsigned int n_particles = 10;
    
    
        /* 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(0, 1);
        alpha_20.set_partial_derivative(2, 0);
    
        /* 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 */
        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_g;
    
    
        /* ISOTROPIC TRAIN SET */
        frac = (d1_e - d1_s) / (train_size - 1);
        for(size_t i = 0; i < train_size; ++i){
            inp = {0.0, 0.0};
            out = {0.0};
            data_vec_g.emplace_back(std::make_pair(inp, out));
        }
        DataSet ds_00(&data_vec_g);
    
        /* ISOTROPIC TRAIN DATA FOR THE FIRST DIRICHLET BC */
        std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_t;
        frac = (d1_e - d1_s) / (train_size - 1);
        for(size_t i = 0; i < train_size; ++i){
            t = i * frac;
            inp = { 0.0, t };
            out = { std::sin( t ) };
    
            data_vec_t.emplace_back(std::make_pair(inp, out));
        }
        DataSet ds_t(&data_vec_t);
    
        /* ISOTROPIC TRAIN DATA FOR THE SECOND DIRICHLET BC */
        std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_x;
        frac = (d1_e - d1_s) / (train_size - 1);
        for(size_t i = 0; i < train_size; ++i){
            x = i * frac;
            inp = { x, 0.0 };
            out = { std::pow(E, -0.707106781 * x) * std::sin( -0.707106781 * x )};
    
            data_vec_x.emplace_back(std::make_pair(inp, out));
        }
        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 );
    
    
    
        /* PARTICLE SWARM TRAINING METHOD SETUP */
    
    
        //must encapsulate each of the partial error functions
        double *domain_bounds = new double[ 6 * n_inner_neurons ];
        for(unsigned int i = 0; i < 3 * n_inner_neurons; ++i){
            domain_bounds[2 * i] = -800.0;
            domain_bounds[2 * i + 1] = 800.0;
        }
    
        double c1 = 0.5, c2 = 0.75, w = 0.8;
    
    
    
        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();
    //    std::vector<double> parameters(3 * n_inner_neurons);//w1, a1, b1, w2, a2, b2, ... , wm, am, bm
    //    std::vector<double> *weight_params = solution->get_parameter_ptr_weights();
    //    std::vector<double> *biases_params = solution->get_parameter_ptr_biases();
    //    for(size_t i = 0; i < n_inner_neurons; ++i){
    //        parameters[3 * i] = weight_params->at(i);
    //        parameters[3 * i + 1] = weight_params->at(i + n_inner_neurons);
    //        parameters[3 * i + 2] = biases_params->at(i);
    //    }
    //
    //    unsigned int n_test_points = 150;
    //    std::vector<double> input(1), output(1);
    //    double x;
    //    for(unsigned int i = 0; i < n_test_points; ++i){
    //
    //        x = i * ((d1_e - d1_s) / (n_test_points - 1)) + d1_s;
    //        input[0] = x;
    //
    //        solution->eval_single(input, output);
    //
    //        std::cout << i + 1 << " " << x << " " << std::pow(E, -2*x) * (3*x + 1)<< " " << output[0] << " " << std::pow(E, -2*x) * (1 - 6*x)<< " " << eval_approx_df(x, n_inner_neurons, parameters) << " " << 4 * std::pow(E, -2*x) * (3*x - 2)<< " " << eval_approx_ddf(x, n_inner_neurons, parameters) << std::endl;
    //    }
    
    
        delete [] domain_bounds;
    }
    
    int main() {
    
        unsigned int n_inner_neurons = 2;
    
        test_odr(n_inner_neurons);
    
    
        return 0;
    }