/**
* Example solving the following ODE:
*
* g(t) = (d^2/d^t)y(t) + 4 (d/dt)y(t) + 4y(t) = 0, for t in [0, 4]
* y(0) = 1
* (d/dt)y(0) = 1
*
* -------------------------------------------
* Analytical solution: e^(-2x) * (3x + 1)
* NN representation: sum over [a_i * (1 + e^(-x * w_i + b_i))^(-1)]
* -------------------------------------------
* Optimal NN setting with biases (2 inner neurons)
* Path 1. w = -1.66009975, b = -0.40767447, a = 2.46457042
* Path 2. w = -4.38622765, b = 2.75707816, a = -8.04752347
* @author Michal Kravčenko
* @date 17.7.18 -
*/
#include <random>
#include <iostream>
#include <chrono>
#include <4neuro.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");
solver.randomize_parameters();
// TODO does not work (poor design of netsum)
// l4n::LevenbergMarquardt leven(10000, 0, 1e-6, 1e-6, 1e-6);
// solver.solve(leven);
l4n::GradientDescent gd(accuracy,
1000,
500000);
solver.solve(gd);
}
void export_solution(size_t n_test_points,
double te,
double ts,
l4n::DESolver& solver,
l4n::MultiIndex& alpha_0,
l4n::MultiIndex& alpha_1,
l4n::MultiIndex& alpha_2,
const std::string prefix) {
l4n::NeuralNetwork* solution = solver.get_solution(alpha_0);
l4n::NeuralNetwork* solution_d = solver.get_solution(alpha_1);
l4n::NeuralNetwork* solution_dd = solver.get_solution(alpha_2);
/* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
/* first boundary condition & its error */
char buff[256];
sprintf(buff,
"%sdata_1d_ode1.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);
std::vector<double> inp(1), out(1);
for (size_t i = 0; i < n_test_points; ++i) {
double x = frac * i + ts;
inp[0] = x;
solution->eval_single(inp,
out);
double F = out[0];
solution_d->eval_single(inp,
out);
double DF = out[0];
solution_dd->eval_single(inp,
out);
double DDF = out[0];
ofs << i + 1 << " " << x << " " << std::pow(l4n::E,
-2 * x) * (3 * x + 1) << " " << F << " "
<< std::pow(l4n::E,
-2 * x) * (1 - 6 * x) << " " << DF << " " << 4 * std::pow(l4n::E,
-2 * x) * (3 * x - 2)
<< " " << DDF << 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%%\r",
final_fn.c_str(),
100.0);
std::cout.flush();
ofs.close();
std::cout
<< "********************************************************************************************************************************************"
<< std::endl;
}
void test_ode(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 and Gradient descent method!" << std::endl;
std::cout
<< "********************************************************************************************************************************************"
<< std::endl;
/* SOLVER SETUP */
size_t n_inputs = 1;
size_t n_equations = 3;
l4n::DESolver solver_01(n_equations,
n_inputs,
n_inner_neurons);
/* SETUP OF THE EQUATIONS */
l4n::MultiIndex alpha_0(n_inputs);
l4n::MultiIndex alpha_1(n_inputs);
l4n::MultiIndex alpha_2(n_inputs);
alpha_2.set_partial_derivative(0,
2);
alpha_1.set_partial_derivative(0,
1);
/* the governing differential equation */
solver_01.add_to_differential_equation(0,
alpha_0,
"4.0");
solver_01.add_to_differential_equation(0,
alpha_1,
"4.0");
solver_01.add_to_differential_equation(0,
alpha_2,
"1.0");
/* dirichlet boundary condition */
solver_01.add_to_differential_equation(1,
alpha_0,
"1.0");
/* neumann boundary condition */
solver_01.add_to_differential_equation(2,
alpha_1,
"1.0");
/* SETUP OF THE TRAINING DATA */
std::vector<double> inp(1), out(1);
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;
std::vector<double> test_points(train_size);
/* ISOTROPIC TRAIN SET */
frac = (d1_e - d1_s) / (train_size - 1);
for (unsigned int i = 0; i < train_size; ++i) {
inp[0] = frac * i;
out[0] = 0.0;
data_vec_g.push_back(std::make_pair(inp,
out));
test_points[i] = inp[0];
}
/* CHEBYSCHEV TRAIN SET */
l4n::DataSet ds_00(&data_vec_g);
/* TRAIN DATA FOR DIRICHLET BC */
std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_y;
inp = {0.0};
out = {1.0};
data_vec_y.emplace_back(std::make_pair(inp,
out));
l4n::DataSet ds_01(&data_vec_y);
/* TRAIN DATA FOR NEUMANN BC */
std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_dy;
inp = {0.0};
out = {1.0};
data_vec_dy.emplace_back(std::make_pair(inp,
out));
l4n::DataSet ds_02(&data_vec_dy);
/* 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_01);
solver_01.set_error_function(2,
l4n::ErrorFunctionType::ErrorFuncMSE,
ds_02);
/* TRAINING METHOD SETUP */
/* optimize_via_particle_swarm( solver_01, alpha_0, max_iters, n_particles );
export_solution( n_test_points, te, ts, solver_01 , alpha_0, alpha_1, alpha_2, "particle_" );*/
auto start = std::chrono::system_clock::now();
optimize_via_gradient_descent(solver_01,
accuracy);
export_solution(n_test_points,
te,
ts,
solver_01,
alpha_0,
alpha_1,
alpha_2,
"gradient_");
auto end = std::chrono::system_clock::now();
std::chrono::duration<double> elapsed_seconds = end - start;
std::cout << "elapsed time: " << elapsed_seconds.count() << std::endl;
}
int main() {
std::cout << "Running lib4neuro Ordinary Differential Equation example 1" << std::endl;
std::cout
<< "********************************************************************************************************************************************"
<< std::endl;
std::cout << " Governing equation: y''(x) + 4y'(x) + 4y(x) = 0.0, for x in [0, 4]" << std::endl;
std::cout << "Dirichlet boundary condition: y(0.0) = 1.0" << std::endl;
std::cout << " Neumann boundary condition: y'(0.0) = 1.0" << 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;
unsigned int n_inner_neurons = 2;
unsigned int train_size = 10;
double accuracy = 1e-1;
double ds = 0.0;
double de = 4.0;
unsigned int test_size = 10;
double ts = ds;
double te = de + 2;
size_t particle_swarm_max_iters = 10;
size_t n_particles = 2;
test_ode(accuracy,
n_inner_neurons,
train_size,
ds,
de,
test_size,
ts,
te,
particle_swarm_max_iters,
n_particles);
return 0;
}