tpl_mincost.H

Go to the documentation of this file.

 
/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
#ifndef TPL_MINCOST_H
#define TPL_MINCOST_H
 
# include <ah-graph-concepts.H>
 
#include <limits>
#include <queue>
#include <functional>
#include <tpl_netcost.H>
#include <tpl_array.H>
#include <tpl_dynBinHeap.H>
#include <ah-errors.H>
 
namespace Aleph {
 
//==============================================================================
// SUCCESSIVE SHORTEST PATHS (SSP) ALGORITHM
//==============================================================================
 
template <typename Flow_Type>
struct SSP_Node_Info
{
  Flow_Type potential{0};      
  Flow_Type distance{0};       
  bool in_tree{false};         
  void *parent_arc{nullptr};   
  void *parent_node{nullptr};  
  bool forward{true};          
 
  void reset()
  {
    distance = std::numeric_limits<Flow_Type>::max();
    in_tree = false;
    parent_arc = nullptr;
    parent_node = nullptr;
  }
};
 
template <class Net>
inline SSP_Node_Info<typename Net::Flow_Type> &ssp_info(typename Net::Node *p) noexcept
{
  return *static_cast<SSP_Node_Info<typename Net::Flow_Type> *>(NODE_COOKIE(p));
}
 
template <class Net>
typename Net::Flow_Type reduced_cost(const Net &net, typename Net::Arc *arc, bool forward)
{
  using Flow_Type = typename Net::Flow_Type;
 
  auto src = net.get_src_node(arc);
  auto tgt = net.get_tgt_node(arc);
 
  Flow_Type pi_src = ssp_info<Net>(src).potential;
  Flow_Type pi_tgt = ssp_info<Net>(tgt).potential;
 
  // Forward arc: src -> tgt with cost arc->cost
  // Reduced cost = cost + π(src) - π(tgt)
  if (forward)
    return arc->cost + pi_src - pi_tgt;
 
  // Backward arc: tgt -> src with cost -arc->cost
  // Reduced cost = -cost + π(tgt) - π(src)
  return -arc->cost + pi_tgt - pi_src;
}
 
template <class Net>
bool ssp_shortest_path(Net &net, typename Net::Node *source, typename Net::Node *sink)
{
  using Node = typename Net::Node;
  using Arc = typename Net::Arc;
  using Flow_Type = typename Net::Flow_Type;
 
  const Flow_Type INF = std::numeric_limits<Flow_Type>::max();
 
  // Reset distances
  for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    ssp_info<Net>(it.get_curr()).reset();
 
  ssp_info<Net>(source).distance = 0;
 
  // Priority queue: (distance, node)
  using PQ_Entry = std::pair<Flow_Type, Node *>;
  DynBinHeap<PQ_Entry> pq;
 
  pq.insert({0, source});
 
  while (not pq.is_empty())
    {
      auto [dist, u] = pq.getMin();
 
      if (dist > ssp_info<Net>(u).distance)
        continue;  // Outdated entry
 
      // Explore all adjacent arcs
      for (typename Net::Node_Arc_Iterator it(u); it.has_curr(); it.next_ne())
        {
          Arc *arc = it.get_curr();
          Node *v = net.get_connected_node(arc, u);
 
          // Determine direction and residual capacity
          bool forward = (net.get_src_node(arc) == u);
          Flow_Type residual = forward ? (arc->cap - arc->flow) : arc->flow;
 
          if (residual <= Flow_Type{0})
            continue;
 
          // Compute reduced cost
          Flow_Type rc = reduced_cost<Net>(net, arc, forward);
 
          // Relaxation (saturate to avoid integer overflow)
          const Flow_Type INF_D = std::numeric_limits<Flow_Type>::max();
          const Flow_Type ud = ssp_info<Net>(u).distance;
          Flow_Type new_dist;
          if constexpr (std::is_integral_v<Flow_Type>)
            new_dist = (rc > 0 and ud > INF_D - rc) ? INF_D : ud + rc;
          else
            new_dist = ud + rc;
          if (new_dist < ssp_info<Net>(v).distance)
            {
              ssp_info<Net>(v).distance = new_dist;
              ssp_info<Net>(v).parent_arc = arc;
              ssp_info<Net>(v).parent_node = u;
              ssp_info<Net>(v).forward = forward;
              pq.insert({new_dist, v});
            }
        }
    }
 
  return ssp_info<Net>(sink).distance < INF;
}
 
template <class Net>
std::pair<typename Net::Flow_Type, typename Net::Flow_Type> successive_shortest_paths(Net &net);
 
template <class Net>
bool ssp_init_potentials(Net &net, typename Net::Node *source)
{
  using Flow_Type = typename Net::Flow_Type;
 
  const Flow_Type INF    = std::numeric_limits<Flow_Type>::max();
  const Flow_Type LOWEST = std::numeric_limits<Flow_Type>::lowest();
  const size_t n = net.vsize();
 
  // Overflow-safe addition/subtraction clamped to [LOWEST, INF].
  auto sat_add = [INF, LOWEST](Flow_Type a, Flow_Type b) noexcept -> Flow_Type {
    if constexpr (std::is_integral_v<Flow_Type>) {
      if (b > 0 and a > INF - b)    return INF;
      if (b < 0 and a < LOWEST - b) return LOWEST;
    }
    return a + b;
  };
  auto sat_sub = [INF, LOWEST](Flow_Type a, Flow_Type b) noexcept -> Flow_Type {
    if constexpr (std::is_integral_v<Flow_Type>) {
      if (b < 0 and a > INF + b)    return INF;
      if (b > 0 and a < LOWEST + b) return LOWEST;
    }
    return a - b;
  };
 
  // Initialize distances
  for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    ssp_info<Net>(it.get_curr()).potential = INF;
  ssp_info<Net>(source).potential = 0;
 
  // Bellman-Ford: relax all edges n-1 times
  for (size_t i = 0; i < n - 1; ++i)
    {
      bool changed = false;
      for (Arc_Iterator<Net> it(net); it.has_curr(); it.next_ne())
        {
          auto arc = it.get_curr();
          auto src = net.get_src_node(arc);
          auto tgt = net.get_tgt_node(arc);
 
          // Forward arc (if has residual capacity)
          if (arc->cap > arc->flow)
            {
              Flow_Type d = ssp_info<Net>(src).potential;
              if (d < INF)
                {
                  Flow_Type new_d = sat_add(d, arc->cost);
                  if (new_d < ssp_info<Net>(tgt).potential)
                    {
                      ssp_info<Net>(tgt).potential = new_d;
                      changed = true;
                    }
                }
            }
 
          // Backward arc (if has flow to cancel)
          if (arc->flow > 0)
            {
              Flow_Type d = ssp_info<Net>(tgt).potential;
              if (d < INF)
                {
                  Flow_Type new_d = sat_sub(d, arc->cost);  // backward: cost is negated
                  if (new_d < ssp_info<Net>(src).potential)
                    {
                      ssp_info<Net>(src).potential = new_d;
                      changed = true;
                    }
                }
            }
        }
      if (not changed)
        break;
    }
 
  // Negative-cycle detection: one extra relaxation pass.
  // If any edge can still be relaxed, a negative cycle is reachable.
  for (Arc_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    {
      auto arc = it.get_curr();
      auto src = net.get_src_node(arc);
      auto tgt = net.get_tgt_node(arc);
 
      if (arc->cap > arc->flow)
        {
          Flow_Type d = ssp_info<Net>(src).potential;
          if (d < INF and sat_add(d, arc->cost) < ssp_info<Net>(tgt).potential)
            return false;
        }
 
      if (arc->flow > 0)
        {
          Flow_Type d = ssp_info<Net>(tgt).potential;
          if (d < INF and sat_sub(d, arc->cost) < ssp_info<Net>(src).potential)
            return false;
        }
    }
 
  // Set unreachable nodes to 0 (they won't affect anything)
  for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    if (ssp_info<Net>(it.get_curr()).potential >= INF)
      ssp_info<Net>(it.get_curr()).potential = 0;
 
  return true;
}
 
template <FlowNetwork Net>
std::pair<typename Net::Flow_Type, typename Net::Flow_Type> successive_shortest_paths(Net &net)
{
  using Node = typename Net::Node;
  using Arc = typename Net::Arc;
  using Flow_Type = typename Net::Flow_Type;
 
  ah_domain_error_if(not(net.is_single_source() and net.is_single_sink()))
    << "Network must have single source and single sink";
 
  Node *source = net.get_source();
  Node *sink = net.get_sink();
 
  // Allocate node info; fixed size, no reallocation, so element addresses are stable
  Array<SSP_Node_Info<Flow_Type>> node_info(net.vsize(), SSP_Node_Info<Flow_Type>{});
  size_t idx = 0;
  for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne(), ++idx)
    NODE_COOKIE(it.get_curr()) = &node_info[idx];
 
  // RAII: clear node cookies on every exit path so callers never see
  // dangling pointers into the local node_info vector.
  struct CookieClearer
  {
    Net &net_;
    ~CookieClearer()
    {
      for (Node_Iterator<Net> it(net_); it.has_curr(); it.next_ne())
        NODE_COOKIE(it.get_curr()) = nullptr;
    }
  } cookie_clearer{net};
 
  // Initialize potentials using Bellman-Ford to handle possible negative costs
  ah_domain_error_if(not ssp_init_potentials(net, source))
    << "successive_shortest_paths: negative-cost cycle detected in residual graph";
 
  Flow_Type total_flow{0};
  Flow_Type total_cost{0};
 
  // Main loop: find shortest augmenting paths
  while (ssp_shortest_path(net, source, sink))
    {
      // Find minimum residual capacity on path
      Flow_Type path_flow = std::numeric_limits<Flow_Type>::max();
      for (Node *v = sink; v != source;)
        {
          Arc *arc = static_cast<Arc *>(ssp_info<Net>(v).parent_arc);
          bool forward = ssp_info<Net>(v).forward;
 
          Flow_Type residual = forward ? (arc->cap - arc->flow) : arc->flow;
          path_flow = std::min(path_flow, residual);
 
          v = static_cast<Node *>(ssp_info<Net>(v).parent_node);
        }
 
      // Augment flow and compute cost
      for (Node *v = sink; v != source;)
        {
          Arc *arc = static_cast<Arc *>(ssp_info<Net>(v).parent_arc);
 
          if (bool forward = ssp_info<Net>(v).forward)
            {
              arc->flow += path_flow;
              total_cost += path_flow * arc->cost;
            }
          else
            {
              arc->flow -= path_flow;
              total_cost -= path_flow * arc->cost;
            }
 
          v = static_cast<Node *>(ssp_info<Net>(v).parent_node);
        }
 
      total_flow += path_flow;
 
      // Update potentials for reachable nodes; track max updated potential
      // and max finite distance for the unreachable-node shift below.
      Flow_Type max_potential = std::numeric_limits<Flow_Type>::lowest();
      Flow_Type max_distance  = Flow_Type{0};
      for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
        if (auto p = it.get_curr(); ssp_info<Net>(p).distance < std::numeric_limits<Flow_Type>::max())
          {
            max_distance  = std::max(max_distance, ssp_info<Net>(p).distance);
            ssp_info<Net>(p).potential += ssp_info<Net>(p).distance;
            max_potential = std::max(max_potential, ssp_info<Net>(p).potential);
          }
 
      // Unreachable nodes are assigned potential = max_potential + max_distance.
      // This ensures that for any residual arc (u→v) with u unreachable and v
      // reachable, the reduced cost remains non-negative: the shift by
      // max_distance dominates the largest finite d(v) seen in Dijkstra.
      // Saturate to avoid integer overflow when both terms are large.
      Flow_Type unreachable_potential;
      if constexpr (std::is_integral_v<Flow_Type>)
        {
          const Flow_Type INF_P = std::numeric_limits<Flow_Type>::max();
          unreachable_potential =
            (max_distance > 0 and max_potential > INF_P - max_distance)
            ? INF_P
            : max_potential + max_distance;
        }
      else
        unreachable_potential = max_potential + max_distance;
 
      for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
        if (auto p = it.get_curr();
            ssp_info<Net>(p).distance >= std::numeric_limits<Flow_Type>::max())
          ssp_info<Net>(p).potential = unreachable_potential;
    }
 
  return {total_flow, total_cost};
}
 
template <class Net>
struct Successive_Shortest_Paths
{
  std::pair<typename Net::Flow_Type, typename Net::Flow_Type> operator()(Net &net) const
  {
    return successive_shortest_paths(net);
  }
};
 
//==============================================================================
// MINIMUM COST FLOW WITH DEMAND (TRANSSHIPMENT)
//==============================================================================
 
template <typename Flow_Type>
struct Transshipment_Node_Info
{
  Flow_Type demand{0};  
};
 
template <typename Flow_Type>
struct TransshipmentResult
{
  bool feasible{false};     
  Flow_Type total_cost{0};  
  Flow_Type total_flow{0};  
};
 
template <class Net, class GetDemand>
TransshipmentResult<typename Net::Flow_Type> solve_transshipment(Net &net, GetDemand get_demand)
{
  using Node = typename Net::Node;
  using Flow_Type = typename Net::Flow_Type;
 
  TransshipmentResult<Flow_Type> result;
 
  // Calculate total supply and demand
  Flow_Type total_supply{0};
  Flow_Type total_demand{0};
 
  for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    {
      Flow_Type d = get_demand(it.get_curr());
      if (d > 0)
        total_supply += d;
      else
        total_demand -= d;
    }
 
  // Check balance
  if (total_supply != total_demand)
    {
      result.feasible = false;
      return result;
    }
 
  // Zero total flow: trivially feasible with no flow needed.
  if (total_supply == Flow_Type{0})
    {
      result.feasible = true;
      return result;
    }
 
  // RAII guard: removes super nodes (and all incident arcs) on every exit path.
  // src_/snk_ are set incrementally so a throw during the second insert_node()
  // still cleans up the first node.
  struct SuperNodeGuard
  {
    Net &net_;
    Node *src_{nullptr};
    Node *snk_{nullptr};
    ~SuperNodeGuard()
    {
      if (snk_)
        net_.remove_node(snk_);
      if (src_)
        net_.remove_node(src_);
    }
  } guard{net};
 
  // Create super-source and super-sink
  guard.src_ = net.insert_node();
  Node *super_source = guard.src_;
  guard.snk_ = net.insert_node();
  Node *super_sink = guard.snk_;
 
  // Connect super-source to supply nodes
  for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    {
      Node *p = it.get_curr();
      if (p == super_source or p == super_sink)
        continue;
 
      Flow_Type d = get_demand(p);
      if (d > 0)
        net.insert_arc(super_source, p, d, Flow_Type{0});  // zero cost
      else if (d < 0)
        net.insert_arc(p, super_sink, -d, Flow_Type{0});  // zero cost
    }
 
  // Source and sink are detected automatically by network topology
 
  // Solve min-cost max-flow
  auto [flow, cost] = successive_shortest_paths(net);
 
  result.total_flow = flow;
  result.total_cost = cost;
  result.feasible = (flow == total_supply);
 
  return result;
}
 
//==============================================================================
// ASSIGNMENT PROBLEM
//==============================================================================
 
template <typename Cost_Type>
struct AssignmentResult
{
  bool feasible{false};
  Cost_Type total_cost{0};
  DynList<std::pair<size_t, size_t>> assignments;  
};
 
template <typename Cost_Type>
AssignmentResult<Cost_Type> solve_assignment(const std::vector<std::vector<Cost_Type>> &costs)
{
  using Net = Net_Cost_Graph<Net_Cost_Node<Empty_Class>, Net_Cost_Arc<Empty_Class, Cost_Type>>;
  using Node = typename Net::Node;
 
  AssignmentResult<Cost_Type> result;
 
  size_t n = costs.size();
  if (n == 0)
    {
      result.feasible = true;
      return result;
    }
 
  // Verify square matrix: a non-square costs matrix is malformed input.
  for (const auto &row : costs)
    ah_invalid_argument_if(row.size() != n)
      << "solve_assignment: costs matrix is not square (" << n << " x " << row.size() << ")";
 
  Net net;
 
  // Create source and sink
  Node *source = net.insert_node();
  Node *sink = net.insert_node();
 
  // Create worker nodes
  auto workers = Array<Node *>::create(n);
  for (size_t i = 0; i < n; ++i)
    {
      workers[i] = net.insert_node();
      net.insert_arc(source, workers[i], Cost_Type{1}, Cost_Type{0});
    }
 
  // Create task nodes
  auto tasks = Array<Node *>::create(n);
  DynMapTree<Node *, size_t> taskIndex;  // Reverse lookup: task node -> index
  for (size_t j = 0; j < n; ++j)
    {
      tasks[j] = net.insert_node();
      taskIndex.insert(tasks[j], j);
      net.insert_arc(tasks[j], sink, Cost_Type{1}, Cost_Type{0});
    }
 
  // Create worker-task arcs with costs
  for (size_t i = 0; i < n; ++i)
    for (size_t j = 0; j < n; ++j)
      net.insert_arc(workers[i], tasks[j], Cost_Type{1}, costs[i][j]);
 
  // Source and sink are detected automatically by network topology
 
  // Solve
  auto [flow, cost] = successive_shortest_paths(net);
 
  result.feasible = (static_cast<size_t>(flow) == n);
  result.total_cost = cost;
 
  // Extract assignments from flow
  if (result.feasible)
    for (size_t i = 0; i < n; ++i)
      for (typename Net::Node_Arc_Iterator it(workers[i]); it.has_curr(); it.next_ne())
        {
          auto arc = it.get_curr();
          if (arc->flow <= 0)
            continue;
          auto tgt = net.get_tgt_node(arc);
 
          // O(log n) lookup instead of O(n) linear search
          if (taskIndex.contains(tgt))
            {
              size_t j = taskIndex.find(tgt);
              result.assignments.append(std::make_pair(i, j));
            }
        }
 
  return result;
}
 
//==============================================================================
// TRANSPORTATION PROBLEM
//==============================================================================
 
template <typename Cost_Type>
struct TransportationResult
{
  bool feasible{false};
  Cost_Type total_cost{0};
  std::vector<std::vector<Cost_Type>> shipments;  
};
 
template <typename Cost_Type>
TransportationResult<Cost_Type> solve_transportation(const std::vector<Cost_Type> &supplies,
                                                     const std::vector<Cost_Type> &demands,
                                                     const std::vector<std::vector<Cost_Type>> &costs)
{
  using Net = Net_Cost_Graph<Net_Cost_Node<Empty_Class>, Net_Cost_Arc<Empty_Class, Cost_Type>>;
  using Node = typename Net::Node;
  using Arc = typename Net::Arc;
 
  TransportationResult<Cost_Type> result;
 
  size_t m = supplies.size();  // sources
  size_t n = demands.size();   // destinations
 
  // Validate costs matrix dimensions
  ah_invalid_argument_if(costs.size() != m)
    << "solve_transportation: costs has " << costs.size() << " rows but " << m << " supply points";
  for (size_t i = 0; i < m; ++i)
    ah_invalid_argument_if(costs[i].size() != n)
      << "solve_transportation: costs[" << i << "] has " << costs[i].size() << " columns but " << n
      << " demand points";
 
  // Validate non-negative supplies and demands
  for (size_t i = 0; i < m; ++i)
    ah_invalid_argument_if(supplies[i] < Cost_Type{0})
      << "solve_transportation: supplies[" << i << "] is negative";
  for (size_t j = 0; j < n; ++j)
    ah_invalid_argument_if(demands[j] < Cost_Type{0})
      << "solve_transportation: demands[" << j << "] is negative";
 
  // Compute totals first — needed both for the early-exit guard and the
  // balance check that follows.
  Cost_Type total_supply = Cost_Type{0};
  Cost_Type total_demand = Cost_Type{0};
  for (auto s : supplies)
    total_supply += s;
  for (auto d : demands)
    total_demand += d;
 
  // If either side is empty no graph can be built.  Return immediately based
  // on whether the totals are balanced (feasible iff both are zero).
  if (m == 0 or n == 0)
    {
      result.feasible = (total_supply == total_demand);
      return result;
    }
 
  // Check supply-demand balance for the non-trivial case.
  if (total_supply != total_demand)
    {
      result.feasible = false;
      return result;
    }
 
  Net net;
 
  // Create source and sink
  Node *source = net.insert_node();
  Node *sink = net.insert_node();
 
  // Create supply nodes
  std::vector<Node *> supply_nodes(m);
  for (size_t i = 0; i < m; ++i)
    {
      supply_nodes[i] = net.insert_node();
      net.insert_arc(source, supply_nodes[i], supplies[i], Cost_Type{0});
    }
 
  // Create demand nodes
  std::vector<Node *> demand_nodes(n);
  for (size_t j = 0; j < n; ++j)
    {
      demand_nodes[j] = net.insert_node();
      net.insert_arc(demand_nodes[j], sink, demands[j], Cost_Type{0});
    }
 
  // Create arcs between supplies and demands
  // Store arcs for extracting solution
  std::vector<std::vector<Arc *>> arcs(m, std::vector<Arc *>(n));
  for (size_t i = 0; i < m; ++i)
    for (size_t j = 0; j < n; ++j)
      arcs[i][j] = net.insert_arc(supply_nodes[i],
                                  demand_nodes[j],
                                  std::min(supplies[i], demands[j]),
                                  costs[i][j]);
 
  // Source and sink are detected automatically by network topology
 
  // Solve
  auto [flow, cost] = successive_shortest_paths(net);
 
  result.feasible = (flow == total_supply);
  result.total_cost = cost;
 
  // Extract shipments
  result.shipments.resize(m, std::vector<Cost_Type>(n, Cost_Type{0}));
  for (size_t i = 0; i < m; ++i)
    for (size_t j = 0; j < n; ++j)
      result.shipments[i][j] = arcs[i][j]->flow;
 
  return result;
}
 
//==============================================================================
// MIN-COST FLOW STATISTICS
//==============================================================================
 
template <typename Flow_Type>
struct MinCostFlowStats
{
  Flow_Type max_flow{0};
  Flow_Type min_cost{0};
  size_t iterations{0};
  size_t augmentations{0};
  double elapsed_ms{0};
 
  void print() const
  {
    std::cout << "=== Min-Cost Flow Statistics ===\n"
              << "Max flow: " << max_flow << "\n"
              << "Min cost: " << min_cost << "\n"
              << "Iterations: " << iterations << "\n"
              << "Augmentations: " << augmentations << "\n"
              << "Time: " << elapsed_ms << " ms\n";
  }
};
 
//==============================================================================
// WARM START FOR NETWORK SIMPLEX
//==============================================================================
 
template <class Net>
DynMapTree<typename Net::Arc *, typename Net::Flow_Type> save_flow_solution(const Net &net)
{
  DynMapTree<typename Net::Arc *, typename Net::Flow_Type> solution;
 
  for (Arc_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    {
      auto arc = it.get_curr();
      solution[arc] = arc->flow;
    }
 
  return solution;
}
 
template <class Net>
void restore_flow_solution(Net &net,
                           const DynMapTree<typename Net::Arc *, typename Net::Flow_Type> &solution)
{
  for (Arc_Iterator<Net> it(net); it.has_curr(); it.next_ne())
    {
      auto arc = it.get_curr();
      if (solution.has(arc))
        arc->flow = solution[arc];
    }
}
 
}  // namespace Aleph
 
#endif  // TPL_MINCOST_H