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 <limits>
#include <queue>
#include <functional>
#include <tpl_netcost.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,
                                      typename Net::Node* from, bool forward)
{
  (void)from;
  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*>;
  std::priority_queue<PQ_Entry, std::vector<PQ_Entry>, std::greater<PQ_Entry>> pq;
 
  pq.push({0, source});
 
  while (not pq.empty())
    {
      auto [dist, u] = pq.top();
      pq.pop();
 
      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, u, forward);
 
          // Relaxation
          Flow_Type new_dist = ssp_info<Net>(u).distance + 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.push({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 size_t n = net.vsize();
  
  // 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 = 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 = d - arc->cost;  // negative cost for backward
                  if (new_d < ssp_info<Net>(src).potential)
                    {
                      ssp_info<Net>(src).potential = new_d;
                      changed = true;
                    }
                }
            }
        }
      if (!changed) break;
    }
  
  // 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;  // No negative cycle check for simplicity
}
 
template <class 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 using std::vector to allow taking address of elements
  std::vector<SSP_Node_Info<Flow_Type>> node_info(net.vsize());
  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];
 
  // Initialize potentials using Bellman-Ford to handle possible negative costs
  ssp_init_potentials(net, source);
 
  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);
          bool forward = ssp_info<Net>(v).forward;
 
          if (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
      Flow_Type max_potential = std::numeric_limits<Flow_Type>::lowest();
      for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
        {
          auto p = it.get_curr();
          if (ssp_info<Net>(p).distance < std::numeric_limits<Flow_Type>::max())
            {
            ssp_info<Net>(p).potential += ssp_info<Net>(p).distance;
              max_potential = std::max(max_potential, ssp_info<Net>(p).potential);
            }
        }
      
      // For unreachable nodes, set potential to max to avoid negative reduced costs
      // when backward arcs to them become available
      for (Node_Iterator<Net> it(net); it.has_curr(); it.next_ne())
        {
          auto p = it.get_curr();
          if (ssp_info<Net>(p).distance >= std::numeric_limits<Flow_Type>::max())
            ssp_info<Net>(p).potential = max_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;
    }
 
  // Create super-source and super-sink
  Node* super_source = net.insert_node();
  Node* super_sink = net.insert_node();
 
  // 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);
 
  // Clean up (remove super nodes)
  // Note: In real implementation, would need to properly remove these
 
  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
  for (const auto& row : costs)
    if (row.size() != n)
      {
        result.feasible = false;
        return result;
      }
 
  Net net;
 
  // Create source and sink
  Node* source = net.insert_node();
  Node* sink = net.insert_node();
 
  // Create worker nodes
  std::vector<Node*> workers(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
  std::vector<Node*> tasks(n);
  for (size_t j = 0; j < n; ++j)
    {
      tasks[j] = net.insert_node();
      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();
              auto tgt = net.get_tgt_node(arc);
 
              // Find which task this is
              for (size_t j = 0; j < n; ++j)
                if (tasks[j] == tgt and arc->flow > 0)
                  {
                    result.assignments.append(std::make_pair(i, j));
                    break;
                  }
            }
        }
    }
 
  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
 
  if (m == 0 or n == 0)
    {
      result.feasible = true;
      return result;
    }
 
  // Check supply-demand balance
  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 (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]) + 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 = (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