Categories
Code Dump

Successive shortest paths in C++

This solves the assignment problem by solving the associated minimum cost flow using the successive shortest paths method.

#include <iostream>
#include <vector>
#include <fstream>
#include <algorithm>
#include <limits>

using std::vector;

using node = int;
using edge = int;
using capacity = int;
using flow = int;
using cost = int;

struct Edge {
  node from;   // from < 2^31
  node to;     // to   < 2^31
  cost c;      // cost < 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;
  }
  cost get_cost(edge e) {
    return e >= 0 ? edges[e].c : - edges[inverse(e)].c;
  }
  capacity 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();
  }
  vector<Edge>& get_edges() {
    return edges;
  }
private:
  vector<Edge> edges;
};

class AssignmentProblem {
public:
  AssignmentProblem(int num_nodes, vector<Edge>& edges)
    : ec(edges), pi(num_nodes + 2, 0), num_nodes{num_nodes}, s{num_nodes}, t{num_nodes + 1} {
    // Since all the edges go from A to B, 
    // we can calculate the initial potential pi in linear time:
    cost max = 0;
    for(size_t i = 0; i < edges.size(); ++i) {
      pi[edges[i].from] = - std::min( - pi[edges[i].from], edges[i].c);
      max = - std::min(edges[i].c, - max);
    }
    pi[s] = max;
    for(int i = 0; i < num_nodes/2; ++i) {
      edges.push_back(Edge {s,i,0,1,0});
    }
    for(int i = num_nodes/2; i < num_nodes; ++i) {
      edges.push_back(Edge {i,t,0,1,0});
    }
    ec = EdgeCollection(edges);
  }

  vector<Edge> succ_shortest_paths() {
    for(int b_of_s = num_nodes/2; b_of_s > 0; --b_of_s) {
      vector<int>  l(num_nodes + 2, std::numeric_limits<int>::max());
      vector<bool> r(num_nodes + 2, false);
      vector<edge> p(num_nodes + 2, -1);
      l[s] = 0;
      node next = s;
      vector<vector<size_t>> incident(num_nodes + 2);
      for(size_t i = 0; i < ec.size(); ++i) {
        if(ec.capacity_left(i) > 0) {
          incident[ec.from(i)].push_back(i);
        }
        if(ec.get_flow(i) > 0) {
          incident[ec.to(i)].push_back(ec.inverse(i));
        }
      }
      while(next != -1) { // Dijkstra
        r[next] = true;

        for(size_t i = 0; i < incident[next].size(); ++i) {
          edge e = incident[next][i]; node w = ec.to(e);
          cost c = l[next] + ec.get_cost(e) + pi[next] - pi[w];
          if(!r[w] and l[w] > c) {
            l[w] = c; p[w] = incident[next][i];
          }
        }

        next = -1; cost min = std::numeric_limits<int>::max();
        for(size_t i = 0; i < r.size(); ++i) {
          if(!r[i] and l[i] < min) {
            min = l[i];
            next = i;
          }
        }
      }

      node w = t;
      while(w != s) {
        edge way = p[w];
        ec.augment(way, 1); // In every iteration gamma will be 1.
        w = ec.from(way);
      }

      for(size_t i = 0; i < pi.size(); ++i) {
        pi[i] += l[i];
      }
    }
    return ec.get_edges();
  }
private:
  EdgeCollection ec;
  vector<int> pi;
  int num_nodes;
  node s, t;
};

void solve(int n, vector<Edge>& edges) {
  AssignmentProblem ap(n, edges);
  vector<Edge> f = ap.succ_shortest_paths();
  int total_cost = 0;
  for(size_t i = 0; i + n < f.size(); ++i) if(f[i].f == 1) total_cost += f[i].c;
  std::cout << total_cost << std::endl;
  for(size_t i = 0; i + n < f.size(); ++i) {
    Edge e = f[i];
    if(e.f == 1) {
      std::cout << e.from << " " << e.to << std::endl;
    }
  }
}

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; cost c;
      while(file >> a >> b >> c) {
        // and then the edges...
        edges.push_back(Edge {a,b,c,1,0});
      }
      solve(n, edges);
    } 
  }
}
Categories
Code Dump

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;
      }
    }
  }
}
Categories
Code Dump

Kruskals algorithm in C++

This is a faster version of the standard algorithm for finding a minimum spanning tree using a Union-Find datastructure.

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

using std::vector;

// A Branching object is a Union-Find-Structure with elements
// of type Branching::node (unsigned int).
// Union and Find are as in the lecture, but there is no MakeSet
// method as we know in advance how many elements we will need.
class Branching {
  public:
    using node = unsigned int;
    // Create a Branching of size n.
    Branching(unsigned int n) {
      size = n;
      rank = vector<node>(n,0); // Initially ranks are 0
      parent = vector<node>(n,0);
      for(node i = 0; i < n; ++i) {
        parent[i] = i; // Initially each element is its own parent.
      }
    }

    // Find with path-compression
    node Find(node x) {
      if(x != parent[x]) {
        parent[x] = Find(parent[x]);
      }
      return parent[x];
    }

    // Union as in lecture
    void Union(node x, node y) {
      node x_bar = Find(x);
      node y_bar = Find(y);
      if(rank[x_bar] > rank[y_bar]) {
        parent[y_bar] = x_bar;
      } else {
        parent[x_bar] = y_bar;
        if(rank[x_bar] == rank[y_bar]) {
          ++rank[y_bar];
        }
      }
    }
  private:
    unsigned int size;
    vector<node> parent;
    vector<node> rank;
};

struct edge {
  Branching::node from; // from < 2^32
  Branching::node to;   // to   < 2^32
  long long cost;       // abs(cost) < 2^32
};

// Kruskals algorithm
vector<edge> kruskal(unsigned int n, vector<edge>& edges) {
  // Lambda functions were introduced in C++11,
  // enable -std=c++11 if this fails to compile.
  std::sort(edges.begin(), edges.end(), [](edge e1, edge e2) {
    return e1.cost < e2.cost;
  });
  Branching nodes(n);
  vector<edge> mst(0);
  for(auto e : edges) {
    // We identify connected components in mst with components
    // in the Branching structure.
    // If the edge connects two disjoint components...
    if(nodes.Find(e.from) != nodes.Find(e.to)) {
      // add the edge to mst and union them in nodes.
      nodes.Union(e.from, e.to);
      mst.push_back(e);
    }
  }
  return mst;
}

// Output the result of Kruskals algorithm.
void output(unsigned int n, vector<edge> mst) {
  if(mst.size() != n - 1) {
    // A minimum spanning tree has exactly n-1 edges
    // if mst has less, the graph can't be connected.
    std::cout << "Der Graph ist nicht zusammenhängend!" << std::endl;
  } else {
    long long sum = 0;
    for(auto e : mst) {
      sum += e.cost; // Sum the costs.
    }
    std::cout << "Es gibt einen MST mit Gewicht " << sum << "." << std::endl;
    for(auto e : mst) {
      // Give back each edge in the order in which they were added to mst.
      std::cout << e.from << " -> " << e.to << " : "  << e.cost << std::endl;
    }
  }
}

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...
      Branching::node a, b; long long c;
      while(file >> a >> b >> c) {
        // and then the edges...
        edges.push_back(edge {a,b,c});
      }
      // and then we compute the MST.
      output(n, kruskal(n, edges));
    }
  }
}