diff --git a/include/4neuro.h b/include/4neuro.h
index 21e4e3935b8d16c61d6045e353ea7eea661b1da0..5776cdd38be57fa333ddc27e5cb5bdfbd5c87e2c 100644
--- a/include/4neuro.h
+++ b/include/4neuro.h
@@ -14,6 +14,7 @@
 #include "../src/NetConnection/ConnectionWeight.h"
 #include "../src/NetConnection/ConnectionWeightIdentity.h"
 #include "../src/Network/NeuralNetwork.h"
+#include "../src/Network/NeuralNetworkSum.h"
 #include "../src/Neuron/Neuron.h"
 #include "../src/Neuron/NeuronBinary.h"
 #include "../src/Neuron/NeuronLinear.h"
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ebdb120b32889a9c529555f821d23e763cbe91fb..ac3b7371a8962969800e44e6849722ae820dc202 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -13,7 +13,7 @@ add_library(4neuro SHARED
         NetConnection/ConnectionWeightIdentity.cpp
         LearningMethods/ParticleSwarm.cpp
         DataSet/DataSet.cpp
-        ErrorFunction/ErrorFunctions.cpp)
+        ErrorFunction/ErrorFunctions.cpp Network/NeuralNetworkSum.cpp Network/NeuralNetworkSum.h)
 
 target_link_libraries(4neuro boost_serialization)
 
@@ -37,6 +37,9 @@ target_link_libraries(net_test_2 4neuro)
 add_executable(net_test_3 net_test_3.cpp)
 target_link_libraries(net_test_3 4neuro)
 
+add_executable(net_test_ode_1 net_test_ode_1.cpp)
+target_link_libraries(net_test_ode_1 4neuro)
+
 ##############
 # UNIT TESTS #
 ##############
diff --git a/src/DataSet/DataSet.h b/src/DataSet/DataSet.h
index 6486f20286713c6f52a9cdbb04d0e8e81680a10c..8fd1d585d8ffb535d3c75c0c6dcae3bfb05416d9 100644
--- a/src/DataSet/DataSet.h
+++ b/src/DataSet/DataSet.h
@@ -1,6 +1,7 @@
 //
 // Created by martin on 7/13/18.
 //
+//TODO generovani dat
 
 #ifndef INC_4NEURO_DATASET_H
 #define INC_4NEURO_DATASET_H
diff --git a/src/ErrorFunction/ErrorFunctions.cpp b/src/ErrorFunction/ErrorFunctions.cpp
index 47187ffa41c2adc9de7d9b8391cbcebc089fcafb..a98e0603295b882c81b5f4b2422ca179eec66a0e 100644
--- a/src/ErrorFunction/ErrorFunctions.cpp
+++ b/src/ErrorFunction/ErrorFunctions.cpp
@@ -53,6 +53,7 @@ double MSE::eval(double *weights) {
 
 MSE_SUM::MSE_SUM() {
     this->summand = nullptr;
+    this->summand_coefficient = nullptr;
     this->dimension = 0;
 }
 
@@ -60,23 +61,31 @@ MSE_SUM::~MSE_SUM(){
     if( this->summand ){
         delete this->summand;
     }
+    if( this->summand_coefficient ){
+        delete this->summand_coefficient;
+    }
 }
 
 double MSE_SUM::eval(double *weights) {
     double output = 0.0;
 
-    for(ErrorFunction *f: *this->summand){
-        output += f->eval( weights );
+    for( int i = 0; i < this->summand->size(); ++i ){
+        output += this->summand->at( i )->eval( weights ) * this->summand_coefficient->at( i );
     }
 
     return output;
 }
 
-void MSE_SUM::add_error_function(ErrorFunction *F) {
+void MSE_SUM::add_error_function(ErrorFunction *F, double alpha) {
     if(!this->summand){
         this->summand = new std::vector<ErrorFunction*>(0);
     }
-    this->summand->push_back(F);
+    this->summand->push_back( F );
+
+    if(!this->summand_coefficient){
+        this->summand_coefficient = new std::vector<double>(0);
+    }
+    this->summand_coefficient->push_back( alpha );
 
     if(F->get_dimension() > this->dimension){
         this->dimension = F->get_dimension();
diff --git a/src/ErrorFunction/ErrorFunctions.h b/src/ErrorFunction/ErrorFunctions.h
index bb639d587db6db95814cae7a6c514a457fb7aa66..48c0731e64c532111620bb7049ed3d9927a0a5ab 100644
--- a/src/ErrorFunction/ErrorFunctions.h
+++ b/src/ErrorFunction/ErrorFunctions.h
@@ -78,7 +78,7 @@ public:
      *
      * @param F
      */
-    void add_error_function(ErrorFunction *F);
+    void add_error_function(ErrorFunction *F, double alpha = 1.0);
 
     /**
      *
@@ -88,6 +88,7 @@ public:
 
 private:
     std::vector<ErrorFunction*>* summand;
+    std::vector<double> *summand_coefficient;
 };
 
 
diff --git a/src/LearningMethods/ParticleSwarm.cpp b/src/LearningMethods/ParticleSwarm.cpp
index 42004f21069a5e1e198debfd1d942068910c9631..77bdff10562285c8be25c2c2d147bff79c7f4c46 100644
--- a/src/LearningMethods/ParticleSwarm.cpp
+++ b/src/LearningMethods/ParticleSwarm.cpp
@@ -46,6 +46,7 @@ Particle::Particle(ErrorFunction* ef, double *domain_bounds) {
         this->optimal_coordinate[i] = this->coordinate[i];
     }
 
+//    printf("coordinate_dim: %d\n", this->coordinate_dim);
     this->optimal_value = this->ef->eval(this->coordinate);
 
 //    this->print_coordinate();
diff --git a/src/NetConnection/Connection.cpp b/src/NetConnection/Connection.cpp
index 6087bfba74a3fe0ad05b7df7165c18d7644a05d0..cd6b07918e0900f18d4658e96ab34c82b572746b 100644
--- a/src/NetConnection/Connection.cpp
+++ b/src/NetConnection/Connection.cpp
@@ -45,7 +45,7 @@ Neuron* Connection::get_neuron_out() {
 //}
 
 void Connection::pass_signal() {
-//    printf("passing signal between neurons %d -> %d, value: %f * %f\n", this->neuron_in, this->neuron_out, this->neuron_in->get_state(), this->con->eval());
+//    printf("  passing signal between neurons %d -> %d, value: %f * %f\n", this->neuron_in, this->neuron_out, this->neuron_in->get_state(), this->con->eval());
     this->neuron_out->adjust_potential(this->neuron_in->get_state() * this->con->eval());
 }
 
diff --git a/src/Network/NeuralNetwork.cpp b/src/Network/NeuralNetwork.cpp
index c3ccd6a109dce494a4ee3385560a9768d56c8655..0cde711c3af3c0e4c862063062d4c30eb71336dd 100644
--- a/src/Network/NeuralNetwork.cpp
+++ b/src/Network/NeuralNetwork.cpp
@@ -274,21 +274,22 @@ int NeuralNetwork::add_neuron(Neuron *n) {
     return n->get_idx();
 }
 
-void NeuralNetwork::add_connection_simple(int n1_idx, int n2_idx) {
-    add_connection_simple(n1_idx, n2_idx, -1);
+size_t NeuralNetwork::add_connection_simple(int n1_idx, int n2_idx) {
+    return add_connection_simple(n1_idx, n2_idx, -1);
 }
 
-void NeuralNetwork::add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx) {
-    add_connection_simple(n1_idx, n2_idx, weight_idx, 1);
+size_t NeuralNetwork::add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx) {
+    return add_connection_simple(n1_idx, n2_idx, weight_idx, 1);
 }
 
-void NeuralNetwork::add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx, double weight_value) {
+size_t NeuralNetwork::add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx, double weight_value) {
     // TODO generate weight_value automatically from normal distribution
 
     if(weight_idx < 0 || weight_idx >= this->connection_weights->size()){
         //this weight is a new one, we add it to the system of weights
         this->connection_weights->push_back(weight_value);
         weight_idx = (int)this->connection_weights->size() - 1;
+        this->n_weights++;
     }
     Neuron *neuron_out = this->neurons->at(n1_idx);
     Neuron *neuron_in = this->neurons->at(n2_idx);
@@ -300,10 +301,42 @@ void NeuralNetwork::add_connection_simple(int n1_idx, int n2_idx, size_t weight_
 
     neuron_out->add_connection_out(u1u2);
     neuron_in->add_connection_in(u1u2);
+
+    return weight_idx;
+}
+
+void NeuralNetwork::add_connection_shared_power1(int n1_idx, int n2_idx, size_t weight_idx) {
+    std::function<double(double *, size_t*, size_t)> weight_function = [](double * weight_array, size_t * index_array, size_t n_params){
+        return weight_array[index_array[0]];
+    };
+
+    size_t* weight_indices = new size_t[1];
+    double* weight_values = new double[1];
+
+    weight_indices[0] = weight_idx;
+    this->add_connection_general( n1_idx, n2_idx, &weight_function, weight_indices, weight_values, 1);
+
+    delete [] weight_indices;
+    delete [] weight_values;
+}
+
+void NeuralNetwork::add_connection_shared_power2(int n1_idx, int n2_idx, size_t weight_idx) {
+    std::function<double(double *, size_t*, size_t)> weight_function = [](double * weight_array, size_t * index_array, size_t n_params){
+        return weight_array[index_array[0]] * weight_array[index_array[0]];
+    };
+
+    size_t* weight_indices = new size_t[1];
+    double* weight_values = new double[1];
+
+    weight_indices[0] = weight_idx;
+    this->add_connection_general( n1_idx, n2_idx, &weight_function, weight_indices, weight_values, 1);
+
+    delete [] weight_indices;
+    delete [] weight_values;
 }
 
-void NeuralNetwork::add_connection_general(int n1_idx, int n2_idx, std::function<double(double *, int*, int)> *f,
-                                           int* weight_indices, double* weight_values, size_t n_weights) {
+void NeuralNetwork::add_connection_general(int n1_idx, int n2_idx, std::function<double(double *, size_t*, size_t)> *f,
+                                           size_t* weight_indices, double* weight_values, size_t n_weights) {
 
     ConnectionWeight *con_weight_u1u2 = new ConnectionWeight(n_weights, this->connection_weights);
     //we analyze weights
@@ -328,6 +361,7 @@ void NeuralNetwork::add_connection_general(int n1_idx, int n2_idx, std::function
 
     neuron_out->add_connection_out(u1u2);
     neuron_in->add_connection_in(u1u2);
+
 }
 
 void NeuralNetwork::determine_inputs_outputs() {
@@ -372,6 +406,8 @@ void NeuralNetwork::set_weight_array(std::vector<double> *weight_ptr) {
     }
     this->connection_weights = weight_ptr;
     this->delete_weights = false;
+
+    this->n_weights = this->connection_weights->size();
 }
 
 void NeuralNetwork::eval_single(std::vector<double> &input, std::vector<double> &output, double * custom_weights) {
@@ -448,8 +484,8 @@ void NeuralNetwork::eval_single(std::vector<double> &input, std::vector<double>
         //we iterate through the active neurons and propagate the signal
         for(i = 0; i < active_set_size[idx1]; ++i){
             active_neuron = this->active_eval_set->at(i + n * idx1);
-            active_neuron->activate();
-//            printf(" active neuron %d, state: %f\n", active_neuron, active_neuron->get_state());
+            active_neuron->activate();//computes the activation function on its potential
+//            printf(" active neuron %d, potential: %f, state: %f\n", active_neuron, active_neuron->get_potential(), active_neuron->get_state());
 
             for(Connection* connection: *(active_neuron->get_connections_out())){
                 connection->pass_signal();
@@ -463,6 +499,7 @@ void NeuralNetwork::eval_single(std::vector<double> &input, std::vector<double>
                 }
             }
         }
+//        printf("-----------\n");
 
         idx1 = idx2;
         idx2 = (idx1 + 1) % 2;
@@ -489,20 +526,25 @@ void NeuralNetwork::copy_weights(double *weights) {
     }
 }
 
+std::vector<double>* NeuralNetwork::get_weights_ptr(){
+    return this->connection_weights;
+}
+
 void NeuralNetwork::randomize_weights() {
 
-    if( this->neurons->size() <= 0 || !this->in_out_determined ){
+    if( this->n_weights <= 0 ){
+
         return;
     }
 
     boost::random::mt19937 gen;
 
     // Init weight guess ("optimal" for logistic activation functions)
-    double r = 4 * sqrt(6./(this->n_inputs + this->n_outputs));
+    double r = 4 * sqrt(6./(this->n_weights));
 
     boost::random::uniform_real_distribution<> dist(-r, r);
 
-    for(unsigned int i = 0; i < this->n_weights; i++) {
+    for(size_t i = 0; i < this->n_weights; i++) {
         this->connection_weights->at(i) = dist(gen);
     }
 }
@@ -552,3 +594,18 @@ void NeuralNetwork::specify_output_neurons(std::vector<size_t> &output_neurons_i
     this->n_outputs = this->output_neurons->size();
 }
 
+void NeuralNetwork::print_weights() {
+    printf("Connection weights: ");
+    if(this->connection_weights){
+        for( size_t i = 0; i < this->connection_weights->size() - 1; ++i){
+            printf("%f, ", this->connection_weights->at(i));
+        }
+        printf("%f", this->connection_weights->at(this->connection_weights->size() - 1));
+    }
+
+    printf("\n");
+}
+
+void NeuralNetwork::print_stats(){
+    printf("Number of neurons: %d, number of weights: %d\n", this->neurons->size(), this->n_weights);
+}
\ No newline at end of file
diff --git a/src/Network/NeuralNetwork.h b/src/Network/NeuralNetwork.h
index ceaa554e182057163534a7ab68eb4dc5953d1d3e..1bfa39a5c47e8bd1110958e1a0291ac85a42e1fe 100644
--- a/src/Network/NeuralNetwork.h
+++ b/src/Network/NeuralNetwork.h
@@ -79,11 +79,7 @@ private:
       */
      void determine_inputs_outputs();
 
-     /**
-      *
-      * @param weight_ptr
-      */
-     void set_weight_array( std::vector<double>* weight_ptr );
+
 
 public:
 
@@ -113,7 +109,7 @@ public:
      * @param[in] input
      * @param[in,out] output
      */
-    void eval_single(std::vector<double> &input, std::vector<double> &output, double *custom_weights = nullptr);
+    virtual void eval_single(std::vector<double> &input, std::vector<double> &output, double *custom_weights = nullptr);
 
 
 
@@ -128,25 +124,46 @@ public:
      *
      * @param n1_idx
      * @param n2_idx
+     * @return
+     */
+    size_t add_connection_simple(int n1_idx, int n2_idx);
+
+    /**
+     *
+     * @param n1_idx
+     * @param n2_idx
+     * @param weight_idx
+     * @return
+     */
+    size_t add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx);
+
+    /**
+     *
+     * @param n1_idx
+     * @param n2_idx
+     * @param weight_idx
+     * @param weight_value
+     * @return
      */
-    void add_connection_simple(int n1_idx, int n2_idx);
+    size_t add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx, double weight_value);
 
     /**
      *
      * @param n1_idx
      * @param n2_idx
      * @param weight_idx
+     * @param power_degree
      */
-    void add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx);
+    void add_connection_shared_power1(int n1_idx, int n2_idx, size_t weight_idx);
 
     /**
      *
-     * @param[in] n1_idx
-     * @param[in] n2_idx
-     * @param[in] weight_idx
-     * @param[in] weight_value
+     * @param n1_idx
+     * @param n2_idx
+     * @param weight_idx
+     * @param power_degree
      */
-    void add_connection_simple(int n1_idx, int n2_idx, size_t weight_idx, double weight_value);
+    void add_connection_shared_power2(int n1_idx, int n2_idx, size_t weight_idx);
 
     /**
      *
@@ -157,8 +174,14 @@ public:
      * @param weight_values
      * @param n_weights
      */
-    void add_connection_general(int n1_idx, int n2_idx, std::function<double(double *, int*, int)> *f,
-                                int* weight_indices, double* weight_values, size_t n_weights);
+    void add_connection_general(int n1_idx, int n2_idx, std::function<double(double *, size_t*, size_t)> *f,
+                                size_t* weight_indices, double* weight_values, size_t n_weights);
+
+    /**
+     *
+     * @param weight_ptr
+     */
+    void set_weight_array( std::vector<double>* weight_ptr );
 
     /**
      *
@@ -166,6 +189,12 @@ public:
      */
     void copy_weights(double *weights);
 
+    /**
+     *
+     * @return
+     */
+    std::vector<double>* get_weights_ptr();
+
     /**
      *
      */
@@ -187,7 +216,7 @@ public:
      *
      * @return
      */
-    size_t get_n_weights();
+    virtual size_t get_n_weights();
 
     /**
      *
@@ -201,6 +230,15 @@ public:
      */
     void specify_output_neurons(std::vector<size_t> &output_neurons_indices);
 
+    /**
+     *
+     */
+    void print_weights();
+
+    /**
+     *
+     */
+    void print_stats();
 };
 
 
diff --git a/src/Network/NeuralNetworkSum.cpp b/src/Network/NeuralNetworkSum.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6bc4be452a1b81b238fc01326c45028626f7ccbd
--- /dev/null
+++ b/src/Network/NeuralNetworkSum.cpp
@@ -0,0 +1,57 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 18.7.18 -
+ */
+
+#include "NeuralNetworkSum.h"
+
+NeuralNetworkSum::NeuralNetworkSum(){
+    this->summand = nullptr;
+    this->summand_coefficient = nullptr;
+}
+
+NeuralNetworkSum::~NeuralNetworkSum() {
+    if( this->summand ){
+        delete this->summand;
+    }
+    if( this->summand_coefficient ){
+        delete this->summand_coefficient;
+    }
+}
+
+void NeuralNetworkSum::add_network(NeuralNetwork *net, double alpha) {
+    if(!this->summand){
+        this->summand = new std::vector<NeuralNetwork*>(0);
+    }
+    this->summand->push_back( net );
+
+    if(!this->summand_coefficient){
+        this->summand_coefficient = new std::vector<double>(0);
+    }
+    this->summand_coefficient->push_back( alpha );
+}
+
+void NeuralNetworkSum::eval_single(std::vector<double> &input, std::vector<double> &output, double *custom_weights) {
+    std::vector<double> mem_output(output.size());
+    std::fill(output.begin(), output.end(), 0.0);
+
+    for(int ni = 0; ni < this->summand->size(); ++ni){
+        this->summand->at(ni)->eval_single(input, mem_output, custom_weights);
+        double alpha = this->summand_coefficient->at(ni);
+        for(int j = 0; j < output.size(); ++j){
+            output[j] += mem_output[j] * alpha;
+        }
+    }
+
+}
+
+size_t NeuralNetworkSum::get_n_weights(){
+    //TODO stupid solution, assumes the networks share weights
+    if(this->summand){
+        return this->summand->at(0)->get_n_weights();
+    }
+
+    return 0;
+}
\ No newline at end of file
diff --git a/src/Network/NeuralNetworkSum.h b/src/Network/NeuralNetworkSum.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c0cba3651a1e87906c1e38948fb9f13c097107b
--- /dev/null
+++ b/src/Network/NeuralNetworkSum.h
@@ -0,0 +1,30 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 18.7.18 -
+ */
+
+#ifndef INC_4NEURO_NEURALNETWORKSUM_H
+#define INC_4NEURO_NEURALNETWORKSUM_H
+
+#include "NeuralNetwork.h"
+
+class NeuralNetworkSum : public NeuralNetwork {
+private:
+    std::vector<NeuralNetwork*> * summand;
+    std::vector<double> * summand_coefficient;
+
+public:
+    NeuralNetworkSum();
+    virtual ~NeuralNetworkSum();
+
+    void add_network( NeuralNetwork *net, double alpha = 1.0 );
+
+    virtual void eval_single(std::vector<double> &input, std::vector<double> &output, double *custom_weights = nullptr) override;
+
+    virtual size_t get_n_weights();
+};
+
+
+#endif //INC_4NEURO_NEURALNETWORKSUM_H
diff --git a/src/Neuron/NeuronLogistic.cpp b/src/Neuron/NeuronLogistic.cpp
index 5113ecd7f2fbc22ab6f8f4308a34ddd20f1b9f8b..dabcb04722378f9c281f5348b68f34333158edae 100644
--- a/src/Neuron/NeuronLogistic.cpp
+++ b/src/Neuron/NeuronLogistic.cpp
@@ -28,8 +28,8 @@ void NeuronLogistic::activate( ) {
     double b = this->activation_function_parameters[1];
     double x = this->potential;
 
-    double ex = std::pow(E, a - x);
-    this->state = std::pow(1.0 + ex, -b);
+    double ex = std::pow(E, b - x);
+    this->state = std::pow(1.0 + ex, -a);
 
 }
 
@@ -40,9 +40,9 @@ double NeuronLogistic::activation_function_get_partial_derivative(int param_idx
     double x = this->potential;
 
     if(param_idx == 0){
-        double ex = std::pow(E, a - x);
+        double ex = std::pow(E, b - x);
 
-        double exa= -std::pow(ex + 1.0, -b);
+        double exa= -std::pow(ex + 1.0, -a);
 
         return exa * std::log(ex + 1.0);
     }
@@ -51,10 +51,10 @@ double NeuronLogistic::activation_function_get_partial_derivative(int param_idx
          * TODO
          * Could be write as activation_function_get_derivative() * -1
          */
-        double ex = std::pow(E, a - x);
-        double ex2 = std::pow(ex + 1.0, -b - 1.0);
+        double ex = std::pow(E, b - x);
+        double ex2 = std::pow(ex + 1.0, -a - 1.0);
 
-        return -b * ex * ex2;
+        return -a * ex * ex2;
     }
 
     return 0.0;
@@ -67,9 +67,9 @@ double NeuronLogistic::activation_function_get_derivative( ) {
     double b = this->activation_function_parameters[1];
     double x = this->potential;
 
-    double ex = std::pow(E, a - x);
-    double ex2 = std::pow(ex + 1.0, -b - 1.0);
+    double ex = std::pow(E, b - x);
+    double ex2 = std::pow(ex + 1.0, -a - 1.0);
 
-    return b * ex * ex2;
+    return a * ex * ex2;
 
 }
\ No newline at end of file
diff --git a/src/net_test_ode_1.cpp b/src/net_test_ode_1.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..70e07407e430c8d0278b1790a2b2cf41c6f901d7
--- /dev/null
+++ b/src/net_test_ode_1.cpp
@@ -0,0 +1,194 @@
+/**
+ * Example solving the following ODE:
+ *
+ * g(t) = (d^2/d^t)y(t) + 4 (d/dt)y(t) + 4y(t) = 0, for t in [0, 4]
+ * y(0) = 1
+ * (d/dt)y(0) = 1
+ *
+ * @author Michal KravÄŤenko
+ * @date 17.7.18 -
+ */
+
+#include <vector>
+#include <utility>
+
+#include "../include/4neuro.h"
+
+int main() {
+
+    int n_inner_neurons = 4;
+
+    std::vector<double> inp, out;
+
+    double d1_s = 0.0, d1_e = 4.0, frac, alpha;
+
+    /* TRAIN DATA FOR g(t) */
+    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_g;
+    unsigned int train_size = 100;
+
+    /* ISOTROPIC TRAIN SET */
+    frac = (d1_e - d1_s) / (train_size - 1);
+    for(unsigned int i = 0; i < train_size; ++i){
+        inp = {frac * i};
+        out = {0.0};
+        data_vec_g.emplace_back(std::make_pair(inp, out));
+    }
+
+    /* CHEBYSCHEV TRAIN SET */
+//    alpha = PI / (train_size - 1);
+//    frac = 0.5 * (d1_e - d1_s);
+//    for(unsigned int i = 0; i < train_size; ++i){
+//        inp = {(std::cos(alpha * i) + 1.0) * frac + d1_s};
+//        out = {0.0};
+//        data_vec_g.emplace_back(std::make_pair(inp, out));
+//    }
+    DataSet ds_g(&data_vec_g);
+
+    /* TRAIN DATA FOR DIRICHLET BC */
+    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_y;
+    inp = {0.0};
+    out = {1.0};
+    data_vec_y.emplace_back(std::make_pair(inp, out));
+    DataSet ds_y(&data_vec_y);
+
+    /* TRAIN DATA FOR NEUMANN BC */
+    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_dy;
+    inp = {0.0};
+    out = {1.0};
+    data_vec_dy.emplace_back(std::make_pair(inp, out));
+    DataSet ds_dy(&data_vec_dy);
+
+    /* DEFINITION OF NNs REPRESENTING y(t) and its derivatives */
+    NeuralNetwork y;
+    NeuralNetwork dy;
+    NeuralNetwork ddy;
+
+    /* Input neurons*/
+    NeuronLinear *i_y = new NeuronLinear(0.0, 1.0);  //f(x) = x
+    NeuronLinear *i_dy = new NeuronLinear(0.0, 1.0);  //f(x) = x
+    NeuronLinear *i_ddy = new NeuronLinear(0.0, 1.0);  //f(x) = x
+
+    /* Output neurons*/
+    NeuronLinear *o_y = new NeuronLinear(0.0, 1.0);  //f(x) = x
+    NeuronLinear *o_dy = new NeuronLinear(0.0, 1.0);  //f(x) = x
+    NeuronLinear *o_ddy = new NeuronLinear(0.0, 1.0);  //f(x) = x
+
+    int idxi = y.add_neuron(i_y);
+    int idxo = y.add_neuron(o_y);
+    int idxinn, idxinn2;
+
+    dy.add_neuron(i_dy);
+    dy.add_neuron(o_dy);
+    ddy.add_neuron(i_ddy);
+    ddy.add_neuron(o_ddy);
+
+    std::vector<double> *weight_ptr;
+
+    weight_ptr = y.get_weights_ptr( );
+
+    dy.set_weight_array( weight_ptr );
+    ddy.set_weight_array( weight_ptr );
+
+    size_t conn_1_idx, conn_2_idx;
+    for( int i = 0; i < n_inner_neurons; ++i ){
+        /* Inner neurons and connections of y(t) */
+        NeuronLogistic *inn_i = new NeuronLogistic(1.0, 0.0); //f(x) = 1.0 / (1.0 + e^(-x))
+        idxinn = y.add_neuron( inn_i );
+
+        conn_1_idx = y.add_connection_simple(idxi, idxinn);
+        conn_2_idx = y.add_connection_simple(idxinn, idxo);
+
+
+        //TODO INCORRECT neurons, must be derivatives of Logistic function
+        /* Inner neurons and connections of dy(t), SHARING WEIGHTS */
+        NeuronLogistic *inn_i_dy = new NeuronLogistic(1.0, 0.0); //f(x) = 1.0 / (1.0 + e^(-x))
+        idxinn = dy.add_neuron( inn_i_dy );
+        dy.add_connection_simple(idxi, idxinn, conn_1_idx);
+        NeuronLinear *identity_dy = new NeuronLinear(0.0, 1.0); //f(x) = x
+        idxinn2 = dy.add_neuron( identity_dy );
+        dy.add_connection_simple(idxinn, idxinn2, conn_2_idx);
+        dy.add_connection_shared_power1(idxinn2, idxo, conn_1_idx);
+
+        /* Inner neurons and connections of ddy(t), SHARING WEIGHTS */
+        NeuronLogistic *inn_i_ddy = new NeuronLogistic(1.0, 0.0); //f(x) = 1.0 / (1.0 + e^(-x))
+        idxinn = ddy.add_neuron( inn_i_ddy );
+        ddy.add_connection_simple(idxi, idxinn, conn_1_idx);
+        NeuronLinear *identity_ddy = new NeuronLinear(0.0, 1.0); //f(x) = x
+        idxinn2 = ddy.add_neuron( identity_ddy );
+        ddy.add_connection_simple(idxinn, idxinn2, conn_2_idx);
+        ddy.add_connection_shared_power2(idxinn2, idxo, conn_1_idx);
+    }
+
+    y.randomize_weights();
+
+    /* SPECIFICATION OF INPUTS */
+    std::vector<size_t> net_input_neurons_indices(1);
+    std::vector<size_t> net_output_neurons_indices(1);
+    net_input_neurons_indices[0] = idxi;
+    net_output_neurons_indices[0] = idxo;
+    y.specify_input_neurons(net_input_neurons_indices);
+    y.specify_output_neurons(net_output_neurons_indices);
+
+    dy.specify_input_neurons(net_input_neurons_indices);
+    dy.specify_output_neurons(net_output_neurons_indices);
+
+    ddy.specify_input_neurons(net_input_neurons_indices);
+    ddy.specify_output_neurons(net_output_neurons_indices);
+
+
+//    y.print_weights();
+//    dy.print_weights();
+//    ddy.print_weights();
+
+//    y.print_stats();
+//    dy.print_stats();
+//    ddy.print_stats();
+
+    /* DEFINITION OF NN SUM */
+    NeuralNetworkSum g;
+    g.add_network(&y, 4.0);
+    g.add_network(&dy, 4.0);
+    g.add_network(&ddy, 1.0);
+
+    /* DEFINITION OF THE PARTIAL ERROR FUNCTIONS */
+    MSE error_y( &y, &ds_y );
+    MSE error_dy( &dy, &ds_dy );
+    MSE error_test_ddy(&ddy, &ds_g);
+    MSE error_g( &g, &ds_g );
+
+    /* DEFINITION OF THE GLOBAL ERROR FUNCTION */
+    MSE_SUM mse_sum;
+    mse_sum.add_error_function( &error_y );
+    mse_sum.add_error_function( &error_dy );
+    mse_sum.add_error_function( &error_g );
+
+//    double weights[2] = {1.0, 1.0};
+//    error_y.eval(weights);
+
+    /* TRAINING METHOD SETUP */
+    unsigned int max_iters = 500;
+
+
+    //must encapsulate each of the partial error functions
+    double *domain_bounds = new double[4 * n_inner_neurons];
+    for(int i = 0; i < 2 * n_inner_neurons; ++i){
+        domain_bounds[2 * i] = -800.0;
+        domain_bounds[2 * i + 1] = 800.0;
+    }
+
+    double c1 = 0.5, c2 = 1.5, w = 0.8;
+
+    unsigned int n_particles = 10;
+
+    ParticleSwarm swarm_01(&mse_sum, domain_bounds, c1, c2, w, n_particles, max_iters);
+//    ParticleSwarm swarm_01(&error_y, domain_bounds, c1, c2, w, n_particles, max_iters);
+//    ParticleSwarm swarm_01(&error_dy, domain_bounds, c1, c2, w, n_particles, max_iters);
+//    ParticleSwarm swarm_01(&error_test_ddy, domain_bounds, c1, c2, w, n_particles, max_iters);
+//    ParticleSwarm swarm_01(&error_g, domain_bounds, c1, c2, w, n_particles, max_iters);
+
+    swarm_01.optimize(0.5, 0.02, 0.9);
+
+
+    delete [] domain_bounds;
+    return 0;
+}
\ No newline at end of file