From 8de7352d3ff2ce8aadd2b73155e0a6e9613f8f81 Mon Sep 17 00:00:00 2001
From: Michal Kravcenko <michal.kravcenko@vsb.cz>
Date: Wed, 5 Sep 2018 16:59:17 +0200
Subject: [PATCH] ADD: added a new example for harmonic oscilator problem FIX:
fixed an issue with ExprtkWrapper being incorrectly destryed before being
used in the NetworkSum class
---
src/General/ExprtkWrapper.cpp | 35 +++--
src/General/ExprtkWrapper.h | 8 +-
src/Network/NeuralNetworkSum.cpp | 16 +-
src/Network/NeuralNetworkSum.h | 2 +-
src/examples/net_test_harmonic_oscilator.cpp | 156 +++++++++++++++++++
5 files changed, 202 insertions(+), 15 deletions(-)
diff --git a/src/General/ExprtkWrapper.cpp b/src/General/ExprtkWrapper.cpp
index 9ec67262..03738f64 100644
--- a/src/General/ExprtkWrapper.cpp
+++ b/src/General/ExprtkWrapper.cpp
@@ -10,26 +10,41 @@
ExprtkWrapper::ExprtkWrapper( std::string expression_string ) {
+ this->expression_str = expression_string;
+
this->symbol_table = new symbol_table_t( );
- symbol_table->add_variable("x", this->t);
- symbol_table->add_variable("y", this->y);
- symbol_table->add_variable("z", this->z);
- symbol_table->add_variable("t", this->t);
+ this->symbol_table->add_variable("x", this->x);
+ this->symbol_table->add_variable("y", this->y);
+ this->symbol_table->add_variable("z", this->z);
+ this->symbol_table->add_variable("t", this->t);
+ this->symbol_table->add_variable("f", this->z);
this->expression = new expression_t( );
this->expression->register_symbol_table( *this->symbol_table );
this->parser = new parser_t( );
- parser->compile(expression_string, *this->expression );
+ parser->compile(this->expression_str, *this->expression );
}
-ExprtkWrapper::ExprtkWrapper( ) {
+ExprtkWrapper::~ExprtkWrapper() {
-}
+ if( this->expression ){
+ delete this->expression;
+ this->expression = nullptr;
+ }
+
+ if( this->symbol_table ){
+ delete this->symbol_table;
+ this->symbol_table = nullptr;
+ }
+
+ if( this->parser ){
+ delete this->parser;
+ this->parser = nullptr;
+ }
-ExprtkWrapper::~ExprtkWrapper() {
}
double ExprtkWrapper::eval(double x1, double x2, double x3, double x4) {
@@ -58,5 +73,7 @@ double ExprtkWrapper::eval(std::vector<double> &p) {
this->t = p[3];
}
- return this->expression->value();
+ double result = this->expression->value();
+
+ return result;
}
\ No newline at end of file
diff --git a/src/General/ExprtkWrapper.h b/src/General/ExprtkWrapper.h
index 2d38fa03..df58142e 100644
--- a/src/General/ExprtkWrapper.h
+++ b/src/General/ExprtkWrapper.h
@@ -74,7 +74,13 @@ private:
/*
* variables
*/
- double x, y, z, t;
+ double x, y, z, t, f;
+
+ /**
+ * referential expression string
+ */
+
+ std::string expression_str;
};
#endif //LIB4NEURO_EXPRTKWRAPPER_H
diff --git a/src/Network/NeuralNetworkSum.cpp b/src/Network/NeuralNetworkSum.cpp
index 8aee1806..2b665fe7 100644
--- a/src/Network/NeuralNetworkSum.cpp
+++ b/src/Network/NeuralNetworkSum.cpp
@@ -15,9 +15,16 @@ NeuralNetworkSum::NeuralNetworkSum(){
NeuralNetworkSum::~NeuralNetworkSum() {
if( this->summand ){
delete this->summand;
+ this->summand = nullptr;
}
if( this->summand_coefficient ){
+
+ for(auto el: *this->summand_coefficient){
+ delete el;
+ }
+
delete this->summand_coefficient;
+ this->summand_coefficient = nullptr;
}
}
@@ -28,9 +35,9 @@ void NeuralNetworkSum::add_network( NeuralNetwork *net, std::string expression_s
this->summand->push_back( net );
if(!this->summand_coefficient){
- this->summand_coefficient = new std::vector<ExprtkWrapper>(0);
+ this->summand_coefficient = new std::vector<ExprtkWrapper*>(0);
}
- this->summand_coefficient->push_back( ExprtkWrapper( expression_string ) );
+ this->summand_coefficient->push_back( new ExprtkWrapper( expression_string ) );
}
void NeuralNetworkSum::eval_single(std::vector<double> &input, std::vector<double> &output, std::vector<double> *custom_weights_and_biases) {
@@ -45,7 +52,7 @@ void NeuralNetworkSum::eval_single(std::vector<double> &input, std::vector<doubl
if( SUM ){
this->summand->at(ni)->eval_single(input, mem_output, custom_weights_and_biases);
- double alpha = this->summand_coefficient->at(ni).eval(input);
+ double alpha = this->summand_coefficient->at(ni)->eval(input);
for(size_t j = 0; j < output.size(); ++j){
output[j] += mem_output[j] * alpha;
@@ -53,7 +60,8 @@ void NeuralNetworkSum::eval_single(std::vector<double> &input, std::vector<doubl
}
else{
//TODO assume the result can be a vector of doubles
- double alpha = this->summand_coefficient->at(ni).eval(input);
+ double alpha = this->summand_coefficient->at(ni)->eval(input);
+
for(size_t j = 0; j < output.size(); ++j){
output[j] += alpha;
}
diff --git a/src/Network/NeuralNetworkSum.h b/src/Network/NeuralNetworkSum.h
index 3c366430..30f272b4 100644
--- a/src/Network/NeuralNetworkSum.h
+++ b/src/Network/NeuralNetworkSum.h
@@ -19,7 +19,7 @@ private:
friend class boost::serialization::access;
std::vector<NeuralNetwork*> * summand;
- std::vector<ExprtkWrapper> * summand_coefficient;
+ std::vector<ExprtkWrapper*> * summand_coefficient;
template <class Archive>
void serialize(Archive & ar, const unsigned int version) {
diff --git a/src/examples/net_test_harmonic_oscilator.cpp b/src/examples/net_test_harmonic_oscilator.cpp
index 69290d67..17fddd01 100644
--- a/src/examples/net_test_harmonic_oscilator.cpp
+++ b/src/examples/net_test_harmonic_oscilator.cpp
@@ -14,7 +14,163 @@
#include "../../include/4neuro.h"
#include "../Solvers/DESolver.h"
+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;
+ DESolver solver( n_equations, n_inputs, n_inner_neurons );
+
+ /* SETUP OF THE EQUATIONS */
+ MultiIndex alpha_0( n_inputs );
+ 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));
+
+ DataSet ds_00(&data_vec_g);
+
+ /* Placing the conditions into the solver */
+ solver.set_error_function( 0, ErrorFunctionType::ErrorFuncMSE, &ds_00 );
+
+ /* PARTICLE SWARM TRAINING METHOD SETUP */
+ size_t total_dim = (2 + n_inputs) * n_inner_neurons;
+
+ std::vector<double> params(total_dim), params_analytical(total_dim);
+ std::random_device seeder;
+ std::mt19937 gen(seeder());
+ std::uniform_real_distribution<double> dist(-10.0, 10.0);
+
+ std::vector<double> input(1);
+ //must encapsulate each of the partial error functions
+ std::vector<double> domain_bounds(2 * total_dim);
+ for(unsigned int i = 0; i < total_dim; ++i){
+ domain_bounds[2 * i] = -10.0;
+ domain_bounds[2 * i + 1] = 10.0;
+ }
+
+
+ double c1 = 1.7, c2 = 1.7, w = 0.7;
+ /* 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.9;
+ solver.solve_via_particle_swarm( &domain_bounds, c1, c2, w, n_particles, max_iters, gamma, epsilon, delta );
+
+ NeuralNetwork *solution = solver.get_solution( alpha_0 );
+ std::vector<double> parameters(total_dim);//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);
+
+ printf("Path %3d. w%d = %15.8f, b%d = %15.8f, a%d = %15.8f\n", (int)(i + 1), (int)(i + 1), parameters[3 * i], (int)(i + 1), parameters[3 * i + 2], (int)(i + 1), parameters[3 * i + 1]);
+ }
+
+ /* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
+ /* first boundary condition & its error */
+ std::ofstream ofs("data_1d_osc.txt", std::ofstream::out);
+ printf("Exporting files 'data_1d_osc.txt': %7.3f%%\r", 0.0);
+ frac = (te - ts) / (n_test_points - 1);
+
+ for(size_t i = 0; i < n_test_points; ++i){
+ double x = frac * i + ts;
+
+ inp[0] = x;
+ solution->eval_single(inp, out);
+ ofs << i + 1 << " " << x << " " << out[0] << " " << std::endl;
+
+ printf("Exporting files 'data_1d_osc.txt': %7.3f%%\r", (100.0 * i) / (n_test_points - 1));
+ std::cout.flush();
+ }
+ printf("Exporting files 'data_1d_osc.txt': %7.3f%%\n", 100.0);
+ std::cout.flush();
+ ofs.close();
+
+ inp[0] = -1.0;
+ solution->eval_single(inp, out);
+ printf("y(-1) = %f\n", out[0]);
+ inp[0] = 0.0;
+ solution->eval_single(inp, out);
+ printf("y( 0) = %f\n", out[0]);
+ inp[0] = 1.0;
+ solution->eval_single(inp, out);
+ printf("y( 1) = %f\n", out[0]);
+ std::cout << "********************************************************************************************************************************************" <<std::endl;
+}
+
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);
+// std::string expression_string = "-x";
+// std::string expression_string_1 = "1.0";
+// ExprtkWrapper f(expression_string);
+// ExprtkWrapper f1(expression_string_1);
+//
+//
+// f1.eval();
+//
+// std::vector<double> inp(1);
+//
+// inp = {150};
+// double result = f.eval(inp);
+//
+// f1.eval();
+// inp = {15};
+// result = f.eval(inp);
return 0;
}
\ No newline at end of file
--
GitLab