net_test_harmonic_oscilator.cpp

/**
 * Example solving the eigenvalue problem:
 *
 *
 *
 * @author Michal Kravčenko
 * @date 3.9.18 -
 */

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

#include <4neuro.h>

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

    char buff[256];
    sprintf(buff,
            "%sdata_1d_osc.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), x;

    std::vector<double> inp(1), out(1);

    for (size_t i = 0; i < n_test_points; ++i) {
        x = frac * i + ts;

        inp[0] = x;
        solution->eval_single(inp,
                              out);
        ofs << i + 1 << " " << x << " " << out[0] << " " << 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%%\n",
           final_fn.c_str(),
           100.0);
    std::cout.flush();
    ofs.close();
}

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 test_harmonic_oscilator_fixed_E(double EE,
                                     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" << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;

    /* SOLVER SETUP */
    size_t        n_inputs    = 1;
    size_t        n_equations = 1;
    l4n::DESolver solver(n_equations,
                         n_inputs,
                         n_inner_neurons);

    /* SETUP OF THE EQUATIONS */
    l4n::MultiIndex alpha_0(n_inputs);
    l4n::MultiIndex alpha_2(n_inputs);
    alpha_2.set_partial_derivative(0,
                                   2);

    /* the governing differential equation */
    char buff[255];
    std::sprintf(buff,
                 "%f",
                 -EE);
    std::string eigenvalue(buff);
    solver.add_to_differential_equation(0,
                                        alpha_2,
                                        "-1.0");
    solver.add_to_differential_equation(0,
                                        alpha_0,
                                        "x^2");
    solver.add_to_differential_equation(0,
                                        alpha_0,
                                        eigenvalue);

    /* 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;


    /* ISOTROPIC TRAIN SET */
    frac = (d1_e - d1_s) / (train_size - 1);
    for (unsigned int i = 0; i < train_size; ++i) {
        inp = {frac * i + d1_s};
        out = {0.0};
        data_vec_g.emplace_back(std::make_pair(inp,
                                               out));
    }
    inp                 = {0.0};
    out                 = {1.0};
    data_vec_g.emplace_back(std::make_pair(inp,
                                           out));

    l4n::DataSet ds_00(&data_vec_g);

    /* Placing the conditions into the solver */
    solver.set_error_function(0,
                              l4n::ErrorFunctionType::ErrorFuncMSE,
                              ds_00);

    /* PARTICLE SWARM TRAINING METHOD SETUP */
    size_t total_dim = (2 + n_inputs) * n_inner_neurons;

    optimize_via_gradient_descent(solver,
                                  accuracy);
    export_solution(n_test_points,
                    te,
                    ts,
                    solver,
                    alpha_0,
                    "gradient_");
}

int main() {
    std::cout << "Running lib4neuro harmonic Oscilator example   1" << std::endl;
    std::cout
        << "********************************************************************************************************************************************"
        << std::endl;
    std::cout << "          Governing equation: -y''(x) + x^2 * y(x) = E * y(x)" << 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;

    double       EE              = -1.0;
    unsigned int n_inner_neurons = 2;
    unsigned int train_size      = 10;
    double       accuracy        = 1e-3;
    double       ds              = -5.0;
    double       de              = 5.0;

    unsigned int test_size = 300;
    double       ts        = -6.0;
    double       te        = 6.0;

    size_t particle_swarm_max_iters = 1000;
    size_t n_particles              = 100;
    test_harmonic_oscilator_fixed_E(EE,
                                    accuracy,
                                    n_inner_neurons,
                                    train_size,
                                    ds,
                                    de,
                                    test_size,
                                    ts,
                                    te,
                                    particle_swarm_max_iters,
                                    n_particles);

    return 0;
}