diff --git a/R/clarke_wright.R b/R/clarke_wright.R
index 69adf14..f493836 100644
--- a/R/clarke_wright.R
+++ b/R/clarke_wright.R
@@ -77,7 +77,7 @@
 #'    * `site` - The site index (i.e. the index of the (1-indexed) `demand` vector)
 #'              associated to the run.
 #'    * `order` - Integer values providing the visiting order within each run.
-#'    * `load` - The departing load on site `site` in units of `demand` per particular run.
+#'    * `load` - The load after visiting site `site` in units of `demand` per particular run.
 #'  
 #' 
 #'  Unless a site demand exceeds the vehicle capacities it is always assigned
diff --git a/src/Makefile.debug b/src/Makefile.debug
index 968c1a6..164cee3 100644
--- a/src/Makefile.debug
+++ b/src/Makefile.debug
@@ -1,7 +1,7 @@
-debug: clarke_wright.cpp fleet.cpp cpp11.cpp distmat.cpp run.cpp router.cpp tsp_greedy.cpp
+debug: clarke_wright.cpp fleet.cpp cpp11.cpp symmat.cpp run.cpp router.cpp tsp_greedy.cpp
 	g++ -std=gnu++17 -I. \
 		-fstack-protector-strong -fstack-clash-protection -Wformat -fcf-protection \
 		-Wdate-time -Wall -pedantic -Wformat-security -g -O0 -fdiagnostics-color=always \
 		-fsanitize=address,undefined \
-		clarke_wright.cpp fleet.cpp distmat.cpp run.cpp router.cpp tsp_greedy.cpp \
+		clarke_wright.cpp fleet.cpp symmat.cpp run.cpp router.cpp tsp_greedy.cpp \
 		-o debug
diff --git a/src/clarke_wright.cpp b/src/clarke_wright.cpp
index 6d7b585..95a6cb2 100644
--- a/src/clarke_wright.cpp
+++ b/src/clarke_wright.cpp
@@ -1,15 +1,17 @@
 #include <iterator>
 #include <vector>
+#include <functional>
 #include "router.h"
+#include "symmat.h"
 #include <cassert>
 
 // selects demand with a given sign, and always returns positive values
-std::tuple<std::vector<int>, std::vector<double>, distmat<double>>
-select_demand(const std::vector<double> &demand, const distmat<double> &distm, int sign)
+std::tuple<std::vector<int>, std::vector<double>, Distmat>
+select_demand(const std::vector<double> &demand, const Distmat &distm, int sign)
 {
   std::vector<int> new_ind;
   std::vector<double> new_demand;
-  distmat<double> new_distm;
+  Distmat new_distm;
 
   new_ind.reserve(demand.size());
   new_demand.reserve(demand.size());
@@ -37,7 +39,7 @@ select_demand(const std::vector<double> &demand, const distmat<double> &distm, i
     new_distm = distm.sub(new_ind_dist);
   }
 
-  return std::tuple<std::vector<int>, std::vector<double>, distmat<double>>(new_ind, new_demand, new_distm);
+  return std::tuple<std::vector<int>, std::vector<double>, Distmat>(new_ind, new_demand, new_distm);
 }
 
 tbls cpp_clarke_wright(
@@ -60,7 +62,7 @@ tbls cpp_clarke_wright(
   }
 
   auto fleet = std::make_shared<Fleet>(n_res, capacities, restricted_vehicles);
-  distmat<double> distm(distances);
+  Distmat distm(distances);
 
   // we can have positive and negative demands
   // positive = sites are sinks and origin is source
@@ -69,112 +71,55 @@ tbls cpp_clarke_wright(
   // Strategy:
   // Consider and solve positive and negative demands separately (we can use the same solver)
   // Combine those solutions afterwards
+  /*
+    bool have_pos = false;
+    bool have_neg = false;
 
-  bool have_pos = false;
-  bool have_neg = false;
-
-  for (size_t i = 0; i < demand.size(); i++)
-  {
-    if (demand[i] > 0)
+    for (size_t i = 0; i < demand.size(); i++)
     {
-      have_pos = true;
-    }
-    if (demand[i] < 0)
-    {
-      have_neg = true;
-    }
-    if (have_pos && have_neg)
-    {
-      break;
+      if (demand[i] > 0)
+      {
+        have_pos = true;
+      }
+      if (demand[i] < 0)
+      {
+        have_neg = true;
+      }
+      if (have_pos && have_neg)
+      {
+        break;
+      }
     }
-  }
-
+  */
   /* Positive */
+  /*
+    std::vector<int> ind_pos;
+    std::vector<double> demand_pos;
+    distmat<double> distances_pos;
 
-  std::vector<int> ind_pos;
-  std::vector<double> demand_pos;
-  distmat<double> distances_pos;
-
-  if (have_pos)
-  {
-    std::tie(ind_pos, demand_pos, distances_pos) =
-        select_demand(demand, distm, 1);
-  }
-
-  Router runm_pos(
-      demand_pos,
-      std::make_unique<distmat<double>>(distances_pos),
-      fleet);
-
-  if (have_pos)
-  {
-    callback(runm_pos);
-    while (runm_pos.relink_best())
-    {
-      callback(runm_pos);
-    };
-    while (runm_pos.opt_vehicles())
-    {
-    };
-
-    if (!have_neg)
+    if (have_pos)
     {
-      return (runm_pos.runs_as_tbls(demand_pos));
+      std::tie(ind_pos, demand_pos, distances_pos) =
+          select_demand(demand, distm, 1);
     }
-  }
-
-  /* Negative */
-
-  std::vector<int> ind_neg;
-  std::vector<double> demand_neg;
-  distmat<double> distances_neg;
-
-  if (have_neg)
-  {
-    std::tie(ind_neg, demand_neg, distances_neg) =
-        select_demand(demand, distm, -1);
-  }
-
-  Router runm_neg(demand_neg,
-                      std::make_unique<distmat<double>>(distances_neg),
-                      fleet);
-
-  if (have_neg)
-  {
-    callback(runm_neg);
-    while (runm_neg.relink_best())
-    {
-      callback(runm_neg);
-    };
-
-    while (runm_neg.opt_vehicles())
-    {
-    };
-
-    if (!have_pos)
-    {
-      for (auto &dmnd : demand_neg)
-      {
-        dmnd = -dmnd;
-      }
-      return (runm_neg.runs_as_tbls(demand_neg));
-    }
-  }
+  */
 
   // if we have both, combine then and optimize again
-  Router runm_all(runm_pos, runm_neg, distm, ind_pos, ind_neg);
+  Router router(
+    std::make_shared<std::vector<double>>(demand),
+    std::make_unique<Distmat>(distm), fleet);
 
-  callback(runm_all);
-  while (runm_all.relink_best([](double l1, double l2)
-                              { return (std::max(l1, l2)); }))
+  callback(router);
+  while (router.relink_best())
   {
-    callback(runm_all);
+    callback(router);
   };
-  while (runm_all.opt_vehicles())
+  while (router.optimize_vehicles())
   {
   };
+  // router.optimize_runs_order();
 
-  return (runm_all.runs_as_tbls(demand));
+  return (router.runs_as_tbls());
 }
 
 #ifndef NDEBUG
@@ -183,10 +128,10 @@ int main()
 {
   tbls cols =
       cpp_clarke_wright(
-          std::vector<double>{-3, -2, 1, 2},
-          std::vector<double>{5, 10, 5, 4, 5, 2, 12, 13, 10, 3},
-          std::vector<int>{100},
-          std::vector<double>{5},
+          std::vector<double>{14.14, 14.37, 7.86},
+          std::vector<double>{8.128214, 7.837354, 3.162988, 12.616152, 6.427298, 10.7306495},
+          std::vector<int>{2, 100},
+          std::vector<double>{33, 44},
           std::vector<int>{},
           std::vector<int>{});
 
@@ -203,4 +148,16 @@ int main()
 
   return 0;
 }
-#endif
\ No newline at end of file
+#endif
+
+
+/*
+> demand
+[1] 14.148060 14.370754  7.861395
+> pos
+     pos_x     pos_y
+1 0.000000  0.000000
+2 6.608953  4.731766
+3 2.834910 -7.306668
+4 0.381919  3.139846
+*/
\ No newline at end of file
diff --git a/src/fleet.cpp b/src/fleet.cpp
index fe51f89..2d42e15 100644
--- a/src/fleet.cpp
+++ b/src/fleet.cpp
@@ -9,7 +9,7 @@ int Fleet::capacity(const int vehicle) const
     return vehicle_caps[vehicle];
 }
 
-void Fleet::release_vehicle(const int vehicle)
+void Fleet::release_vehicle(const VehicleTypeID vehicle)
 {
     if (vehicle_avail[vehicle] < INT_MAX)
     {
@@ -17,13 +17,8 @@ void Fleet::release_vehicle(const int vehicle)
     }
 }
 
-void Fleet::reserve_vehicle(const int vehicle)
+void Fleet::reserve_vehicle(const VehicleTypeID vehicle)
 {
-    if (vehicle == -1)
-    {
-        throw std::runtime_error("blah.");
-    }
-
     if (vehicle_avail[vehicle] <= 0)
     {
         throw std::runtime_error("No available vehicles of this type to pop.");
@@ -36,22 +31,21 @@ void Fleet::reserve_vehicle(const int vehicle)
 }
 
 template <typename ForwardIt>
-bool Fleet::is_vehicle_restricted(const int vehicle, ForwardIt sites) const
+bool Fleet::is_vehicle_restricted(const VehicleTypeID vehicle, ForwardIt sites) const
 {
-    bool restricted = false;
-
     for (const auto site : sites)
     {
-        restricted = restricted ||
-                     (this->restricted_vehicles[site].find(vehicle) !=
-                      this->restricted_vehicles[site].end());
+        if (this->restricted_vehicles[site].find(vehicle) != this->restricted_vehicles[site].end())
+        {
+            return true;
+        }
     }
 
-    return restricted;
+    return false;
 }
 
 template <typename ForwardIt>
-int Fleet::find_fitting_vehicle(ForwardIt sites, const double max_load, bool max_load_soft) const
+std::optional<VehicleTypeID> Fleet::find_fitting_vehicle(ForwardIt sites, const double max_load, bool max_load_soft) const
 {
     for (size_t vehicle = 0; vehicle < vehicle_caps.size(); vehicle++)
     {
@@ -59,7 +53,7 @@ int Fleet::find_fitting_vehicle(ForwardIt sites, const double max_load, bool max
             max_load <= vehicle_caps[vehicle] &&
             !is_vehicle_restricted(vehicle, sites))
         {
-            return vehicle;
+            return (VehicleTypeID)vehicle;
         }
     }
 
@@ -71,13 +65,13 @@ int Fleet::find_fitting_vehicle(ForwardIt sites, const double max_load, bool max
             if (vehicle_avail[vehicle] >= 1 &&
                 !is_vehicle_restricted(vehicle, sites))
             {
-                return vehicle;
+                return (VehicleTypeID)vehicle;
             }
         }
     }
 
-    return -1;
+    return std::nullopt;
 }
 
-template int Fleet::find_fitting_vehicle(union_view<int, std::unordered_set>, const double, bool) const;
-template int Fleet::find_fitting_vehicle(std::unordered_set<int>, const double, bool) const;
+template std::optional<VehicleTypeID> Fleet::find_fitting_vehicle(union_view<int, std::list>, const double, bool) const;
+template std::optional<VehicleTypeID> Fleet::find_fitting_vehicle(std::list<int>, const double, bool) const;
diff --git a/src/fleet.h b/src/fleet.h
index 36a5185..3327681 100644
--- a/src/fleet.h
+++ b/src/fleet.h
@@ -2,7 +2,11 @@
 #define FLEET_H
 
 #include <vector>
+#include <list>
 #include <unordered_set>
+#include <optional>
+
+typedef int VehicleTypeID;
 
 class Fleet
 {
@@ -12,14 +16,16 @@ class Fleet
       : vehicle_avail(vehicle_avail),
         vehicle_caps(vehicle_caps),
         restricted_vehicles(restricted_vehicles) {};
-  void reserve_vehicle(const int vehicle);
-  void release_vehicle(const int vehicle);
+  void reserve_vehicle(const VehicleTypeID vehicle);
+  void release_vehicle(const VehicleTypeID vehicle);
 
   template <typename ForwardIt>
-  int find_fitting_vehicle(ForwardIt sites, const double max_load,
-                           bool max_load_soft) const;
+  std::optional<VehicleTypeID> find_fitting_vehicle(
+      ForwardIt sites,
+      const double max_load,
+      bool max_load_soft) const;
 
-  int capacity(const int vehicle) const;
+  int capacity(const VehicleTypeID vehicle) const;
 
 private:
   // Number and capacity of particular vehicle types
@@ -30,7 +36,7 @@ class Fleet
   std::vector<std::unordered_set<int>> restricted_vehicles;
 
   template <typename ForwardIt>
-  bool is_vehicle_restricted(int vehicle, ForwardIt sites) const;
+  bool is_vehicle_restricted(const VehicleTypeID vehicle, ForwardIt sites) const;
 };
 
 #endif
\ No newline at end of file
diff --git a/src/r_clarke_wright.cpp b/src/r_clarke_wright.cpp
index 7a3ea52..8a03c0c 100644
--- a/src/r_clarke_wright.cpp
+++ b/src/r_clarke_wright.cpp
@@ -1,7 +1,8 @@
 #include <cpp11.hpp>
 #include <iterator>
 #include <vector>
-#include "runmanager.h"
+#include "router.h"
+#include <functional>
 
 using namespace cpp11;
 
@@ -38,7 +39,7 @@ tbls cpp_clarke_wright(
     const std::vector<double> &capacities,
     const std::vector<int> &restr_sites,
     const std::vector<int> &restr_vehicles,
-    std::function<void(RunManager &)> callback = [](RunManager &) {});
+    std::function<void(Router &)> callback = [](Router &) {});
 
 [[cpp11::register]]
 cpp11::writable::list r_cpp_clarke_wright(
@@ -77,8 +78,8 @@ list cpp_clarke_wright_stepwise(
       capacities,
       restr_sites,
       restr_vehicles,
-      [&steps, &demand](RunManager &runm)
-      { steps.push_back(tbls_to_dfs(runm.runs_as_tbls(demand))); });
+      [&steps, &demand](Router &router)
+      { steps.push_back(tbls_to_dfs(router.runs_as_tbls())); });
 
   return steps;
 }
diff --git a/src/router.cpp b/src/router.cpp
index 4d093dd..a68df84 100644
--- a/src/router.cpp
+++ b/src/router.cpp
@@ -8,222 +8,112 @@
 #include <vector>
 #include <memory>
 #include <cassert>
+#include <optional>
 
-Router::Router(const std::vector<double> &demand,
-               std::unique_ptr<distmat<double>> distances,
-               std::shared_ptr<Fleet> fleet)
-    : fleet(fleet),
-      distances(std::move(distances)),
-      sites_relinked(demand.size(), 0),
-      runs(demand.size())
+// creates potentially multiple singleton runs but returns only a single "run".
+run Router::create_initial_runs(Site s, double demand, std::shared_ptr<Fleet> fleet,
+                                std::shared_ptr<Distmat> distances)
 {
-  this->savings = calc_savings(*(this->distances));
-
-  runs.reserve(demand.size()); // just in case
-  for (size_t i = 0; i < demand.size(); i++)
+  VehicleTypeID vehicle;
+  try
   {
-    if (demand[i] <= 0)
-    {
-      throw std::runtime_error("Router: demand has to be strictly positive");
-    }
-
-    runs[i] = std::make_shared<run>(i, demand[i]);
+    vehicle = fleet->find_fitting_vehicle(std::list<Site>{s}, demand, true).value();
   }
-
-  fixed_singleton_runs = std::vector<run>();
-
-  // first vehicle assignments (iterate over runs)
-  // initial runs have only a single site
-  for (auto &run : runs)
+  catch (const std::bad_optional_access &e)
   {
-    int vehicle = fleet->find_fitting_vehicle(run->sites(),
-                                              run->max_load,
-                                              true);
-    fleet->reserve_vehicle(vehicle);
-
-    // special treatment for the case when demand is higher than capacity
-    while (run->max_load > fleet->capacity(vehicle))
-    {
-      run->max_load -= fleet->capacity(vehicle);
-      this->fixed_singleton_runs.emplace_back(
-          *(run->sites().begin()), fleet->capacity(vehicle), vehicle);
-
-      vehicle = fleet->find_fitting_vehicle(run->sites(),
-                                            run->max_load,
-                                            true);
-
-      if (vehicle == -1)
-      {
-        throw std::runtime_error(
-            "Not enough vehicles available to fulfill all demands trivially."
-            " Solver cannot proceed in that case.");
-      }
-
-      fleet->reserve_vehicle(vehicle);
-    }
-
-    // only add the last one to the state
-    run->vehicle = vehicle;
+    throw std::runtime_error(
+        "Not enough vehicles available to fulfill all demands trivially."
+        " Solver cannot proceed in that case.");
   }
-}
+  fleet->reserve_vehicle(vehicle);
+  int capacity = fleet->capacity(vehicle);
 
-int unique_count(union_view<int, std::vector> uv)
-{
-  std::unordered_set<int> s;
-  for (auto it = uv.begin(); it != uv.end(); ++it)
+  if (demand > capacity)
   {
-    s.insert(*it);
+    this->fixed_singleton_runs.emplace_back(s, capacity, vehicle, distances);
+    return (create_initial_runs(s, demand - capacity, fleet, distances));
   }
-  return s.size();
-}
-
-Router::Router(const Router &runm1, const Router &runm2,
-               const distmat<double> &new_distances,
-               const std::vector<int> &site_ind_map1,
-               const std::vector<int> &site_ind_map2)
-    : fleet(runm1.fleet),
-      distances(std::make_unique<distmat<double>>(new_distances))
-{
-  if (runm1.fleet != runm2.fleet)
-  {
-    throw std::runtime_error("Router: cannot combine two Routers with different fleets");
-  }
-
-  // create inverse maps
-  // std::map<int, int> inv_site_ind_map1;
-  // for (size_t i = 0; i < site_ind_map1.size(); i++) {
-  //   inv_site_ind_map1[site_ind_map1[i]] = i;
-  // }
-  // std::map<int, int> inv_site_ind_map2;
-  // for (size_t i = 0; i < site_ind_map2.size(); i++) {
-  //   inv_site_ind_map2[site_ind_map2[i]] = i;
-  // }
-
-  // we have to mainly take care of identifying sites of run1 and run2 correctly
-  size_t site_size = unique_count(union_view(site_ind_map1, site_ind_map2));
-
-  this->sites_relinked = std::vector<int>(site_size);
-  for (size_t i = 0; i < runm1.sites_relinked.size(); i++)
-  {
-    this->sites_relinked[site_ind_map1[i]] = runm1.sites_relinked[i];
-  }
-  for (size_t i = 0; i < runm2.sites_relinked.size(); i++)
+  else
   {
-    this->sites_relinked[site_ind_map2[i]] = runm2.sites_relinked[i];
+    return run(s, demand, vehicle, distances);
   }
+}
 
-  this->runs = std::vector<std::shared_ptr<run>>(site_size);
-  for (const auto &rptr : runm1.runs)
-  {
-    std::unordered_set<int> new_sites;
-    for (auto old_site : rptr->sites())
-    {
-      new_sites.insert(site_ind_map1[old_site]);
-    }
-    auto new_run = std::make_shared<run>(new_sites, rptr->max_load, rptr->vehicle);
-    for (const auto site : new_run->sites())
-    {
-      this->runs[site] = new_run;
-    }
-  }
-  for (const auto &rptr : runm2.runs)
-  {
-    std::unordered_set<int> new_sites;
-    for (auto old_site : rptr->sites())
-    {
-      new_sites.insert(site_ind_map2[old_site]);
-    }
-    auto new_run = std::make_shared<run>(new_sites, rptr->max_load, rptr->vehicle);
-    for (const auto site : new_run->sites())
-    {
-      this->runs[site] = new_run;
-    }
-  }
+Router::Router(const std::shared_ptr<std::vector<double>> demand,
+               const std::unique_ptr<Distmat> distances,
+               std::shared_ptr<Fleet> fleet)
+    : fleet(fleet),
+      distances(std::make_shared<Distmat>(std::move(*distances))),
+      demand(demand),
+      sites_start(demand->size(), true),
+      sites_end(demand->size(), true),
+      runs(demand->size())
+{
+  this->savings = calc_savings(*(this->distances));
 
-  this->fixed_singleton_runs = std::vector<run>();
-  this->fixed_singleton_runs.reserve(runm1.fixed_singleton_runs.size() + runm2.fixed_singleton_runs.size());
-  for (const auto &srun : runm1.fixed_singleton_runs)
-  {
-    std::unordered_set<int> new_sites;
-    for (auto site : srun.sites())
-    {
-      new_sites.insert(site_ind_map1[site]);
-    }
-    this->fixed_singleton_runs.emplace_back(new_sites, srun.max_load, srun.vehicle);
-  }
-  for (const auto &srun : runm2.fixed_singleton_runs)
-  {
-    std::unordered_set<int> new_sites;
-    for (auto site : srun.sites())
-    {
-      new_sites.insert(site_ind_map2[site]);
-    }
-    this->fixed_singleton_runs.emplace_back(new_sites, srun.max_load, srun.vehicle);
-  }
+  fixed_singleton_runs = std::vector<run>();
 
-  // we only want to consider connecting runs from one set to the other
-  this->consider_optim1 = std::vector<int>(site_size, false);
-  this->consider_optim2 = std::vector<int>(site_size, false);
-  for (const auto site : site_ind_map1)
+  runs.reserve(demand->size()); // just in case
+  for (size_t i = 0; i < demand->size(); i++)
   {
-    this->consider_optim1[site] = true;
+    runs[i] = std::make_shared<run>(
+        create_initial_runs(i, (*demand)[i], fleet, this->distances));
   }
-  for (const auto site : site_ind_map2)
-  {
-    this->consider_optim2[site] = true;
-  }
-
-  // recalculate savings (there are more efficient ways)
-  this->savings = calc_savings(*(this->distances));
 }
 
-void Router::combine_runs(const int a, const int b, const int new_vehicle, binop_dbl combine_load)
+void Router::combine_runs(const Site a, const Site b, const VehicleTypeID new_vehicle)
 {
   assert(a != b);
   assert(runs[a] != runs[b]);
-  assert(links_to_origin(a) && links_to_origin(b));
+  assert(end_of_run(a));
+  assert(start_of_run(b));
 
   /*
   In the original algorithm we are only supposed to attempt to combine vertices that
   are connected to the origin.
   We used to keep track of the whole graph to determine that, but that was
   not necessary if the initial state consists of singleton runs all connected
-  to the origin, i.e. the form a cycle (ORIGIN - v - ORIGIN).
+  to the origin, i.e. the form a run_ptrle (ORIGIN - v - ORIGIN).
   When combining two vertices v and b for the first time, v is adjacent to ORIGIN
   and b. When combined a second time (to another vertex c), v is adjacent to b and c.
   -> A vertex can be combined at most two times, then it will not be connected to
      the origin anymore.
   */
 
-  sites_relinked[a] += 1;
-  sites_relinked[b] += 1;
-  runs[a]->combine(*runs[b], new_vehicle, combine_load);
+  runs[a]->combine(*runs[b], new_vehicle);
+
+  sites_end[a] = false;
+  sites_start[b] = false;
 
   // all vertices in the runs are affected,
   // we need to reset the pointer of the ones that b pointed
-  // to to point to the same cycle
+  // to to point to the same run_ptrle
   for (const auto site : runs[a]->sites())
   {
     runs[site] = runs[a];
   }
 }
 
-bool Router::links_to_origin(const int a) const
+bool Router::end_of_run(const Site a) const
+{
+  return sites_end[a];
+}
+
+bool Router::start_of_run(const Site a) const
 {
-  return (sites_relinked[a] < 2);
+  return sites_start[a];
 }
 
-bool Router::edges_share_run(const int a, const int b) const
+bool Router::sites_share_run(const Site a, const Site b) const
 {
   return (runs[a] == runs[b]);
 }
 
-// we create a symmat that is one size smaller than the distances
+// we create a Distmat that is one size smaller than the distances
 // (only calculate for sites)
-distmat<double> Router::calc_savings(const distmat<double> &d) const
+Distmat Router::calc_savings(const Distmat &d) const
 {
-  distmat<double> savings(d.size() - 1, 0);
-
+  Distmat savings(d.size() - 1, 0);
   for (int i = 1; i < savings.size(); i++)
   {
     for (int j = 0; j < i; j++)
@@ -235,47 +125,79 @@ distmat<double> Router::calc_savings(const distmat<double> &d) const
   return savings;
 }
 
+std::optional<std::tuple<Site, Site>> Router::run_merge_order(const Site i, const Site j) const
+{
+  // if we have only attach a singleton run, we attach the new one at the end
+  // if its associated to negative demand, and in the beginning, if it's positive
+  bool i_is_singleton = end_of_run(i) && start_of_run(i);
+  bool j_is_singleton = end_of_run(j) && start_of_run(j);
+  if (
+      (i_is_singleton && ((*demand)[i] > 0) && start_of_run(j)) ||
+      (j_is_singleton && ((*demand)[j] > 0) && end_of_run(i)) ||
+      (end_of_run(i) && start_of_run(j)))
+  {
+    return std::make_tuple(i, j);
+  }
+
+  if (
+      (i_is_singleton && ((*demand)[i] < 0) && end_of_run(j)) ||
+      (j_is_singleton && ((*demand)[j] < 0) && start_of_run(i)) ||
+      (end_of_run(j) && start_of_run(i)))
+  {
+    return std::make_tuple(j, i);
+  }
+
+  return std::nullopt;
+}
+
 // returns Site 1, Site 2, Used vehicle
-std::tuple<int, int, int> Router::best_link(binop_dbl combine_load) const
+std::optional<std::tuple<Site, Site, VehicleTypeID>> Router::best_link() const
 {
-  std::tuple<int, int, int> best_link = {-1, -1, -1};
+  std::optional<std::tuple<Site, Site, VehicleTypeID>> best_link(std::nullopt);
   double max_val = 0;
 
-  for (int i = 1; i < savings.size(); i++)
+  for (Site i = 1; i < savings.size(); i++)
   {
-    for (int j = 0; j < i; j++)
+    for (Site j = 0; j < i; j++)
     {
       // printf("---\n");
       // printf("Link (%d,%d)\n", i, j);
       // printf("orig1 %d\n", runm.links_to_origin(i));
       // printf("orig2 %d\n", runm.links_to_origin(j));
       // printf("selected vehicle %d\n", select_vehicle(vehicle_avail, vehicle_caps, site_vehicle, load, restricted_vehicles, runm, i, j));
-      // printf("share cycle %d\n", runm.edges_share_cycle(i, j));
+      // printf("share run_ptrle %d\n", runm.sites_share_run_ptrle(i, j));
 
-      int selected_vehicle;
+      std::optional<VehicleTypeID> selected_vehicle;
       double saving;
 
-      if (is_considered(i, j) &&
-          !edges_share_run(i, j) &&
+      if (!sites_share_run(i, j) &&
           ((saving = savings.get(i, j)) > max_val) &&
-          links_to_origin(i) && links_to_origin(j))
+          ((end_of_run(i) && start_of_run(j)) || (end_of_run(j) && start_of_run(i))))
       {
-        fleet->release_vehicle(runs[i]->vehicle);
-        fleet->release_vehicle(runs[j]->vehicle);
+        // we need to decide how to merge runs together
+        std::optional<std::tuple<Site, Site>> mo = run_merge_order(i, j);
+
+        if (mo)
+        {
+          auto [end, start] = mo.value();
 
-        selected_vehicle =
-            fleet->find_fitting_vehicle(
-                union_view(runs[i]->sites(), runs[j]->sites()),
-                combine_load(runs[i]->max_load, runs[j]->max_load),
-                false);
+          fleet->release_vehicle(runs[end]->vehicle());
+          fleet->release_vehicle(runs[start]->vehicle());
 
-        fleet->reserve_vehicle(runs[i]->vehicle);
-        fleet->reserve_vehicle(runs[j]->vehicle);
+          selected_vehicle =
+              fleet->find_fitting_vehicle(
+                  union_view(runs[end]->sites(), runs[start]->sites()),
+                  runs[end]->combined_max_load(*runs[start]),
+                  false);
 
-        if (selected_vehicle != -1)
-        {
-          max_val = saving;
-          best_link = {i, j, selected_vehicle};
+          fleet->reserve_vehicle(runs[end]->vehicle());
+          fleet->reserve_vehicle(runs[start]->vehicle());
+
+          if (selected_vehicle)
+          {
+            max_val = saving;
+            best_link = {end, start, selected_vehicle.value()};
+          }
         }
       }
     }
@@ -286,51 +208,30 @@ std::tuple<int, int, int> Router::best_link(binop_dbl combine_load) const
 
 // TRUE if something got relinked,
 // FALSE if nothing got relinked (i.e. the procedure stabilized)
-bool Router::relink_best(binop_dbl combine_load)
+bool Router::relink_best()
 {
-  int a;
-  int b;
-  int vehicle;
-  std::tie(a, b, vehicle) = best_link(combine_load);
-
-  // printf("---\n");
-  // printf("Best Link (%d,%d)\n", a, b);
-  // printf("orig1 %d\n", runm.links_to_origin(a));
-  // printf("orig2 %d\n", runm.links_to_origin(b));
-  // printf("selected vehicle %d\n", vehicle);
-  // printf("share cycle %d\n", runm.edges_share_cycle(a, b));
-
-  if (!((a == b) && (a == -1)))
+  std::optional<std::tuple<Site, Site, VehicleTypeID>> best_link = this->best_link();
+  if (!best_link)
   {
-    // return two vehicles
-    fleet->release_vehicle(this->runs[a]->vehicle);
-    fleet->release_vehicle(this->runs[b]->vehicle);
-    fleet->reserve_vehicle(vehicle);
+    return false;
+  }
 
-    combine_runs(a, b, vehicle, combine_load);
+  Site a;
+  Site b;
+  VehicleTypeID vehicle;
+  std::tie(a, b, vehicle) = best_link.value();
 
-    // we use consider_optim1 only when we combine positive and negative demand runs
-    // so we can simply check the existance of this variable
-    if (consider_optim1.size() > 0)
-    {
-      // we want to combine pos and negative tours only once, so we remove the relevant sites here.
-      for (const int site : runs[a]->sites())
-      {
-        // TODO: this is an abuse of sites_relinked, where we use the fact that elements with sites_relinked = 2
-        // are not considered.
-        sites_relinked[site] = 2;
-      }
-    }
+  // return two vehicles
+  fleet->release_vehicle(this->runs[a]->vehicle());
+  fleet->release_vehicle(this->runs[b]->vehicle());
+  fleet->reserve_vehicle(vehicle);
 
-    return true;
-  }
-  else
-  {
-    return false;
-  }
+  combine_runs(a, b, vehicle);
+
+  return true;
 }
 
-bool Router::opt_vehicles()
+bool Router::optimize_vehicles()
 {
   /*
   It might seem more efficient to release all vehicles first and then determine them from
@@ -343,24 +244,10 @@ bool Router::opt_vehicles()
   bool changed = false;
 
   // then reassign fitting vehicles
-  for (auto &run : runs)
+  for (auto run : runs)
   {
-    int old_vehicle = run->vehicle;
-
-    // we are guaranteed to find at least that vehicle one again
-    fleet->release_vehicle(old_vehicle);
-
-    int vehicle =
-        fleet->find_fitting_vehicle(
-            run->sites(),
-            run->max_load,
-            false);
-
-    fleet->reserve_vehicle(vehicle);
-
-    if (old_vehicle != vehicle)
+    if (run->reassign_vehicle(*fleet))
     {
-      run->vehicle = vehicle;
       changed = true;
     }
   }
@@ -368,34 +255,22 @@ bool Router::opt_vehicles()
   return changed;
 }
 
-bool Router::is_considered(const int site1, const int site2) const
+void Router::optimize_runs_order()
 {
-  return (consider_optim1.size() == 0 ||
-          (consider_optim1[site1] && consider_optim2[site2]) || (consider_optim1[site2] && consider_optim2[site1]));
-}
-
-double run_distance(const std::vector<int> ordered_sites,
-                    const distmat<double> &d)
-{
-  auto it = ordered_sites.begin();
-  double distance = d.get(0, 1 + *it);
-  for (; it < (ordered_sites.end() - 1); it++)
+  for (auto run : runs)
   {
-    distance += d.get(1 + *it, 1 + *(it + 1));
+    run->optimize_route_order();
   }
-  distance += d.get(0, 1 + *it);
-
-  return distance;
 }
 
-tbls Router::runs_as_tbls(const std::vector<double> &demand) const
+tbls Router::runs_as_tbls() const
 {
   typedef std::shared_ptr<run> T;
 
   size_t col_size = runs.size() + fixed_singleton_runs.size();
 
   std::map<T, int> visited_runs;
-  std::map<int, std::vector<int>> orders;
+  std::map<int, std::list<Site>> orders;
   std::map<int, double> run_dists;
 
   tbl_run_site run_site_cols = {
@@ -410,26 +285,26 @@ tbls Router::runs_as_tbls(const std::vector<double> &demand) const
   size_t i = 0;
   for (; i < runs.size(); i++)
   {
-    std::vector<int> order;
+    std::list<Site> order;
     double run_dist;
-    T cyc = runs[i];
+    std::shared_ptr<run> run_ptr = runs[i];
 
     std::get<1>(run_site_cols)[i] = i;
 
-    // check if we have seen cyc before
-    if (visited_runs.count(cyc) > 0)
+    // check if we have seen run_ptr before
+    if (visited_runs.count(run_ptr) > 0)
     {
-      order = orders[visited_runs[cyc]];
-      run_dist = run_dists[visited_runs[cyc]];
+      order = orders[visited_runs[run_ptr]];
+      run_dist = run_dists[visited_runs[run_ptr]];
 
-      std::get<0>(run_site_cols)[i] = visited_runs[cyc];
+      std::get<0>(run_site_cols)[i] = visited_runs[run_ptr];
     }
     else // if we did not see it before
     {
-      visited_runs.insert({cyc, run_id});
+      visited_runs.insert({run_ptr, run_id});
       // we reorder each run again (by solving the TSP)
-      order = cyc->ordered_sites(*(this->distances), consider_optim1, consider_optim2);
-      run_dist = run_distance(order, *(this->distances));
+      order = run_ptr->sites();
+      run_dist = run_ptr->distance();
 
       orders.insert({run_id, order});
       run_dists.insert({run_id, run_dist});
@@ -454,29 +329,13 @@ tbls Router::runs_as_tbls(const std::vector<double> &demand) const
   for (const auto &[run, run_id] : visited_runs)
   {
     std::get<0>(run_cols)[run_id] = run_id;
-    std::get<1>(run_cols)[run_id] = run->vehicle;
-    std::get<2>(run_cols)[run_id] = run->max_load;
+    std::get<1>(run_cols)[run_id] = run->vehicle();
+    std::get<2>(run_cols)[run_id] = run->max_load();
     std::get<3>(run_cols)[run_id] = run_dists[run_id];
 
-    // we also perform the load calculation here per run
-    double sdemand = 0;
-    double min_sdemand = 0;
-    for (auto site : orders[run_id])
-    {
-      sdemand -= demand[site];
-      // first write in the demand wrong by an offset
-      std::get<3>(run_site_cols)[site] = sdemand;
-
-      if (sdemand < min_sdemand)
-      {
-        min_sdemand = sdemand;
-      }
-    }
-
-    // now update everything by that offset
-    for (auto site : orders[run_id])
+    for (const auto [site, load] : run->load_after_visit(*demand))
     {
-      std::get<3>(run_site_cols)[site] -= min_sdemand;
+      std::get<3>(run_site_cols)[site] = load;
     }
   }
 
@@ -487,11 +346,11 @@ tbls Router::runs_as_tbls(const std::vector<double> &demand) const
     std::get<0>(run_site_cols)[i] = run_id;
     std::get<1>(run_site_cols)[i] = site;
     std::get<2>(run_site_cols)[i] = 0;
-    std::get<3>(run_site_cols)[i] = run.max_load;
+    std::get<3>(run_site_cols)[i] = run.max_load();
 
     std::get<0>(run_cols)[i] = run_id;
-    std::get<1>(run_cols)[i] = run.vehicle;
-    std::get<2>(run_cols)[i] = run.max_load;
+    std::get<1>(run_cols)[i] = run.vehicle();
+    std::get<2>(run_cols)[i] = run.max_load();
     std::get<3>(run_cols)[i] = 2 * distances->get(0, 1 + site);
 
     run_id++;
@@ -500,7 +359,7 @@ tbls Router::runs_as_tbls(const std::vector<double> &demand) const
 
   tbl_site site_cols = {
       std::vector<int>(runs.size()),
-      std::vector<double>(runs.size()) = demand};
+      std::vector<double>(runs.size()) = *demand};
 
   for (i = 0; i < runs.size(); i++)
   {
diff --git a/src/router.h b/src/router.h
index b233be9..5e4b001 100644
--- a/src/router.h
+++ b/src/router.h
@@ -7,6 +7,8 @@
 
 #include "run.h"
 #include "fleet.h"
+#include "site.h"
+#include "symmat.h"
 
 using tbl_run =
     std::tuple<              // Run
@@ -39,63 +41,59 @@ using tbls =
 class Router
 {
 public:
-  // creates one singleton runs for each site with the given demand and already assigns
-  // vehicles from a fleet
-  Router(const std::vector<double> &demand,
-         std::unique_ptr<distmat<double>> distances,
-         std::shared_ptr<Fleet> fleet);
-
-  // creates a new Router by combining two existing ones
-  // note that their fleets have to be identical for that to make sense
-  Router(const Router &runm1, const Router &runm2,
-         const distmat<double> &new_distances,
-         const std::vector<int> &site_ind_map1,
-         const std::vector<int> &site_ind_map2);
-
-  bool relink_best(binop_dbl combine_load = [](double l1, double l2)
-                   { return (l1 + l2); });
-
-  // After we have the final routes, we might still be able to assign
-  // better vehicles for each route
-  // (we might have released some high-priority vehicles on the way which
-  // are now unused)
-  bool opt_vehicles();
-
-  // returns the current runs as column vectors for
-  // data frame creation
-  tbls runs_as_tbls(const std::vector<double> &demand) const;
-
-  std::shared_ptr<Fleet> fleet;
-  const std::unique_ptr<distmat<double>> distances;
+    // creates one singleton runs for each site with the given demand and already assigns
+    // vehicles from a fleet
+    Router(const std::shared_ptr<std::vector<double>> demand,
+           const std::unique_ptr<Distmat> distances,
+           std::shared_ptr<Fleet> fleet);
+
+    run create_initial_runs(Site s, double demand, std::shared_ptr<Fleet> fleet,
+                            const std::shared_ptr<Distmat> distances);
+    // creates a new Router by combining two existing ones
+    // note that their fleets have to be identical for that to make sense
+
+    bool relink_best();
+
+    // After we have the final routes, we might still be able to assign
+    // better vehicles for each route
+    // (we might have released some high-priority vehicles on the way which
+    // are now unused)
+    bool optimize_vehicles();
+    void optimize_runs_order();
+
+    // returns the current runs as column vectors for
+    // data frame creation
+    tbls runs_as_tbls() const;
 
 private:
-  // combines the two runs traversing site a and site b with the new vehicle new_vehicle.
-  void combine_runs(const int a, const int b, const int new_vehicle, binop_dbl combine_load);
+    // combines the two runs traversing site a and site b with the new vehicle new_vehicle.
+    void combine_runs(const Site a, const Site b, const VehicleTypeID new_vehicle);
 
-  // is site a directly linked to the origin via its traversing run?
-  bool links_to_origin(const int a) const;
+    // is site a directly linked to the origin via its traversing run?
+    bool end_of_run(const Site a) const;
+    bool start_of_run(const Site a) const;
+    std::optional<std::tuple<Site, Site>> run_merge_order(const Site a, const Site b) const;
 
-  // are the sites a and site b traversed by the same run?
-  bool edges_share_run(const int a, const int b) const;
+    // are the sites a and site b traversed by the same run?
+    bool sites_share_run(const Site a, const Site b) const;
 
-  distmat<double> calc_savings(const distmat<double> &d) const;
+    Distmat calc_savings(const Distmat &d) const;
 
-  std::tuple<int, int, int> best_link(binop_dbl combine_load) const;
+    std::optional<std::tuple<Site, Site, VehicleTypeID>> best_link() const;
 
-  distmat<double> savings;
-  std::vector<int> sites_relinked;
-  std::vector<run> fixed_singleton_runs; // those runs are not dynamic, i.e. they won't be changed
+    std::shared_ptr<Fleet> fleet;
+    const std::shared_ptr<Distmat> distances;
+    const std::shared_ptr<std::vector<double>> demand;
 
-  // a vector of runs (of length of the sites): each site has a reference to
-  // the runs it belongs to (which in turn has all the other references)
-  std::vector<std::shared_ptr<run>> runs;
+    Distmat savings;
+    // std::vector<int> sites_relinked;
+    std::vector<bool> sites_start;
+    std::vector<bool> sites_end;
+    std::vector<run> fixed_singleton_runs; // those runs are not dynamic, i.e. they won't be changed
 
-  bool is_considered(const int site1, const int site2) const;
-
-  // site-indexed: consider only site1-site2 combination for optimization
-  // if vectors are not empty
-  std::vector<int> consider_optim1;
-  std::vector<int> consider_optim2;
+    // a vector of runs (of length of the sites): each site has a reference to
+    // the runs it belongs to (which in turn has all the other references)
+    std::vector<std::shared_ptr<run>> runs;
 };
 
 #endif
diff --git a/src/run.cpp b/src/run.cpp
index d03acca..c47bf6b 100644
--- a/src/run.cpp
+++ b/src/run.cpp
@@ -1,55 +1,83 @@
 #include "run.h"
-#include "tsp_greedy.h"
+#include <stdexcept>
+#include <cassert>
 
-void run::combine(run &other_run, int new_vehicle, binop_dbl combine_load)
+void run::combine(run &other_run, VehicleTypeID new_vehicle)
 {
-    auto other_sites = other_run.sites();
-    _sites.insert(other_sites.begin(), other_sites.end());
+    assert(this->_distances == other_run._distances);
 
-    // new max load
-    double new_max_load = combine_load(this->max_load, other_run.max_load);
-    this->max_load = new_max_load;
-    other_run.max_load = new_max_load;
+    _sites.splice(_sites.end(), other_run._sites);
 
-    this->vehicle = new_vehicle;
+    this->_vehicle = new_vehicle;
+    this->_max_load = this->combined_max_load(other_run);
+    // needs to be after combined_max_load(), as it uses initial_load()
+    this->_initial_load += other_run._initial_load;
+    this->_final_load += other_run._final_load;
+
+    this->_distance =
+        this->_distance + other_run._distance -
+        this->_distances->get(0, 1 + *this->_sites.begin()) -
+        this->_distances->get(0, 1 + *other_run._sites.rbegin()) +
+        this->_distances->get(1 + *this->_sites.rbegin(),
+                              1 + *other_run._sites.begin());
 }
 
-std::vector<int> run::ordered_sites(const distmat<double> &distances) const
+double run::combined_max_load(const run &other_run) const
 {
-    return tsp_greedy(_sites, distances);
+    // note that this is not symmetric in the order of runs
+    return std::max(
+        this->_max_load + other_run._initial_load,
+        other_run._max_load + this->_final_load);
 }
 
-std::vector<int> run::ordered_sites(
-    const distmat<double> &distances,
-    const std::vector<int> &first,
-    const std::vector<int> &last) const
+std::map<Site, double> run::load_after_visit(const std::vector<double> &demand) const
 {
-    if (first.size() == 0 || last.size() == 0)
+    std::map<Site, double> loads;
+
+    double cur_load = _initial_load;
+    for (Site site : _sites)
     {
-        return ordered_sites(distances);
+        cur_load -= demand[site];
+        loads[site] = cur_load;
     }
 
-    // solve two tsp problems
-    std::unordered_set<int> sites_first;
-    std::unordered_set<int> sites_last;
-    for (const auto site : _sites)
+    return loads;
+}
+
+bool run::reassign_vehicle(Fleet &fleet)
+{
+    VehicleTypeID old_vehicle = this->_vehicle;
+
+    // we are guaranteed to find at least that vehicle once again
+    // so the usage of find_fitting_vehicle(...).value() is safe here
+    fleet.release_vehicle(old_vehicle);
+
+    VehicleTypeID vehicle =
+        fleet.find_fitting_vehicle(
+                 this->_sites,
+                 this->_max_load,
+                 false)
+            .value();
+
+    fleet.reserve_vehicle(vehicle);
+
+    if (old_vehicle != vehicle)
     {
-        if (first[site])
-        {
-            sites_first.insert(site);
-        }
-        if (last[site])
-        {
-            sites_last.insert(site);
-        }
+        this->_vehicle = vehicle;
+        return true;
     }
+    return false;
+}
 
-    std::vector<int> ordered_sites_first = tsp_greedy(sites_first, distances);
-    std::vector<int> ordered_sites_last = tsp_greedy(sites_last, distances);
-
-    ordered_sites_first.insert(
-        ordered_sites_first.end(),
-        ordered_sites_last.rbegin(), ordered_sites_last.rend());
+void run::optimize_route_order()
+{
+    std::list<Site> tsp_order;
+    double dist;
+    std::tie(tsp_order, dist) = tsp_greedy(_sites, *_distances);
 
-    return ordered_sites_first;
+    if (dist < this->_distance)
+    {
+        this->_sites = tsp_order;
+        this->_distance = dist;
+    }
 }
\ No newline at end of file
diff --git a/src/run.h b/src/run.h
index 6c36223..ba565be 100644
--- a/src/run.h
+++ b/src/run.h
@@ -1,45 +1,64 @@
 #ifndef RUN_H
 #define RUN_H
 
-#include <unordered_set>
-#include <functional>
-#include "distmat.h"
-
-using binop_dbl = std::function<double(const double, const double)>;
+#include <list>
+#include <map>
+#include <optional>
+#include <cmath>
+#include <memory>
+#include "site.h"
+#include "fleet.h"
+#include "symmat.h"
+#include "tsp_greedy.h"
 
 class run
 {
 public:
-  double max_load;
-  int vehicle;
-
-  run(int site, double max_load) // initialize a run with a single site
-      : max_load(max_load),
-        vehicle(-1),
-        _sites(std::unordered_set<int>{site}) {};
-  run(int site, double max_load, int vehicle)
-      : max_load(max_load),
-        vehicle(vehicle),
-        _sites(std::unordered_set<int>{site}) {};
-  run(std::unordered_set<int> &sites, double max_load, int vehicle)
-      : max_load(max_load),
-        vehicle(vehicle),
-        _sites(sites) {};
-  void combine(run &other_run, int new_vehicle, binop_dbl combine_load);
-  const std::unordered_set<int> &sites() const
+  run(Site site, double site_demand, VehicleTypeID vehicle,
+      std::shared_ptr<Distmat> distances)
+      : _initial_load(site_demand > 0 ? site_demand : 0),
+        _final_load(site_demand > 0 ? 0 : -site_demand),
+        _max_load(std::abs(site_demand)),
+        _sites(std::list<Site>{site}),
+        _vehicle(vehicle),
+        _distances(distances),
+        _distance(2 * distances->get(0, 1 + site)) {};
+
+  double combined_max_load(const run &other_run) const;
+  std::map<Site, double> load_after_visit(const std::vector<double> &demand) const;
+
+  void combine(run &other_run, VehicleTypeID new_vehicle);
+  bool reassign_vehicle(Fleet &fleet);
+  void optimize_route_order();
+
+  double distance() const
+  {
+    return _distance;
+  }
+  const std::list<Site> &sites() const
   {
     return _sites;
   }
-  std::vector<int> ordered_sites(const distmat<double> &distances) const;
-
-  // special version that has predefined order requirements:
-  // all "first" sites must come before "last" sites
-  std::vector<int> ordered_sites(const distmat<double> &distances,
-                                 const std::vector<int> &first,
-                                 const std::vector<int> &last) const;
+  double max_load() const
+  {
+    return _max_load;
+  };
+  const VehicleTypeID vehicle() const
+  {
+    return _vehicle;
+  }
 
 private:
-  std::unordered_set<int> _sites;
+  double _initial_load; // load when leaving the origin
+  double _final_load;   // load when arriving at the origin
+  double _max_load;
+
+  std::list<Site> _sites;
+
+  VehicleTypeID _vehicle;
+
+  std::shared_ptr<Distmat> _distances;
+  double _distance;
 };
 
 #endif
\ No newline at end of file
diff --git a/src/site.h b/src/site.h
new file mode 100644
index 0000000..feed24f
--- /dev/null
+++ b/src/site.h
@@ -0,0 +1,6 @@
+#ifndef SITE_H
+#define SITE_H
+
+typedef int Site;
+
+#endif
\ No newline at end of file
diff --git a/src/distmat.cpp b/src/symmat.cpp
similarity index 78%
rename from src/distmat.cpp
rename to src/symmat.cpp
index 634df79..660de12 100644
--- a/src/distmat.cpp
+++ b/src/symmat.cpp
@@ -1,23 +1,23 @@
-#include "distmat.h"
+#include "symmat.h"
 #include <cmath>
 #include <stdexcept>
 
 template <class T>
-distmat<T>::distmat(int size, T def)
+Symmat<T>::Symmat(int size, T def)
 {
   m_size = size;
   data = std::vector<T>(size * (size - 1) / 2, def);
 }
 
 template <class T>
-distmat<T>::distmat(const std::vector<T> &vec)
+Symmat<T>::Symmat(const std::vector<T> &vec)
 {
   m_size = std::round((1 + sqrt(1 + 8 * vec.size())) / 2);
   data = vec;
 }
 
 template <class T>
-T &distmat<T>::acc(const int i, const int j)
+T &Symmat<T>::acc(const int i, const int j)
 {
   // for i < j < n we have
   // n * i - i*(i+1)/2 + (j-i) - 1
@@ -42,7 +42,7 @@ T &distmat<T>::acc(const int i, const int j)
 }
 
 template <class T>
-T distmat<T>::get(const int i, const int j) const
+T Symmat<T>::get(const int i, const int j) const
 {
   if (i >= m_size)
   {
@@ -64,15 +64,15 @@ T distmat<T>::get(const int i, const int j) const
 }
 
 template <class T>
-int distmat<T>::size() const
+int Symmat<T>::size() const
 {
   return m_size;
 }
 
 template <class T>
-distmat<T> distmat<T>::sub(std::vector<int> &subvec) const
+Symmat<T> Symmat<T>::sub(std::vector<int> &subvec) const
 {
-  distmat<T> submat(subvec.size(), 0);
+  Symmat<T> submat(subvec.size(), 0);
 
   for (size_t i = 0; i < subvec.size(); i++)
   {
@@ -85,4 +85,4 @@ distmat<T> distmat<T>::sub(std::vector<int> &subvec) const
   return submat;
 }
 
-template class distmat<double>;
+template class Symmat<double>;
diff --git a/src/distmat.h b/src/symmat.h
similarity index 62%
rename from src/distmat.h
rename to src/symmat.h
index e4f1c0e..21157e3 100644
--- a/src/distmat.h
+++ b/src/symmat.h
@@ -4,19 +4,21 @@
 #include <vector>
 
 template <class T>
-class distmat
+class Symmat
 {
 public:
-  distmat(int n = 0, T def = T());
-  distmat(const std::vector<T> &vec);
+  Symmat(int n = 0, T def = T());
+  Symmat(const std::vector<T> &vec);
   T &acc(const int i, const int j);
   T get(const int i, const int j) const;
   int size() const;
-  distmat<T> sub(std::vector<int> &subvec) const;
+  Symmat<T> sub(std::vector<int> &subvec) const;
 
 private:
   std::vector<T> data;
   int m_size; // m_size * m_size elements in matrix
 };
 
+typedef Symmat<double> Distmat;
+
 #endif
diff --git a/src/tsp_greedy.cpp b/src/tsp_greedy.cpp
index 606c04e..c5175db 100644
--- a/src/tsp_greedy.cpp
+++ b/src/tsp_greedy.cpp
@@ -1,45 +1,50 @@
-#ifndef TSP_GREEDY
-#define TSP_GREEDY
-
 #include <vector>
 #include <unordered_set>
 #include <limits>
 #include <algorithm>
-#include "distmat.h"
+#include <list>
+#include "tsp_greedy.h"
+#include "site.h"
 
-std::vector<int> tsp_greedy(const std::unordered_set<int> sites,
-                            const distmat<double> &distances)
+template <typename Container>
+std::tuple<Container, double> tsp_greedy(const Container &sites, const Distmat &distances)
 {
-  std::vector<int> run;
-  run.reserve(sites.size());
+  Container run;
+  // run.reserve(sites.size());
+
+  double total_dist = 0;
 
   int ref_site = -1;
   int next_site = ref_site;
 
+  std::unordered_set<Site> not_visited;
+  for (const auto site : sites)
+  {
+    not_visited.insert(site);
+  }
+
   int cont = true;
   do
   {
     double min_dist = std::numeric_limits<double>::max();
     // get the nearest site from origin
-    for (const auto site : sites)
+    for (const auto site : not_visited)
     {
-      // if it's not already in the list check the distance
-      if (std::find(run.begin(), run.end(), site) == run.end())
-      {
-        double dist = distances.get(ref_site + 1, site + 1);
+      double dist = distances.get(ref_site + 1, site + 1);
 
-        if (dist < min_dist)
-        {
-          min_dist = dist;
-          next_site = site;
-        }
+      if (dist < min_dist)
+      {
+        min_dist = dist;
+        next_site = site;
       }
     }
 
     if (next_site != ref_site)
     {
       run.push_back(next_site);
+      not_visited.erase(next_site);
       ref_site = next_site;
+      total_dist += min_dist;
     }
     else
     {
@@ -48,7 +53,9 @@ std::vector<int> tsp_greedy(const std::unordered_set<int> sites,
 
   } while (cont);
 
-  return run;
+  total_dist += distances.get(0, next_site + 1);
+
+  return std::make_tuple(run, total_dist);
 }
 
-#endif
\ No newline at end of file
+template std::tuple<std::list<Site>, double> tsp_greedy(const std::list<Site> &sites, const Distmat &distances);
\ No newline at end of file
diff --git a/src/tsp_greedy.h b/src/tsp_greedy.h
index 9b03dae..ff8c689 100644
--- a/src/tsp_greedy.h
+++ b/src/tsp_greedy.h
@@ -1,6 +1,11 @@
+#ifndef TSP_GREEDY
+#define TSP_GREEDY
+
 #include <vector>
 #include <unordered_set>
-#include "distmat.h"
+#include "symmat.h"
+
+template <typename Container>
+std::tuple<Container, double> tsp_greedy(const Container &sites, const Distmat &distances);
 
-std::vector<int> tsp_greedy(const std::unordered_set<int> sites,
-                            const distmat<double> &distances);
+#endif
\ No newline at end of file
diff --git a/tests/testthat/test-clarke_wright.R b/tests/testthat/test-clarke_wright.R
index b399071..b647f0c 100644
--- a/tests/testthat/test-clarke_wright.R
+++ b/tests/testthat/test-clarke_wright.R
@@ -164,7 +164,7 @@ test_that("Vehicles are not assigned to restricted sites", {
         )
 
       runs <- res$visits[res$visits$site == 0, "run"]
-      
+
       expect_false(
         0 %in% res$runs[res$runs$run %in% runs, "vehicle"]
       )
@@ -341,22 +341,18 @@ test_that("Truck loads are always within physical boundaries at any point on the
           demand_net$distances,
           data.frame(n = c(NA_integer_, 3L), caps = c(60, max_cap))
         )
-      
+
       expect_all_true(res$visits$load >= 0)
       expect_all_true(res$visits$load <= max_cap)
     }
   )
 })
 
-test_that("Max of positive demand sum and negative demand sum
-          for each run equals the max_load", {
+test_that("Single-signed demand: Max of demand sum for each run equals max_load", {
   skip_if_not_installed("hedgehog")
 
-  # requirement for that to be true:
-  # * only count positive demands
-  # * demand is always <= vehicle capacity
   hedgehog::forall(
-    gen.demand_net(max_sites = 10L),
+    gen.demand_net(max_sites = 10L, min_demand = 1L),
     function(demand_net) {
       res <-
         clarke_wright(
@@ -371,19 +367,42 @@ test_that("Max of positive demand sum and negative demand sum
         by = "site"
       )
 
-      pos_load <-
+      load <-
         as.numeric(
-          by(res1, res1$run, function(x) sum(pmax(x$demand, 0)))
+          by(res1, res1$run, function(x) sum(x$demand))
+        )
+
+      expect_equal(
+        res$runs$max_load,
+        abs(load)
+      )
+    }
+  )
+
+  hedgehog::forall(
+    gen.demand_net(max_sites = 10L, max_demand = -1L),
+    function(demand_net) {
+      res <-
+        clarke_wright(
+          demand_net$demand,
+          demand_net$distances,
+          data.frame(n = c(NA_integer_, 3L), caps = c(60, 120))
         )
 
-      neg_load <-
+      res1 <- merge(
+        res$visits,
+        res$sites,
+        by = "site"
+      )
+
+      load <-
         as.numeric(
-          by(res1, res1$run, function(x) sum(pmin(x$demand, 0)))
+          by(res1, res1$run, function(x) sum(x$demand))
         )
 
       expect_equal(
         res$runs$max_load,
-        pmax(pos_load, -neg_load)
+        abs(load)
       )
     }
   )
@@ -407,7 +426,7 @@ test_that("Example scenario with negative demand yields a single run", {
 
 
 test_that("Last result of clarke_wright stepwise should be the original result", {
-skip_if_not_installed("hedgehog")
+  skip_if_not_installed("hedgehog")
 
   # requirement for that to be true:
   # * only count positive demands