Ticket #296: ticket296.patch

lemon/multicommodity_flow.h

@@ -2 +2 @@
  *
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2003-2009
+ * Copyright (C) 2003-2021
  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
  * (Egervary Research Group on Combinatorial Optimization, EGRES).
  *
@@ -14 +14 @@
  * express or implied, and with no claim as to its suitability for any
  * purpose.
  *
+ * Additional Copyright file multicommodity_flow.h: Lancaster University (2021)
  */
 
-#ifndef LEMON_MULTICOMMODITY_FLOW_H
-#define LEMON_MULTICOMMODITY_FLOW_H
+#ifndef LEMON_MULTICOMMODITY_FLOW_H_
+#define LEMON_MULTICOMMODITY_FLOW_H_
 
 /// \ingroup approx
 ///
 /// \file
-/// \brief Approximation algorithms for solving multicommodity flow problems.
+/// \brief Fleischer's algorithms for finding approximate multicommodity flows.
 
+#include <algorithm>  // min
+#include <limits>     // numeric_limits
+#include <random>     // random_device, mt19937
+#include <set>
 #include <vector>
-#include <lemon/bfs.h>
-#include <lemon/dijkstra.h>
-#include <lemon/preflow.h>
-#include <lemon/adaptors.h>
-#include <lemon/tolerance.h>
-#include <lemon/path.h>
-#include <lemon/math.h>
 
-//
-// TODO:
-// 
-//   - The algorithms should provide primal and dual solution as well:
-//       primal objective value, primal solution (flow), primal solution (flow) for each commodity
-//       dual objective value, dual solution
-//     Note that these algorithms are approximation methods, so the found primal and dual obj. values
-//     are not necessarily optimal, ie. the found primal obj. value could be smaller than the found
-//     dual objective value and their ratio provides a proof for the approximation factor of both
-//     the primal and the dual solution.
-//   
-//   - compute an upper bound for the optimal primal obj. value, e.g. using
-//       x_i = min ( out_cap(s_i), in_cap(t_i) )  for each commodity i
-//     the sum of x_i or the minimum of x_i/d_i ratios can be utilized.
-//     
-//   - give back as good approximation guarantee r as possible for which:
-//     dual = r*primal
-//
-//   - revise the implementations and their API
-//
+#include "lemon/adaptors.h"
+#include "lemon/bfs.h"
+#include "lemon/concepts/maps.h"
+#include "lemon/dijkstra.h"
+#include "lemon/maps.h"
+#include "lemon/math.h"
+#include "lemon/path.h"
+#include "lemon/preflow.h"
+#include "lemon/static_graph.h"
+#include "lemon/tolerance.h"
 
 namespace lemon {
 
-  /// \addtogroup approx
+/// \brief Default traits class of \ref ApproxMCF class.
+///
+/// Default traits class of \ref ApproxMCF class.
+/// \tparam GR The type of the digraph.
+/// \tparam CAPType The type of the capacity values.
+/// \tparam CAPMap The type of the capacity map.
+/// \tparam V The type for the flows and costs.
+/// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
+#ifdef DOXYGEN
+template<typename GR, typename CAPType, typename CAPMap, typename V>
+#else
+template<
+    typename GR,
+    typename CAPMap = typename GR::template ArcMap<double>,
+    typename V      = double>
+#endif
+struct ApproxMCNFDefaultTraits {
+  /// The type of the digraph
+  typedef GR Digraph;
+  /// The type of the capacity map
+  typedef CAPMap CapacityMap;
+  /// The type of the capacity map values. By default is double.
+  typedef typename CapacityMap::Value CapacityValue;
+
+  /// \brief The type used for flow values and costs.
+  ///
+  /// \c CapacityValue must be convertible to \c Value.
+  /// By default is \c double.
+  typedef V                                        Value;
+  typedef typename Digraph::template ArcMap<Value> ValueMap;
+
+  /// The tolerance type used for internal computations
+  typedef lemon::Tolerance<Value> Tolerance;
+
+  /// \brief The path type of the decomposed flows
+  ///
+  /// It must conform to the \ref lemon::concepts::Path "Path" concept
+  typedef lemon::Path<Digraph>                 Path;
+  typedef typename lemon::Path<Digraph>::ArcIt PathArcIt;
+
+  /// \brief Instantiates a ValueMap.
+  ///
+  /// This function instantiates a \ref ValueMap.
+  /// \param digraph The digraph for which we would like to define
+  /// the large cost map. This can be used for flows or otherwise.
+  static ValueMap* createValueMapPtr(
+      const Digraph& digraph,
+      const Value&   value = 0) {
+    return new ValueMap(digraph, value);
+  }
+  // Same as above just const
+  static const ValueMap* createValueMapConstPtr(
+      const Digraph& digraph,
+      const Value&   value = 0) {
+    return new ValueMap(digraph, value);
+  }
+};
+
+/// \addtogroup approx
+/// @{
+/// \brief Fleischer's algorithms for the multicommodity flow problem.
+///
+/// Default traits class of \ref ApproxMCF class.
+/// \tparam GR The type of the digraph.
+/// \tparam TR The type of the traits class.
+#ifdef DOXYGEN
+template<typename GR, typename TR>
+#else
+template<
+    typename GR,
+    typename CM = typename GR::template ArcMap<double>,
+    typename TR = ApproxMCNFDefaultTraits<GR, CM>>
+#endif
+class ApproxMCF {
+ public:
+  /// The type of the digraph
+  typedef typename TR::Digraph Digraph;
+  /// The type of the capacity map
+  typedef typename TR::CapacityMap CapacityMap;
+  /// The type of the capacity values
+  typedef typename TR::CapacityValue CapacityValue;
+
+  /// \brief The large cost type
+  ///
+  /// The large cost type used for internal computations.
+  /// By default is \c double.
+  typedef typename TR::Value Value;
+  /// Arc Map with values of type \ref Value
+  typedef typename TR::ValueMap ValueMap;
+
+  /// The tolerance type
+  typedef typename TR::Tolerance Tolerance;
+
+  /// \brief The path type of the found cycles
+  ///
+  /// The path type of the found cycles.
+  /// Using the \ref lemon::ApproxMCNFDefaultTraits "default traits class",
+  /// it is \ref lemon::Path "Path<Digraph>".
+  typedef typename TR::Path Path;
+  /// The arc iterator for paths.
+  typedef typename TR::PathArcIt PathArcIt;
+
+  /// The \ref lemon::ApproxMCNFDefaultTraits "traits class" of the
+  /// algorithm
+  typedef TR Traits;
+
+  /// Types of approximation type.
+  enum ApproxType {
+    /// Maximum flow using Fleischer's algorithm.
+    FLEISCHER_MAX,
+    /// Maximum concurrent flow using Fleischer's algorithm.
+    FLEISCHER_MAX_CONCURRENT,
+    /// Minimum cost concurrent flow using binary search with cost-bounded
+    /// \ref FLEISCHER_MAX_CONCURRENT
+    BINARY_SEARCH_MIN_COST
+  };
+
+ private:
+  struct PathData {
+    Path path;
+    // The amount of flow through the path
+    Value value;
+    // The minimum capacity along all edges in the path
+    CapacityValue min_cap;
+    // The actual cost of the path
+    Value cost;
+    int   hash;
+
+    PathData() {}
+    PathData(const Path& p, Value v, CapacityValue m, Value c, int h = 0)
+        : path(p), value(v), min_cap(m), cost(c), hash(h) {}
+  };
+
+  int hash(const Path& path) const {
+    int h = path.length();
+    for (PathArcIt a(path); a != INVALID; ++a)
+      h = (h << 4) ^ (h >> 28) ^ _graph.id(a);
+    return h;
+  }
+
+  // Matrix indexed by commodity with all the paths found. i.e. [i][p] where i
+  // is the commodity index and p is the path index.
+  typedef std::vector<std::vector<PathData>> PathMatrix;
+  // Vector indexed by commodity with \ref ValueMap pointers.
+  typedef std::vector<ValueMap*>       ValueMapVector;
+  typedef std::vector<const ValueMap*> ConstValueMapVector;
+
+  TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+  const Digraph& _graph;
+  const int      _num_nodes;
+  const int      _num_arcs;
+  // Input capacity map
+  const CapacityMap& _cap;
+
+  // The number of commodities
+  int _num_commodities;
+  // The source nodes for each commodity
+  std::vector<Node> _source;
+  // The target nodes for each commodity
+  std::vector<Node> _target;
+  // The demands for each commodity (may be 0, or undefined).
+  std::vector<CapacityValue> _demand;
+  // The cost for each commodity (may be 0, or undefined)
+  // See \ref ValueMapVector.
+  ConstValueMapVector _cost;
+
+  /**
+   * Algorithm parameters
+   */
+  // type of approximation algorithm. See \ref ApproxType.
+  ApproxType _approx_type;
+  // Stopping criterion
+  Value _epsilon;
+  // Whether the internal heuristics will be applied at the end of the
+  // algorithms
+  bool _heuristic;
+  // Whether there is feasible unsplitable flow
+  bool _feasible_unsplit;
+  // Objective function of the primal problem
+  Value _obj;
+  // Objective function of the dual problem
+  Value _dual_obj;
+  bool  _local_cost;
+  //
+  Value _delta;
+  // Expected maximum number of iterations
+  Value _max_iter;
+  //
+  Value _limit;
+  // Number of iterations
+  int _iter;
+  // Tolerance function. See \ref lemon::Tolerance.
+  Tolerance _tol;
+
+  // For the BINARY_SEARCH_MIN_COST case
+  // Cost of the final solution.
+  Value _total_cost;
+  // Cost bounding factor.
+  Value _phi;
+  // Cost bound.
+  Value _cost_bound;
+
+  // The aggregated flow per arc over all commodities
+  ValueMap _total_flow;
+  // The value of the length per arc (the dual solution).
+  ValueMap _len;
+  // The flow for each commodity. See \ref ValueMapVector.
+  ValueMapVector _flow;
+  // The paths found for different commodities. See \ref PathMatrix.
+  PathMatrix _paths;
+
+  // Solution Attributes
+  std::vector<Path>  _final_paths;
+  std::vector<Value> _final_flows;
+  std::vector<Value> _final_costs;
+
+ public:
+  /// \brief Default constructor
+  ApproxMCF(const Digraph& gr, const CapacityMap& cp)
+      : _graph(gr),
+        _num_nodes(countNodes(_graph)),
+        _num_arcs(countArcs(_graph)),
+        _cap(cp),
+        _num_commodities(0),
+        _feasible_unsplit(true),
+        _obj(0),
+        _dual_obj(0),
+        _local_cost(true),
+        _iter(0),
+        _cost_bound(0),
+        _total_flow(_graph),
+        _len(_graph) {}
+
+  ~ApproxMCF() {
+    for (int i = 0; i < _num_commodities; ++i) {
+      delete _flow[i];
+      _flow[i] = nullptr;
+      if (_local_cost) {
+        delete _cost[i];
+        _cost[i] = nullptr;
+      }
+    }
+    _flow.clear();
+    _cost.clear();
+  }
+
+  /// \name Parameters
+  /// The parameters of the algorithm can be specified using these
+  /// functions.
+
   /// @{
 
-  /// \brief Implementation of an approximation algorithm for finding
-  /// a maximum multicommodity flow.
+  /// \brief Set single commodity parameters.
   ///
-  /// \ref MaxMulticommodityFlow implements Lisa Fleischer's algorithm
-  /// for finding a maximum multicommodity flow.
+  /// \param s The source node.
+  /// \param t The target node.
+  /// \param d The required amount of flow from node \c s to node \c t
+  /// (i.e. the supply of \c s and the demand of \c t).
+  /// \param cost_map The cost arc map.
   ///
-  /// \tparam Digraph The type of the digraph the algorithm runs on.
-  /// \tparam CAP The type of the capacity map.
-  /// The default map type is \ref concepts::Digraph::ArcMap
-  /// "GR::ArcMap<double>".
-  /// \tparam V The value type used in the algorithm.
-  /// The default value type is \c CAP::Value.
-#ifdef DOXYGEN
-  template < typename GR,
-             typename CAP,
-             typename V >
-#else
-  template < typename GR,
-             typename CAP = typename GR::template ArcMap<double>,
-             typename V = typename CAP::Value >
-#endif
-  class MaxMulticommodityFlow
-  {
-    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
-
-    typedef SimplePath<GR> SPath;
-    typedef typename SPath::ArcIt PathArcIt;
-
-  public:
-
-    ///  The type of the digraph the algorithm runs on.
-    typedef GR Digraph;
-    /// The type of the capacity map.
-    typedef CAP CapacityMap;
-    /// The value type of the algorithm.
-    typedef V Value;
-
-#ifdef DOXYGEN
-    /// The type of the flow map (primal solution).
-    typedef Digraph::ArcMap<Value> FlowMap;
-    /// The type of the length map (dual soluiton).
-    typedef Digraph::ArcMap<Value> LengthMap;
-#else
-    /// The type of the flow map (primal solution).
-    typedef typename Digraph::template ArcMap<Value> FlowMap;
-    /// The type of the length map (dual soluiton).
-    typedef typename Digraph::template ArcMap<Value> LengthMap;
-#endif
-
-  private:
-
-    // The digraph the algorithm runs on
-    const Digraph &_g;
-    // The capacity map
-    const CapacityMap &_cap;
-
-    // The number of commdities
-    int _k;
-    // The source nodes for each commodity
-    std::vector<Node> _s;
-    // The sink nodes for each commodity
-    std::vector<Node> _t;
-    // The weight of each commodity
-    // std::vector<Weight> _w;
+  /// \return <tt>(*this)</tt>
+  ApproxMCF& addCommodity(
+      const Node&          s,
+      const Node&          t,
+      const CapacityValue& d,
+      const ValueMap&      cost_map) {
+    _source.push_back(s);
+    _target.push_back(t);
+    if (d > 0)
+      _demand.push_back(d);
+    _local_cost = false;
+    _cost.push_back(&cost_map);
+    ++_num_commodities;
+    return *this;
+  }
 
-    // The flow map (the primal solution).
-    FlowMap _flow;
-    // The length map (the dual solution).
-    LengthMap _len;
-
-    struct PathData {
-      SPath path;
-      Value value;
-      Value min_cap;
-      int hash;
-
-      PathData() {}
-      PathData(const SPath& p, Value v, Value mc, int h = 0) :
-        path(p), value(v), min_cap(mc), hash(h) {}
-    };
-
-  public:
-
-    /// \brief Constructor.
-    ///
-    /// Constructor.
-    ///
-    /// \param digraph The digraph the algorithm runs on.
-    /// \param capacity The capacities of the edges.
-    MaxMulticommodityFlow( const Digraph &digraph,
-                           const CapacityMap &capacity ) :
-      _g(digraph), _cap(capacity), _k(0),
-      _flow(_g), _len(_g)
-    {}
+  /// @overload no cost map, optional demand
+  ApproxMCF&
+  addCommodity(const Node& s, const Node& t, const CapacityValue& d = 0) {
+    const ValueMap* cost_map_ptr = Traits::createValueMapConstPtr(_graph);
+    addCommodity(s, t, d, *cost_map_ptr);
+    _local_cost = true;
+    return *this;
+  }
 
-    /// \brief Constructor.
-    ///
-    /// Constructor.
-    ///
-    /// \param gr The digraph the algorithm runs on.
-    /// \param capacity The capacities of the edges.
-    /// \param k The number of commodities.
-    /// \param sources A vector containing the source nodes
-    /// (its size must be at least \c k).
-    /// \param sinks A vector containing the sink nodes
-    /// (its size must be at least \c k).
-    MaxMulticommodityFlow( const Digraph &gr,
-                           const CapacityMap &capacity,
-                           int k,
-                           const std::vector<Node>& sources,
-                           const std::vector<Node>& sinks ) :
-      _g(gr), _cap(capacity),
-      _k(k), _s(sources), _t(sinks),
-      _flow(_g), _len(_g)
-    {
-      _s.resize(k);
-      _t.resize(k);
-    }
+  /// @}
 
-    /// Destructor.
-    ~MaxMulticommodityFlow() {}
+  /// \name Execution control
 
-    /// \brief Add a new commodity.
-    ///
-    /// This function adds a new commodity with the given source and
-    /// sink nodes.
-    void addCommodity(const Node& s, const Node& t) {
-      ++_k;
-      _s.push_back(s);
-      _t.push_back(t);
-    }
-
-    /// \name Execution control
-
-    /// @{
+  /// @{
 
-    /// \brief Run the algorithm for finding a maximum multicommodity flow.
-    ///
-    /// This function runs the algorithm for finding a maximum multicommodity
-    /// flow.
-    /// \param r The approximation factor. It must be greater than 1.
-    /// \return The effective approximation factor (at least 1), i.e. the
-    /// ratio of the values of the found dual and primal solutions.
-    /// Therefore both the primal and dual solutions are r-approximate.
-    ///
-    /// \note The smaller parameter \c r is, the slower the algorithm runs,
-    /// but it might find better solution. According to our tests, using
-    /// the default value the algorithm performs well and it usually manage
-    /// to find the exact optimum due to heuristic improvements.
-    Value run(Value r = 2.0, Value opt = 0) {
-
-    #define MMF_FLEISCHER_IMPR
-    #define MMF_HEURISTIC_IMPR
-    //#define MMF_EXP_UPDATE
-    //#define MMF_DEBUG_OUTPUT
-
-      // Initialize epsilon and delta parameters
-      Tolerance<Value> ftol, ctol;
-      if (!ftol.less(1.0, r)) return 0;
-#ifdef MMF_FLEISCHER_IMPR
-  #ifdef MMF_EXP_UPDATE
-      Value epsilon = (r + 1 - 2 * pow(r, 0.5)) / (r - 1);
-  #else
-      Value epsilon = (2 * r + 1 - pow(8 * r + 1, 0.5)) / (2 * r);
-  #endif
-#else
-  #ifdef MMF_EXP_UPDATE
-      Value epsilon = (2 * r + 1 - pow(8 * r + 1, 0.5)) / (2 * r);
-  #else
-      Value epsilon = 1 - pow(r, -0.5);
-  #endif
-#endif
-      
-#ifdef MMF_DEBUG_OUTPUT      
-      std::cout << "r = " << r << "; epsilon = " << epsilon << "\n";
-#endif
-      
-      int node_num = countNodes(_g);
-      Value delta = pow(1 + epsilon, 1 - 1/epsilon) /
-                    pow(node_num, 1/epsilon);
-      ctol.epsilon(delta * 0.0001);
-
-      // Initialize flow values and lengths
-      for (ArcIt a(_g); a != INVALID; ++a) {
-        _flow[a] = 0;
-        _len[a] = delta;
-      }
-
-      // Check commodities
-      for (int i = 0; i < _k; ++i) {
-        if (_s[i] == _t[i] || !bfs(_g).run(_s[i], _t[i])) {
-          _s[i] = _s[_k - 1];
-          _t[i] = _t[_k - 1];
-          --_k;
-        }
-      }
-      if (_k <= 0) return 0;
-
+  /// \brief Run the algorithm for finding a multicommodity flow.
+  ///
+  /// This function runs the algorithm for finding a maximum/minimum cost
+  /// multicommodity flow.
+  ///
+  /// \param at  The type of approximation algorithm to use. @see \ref
+  /// ApproxMCF::ApproxType.
+  /// \param heuristic Whether the internal heuristics will be applied at the
+  // end of Fleischer's algorithm.
+  /// \param epsilon The stopping criterion. Must be positive and larger than
+  /// 0. The smaller the value, the more accurate the result, but the longer the
+  /// execution time. By default is 0.1.
+  void run(
+      const ApproxType& at        = FLEISCHER_MAX,
+      const bool&       heuristic = true,
+      const Value&      epsilon   = 0.1) {
+    _approx_type = at;
+    _epsilon     = epsilon;
+    _heuristic   = heuristic;
+    check();
+    switch (_approx_type) {
+      case FLEISCHER_MAX:
+        runFleischerMax();
+        break;
+      case BINARY_SEARCH_MIN_COST:
+        runBinarySearchMinCost();
+        break;
+      case FLEISCHER_MAX_CONCURRENT:
+        runFleischerMaxConcurrent();
+        break;
+    }
+  }
 
-      // Successive improvement along shortest paths
-      std::vector<std::vector<PathData> > paths(_k);
-#ifdef MMF_FLEISCHER_IMPR
-      Value limit = delta * (1 + epsilon);
-      int i = 0;
-#else
-      std::vector<SPath> curr_paths(_k);
-#endif
-      SPath curr_path;
-      Value curr_dist;
-      int iter = 0;
-      while (true) {
-#ifdef MMF_FLEISCHER_IMPR
-        // Find the shortest path for the i-th commodity
-        // (w.r.t the current length function)
-        dijkstra(_g, _len).path(curr_path).dist(curr_dist)
-          .run(_s[i], _t[i]);
-        
-        // Check if the path can be used and increase i or the
-        // length limit if necessary
-        if (!ctol.less(curr_dist, limit)) {
-          if (++i == _k) {
-            if (!ctol.less(limit, 1.0)) break;
-            i = 0;
-            limit *= (1 + epsilon);
-            if (ctol.less(1.0, limit)) limit = 1.0;
-          }
-          continue;
-        }
-#else
-        // Find the shortest path
-        Value min_dist = 1.0;
-        int i = -1;
-        for (int j = 0; j < _k; ++j) {
-          if ( !dijkstra(_g, _len).path(curr_paths[j]).dist(curr_dist).
-                run(_s[j], _t[j]) ) {
-            std::cerr << "\n\nERROR!!! "<<j<<":"<<_k<<"\n\n";
-          }
-          if (curr_dist < min_dist) {
-            min_dist = curr_dist;
-            i = j;
-          }
-        }
-        if (i < 0) break;
-        curr_path = curr_paths[i];
-#endif
+  /// @}
+
+ private:
+  // Initialize flow values and lengths
+  void init() {
+    _iter     = 0;
+    _obj      = 0;
+    _dual_obj = 0;
+    // Allocate paths
+    _paths.resize(_num_commodities);
+    // Solution paths
+    _final_paths.resize(_num_commodities);
+    _final_flows.resize(_num_commodities);
+    _final_costs.resize(_num_commodities);
 
-        ++iter;
-
-        // Check if the path is used before
-        int h = hash(curr_path);
-        int j;
-        for (j = 0; j < int(paths[i].size()); ++j) {
-          bool b = (paths[i][j].hash == h) &&
-                   (paths[i][j].path.length() == curr_path.length());
-          if (b) {
-            for (PathArcIt a1(paths[i][j].path), a2(curr_path);
-                 b && a1 != INVALID && a2 != INVALID; ++a1, ++a2)
-              b = _g.id(a1) == _g.id(a2);
-            if (b) break;
-          }
-        }
+    for (int i = 0; i < _num_commodities; ++i)
+      _paths[i].clear();
+    _delta =
+        pow(1 + _epsilon, 1 - 1 / _epsilon) / pow(_num_nodes, 1 / _epsilon);
+    _tol.epsilon(_delta * 1e-5);
+    _limit = _delta * (1 + _epsilon);
 
-        // Set/modify the flow value for the found path
-        Value gamma;
-        if (j == int(paths[i].size())) {
-          gamma = std::numeric_limits<Value>::infinity();
-          for (PathArcIt a(curr_path); a != INVALID; ++a) {
-            if (_cap[a] < gamma) gamma = _cap[a];
-          }
-          paths[i].push_back(PathData(curr_path, gamma, gamma, h));
-        } else {
-          gamma = paths[i][j].min_cap;
-          paths[i][j].value += gamma;
-        }
-
-        // Modify the arc lengths along the current path
-#ifdef MMF_EXP_UPDATE
-        for (PathArcIt a(curr_path); a != INVALID; ++a)
-          _len[a] *= exp(epsilon * gamma / _cap[a]);
-#else
-        for (PathArcIt a(curr_path); a != INVALID; ++a)
-          _len[a] *= 1 + epsilon * gamma / _cap[a];
-#endif
-      }
-
-      // Setting flow values
-      for (int i = 0; i < _k; ++i) {
-        for (int j = 0; j < int(paths[i].size()); ++j) {
-          Value val = paths[i][j].value;
-          for (PathArcIt a(paths[i][j].path); a != INVALID; ++a)
-            _flow[a] += val;
+    // Allocate/reset flow per commodity and arc
+    _flow.resize(_num_commodities);
+    for (int i = 0; i < _num_commodities; ++i)
+      if (!_flow[i]) {
+        _flow[i] = Traits::createValueMapPtr(_graph);
+      } else {
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          (*_flow[i])[a] = 0;
         }
       }
 
-      // Scaling the solution
-      Value lambda = 0;
-      for (ArcIt a(_g); a != INVALID; ++a) {
-        if (_flow[a] / _cap[a] > lambda) lambda = _flow[a] / _cap[a];
+    // Initialise length and flow map
+    for (ArcIt a(_graph); a != INVALID; ++a) {
+      _len[a]        = _delta;
+      _total_flow[a] = 0;
+    }
+
+    // Initialise phi for cost bounding
+    if (_approx_type == BINARY_SEARCH_MIN_COST && _cost_bound > 0)
+      _phi = _delta / _cost_bound;
+    else
+      _phi = 0;
+  }
+
+  // Preprocessing checks to ensure that:
+  //
+  // 1. Commodities are well-defined
+  // 2. Input sizes are consistent
+  // 3. Concurrent case: Feasible split flow exists for each commodity that
+  // satisfies demand. In the case this doesn't pass, will throw an exception.
+  //
+  // TODO(): Group commodities by source, so we only have to solve the shortest
+  // path once for each group (and length values).
+  void check() {
+    LEMON_ASSERT(_tol.positive(_epsilon), "Epsilon must be greater than 0");
+
+    // Remove commodities with source = sink or unreachable
+    for (int i = 0; i < _num_commodities; ++i) {
+      if (_source[i] == _target[i] ||
+          !bfs(_graph).run(_source[i], _target[i])) {
+        _source[i] = _source[_num_commodities - 1];
+        _target[i] = _target[_num_commodities - 1];
+        --_num_commodities;
+      }
+    }
+
+    const int& s_size = _source.size();
+    const int& t_size = _target.size();
+    const int& d_size = _demand.size();
+
+    switch (_approx_type) {
+      case FLEISCHER_MAX: {
+        LEMON_ASSERT(_tol.positive(_num_commodities), "No commodities");
+        // Check consistency of vector sizes
+        LEMON_ASSERT(
+            (_num_commodities == s_size && s_size == t_size &&
+             !_tol.positive(d_size)),
+            "All vectors should be of size equal to the number of commodities, "
+            "and have no demand.");
+        break;
       }
-      if (ftol.positive(lambda)) {
-        lambda = 1 / lambda;
-        for (ArcIt a(_g); a != INVALID; ++a) _flow[a] *= lambda;
-        for (int i = 0; i < _k; ++i) {
-          for (int j = 0; j < int(paths[i].size()); ++j)
-            paths[i][j].value *= lambda;
+      case FLEISCHER_MAX_CONCURRENT:
+      case BINARY_SEARCH_MIN_COST: {
+        // Check consistency of vector sizes
+        LEMON_ASSERT(
+            (_num_commodities == s_size && s_size == t_size &&
+             s_size == d_size),
+            "All vectors should be of size equal to the number of commodities");
+        // Remove commodities with 0 demand.
+        for (int i = 0; i < _num_commodities; ++i) {
+          if (!_tol.positive(_demand[i])) {
+            _source[i] = _source[_num_commodities - 1];
+            _target[i] = _target[_num_commodities - 1];
+            --_num_commodities;
+          }
+        }
+        LEMON_ASSERT(_tol.positive(_num_commodities), "No commodities");
+
+        // Check if feasible routing of demand exists for each commodity
+        for (int i = 0; i < _num_commodities; ++i) {
+          // Only check commodity if the demand in non-zero
+          // Check if feasible split flow exists
+          // Total capacity outgoing at the source
+          CapacityValue total_cap_source = 0;
+          CapacityValue max_cap_source   = 0;
+          for (OutArcIt a(_graph, _source[i]); a != INVALID; ++a) {
+            total_cap_source += _cap[a];
+            if (max_cap_source < _cap[a])
+              max_cap_source = _cap[a];
+          }
+
+          // Total capacity incoming at the target
+          CapacityValue total_cap_target = 0;
+          CapacityValue max_cap_target   = 0;
+          for (InArcIt a(_graph, _target[i]); a != INVALID; ++a) {
+            total_cap_target += _cap[a];
+            if (max_cap_target < _cap[a])
+              max_cap_target = _cap[a];
+          }
+
+          // Check if enough capacity to satisfy the demand
+          LEMON_ASSERT(
+              (total_cap_source >= _demand[i] &&
+               total_cap_target >= _demand[i]),
+              "No feasible (split) flow to satisfy demand for commodity" +
+                  std::to_string(i));
+
+          // Check if feasible unsplittable flow exists that satisfies demand
+          if (max_cap_source < _demand[i] || max_cap_target < _demand[i])
+            _feasible_unsplit = false;
         }
       }
+    }
+  }
 
-#ifdef MMF_DEBUG_OUTPUT
-      Value value1 = 0;
-      for (int i = 0; i < _k; ++i) {
-        for (InArcIt a(_g, _t[i]); a != INVALID; ++a) value1 += _flow[a];
-        for (OutArcIt a(_g, _t[i]); a != INVALID; ++a) value1 -= _flow[a];
+  // Run standard max concurrent algorithm iteratively updating the budget, or
+  // cost bound, a la binary search to obtain an approximation to the
+  // multicommodity flow minimum cost
+  void runBinarySearchMinCost() {
+    // Run standard max concurrent
+    runFleischerMaxConcurrent();
+
+    Value lower_bound = 0, upper_bound = getTotalCost();
+
+    // Best solution attributes
+    Value          curr_best_cost = std::numeric_limits<Value>::infinity();
+    ValueMapVector curr_best_flow;
+
+    while ((upper_bound - lower_bound) / upper_bound > _epsilon &&
+           upper_bound >= lower_bound) {
+      // Run Max Concurrent algorithm with cost bound
+      _cost_bound = (upper_bound + lower_bound) / 2;
+      runFleischerMaxConcurrent();
+      // Update upper/lower bounds
+      if (_total_cost < _cost_bound)
+        upper_bound = _cost_bound;
+      else
+        lower_bound = _cost_bound;
+      // Save current solution if appropriate
+      if (_total_cost < curr_best_cost) {
+        curr_best_cost = _total_cost;
+        curr_best_flow = _flow;
       }
-      Value max_opt = value1 / (1 - epsilon);
+    }
+    // Update flow
+    _total_cost = curr_best_cost;
+    _flow       = curr_best_flow;
+  }
 
-      std::cout << "\n";
-      std::cout << "Computation finished. Iterations: " << iter << "\n";
-      std::cout << "Solution value:   " << value1;
-      if (opt > 0) std::cout << "  \t" << value1 / opt * 100 << "%";
-      std::cout << std::endl;
-#else
-      opt = opt; // avoid warning
-#endif
+  void runFleischerMaxConcurrent() {
+    // Successive improvement along shortest paths
+    init();
+    Path  curr_path;
+    Value curr_dist;
+    while (true) {
+      // Find the shortest path for the i-th commodity
+      // (w.r.t the current length function)
+      for (int i = 0; i < _num_commodities; ++i) {
+        // Route demand
+        CapacityValue d_i = _demand[i];
+        while (true) {  // while (phase)
+          dijkstra(_graph, _len)
+              .path(curr_path)
+              .dist(curr_dist)
+              .run(_source[i], _target[i]);
+          // Get
+          const CapacityValue& routed_flow = updatePaths(curr_path, i, d_i);
+          // Update demand
+          d_i -= routed_flow;
+          updateLen(curr_path, routed_flow, i);
+          updateDualObjective();
+          if (d_i <= 0 || !_tol.less(_dual_obj, 1))
+            break;
+        }  // end while (phase)
+      }
+      ++_iter;
+      updateDualObjective();
+      if (!_tol.less(_dual_obj, 1))
+        break;
+    }  // end while (main)
+
+    // Clean up
+    scaleFinalFlowsMaxConcurrent();
+
+    if (_heuristic)
+      randomisedRounding();
+
+    setFlowPerCommodity();
+    setFinalCost();
+    saveFinalPaths();
+  }
+
+  void runFleischerMax() {
+    // Successive improvement along shortest paths
+    int   i = 0;
+    Path  curr_path;
+    Value curr_dist;
+    init();
+    while (true) {  // start while
+      dijkstra(_graph, _len)
+          .path(curr_path)
+          .dist(curr_dist)
+          .run(_source[i], _target[i]);
 
-#ifdef MMF_DEBUG_OUTPUT
-      Value tmp = 0;
-      for (ArcIt a(_g); a != INVALID; ++a)
-        tmp += _cap[a] * _len[a];
-      std::cout << "Final dual: " << tmp << "\n";
-#endif
+      // Check if the path can be used and increase i or the length limit if
+      // necessary
+      if (!_tol.less(curr_dist, _limit)) {
+        if (++i == _num_commodities) {
+          if (!_tol.less(_limit, 1.0))
+            break;
+          i = 0;
+          _limit *= (1 + _epsilon);
+          if (_tol.less(1.0, _limit))
+            _limit = 1.0;
+        }
+        continue;
+      }
+      const CapacityValue& routed_flow = updatePaths(curr_path, i);
+      // Update dual values
+      updateLen(curr_path, routed_flow, i);
+      ++_iter;
+    }  // end while
+
+    scaleFinalFlowsMaxFlow();
+
+    if (_heuristic) {
+      runHeurMax1();
+      runHeurMax2();
+    }
+
+    setFlowPerCommodity();
+
+    // Get dual objective
+    _dual_obj = getDualObjective();
+    // Get objective value
+    _obj = 0;
+    for (int i = 0; i < _num_commodities; ++i) {
+      for (InArcIt a(_graph, _target[i]); a != INVALID; ++a)
+        _obj += (*_flow[i])[a];
+      for (OutArcIt a(_graph, _target[i]); a != INVALID; ++a)
+        _obj -= (*_flow[i])[a];
+    }
+  }
 
-#ifdef MMF_HEURISTIC_IMPR
-      // Heuristic improvement 1: improve the solution along one of the
-      // found paths if it is possible
-      int opi = 0;
-      for (int i = 0; i < _k; ++i) {
-        for (int j = 0; j < int(paths[i].size()); ++j) {
-          Value d = std::numeric_limits<Value>::infinity();
-          for (PathArcIt a(paths[i][j].path); a != INVALID; ++a)
-            if (_cap[a] - _flow[a] < d) d = _cap[a] - _flow[a];
-          if (ftol.positive(d)) {
-            paths[i][j].value += d;
-            for (PathArcIt a(paths[i][j].path); a != INVALID; ++a)
-              _flow[a] += d;
-            ++opi;
+  // Randomised rounding - if flow is split over multiple paths, randomised
+  // rounding randomly rounds up the largest proportion of flow/cost.
+  //
+  // Only runs if flow can be routed without splitting.
+  // If costs are given, to make expensive paths less likely to occur, we use,
+  // for each path j a value[j] = flow / cost for the discrete distribution.
+  // This distribution samples path j with a probability value[j] / sum values
+  void randomisedRounding() {
+    // Only works if flow can be routed without splitting.
+    if (!_feasible_unsplit)
+      return;
+    for (int i = 0; i < _num_commodities; ++i) {
+      const int&         p_size       = static_cast<int>(_paths[i].size());
+      Value              total_flow_i = 0;
+      std::vector<Value> values(p_size, 0);
+
+      for (int j = 0; j < p_size; ++j) {
+        const Value& val  = _paths[i][j].value;
+        const Value& cost = _paths[i][j].cost;
+        total_flow_i += val;
+        values[j] = val;
+        if (_tol.positive(cost))
+          values[j] /= cost;
+      }
+      // Get random index wrt values[j] / sum of values
+      int                          chosen_p_idx;
+      std::random_device           rd;
+      std::mt19937                 gen(rd());
+      std::discrete_distribution<> d(values.begin(), values.end());
+      chosen_p_idx = d(gen);
+      // Check routing all demand through chosen_p_idx remains capacity feasible
+      bool capacity_feasible = true;
+      // cost is all demand routed through here
+      for (PathArcIt a(_paths[i][chosen_p_idx].path); a != INVALID; ++a) {
+        if (_total_flow[a] - _paths[i][chosen_p_idx].value + _demand[i] >
+            _cap[a]) {
+          capacity_feasible = false;
+          break;
+        }
+      }
+      if (capacity_feasible) {
+        for (int j = 0; j < p_size; ++j) {
+          // Update total flows and path flow
+          if (j == chosen_p_idx) {  // path j routes demand
+            for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) {
+              _total_flow[a] = _total_flow[a] - _paths[i][j].value + _demand[i];
+            }
+            _paths[i][j].value = _demand[i];
+          } else {  // path j flow is 0
+            for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) {
+              _total_flow[a] = _total_flow[a] - _paths[i][j].value;
+            }
+            _paths[i][j].value = 0;
           }
         }
       }
+    }
+  }
 
-      // Heuristic improvement 2: improve the solution along new paths
-      // if it is possible
-      int npi = 0;
-      BoolArcMap filter(_g);
-      for (ArcIt a(_g); a != INVALID; ++a)
-        filter[a] = ftol.positive(_cap[a] - _flow[a]);
-      FilterArcs<const Digraph, BoolArcMap> sub_graph(_g, filter);
-      SPath p;
-      for (int i = 0; i < _k; ++i) {
-        while (bfs(sub_graph).path(p).run(_s[i], _t[i])) {
-          Value d = std::numeric_limits<Value>::infinity();
-          for (PathArcIt a(p); a != INVALID; ++a)
-            if (_cap[a] - _flow[a] < d) d = _cap[a] - _flow[a];
-          paths[i].push_back(PathData(p, d, d));
-          for (PathArcIt a(p); a != INVALID; ++a) {
-            _flow[a] += d;
-            filter[a] = ftol.positive(_cap[a] - _flow[a]);
-          }
-          ++npi;
+  // Heuristic improvement 1 (for FLEISCHER_MAX only).
+  //
+  // Attempts to improve the solution along one of the found paths if it is
+  // possible
+  void runHeurMax1() {
+    for (int i = 0; i < _num_commodities; ++i) {
+      for (int j = 0; j < static_cast<int>(_paths[i].size()); ++j) {
+        CapacityValue d = std::numeric_limits<CapacityValue>::infinity();
+        for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a)
+          if (_cap[a] - _total_flow[a] < d)
+            d = _cap[a] - _total_flow[a];
+        if (_tol.positive(d)) {
+          _paths[i][j].value += d;
+          for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a)
+            _total_flow[a] += d;
+        }
+      }
+    }
+  }
+
+  // Heuristic improvement 2 (for FLEISCHER_MAX only).
+  //
+  // Attempts to improve the solution along new paths
+  void runHeurMax2() {
+    BoolArcMap filter(_graph);
+    for (ArcIt a(_graph); a != INVALID; ++a)
+      filter[a] = _tol.positive(_cap[a] - _total_flow[a]);
+    FilterArcs<const Digraph, BoolArcMap> sub_graph(_graph, filter);
+    Path                                  p;
+    for (int i = 0; i < _num_commodities; ++i) {
+      while (bfs(sub_graph).path(p).run(_source[i], _target[i])) {
+        CapacityValue d = std::numeric_limits<CapacityValue>::infinity();
+        for (PathArcIt a(p); a != INVALID; ++a)
+          if (_cap[a] - _total_flow[a] < d)
+            d = _cap[a] - _total_flow[a];
+        _paths[i].push_back(PathData(p, d, d, 0, 0));
+        for (PathArcIt a(p); a != INVALID; ++a) {
+          _total_flow[a] += d;
+          filter[a] = _tol.positive(_cap[a] - _total_flow[a]);
         }
       }
-#endif
-
-#ifdef MMF_DEBUG_OUTPUT
-      Value value2 = 0;
-      for (int i = 0; i < _k; ++i) {
-        for (InArcIt a(_g, _t[i]); a != INVALID; ++a) value2 += _flow[a];
-        for (OutArcIt a(_g, _t[i]); a != INVALID; ++a) value2 -= _flow[a];
-      }
-      epsilon = 1 - value2 / max_opt;
+    }
+  }
 
-  #ifdef MMF_HEURISTIC_IMPR
-      std::cout << "\nHeuristics applied. Old impr: " << opi << "  new impr: " << npi << "\n";
-  #endif
-      std::cout << "Solution value:   " << value2;
-      if (opt > 0) std::cout << "  \t" << value2 / opt * 100 << "%";
-      std::cout << std::endl << std::endl;
-#endif
+  // Update paths for current commodity and return flow routed
+  Value updatePaths(
+      const Path&          curr_path,
+      const int&           i,
+      const CapacityValue& d_i =
+          std::numeric_limits<CapacityValue>::infinity()) {
+    int h = hash(curr_path);
+    int j;  // path index
+    for (j = 0; j < static_cast<int>(_paths[i].size()); ++j) {
+      bool b = (_paths[i][j].hash == h) &&
+               (_paths[i][j].path.length() == curr_path.length());
+      if (b) {
+        for (PathArcIt a1(_paths[i][j].path), a2(curr_path);
+             b && a1 != INVALID && a2 != INVALID;
+             ++a1, ++a2)
+          b = _graph.id(a1) == _graph.id(a2);
+        if (b)
+          break;
+      }
+    }
+    const bool path_seen = !(j == static_cast<int>(_paths[i].size()));
 
-      // TODO: better return value
-      return r;
+    // Set/modify the flow value for the found path
+    Value         routed_flow = 0;
+    CapacityValue min_cap;
+    if (!path_seen) {
+      // Add new path
+      Value path_cost = 0;
+      min_cap         = std::numeric_limits<CapacityValue>::infinity();
+      for (PathArcIt a(curr_path); a != INVALID; ++a) {
+        path_cost += (*_cost[i])[a];
+        if (_cap[a] < min_cap)
+          min_cap = _cap[a];
+      }
+      routed_flow = std::min(min_cap, d_i);
+      // Save path
+      _paths[i].push_back(
+          PathData(curr_path, routed_flow, min_cap, path_cost, h));
+    } else {
+      // Update path flow
+      min_cap     = _paths[i][j].min_cap;
+      routed_flow = std::min(min_cap, d_i);
+      _paths[i][j].value += routed_flow;
     }
+    return routed_flow;
+  }
 
-    /// @}
+  // Modify the arc lengths along the current path
+  void updateLen(const Path& p, const Value& f, const int& i) {
+    Value path_cost = 0;
+    for (PathArcIt a(p); a != INVALID; ++a) {
+      const Value& c_a = (*_cost[i])[a];
+      _len[a] *= 1 + _epsilon * f / _cap[a] + c_a * _phi;
+      path_cost += c_a;
+    }
+    // Update phi if needed
+    if (_tol.positive(_phi))
+      _phi *= (1 + (_epsilon * path_cost * f) / _cost_bound);
+  }
 
-  private:
+  void updateDualObjective() {
+    _dual_obj = 0;
+    for (ArcIt a(_graph); a != INVALID; ++a)
+      _dual_obj += _cap[a] * _len[a];
+  }
 
-    int hash(const SPath& path) {
-      int h = path.length();
-      for (PathArcIt a(path); a != INVALID; ++a)
-        h = (h<<4) ^ (h>>28) ^ _g.id(a);
-      return h;
+  void setFinalCost() {
+    _total_cost = 0;
+    for (int i = 0; i < _num_commodities; ++i) {
+      for (int j = 0; j < static_cast<int>((_paths[i].size())); ++j) {
+        for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a)
+          _total_cost += (*_cost[i])[a] * (*_flow[i])[a];
+      }
     }
+  }
 
-  public:
-
-    /// \name Query Functions
-    /// The result of the algorithm can be obtained using these
-    /// functions.\n
-    /// \ref run() must be called before using them.
-
-    /// @{
+  // Scale flows in the \ref FLEISCHER_MAX_CONCURRENT.
+  // Flows are scaled by the number of iterations.
+  void scaleFinalFlowsMaxConcurrent() {
+    int scaler = _iter;
+    for (int i = 0; i < _num_commodities; ++i) {
+      for (int j = 0; j < static_cast<int>(_paths[i].size()); ++j) {
+        _paths[i][j].value /= scaler;
+        for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) {
+          _total_flow[a] += _paths[i][j].value;
+        }
+      }
+    }
+    _obj = (_iter) / scaler;
+  }
 
-    /// \brief Return the flow value on the given arc.
-    ///
-    /// This function returns the flow value on the given arc.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Value flow(const Arc& arc) const {
-      return _flow[arc];
+  // Update total flow map
+  void updateTotalFlow() {
+    for (int i = 0; i < _num_commodities; ++i) {
+      for (int j = 0; j < static_cast<int>(_paths[i].size()); ++j) {
+        CapacityValue val = _paths[i][j].value;
+        for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a)
+          _total_flow[a] += val;
+      }
     }
+  }
 
-    /// \brief Return a const reference to an arc map storing the
-    /// found flow.
-    ///
-    /// This function returns a const reference to an arc map storing
-    /// the found flow.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    const FlowMap& flowMap() const {
-      return _flow;
+  // Scale flows per commodity for the FLEISCHER_MAX case.
+  // scaler = max {total_flow_a / cap_a} such that all flows are capacity
+  // feasible
+  void scaleFinalFlowsMaxFlow() {
+    updateTotalFlow();
+    Value scaler = 0;
+    for (ArcIt a(_graph); a != INVALID; ++a) {
+      if (_total_flow[a] / _cap[a] > scaler)
+        scaler = _total_flow[a] / _cap[a];
+    }
+    if (_tol.positive(scaler)) {
+      for (ArcIt a(_graph); a != INVALID; ++a) {
+        _total_flow[a] /= scaler;
+      }
+      for (int i = 0; i < _num_commodities; ++i) {
+        for (int j = 0; j < static_cast<int>(_paths[i].size()); ++j)
+          _paths[i][j].value /= scaler;
+      }
     }
+  }
+
+  // Set flow decomposed solution for exposure
+  void setFlowPerCommodity() {
+    // Setting flow values using saved paths
+    for (int i = 0; i < _num_commodities; ++i) {
+      for (int j = 0; j < static_cast<int>(_paths[i].size()); ++j) {
+        Value val = _paths[i][j].value;
+        for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a)
+          (*_flow[i])[a] += val;
+      }
+    }
+  }
+
+  void saveFinalPaths() {
+    for (int i = 0; i < _num_commodities; ++i) {
+      const int& p_size     = static_cast<int>(_paths[i].size());
+      bool       saved_path = false;
+      Value      path_flow  = 0;
+      Value      path_cost  = std::numeric_limits<Value>::infinity();
+      Path       path;
 
-    /// \brief Return the length of the given arc (the dual value).
-    ///
-    /// This function returns the length of the given arc (the dual value).
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Value length(const Arc& arc) const {
-      return _len[arc];
+      for (int j = 0; j < p_size; ++j) {
+        const Value& v = _paths[i][j].value;
+        const Value& c = _paths[i][j].cost;
+        if (_tol.positive(v)) {
+          if (saved_path)
+            std::cout << "[WARNING] multiple paths with positive flow"
+                      << std::endl;
+          // Keep path with largest flow
+          if (path_flow < v && c < path_cost) {
+            path       = _paths[i][j].path;
+            path_flow  = v;
+            path_cost  = c;
+            saved_path = true;
+          }
+        }
+      }
+      _final_paths[i] = path;
+      _final_flows[i] = path_flow;
+      _final_costs[i] = path_cost;
     }
+  }
+
+ public:
+  const Value& flow(const int& k, const Arc& arc) const {
+    return (*_flow[k])[arc];
+  }
 
-    /// \brief Return a const reference to an arc map storing the
-    /// found lengths (the dual solution).
-    ///
-    /// This function returns a const reference to an arc map storing
-    /// the found lengths (the dual solution).
-    ///
-    /// \pre \ref run() must be called before using this function.
-    const LengthMap& lengthMap() const {
-      return _len;
-    }
+  const Value& flow(const int& k, const int& arc_id) const {
+    return (*_flow[k])[_graph.arcFromId(arc_id)];
+  }
+
+  const Value& length(const Arc& arc) const { return _len[arc]; }
+
+  const ValueMap& lengthMap() const { return _len; }
+
+  const Value& getObjective() const { return _obj; }
+
+  const Value& getDualObjective() const { return _dual_obj; }
+
+  const Value& getTotalCost() const { return _total_cost; }
 
-    /// @}
+  const ValueMap& getFlowMap(const int& k) const { return *_flow[k]; }
 
-  }; //class MaxMulticommodityFlow
+  const Path& getPath(const int& k) const { return _final_paths[k]; }
+
+  const Value& getPathCost(const int& k) const { return _final_costs[k]; }
 
-  ///@}
+  const Value& getPathFlow(const int& k) const { return _final_flows[k]; }
+};  // namespace lemon
+/// @}
 
-} //namespace lemon
+}  // namespace lemon
 
-#endif //LEMON_MULTICOMMODITY_FLOW_H
+#endif  // LEMON_MULTICOMMODITY_FLOW_H_

new file test/multicommodity_flow_test.cc

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2021
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <lemon/concepts/digraph.h>
+#include <lemon/concepts/maps.h>
+#include <lemon/lgf_reader.h>
+#include <lemon/multicommodity_flow.h>
+#include <lemon/smart_graph.h>
+
+#include <iostream>
+#include <sstream>
+
+#include "test/test_tools.h"
+
+char test_lgf[] =
+    "@nodes\n"
+    "label\n"
+    "0\n"
+    "1\n"
+    "2\n"
+    "3\n"
+    "@arcs\n"
+    "    label  cap   cost  solution_1  solution_2 solution_1_max  "
+    "solution_2_max\n"
+    "0 1   0    10.0  1.0      10.0        0.0         10.0    0.0\n"
+    "0 2   1    10.0  1.0       0.0        10.0        8.3     1.7\n"
+    "1 3   2    10.0  1.0      10.0        0.0         10.0    0.0\n"
+    "2 3   3    10.0  1000.0   0.0         0.0         8.3     0.0\n"
+    "@attributes\n"
+    "source_1 0\n"
+    "target_1 3\n"
+    "source_2 0\n"
+    "target_2 2\n";
+
+using namespace lemon;
+
+// Returns string for test-naming
+template<typename GR, typename A = typename GR::Arc>
+std::string getEdgeName(const GR& graph, const A& edge) {
+  return std::to_string(graph.id(graph.source(edge))) + "->" +
+         std::to_string(graph.id(graph.target(edge)));
+}
+
+template<class Digraph>
+void checkApproxMCF() {
+  TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+  Digraph                                           graph;
+  Node                                              s1, t1, s2, t2;
+  typedef typename Digraph::template ArcMap<double> DoubleMap;
+  DoubleMap cap(graph), cost(graph), sol1(graph), sol2(graph), sol1_max(graph),
+      sol2_max(graph);
+  Tolerance<double> tol;
+  tol.epsilon(1e-1);  // 1 decimal place
+
+  std::istringstream input(test_lgf);
+  DigraphReader<Digraph>(graph, input)
+      .arcMap("cap", cap)
+      .arcMap("cost", cost)
+      .arcMap("solution_1", sol1)  // flow in the solution for commodity 1
+      .arcMap("solution_2", sol2)  // flow in the solution for commodity 1
+      .arcMap(
+          "solution_1_max",
+          sol1_max)  // flow in the solution for commodity (max-flow) 1
+      .arcMap(
+          "solution_2_max",
+          sol2_max)          // flow in the solution for commodity (max-flow) 1
+      .node("source_1", s1)  // source node commodity 1
+      .node("target_1", t1)  // target node commodity 1
+      .node("source_2", s2)  // source node commodity 2
+      .node("target_2", t2)  // target node commodity 2
+      .run();
+
+  // Max
+  {
+    ApproxMCF<Digraph> max(graph, cap);
+    max.addCommodity(s1, t1).addCommodity(s2, t2).run();
+    for (ArcIt e(graph); e != INVALID; ++e) {
+      const std::string test_name =
+          "[FLEISCHER_MAX] flow on " + getEdgeName<Digraph>(graph, e);
+      // Check flow for commodity 1 (index 0) to 1 decimal place
+      check(
+          !tol.different(max.flow(0, e), sol1_max[e]),
+          test_name + " commodity 1");
+      // Check flow for commodity 2 (index 1) to 1 decimal place
+      check(
+          !tol.different(max.flow(1, e), sol2_max[e]),
+          test_name + " commodity 2");
+    }
+  }
+
+  // Max Concurrent
+  {
+    ApproxMCF<Digraph> max(graph, cap);
+    max.addCommodity(s1, t1, 10.0)
+        .addCommodity(s2, t2, 10.0)
+        .run(ApproxMCF<Digraph>::FLEISCHER_MAX_CONCURRENT);
+    for (ArcIt e(graph); e != INVALID; ++e) {
+      const std::string test_name = "[FLEISCHER_MAX_CONCURRENT] flow on " +
+                                    getEdgeName<Digraph>(graph, e);
+      // Check commodity 1 (index 0)
+      check(max.flow(0, e) == sol1[e], test_name + " commodity 1");
+      // Check commodity 2 (index 1)
+      check(max.flow(1, e) == sol2[e], test_name + " commodity 2");
+    }
+  }
+
+  // Min cost
+  {
+    ApproxMCF<Digraph> max(graph, cap);
+    max.addCommodity(s1, t1, 10.0, cost)
+        .addCommodity(s2, t2, 10.0, cost)
+        .run(ApproxMCF<Digraph>::BINARY_SEARCH_MIN_COST);
+    for (ArcIt e(graph); e != INVALID; ++e) {
+      const std::string test_name =
+          "[BINARY_SEARCH_MIN_COST] flow on " + getEdgeName<Digraph>(graph, e);
+      // Check commodity 1 (index 0)
+      check(max.flow(0, e) == sol1[e], test_name + " commodity 1");
+      // Check commodity 2 (index 1)
+      check(max.flow(1, e) == sol2[e], test_name + " commodity 2");
+    }
+    // Check total cost
+    check(max.getTotalCost() == 30.0, "min-cost failed");
+  }
+
+  // Max Concurrent - 3 commodities
+  {
+    ApproxMCF<Digraph> max(graph, cap);
+    max.addCommodity(s1, t1, 5.0, cost)
+        .addCommodity(s1, t2, 10.0, cost)
+        .addCommodity(s1, t1, 5.0, cost)
+        .run(ApproxMCF<Digraph>::FLEISCHER_MAX_CONCURRENT);
+
+    const DoubleMap& flow = max.getFlowMap(0);
+    check(
+        max.getTotalCost() == 30.0,
+        "3 commodities FLEISCHER_MAX_CONCURRENT failed");
+  }
+}
+
+template<class Digraph>
+void checkApproxMCFDisconnected() {
+  TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+  Digraph                                           graph;
+  typedef typename Digraph::template ArcMap<double> DoubleMap;
+  DoubleMap cap(graph), cost(graph), sol1(graph), sol2(graph);
+
+  Node n0 = graph.addNode();
+  Node n1 = graph.addNode();
+  Node n2 = graph.addNode();
+  Node n3 = graph.addNode();
+
+  Arc a;
+  a       = graph.addArc(n0, n1);
+  cap[a]  = 10.0;
+  cost[a] = 1.0;
+  sol1[a] = 10.0;
+  sol2[a] = 0.0;
+
+  a       = graph.addArc(n2, n3);
+  cap[a]  = 10.0;
+  cost[a] = 1.0;
+  sol1[a] = 0.0;
+  sol2[a] = 10.0;
+
+  // Max
+  {
+    ApproxMCF<Digraph> max(graph, cap);
+    max.addCommodity(n0, n1).addCommodity(n2, n3).run();
+    for (ArcIt e(graph); e != INVALID; ++e) {
+      const std::string test_name =
+          "[FLEISCHER_MAX] flow on " + getEdgeName<Digraph>(graph, e);
+      // Check commodity 1
+      check(max.flow(0, e) == sol1[e], test_name + " commodity 1");
+      // Check commodity 2
+      check(max.flow(1, e) == sol2[e], test_name + " commodity 2");
+    }
+  }
+  // Max Concurrent
+  {
+    ApproxMCF<Digraph> max(graph, cap);
+    max.addCommodity(n0, n1, 10.0, cost)
+        .addCommodity(n2, n3, 10.0, cost)
+        .run(ApproxMCF<Digraph>::FLEISCHER_MAX_CONCURRENT);
+    for (ArcIt e(graph); e != INVALID; ++e) {
+      const std::string test_name = "[FLEISCHER_MAX_CONCURRENT] flow on " +
+                                    getEdgeName<Digraph>(graph, e);
+      // Check commodity 1
+      check(max.flow(0, e) == sol1[e], test_name + " commodity 1");
+      // Check commodity 2
+      check(max.flow(1, e) == sol2[e], test_name + " commodity 2");
+    }
+  }
+}
+
+int main() {
+  checkApproxMCF<SmartDigraph>();
+  checkApproxMCF<ListDigraph>();
+  checkApproxMCFDisconnected<SmartDigraph>();
+  checkApproxMCFDisconnected<ListDigraph>();
+}