Algoteka
Log in to submit your own samples!

Primal-Dual Min-Cost Max-Flow Algorithm

Problem by oml1111
# Tech tags Title Creator Created date
1 0
2022-10-15 20:39
View all samples for this language (1 verified and 0 unverified)

O(F n^2m) Readable Solution | C++ |

By oml1111 |
0 likes

An \(O(F n^2m)\) implementation of a variation of the Primal-Dual algorithm (\(F\) is the maximum possible flow, \(n\) is the number of vertices, \(m\) is the number of edges). Usually performs much better in practice. Performs \(O(F)\) augmenting loops in which it uses a modification of the Bellman-Ford algorithm to compute the potential for each node and uses a modification of the Dinic's algorithm to saturate paths with cost equal to the potential of the sink node.

Caveat: The algorithm will fail if there exists a negative-cost cycle

Code

#include<vector>
#include<deque>
#include<limits>
#include<functional>
#include<algorithm>
#include<iostream>


struct PrimalDualFlowNetwork {
    struct Arc {
        int start_;
        int end_;
        int flow_;
        int capacity_;
        int cost_;
        
        int get_dest(int from) {  // Return the destination node if we were to traverse the arc from node `from`
            if(from == start_) {
                return end_;
            } else {
                return start_;
            }
        }
        
        void add_flow(int from, int to_add) {  // Adds flow from originating vertex `from`
            if(from == start_) {
                flow_ += to_add;
            } else {
                flow_ -= to_add;
            }
        }
        
        int get_capacity(int from) {  // Gets the capacity of the edge if the originating vertex is `from`
            if(from == start_) {
                return capacity_ - flow_;
            } else {
                return flow_;
            }
        }
        
        int get_cost_from(int from) {
            if(from == start_) {
                return cost_;
            } else {
                return -cost_;
            }
        }
    };
    
    struct Node {
        int index_;
        std::vector<Arc*> connected_arcs_;
    };
    
    std::vector<Node> nodes_;
    std::deque<Arc> arcs_;
    
    void add_node() {
        nodes_.push_back({int(nodes_.size()), {}});
    }
    
    void add_arc(int start, int end, int flow, int capacity, int cost) {
        arcs_.push_back({start, end, flow, capacity, cost});
        nodes_[start].connected_arcs_.push_back(&arcs_.back());
        nodes_[end].connected_arcs_.push_back(&arcs_.back());
    }
    
    // Dinic that fills the zero-cost path
    int dinic_with_potentials(int source_i, int sink_i, std::vector<int> &potentials) {
        int result = 0;
        
        auto get_modified_cost = [&](int node_i, Arc* arc) {
            if(node_i == arc->start_)
                return potentials[arc->start_] - potentials[arc->end_] + arc->cost_;
            return potentials[arc->end_] - potentials[arc->start_] - arc->cost_;
        };
        
        while(1) {
            // First divide into layers with a breadth-first search
            std::vector<int> level(nodes_.size(), -1);
            {
                std::deque<std::tuple<int, int> > queue(1, std::make_tuple(0, source_i) );
                level[source_i] = 0;
                
                while(queue.size() > 0) {
                    int cur_level;
                    int cur_node_i;
                    std::tie(cur_level, cur_node_i) = queue.front();  // Extract the values from queue
                    queue.pop_front();
                    
                    for(Arc* arc : nodes_[cur_node_i].connected_arcs_) {
                        if(arc->get_capacity(cur_node_i) > 0 && get_modified_cost(cur_node_i, arc) == 0) {
                            int next_i = arc->get_dest(cur_node_i);
                            if(level[next_i] == -1) {
                                level[next_i] = cur_level+1;
                                queue.push_back(std::make_tuple(cur_level+1, next_i));
                            }
                        }
                    }
                }
                
                if(level[sink_i] == -1) {
                    return result;  // Path to sink no longer exists, so return
                }
            }
            
            // Now perform DFS to push through the blocking flow
            // Since due to leveling we never push back flow, we can speed up DFS by remembering arc iteration indices
            std::vector<int> arc_index(nodes_.size(), 0);
            
            std::function<int (int, int)> dfs = [&](int cur_i, int capacity) {
                if(cur_i == sink_i) {
                    return capacity;
                }
                
                for(int& i = arc_index[cur_i]; i < (int)nodes_[cur_i].connected_arcs_.size(); i++) {
                    Arc* arc = nodes_[cur_i].connected_arcs_[i];
                    int next_i = arc->get_dest(cur_i);
                    if(arc->get_capacity(cur_i) > 0 && get_modified_cost(cur_i, arc) == 0 && level[cur_i] < level[next_i]) {
                        int pushed_flow = dfs(next_i, std::min(capacity, arc->get_capacity(cur_i)));
                        if(pushed_flow > 0) {
                            arc->add_flow(cur_i, pushed_flow);
                            return pushed_flow;
                        }
                    }
                }
                return 0;
            };
            int pushed_flow;
            do {
                pushed_flow = dfs(source_i, std::numeric_limits<int>::max());
                result += pushed_flow;
            } while(pushed_flow > 0);
        }
    }
    
    // Primal-Dual min-cost max-flow algorithm
    // If there is an negative cost cycle initially, then it goes into infinite loop
    int min_cost_max_flow(int source_i, int sink_i) {
        int result = 0;
        
        while(1) {
            // First calculate the potentials with Bellman–Ford derivative
            // It starts from a single vertex and is optimized to only operate on active vertices on each layer
            // Thus, it works more like a BFS derivative
            std::vector<int> potentials(nodes_.size(), std::numeric_limits<int>::max());
            {
                    
                std::deque<std::pair<int, int> > front;
                front.push_back({0, source_i});
                
                while(front.size() > 0) {
                    int potential;
                    int cur_i;
                    std::tie(potential, cur_i) = front.front();
                    front.pop_front();
                    
                    if(potential >= potentials[cur_i]) {
                        continue;
                    }
                    potentials[cur_i] = potential;
                    
                    for(Arc* arc : nodes_[cur_i].connected_arcs_) if(arc->get_capacity(cur_i) > 0) {
                        // Traverse the arc if there is some remaining capacity
                        front.push_back({potential + arc->get_cost_from(cur_i), arc->get_dest(cur_i)});
                    }
                }
            }
            
            // Next find the zero cost when modifying with potentials
            {
                int flow_pushed = dinic_with_potentials(source_i, sink_i, potentials);
                result += flow_pushed * (potentials[sink_i] - potentials[source_i]);  // Add the cost of the pushed flow
                if(flow_pushed == 0) {
                    return result;
                }
            }
        }
    }
};


int main() {
    PrimalDualFlowNetwork graph;
    
    for(int i=0; i<6; i++) graph.add_node();  // Initialize 6 nodes
    // Add edges of our graph
    graph.add_arc(0, 1, 0, 4, 1);
    graph.add_arc(0, 3, 0, 2, 5);
    graph.add_arc(1, 2, 0, 2, 1);
    graph.add_arc(1, 3, 0, 6, 1);
    graph.add_arc(2, 1, 0, 2, 1);
    graph.add_arc(2, 5, 0, 4, 0);
    graph.add_arc(3, 4, 0, 8, 1);
    graph.add_arc(4, 2, 0, 6, -3);
    graph.add_arc(4, 5, 0, 4, 1);
    
    int result = graph.min_cost_max_flow(0, 5);  // Find the min-cost max-flow
    std::cout << result << '\n';
    
    return 0;
}

Further reading

Sample problem:
Minimum Cost Maximum Flow - hackerearth.com

References

objects
std::cout en.cppreference.com cplusplus.com
classes
std::deque en.cppreference.com cplusplus.com
std::function en.cppreference.com cplusplus.com
std::numeric_limits en.cppreference.com cplusplus.com
std::tuple en.cppreference.com cplusplus.com
std::vector en.cppreference.com cplusplus.com
functions
std::basic_ostream::operator<< en.cppreference.com cplusplus.com
std::deque::front en.cppreference.com cplusplus.com
std::deque::pop_front en.cppreference.com cplusplus.com
std::deque::push_back en.cppreference.com cplusplus.com
std::function::operator() en.cppreference.com cplusplus.com
std::make_tuple en.cppreference.com cplusplus.com
std::min en.cppreference.com cplusplus.com
std::numeric_limits::max en.cppreference.com cplusplus.com
std::reverse en.cppreference.com cplusplus.com
std::tie en.cppreference.com cplusplus.com
std::vector::operator[] en.cppreference.com cplusplus.com
std::vector::resize en.cppreference.com cplusplus.com
std::vector::size en.cppreference.com cplusplus.com
std::vector::vector en.cppreference.com cplusplus.com

Problem Description

Implement the primal-dual algorithm for finding minimum cost max-flow. Construct a flow network and output the value of the min-cost max-flow on it.

Minimum Cost Flow Part Two: Primal-Dual Algorithm - topcoder.com
Flow network - wikipedia.org

View sample discussion (0 comments)
View problem discussion (0 comments)