diff --git a/.gitmodules b/.gitmodules
index 25716f78..427f4910 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -10,3 +10,6 @@
 [submodule "lib/cxxopts"]
 	path = lib/cxxopts
 	url = https://github.com/tud-zih-energy/FIRESTARTER.git
+[submodule "lib/SOT"]
+	path = lib/SOT
+	url = https://github.com/dme65/SOT.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ba6bf49d..f891110d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -73,15 +73,11 @@ if(NOT DEFINED ASMJIT_STATIC)
 	set(ASMJIT_STATIC TRUE)
 endif()
 
-add_subdirectory(lib/asmjit)
-add_subdirectory(lib/nitro)
+add_subdirectory(lib)
 
 include_directories(include)
 include_directories(lib/cxxopts/include)
-
-set(JSON_BuildTests OFF CACHE INTERNAL "")
-set(JSON_Install OFF CACHE INTERNAL "")
-add_subdirectory(lib/json)
+include_directories(lib/SOT/include)
 
 set(THREADS_PREFER_PTHREAD_FLAG ON)
 find_package(Threads REQUIRED)
diff --git a/derivation.nix b/derivation.nix
index 6c1215d1..cd7ab8eb 100644
--- a/derivation.nix
+++ b/derivation.nix
@@ -6,6 +6,7 @@
 , git
 , pkgconfig
 , cudatoolkit
+, armadillo
 , withCuda ? false
 , linuxPackages
 }:
@@ -53,14 +54,16 @@ stdenv.mkDerivation rec {
   nativeBuildInputs = [ cmake git pkgconfig ];
 
   buildInputs = if withCuda then
-    [ glibc_multi cudatoolkit linuxPackages.nvidia_x11 hwloc ]
+    [ glibc_multi cudatoolkit linuxPackages.nvidia_x11 hwloc armadillo ]
     else
-    [ glibc.static hwloc ];
+    [ glibc_multi hwloc armadillo ];
+#    [ glibc.static hwloc armadillo ];
 
   cmakeFlags = [
     "-DFIRESTARTER_BUILD_HWLOC=OFF"
     "-DCMAKE_C_COMPILER_WORKS=1"
     "-DCMAKE_CXX_COMPILER_WORKS=1"
+    "-DFIRESTARTER_LINK_STATIC=OFF"
   ] ++ lib.optionals withCuda [
    "-DFIRESTARTER_BUILD_TYPE=FIRESTARTER_CUDA"
   ];
diff --git a/include/firestarter/Firestarter.hpp b/include/firestarter/Firestarter.hpp
index 9d38799a..a34a05fa 100644
--- a/include/firestarter/Firestarter.hpp
+++ b/include/firestarter/Firestarter.hpp
@@ -83,8 +83,9 @@ class Firestarter {
               std::string const &optimizationAlgorithm,
               std::vector<std::string> const &optimizationMetrics,
               std::chrono::seconds const &evaluationDuration,
-              unsigned individuals, std::string const &optimizeOutfile,
-              unsigned generations, double nsga2_cr, double nsga2_m);
+              std::string const &optimizeOutfile, unsigned nsga2_individuals,
+              unsigned nsga2_generations, double nsga2_cr, double nsga2_m,
+              unsigned samo_is_individuals, unsigned samo_is_maxEvaluations);
 
   ~Firestarter();
 
@@ -113,11 +114,14 @@ class Firestarter {
   const std::string _optimizationAlgorithm;
   const std::vector<std::string> _optimizationMetrics;
   const std::chrono::seconds _evaluationDuration;
-  const unsigned _individuals;
   const std::string _optimizeOutfile;
-  const unsigned _generations;
+  unsigned _individuals;
+  const unsigned _nsga2_individuals;
+  const unsigned _nsga2_generations;
   const double _nsga2_cr;
   const double _nsga2_m;
+  const unsigned _samo_is_individuals;
+  const unsigned _samo_is_maxEvaluations;
 
 #ifndef FIRESTARTER_BUILD_CUDA_ONLY
 #if defined(__i386__) || defined(_M_IX86) || defined(__x86_64__) ||            \
diff --git a/include/firestarter/Optimizer/Algorithm/SAMO_IS.hpp b/include/firestarter/Optimizer/Algorithm/SAMO_IS.hpp
new file mode 100644
index 00000000..6c0bfff8
--- /dev/null
+++ b/include/firestarter/Optimizer/Algorithm/SAMO_IS.hpp
@@ -0,0 +1,27 @@
+#pragma once
+
+#include <firestarter/Optimizer/Algorithm.hpp>
+
+namespace firestarter::optimizer::algorithm {
+
+class SAMO_IS : public Algorithm {
+public:
+  SAMO_IS(unsigned maxEvaluations, unsigned nsga2_individuals,
+          unsigned nsga2_generations, double nsga2_cr, double nsga2_m);
+  ~SAMO_IS() {}
+
+  void checkPopulation(firestarter::optimizer::Population const &pop,
+                       std::size_t populationSize) override;
+
+  firestarter::optimizer::Population
+  evolve(firestarter::optimizer::Population &pop) override;
+
+private:
+  unsigned _maxEvaluations;
+  unsigned _nsga2_individuals;
+  unsigned _nsga2_generations;
+  double _nsga2_cr;
+  double _nsga2_m;
+};
+
+} // namespace firestarter::optimizer::algorithm
diff --git a/include/firestarter/Optimizer/History.hpp b/include/firestarter/Optimizer/History.hpp
index 7dfb3887..e85fe562 100644
--- a/include/firestarter/Optimizer/History.hpp
+++ b/include/firestarter/Optimizer/History.hpp
@@ -74,6 +74,16 @@ struct History {
       _f = {};
 
 public:
+  inline static std::size_t size() { return _x.size(); }
+
+  inline static std::vector<Individual> const &x() { return _x; }
+
+  inline static std::vector<
+      std::map<std::string, firestarter::measurement::Summary>> const &
+  f() {
+    return _f;
+  }
+
   inline static void append(
       std::vector<unsigned> const &ind,
       std::map<std::string, firestarter::measurement::Summary> const &metric) {
diff --git a/include/firestarter/Optimizer/Population.hpp b/include/firestarter/Optimizer/Population.hpp
index 7a9ff579..fef199a7 100644
--- a/include/firestarter/Optimizer/Population.hpp
+++ b/include/firestarter/Optimizer/Population.hpp
@@ -40,16 +40,18 @@ class Population {
   // Construct a population from a problem.
   Population() = default;
 
-  Population(std::shared_ptr<Problem> &&problem)
-      : _problem(std::move(problem)), gen(rd()) {}
+  Population(std::shared_ptr<Problem> &&problem, bool saveHistory = true)
+      : _problem(std::move(problem)), _saveHistory(saveHistory), gen(rd()) {}
 
   Population(Population &pop)
-      : _problem(pop._problem), _x(pop._x), _f(pop._f), gen(rd()) {}
+      : _problem(pop._problem), _x(pop._x), _f(pop._f),
+        _saveHistory(pop._saveHistory), gen(rd()) {}
 
   Population &operator=(Population const &pop) {
     _problem = std::move(pop._problem);
     _x = pop._x;
     _f = pop._f;
+    _saveHistory = pop._saveHistory;
     gen = pop.gen;
 
     return *this;
@@ -89,6 +91,8 @@ class Population {
   std::vector<Individual> _x;
   std::vector<std::vector<double>> _f;
 
+  bool _saveHistory;
+
   std::random_device rd;
   std::mt19937 gen;
 };
diff --git a/include/firestarter/Optimizer/Problem.hpp b/include/firestarter/Optimizer/Problem.hpp
index f88b0bc3..cb8745b0 100644
--- a/include/firestarter/Optimizer/Problem.hpp
+++ b/include/firestarter/Optimizer/Problem.hpp
@@ -42,7 +42,27 @@ class Problem {
 
   virtual std::vector<double>
   fitness(std::map<std::string, firestarter::measurement::Summary> const
-              &summaries) = 0;
+              &summaries) {
+    std::vector<double> values = {};
+
+    for (auto const &metricName : this->metrics()) {
+      auto findName = [metricName](auto const &summary) {
+        return metricName.compare(summary.first) == 0;
+      };
+
+      auto it = std::find_if(summaries.begin(), summaries.end(), findName);
+
+      if (it == summaries.end()) {
+        continue;
+      }
+
+      // round to two decimal places after the comma
+      auto value = std::round(it->second.average * 100.0) / 100.0;
+      values.push_back(value);
+    }
+
+    return values;
+  }
 
   // get the bounds of the problem
   virtual std::vector<std::tuple<unsigned, unsigned>> getBounds() const = 0;
@@ -59,6 +79,8 @@ class Problem {
   // get the number of fitness evaluations
   unsigned long long getFevals() const { return _fevals; };
 
+  virtual std::vector<std::string> metrics() const = 0;
+
 protected:
   // number of fitness evaluations
   unsigned long long _fevals;
diff --git a/include/firestarter/Optimizer/Problem/CLIArgumentProblem.hpp b/include/firestarter/Optimizer/Problem/CLIArgumentProblem.hpp
index 9e33b7c4..9df7aa39 100644
--- a/include/firestarter/Optimizer/Problem/CLIArgumentProblem.hpp
+++ b/include/firestarter/Optimizer/Problem/CLIArgumentProblem.hpp
@@ -87,38 +87,6 @@ class CLIArgumentProblem final : public firestarter::optimizer::Problem {
     return _measurementWorker->getValues(_startDelta, _stopDelta);
   }
 
-  std::vector<double> fitness(
-      std::map<std::string, firestarter::measurement::Summary> const &summaries)
-      override {
-    std::vector<double> values = {};
-
-    for (auto const &metricName : _metrics) {
-      auto findName = [metricName](auto const &summary) {
-        auto invertedName = "-" + summary.first;
-        return metricName.compare(summary.first) == 0 ||
-               metricName.compare(invertedName) == 0;
-      };
-
-      auto it = std::find_if(summaries.begin(), summaries.end(), findName);
-
-      if (it == summaries.end()) {
-        continue;
-      }
-
-      // round to two decimal places after the comma
-      auto value = std::round(it->second.average * 100.0) / 100.0;
-
-      // invert metric
-      if (metricName[0] == '-') {
-        value *= -1.0;
-      }
-
-      values.push_back(value);
-    }
-
-    return values;
-  }
-
   // get the bounds of the problem
   std::vector<std::tuple<unsigned, unsigned>> getBounds() const override {
     std::vector<std::tuple<unsigned, unsigned>> vec(
@@ -130,6 +98,8 @@ class CLIArgumentProblem final : public firestarter::optimizer::Problem {
   // get the number of objectives.
   std::size_t getNobjs() const override { return _metrics.size(); }
 
+  std::vector<std::string> metrics() const override { return _metrics; }
+
 private:
   std::function<void(std::vector<std::pair<std::string, unsigned>> const &)>
       _changePayloadFunction;
diff --git a/include/firestarter/Optimizer/Problem/SurrogateProblem.hpp b/include/firestarter/Optimizer/Problem/SurrogateProblem.hpp
new file mode 100644
index 00000000..bde6153f
--- /dev/null
+++ b/include/firestarter/Optimizer/Problem/SurrogateProblem.hpp
@@ -0,0 +1,135 @@
+#pragma once
+
+#include <firestarter/Optimizer/Problem.hpp>
+#include <firestarter/Optimizer/Surrogate/SurrogateSelector.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <functional>
+#include <memory>
+#include <thread>
+#include <tuple>
+#include <utility>
+
+namespace firestarter::optimizer::problem {
+
+class SurrogateProblem final : public firestarter::optimizer::Problem {
+public:
+  SurrogateProblem(std::vector<std::string> const &metrics,
+                   std::vector<std::tuple<unsigned, unsigned>> const &bounds)
+      : _metrics(metrics), _bounds(bounds) {
+    assert(_metrics.size() != 0);
+
+    auto const &hist_x = History::x();
+    auto const dims = hist_x[0].size();
+
+    arma::vec boundsLow(dims);
+    arma::vec boundsUp(dims);
+    arma::mat x(dims, hist_x.size());
+
+    {
+      for (std::size_t i = 0; i < dims; ++i) {
+        boundsLow(i) = std::get<0>(bounds[i]);
+        boundsUp(i) = std::get<1>(bounds[i]);
+      }
+
+      for (std::size_t i = 0; i < hist_x.size(); ++i) {
+        x.col(i) = arma::conv_to<arma::vec>::from(hist_x[i]);
+      }
+    }
+
+    for (auto const &metric : metrics) {
+      std::string strippedMetricName;
+      arma::vec y(hist_x.size());
+
+      if (metric[0] == '-') {
+        strippedMetricName = metric.substr(1);
+      } else {
+        strippedMetricName = metric;
+      }
+
+      for (std::size_t i = 0; i < hist_x.size(); ++i) {
+        // fill y with 3 digits precision
+        y(i) = (double)((int)(History::find(hist_x[i])
+                                  .value()[strippedMetricName]
+                                  .average *
+                              1000)) /
+               1000.0;
+      }
+
+      auto model = std::make_unique<
+          firestarter::optimizer::surrogate::SurrogateSelector>(boundsLow,
+                                                                boundsUp, x, y);
+      log::info() << "Using surrogate model " << model->name() << " for metric "
+                  << metric;
+      _models.push_back(std::move(model));
+    }
+  }
+
+  ~SurrogateProblem() {}
+
+  // return all available metrics for the individual
+  std::map<std::string, firestarter::measurement::Summary>
+  metrics(std::vector<unsigned> const &individual) override {
+    std::map<std::string, firestarter::measurement::Summary> metrics = {};
+    for (std::size_t i = 0; i < _metrics.size(); ++i) {
+      auto name = _metrics[i];
+      auto value = _models[i]->eval(arma::conv_to<arma::vec>::from(individual));
+      firestarter::measurement::Summary summary;
+      summary.average = value;
+      metrics[name] = summary;
+    }
+    return metrics;
+  }
+
+  std::vector<double> fitness(
+      std::map<std::string, firestarter::measurement::Summary> const &summaries)
+      override {
+    std::vector<double> values = {};
+
+    for (auto const &metricName : _metrics) {
+      auto findName = [metricName](auto const &summary) {
+        auto invertedName = "-" + summary.first;
+        return metricName.compare(summary.first) == 0 ||
+               metricName.compare(invertedName) == 0;
+      };
+
+      auto it = std::find_if(summaries.begin(), summaries.end(), findName);
+
+      if (it == summaries.end()) {
+        continue;
+      }
+
+      // round to two decimal places after the comma
+      auto value = std::round(it->second.average * 100.0) / 100.0;
+
+      // invert metric
+      if (metricName[0] == '-') {
+        value *= -1.0;
+      }
+
+      values.push_back(value);
+    }
+
+    return values;
+  }
+
+  // get the bounds of the problem
+  std::vector<std::tuple<unsigned, unsigned>> getBounds() const override {
+    return _bounds;
+  }
+
+  // get the number of objectives.
+  std::size_t getNobjs() const override { return _metrics.size(); }
+
+  std::vector<std::string> metrics() const override { return _metrics; }
+
+private:
+  std::vector<std::string> _metrics;
+  std::vector<std::tuple<unsigned, unsigned>> _bounds;
+  std::vector<
+      std::unique_ptr<firestarter::optimizer::surrogate::SurrogateSelector>>
+      _models;
+};
+
+} // namespace firestarter::optimizer::problem
diff --git a/include/firestarter/Optimizer/Surrogate/Model/SOT.hpp b/include/firestarter/Optimizer/Surrogate/Model/SOT.hpp
new file mode 100644
index 00000000..96ee8e7f
--- /dev/null
+++ b/include/firestarter/Optimizer/Surrogate/Model/SOT.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+#include <firestarter/Logging/Log.hpp>
+#include <firestarter/Optimizer/Surrogate/SurrogateModel.hpp>
+
+#include <sot.h>
+
+namespace firestarter::optimizer::surrogate::model {
+
+template <class T> class SOT : public SurrogateModel {
+  static_assert(std::is_base_of<sot::Surrogate, T>::value,
+                "T must extend sot::Surrogate");
+
+public:
+  // train and select the fitting surrogate model
+  SOT(arma::vec const &boundsLow, arma::vec const &boundsUp, arma::mat const &x,
+      arma::vec const &y) {
+    assert(boundsLow.n_elem == boundsUp.n_elem);
+    assert(boundsLow.n_elem == x.n_rows);
+    assert(x.n_cols == y.n_elem);
+
+    _name = typeid(T).name();
+    _model =
+        std::make_unique<T>(y.n_elem, boundsLow.n_elem, boundsLow, boundsUp);
+
+    _model->addPoints(x, y);
+
+    _model->fit();
+  }
+
+  ~SOT() {}
+
+  // eval the selected surrogate model
+  double eval(arma::vec const &x) override { return _model->eval(x); }
+
+private:
+  std::unique_ptr<T> _model;
+};
+
+} // namespace firestarter::optimizer::surrogate::model
diff --git a/include/firestarter/Optimizer/Surrogate/SurrogateModel.hpp b/include/firestarter/Optimizer/Surrogate/SurrogateModel.hpp
new file mode 100644
index 00000000..15f08ad1
--- /dev/null
+++ b/include/firestarter/Optimizer/Surrogate/SurrogateModel.hpp
@@ -0,0 +1,24 @@
+#pragma once
+
+#include <armadillo>
+
+namespace firestarter::optimizer::surrogate {
+
+class SurrogateModel {
+public:
+  // train and select the fitting surrogate model
+  SurrogateModel() {}
+
+  virtual ~SurrogateModel() {}
+
+  // get the name of the selected surrogate model
+  std::string const &name() const { return _name; }
+
+  // eval the selected surrogate model
+  virtual double eval(arma::vec const &x) = 0;
+
+protected:
+  std::string _name;
+};
+
+} // namespace firestarter::optimizer::surrogate
diff --git a/include/firestarter/Optimizer/Surrogate/SurrogateSelector.hpp b/include/firestarter/Optimizer/Surrogate/SurrogateSelector.hpp
new file mode 100644
index 00000000..b61780b4
--- /dev/null
+++ b/include/firestarter/Optimizer/Surrogate/SurrogateSelector.hpp
@@ -0,0 +1,101 @@
+#pragma once
+
+#include <firestarter/Optimizer/Surrogate/Model/SOT.hpp>
+#include <firestarter/Optimizer/Surrogate/SurrogateModel.hpp>
+
+#include <armadillo>
+
+#define REGISTER(NAME...)                                                      \
+  [](arma::vec const &boundsLow, arma::vec const &boundsUp,                    \
+     arma::mat const &x, arma::vec const &y) -> SurrogateModel * {             \
+    return new model::NAME(boundsLow, boundsUp, x, y);                         \
+  }
+
+namespace firestarter::optimizer::surrogate {
+
+class SurrogateSelector {
+public:
+  // train and select the fitting surrogate model
+  SurrogateSelector(arma::vec const &boundsLow, arma::vec const &boundsUp,
+                    arma::mat const &x, arma::vec const &y) {
+    arma::mat x_train, x_eval;
+    arma::vec y_train, y_eval_true;
+
+    // split the dataset into train and validation
+    arma::vec train_idx;
+    arma::vec eval_idx;
+    {
+      // random shuffle of indicies from 0 to y.n_elem - 1
+      arma::uvec shuffle = arma::randperm(y.n_elem);
+      // get first 90% for train and last 10% for eval
+      auto split_point = std::ceil((double)shuffle.n_elem * 0.9);
+
+      arma::uvec train_idx = shuffle.head(split_point);
+      arma::uvec eval_idx = shuffle.tail(shuffle.n_elem - split_point);
+
+      if (eval_idx.n_elem == 0) {
+        assert((false, "No elements left for evaluation of surrogate models in "
+                       "__FILE__ __LINE__"));
+      }
+
+      x_train = x.cols(train_idx);
+      x_eval = x.cols(eval_idx);
+      y_train = y(train_idx);
+      y_eval_true = y(eval_idx);
+    }
+
+    // train the models
+    std::vector<std::unique_ptr<SurrogateModel>> models;
+    for (auto const &ctor : _modelsCtor) {
+      models.push_back(std::move(std::unique_ptr<SurrogateModel>(
+          ctor(boundsLow, boundsUp, x_train, y_train))));
+    }
+
+    // vector of mean squared error
+    arma::vec mse(models.size());
+
+    for (std::size_t i = 0; i < models.size(); ++i) {
+      auto &model = models[i];
+      arma::vec y_eval(y_eval_true.n_elem);
+
+      for (std::size_t j = 0; j < x_eval.n_cols; ++j) {
+        arma::vec x = x_eval.col(j);
+        y_eval(j) = model->eval(x);
+      }
+
+      arma::vec diff = y_eval - y_eval_true;
+      mse(i) = sum(diff % diff) / (double)y_eval.n_elem;
+
+      log::info() << "mean squared error of " << model->name() << " = "
+                  << mse(i);
+    }
+
+    auto idx = arma::index_min(mse);
+
+    _model = std::move(models[idx]);
+  }
+
+  ~SurrogateSelector() {}
+
+  // get the name of the selected surrogate model
+  std::string const &name() const { return _model->name(); }
+
+  // eval the selected surrogate model
+  double eval(arma::vec const &x) { return _model->eval(x); }
+
+private:
+  // list of surrogate ctors
+  const std::vector<std::function<SurrogateModel *(
+      arma::vec const &boundsLow, arma::vec const &boundsUp, arma::mat const &x,
+      arma::vec const &y)>>
+      _modelsCtor = {
+          REGISTER(SOT<sot::CubicRBF>), REGISTER(SOT<sot::TpsRBF>),
+          REGISTER(
+              SOT<sot::RBFInterpolant<sot::LinearKernel, sot::ConstantTail>>)};
+
+#undef REGISTER
+
+  std::unique_ptr<SurrogateModel> _model;
+};
+
+} // namespace firestarter::optimizer::surrogate
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
new file mode 100644
index 00000000..16bccae4
--- /dev/null
+++ b/lib/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_subdirectory(asmjit)
+
+add_subdirectory(nitro)
+
+set(JSON_BuildTests OFF CACHE INTERNAL "")
+set(JSON_Install OFF CACHE INTERNAL "")
+add_subdirectory(json)
diff --git a/lib/SOT b/lib/SOT
new file mode 160000
index 00000000..269031ab
--- /dev/null
+++ b/lib/SOT
@@ -0,0 +1 @@
+Subproject commit 269031ab32680e1d59dc3480d0602f5bc1969244
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ee7e92a1..6fdb4fd0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -38,6 +38,7 @@ if(CMAKE_SYSTEM_NAME STREQUAL "Linux")
 		firestarter/Optimizer/OptimizerWorker.cpp
 		firestarter/Optimizer/Util/MultiObjective.cpp
 		firestarter/Optimizer/Algorithm/NSGA2.cpp
+		firestarter/Optimizer/Algorithm/SAMO_IS.cpp
 		)
 endif()
 
@@ -52,6 +53,8 @@ if ("${FIRESTARTER_BUILD_TYPE}" STREQUAL "FIRESTARTER_CUDA")
 	find_package(CUDA REQUIRED)
 	include_directories(${CUDA_INCLUDE_DIRS})
 
+	find_package(Armadillo REQUIRED)
+
 	SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DFIRESTARTER_BUILD_CUDA")
 
 	add_executable(FIRESTARTER_CUDA
@@ -68,6 +71,7 @@ if ("${FIRESTARTER_BUILD_TYPE}" STREQUAL "FIRESTARTER_CUDA")
 
 	target_link_libraries(FIRESTARTER_CUDA
 		hwloc
+		${ARMADILLO_LIBRARIES}
 		AsmJit::AsmJit
 		Nitro::log
 		nlohmann_json::nlohmann_json
@@ -124,6 +128,8 @@ elseif(${FIRESTARTER_BUILD_TYPE} STREQUAL "FIRESTARTER")
 		)
 	target_compile_features(FIRESTARTER PRIVATE cxx_std_17)
 
+	find_package(Armadillo REQUIRED)
+
 	# static linking is not supported on Darwin, see Apple Technical QA1118
 	if(CMAKE_SYSTEM_NAME STREQUAL "Darwin")
 		find_library( COREFOUNDATION_LIBRARY CoreFoundation )
@@ -165,6 +171,7 @@ elseif(${FIRESTARTER_BUILD_TYPE} STREQUAL "FIRESTARTER")
 
 	target_link_libraries(FIRESTARTER
 		hwloc
+		${ARMADILLO_LIBRARIES}
 		AsmJit::AsmJit
 		Nitro::log
 		nlohmann_json::nlohmann_json
diff --git a/src/firestarter/Firestarter.cpp b/src/firestarter/Firestarter.cpp
index 5b681d6a..efbf86fb 100644
--- a/src/firestarter/Firestarter.cpp
+++ b/src/firestarter/Firestarter.cpp
@@ -24,6 +24,7 @@
 #ifndef FIRESTARTER_BUILD_CUDA_ONLY
 #if defined(linux) || defined(__linux__)
 #include <firestarter/Optimizer/Algorithm/NSGA2.hpp>
+#include <firestarter/Optimizer/Algorithm/SAMO_IS.hpp>
 #include <firestarter/Optimizer/History.hpp>
 #include <firestarter/Optimizer/Problem/CLIArgumentProblem.hpp>
 extern "C" {
@@ -61,9 +62,10 @@ Firestarter::Firestarter(
     std::chrono::seconds const &preheat,
     std::string const &optimizationAlgorithm,
     std::vector<std::string> const &optimizationMetrics,
-    std::chrono::seconds const &evaluationDuration, unsigned individuals,
-    std::string const &optimizeOutfile, unsigned generations, double nsga2_cr,
-    double nsga2_m)
+    std::chrono::seconds const &evaluationDuration,
+    std::string const &optimizeOutfile, unsigned nsga2_individuals,
+    unsigned nsga2_generations, double nsga2_cr, double nsga2_m,
+    unsigned samo_is_individuals, unsigned samo_is_maxEvaluations)
     : _argc(argc), _argv(argv), _timeout(timeout), _loadPercent(loadPercent),
       _period(period), _dumpRegisters(dumpRegisters),
       _dumpRegistersTimeDelta(dumpRegistersTimeDelta),
@@ -74,9 +76,11 @@ Firestarter::Firestarter(
       _stopDelta(stopDelta), _measurement(measurement), _optimize(optimize),
       _preheat(preheat), _optimizationAlgorithm(optimizationAlgorithm),
       _optimizationMetrics(optimizationMetrics),
-      _evaluationDuration(evaluationDuration), _individuals(individuals),
-      _optimizeOutfile(optimizeOutfile), _generations(generations),
-      _nsga2_cr(nsga2_cr), _nsga2_m(nsga2_m) {
+      _evaluationDuration(evaluationDuration),
+      _optimizeOutfile(optimizeOutfile), _nsga2_individuals(nsga2_individuals),
+      _nsga2_generations(nsga2_generations), _nsga2_cr(nsga2_cr),
+      _nsga2_m(nsga2_m), _samo_is_individuals(samo_is_individuals),
+      _samo_is_maxEvaluations(samo_is_maxEvaluations) {
   int returnCode;
 
   _load = (_period * _loadPercent) / 100;
@@ -284,7 +288,13 @@ Firestarter::Firestarter(
 
     if (_optimizationAlgorithm == "NSGA2") {
       _algorithm = std::make_unique<firestarter::optimizer::algorithm::NSGA2>(
-          _generations, _nsga2_cr, _nsga2_m);
+          _nsga2_generations, _nsga2_cr, _nsga2_m);
+      _individuals = _nsga2_individuals;
+    } else if (_optimizationAlgorithm == "SAMO-IS") {
+      _algorithm = std::make_unique<firestarter::optimizer::algorithm::SAMO_IS>(
+          _samo_is_maxEvaluations, _nsga2_individuals, _nsga2_generations,
+          _nsga2_cr, _nsga2_m);
+      _individuals = _samo_is_individuals;
     } else {
       throw std::invalid_argument("Algorithm " + _optimizationAlgorithm +
                                   " unknown.");
diff --git a/src/firestarter/Main.cpp b/src/firestarter/Main.cpp
index 71174b43..0275e976 100644
--- a/src/firestarter/Main.cpp
+++ b/src/firestarter/Main.cpp
@@ -78,11 +78,13 @@ struct Config {
   std::string optimizationAlgorithm;
   std::vector<std::string> optimizationMetrics;
   std::chrono::seconds evaluationDuration;
-  unsigned individuals;
   std::string optimizeOutfile = "";
-  unsigned generations;
+  unsigned nsga2_individuals;
+  unsigned nsga2_generations;
   double nsga2_cr;
   double nsga2_m;
+  unsigned samo_is_individuals;
+  unsigned samo_is_maxEvaluations;
 
   Config(int argc, const char **argv);
 };
@@ -249,20 +251,24 @@ Config::Config(int argc, const char **argv) {
       cxxopts::value<unsigned>()->default_value("240"), "N");
 
   parser.add_options("optimization")
-    ("optimize", "Run the optimization with one of these algorithms: NSGA2.\nCannot be combined with --measurement.",
+    ("optimize", "Run the optimization with one of these algorithms: NSGA2, SAMO-IS.\nCannot be combined with --measurement.",
       cxxopts::value<std::string>())
     ("optimize-outfile", "Dump the output of the optimization into this\nfile, default: $PWD/$HOSTNAME_$DATE.json",
       cxxopts::value<std::string>())
     ("optimization-metric", "Use a metric for optimization. Metrics listed\nwith cli argument --list-metrics or specified\nwith --metric-from-stdin are valid.",
       cxxopts::value<std::vector<std::string>>())
-    ("individuals", "Number of individuals for the population. For\nNSGA2 specify at least 5 and a multiple of 4,\ndefault: 20",
+    ("nsga2-individuals", "Number of individuals for the population.\nSpecify at least 5 and a multiple of 4,\ndefault: 20",
       cxxopts::value<unsigned>()->default_value("20"))
-    ("generations", "Number of generations, default: 20",
+    ("nsga2-generations", "Number of generations, default: 20",
       cxxopts::value<unsigned>()->default_value("20"))
     ("nsga2-cr", "Crossover probability. Must be in range [0,1[\ndefault: 0.6",
       cxxopts::value<double>()->default_value("0.6"))
     ("nsga2-m", "Mutation probability. Must be in range [0,1]\ndefault: 0.4",
-      cxxopts::value<double>()->default_value("0.4"));
+      cxxopts::value<double>()->default_value("0.4"))
+    ("samo-is-individuals", "Number of individuals for the population.",
+      cxxopts::value<unsigned>()->default_value("20"))
+    ("samo-is-max-evaluations", "Maximum number of evaluations, default: 800",
+      cxxopts::value<unsigned>()->default_value("800"));
 #endif
   // clang-format on
 
@@ -437,16 +443,21 @@ Config::Config(int argc, const char **argv) {
       evaluationDuration = timeout;
       // this will deactivate the watchdog worker
       timeout = std::chrono::seconds::zero();
-      individuals = options["individuals"].as<unsigned>();
       if (options.count("optimize-outfile")) {
         optimizeOutfile = options["optimize-outfile"].as<std::string>();
       }
-      generations = options["generations"].as<unsigned>();
+      nsga2_individuals = options["nsga2-individuals"].as<unsigned>();
+      nsga2_generations = options["nsga2-generations"].as<unsigned>();
       nsga2_cr = options["nsga2-cr"].as<double>();
       nsga2_m = options["nsga2-m"].as<double>();
+      samo_is_individuals = options["samo-is-individuals"].as<unsigned>();
+      samo_is_maxEvaluations =
+          options["samo-is-max-evaluations"].as<unsigned>();
 
-      if (optimizationAlgorithm != "NSGA2") {
-        throw std::invalid_argument("Option --optimize must be any of: NSGA2");
+      if (optimizationAlgorithm != "NSGA2" &&
+          optimizationAlgorithm != "SAMO-IS") {
+        throw std::invalid_argument(
+            "Option --optimize must be any of: NSGA2, SAMO-IS");
       }
     }
 #endif
@@ -481,8 +492,9 @@ int main(int argc, const char **argv) {
         cfg.gpuUseDouble, cfg.listMetrics, cfg.measurement, cfg.startDelta,
         cfg.stopDelta, cfg.measurementInterval, cfg.metricPaths,
         cfg.stdinMetrics, cfg.optimize, cfg.preheat, cfg.optimizationAlgorithm,
-        cfg.optimizationMetrics, cfg.evaluationDuration, cfg.individuals,
-        cfg.optimizeOutfile, cfg.generations, cfg.nsga2_cr, cfg.nsga2_m);
+        cfg.optimizationMetrics, cfg.evaluationDuration, cfg.optimizeOutfile,
+        cfg.nsga2_individuals, cfg.nsga2_generations, cfg.nsga2_cr, cfg.nsga2_m,
+        cfg.samo_is_individuals, cfg.samo_is_maxEvaluations);
 
     firestarter.mainThread();
 
diff --git a/src/firestarter/Optimizer/Algorithm/NSGA2.cpp b/src/firestarter/Optimizer/Algorithm/NSGA2.cpp
index ad010122..5d5ac37b 100644
--- a/src/firestarter/Optimizer/Algorithm/NSGA2.cpp
+++ b/src/firestarter/Optimizer/Algorithm/NSGA2.cpp
@@ -86,7 +86,8 @@ NSGA2::evolve(firestarter::optimizer::Population &pop) {
   {
     std::stringstream ss;
 
-    ss << std::endl << std::setw(7) << "Gen:" << std::setw(15) << "Fevals:";
+    ss << std::endl
+       << std::setw(7) << "NSGA-II Gen:" << std::setw(15) << "Fevals:";
     for (decltype(prob.getNobjs()) i = 0; i < prob.getNobjs(); ++i) {
       ss << std::setw(15) << "ideal" << std::to_string(i + 1u) << ":";
     }
@@ -99,7 +100,8 @@ NSGA2::evolve(firestarter::optimizer::Population &pop) {
       std::vector<double> idealPoint = util::ideal(pop.f());
       std::stringstream ss;
 
-      ss << std::setw(7) << gen << std::setw(15) << prob.getFevals() - fevals0;
+      ss << std::setw(7 + 8) << gen << std::setw(15)
+         << prob.getFevals() - fevals0;
       for (decltype(idealPoint.size()) i = 0; i < idealPoint.size(); ++i) {
         ss << std::setw(15) << idealPoint[i];
       }
diff --git a/src/firestarter/Optimizer/Algorithm/SAMO_IS.cpp b/src/firestarter/Optimizer/Algorithm/SAMO_IS.cpp
new file mode 100644
index 00000000..d4bb9d7a
--- /dev/null
+++ b/src/firestarter/Optimizer/Algorithm/SAMO_IS.cpp
@@ -0,0 +1,148 @@
+#include <firestarter/Logging/Log.hpp>
+#include <firestarter/Optimizer/Algorithm/NSGA2.hpp>
+#include <firestarter/Optimizer/Algorithm/SAMO_IS.hpp>
+#include <firestarter/Optimizer/Individual.hpp>
+#include <firestarter/Optimizer/Problem/SurrogateProblem.hpp>
+#include <firestarter/Optimizer/Util/MultiObjective.hpp>
+
+#include <armadillo>
+
+#include <algorithm>
+#include <set>
+#include <stdexcept>
+
+using namespace firestarter::optimizer::algorithm;
+
+SAMO_IS::SAMO_IS(unsigned maxEvaluations, unsigned nsga2_individuals,
+                 unsigned nsga2_generations, double nsga2_cr, double nsga2_m)
+    : Algorithm(), _maxEvaluations(maxEvaluations),
+      _nsga2_individuals(nsga2_individuals),
+      _nsga2_generations(nsga2_generations), _nsga2_cr(nsga2_cr),
+      _nsga2_m(nsga2_m) {
+  if (nsga2_cr >= 1. || nsga2_cr < 0.) {
+    throw std::invalid_argument("The crossover probability must be in the "
+                                "[0,1[ range, while a value of " +
+                                std::to_string(nsga2_cr) + " was detected");
+  }
+  if (nsga2_m < 0. || nsga2_m > 1.) {
+    throw std::invalid_argument("The mutation probability must be in the [0,1] "
+                                "range, while a value of " +
+                                std::to_string(nsga2_m) + " was detected");
+  }
+}
+
+void SAMO_IS::checkPopulation(firestarter::optimizer::Population const &pop,
+                              std::size_t populationSize) {
+  const auto &prob = pop.problem();
+
+  if (!prob.isMO()) {
+    throw std::invalid_argument("SAMO-IS is a multiobjective algorithms, while "
+                                "number of objectives is " +
+                                std::to_string(prob.getNobjs()));
+  }
+
+  auto minPopulationSize = 2 * (prob.getDims() + 1);
+
+  if (populationSize < minPopulationSize) {
+    throw std::invalid_argument(
+        "for SAMO-IS the population size must be "
+        "greater equal 2 times the number of (variables plus one)."
+        "The population size must be greater equal " +
+        std::to_string(minPopulationSize) + ", while it is " +
+        std::to_string(populationSize));
+  }
+}
+
+firestarter::optimizer::Population
+SAMO_IS::evolve(firestarter::optimizer::Population &pop) {
+  const auto &prob = pop.problem();
+  const auto bounds = prob.getBounds();
+  auto NP = pop.size();
+  auto fevals0 = prob.getFevals();
+
+  this->checkPopulation(
+      const_cast<firestarter::optimizer::Population const &>(pop), NP);
+
+  std::random_device rd;
+  std::mt19937 rng(rd());
+
+  std::pair<Individual, Individual> children;
+
+  for (unsigned gen = 1u; History::x().size() < _maxEvaluations; ++gen) {
+
+    {
+      std::stringstream ss;
+
+      ss << std::endl
+         << std::setw(7) << "SAMO-IS Gen:" << std::setw(15) << "Fevals:";
+      for (decltype(prob.getNobjs()) i = 0; i < prob.getNobjs(); ++i) {
+        ss << std::setw(15) << "ideal" << std::to_string(i + 1u) << ":";
+      }
+      ss << std::endl;
+
+      // Print the logs
+      std::vector<double> idealPoint = util::ideal(pop.f());
+
+      ss << std::setw(7 + 8) << gen << std::setw(15)
+         << prob.getFevals() - fevals0;
+      for (decltype(idealPoint.size()) i = 0; i < idealPoint.size(); ++i) {
+        ss << std::setw(15) << idealPoint[i];
+      }
+
+      firestarter::log::info() << ss.str();
+    }
+
+    // At each generation we make a copy of the population into popnew
+    firestarter::optimizer::Population popnew(pop);
+
+    // run nsga-ii with the surrogate and the limit for 20 generations with an
+    // population size of 100
+    auto surrogateProblem =
+        std::make_shared<firestarter::optimizer::problem::SurrogateProblem>(
+            prob.metrics(), bounds);
+    firestarter::optimizer::Population surrogatePopulation(
+        std::move(surrogateProblem), false);
+
+    surrogatePopulation.generateInitialPopulation(_nsga2_individuals);
+    NSGA2 nsga2(_nsga2_generations, _nsga2_cr, _nsga2_m);
+
+    auto solutions = nsga2.evolve(surrogatePopulation);
+
+    auto fnds_res = util::fast_non_dominated_sorting(solutions.f());
+    auto ndf = std::get<0>(fnds_res);
+
+    std::vector<Individual> C_ND_x(pop.x());
+    std::vector<std::vector<double>> C_ND_f(pop.f());
+
+    // add the non dominated solutions to C_NDset (_x and _f)
+    for (auto const &idx : ndf[0]) {
+      C_ND_x.push_back(solutions.x()[idx]);
+      C_ND_f.push_back(solutions.f()[idx]);
+    }
+
+    {
+
+      auto max_num_solutions = std::min(NP, C_ND_x.size());
+
+      auto best_idx = util::select_best_N_mo(C_ND_f, max_num_solutions);
+
+      // only add new solutions
+      for (auto const &idx : best_idx) {
+        if (!History::find(C_ND_x[idx]).has_value()) {
+          popnew.append(C_ND_x[idx]);
+        }
+      }
+    }
+
+    {
+      auto best_idx = util::select_best_N_mo(popnew.f(), NP);
+      // We insert into the population
+      for (decltype(NP) i = 0;
+           (i < NP) && (History::x().size() < _maxEvaluations); ++i) {
+        pop.insert(i, popnew.x()[best_idx[i]], popnew.f()[best_idx[i]]);
+      }
+    }
+  }
+
+  return pop;
+}
diff --git a/src/firestarter/Optimizer/OptimizerWorker.cpp b/src/firestarter/Optimizer/OptimizerWorker.cpp
index 48819fd5..1a5b8c75 100644
--- a/src/firestarter/Optimizer/OptimizerWorker.cpp
+++ b/src/firestarter/Optimizer/OptimizerWorker.cpp
@@ -62,7 +62,8 @@ void *OptimizerWorker::optimizerThread(void *optimizerWorker) {
   std::this_thread::sleep_for(_this->_preheat);
 
   // For NSGA2 we start with a initial population
-  if (_this->_optimizationAlgorithm == "NSGA2") {
+  if (_this->_optimizationAlgorithm == "NSGA2" ||
+      _this->_optimizationAlgorithm == "SAMO-IS") {
     _this->_population.generateInitialPopulation(_this->_individuals);
   }
 
diff --git a/src/firestarter/Optimizer/Population.cpp b/src/firestarter/Optimizer/Population.cpp
index 9dc1b3a6..981cda24 100644
--- a/src/firestarter/Optimizer/Population.cpp
+++ b/src/firestarter/Optimizer/Population.cpp
@@ -63,19 +63,23 @@ void Population::append(Individual const &ind) {
   std::map<std::string, firestarter::measurement::Summary> metrics;
 
   // check if we already evaluated this individual
-  auto optional_metric = History::find(ind);
-  if (optional_metric.has_value()) {
-    metrics = optional_metric.value();
-  } else {
-    metrics = this->_problem->metrics(ind);
-  }
+  if (_saveHistory) {
+    auto optional_metric = History::find(ind);
+    if (optional_metric.has_value()) {
+      metrics = optional_metric.value();
+    } else {
+      metrics = this->_problem->metrics(ind);
+    }
 
-  auto fitness = this->_problem->fitness(metrics);
+    auto fitness = this->_problem->fitness(metrics);
 
-  this->append(ind, fitness);
+    this->append(ind, fitness);
 
-  if (!optional_metric.has_value()) {
-    History::append(ind, metrics);
+    if (!optional_metric.has_value()) {
+      History::append(ind, metrics);
+    }
+  } else {
+    this->append(ind, this->_problem->fitness(this->_problem->metrics(ind)));
   }
 }
 
@@ -85,7 +89,7 @@ void Population::append(Individual const &ind, std::vector<double> const &fit) {
   for (auto const &v : fit) {
     ss << v << " ";
   }
-  firestarter::log::trace() << ss.str();
+  // firestarter::log::trace() << ss.str();
 
   assert(this->problem().getNobjs() == fit.size());
   assert(this->problem().getDims() == ind.size());
@@ -107,7 +111,8 @@ Individual Population::getRandomIndividual() {
   auto dims = this->problem().getDims();
   auto const bounds = this->problem().getBounds();
 
-  firestarter::log::trace() << "Generating random individual of size: " << dims;
+  // firestarter::log::trace() << "Generating random individual of size: " <<
+  // dims;
 
   Individual out(dims);
 
@@ -117,8 +122,8 @@ Individual Population::getRandomIndividual() {
 
     out[i] = std::uniform_int_distribution<unsigned>(lb, ub)(this->gen);
 
-    firestarter::log::trace()
-        << "  - " << i << ": [" << lb << "," << ub << "]: " << out[i];
+    // firestarter::log::trace()
+    //    << "  - " << i << ": [" << lb << "," << ub << "]: " << out[i];
   }
 
   return out;