All Pairs Shortest Paths

apsp.cpp

#include <iostream>
#include <vector>
#include <random>
#include <numeric>
#include <algorithm>
#include <cassert>
#include <tuple>
#include <utility>
#include <map>
#include <iomanip>
#include <queue>
#include <limits>

using namespace std;

namespace hidden_path_algorithm {
	template<class T=int>
	struct edge {
		using dist_type = T;
		using self_type = edge<T>;
		size_t from;
		size_t to;
		dist_type dist;
		size_t next;

		bool operator < (const self_type& x) const {
			return this->dist < x.dist;
		}
		bool operator > (const self_type& x) const {
			return this->dist > x.dist;
		}
	};

	template<class T=int>
	class graph {
		using dist_type = T;
		using path = edge<T>;
		const dist_type inf = numeric_limits<dist_type>::max() / 2;

		size_t n_;
		vector<path> e_;
		bool is_apsp_calculated_;

	public:
		vector<vector<path>> shortest_path_;

		graph(size_t n): n_(n), e_(), is_apsp_calculated_(false) {}

		void add_edge(size_t from, size_t to, dist_type dist){
			e_.push_back(path{
				from,
				to,
				dist,
				to
			});
		}

		vector<vector<path>>& all_pairs_shortest_paths(){
			assert(is_apsp_calculated_ == false);
			is_apsp_calculated_ = true;

			shortest_path_ = vector<vector<path>>(n_, vector<path>(n_));
			for(size_t i=0; i<n_; i++) for(size_t j=0; j<n_; j++){
				if(i==j){
					shortest_path_[i][j] = path{
						i,i,0,i
					};
				}else{
					shortest_path_[i][j] = path{
						i,j,inf,j
					};
				}
			}
			for(const auto& p: e_){
				if(p < shortest_path_[p.from][p.to]){
					shortest_path_[p.from][p.to] = p;
				}
			}

			vector<vector<path>> optimal_into_edges(n_);
			vector<vector<path>> discovered_paths(n_);
			for(int i=0; i<(int)n_; i++) discovered_paths[i].reserve(n_);
			priority_queue<path, vector<path>, greater<path>> pq(begin(e_), end(e_));

			auto update = [&](const path& p, const path& q){
				path next_path{
					p.from,
					q.to,
					p.dist + q.dist,
					p.next,
				};
				if(next_path < shortest_path_[next_path.from][next_path.to]){
					shortest_path_[next_path.from][next_path.to] = next_path;
					pq.push(next_path);
				}
			};

			while(!pq.empty()){
				path p = pq.top();
				pq.pop();
				if(shortest_path_[p.from][p.to] < p) continue;

				discovered_paths[p.from].push_back(p);
				if(p.next == p.to){ // is edge
					optimal_into_edges[p.to].push_back(p);

					for(const path& q: discovered_paths[p.to]){
						update(p, q);
					}
				}
				for(const path& r: optimal_into_edges[p.from]){
					update(r, p);
				}
			}

			return shortest_path_;
		}

		dist_type get_dist(size_t s, size_t t) const {
			assert(is_apsp_calculated_ == true);
			dist_type dist = shortest_path_[s][t].dist;
			return dist;
		}

		pair<dist_type, vector<size_t>> get_path(size_t s, size_t t) const {
			assert(is_apsp_calculated_ == true);
			dist_type dist = shortest_path_[s][t].dist;
			if(dist == inf) return {dist, {}};

			vector<size_t> path{s};
			while(s != t){
				s = shortest_path_[s][t].next;
				path.push_back(s);
			}
			return {dist, path};
		}
	};
}


namespace SHORT_Algorith {
	template<class T=int>
	struct edge {
		using dist_type = T;
		using self_type = edge<T>;
		size_t from;
		size_t to;
		dist_type dist;
		size_t next;

		bool operator < (const self_type& x) const {
			return this->dist < x.dist;
		}
		bool operator > (const self_type& x) const {
			return this->dist > x.dist;
		}
	};

	template<class C=int>
	class graph {
		using dist_type = C;
		using path = edge<C>;
		const dist_type inf = numeric_limits<dist_type>::max() / 2;

		size_t n_;
		vector<path> e_;
		bool is_apsp_calculated_;

	public:
		vector<vector<path>> shortest_path_;
		vector<vector<path>> essential_subgraph_;

		graph(size_t n): n_(n), e_(), is_apsp_calculated_(false) {}

		void add_edge(size_t from, size_t to, dist_type dist){
			e_.push_back(path{
				from,
				to,
				dist,
				to,
			});
		}

		vector<vector<path>>& all_pairs_shortest_paths(){
			assert(is_apsp_calculated_ == false);
			is_apsp_calculated_ = true;

			shortest_path_ = vector<vector<path>>(n_, vector<path>(n_));
			for(size_t i=0; i<n_; i++) for(size_t j=0; j<n_; j++){
				if(i==j){
					shortest_path_[i][j] = path{ i,i,0,i };
				}else{
					shortest_path_[i][j] = path{ i,j,inf,j };
				}
			}

			vector<vector<int>> T(n_, vector<int>(n_, -1));
			vector<vector<size_t>> T_nodes(n_);
			for(int i=0; i<(int)n_; i++){
				T[i][i] = 0;
				T_nodes[i] = {(size_t)i};
				T_nodes[i].reserve(n_);
			}
			using PQ = priority_queue<path, vector<path>, greater<path>>;
			vector<PQ> fringe(n_);
			vector<vector<path>> H(n_);

			PQ pq(begin(e_), end(e_));

			auto dijkstra_process = [&](size_t v, size_t z, path e){
				auto d = e.dist + shortest_path_[v][z].dist;
				path next_path{v, e.to, d, shortest_path_[v][z].next}; 
				if(v == z) next_path = e;
				if(shortest_path_[v][e.to] > next_path){
					fringe[v].push(next_path);
				}
			};

			auto exists_altenate_path = [&](path p, bool is_postprocess = false){
				const bool ACCEPT = false;
				const bool REJECT = true;
				size_t v = p.from;
				if(!is_postprocess && shortest_path_[p.from][p.to].dist <= p.dist) return REJECT;
				for(int z: T_nodes[v]){
					for(auto& i=T[v][z]; i<(int)H[z].size(); i++){
						auto& e = H[z][i];
						dijkstra_process(v, z, e);
					}
				}
				while(true){
					if(fringe[v].empty()) return ACCEPT;
					auto e = fringe[v].top();
					if(!is_postprocess && e > p) return ACCEPT;
					fringe[v].pop();
					size_t z = e.to;
					if(T[v][z] >= 0) continue;
					shortest_path_[v][z] = e;
					T_nodes[v].push_back(z);
					T[v][z] = 0;
					for(auto& i=T[v][z]; i<(int)H[z].size(); i++){
						auto& e = H[z][i];
						dijkstra_process(v, z, e);
					}

					if(!is_postprocess && z == p.to) return REJECT;
				}
			};

			while(!pq.empty()){
				path p = pq.top();
				pq.pop();

				if(exists_altenate_path(p) == false) {
					H[p.from].push_back(p);
				}
			}

			// postprocessing
			for(size_t i=0; i<n_; i++){
				exists_altenate_path(path{i, i, -1}, true);
			}

			essential_subgraph_ = H;
			return shortest_path_;
		}

		dist_type get_dist(size_t s, size_t t) const {
			assert(is_apsp_calculated_ == true);
			dist_type dist = shortest_path_[s][t].dist;
			return dist;
		}
	};

}

namespace SKG{
struct AliasMethod{
  int n;
  vector<double> values;
  vector<double> thresholds;
  vector<int> low;
  vector<int> high;
  AliasMethod(){}

  AliasMethod(vector<double> v)
   : n(v.size()), values(v), thresholds(v.size()), low(v.size()), high(v.size())
  {
    double s = accumulate(values.begin(), values.end(), 0.0);
    s = 1.0 / s;
    for(auto& x: values){
      x = x * s * n;
    }

    vector<int> small, big, one;
    for(int i=0; i<n; i++){
      if(values[i] > 1.0){
        big.push_back(i);
      }else if(values[i] < 1.0){
        small.push_back(i);
      }else{
        one.push_back(i);
      }
    }

    for(int i=0; i<n; i++){
      if(one.size()){
        int a = one.back();
        one.pop_back();
        thresholds[i] = 1.0;
        low[i] = a;
        high[i] = a;
        values[a] = 0.0;
      }else if(small.size() && big.size()){
        int a = small.back();
        small.pop_back();
        int b = big.back();
        big.pop_back();

        thresholds[i] = values[a];
        low[i] = a;
        high[i] = b;
        values[b] -= 1.0 - values[a];
        values[a] = 0.0;
        if(values[b] > 1.0){
          big.push_back(b);
        }else if(values[b] == 1.0){
          one.push_back(a);
        }else if(values[b] > 0.0){
          small.push_back(b);
        }
      }else{
        swap(one, big.size() ? big : small);
        i--;
      }
    }
  }

  int get(double p){
    assert(0.0 <= p && p < 1.0);
    p *= n;
    int i = p;
    p -= i;
    assert(0.0 <= p && p < 1.0);
    if(thresholds[i] > p){
      return low[i];
    }else{
      return high[i];
    }
  }
};


template<class T>
struct edge {
  int from;
  int to;
  T weight;
};


// Stochastic Kronecker Graph Generation
class SKG {
  int n;
  AliasMethod a;
  mt19937 mt;
  uniform_real_distribution<double> dst;

  static vector<double> flatten_weight_mat(const vector<vector<double>>& w){
    int n = w.size();
    assert((int)w[0].size() == n);

    vector<double> weight(n*n);
    for(int i=0; i<n; i++) for(int j=0; j<n; j++){
      weight[i*n + j] = w[i][j];
    }
    return weight;
  }

  tuple<int, int, double> next_edge(int K){
    K--;
    int u=0, v=0;
    double p = dst(mt);
    int w = a.get(p);
    u = w/n;
    v = w%n;

    if(K > 0){ 
      auto tmp = next_edge(K);
      u = u + get<0>(tmp) * n;
      v = v + get<1>(tmp) * n;
      p = p * get<2>(tmp);
    }
    return make_tuple(u, v, p);
  }

public:
  SKG(const vector<vector<double>>& weight, long long seed):
    n(weight.size()),
    a(flatten_weight_mat(weight)),
    mt(seed),
    dst(0.0, 1.0)
  { }

  // generate graph where |G| = n^k
  vector<edge<double>> generate(
    int K, 
    int edge_factor
  ){
    long long m = 1;
    for(int k=0; k<K; k++, m=m*n){}
    vector<edge<double>> res;

  	vector<vector<bool>> used(m, vector<bool>(m));

    int num_vertex = m;
    cerr << "num_vertex = " << num_vertex << endl;
    int num_target_edges = m * edge_factor;
    cerr << "num_target_edges = " << num_target_edges << endl;

    for(int itr=0; itr < num_target_edges; itr++){
      int u, v;
      double w;
      tie(u, v, w) = next_edge(K);

      if(u==v) continue;
      // if(u>v) swap(u,v);
      if(used[u][v]) continue;
      used[u][v] = true;
      res.push_back(edge<double>{u, v, w});
    }
    return res;
  }
};

}


#include <sstream>

class Graphviz {
public:
	stringstream data;
	Graphviz():data(){}

	void add_edge(int from, int to, int weight, bool visible=true){
		if(visible){
	    data << from << " -> " << to << " [weight=" << weight << "];\n";
	  }else{
	    data << from << " -> " << to << " [weight=" << weight << ", style=invis ];\n";
	  }
	}

	void write(ostream& os){
		os << "digraph {\n";
		os << "node[shape=circle];\n";
		os << data.str();
		os << "}" << endl;
	}
};

#include <random>
#include <chrono>
#include <fstream>
#include <set>

int main(int argc, char* argv[]){
  long long seed = chrono::duration_cast<chrono::nanoseconds>(chrono::steady_clock::now().time_since_epoch()).count();
  long long depth = 2;
  long long edge_factor = 2;

  bool output_graph = false;
  for(int i=1; i<argc; i++){
    if(argv[i] == "-seed"s){
      seed = stoll(argv[++i]);
    }else if(argv[i] == "-depth"s){
      depth = stoll(argv[++i]);
    }else if(argv[i] == "-e"s){
      edge_factor = stoll(argv[++i]);
    }else if(argv[i] == "-o"s){
    	output_graph = true;
    	cerr << "OUTPUT ENABLED" << endl;
    }
  }

  vector<vector<double>> w = {
    {0.57, 0.19},
    {0.19, 0.05},
  };

	SKG::SKG skg(w, seed);

	int n = 1;
	for(int i=0; i<depth; i++, n*=w.size());

	auto generated_edges = skg.generate(depth, edge_factor);

	int inf = 1<<28;
	vector<vector<int>> dist(n, vector<int>(n, inf));
	for(int i=0; i<n; i++){
		dist[i][i] = 0;
	}
	for(auto e: generated_edges){
		dist[e.from][e.to] = min(dist[e.from][e.to], (int)(e.weight * 100000));
		// dist[e.to][e.from] = min(dist[e.to][e.from], (int)(e.weight * 100000));
	}

	// wf
	auto dist_wf = dist;
	{
		auto s0 = chrono::steady_clock::now();
		for(int k=0; k<n; k++) for(int i=0; i<n; i++) for(int j=0; j<n; j++){
			dist_wf[i][j] = min(dist_wf[i][j], dist_wf[i][k]+dist_wf[k][j]);
		}
		auto s1 = chrono::steady_clock::now();
		double elapsed_s = chrono::duration_cast<chrono::milliseconds>(s1 - s0).count();
		cerr << "WARSHALL_FLOYD: " << elapsed_s << "[ms]" << endl;
	}
	{
		SHORT_Algorith::graph<int> g(n);
		for(auto e: generated_edges){
			g.add_edge(e.from, e.to, (int)(e.weight * 100000));
			// g.add_edge(e.to, e.from, (int)(e.weight * 100000));
		}

		auto t0 = chrono::steady_clock::now();
		g.all_pairs_shortest_paths();
		auto t1 = chrono::steady_clock::now();
		double elapsed_t = chrono::duration_cast<chrono::milliseconds>(t1 - t0).count();
		cerr << "SHORT         : " << elapsed_t << "[ms]" << endl;

		for(int i=0; i<n; i++) for(int j=0; j<n; j++){
			if(dist_wf[i][j] == inf){
				assert(g.get_dist(i,j) >= inf);
			}else{
				assert(dist_wf[i][j] == g.get_dist(i,j));
			}
		}

		if(output_graph){
			ofstream ofs("graph.dot");
			Graphviz gv;
			for(auto e: generated_edges)
				gv.add_edge(e.from, e.to, e.weight * 100000);
			gv.write(ofs);
			ofs.close();
		}
		if(output_graph){
			ofstream ofs("essential_subgraph.dot");
			Graphviz gv;
			set<pair<int,int>> essential;
			for(int i=0; i<n; i++) for(auto e: g.essential_subgraph_[i])
				essential.insert({e.from, e.to});
			for(auto e: generated_edges)
				gv.add_edge(e.from, e.to, e.weight * 100000, essential.count({e.from, e.to}));
			gv.write(ofs);
			ofs.close();
		}
	}

	{
		hidden_path_algorithm::graph<int> g(n);
		for(auto e: generated_edges){
			g.add_edge(e.from, e.to, (int)(e.weight * 100000));
			// g.add_edge(e.to, e.from, (int)(e.weight * 100000));
		}

		auto t0 = chrono::steady_clock::now();
		g.all_pairs_shortest_paths();
		auto t1 = chrono::steady_clock::now();
		double elapsed_t = chrono::duration_cast<chrono::milliseconds>(t1 - t0).count();
		cerr << "HIDDEN_PATH   : " << elapsed_t << "[ms]" << endl;

		for(int i=0; i<n; i++) for(int j=0; j<n; j++){
			if(dist_wf[i][j] == inf){
				assert(g.get_dist(i,j) >= inf);
			}else{
				assert(dist_wf[i][j] == g.get_dist(i,j));
			}
		}
	}

	// dijkstra
	{
		auto dist_d = dist;
		vector<vector<pair<int,int>>> g(n);
		for(auto e: generated_edges){
			g[e.from].push_back({(int)(e.weight * 100000), e.to});
			// g[e.to].push_back({(int)(e.weight * 100000), e.from});
		}
		auto r0 = chrono::steady_clock::now();
		for(int s=0; s<n; s++){
			priority_queue<pair<int,int>, vector<pair<int,int>>, greater<pair<int,int>>> pq;
			auto& d = dist_d[s];
			d = vector<int>(n, inf);
			d[s] = 0;
			for(auto& e: g[s]){
				d[e.second] = min(d[e.second], e.first);
				pq.push(e);
			}
			while(!pq.empty()){
				auto p = pq.top();
				pq.pop();
				if(d[p.second] < p.first) continue;
				for(auto& e: g[p.second]){
					int dd = p.first + e.first;
					if(dd < d[e.second]){
						d[e.second] = dd;
						pq.push({dd, e.second});
					}
				}
			}
		}
		auto r1 = chrono::steady_clock::now();
		double elapsed_r = chrono::duration_cast<chrono::milliseconds>(r1 - r0).count();
		cerr << "DIJKSTRA      : " << elapsed_r << "[ms]" << endl;

		for(int i=0; i<n; i++) for(int j=0; j<n; j++){
			if(dist_wf[i][j] == inf){
				assert(dist_d[i][j] >= inf);
			}else{
				assert(dist_wf[i][j] == dist_d[i][j]);
			}
		}
	}


	return 0;
}