Johnson's algorithm for all-pairs shortest paths with non-integer (double) weights allowed.

johnson.cpp

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <climits>
#include <float.h>
#include <exception>
#include <set>

// This implementation of "Johnson's algorithm" 
// is based on Ashley Holman's https://gist.github.com/ashleyholman/6793360
// although there was no apparent copyright.
// 
// I (Brian Gauch) modified it so that weights are doubles.
// More compiler-friendly without "using" keyword.
// Also prints out all shortest paths at the end.
// There was no apparent copyright on the original implementation
// by Ashley Homlan, although it was shared publicly.
// My contribution is under the MIT license (below), 
// if possible given that is based on previous work.

/*
Copyright 2017 Brian Gauch

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.
*/

using namespace std;

struct Edge {
  int head;
  double cost;
};

//using Graph = vector<vector<Edge>>;
//using SingleSP = vector<long>;
//using AllSP = vector<vector<long>>;
const double INF = DBL_MAX;

// first line is in the format <N> <M> where N is the number of vertices, M is the number of edges that follow
// each following line represents a directed edge in the form <Source> <Dest> <Cost>
vector<vector<Edge>> loadgraph(istream& is) {
  vector<vector<Edge>> g;

  int n, m;
  is >> n >> m;
  cout << "Graph has " << n << " vertices and " << m << " edges." << endl;
  g.resize(n+1);
  while (is) 
  {
    int v1, v2;
    double c;
    is >> v1 >> v2 >> c;
    //cout << "got cost:" << c << endl;
    if (is) 
    {
      g[v1].push_back({v2, c});
    }
  }

  return g;
}

vector<vector<Edge>> addZeroEdge(vector<vector<Edge>> g) {
  // add a zero-cost edge from vertex 0 to all other edges
  for (int i = 1; i < g.size(); i++) {
    g[0].push_back({i, 0.0});
  }

  return g;
}

vector<double> bellmanford(vector<vector<Edge>> &g, int s) {
  vector<vector<double>> memo(g.size()+2, vector<double>(g.size(), INF));

  // initialise base case
  memo[0][s] = 0.0;

  for (int i = 1; i < memo.size(); i++) {
    // compute shortest paths from s to all vertices, with max hop-count i
    for (int n = 0; n < g.size(); n++) {
      if (memo[i-1][n] < memo[i][n]) {
        memo[i][n] = memo[i-1][n];
      }
      for (auto& e : g[n]) {
        if (memo[i-1][n] != INF) {
          if (memo[i-1][n] + e.cost < memo[i][e.head]) {
            memo[i][e.head] = memo[i-1][n] + e.cost;
          }
        }
      }
    }
  }

  // check if the last iteration differed from the 2nd-last
  for (int j = 0; j < g.size(); j++) {
    if (memo[g.size()+1][j] != memo[g.size()][j]) {
      throw string{"negative cycle found"};
    }
  }

  return memo[g.size()];
}

vector<double> djikstra(const vector<vector<Edge>>& g, int s) {
  vector<double> dist(g.size(), INF);
  set<pair<double,int>> frontier;

  frontier.insert({0.0,s});

  while (!frontier.empty()) {
    pair<double,int> p = *frontier.begin();
    frontier.erase(frontier.begin());

    double d = p.first;
    //cout << "d=" << d << endl;
    int n = p.second;

    // this is our shortest path to n
    dist[n] = d;

    // now look at all edges out from n to update the frontier
    for (auto e : g[n]) {
      // update this node in the frontier if we have a shorter path
      if (dist[n]+e.cost < dist[e.head]) {
        if (dist[e.head] != INF) {
          // we've seen this node before, so erase it from the set in order to update it
          frontier.erase(frontier.find({dist[e.head], e.head}));
        }
        frontier.insert({dist[n]+e.cost, e.head});
        dist[e.head] = dist[n]+e.cost;
      }
    }
  }

  return dist;
}

vector<vector<double>> johnson(vector<vector<Edge>> &g) {
  // now build "g prime" which is g with a new edge added from vertex 0 to all other edges, with cost 0
  vector<vector<Edge>> gprime = addZeroEdge(g);

  // now run Bellman-Ford to get single-source shortest paths from s to all other vertices
  vector<double> ssp;
  try {
    ssp = bellmanford(gprime, 0);
  } catch (string e) {
    cout << "Negative cycles found in graph.  Cannot compute shortest paths." << endl;
    throw e;
  }

  // no re-weight each edge (u,v) in g to be: cost + ssp[u] - ssp[v].
  for (int i = 1; i < g.size(); i++) {
    for (auto &e : g[i]) {
      e.cost = e.cost + ssp[i] - ssp[e.head];
    }
  }

  // now that we've re-weighted our graph, we can invoke N iterations of Djikstra to find
  // all pairs shortest paths
  vector<vector<double>> allsp(g.size());
  for (int i = 1; i < g.size(); i++) {
    allsp[i] = djikstra(g, i);
  }

  // now re-weight the path costs back to their original weightings
  for (int u = 1; u < g.size(); u++) {
    for (int v = 1; v < g.size(); v++) {
      if (allsp[u][v] != INF) {
        allsp[u][v] += ssp[v] - ssp[u];
      }
    }
  }

  return allsp;
}

int main (int argc, char *argv[]) {
  if (argc < 2) {
    cerr << "Usage: " << argv[0] << " <infile>" << endl;
    return 1;
  }

  ifstream is {argv[1]};
  if (!is) {
    cerr << "Couldn't open input file: " << argv[1] << endl;
    return 1;
  }

  vector<vector<Edge>> g = loadgraph(is);

  // run Johnson's algorithm to get all pairs shortest paths
  vector<vector<double>> asp = johnson(g);
  // display all shortest paths
  for (int i = 1; i < g.size(); i++)
  {
    for (int j = 1; j < g.size(); j++)
    {
      cout << i << "==>" << j << ":" << std::setprecision (5) << asp[i][j] << endl;
    }
  }

  // find the "shortest shortest path", ie. the path with lowest cost,
  // amongst all shortest path pairs.
  double shortest = INF;
  for (int i = 1; i < g.size(); i++) {
    for (int j = 1; j < g.size(); j++) {
      if (asp[i][j] < shortest) {
        shortest = asp[i][j];
      }
    }
  }

  cout << "Shortest shortest path = " << shortest << endl;
}