diff --git a/include/4neuro.h b/include/4neuro.h
index 09f1b997e26a83877b64a7cedb337f629e9640df..2ccc19415ff0a5eb8910422951ce84c14994c881 100644
--- a/include/4neuro.h
+++ b/include/4neuro.h
@@ -32,6 +32,7 @@
 #include "../src/CrossValidator/CrossValidator.h"
 #include "../src/constants.h"
 #include "../src/Coordinates/coordinates.h"
+#include "../src/LearningMethods/NelderMead.h"
 
 // Abbreviate lib4neuro namespace to l4n
 namespace l4n = lib4neuro;
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index c9d0642ef2ebee2d745ad5676726eaf5573e915f..656da507cf722c9bc850d5d2c7bb8801b22fa068 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -119,6 +119,7 @@ IF("${BUILD_LIB}" STREQUAL "yes")
         Reader/XYZReader.cpp
         Reader/Reader.cpp
         Coordinates/coordinates.cpp
+        LearningMethods/NelderMead.cpp
     )
 
     # Detect Threading library
diff --git a/src/LearningMethods/GradientDescentBB.cpp b/src/LearningMethods/GradientDescentBB.cpp
index 971d15efd24bb78dd2a4f209f12933e793189ce8..b1e48e4a6212f33623036d16b3b437c23a6a213e 100644
--- a/src/LearningMethods/GradientDescentBB.cpp
+++ b/src/LearningMethods/GradientDescentBB.cpp
@@ -19,9 +19,7 @@ namespace lib4neuro {
         this->batch             = batch;
     }
 
-    GradientDescentBB::~GradientDescentBB() {
-    }
-
+    GradientDescentBB::~GradientDescentBB() {}
 
     void GradientDescentBB::optimize(lib4neuro::ErrorFunction& ef,
                                      std::ofstream* ofs) {
diff --git a/src/LearningMethods/NelderMead.cpp b/src/LearningMethods/NelderMead.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e0dc15ce4440d282b42caab2d4a4bac75f5adcce
--- /dev/null
+++ b/src/LearningMethods/NelderMead.cpp
@@ -0,0 +1,208 @@
+//
+// Created by martin on 05.09.19.
+//
+
+#include <armadillo>
+#include <utility>
+#include <random>
+
+#include "../message.h"
+
+#include "NelderMead.h"
+
+lib4neuro::NelderMead::NelderMead(double reflectionCoefficient,
+                                  double expansionCoefficient,
+                                  double contractionCoefficient,
+                                  double shrinkCoefficient,
+                                  double targetStandardDeviation,
+                                  double targetMaximalEdgeSize,
+                                  long long int maximumNiters,
+//                                  std::vector<double> initial_position,
+                                  double initial_edge_size)
+    : maximum_niters(maximumNiters), reflection_coefficient(reflectionCoefficient),
+      expansion_coefficient(expansionCoefficient), contraction_coefficient(contractionCoefficient),
+      shrink_coefficient(shrinkCoefficient), target_standard_deviation(targetStandardDeviation),
+      target_maximal_edge_size(targetMaximalEdgeSize),
+      initial_edge_size(initial_edge_size) {}
+
+lib4neuro::NelderMead::NelderMead(double targetStandardDeviation,
+                                  double targetMaximalEdgeSize,
+                                  long long int maximumNiters,
+//                                  std::vector<double> initial_position,
+                                  double initial_edge_size)
+                                  : NelderMead(1,
+                                               2,
+                                               0.5,
+                                               0.5,
+                                               targetStandardDeviation,
+                                               targetMaximalEdgeSize,
+                                               maximumNiters,
+//                                               initial_position,
+                                               initial_edge_size) {}
+
+lib4neuro::NelderMead::NelderMead(double initial_edge_size) : NelderMead(1,
+                                                                         2,
+                                                                         0.5,
+                                                                         0.5,
+                                                                         10e-7,
+                                                                         10e-4,
+                                                                         5000,
+//                                                                         initial_position,
+                                                                         initial_edge_size) {}
+
+long long int lib4neuro::NelderMead::getMaximumNiters() const {
+    return maximum_niters;
+}
+
+double lib4neuro::NelderMead::getReflectionCoefficient() const {
+    return reflection_coefficient;
+}
+
+double lib4neuro::NelderMead::getExpansionCoefficient() const {
+    return expansion_coefficient;
+}
+
+double lib4neuro::NelderMead::getContractionCoefficient() const {
+    return contraction_coefficient;
+}
+
+double lib4neuro::NelderMead::getShrinkCoefficient() const {
+    return shrink_coefficient;
+}
+
+double lib4neuro::NelderMead::getTargetStandardDeviation() const {
+    return target_standard_deviation;
+}
+
+double lib4neuro::NelderMead::getTargetMaximalEdgeSize() const {
+    return target_maximal_edge_size;
+}
+
+lib4neuro::NelderMead::~NelderMead() {}
+
+void lib4neuro::NelderMead::optimize(lib4neuro::ErrorFunction& ef,
+                                     std::ofstream* ofs) {
+    COUT_DEBUG("NelderMead optimizer starting...");
+    COUT_DEBUG("Initial network error: " + std::to_string(ef.eval()));
+    size_t problem_dim = ef.get_dimension();
+    unsigned int best = 0;
+    unsigned int worst = problem_dim;
+    unsigned int second_worst = problem_dim-1;
+
+    /* Initialize simplex */
+    /* Every column corresponds to a coordinates of one point in a simplex */
+    arma::Mat<double> simplex(problem_dim, problem_dim + 1);
+
+    /* Initialize first point coordinates with network's parameters */
+//    this->initial_position = ef.get_parameters();
+
+    arma::Col<double> column(ef.get_parameters());
+    simplex.col(0) = column;
+
+    for(size_t i = 1; i <= problem_dim; i++) {
+        simplex.col(i) = simplex.col(i-1);
+        simplex.col(i)(i-1) += this->initial_edge_size;
+    }
+
+    /* Function values at simplex's points */
+    /* Single output is a single column */
+    arma::Col<double> simplex_function_values(problem_dim+1);
+    std::vector<double> std_vec_tmp;
+
+    /* Indices of a sorted simplex points, according to their function values */
+    arma::Col<arma::uword> simplex_indices;
+
+    std::cout << "init error: " << ef.eval() << std::endl;
+
+    std::vector<double> init_net_parameters = ef.get_parameters();
+
+    this->optimal_parameters = ef.get_parameters();
+
+    for(unsigned int iter_ind = 0; iter_ind < this->maximum_niters; iter_ind++) {
+//        ef.set_parameters(this->optimal_parameters);
+
+        std::cout << "iter ind " << iter_ind << " iter: " << ef.eval() << std::endl;
+        /* Eval function */
+        for(size_t j = 0; j <= problem_dim; j++) {
+            std_vec_tmp = arma::conv_to<std::vector<double>>::from(simplex.col(j));
+            simplex_function_values(j) = ef.eval(&std_vec_tmp);  //TODO do in a more efficient way
+            ef.set_parameters(this->optimal_parameters);
+        }
+
+        std::cout << ef.eval() << std::endl;
+
+        /* Sort simplex points according to the error function output */
+        simplex_indices = arma::sort_index(simplex_function_values);
+        simplex = simplex.cols(simplex_indices);
+        simplex_function_values = arma::sort(simplex_function_values);
+
+        this->optimal_parameters = arma::conv_to<std::vector<double>>::from(simplex.col(best));
+        ef.set_parameters(this->optimal_parameters);
+
+        /* Termination condition */
+
+
+        /* Compute centroid */
+        arma::Col<double> centroid(problem_dim);
+        for(size_t j = 0; j < problem_dim; j++) {
+            centroid(j) = arma::mean(simplex.row(j).subvec(0, problem_dim-1));
+        }
+
+        /* Reflection */
+        arma::Col<double> xr = centroid + this->reflection_coefficient*(centroid - simplex.col(worst));  // Reflection point
+        std_vec_tmp = arma::conv_to<std::vector<double>>::from(xr);
+        double xr_val = ef.eval(&std_vec_tmp);
+        ef.set_parameters(this->optimal_parameters);
+
+        if(xr_val < simplex_function_values(best) ) {
+            /* Expansion */
+            arma::Col<double> xe = centroid + this->expansion_coefficient*(xr - centroid);
+            std_vec_tmp = arma::conv_to<std::vector<double>>::from(xe);
+            double xe_val = ef.eval(&std_vec_tmp);
+            ef.set_parameters(this->optimal_parameters);
+
+            if(xe_val < xr_val) {
+//                std::cout << simplex_function_values(problem_dim) << " " << xe_val << std::endl;
+                simplex.col(worst) = xe;
+                continue;
+            } else {
+//                std::cout << simplex_function_values(problem_dim) << " " << xr_val << std::endl;
+                simplex.col(worst) = xr;
+                continue;
+            }
+
+        } else if(xr_val < simplex_function_values(second_worst)) {
+            /* Reflection */
+            simplex.col(worst) = xr;
+            continue;
+
+        } else {
+            /* Contraction */
+            arma::Col<double> xc = centroid + this->contraction_coefficient*(simplex.col(worst) - centroid);
+            std_vec_tmp = arma::conv_to<std::vector<double>>::from(xc);
+            double xc_value = ef.eval(&std_vec_tmp);
+            ef.set_parameters(this->optimal_parameters);
+
+            if(xc_value < simplex_function_values(worst)) {
+                simplex.col(worst) = xc;
+                continue;
+            } else {
+                /* Shrink */
+                for(size_t j = 0; j <= problem_dim; j++) {
+                    if(j == best) {
+                        continue;
+                    }
+
+                    simplex.col(j) = centroid + this->shrink_coefficient*(simplex.col(j) - simplex.col(best));
+                }
+            }
+        }
+    }
+
+    std::cout << "end" << std::endl;
+
+//    std_vec_tmp = arma::conv_to<std::vector<double>>::from(simplex.col(best));
+//    ef.set_parameters(std_vec_tmp);
+    this->optimal_parameters = arma::conv_to<std::vector<double>>::from(simplex.col(best));
+    ef.set_parameters(this->optimal_parameters);
+}
diff --git a/src/LearningMethods/NelderMead.h b/src/LearningMethods/NelderMead.h
new file mode 100644
index 0000000000000000000000000000000000000000..830bb8ec88d0d3e628ef9af5b8493b99a79bdf7b
--- /dev/null
+++ b/src/LearningMethods/NelderMead.h
@@ -0,0 +1,90 @@
+//
+// Created by martin on 05.09.19.
+//
+
+#ifndef LIB4NEURO_NELDERMEAD_H
+#define LIB4NEURO_NELDERMEAD_H
+
+#include <vector>
+
+#include "LearningMethod.h"
+
+
+namespace lib4neuro {
+    class NelderMead : public LearningMethod {
+    protected:
+        /**
+         * Maximal number of iterations - optimization will stop after that, even if not converged
+         */
+        unsigned int maximum_niters;
+
+        double reflection_coefficient;
+
+        double expansion_coefficient;
+
+        double contraction_coefficient;
+
+        double shrink_coefficient;
+
+        /**
+         * Threshold for a standard deviation of simplex function values - serving as a part of a terminating condition
+         */
+        double target_standard_deviation;
+
+        /**
+         * Threshold for a maximal edge size of a simplex - serving as a part of a terminating condition
+         */
+        double target_maximal_edge_size;
+
+//        /**
+//         * Initial coordinates of the first point in the simplex
+//         */
+//        std::vector<double> initial_position;
+
+        /**
+         * Initial size of the simplexes' edges
+         */
+         double initial_edge_size;
+
+    public:
+        LIB4NEURO_API NelderMead(double reflectionCoefficient,
+                                 double expansionCoefficient,
+                                 double contractionCoefficient,
+                                 double shrinkCoefficient,
+                                 double targetStandardDeviation,
+                                 double targetMaximalEdgeSize,
+                                 long long int maximumNiters,
+//                                 std::vector<double> initial_position,
+                                 double initial_edge_size);
+
+        LIB4NEURO_API NelderMead(double targetStandardDeviation,
+                                 double targetMaximalEdgeSize,
+                                 long long int maximumNiters,
+//                                 std::vector<double> initial_position,
+                                 double initial_edge_size);
+
+        LIB4NEURO_API explicit NelderMead(double initial_edge_size);
+
+        LIB4NEURO_API virtual ~NelderMead();
+
+        // TODO rewrite method names to snake case
+        LIB4NEURO_API [[nodiscard]] long long int getMaximumNiters() const;
+
+        LIB4NEURO_API [[nodiscard]] double getReflectionCoefficient() const;
+
+        LIB4NEURO_API [[nodiscard]] double getExpansionCoefficient() const;
+
+        LIB4NEURO_API [[nodiscard]] double getContractionCoefficient() const;
+
+        LIB4NEURO_API [[nodiscard]] double getShrinkCoefficient() const;
+
+        LIB4NEURO_API [[nodiscard]] double getTargetStandardDeviation() const;
+
+        LIB4NEURO_API [[nodiscard]] double getTargetMaximalEdgeSize() const;
+
+        LIB4NEURO_API void optimize(lib4neuro::ErrorFunction& ef,
+                                    std::ofstream* ofs = nullptr) override;
+    };
+}
+
+#endif //LIB4NEURO_NELDERMEAD_H
diff --git a/src/LearningMethods/ParticleSwarm.cpp b/src/LearningMethods/ParticleSwarm.cpp
index c8528554ffd7ffc7747aa92f6d4c24e699071b68..7ba249b27dc759c361514a24067df32a0afb9e87 100644
--- a/src/LearningMethods/ParticleSwarm.cpp
+++ b/src/LearningMethods/ParticleSwarm.cpp
@@ -292,8 +292,6 @@ namespace lib4neuro {
      */
     void ParticleSwarm::optimize(lib4neuro::ErrorFunction& ef,
                                  std::ofstream* ofs) {
-        //TODO add output to the 'ofs'
-
         COUT_INFO("Finding optima via Globalized Particle Swarm method..." << std::endl);
         if (ofs && ofs->is_open()) {
             *ofs << "Finding optima via Globalized Particle Swarm method..." << std::endl;
diff --git a/src/examples/dev_sandbox.cpp b/src/examples/dev_sandbox.cpp
index 624aab954b3ab053912a3b17854af10770b61d0c..36d919476430c5d970fb47c0e0745de88c50ce99 100644
--- a/src/examples/dev_sandbox.cpp
+++ b/src/examples/dev_sandbox.cpp
@@ -1,6 +1,8 @@
 //
 // Created by martin on 20.08.19.
 //
+
+#include <armadillo>
 #include <exception>
 
 #include <4neuro.h>
@@ -107,6 +109,21 @@ double optimize_via_LBMQ(l4n::NeuralNetwork& net,
 	return err;
 }
 
+double optimize_via_NelderMead(l4n::NeuralNetwork& net, l4n::ErrorFunction& ef) {
+    l4n::NelderMead nm(10e-4, 10e-7, 10, 100);
+
+    nm.optimize(ef);
+    net.copy_parameter_space(nm.get_parameters());
+
+    /* ERROR CALCULATION */
+    double err = ef.eval(nullptr);
+     std::cout << "Run finished! Error of the network[Nelder-Mead]: " << err << std::endl;
+
+    /* Just for validation test purposes - NOT necessary for the example to work! */
+    return err;
+
+}
+
 int main() {
 
     try{
@@ -142,7 +159,7 @@ int main() {
         elements[l4n::ELEMENT_SYMBOL::He] = &helium;
 
         /* Read data */
-        l4n::XYZReader reader("/home/bes0030/HE4+T0.xyz", true);
+        l4n::XYZReader reader("/home/martin/Desktop/HE4+T0.xyz", true);
         reader.read();
 
         std::cout << "Finished reading data" << std::endl;
@@ -160,12 +177,15 @@ int main() {
 
         l4n::MSE mse(&net, ds.get());
 
-
         net.randomize_parameters();
-         optimize_via_particle_swarm(net, mse);
+//         optimize_via_particle_swarm(net, mse);
+
+
+        optimize_via_NelderMead(net, mse);
+
 //        double err1 = optimize_via_LBMQ(net, mse);
         double err2 = optimize_via_gradient_descent(net, mse);
-		
+
         std::cout << "Weights: " << net.get_min_max_weight().first << " " << net.get_min_max_weight().second << std::endl;
 
         /* Print fit comparison with real data */
@@ -189,5 +209,12 @@ int main() {
         exit(EXIT_FAILURE);
     }
 
+//    arma::Mat<double> m = {{1,2,3}, {4,5,6}, {7,8,9}};
+//    arma::Col<double> v = arma::conv_to<std::vector<double>(m);
+
+
+
+//    std::cout << arma::stddev(m) << std::endl;
+
     return 0;
 }