Categories

# Push-relabel in C++

In a course on discrete maths I took in my second year, we had to implement the push-relabel algorithm in C++. This is a modified solution that has a worse theoretical, but better practical running time.

``````#include <iostream>
#include <vector>
#include <fstream>
#include <algorithm>
#include <map>

using std::vector;
using node = int;
using edge = int;
using capacity = int;
using flow = int;

struct Edge {
node from;   // from < 2^31
node to;     // to   < 2^31
capacity u;  // capacity < 2^31
flow f;      // flow <= capacity
};

// In an EdgeCollection we store both edges and their inverses.
// A edge is represented by an integer [0,edges.size) and its inverse edge
// by the negative minus one [-edges.size, 0).
class EdgeCollection {
public:
EdgeCollection(vector<Edge>& es): edges{es} {}
edge inverse(edge e) {
return -e-1;
}
node from(edge e) {
return e >= 0 ? edges[e].from : edges[inverse(e)].to;
}
node to(edge e) {
return e >= 0 ? edges[e].to : edges[inverse(e)].from;
}
node capacity_left(edge e) {
return e >= 0 ? edges[e].u - edges[e].f : get_flow(inverse(e));
}
flow get_flow(edge e) {
return e >= 0 ? edges[e].f : capacity_left(inverse(e));
}
void augment(edge e, flow f) {
if(e >= 0)    edges[e].f += f;
else edges[inverse(e)].f -= f;
}
std::size_t size() {
return edges.size();
}
private:
vector<Edge> edges;
};

// GoldbergTarjan contains the main push_relabel algorithm.
// We maintain the following invariants:
//  - incident[v] contains a list of edges starting in v,
//    inverse or otherwise, without duplicates.
//  - ex[v] stores the excess at v.
//  - phi[v] stores the distance phi(v).
//  - L[i] stores the active nodes v with phi[v] == i
//    without duplicates.
//  - A[v] stores a vector of allowed edges starting in v,
//    inverse or otherwise, without duplicates.
//    We may keep edges that were allowed once but aren't anymore.
class GoldbergTarjan {
public:
GoldbergTarjan(int num_nodes, EdgeCollection ec, node source, node t_node)
:phi(num_nodes,0), s{source}, t{t_node}, incoming(num_nodes),
outgoing(num_nodes), A(num_nodes), edges{ec}, ex(num_nodes, 0) {
for(std::size_t i = 0; i < edges.size(); ++i) {
if(edges.from(i) == s) {
edges.augment(i, edges.capacity_left(i));
if(edges.to(i) != t and edges.get_flow(i) > 0)
L.insert({0, edges.to(i)});
ex[edges.to(i)] += edges.get_flow(i);
ex[s] -= edges.get_flow(i);
}
outgoing[edges.from(i)].push_back(i);
incoming[edges.to(i)].push_back(edges.inverse(i));
}
phi[s] = num_nodes; // phi[v] == 0 else (see initialization)
}

EdgeCollection push_relabel() {
node v = get_maximum_active_node();
bool still_active = false;
while(v != -1) {
if(A[v].size() != 0) {
still_active = push(A[v].back());
} else {
relabel(v);
// v is active and since it was the maximum active node before
// it will be now (phi[v] never decreases).
still_active = true;
}
v = still_active ? v : get_maximum_active_node();
}
return edges;
}

flow value() {
return -ex[s];
}
private:
vector<unsigned> phi;
node s,t;
vector<vector<edge>> incoming;
vector<vector<edge>> outgoing;
std::multimap<int, node> L;
vector<vector<edge>> A;
EdgeCollection edges;
vector<int> ex;

bool push(edge e) { // returns true iff ex[v] remains > 0
node v = edges.from(e);
node w = edges.to(e);
if(phi[v] != phi[w] + 1) {
A[v].pop_back(); // We remove not-any-more-allowed edges
return true;
}
bool still_active = true;
flow gamma = std::min(ex[v], edges.capacity_left(e));
edges.augment(e, gamma);
ex[v] -= gamma;
ex[w] += gamma;
if(ex[v] == 0) {
L.erase(--L.end());
still_active = false;
}
if(edges.capacity_left(e) == 0)          A[v].pop_back();
if(ex[w] == gamma and w != s and w != t) L.insert({phi[w], w});
return still_active;
}

void relabel(node v) {
unsigned m = 2*phi.size();
L.erase(--L.end());
// Splitting outgoing and incoming is better for branch-prediction.
std::vector<vector<edge>*> incident({&(outgoing[v]), &(incoming[v])});
for(auto* vec : incident)
for(auto e : *vec)
if(edges.capacity_left(e) > 0)
m = std::min(m, phi[edges.to(e)] + 1);
phi[v] = m;
L.insert({phi[v],v});
for(auto* vec : incident)
for(auto e : *vec)
if(phi[v] == phi[edges.to(e)] + 1 and edges.capacity_left(e) > 0)
A[v].push_back(e);
}
node get_maximum_active_node() {
if(L.size() != 0) return L.rbegin()->second;
return -1;
}
};

int main(int argc, char *argv[]) {
if(argc != 2) {
// We expect exactly one argument...
std::cout << "Aufruf: " << argv[0] << " <filename>" << std::endl;
} else {
// which should be a filename...
std::ifstream file(argv[1]);
unsigned int n = 0;
vector<Edge> edges(0);
if(!file.is_open()) {
// of a non-busy file.
std::cout << "Konnte Datei nicht öffnen." << std::endl;
} else {
file >> n; // First we read the number of nodes...
node a, b; capacity c;
while(file >> a >> b >> c) {
// and then the edges...
edges.push_back(Edge {a,b,c});
}
}
GoldbergTarjan gt(n, EdgeCollection(edges), 0, 1);
EdgeCollection result = gt.push_relabel();
std::cout << gt.value() << std::endl;
for(std::size_t i = 0; i < result.size(); ++i) {
if(result.get_flow(i) > 0) {
std::cout << i << " " << result.get_flow(i) << std::endl;
}
}
}
}``````