diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e216e3a62e15b23f69451e52694f6be30f760292..10f18a2f0a957d755f80a531fa4bf62c2cc0be2e 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -19,11 +19,11 @@ if ("${BUILD_LIB}" STREQUAL "yes")
 		NetConnection/ConnectionFunctionGeneral.cpp        
 		NetConnection/ConnectionFunctionIdentity.cpp        
 		LearningMethods/ParticleSwarm.cpp        
-		DataSet/DataSet.cpp        
+		DataSet/DataSet.cpp
 		ErrorFunction/ErrorFunctions.cpp        
 		Solvers/DESolver.cpp
 		General/ExprtkWrapper.cpp
-		LearningMethods/ILearningMethods.h
+		LearningMethods/ILearningMethods.cpp
 	)
 
 	if(WIN32)
diff --git a/src/LearningMethods/ILearningMethods.cpp b/src/LearningMethods/ILearningMethods.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6aa47daf102bea9655ef0dc8f33c9dce54073092
--- /dev/null
+++ b/src/LearningMethods/ILearningMethods.cpp
@@ -0,0 +1,9 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 10.9.18 -
+ */
+
+#include "ILearningMethods.h"
+
diff --git a/src/LearningMethods/ILearningMethods.h b/src/LearningMethods/ILearningMethods.h
index 31578344de9946b76df7aef97b8640add6d8e872..c49891fd0cc07ccefcf481d24e1058f42d17185a 100644
--- a/src/LearningMethods/ILearningMethods.h
+++ b/src/LearningMethods/ILearningMethods.h
@@ -12,16 +12,23 @@
 #include "../ErrorFunction/ErrorFunctions.h"
 
 class ILearningMethods {
+private:
+
+    /**
+     *
+     */
+    ErrorFunction *ef = nullptr;
+
 public:
     /*
      * Runs the method specific learning algorithm minimizing the given error function
      */
-    virtual void optimize( ErrorFunction &error ) = 0;
+    virtual void optimize( ErrorFunction &ef ) = 0;
 
     /*
      * Updates the optimal weight&bias settings in the passed vector
      */
-    virtual void get_optimal_configuration( std::vector<double> &config ) = 0;
+    virtual std::vector<double>* get_parameters( ) = 0;
 };
 
 
diff --git a/src/LearningMethods/ParticleSwarm.cpp b/src/LearningMethods/ParticleSwarm.cpp
index 85a8bd280ca71789cfacd27ef4e270f3d8ca3507..763932404ea155386244b962e30db4173006ad32 100644
--- a/src/LearningMethods/ParticleSwarm.cpp
+++ b/src/LearningMethods/ParticleSwarm.cpp
@@ -47,8 +47,9 @@ void Particle::randomize_velocity() {
     }
 }
 
-Particle::Particle(ErrorFunction* ef, double *domain_bounds) {
-    //TODO better generating of random numbers
+Particle::Particle(ErrorFunction *ef, std::vector<double> *domain_bounds) {
+
+    this->ef = ef;
     this->domain_bounds = domain_bounds;
     this->coordinate_dim = ef->get_dimension();
     this->ef = ef;
@@ -62,8 +63,6 @@ Particle::Particle(ErrorFunction* ef, double *domain_bounds) {
     this->randomize_parameters();
     this->randomize_coordinates();
 
-
-
     for(unsigned int i = 0; i < this->coordinate_dim; ++i){
         (*this->optimal_coordinate)[i] = (*this->coordinate)[i];
     }
@@ -120,7 +119,7 @@ double Particle::change_coordinate(double w, double c1, double c2, std::vector<d
     std::vector<double> *random_global_best;
     std::random_device rand_dev;
     std::mt19937 engine{rand_dev()};
-    std::uniform_int_distribution<int> dist(0, global_min_vec.size() - 1);
+    std::uniform_int_distribution<size_t> dist(0, global_min_vec.size() - 1);
     random_global_best = &global_min_vec[dist(engine)];
 
     // TODO use std::sample to choose random vector
@@ -132,12 +131,12 @@ double Particle::change_coordinate(double w, double c1, double c2, std::vector<d
                   + c2 * this->r2 * (glob_min_coord[i] - (*this->coordinate)[i])
                   + (c1+c2)/2 * this->r3 * ((*random_global_best)[i] - (*this->coordinate)[i]);
 
-        if ((*this->coordinate)[i] + vel_mem > this->domain_bounds[2 * i + 1]) {
+        if ((*this->coordinate)[i] + vel_mem > this->domain_bounds->at(2 * i + 1)) {
             this->randomize_velocity();
             this->randomize_parameters();
             this->randomize_coordinates();
             break;
-        } else if ((*this->coordinate)[i] + vel_mem < this->domain_bounds[2 * i]) {
+        } else if ((*this->coordinate)[i] + vel_mem < this->domain_bounds->at(2 * i)) {
             this->randomize_velocity();
             this->randomize_parameters();
             this->randomize_coordinates();
@@ -178,17 +177,25 @@ void Particle::print_coordinate() {
     printf("%10.8f\n", (*this->coordinate)[this->coordinate_dim - 1]);
 }
 
-ParticleSwarm::ParticleSwarm(ErrorFunction* ef, std::vector<double> *domain_bounds,
-                             double c1, double c2, double w, size_t n_particles, size_t iter_max) {
+ParticleSwarm::ParticleSwarm(
+        std::vector<double> *domain_bounds,
+        double c1,
+        double c2,
+        double w,
+        double gamma,
+        double epsilon,
+        double delta,
+        size_t n_particles,
+        size_t iter_max
+    ) {
+
     srand(time(NULL));
 
-    if(domain_bounds->size() < 2 * ef->get_dimension()){
-        std::cerr << "The supplied domain bounds dimension is too low! It should be at least " << 2 * ef->get_dimension() << "\n" << std::endl;
-    }
 
-    this->f = ef;
 
-    this->func_dim = ef->get_dimension();
+    if(epsilon < 0 || gamma < 0 || delta < 0) {
+        throw std::invalid_argument("Parameters 'gamma', 'epsilon' and 'delta' must be greater than or equal to zero!");
+    }
 
     this->c1 = c1;
 
@@ -198,16 +205,21 @@ ParticleSwarm::ParticleSwarm(ErrorFunction* ef, std::vector<double> *domain_boun
 
     this->w = w;
 
+    this->gamma = gamma;
+
+    this->epsilon = epsilon;
+
+    this->delta = delta;
+
     this->n_particles = n_particles;
 
     this->iter_max = iter_max;
 
-    this->particle_swarm = new Particle*[this->n_particles];
-    this->domain_bounds = &(domain_bounds->at(0));
+    this->particle_swarm = new std::vector<Particle*>( this->n_particles );
 
-    for( size_t pi = 0; pi < this->n_particles; ++pi ){
-        this->particle_swarm[pi] = new Particle(ef, this->domain_bounds);
-    }
+    this->domain_bounds = domain_bounds;
+
+    std::fill( this->particle_swarm->begin(), this->particle_swarm->end(), nullptr );
 
 }
 
@@ -215,10 +227,11 @@ ParticleSwarm::~ParticleSwarm() {
 
     if( this->particle_swarm ){
         for( size_t i = 0; i < this->n_particles; ++i ){
-            delete this->particle_swarm[i];
+            delete this->particle_swarm->at( i );
         }
 
-        delete [] this->particle_swarm;
+        delete this->particle_swarm;
+        this->particle_swarm = nullptr;
     }
 
     if(this->p_min_glob){
@@ -234,13 +247,28 @@ ParticleSwarm::~ParticleSwarm() {
  * @param epsilon
  * @param delta
  */
-void ParticleSwarm::optimize( double gamma, double epsilon, double delta) {
-    if(epsilon < 0 || gamma < 0 || delta < 0) {
-        throw std::invalid_argument("Parameters 'gamma', 'epsilon' and 'delta' must be greater than or equal to zero!");
+void ParticleSwarm::optimize( ErrorFunction &ef ) {
+
+    if(this->domain_bounds->size() < 2 * ef.get_dimension()){
+        std::cerr << "The supplied domain bounds dimension is too low! It should be at least " << 2 * ef.get_dimension() << "\n" << std::endl;
     }
+
+    this->func_dim = ef.get_dimension();
+
+    /* initialize the particles */
+    for( size_t pi = 0; pi < this->n_particles; ++pi ){
+        if(this->particle_swarm->at( pi )){
+            delete this->particle_swarm->at( pi );
+        }
+        this->particle_swarm->at( pi ) = new Particle( &ef, this->domain_bounds );
+    }
+
     if(!this->p_min_glob){
         this->p_min_glob = new std::vector<double>(this->func_dim);
     }
+    else{
+        this->p_min_glob->resize(this->func_dim);
+    }
 
     size_t outer_it = 0;
     Particle *particle;
@@ -297,7 +325,7 @@ void ParticleSwarm::optimize( double gamma, double epsilon, double delta) {
             }
 
             for(size_t pi=0; pi < this->n_particles; pi++) {
-                particle = this->particle_swarm[pi];
+                particle = this->particle_swarm->at(pi);
                 tmp_velocity = particle->change_coordinate( this->w, this->c1, this->c2, *this->p_min_glob, global_best_vec);
 //                particle->print_coordinate();
 
@@ -313,7 +341,7 @@ void ParticleSwarm::optimize( double gamma, double epsilon, double delta) {
                 // only for info purposes
                 euclidean_dist += this->get_euclidean_distance(particle->get_coordinate(), centroid);
 
-                if(this->get_euclidean_distance(particle->get_coordinate(), centroid) < epsilon) {
+                if(this->get_euclidean_distance(particle->get_coordinate(), centroid) < this->epsilon) {
                     cluster.insert(particle);
                 }
             }
@@ -338,7 +366,7 @@ void ParticleSwarm::optimize( double gamma, double epsilon, double delta) {
 //        }
 
         /* Check if the particles are near to each other AND the maximal velocity is less than 'gamma' */
-        if( cluster.size() > delta * this->n_particles && prev_max_vel_step < gamma * max_vel_step ) {
+        if( cluster.size() > this->delta * this->n_particles && prev_max_vel_step < this->gamma * max_vel_step ) {
             break;
         }
 
@@ -368,18 +396,18 @@ Particle* ParticleSwarm::determine_optimal_coordinate_and_value(std::vector<doub
 
     Particle* p;
 
-    val = this->particle_swarm[0]->get_optimal_value( );
-    this->particle_swarm[0]->get_optimal_coordinate(coord);
-    p = this->particle_swarm[0];
+    val = this->particle_swarm->at(0)->get_optimal_value( );
+    this->particle_swarm->at(0)->get_optimal_coordinate(coord);
+    p = this->particle_swarm->at(0);
 
     for(size_t i = 1; i < this->n_particles; ++i){
 
-        double val_m = this->particle_swarm[i]->get_optimal_value( );
+        double val_m = this->particle_swarm->at(i)->get_optimal_value( );
 
         if(val_m < val){
             val = val_m;
-            this->particle_swarm[i]->get_optimal_coordinate(coord);
-            p = this->particle_swarm[i];
+            this->particle_swarm->at(i)->get_optimal_coordinate(coord);
+            p = this->particle_swarm->at(i);
         }
     }
 
@@ -391,7 +419,7 @@ std::vector<double>* ParticleSwarm::get_centroid_coordinates() {
     std::vector<double>* tmp;
 
     for (size_t pi = 0; pi < this->n_particles; pi++) {
-        tmp = this->particle_swarm[pi]->get_coordinate();
+        tmp = this->particle_swarm->at(pi)->get_coordinate();
 
         for (size_t di = 0; di < this->func_dim; di++) {
             (*coords)[di] += (*tmp)[di];
@@ -415,6 +443,6 @@ double ParticleSwarm::get_euclidean_distance(std::vector<double>* a, std::vector
     return std::sqrt(dist);
 }
 
-std::vector<double>* ParticleSwarm::get_solution() {
+std::vector<double>* ParticleSwarm::get_parameters( ) {
     return this->p_min_glob;
 }
\ No newline at end of file
diff --git a/src/LearningMethods/ParticleSwarm.h b/src/LearningMethods/ParticleSwarm.h
index e0e87f9131fcc1fa751e06dfded6660a9ad96369..aa33816745689b0996838157dd03deb19c882e31 100644
--- a/src/LearningMethods/ParticleSwarm.h
+++ b/src/LearningMethods/ParticleSwarm.h
@@ -22,6 +22,7 @@
 #include "Network/NeuralNetwork.h"
 #include "DataSet/DataSet.h"
 #include "ErrorFunction/ErrorFunctions.h"
+#include "ILearningMethods.h"
 
 
 class Particle{
@@ -42,7 +43,7 @@ private:
 
     ErrorFunction* ef;
 
-    double *domain_bounds;
+    std::vector<double> *domain_bounds;
 
 
     void randomize_coordinates();
@@ -62,7 +63,7 @@ public:
      *
      * @param f_dim
      */
-    LIB4NEURO_API Particle(ErrorFunction* ef, double *domain_bounds);
+    LIB4NEURO_API Particle( ErrorFunction *ef, std::vector<double> *domain_bounds );
     LIB4NEURO_API ~Particle( );
 
     /**
@@ -101,14 +102,14 @@ public:
 };
 
 
-class ParticleSwarm {
+class ParticleSwarm: public ILearningMethods  {
 
 private:
 
     /**
      *
      */
-    Particle** particle_swarm = nullptr;
+    std::vector<Particle*> *particle_swarm = nullptr;
 
     /**
      *
@@ -129,9 +130,15 @@ private:
 
     double w;
 
+    double gamma;
+
+    double epsilon;
+
+    double delta;
+
     double global_optimal_value;
 
-    double *domain_bounds;
+    std::vector<double> *domain_bounds = nullptr;
 
     std::vector<double> *p_min_glob = nullptr;
 
@@ -173,7 +180,17 @@ public:
      * @param iter_max
      */
      //TODO make domain_bounds constant
-    LIB4NEURO_API ParticleSwarm( ErrorFunction* ef, std::vector<double> *domain_bounds, double c1 = 1.711897, double c2 = 1.711897, double w = 0.711897, size_t n_particles = 50, size_t iter_max = 1000 );
+    LIB4NEURO_API ParticleSwarm(
+            std::vector<double> *domain_bounds,
+            double c1 = 1.711897,
+            double c2 = 1.711897,
+            double w = 0.711897,
+            double gamma = 0.5,
+            double epsilon = 0.02,
+            double delta = 0.7,
+            size_t n_particles = 50,
+            size_t iter_max = 1000
+     );
 
     /**
      *
@@ -187,13 +204,13 @@ public:
      * @param epsilon
      * @param delta
      */
-    LIB4NEURO_API void optimize( double gamma, double epsilon, double delta=0.7 );
+    LIB4NEURO_API void optimize( ErrorFunction &ef ) override;
 
     /**
      *
      * @return
      */
-    LIB4NEURO_API std::vector<double>* get_solution();
+    LIB4NEURO_API std::vector<double>* get_parameters( ) override;
 
 
 };
diff --git a/src/Solvers/DESolver.cpp b/src/Solvers/DESolver.cpp
index f84fd9c31dc1a03214fa00ffd0e89c66065f5a49..7889202723aa495cd16ad20f08146a6f7c78fa9f 100644
--- a/src/Solvers/DESolver.cpp
+++ b/src/Solvers/DESolver.cpp
@@ -332,9 +332,7 @@ void DESolver::set_error_function(size_t equation_idx, ErrorFunctionType F, Data
 }
 
 //TODO instead use general method with Optimizer as its argument (create hierarchy of optimizers)
-void DESolver::solve_via_particle_swarm(std::vector<double> *domain_bounds, double c1, double c2, double w,
-                                          size_t n_particles, size_t max_iters, double gamma,
-                                          double epsilon, double delta) {
+void DESolver::solve( ILearningMethods &learning_method ) {
 
     NeuralNetwork *nn;
     DataSet *ds;
@@ -358,14 +356,12 @@ void DESolver::solve_via_particle_swarm(std::vector<double> *domain_bounds, doub
         }
     }
 
-    ParticleSwarm swarm_01(&total_error, domain_bounds, c1, c2, w, n_particles, max_iters);
-
     this->solution->randomize_weights();
     this->solution->randomize_biases();
 
-    swarm_01.optimize(gamma, epsilon, delta);
+    learning_method.optimize( total_error );
 
-    this->solution->copy_parameter_space(swarm_01.get_solution());
+    this->solution->copy_parameter_space( learning_method.get_parameters( ) );
 }
 
 NeuralNetwork* DESolver::get_solution( MultiIndex &alpha ) {
diff --git a/src/Solvers/DESolver.h b/src/Solvers/DESolver.h
index ecc0ad194e01b95987f43ab6f5724b501129ea86..68a3f978c09143b05feaa7c8287dbb5ff4173961 100644
--- a/src/Solvers/DESolver.h
+++ b/src/Solvers/DESolver.h
@@ -144,17 +144,7 @@ public:
 
 
 
-    LIB4NEURO_API void solve_via_particle_swarm(
-            std::vector<double> *domain_bounds,
-            double c1,
-            double c2,
-            double w,
-            size_t n_particles,
-            size_t max_iters,
-            double gamma,
-            double epsilon,
-            double delta
-            );
+    LIB4NEURO_API void solve( ILearningMethods &learning_method );
 
     /**
      * returns the pointer to the object representing the given partial derivative of the solution
diff --git a/src/examples/net_test_1.cpp b/src/examples/net_test_1.cpp
index 0c3cff4c07aa9c9c5ecd21a77179207d7d2bde3d..8d032d8b8a8c8f7c2b8c1c42a4fedbf3c8b3f89e 100644
--- a/src/examples/net_test_1.cpp
+++ b/src/examples/net_test_1.cpp
@@ -71,8 +71,18 @@ int main() {
     MSE mse(&net, &ds);
 
     /* TRAINING METHOD SETUP */
-    std::vector<double> domain_bounds = {-10.0, 10.0, -10.0, 10.0, -10.0, 10.0};
-    ParticleSwarm swarm_01(&mse, &domain_bounds);
+    std::vector<double> domain_bounds(2 * (net.get_n_weights() + net.get_n_biases()));
+
+    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.7;
+    size_t n_particles = 50;
+    size_t iter_max = 1000;
 
     /* if the maximal velocity from the previous step is less than 'gamma' times the current maximal velocity, then one
      * terminating criterion is met */
@@ -82,9 +92,21 @@ int main() {
      * terminating criterion is met ('n' is the total number of particles) */
     double epsilon = 0.02;
     double delta = 0.7;
-    swarm_01.optimize(gamma, epsilon, delta);
 
-    std::vector<double> *parameters = swarm_01.get_solution();
+    ParticleSwarm swarm_01(
+            &domain_bounds,
+            c1,
+            c2,
+            w,
+            gamma,
+            epsilon,
+            delta,
+            n_particles,
+            iter_max
+    );
+    swarm_01.optimize( mse );
+
+    std::vector<double> *parameters = swarm_01.get_parameters();
     net.copy_parameter_space(parameters);
 
     printf("w1 = %10.7f\n", parameters->at( 0 ));
diff --git a/src/examples/net_test_2.cpp b/src/examples/net_test_2.cpp
index c8b84f1c54bbb960b8c6d39ec12d5d66c9be2343..709b2af4c18afb32278afd51d011f183c36e0d3c 100644
--- a/src/examples/net_test_2.cpp
+++ b/src/examples/net_test_2.cpp
@@ -43,7 +43,6 @@ int main() {
     NeuronLinear *i1 = new NeuronLinear( );  //f(x) = x
     NeuronLinear *i2 = new NeuronLinear( );  //f(x) = x
 
-    double b = 1;//bias
     NeuronLinear *i3 = new NeuronLinear( );  //f(x) = x
 
     /* Output neurons */
@@ -91,9 +90,18 @@ int main() {
 //    printf("evaluation of error at point (%f, %f) => %f\n", weights[0], weights[1], mse.eval(weights));
 
     /* TRAINING METHOD SETUP */
-    //must encapsulate each of the partial error functions
-    std::vector<double> domain_bounds = {-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0};
-    ParticleSwarm swarm_01(&mse, &domain_bounds);
+    std::vector<double> domain_bounds(2 * (net.get_n_weights() + net.get_n_biases()));
+
+    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.7;
+    size_t n_particles = 50;
+    size_t iter_max = 1000;
 
     /* if the maximal velocity from the previous step is less than 'gamma' times the current maximal velocity, then one
      * terminating criterion is met */
@@ -102,10 +110,22 @@ int main() {
     /* 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;
-    swarm_01.optimize(gamma, epsilon, delta);
-
-    std::vector<double> *parameters = swarm_01.get_solution();
+    double delta = 0.7;
+
+    ParticleSwarm swarm_01(
+            &domain_bounds,
+            c1,
+            c2,
+            w,
+            gamma,
+            epsilon,
+            delta,
+            n_particles,
+            iter_max
+    );
+    swarm_01.optimize( mse );
+
+    std::vector<double> *parameters = swarm_01.get_parameters();
     net.copy_parameter_space(parameters);
 
     printf("w1 = %10.7f\n", parameters->at( 0 ));
diff --git a/src/examples/net_test_3.cpp b/src/examples/net_test_3.cpp
index 9354fe799c74e89c380e2b23475c7564864d7264..8219524159c6602acfead0d85b58d53578c79e3b 100644
--- a/src/examples/net_test_3.cpp
+++ b/src/examples/net_test_3.cpp
@@ -82,6 +82,7 @@ int main() {
 
     /* CONSTRUCTION OF SUBNETWORKS */
     //TODO subnetworks retain the number of weights, could be optimized to include only the used weights
+    //TODO this example is not working properly, subnet method is not implemented
     std::vector<size_t> subnet_01_input_neurons, subnet_01_output_neurons;
     std::vector<size_t> subnet_02_input_neurons, subnet_02_output_neurons;
 
@@ -104,8 +105,18 @@ int main() {
         mse_sum.add_error_function( &mse_02 );
 
         /* TRAINING METHOD SETUP */
-        std::vector<double> domain_bounds = {-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0,-10.0, 10.0, -10.0, 10.0};
-        ParticleSwarm swarm_01(&mse_sum, &domain_bounds);
+        std::vector<double> domain_bounds(2 * (net.get_n_weights() + net.get_n_biases()));
+
+        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.7;
+        size_t n_particles = 50;
+        size_t iter_max = 1000;
 
         /* if the maximal velocity from the previous step is less than 'gamma' times the current maximal velocity, then one
          * terminating criterion is met */
@@ -114,8 +125,20 @@ int main() {
         /* 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;
-        swarm_01.optimize(gamma, epsilon, delta);
+        double delta = 0.7;
+
+        ParticleSwarm swarm_01(
+                &domain_bounds,
+                c1,
+                c2,
+                w,
+                gamma,
+                epsilon,
+                delta,
+                n_particles,
+                iter_max
+        );
+        swarm_01.optimize( mse_sum );
 
 
     }
diff --git a/src/examples/net_test_harmonic_oscilator.cpp b/src/examples/net_test_harmonic_oscilator.cpp
index 17fddd014d734cd2e8d419b34ee781817e06f8cb..6eecc425d72ebb638b2051defaa52a134c371698 100644
--- a/src/examples/net_test_harmonic_oscilator.cpp
+++ b/src/examples/net_test_harmonic_oscilator.cpp
@@ -61,24 +61,19 @@ void test_harmonic_oscilator_fixed_E(double EE, double accuracy, size_t n_inner_
     /* 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;
+    /* TRAINING METHOD SETUP */
+    size_t total_dim = solver.get_solution( alpha_0 )->get_n_biases() + solver.get_solution( alpha_0 )->get_n_weights();
+    std::vector<double> domain_bounds( 2 * total_dim );
+
+    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;
 
-    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;
@@ -86,8 +81,20 @@ void test_harmonic_oscilator_fixed_E(double EE, double accuracy, size_t n_inner_
     /* 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 );
+    double delta = 0.7;
+
+    ParticleSwarm swarm(
+            &domain_bounds,
+            c1,
+            c2,
+            w,
+            gamma,
+            epsilon,
+            delta,
+            n_particles,
+            max_iters
+    );
+    solver.solve( swarm );
 
     NeuralNetwork *solution = solver.get_solution( alpha_0 );
     std::vector<double> parameters(total_dim);//w1, a1, b1, w2, a2, b2, ... , wm, am, bm
diff --git a/src/examples/net_test_ode_1.cpp b/src/examples/net_test_ode_1.cpp
index f7019f1dd4aaaddef52b5785cd75651307dfbc6c..1e857ba22bf9d1d3b5137613acd11b8ea6a3ab36 100644
--- a/src/examples/net_test_ode_1.cpp
+++ b/src/examples/net_test_ode_1.cpp
@@ -315,168 +315,168 @@ double calculate_gradient( std::vector<double> &data_points, size_t n_inner_neur
 }
 
 
-void test_analytical_gradient_y(std::vector<double> &guess, double accuracy, size_t n_inner_neurons, size_t train_size, double d1_s, double d1_e,
-                                size_t test_size, double ts, double te) {
-
-    std::cout << "Finding a solution via a Gradient Descent method with adaptive step-length" << std::endl;
-    std::cout << "********************************************************************************************************************************************" <<std::endl;
-
-    /* SETUP OF THE TRAINING DATA */
-    std::vector<double> inp, out;
-
-    double frac, alpha, x;
-    double grad_norm = accuracy * 10.0, mem, ai, bi, wi, error, derror, approx, xj, gamma, total_error, sk, sy, sx, sg, beta;
-    double grad_norm_prev = grad_norm;
-    size_t i, j, iter_idx = 0;
-
-    /* TRAIN DATA FOR THE GOVERNING DE */
-    std::vector<double> data_points(train_size);
-
-    /* ISOTROPIC TRAIN SET */
-    frac = (d1_e - d1_s) / (train_size - 1);
-    for(unsigned int i = 0; i < train_size; ++i){
-        data_points[i] = frac * i;
-    }
-
-//    /* CHEBYSCHEV TRAIN SET */
-//    alpha = PI / (train_size );
-//    frac = 0.5 * (d1_e - d1_s);
-//    for(i = 0; i < train_size; ++i){
-//        x = (std::cos(PI - alpha * i) + 1.0) * frac + d1_s;
-//        data_points[i] = x;
+//void test_analytical_gradient_y(std::vector<double> &guess, double accuracy, size_t n_inner_neurons, size_t train_size, double d1_s, double d1_e,
+//                                size_t test_size, double ts, double te) {
+//
+//    std::cout << "Finding a solution via a Gradient Descent method with adaptive step-length" << std::endl;
+//    std::cout << "********************************************************************************************************************************************" <<std::endl;
+//
+//    /* SETUP OF THE TRAINING DATA */
+//    std::vector<double> inp, out;
+//
+//    double frac, alpha, x;
+//    double grad_norm = accuracy * 10.0, mem, ai, bi, wi, error, derror, approx, xj, gamma, total_error, sk, sy, sx, sg, beta;
+//    double grad_norm_prev = grad_norm;
+//    size_t i, j, iter_idx = 0;
+//
+//    /* TRAIN DATA FOR THE GOVERNING DE */
+//    std::vector<double> data_points(train_size);
+//
+//    /* ISOTROPIC TRAIN SET */
+//    frac = (d1_e - d1_s) / (train_size - 1);
+//    for(unsigned int i = 0; i < train_size; ++i){
+//        data_points[i] = frac * i;
 //    }
-
-//    DataSet ds(0.0, 4.0, train_size, 0.0);
-
-    std::vector<double> *gradient_current = new std::vector<double>(3 * n_inner_neurons);
-    std::vector<double> *gradient_prev = new std::vector<double>(3 * n_inner_neurons);
-    std::vector<double> *params_current = new std::vector<double>(guess);
-    std::vector<double> *params_prev = new std::vector<double>(guess);
-    std::vector<double> *conjugate_direction_current = new std::vector<double>(3 * n_inner_neurons);
-    std::vector<double> *conjugate_direction_prev = new std::vector<double>(3 * n_inner_neurons);
-
-    std::vector<double> *ptr_mem;
-
-
-    std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
-    std::fill(gradient_prev->begin(), gradient_prev->end(), 0.0);
-    std::fill(conjugate_direction_current->begin(), conjugate_direction_current->end(), 0.0);
-    std::fill(conjugate_direction_prev->begin(), conjugate_direction_prev->end(), 0.0);
-
-
-
+//
+////    /* CHEBYSCHEV TRAIN SET */
+////    alpha = PI / (train_size );
+////    frac = 0.5 * (d1_e - d1_s);
+////    for(i = 0; i < train_size; ++i){
+////        x = (std::cos(PI - alpha * i) + 1.0) * frac + d1_s;
+////        data_points[i] = x;
+////    }
+//
+////    DataSet ds(0.0, 4.0, train_size, 0.0);
+//
+//    std::vector<double> *gradient_current = new std::vector<double>(3 * n_inner_neurons);
+//    std::vector<double> *gradient_prev = new std::vector<double>(3 * n_inner_neurons);
+//    std::vector<double> *params_current = new std::vector<double>(guess);
+//    std::vector<double> *params_prev = new std::vector<double>(guess);
+//    std::vector<double> *conjugate_direction_current = new std::vector<double>(3 * n_inner_neurons);
+//    std::vector<double> *conjugate_direction_prev = new std::vector<double>(3 * n_inner_neurons);
+//
+//    std::vector<double> *ptr_mem;
+//
+//
+//    std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
+//    std::fill(gradient_prev->begin(), gradient_prev->end(), 0.0);
+//    std::fill(conjugate_direction_current->begin(), conjugate_direction_current->end(), 0.0);
+//    std::fill(conjugate_direction_prev->begin(), conjugate_direction_prev->end(), 0.0);
+//
+//
+//
+////    for (i = 0; i < n_inner_neurons; ++i) {
+////        wi = (*params_current)[3 * i];
+////        ai = (*params_current)[3 * i + 1];
+////        bi = (*params_current)[3 * i + 2];
+////
+////        printf("Path %3d. w = %15.8f, b = %15.8f, a = %15.8f\n", (int)(i + 1), wi, bi, ai);
+////    }
+//
+//    gamma = 1.0;
+//    double prev_val, val = 0.0, c = 2.0;
+//    while( grad_norm > accuracy) {
+//        iter_idx++;
+//        prev_val = val;
+//        grad_norm_prev = grad_norm;
+//
+//        /* reset of the current gradient */
+//        std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
+//        val = calculate_gradient( data_points, n_inner_neurons, params_current, gradient_current );
+//
+//        grad_norm = 0.0;
+//        for(auto v: *gradient_current){
+//            grad_norm += v * v;
+//        }
+//        grad_norm = std::sqrt(grad_norm);
+//
+//        /* Update of the parameters */
+//        /* step length calculation */
+//        if(iter_idx < 10 || iter_idx % 100 == 0){
+//            /* fixed step length */
+//            gamma = 0.1 * accuracy;
+//        }
+//        else{
+//
+//
+//
+//            /* norm of the gradient calculation */
+//
+//            sk = 0.0;
+//            for(i = 0; i < gradient_current->size(); ++i){
+//                sx = (*gradient_current)[i] - (*gradient_prev)[i];
+//                sk += sx * sx;
+//            }
+//            sk = std::sqrt(sk);
+//
+//            /* angle between two consecutive gradients */
+//            double beta = 0.0, sx = 0.0;
+//            for(i = 0; i < gradient_current->size(); ++i){
+//                sx += (gradient_current->at( i ) * gradient_prev->at( i ));
+//            }
+//            sx /= grad_norm * grad_norm_prev;
+//            beta = std::sqrt(std::acos( sx ) / PI);
+//
+//
+////            eval_step_size_simple( gamma, val, prev_val, sk, grad_norm, grad_norm_prev );
+//            eval_step_size_mk( gamma, beta, c, grad_norm_prev, grad_norm, val, prev_val );
+//        }
+//
+//
+//
+//
+//
+//
+//        for(i = 0; i < gradient_current->size(); ++i){
+//            (*params_prev)[i] = (*params_current)[i] - gamma * (*gradient_current)[i];
+//        }
+//
+//        /* switcheroo */
+//        ptr_mem = gradient_prev;
+//        gradient_prev = gradient_current;
+//        gradient_current = ptr_mem;
+//
+//        ptr_mem = params_prev;
+//        params_prev = params_current;
+//        params_current = ptr_mem;
+//
+//
+//        if(iter_idx % 1 == 0){
+//            printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r", (int)iter_idx, gamma, c, grad_norm, val);
+//            std::cout.flush();
+//        }
+//    }
+//    printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r\n", (int)iter_idx, gamma, c, grad_norm, val);
+//    std::cout.flush();
+//
 //    for (i = 0; i < n_inner_neurons; ++i) {
 //        wi = (*params_current)[3 * i];
 //        ai = (*params_current)[3 * i + 1];
 //        bi = (*params_current)[3 * i + 2];
 //
-//        printf("Path %3d. w = %15.8f, b = %15.8f, a = %15.8f\n", (int)(i + 1), wi, bi, ai);
+//        printf("Path %3d. w%d = %15.8f, b%d = %15.8f, a%d = %15.8f\n", (int)(i + 1), (int)(i + 1), wi, (int)(i + 1), bi, (int)(i + 1), ai);
 //    }
-
-    gamma = 1.0;
-    double prev_val, val = 0.0, c = 2.0;
-    while( grad_norm > accuracy) {
-        iter_idx++;
-        prev_val = val;
-        grad_norm_prev = grad_norm;
-
-        /* reset of the current gradient */
-        std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
-        val = calculate_gradient( data_points, n_inner_neurons, params_current, gradient_current );
-
-        grad_norm = 0.0;
-        for(auto v: *gradient_current){
-            grad_norm += v * v;
-        }
-        grad_norm = std::sqrt(grad_norm);
-
-        /* Update of the parameters */
-        /* step length calculation */
-        if(iter_idx < 10 || iter_idx % 100 == 0){
-            /* fixed step length */
-            gamma = 0.1 * accuracy;
-        }
-        else{
-
-
-
-            /* norm of the gradient calculation */
-
-            sk = 0.0;
-            for(i = 0; i < gradient_current->size(); ++i){
-                sx = (*gradient_current)[i] - (*gradient_prev)[i];
-                sk += sx * sx;
-            }
-            sk = std::sqrt(sk);
-
-            /* angle between two consecutive gradients */
-            double beta = 0.0, sx = 0.0;
-            for(i = 0; i < gradient_current->size(); ++i){
-                sx += (gradient_current->at( i ) * gradient_prev->at( i ));
-            }
-            sx /= grad_norm * grad_norm_prev;
-            beta = std::sqrt(std::acos( sx ) / PI);
-
-
-//            eval_step_size_simple( gamma, val, prev_val, sk, grad_norm, grad_norm_prev );
-            eval_step_size_mk( gamma, beta, c, grad_norm_prev, grad_norm, val, prev_val );
-        }
-
-
-
-
-
-
-        for(i = 0; i < gradient_current->size(); ++i){
-            (*params_prev)[i] = (*params_current)[i] - gamma * (*gradient_current)[i];
-        }
-
-        /* switcheroo */
-        ptr_mem = gradient_prev;
-        gradient_prev = gradient_current;
-        gradient_current = ptr_mem;
-
-        ptr_mem = params_prev;
-        params_prev = params_current;
-        params_current = ptr_mem;
-
-
-        if(iter_idx % 1 == 0){
-            printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r", (int)iter_idx, gamma, c, grad_norm, val);
-            std::cout.flush();
-        }
-    }
-    printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r\n", (int)iter_idx, gamma, c, grad_norm, val);
-    std::cout.flush();
-
-    for (i = 0; i < n_inner_neurons; ++i) {
-        wi = (*params_current)[3 * i];
-        ai = (*params_current)[3 * i + 1];
-        bi = (*params_current)[3 * i + 2];
-
-        printf("Path %3d. w%d = %15.8f, b%d = %15.8f, a%d = %15.8f\n", (int)(i + 1), (int)(i + 1), wi, (int)(i + 1), bi, (int)(i + 1), ai);
-    }
-    std::cout << "********************************************************************************************************************************************" <<std::endl;
-
-
-//    data_export_gradient(params_current);
-//    if(total_error < 1e-3 || true){
-//        /* ISOTROPIC TEST SET */
-//        frac = (te - ts) / (test_size - 1);
-//        for(j = 0; j < test_size; ++j){
-//            xj = frac * j + ts;
+//    std::cout << "********************************************************************************************************************************************" <<std::endl;
 //
-//            std::cout << j + 1 << " " << xj << " " << eval_f(xj) << " " << eval_approx_f(xj, n_inner_neurons, *params_current) << " " << eval_df(xj) << " " << eval_approx_df(xj, n_inner_neurons, *params_current) << " " << eval_ddf(xj) << " " << eval_approx_ddf(xj, n_inner_neurons, *params_current) << std::endl;
-//        }
-//    }
-
-    delete gradient_current;
-    delete gradient_prev;
-    delete params_current;
-    delete params_prev;
-    delete conjugate_direction_current;
-    delete conjugate_direction_prev;
-}
+//
+////    data_export_gradient(params_current);
+////    if(total_error < 1e-3 || true){
+////        /* ISOTROPIC TEST SET */
+////        frac = (te - ts) / (test_size - 1);
+////        for(j = 0; j < test_size; ++j){
+////            xj = frac * j + ts;
+////
+////            std::cout << j + 1 << " " << xj << " " << eval_f(xj) << " " << eval_approx_f(xj, n_inner_neurons, *params_current) << " " << eval_df(xj) << " " << eval_approx_df(xj, n_inner_neurons, *params_current) << " " << eval_ddf(xj) << " " << eval_approx_ddf(xj, n_inner_neurons, *params_current) << std::endl;
+////        }
+////    }
+//
+//    delete gradient_current;
+//    delete gradient_prev;
+//    delete params_current;
+//    delete params_prev;
+//    delete conjugate_direction_current;
+//    delete conjugate_direction_prev;
+//}
 
 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){
 
@@ -606,17 +606,18 @@ void test_ode(double accuracy, size_t n_inner_neurons, size_t train_size, double
 //        printf("\tRepresentation test %6d, error of eq1: %10.8f, error of eq2: %10.8f, error of eq3: %10.8f, total error: %10.8f\n", (int)testi, std::sqrt(test_error_eq1), std::sqrt(test_error_eq2), std::sqrt(test_error_eq3), (total_error_analytical - total_error_de_solver) * (total_error_analytical - total_error_de_solver));
 //    }
 
-    /* PARTICLE SWARM TRAINING METHOD SETUP */
+    /* TRAINING METHOD SETUP */
+    std::vector<double> domain_bounds(2 * (solver_01.get_solution( alpha_0 )->get_n_biases() + solver_01.get_solution( alpha_0 )->get_n_weights()));
 
-    //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;
+    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;
 
-    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;
@@ -624,8 +625,20 @@ void test_ode(double accuracy, size_t n_inner_neurons, size_t train_size, double
     /* 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_01.solve_via_particle_swarm( &domain_bounds, c1, c2, w, n_particles, max_iters, gamma, epsilon, delta );
+    double delta = 0.7;
+
+    ParticleSwarm swarm_01(
+            &domain_bounds,
+            c1,
+            c2,
+            w,
+            gamma,
+            epsilon,
+            delta,
+            n_particles,
+            max_iters
+    );
+    solver_01.solve( swarm_01 );
 
     NeuralNetwork *solution = solver_01.get_solution( alpha_0 );
     NeuralNetwork *solution_d = solver_01.get_solution( alpha_1 );
diff --git a/src/examples/net_test_pde_1.cpp b/src/examples/net_test_pde_1.cpp
index f27a9b23497deb74aea9f5fbe1eda9c241197a12..2bce443169faa603c1a9505dc2b3a78366049dc5 100644
--- a/src/examples/net_test_pde_1.cpp
+++ b/src/examples/net_test_pde_1.cpp
@@ -60,24 +60,24 @@ double eval_approx_yt(double x, double t, size_t n_inner_neurons, std::vector<do
     }
     return value;
 }
-//yx(x, t) = ...  (ai * wxi * e^(bi - t * wti - wxi * x))/(e^(bi - t * wti - wxi * x) + 1)^2
-double eval_approx_yx(double x, double t, size_t n_inner_neurons, std::vector<double> &parameters){
-    double value= 0.0, wxi, wti, ai, bi, ei, ei1;
-
-    for(size_t i = 0; i < n_inner_neurons; ++i){
-
-        wxi = parameters[4 * i + 0];
-        wti = parameters[4 * i + 1];
-        ai  = parameters[4 * i + 2];
-        bi  = parameters[4 * i + 3];
-
-        ei = std::pow(E, bi - wxi * x - wti * t);
-        ei1 = ei + 1.0;
-
-        value += (ai * wxi * ei1)/(ei1 * ei1);
-    }
-    return value;
-}
+////yx(x, t) = ...  (ai * wxi * e^(bi - t * wti - wxi * x))/(e^(bi - t * wti - wxi * x) + 1)^2
+//double eval_approx_yx(double x, double t, size_t n_inner_neurons, std::vector<double> &parameters){
+//    double value= 0.0, wxi, wti, ai, bi, ei, ei1;
+//
+//    for(size_t i = 0; i < n_inner_neurons; ++i){
+//
+//        wxi = parameters[4 * i + 0];
+//        wti = parameters[4 * i + 1];
+//        ai  = parameters[4 * i + 2];
+//        bi  = parameters[4 * i + 3];
+//
+//        ei = std::pow(E, bi - wxi * x - wti * t);
+//        ei1 = ei + 1.0;
+//
+//        value += (ai * wxi * ei1)/(ei1 * ei1);
+//    }
+//    return value;
+//}
 //yxx(x, t) = ...  (ai * wxi * e^(bi - t * wti - wxi * x))/(e^(bi - t * wti - wxi * x) + 1)^2
 double eval_approx_yxx(double x, double t, size_t n_inner_neurons, std::vector<double> &parameters){
     double value= 0.0, wxi, wti, ai, bi, ei, ei1;
@@ -95,539 +95,539 @@ double eval_approx_yxx(double x, double t, size_t n_inner_neurons, std::vector<d
     }
     return value;
 }
-
-double eval_approx_da_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, bi, ei, ei1;
-
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return 1.0 / ei1;
-}
-
-double eval_approx_dwx_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return  (ai * x * ei)/(ei1 * ei1);
-}
-
-double eval_approx_dwt_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return  (ai * t * ei)/(ei1 * ei1);
-}
-
-double eval_approx_db_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, bi, ei, ai, ei1;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return -(ai * ei)/(ei1 * ei1);
-}
-
-double eval_approx_da_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, bi, ei, ei1;
-
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return (wti * ei)/(ei1 * ei1);
-}
-
-double eval_approx_dwx_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return (2 * ai * wti * x * ei * ei)/(ei1 * ei1 * ei1) - (ai * wti * x * ei)/(ei1 * ei1);
-}
-
-double eval_approx_dwt_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return  -(ai * t * wti * ei) / (ei1 * ei1) + (2 * ai * t * wti * ei * ei)/(ei1 * ei1 * ei1) + (ai * ei)/(ei1 * ei1);
-}
-
-double eval_approx_db_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ei1 = ei + 1.0;
-
-    return (ai * wti * ei) / (ei1 * ei1) - (2 * ai * wti * ei * ei) / (ei1 * ei1 * ei1);
-}
-
-double eval_approx_da_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ebp= std::pow(E, bi + wxi * x + wti * t);
-    eb = std::pow(E, bi);
-    etx = std::pow(E, wxi * x + wti * t);
-    ei1 = eb + etx;
-
-    return -(wxi * wxi * ebp * (etx - eb))/(ei1 * ei1 * ei1);
-}
-
-double eval_approx_dwx_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ebp= std::pow(E, bi + wxi * x + wti * t);
-    eb = std::pow(E, bi);
-    etx = std::pow(E, wxi * x + wti * t);
-    ei1 = eb + etx;
-
-    return (ai * wxi * wxi * x * ei) / ((ei + 1) * (ei + 1)) - (6 * ai * wxi * wxi * x * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1)) + (6 * ai * wxi *wxi * x * ei * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1) * (ei + 1)) - (2 * ai * wxi * ei) / ((ei + 1) * (ei + 1)) + (4 * ai * wxi * ei * ei)/((ei + 1) * (ei + 1) * (ei + 1));
-}
-
-double eval_approx_dwt_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ebp= std::pow(E, bi + wxi * x + wti * t);
-    eb = std::pow(E, bi);
-    etx = std::pow(E, wxi * x + wti * t);
-    ei1 = eb + etx;
-
-    return (ai * t * wxi * wxi * ei) / ((ei + 1) * (ei + 1)) - (6 * ai * t * wxi * wxi * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1)) + (6 * ai * t * wxi * wxi * ei * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1) * (ei + 1));
-}
-
-double eval_approx_db_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
-    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
-    wxi = parameters[4 * neuron_idx + 0];
-    wti = parameters[4 * neuron_idx + 1];
-    ai =  parameters[4 * neuron_idx + 2];
-    bi =  parameters[4 * neuron_idx + 3];
-
-    ei = std::pow(E, bi - wxi * x - wti * t);
-    ebp= std::pow(E, bi + wxi * x + wti * t);
-    eb = std::pow(E, bi);
-    etx = std::pow(E, wxi * x + wti * t);
-    ei1 = eb + etx;
-
-    return (ai * wxi * wxi * eb * ebp) / (ei1 * ei1 * ei1) - (ai * wxi * wxi * ebp * (etx - eb)) / (ei1 * ei1 * ei1) + (3 * ai * wxi * wxi * eb * ebp * (etx - eb)) / (ei1 * ei1 * ei1 * ei1);
-}
-
-void eval_step_size_simple(double &gamma, double val, double prev_val, double sk, double grad_norm, double grad_norm_prev){
-
-    if(val > prev_val){
-        gamma *= 0.99999;
-    }
-
-    if(sk <= 1e-3 || grad_norm < grad_norm_prev){
-        /* movement on a line */
-        /* new slope is less steep, speed up */
-        gamma *= 1.0005;
-    }
-    else if(grad_norm > grad_norm_prev){
-        /* new slope is more steep, slow down*/
-        gamma /= 1.0005;
-    }
-    else{
-        gamma /= 1.005;
-    }
-//        gamma *= 0.999999;
-}
-
-void eval_step_size_mk( double &gamma, double beta, double &c, double grad_norm_prev, double grad_norm, double fi, double fim ){
-
-    if( fi > fim )
-    {
-        c /= 1.0000005;
-    }
-    else if( fi < fim )
-    {
-        c *= 1.0000005;
-    }
-
-    gamma *= std::pow( c, 1.0 - 2.0 * beta) * std::pow( grad_norm_prev / grad_norm, 1.0 / c );
-
-
-}
-
-double calculate_gradient( std::vector<double> &data_points, size_t n_inner_neurons, std::vector<double> *parameters, std::vector<double> *gradient ){
-    size_t i, j, k;
-    double x, t, mem, derror, total_error, approx;
-
-    size_t train_size = data_points.size();
-
-    /* error of boundary condition: y(0, t) = sin(t) => e1 = 1/n * (sin(t) - y(0, t))^2 */
-    for(i = 0; i < train_size; ++i){
-
-        t = data_points[i];
-        mem = (std::sin(t) - eval_approx_y(0.0, t, n_inner_neurons, *parameters));
-        derror = 2.0 * mem / train_size;
-
-        for(j = 0; j < n_inner_neurons; ++j){
-            (*gradient)[4 * j + 0] -= derror * eval_approx_dwx_y(0, t, j, *parameters);
-            (*gradient)[4 * j + 1] -= derror * eval_approx_dwt_y(0, t, j, *parameters);
-            (*gradient)[4 * j + 2] -= derror *  eval_approx_da_y(0, t, j, *parameters);
-            (*gradient)[4 * j + 3] -= derror *  eval_approx_db_y(0, t, j, *parameters);
-        }
-
-        total_error += mem * mem / train_size;
-    }
-
-
-
-    /* error boundary condition: y(x, 0) = e^(-(0.5)^(0.5)x) * sin(-(0.5)^(0.5)x) => e2 = 1/n * (e^(-(0.5)^(0.5)x) * sin(-(0.5)^(0.5)x) - y(x, 0))^2 */
-    for(i = 0; i < train_size; ++i){
-
-        x = data_points[i];
-        mem = (std::pow(E, -0.707106781 * x) * std::sin( -0.707106781 * x ) - eval_approx_y(x, 0.0, n_inner_neurons, *parameters));
-        derror = 2.0 * mem / train_size;
-
-        for(j = 0; j < n_inner_neurons; ++j){
-            (*gradient)[4 * j + 0] -= derror * eval_approx_dwx_y(x, 0, j, *parameters);
-            (*gradient)[4 * j + 1] -= derror * eval_approx_dwt_y(x, 0, j, *parameters);
-            (*gradient)[4 * j + 2] -= derror *  eval_approx_da_y(x, 0, j, *parameters);
-            (*gradient)[4 * j + 3] -= derror *  eval_approx_db_y(x, 0, j, *parameters);
-        }
-
-        total_error += mem * mem / train_size;
-    }
-
-    /* error of the governing equation: y_xx - y_t = 0 => e3 = 1/n^2 * (0 - y_xx + y_t)^2 */
-    for(i = 0; i < data_points.size(); ++i){
-        x = data_points[i];
-        for(j = 0; j < data_points.size(); ++j){
-            t = data_points[j];
-
-            approx = eval_approx_yxx(x, t, n_inner_neurons, *parameters) - eval_approx_yt(x, t, n_inner_neurons, *parameters);
-            mem = 0.0 - approx;
-            derror = 2.0 * mem / (train_size * train_size);
-            for(k = 0; k < n_inner_neurons; ++k){
-                (*gradient)[4 * k + 0] -= derror * (eval_approx_dwx_yxx(x, t, k, *parameters) - eval_approx_dwx_yt(x, t, k, *parameters));
-                (*gradient)[4 * k + 1] -= derror * (eval_approx_dwt_yxx(x, t, k, *parameters) - eval_approx_dwt_yt(x, t, k, *parameters));
-                (*gradient)[4 * k + 2] -= derror * ( eval_approx_da_yxx(x, t, k, *parameters) -  eval_approx_da_yt(x, t, k, *parameters));
-                (*gradient)[4 * k + 3] -= derror * ( eval_approx_db_yxx(x, t, k, *parameters) -  eval_approx_db_yt(x, t, k, *parameters));
-            }
-            total_error += mem * mem / (train_size * train_size);
-        }
-    }
-    return total_error;
-}
-
-void solve_example_gradient(std::vector<double> &guess, double accuracy, size_t n_inner_neurons, size_t train_size, double ds, double de, size_t n_test_points, double ts, double te){
-    std::cout << "Finding a solution via a Gradient Descent method with adaptive step-length" << std::endl;
-    std::cout << "********************************************************************************************************************************************" <<std::endl;
-    /* SETUP OF THE TRAINING DATA */
-    std::vector<double> inp, out;
-
-    double frac, alpha, x;
-    double grad_norm = accuracy * 10.0, mem, ai, bi, wxi, wti, error, derror, approx, t, gamma, total_error, sk, sy, sx, sg, beta;
-    double grad_norm_prev = grad_norm;
-    size_t i, j, k, iter_idx = 0;
-
-    /* TRAIN DATA FOR THE GOVERNING DE */
-    std::vector<double> data_points(train_size);
-
-    /* ISOTROPIC TRAIN SET */
-    frac = (de - ds) / (train_size - 1);
-    for(i = 0; i < train_size; ++i){
-        data_points[i] = frac * i + ds;
-//        std::cout << data_points[i] << std::endl;
-    }
-
-//    /* CHEBYSCHEV TRAIN SET */
-//    alpha = PI / (train_size );
-//    frac = 0.5 * (de - ds);
+//
+//double eval_approx_da_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, bi, ei, ei1;
+//
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return 1.0 / ei1;
+//}
+//
+//double eval_approx_dwx_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return  (ai * x * ei)/(ei1 * ei1);
+//}
+//
+//double eval_approx_dwt_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return  (ai * t * ei)/(ei1 * ei1);
+//}
+//
+//double eval_approx_db_y(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, bi, ei, ai, ei1;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return -(ai * ei)/(ei1 * ei1);
+//}
+//
+//double eval_approx_da_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, bi, ei, ei1;
+//
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return (wti * ei)/(ei1 * ei1);
+//}
+//
+//double eval_approx_dwx_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return (2 * ai * wti * x * ei * ei)/(ei1 * ei1 * ei1) - (ai * wti * x * ei)/(ei1 * ei1);
+//}
+//
+//double eval_approx_dwt_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return  -(ai * t * wti * ei) / (ei1 * ei1) + (2 * ai * t * wti * ei * ei)/(ei1 * ei1 * ei1) + (ai * ei)/(ei1 * ei1);
+//}
+//
+//double eval_approx_db_yt(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ei1 = ei + 1.0;
+//
+//    return (ai * wti * ei) / (ei1 * ei1) - (2 * ai * wti * ei * ei) / (ei1 * ei1 * ei1);
+//}
+//
+//double eval_approx_da_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ebp= std::pow(E, bi + wxi * x + wti * t);
+//    eb = std::pow(E, bi);
+//    etx = std::pow(E, wxi * x + wti * t);
+//    ei1 = eb + etx;
+//
+//    return -(wxi * wxi * ebp * (etx - eb))/(ei1 * ei1 * ei1);
+//}
+//
+//double eval_approx_dwx_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ebp= std::pow(E, bi + wxi * x + wti * t);
+//    eb = std::pow(E, bi);
+//    etx = std::pow(E, wxi * x + wti * t);
+//    ei1 = eb + etx;
+//
+//    return (ai * wxi * wxi * x * ei) / ((ei + 1) * (ei + 1)) - (6 * ai * wxi * wxi * x * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1)) + (6 * ai * wxi *wxi * x * ei * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1) * (ei + 1)) - (2 * ai * wxi * ei) / ((ei + 1) * (ei + 1)) + (4 * ai * wxi * ei * ei)/((ei + 1) * (ei + 1) * (ei + 1));
+//}
+//
+//double eval_approx_dwt_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ebp= std::pow(E, bi + wxi * x + wti * t);
+//    eb = std::pow(E, bi);
+//    etx = std::pow(E, wxi * x + wti * t);
+//    ei1 = eb + etx;
+//
+//    return (ai * t * wxi * wxi * ei) / ((ei + 1) * (ei + 1)) - (6 * ai * t * wxi * wxi * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1)) + (6 * ai * t * wxi * wxi * ei * ei * ei) / ((ei + 1) * (ei + 1) * (ei + 1) * (ei + 1));
+//}
+//
+//double eval_approx_db_yxx(double x, double t, size_t neuron_idx, std::vector<double> &parameters){
+//    double wxi, wti, ai, bi, ei, ei1, ebp, eb, etx;
+//    wxi = parameters[4 * neuron_idx + 0];
+//    wti = parameters[4 * neuron_idx + 1];
+//    ai =  parameters[4 * neuron_idx + 2];
+//    bi =  parameters[4 * neuron_idx + 3];
+//
+//    ei = std::pow(E, bi - wxi * x - wti * t);
+//    ebp= std::pow(E, bi + wxi * x + wti * t);
+//    eb = std::pow(E, bi);
+//    etx = std::pow(E, wxi * x + wti * t);
+//    ei1 = eb + etx;
+//
+//    return (ai * wxi * wxi * eb * ebp) / (ei1 * ei1 * ei1) - (ai * wxi * wxi * ebp * (etx - eb)) / (ei1 * ei1 * ei1) + (3 * ai * wxi * wxi * eb * ebp * (etx - eb)) / (ei1 * ei1 * ei1 * ei1);
+//}
+//
+//void eval_step_size_simple(double &gamma, double val, double prev_val, double sk, double grad_norm, double grad_norm_prev){
+//
+//    if(val > prev_val){
+//        gamma *= 0.99999;
+//    }
+//
+//    if(sk <= 1e-3 || grad_norm < grad_norm_prev){
+//        /* movement on a line */
+//        /* new slope is less steep, speed up */
+//        gamma *= 1.0005;
+//    }
+//    else if(grad_norm > grad_norm_prev){
+//        /* new slope is more steep, slow down*/
+//        gamma /= 1.0005;
+//    }
+//    else{
+//        gamma /= 1.005;
+//    }
+////        gamma *= 0.999999;
+//}
+//
+//void eval_step_size_mk( double &gamma, double beta, double &c, double grad_norm_prev, double grad_norm, double fi, double fim ){
+//
+//    if( fi > fim )
+//    {
+//        c /= 1.0000005;
+//    }
+//    else if( fi < fim )
+//    {
+//        c *= 1.0000005;
+//    }
+//
+//    gamma *= std::pow( c, 1.0 - 2.0 * beta) * std::pow( grad_norm_prev / grad_norm, 1.0 / c );
+//
+//
+//}
+//
+//double calculate_gradient( std::vector<double> &data_points, size_t n_inner_neurons, std::vector<double> *parameters, std::vector<double> *gradient ){
+//    size_t i, j, k;
+//    double x, t, mem, derror, total_error, approx;
+//
+//    size_t train_size = data_points.size();
+//
+//    /* error of boundary condition: y(0, t) = sin(t) => e1 = 1/n * (sin(t) - y(0, t))^2 */
 //    for(i = 0; i < train_size; ++i){
-//        x = (std::cos(PI - alpha * i) + 1.0) * frac + ds;
-//        data_points[i] = x;
+//
+//        t = data_points[i];
+//        mem = (std::sin(t) - eval_approx_y(0.0, t, n_inner_neurons, *parameters));
+//        derror = 2.0 * mem / train_size;
+//
+//        for(j = 0; j < n_inner_neurons; ++j){
+//            (*gradient)[4 * j + 0] -= derror * eval_approx_dwx_y(0, t, j, *parameters);
+//            (*gradient)[4 * j + 1] -= derror * eval_approx_dwt_y(0, t, j, *parameters);
+//            (*gradient)[4 * j + 2] -= derror *  eval_approx_da_y(0, t, j, *parameters);
+//            (*gradient)[4 * j + 3] -= derror *  eval_approx_db_y(0, t, j, *parameters);
+//        }
+//
+//        total_error += mem * mem / train_size;
 //    }
-
-    std::vector<double> *gradient_current = new std::vector<double>(4 * n_inner_neurons);
-    std::vector<double> *gradient_prev = new std::vector<double>(4 * n_inner_neurons);
-    std::vector<double> *params_current = new std::vector<double>(guess);
-    std::vector<double> *params_prev = new std::vector<double>(guess);
-
-    std::vector<double> *ptr_mem;
-
-    std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
-    std::fill(gradient_prev->begin(), gradient_prev->end(), 0.0);
-
-//    for (i = 0; i < n_inner_neurons; ++i) {
-//        wxi = (*params_current)[4 * i + 0];
-//        wti = (*params_current)[4 * i + 1];
-//        ai  = (*params_current)[4 * i + 2];
-//        bi  = (*params_current)[4 * i + 3];
 //
-//        printf("Path %3d. wx = %15.8f, wt = %15.8f, b = %15.8f, a = %15.8f\n", (int)(i + 1), wxi, wti, bi, ai);
+//
+//
+//    /* error boundary condition: y(x, 0) = e^(-(0.5)^(0.5)x) * sin(-(0.5)^(0.5)x) => e2 = 1/n * (e^(-(0.5)^(0.5)x) * sin(-(0.5)^(0.5)x) - y(x, 0))^2 */
+//    for(i = 0; i < train_size; ++i){
+//
+//        x = data_points[i];
+//        mem = (std::pow(E, -0.707106781 * x) * std::sin( -0.707106781 * x ) - eval_approx_y(x, 0.0, n_inner_neurons, *parameters));
+//        derror = 2.0 * mem / train_size;
+//
+//        for(j = 0; j < n_inner_neurons; ++j){
+//            (*gradient)[4 * j + 0] -= derror * eval_approx_dwx_y(x, 0, j, *parameters);
+//            (*gradient)[4 * j + 1] -= derror * eval_approx_dwt_y(x, 0, j, *parameters);
+//            (*gradient)[4 * j + 2] -= derror *  eval_approx_da_y(x, 0, j, *parameters);
+//            (*gradient)[4 * j + 3] -= derror *  eval_approx_db_y(x, 0, j, *parameters);
+//        }
+//
+//        total_error += mem * mem / train_size;
 //    }
-
-    gamma = 1.0;
-    double val = 0.0, prev_val, c = 2.0;
-    prev_val = 0.0;
-    while( grad_norm > accuracy) {
-        iter_idx++;
-        prev_val = val;
-        grad_norm_prev = grad_norm;
-
-        /* reset of the current gradient */
-        std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
-        val = calculate_gradient( data_points, n_inner_neurons, params_current, gradient_current );
-
-        grad_norm = 0.0;
-        for(auto v: *gradient_current){
-            grad_norm += v * v;
-        }
-        grad_norm = std::sqrt(grad_norm);
-
-        /* Update of the parameters */
-        /* step length calculation */
-        if(iter_idx < 10 || iter_idx % 100 == 0){
-            /* fixed step length */
-            gamma = 0.1 * accuracy;
-        }
-        else{
-
-
-
-            /* norm of the gradient calculation */
-
-            sk = 0.0;
-            for(i = 0; i < gradient_current->size(); ++i){
-                sx = (*gradient_current)[i] - (*gradient_prev)[i];
-                sk += sx * sx;
-            }
-            sk = std::sqrt(sk);
-
-            /* angle between two consecutive gradients */
-            double beta = 0.0, sx = 0.0;
-            for(i = 0; i < gradient_current->size(); ++i){
-                sx += (gradient_current->at( i ) * gradient_prev->at( i ));
-            }
-            sx /= grad_norm * grad_norm_prev;
-            beta = std::sqrt(std::acos( sx ) / PI);
-
-
-//            eval_step_size_simple( gamma, val, prev_val, sk, grad_norm, grad_norm_prev );
-            eval_step_size_mk( gamma, beta, c, grad_norm_prev, grad_norm, val, prev_val );
-        }
-
-
-
-
-
-
-        for(i = 0; i < gradient_current->size(); ++i){
-            (*params_prev)[i] = (*params_current)[i] - gamma * (*gradient_current)[i];
-        }
-
-        /* switcheroo */
-        ptr_mem = gradient_prev;
-        gradient_prev = gradient_current;
-        gradient_current = ptr_mem;
-
-        ptr_mem = params_prev;
-        params_prev = params_current;
-        params_current = ptr_mem;
-
-
-        if(iter_idx % 1 == 0){
-            printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r", (int)iter_idx, gamma, c, grad_norm, val);
-            std::cout.flush();
-        }
-    }
-    printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r\n", (int)iter_idx, gamma, c, grad_norm, val);
-    std::cout.flush();
-
-    for (i = 0; i < n_inner_neurons; ++i) {
-        wxi = (*params_current)[4 * i + 0];
-        wti = (*params_current)[4 * i + 1];
-        ai  = (*params_current)[4 * i + 2];
-        bi  = (*params_current)[4 * i + 3];
-
-        printf("Path %3d. wx%d = %15.8f, wt%d = %15.8f, b%d = %15.8f, a%d = %15.8f\n", (int)(i + 1), (int)(i + 1), wxi , (int)(i + 1), wti, (int)(i + 1), bi, (int)(i + 1), ai);
-    }
-    std::cout << "********************************************************************************************************************************************" <<std::endl;
-
+//
+//    /* error of the governing equation: y_xx - y_t = 0 => e3 = 1/n^2 * (0 - y_xx + y_t)^2 */
+//    for(i = 0; i < data_points.size(); ++i){
+//        x = data_points[i];
+//        for(j = 0; j < data_points.size(); ++j){
+//            t = data_points[j];
+//
+//            approx = eval_approx_yxx(x, t, n_inner_neurons, *parameters) - eval_approx_yt(x, t, n_inner_neurons, *parameters);
+//            mem = 0.0 - approx;
+//            derror = 2.0 * mem / (train_size * train_size);
+//            for(k = 0; k < n_inner_neurons; ++k){
+//                (*gradient)[4 * k + 0] -= derror * (eval_approx_dwx_yxx(x, t, k, *parameters) - eval_approx_dwx_yt(x, t, k, *parameters));
+//                (*gradient)[4 * k + 1] -= derror * (eval_approx_dwt_yxx(x, t, k, *parameters) - eval_approx_dwt_yt(x, t, k, *parameters));
+//                (*gradient)[4 * k + 2] -= derror * ( eval_approx_da_yxx(x, t, k, *parameters) -  eval_approx_da_yt(x, t, k, *parameters));
+//                (*gradient)[4 * k + 3] -= derror * ( eval_approx_db_yxx(x, t, k, *parameters) -  eval_approx_db_yt(x, t, k, *parameters));
+//            }
+//            total_error += mem * mem / (train_size * train_size);
+//        }
+//    }
+//    return total_error;
+//}
+//
+//void solve_example_gradient(std::vector<double> &guess, double accuracy, size_t n_inner_neurons, size_t train_size, double ds, double de, size_t n_test_points, double ts, double te){
+//    std::cout << "Finding a solution via a Gradient Descent method with adaptive step-length" << std::endl;
+//    std::cout << "********************************************************************************************************************************************" <<std::endl;
+//    /* SETUP OF THE TRAINING DATA */
+//    std::vector<double> inp, out;
+//
+//    double frac, alpha, x;
+//    double grad_norm = accuracy * 10.0, mem, ai, bi, wxi, wti, error, derror, approx, t, gamma, total_error, sk, sy, sx, sg, beta;
+//    double grad_norm_prev = grad_norm;
+//    size_t i, j, k, iter_idx = 0;
+//
+//    /* TRAIN DATA FOR THE GOVERNING DE */
+//    std::vector<double> data_points(train_size);
+//
+//    /* ISOTROPIC TRAIN SET */
+//    frac = (de - ds) / (train_size - 1);
+//    for(i = 0; i < train_size; ++i){
+//        data_points[i] = frac * i + ds;
+////        std::cout << data_points[i] << std::endl;
+//    }
+//
+////    /* CHEBYSCHEV TRAIN SET */
+////    alpha = PI / (train_size );
+////    frac = 0.5 * (de - ds);
+////    for(i = 0; i < train_size; ++i){
+////        x = (std::cos(PI - alpha * i) + 1.0) * frac + ds;
+////        data_points[i] = x;
+////    }
+//
+//    std::vector<double> *gradient_current = new std::vector<double>(4 * n_inner_neurons);
+//    std::vector<double> *gradient_prev = new std::vector<double>(4 * n_inner_neurons);
+//    std::vector<double> *params_current = new std::vector<double>(guess);
+//    std::vector<double> *params_prev = new std::vector<double>(guess);
+//
+//    std::vector<double> *ptr_mem;
+//
+//    std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
+//    std::fill(gradient_prev->begin(), gradient_prev->end(), 0.0);
+//
+////    for (i = 0; i < n_inner_neurons; ++i) {
+////        wxi = (*params_current)[4 * i + 0];
+////        wti = (*params_current)[4 * i + 1];
+////        ai  = (*params_current)[4 * i + 2];
+////        bi  = (*params_current)[4 * i + 3];
+////
+////        printf("Path %3d. wx = %15.8f, wt = %15.8f, b = %15.8f, a = %15.8f\n", (int)(i + 1), wxi, wti, bi, ai);
+////    }
+//
+//    gamma = 1.0;
+//    double val = 0.0, prev_val, c = 2.0;
+//    prev_val = 0.0;
+//    while( grad_norm > accuracy) {
+//        iter_idx++;
+//        prev_val = val;
+//        grad_norm_prev = grad_norm;
+//
+//        /* reset of the current gradient */
+//        std::fill(gradient_current->begin(), gradient_current->end(), 0.0);
+//        val = calculate_gradient( data_points, n_inner_neurons, params_current, gradient_current );
+//
+//        grad_norm = 0.0;
+//        for(auto v: *gradient_current){
+//            grad_norm += v * v;
+//        }
+//        grad_norm = std::sqrt(grad_norm);
+//
+//        /* Update of the parameters */
+//        /* step length calculation */
+//        if(iter_idx < 10 || iter_idx % 100 == 0){
+//            /* fixed step length */
+//            gamma = 0.1 * accuracy;
+//        }
+//        else{
+//
+//
+//
+//            /* norm of the gradient calculation */
+//
+//            sk = 0.0;
+//            for(i = 0; i < gradient_current->size(); ++i){
+//                sx = (*gradient_current)[i] - (*gradient_prev)[i];
+//                sk += sx * sx;
+//            }
+//            sk = std::sqrt(sk);
+//
+//            /* angle between two consecutive gradients */
+//            double beta = 0.0, sx = 0.0;
+//            for(i = 0; i < gradient_current->size(); ++i){
+//                sx += (gradient_current->at( i ) * gradient_prev->at( i ));
+//            }
+//            sx /= grad_norm * grad_norm_prev;
+//            beta = std::sqrt(std::acos( sx ) / PI);
+//
+//
+////            eval_step_size_simple( gamma, val, prev_val, sk, grad_norm, grad_norm_prev );
+//            eval_step_size_mk( gamma, beta, c, grad_norm_prev, grad_norm, val, prev_val );
+//        }
+//
+//
+//
+//
+//
+//
+//        for(i = 0; i < gradient_current->size(); ++i){
+//            (*params_prev)[i] = (*params_current)[i] - gamma * (*gradient_current)[i];
+//        }
+//
+//        /* switcheroo */
+//        ptr_mem = gradient_prev;
+//        gradient_prev = gradient_current;
+//        gradient_current = ptr_mem;
+//
+//        ptr_mem = params_prev;
+//        params_prev = params_current;
+//        params_current = ptr_mem;
+//
+//
+//        if(iter_idx % 1 == 0){
+//            printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r", (int)iter_idx, gamma, c, grad_norm, val);
+//            std::cout.flush();
+//        }
+//    }
+//    printf("Iteration %12d. Step size: %15.8f, C: %15.8f, Gradient norm: %15.8f. Total error: %10.8f\r\n", (int)iter_idx, gamma, c, grad_norm, val);
+//    std::cout.flush();
+//
 //    for (i = 0; i < n_inner_neurons; ++i) {
 //        wxi = (*params_current)[4 * i + 0];
 //        wti = (*params_current)[4 * i + 1];
 //        ai  = (*params_current)[4 * i + 2];
 //        bi  = (*params_current)[4 * i + 3];
 //
-//        printf("%f/(1+e^(%f-%fx-%ft)) + ", (int)(i + 1), ai, bi, wxi, wti);
+//        printf("Path %3d. wx%d = %15.8f, wt%d = %15.8f, b%d = %15.8f, a%d = %15.8f\n", (int)(i + 1), (int)(i + 1), wxi , (int)(i + 1), wti, (int)(i + 1), bi, (int)(i + 1), ai);
 //    }
-//    printf("\n");
-
-
-    /* SOLUTION EXPORT */
-    printf("Exporting file 'data_2d_pde1_y.txt' : %7.3f%%\r", 0.0);
-    std::cout.flush();
-
-    std::vector<double> input, output(1);
-    std::ofstream ofs_y("data_2d_pde1_y.txt", std::ofstream::out);
-    frac = (te - ts) / (n_test_points - 1);
-    for(i = 0; i < n_test_points; ++i){
-        x = i * frac + ts;
-        for(j = 0; j < n_test_points; ++j){
-            t = j * frac + ts;
-            ofs_y << x << " " << t << " " << eval_approx_y(x, t, n_inner_neurons, *params_current) << std::endl;
-            printf("Exporting file 'data_2d_pde1_y.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
-            std::cout.flush();
-        }
-    }
-    printf("Exporting file 'data_2d_pde1_y.txt' : %7.3f%%\n", 100.0);
-    std::cout.flush();
-    ofs_y.close();
-
-    printf("Exporting file 'data_2d_pde1_yt.txt' : %7.3f%%\r", 0.0);
-    std::cout.flush();
-    std::ofstream ofs_t("data_2d_pde1_yt.txt", std::ofstream::out);
-    frac = (te - ts) / (n_test_points - 1);
-    for(i = 0; i < n_test_points; ++i){
-        x = i * frac + ts;
-        for(j = 0; j < n_test_points; ++j){
-            t = j * frac + ts;
-            ofs_t << x << " " << t << " " << eval_approx_yt(x, t, n_inner_neurons, *params_current) << std::endl;
-            printf("Exporting file 'data_2d_pde1_yt.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
-            std::cout.flush();
-        }
-    }
-    printf("Exporting file 'data_2d_pde1_yt.txt' : %7.3f%%\n", 100.0);
-    std::cout.flush();
-    ofs_t.close();
-
-    printf("Exporting file 'data_2d_pde1_yx.txt' : %7.3f%%\r", 0.0);
-    std::cout.flush();
-    std::ofstream ofs_x("data_2d_pde1_yx.txt", std::ofstream::out);
-    frac = (te - ts) / (n_test_points - 1);
-    for(i = 0; i < n_test_points; ++i){
-        x = i * frac + ts;
-        for(j = 0; j < n_test_points; ++j){
-            t = j * frac + ts;
-            ofs_x << x << " " << t << " " << eval_approx_yx(x, t, n_inner_neurons, *params_current) << std::endl;
-            printf("Exporting file 'data_2d_pde1_yx.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
-            std::cout.flush();
-        }
-    }
-    printf("Exporting file 'data_2d_pde1_yx.txt' : %7.3f%%\n", 100.0);
-    std::cout.flush();
-    ofs_x.close();
-
-    printf("Exporting file 'data_2d_pde1_yxx.txt' : %7.3f%%\r", 0.0);
-    std::cout.flush();
-    std::ofstream ofs_xx("data_2d_pde1_yxx.txt", std::ofstream::out);
-    frac = (te - ts) / (n_test_points - 1);
-    for(i = 0; i < n_test_points; ++i){
-        x = i * frac + ts;
-        for(j = 0; j < n_test_points; ++j){
-            t = j * frac + ts;
-            ofs_xx << x << " " << t << " " << eval_approx_yxx(x, t, n_inner_neurons, *params_current) << std::endl;
-            printf("Exporting file 'data_2d_pde1_yxx.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
-            std::cout.flush();
-        }
-    }
-    printf("Exporting file 'data_2d_pde1_yxx.txt' : %7.3f%%\n", 100.0);
-    std::cout.flush();
-    ofs_xx.close();
-
-    /* governing equation error */
-    std::ofstream ofs_error("data_2d_pde1_first_equation_error.txt", std::ofstream::out);
-    printf("Exporting file 'data_2d_pde1_first_equation_error.txt' : %7.3f%%\r", 0.0);
-    for(i = 0; i < n_test_points; ++i){
-        x = i * frac + ts;
-        for(j = 0; j < n_test_points; ++j){
-            t = j * frac + ts;
-            ofs_error << x << " " << t << " " << std::fabs(eval_approx_yxx(x, t, n_inner_neurons, *params_current) - eval_approx_yt(x, t, n_inner_neurons, *params_current)) << std::endl;
-            printf("Exporting file 'data_2d_pde1_first_equation_error.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
-            std::cout.flush();
-        }
-    }
-    printf("Exporting file 'data_2d_pde1_first_equation_error.txt' : %7.3f%%\n", 100.0);
-    std::cout.flush();
-    ofs_error.close();
-
-    /* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
-    /* first boundary condition & its error */
-    std::ofstream ofs_bc_t("data_1d_pde1_yt.txt", std::ofstream::out);
-    std::ofstream ofs_bc_x("data_1d_pde1_yx.txt", std::ofstream::out);
-    printf("Exporting files 'data_1d_pde1_yt.txt' and 'data_1d_pde1_yx.txt' : %7.3f%%\r", 0.0);
-    for(i = 0; i < n_test_points; ++i){
-        x = frac * i + ts;
-        t = frac * i + ts;
-
-        double yt = std::sin(t);
-        double yx = std::pow(E, -0.707106781 * x) * std::sin( -0.707106781 * x );
-
-        double evalt = eval_approx_y(0, t, n_inner_neurons, *params_current);
-        double evalx = eval_approx_y(x, 0, n_inner_neurons, *params_current);
-
-        ofs_bc_t << i + 1 << " " << t << " " << yt << " " << evalt << " " << std::fabs(evalt - yt) << std::endl;
-        ofs_bc_x << i + 1 << " " << x << " " << yx << " " << evalx << " " << std::fabs(evalx - yx) << std::endl;
-
-        printf("Exporting files 'data_1d_pde1_yt.txt' and 'data_1d_pde1_yx.txt' : %7.3f%%\r", (100.0 * i) / (n_test_points - 1));
-        std::cout.flush();
-    }
-    printf("Exporting files 'data_1d_pde1_yt.txt' and 'data_1d_pde1_yx.txt' : %7.3f%%\r", 100.0);
-    std::cout.flush();
-    ofs_bc_t.close();
-    ofs_bc_x.close();
-
-    delete gradient_current;
-    delete gradient_prev;
-    delete params_current;
-    delete params_prev;
-}
+//    std::cout << "********************************************************************************************************************************************" <<std::endl;
+//
+////    for (i = 0; i < n_inner_neurons; ++i) {
+////        wxi = (*params_current)[4 * i + 0];
+////        wti = (*params_current)[4 * i + 1];
+////        ai  = (*params_current)[4 * i + 2];
+////        bi  = (*params_current)[4 * i + 3];
+////
+////        printf("%f/(1+e^(%f-%fx-%ft)) + ", (int)(i + 1), ai, bi, wxi, wti);
+////    }
+////    printf("\n");
+//
+//
+//    /* SOLUTION EXPORT */
+//    printf("Exporting file 'data_2d_pde1_y.txt' : %7.3f%%\r", 0.0);
+//    std::cout.flush();
+//
+//    std::vector<double> input, output(1);
+//    std::ofstream ofs_y("data_2d_pde1_y.txt", std::ofstream::out);
+//    frac = (te - ts) / (n_test_points - 1);
+//    for(i = 0; i < n_test_points; ++i){
+//        x = i * frac + ts;
+//        for(j = 0; j < n_test_points; ++j){
+//            t = j * frac + ts;
+//            ofs_y << x << " " << t << " " << eval_approx_y(x, t, n_inner_neurons, *params_current) << std::endl;
+//            printf("Exporting file 'data_2d_pde1_y.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
+//            std::cout.flush();
+//        }
+//    }
+//    printf("Exporting file 'data_2d_pde1_y.txt' : %7.3f%%\n", 100.0);
+//    std::cout.flush();
+//    ofs_y.close();
+//
+//    printf("Exporting file 'data_2d_pde1_yt.txt' : %7.3f%%\r", 0.0);
+//    std::cout.flush();
+//    std::ofstream ofs_t("data_2d_pde1_yt.txt", std::ofstream::out);
+//    frac = (te - ts) / (n_test_points - 1);
+//    for(i = 0; i < n_test_points; ++i){
+//        x = i * frac + ts;
+//        for(j = 0; j < n_test_points; ++j){
+//            t = j * frac + ts;
+//            ofs_t << x << " " << t << " " << eval_approx_yt(x, t, n_inner_neurons, *params_current) << std::endl;
+//            printf("Exporting file 'data_2d_pde1_yt.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
+//            std::cout.flush();
+//        }
+//    }
+//    printf("Exporting file 'data_2d_pde1_yt.txt' : %7.3f%%\n", 100.0);
+//    std::cout.flush();
+//    ofs_t.close();
+//
+//    printf("Exporting file 'data_2d_pde1_yx.txt' : %7.3f%%\r", 0.0);
+//    std::cout.flush();
+//    std::ofstream ofs_x("data_2d_pde1_yx.txt", std::ofstream::out);
+//    frac = (te - ts) / (n_test_points - 1);
+//    for(i = 0; i < n_test_points; ++i){
+//        x = i * frac + ts;
+//        for(j = 0; j < n_test_points; ++j){
+//            t = j * frac + ts;
+//            ofs_x << x << " " << t << " " << eval_approx_yx(x, t, n_inner_neurons, *params_current) << std::endl;
+//            printf("Exporting file 'data_2d_pde1_yx.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
+//            std::cout.flush();
+//        }
+//    }
+//    printf("Exporting file 'data_2d_pde1_yx.txt' : %7.3f%%\n", 100.0);
+//    std::cout.flush();
+//    ofs_x.close();
+//
+//    printf("Exporting file 'data_2d_pde1_yxx.txt' : %7.3f%%\r", 0.0);
+//    std::cout.flush();
+//    std::ofstream ofs_xx("data_2d_pde1_yxx.txt", std::ofstream::out);
+//    frac = (te - ts) / (n_test_points - 1);
+//    for(i = 0; i < n_test_points; ++i){
+//        x = i * frac + ts;
+//        for(j = 0; j < n_test_points; ++j){
+//            t = j * frac + ts;
+//            ofs_xx << x << " " << t << " " << eval_approx_yxx(x, t, n_inner_neurons, *params_current) << std::endl;
+//            printf("Exporting file 'data_2d_pde1_yxx.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
+//            std::cout.flush();
+//        }
+//    }
+//    printf("Exporting file 'data_2d_pde1_yxx.txt' : %7.3f%%\n", 100.0);
+//    std::cout.flush();
+//    ofs_xx.close();
+//
+//    /* governing equation error */
+//    std::ofstream ofs_error("data_2d_pde1_first_equation_error.txt", std::ofstream::out);
+//    printf("Exporting file 'data_2d_pde1_first_equation_error.txt' : %7.3f%%\r", 0.0);
+//    for(i = 0; i < n_test_points; ++i){
+//        x = i * frac + ts;
+//        for(j = 0; j < n_test_points; ++j){
+//            t = j * frac + ts;
+//            ofs_error << x << " " << t << " " << std::fabs(eval_approx_yxx(x, t, n_inner_neurons, *params_current) - eval_approx_yt(x, t, n_inner_neurons, *params_current)) << std::endl;
+//            printf("Exporting file 'data_2d_pde1_first_equation_error.txt' : %7.3f%%\r", (100.0 * (j + i * n_test_points)) / (n_test_points * n_test_points - 1));
+//            std::cout.flush();
+//        }
+//    }
+//    printf("Exporting file 'data_2d_pde1_first_equation_error.txt' : %7.3f%%\n", 100.0);
+//    std::cout.flush();
+//    ofs_error.close();
+//
+//    /* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
+//    /* first boundary condition & its error */
+//    std::ofstream ofs_bc_t("data_1d_pde1_yt.txt", std::ofstream::out);
+//    std::ofstream ofs_bc_x("data_1d_pde1_yx.txt", std::ofstream::out);
+//    printf("Exporting files 'data_1d_pde1_yt.txt' and 'data_1d_pde1_yx.txt' : %7.3f%%\r", 0.0);
+//    for(i = 0; i < n_test_points; ++i){
+//        x = frac * i + ts;
+//        t = frac * i + ts;
+//
+//        double yt = std::sin(t);
+//        double yx = std::pow(E, -0.707106781 * x) * std::sin( -0.707106781 * x );
+//
+//        double evalt = eval_approx_y(0, t, n_inner_neurons, *params_current);
+//        double evalx = eval_approx_y(x, 0, n_inner_neurons, *params_current);
+//
+//        ofs_bc_t << i + 1 << " " << t << " " << yt << " " << evalt << " " << std::fabs(evalt - yt) << std::endl;
+//        ofs_bc_x << i + 1 << " " << x << " " << yx << " " << evalx << " " << std::fabs(evalx - yx) << std::endl;
+//
+//        printf("Exporting files 'data_1d_pde1_yt.txt' and 'data_1d_pde1_yx.txt' : %7.3f%%\r", (100.0 * i) / (n_test_points - 1));
+//        std::cout.flush();
+//    }
+//    printf("Exporting files 'data_1d_pde1_yt.txt' and 'data_1d_pde1_yx.txt' : %7.3f%%\r", 100.0);
+//    std::cout.flush();
+//    ofs_bc_t.close();
+//    ofs_bc_x.close();
+//
+//    delete gradient_current;
+//    delete gradient_prev;
+//    delete params_current;
+//    delete params_prev;
+//}
 
 void solve_example_particle_swarm(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;
@@ -702,22 +702,43 @@ void solve_example_particle_swarm(double accuracy, size_t n_inner_neurons, size_
     solver_01.set_error_function( 1, ErrorFunctionType::ErrorFuncMSE, &ds_t );
     solver_01.set_error_function( 2, ErrorFunctionType::ErrorFuncMSE, &ds_x );
 
-    size_t total_dim = (2 + n_inputs) * n_inner_neurons;
-
-    /* PARTICLE SWARM TRAINING METHOD SETUP */
-
+    size_t total_dim = (solver_01.get_solution( alpha_00 )->get_n_biases() + solver_01.get_solution( alpha_00 )->get_n_weights());
 
-    //must encapsulate each of the partial error functions
+    /* TRAINING METHOD SETUP */
     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;
+
+    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, c2 = 1.7, w = 0.7;
-    double gamma = 0.5, epsilon = 0.02, delta = 0.9;
+    double c1 = 1.7;
+    double c2 = 1.7;
+    double 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.7;
+
+    ParticleSwarm swarm_01(
+            &domain_bounds,
+            c1,
+            c2,
+            w,
+            gamma,
+            epsilon,
+            delta,
+            n_particles,
+            max_iters
+    );
+    solver_01.solve( swarm_01 );
+
 
-    solver_01.solve_via_particle_swarm( &domain_bounds, c1, c2, w, n_particles, max_iters, gamma, epsilon, delta );
     size_t i, j;
     std::vector<double> *w1_ptr = solver_01.get_solution( alpha_00 )->get_parameter_ptr_weights();
     std::vector<double> *w2_ptr = solver_01.get_solution( alpha_00 )->get_parameter_ptr_biases();
diff --git a/src/examples/network_serialization.cpp b/src/examples/network_serialization.cpp
index 71944de610cfee28938526a019c4f9b5c45f731f..caf4b1b0102cce22b1080eaeb79c99e2f03d79c0 100644
--- a/src/examples/network_serialization.cpp
+++ b/src/examples/network_serialization.cpp
@@ -78,8 +78,19 @@ int main() {
     MSE mse(&net, &ds);
 
     /* TRAINING METHOD SETUP */
-    std::vector<double> domain_bounds = {-10.0, 10.0, -10.0, 10.0, -10.0, 10.0};
-    ParticleSwarm swarm_01(&mse, &domain_bounds);
+    std::vector<double> domain_bounds(2 * (net.get_n_weights() + net.get_n_biases()));
+
+    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.7;
+    size_t n_particles = 50;
+    size_t iter_max = 1000;
+
     /* 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;
@@ -87,10 +98,22 @@ int main() {
     /* 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;
-    swarm_01.optimize(gamma, epsilon, delta);
-
-    std::vector<double> *parameters = swarm_01.get_solution();
+    double delta = 0.7;
+
+    ParticleSwarm swarm_01(
+            &domain_bounds,
+            c1,
+            c2,
+            w,
+            gamma,
+            epsilon,
+            delta,
+            n_particles,
+            iter_max
+    );
+    swarm_01.optimize( mse );
+
+    std::vector<double> *parameters = swarm_01.get_parameters();
     net.copy_parameter_space(parameters);
 
     printf("w1 = %10.7f\n", parameters->at( 0 ));
diff --git a/src/examples/seminar.cpp b/src/examples/seminar.cpp
index 800e27df1bf4f6c70ea8f3ca92237dbfb45f55ec..6369f099c4b0f688e2c56abf022f30495228a5be 100644
--- a/src/examples/seminar.cpp
+++ b/src/examples/seminar.cpp
@@ -76,14 +76,18 @@ int main() {
 
 
     /* TRAINING METHOD SETUP */
-    std::vector<double> domain_bounds(2 * (2 + 6));
+    std::vector<double> domain_bounds(2 * (XOR.get_n_weights() + XOR.get_n_biases()));
 
     for(size_t i = 0; i < domain_bounds.size() / 2; ++i){
         domain_bounds[2 * i] = -10;
         domain_bounds[2 * i + 1] = 10;
     }
 
-    ParticleSwarm swarm_01(&mse, &domain_bounds);
+    double c1 = 1.7;
+    double c2 = 1.7;
+    double w = 0.7;
+    size_t n_particles = 50;
+    size_t iter_max = 1000;
 
     /* if the maximal velocity from the previous step is less than 'gamma' times the current maximal velocity, then one
      * terminating criterion is met */
@@ -93,9 +97,21 @@ int main() {
      * terminating criterion is met ('n' is the total number of particles) */
     double epsilon = 0.02;
     double delta = 0.7;
-    swarm_01.optimize(gamma, epsilon, delta);
 
-    std::vector<double> *parameters = swarm_01.get_solution();
+    ParticleSwarm swarm_01(
+            &domain_bounds,
+            c1,
+            c2,
+            w,
+            gamma,
+            epsilon,
+            delta,
+            n_particles,
+            iter_max
+    );
+    swarm_01.optimize( mse );
+
+    std::vector<double> *parameters = swarm_01.get_parameters( );
     XOR.copy_parameter_space(parameters);
 
     /* ERROR CALCULATION */
diff --git a/src/tests/NeuralNetwork_test.cpp b/src/tests/NeuralNetwork_test.cpp
index 4980fb3d46c0bef7e5779bcb1e66918c292f7724..f0c6df790656d545c91b41556dbd806a03d53a77 100644
--- a/src/tests/NeuralNetwork_test.cpp
+++ b/src/tests/NeuralNetwork_test.cpp
@@ -59,7 +59,6 @@ BOOST_AUTO_TEST_SUITE(NeuralNetwork_test)
         Neuron *n1 = new NeuronLinear();
         Neuron *n2 = new NeuronLinear();
         Neuron *n3 = new NeuronLinear();
-        Neuron *n4 = new NeuronLinear();
 
         NeuralNetwork network;
 
@@ -68,7 +67,7 @@ BOOST_AUTO_TEST_SUITE(NeuralNetwork_test)
 
         BOOST_CHECK_EQUAL(1, network.add_neuron(n2, BIAS_TYPE::NEXT_BIAS));
 
-        BOOST_CHECK_EQUAL(2, network.add_neuron(n4, BIAS_TYPE::NO_BIAS));
+        BOOST_CHECK_EQUAL(2, network.add_neuron(n3, BIAS_TYPE::NO_BIAS));
 
         BOOST_CHECK_EQUAL(3, network.get_n_neurons());
 
diff --git a/src/tests/ParticleSwarm_test.cpp b/src/tests/ParticleSwarm_test.cpp
index f61294f4b175a224f98c49c79844cea87266a202..e61e235f4928b0ddfdf89590c901a79e51d3062d 100644
--- a/src/tests/ParticleSwarm_test.cpp
+++ b/src/tests/ParticleSwarm_test.cpp
@@ -40,19 +40,20 @@ BOOST_AUTO_TEST_SUITE(ParticleSwarm_test)
         std::vector<double> domain_bound;
         domain_bound.push_back(5);
         NeuralNetwork network;
-        std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec;
-        std::vector<double> inp, out;
-
-        for (int i = 0; i < 3; i++) {
-            inp.push_back(i);
-            out.push_back(i + 4);
-        }
-
-        data_vec.emplace_back(std::make_pair(inp, out));
-
-        DataSet dataSet(&data_vec);
-        ErrorFunction *error = new MSE(&network, &dataSet);
-        BOOST_CHECK_NO_THROW(ParticleSwarm swarm(error, &domain_bound, 0, 1, 1, 0, 20));
+//        std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec;
+//        std::vector<double> inp, out;
+//
+//        for (int i = 0; i < 3; i++) {
+//            inp.push_back(i);
+//            out.push_back(i + 4);
+//        }
+//
+//        data_vec.emplace_back(std::make_pair(inp, out));
+//
+//        DataSet dataSet(&data_vec);
+//        ErrorFunction *error = new MSE(&network, &dataSet);
+
+        BOOST_CHECK_NO_THROW(ParticleSwarm swarm(&domain_bound, 0, 1, 1, 0.5, 0.05, 0.5, 0, 20));
     }
 
     BOOST_AUTO_TEST_CASE(ParticleSwarm_optimalize_test){
@@ -71,10 +72,15 @@ BOOST_AUTO_TEST_SUITE(ParticleSwarm_test)
 
         DataSet dataSet(&data_vec);
         ErrorFunction *error = new MSE(&network, &dataSet);
-        ParticleSwarm swarm(error,&domain_bound, 0, 1, 1, 0, 20);
-        BOOST_CHECK_THROW(swarm.optimize(-1,1,1), std::invalid_argument) ;
-        BOOST_CHECK_THROW(swarm.optimize(1,-1,1), std::invalid_argument) ;
-        BOOST_CHECK_THROW(swarm.optimize(1,1,-1), std::invalid_argument) ;
+
+        ParticleSwarm swarm(&domain_bound, 0, 1, 1, -1,1,1, 0, 20);
+        BOOST_CHECK_THROW(swarm.optimize( *error ), std::invalid_argument) ;
+
+        ParticleSwarm swarm2(&domain_bound, 0, 1, 1, 1,-1,1, 0, 20);
+        BOOST_CHECK_THROW(swarm2.optimize( *error ), std::invalid_argument) ;
+
+        ParticleSwarm swarm3(&domain_bound, 0, 1, 1, 1,1,-1, 0, 20);
+        BOOST_CHECK_THROW(swarm3.optimize( *error ), std::invalid_argument) ;
 
     }
 
diff --git a/src/tests/Particle_test.cpp b/src/tests/Particle_test.cpp
index 590b7b1b1586d9f4307b18b3f18598bc5e31a9c4..3d1dc169d83b7962a9756bbe1ee0d76655d423d6 100644
--- a/src/tests/Particle_test.cpp
+++ b/src/tests/Particle_test.cpp
@@ -31,7 +31,7 @@
  BOOST_AUTO_TEST_SUITE(Particle_test)
 
     BOOST_AUTO_TEST_CASE(Particle_construction_test) {
-        double domain_bound[5] = {1, 2, 3, 4, 5};
+        std::vector<double> domain_bound = {1, 2, 3, 4, 5};
         Neuron *n1 = new NeuronLinear();
         Neuron *n2 = new NeuronLinear();
         NeuralNetwork network;
@@ -60,13 +60,13 @@
 
         DataSet dataSet(&data_vec);
         ErrorFunction *error = new MSE(&network, &dataSet);
-        Particle particle(error, &domain_bound[0]);
-        BOOST_CHECK_NO_THROW(Particle particle(error, &domain_bound[0]));
+        Particle particle(error, &domain_bound);
+        BOOST_CHECK_NO_THROW(Particle particle(error, &domain_bound));
         //  particle.get_coordinate();
     }
 
     BOOST_AUTO_TEST_CASE(Particle_get_coordinate_test) {
-        double domain_bound[5] = {1, 2, 3, 4, 5};
+         std::vector<double> domain_bound = {1, 2, 3, 4, 5};
         Neuron *n1 = new NeuronLinear();
         Neuron *n2 = new NeuronLinear();
         NeuralNetwork network;
@@ -95,8 +95,8 @@
 
         DataSet dataSet(&data_vec);
         ErrorFunction *error = new MSE(&network, &dataSet);
-        Particle particle1(error, &domain_bound[0]);
-        Particle particle2(error, &domain_bound[0]);
+        Particle particle1( error, &domain_bound );
+        Particle particle2( error, &domain_bound );
 
         BOOST_CHECK(*particle1.get_coordinate() != *particle2.get_coordinate());
     }