net_test_pde_1.cpp 4.96 KB
Newer Older
Michal Kravcenko's avatar
Michal Kravcenko committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
/**
 * 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>

#include "../include/4neuro.h"
#include "Solvers/DESolver.h"


void test_odr(size_t n_inner_neurons){

    /* solution properties */
    size_t train_size = 10;
    double d1_s = 0.0, d1_e = 1.0;

    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;
}