ppsp_test.cc

#include <cstdio>
#include <fstream>
#include <functional>
#include <numeric>
#include <sstream>

#include "connectivity.h"
#include "hashbag.h"
#include "parlay/primitives.h"
#include "parlay/sequence.h"

namespace sssp {
#include "sssp.h"
}
namespace ppsp {
#include "ppsp.h"
}
namespace astar {
#include "astar.h"
}

#define DISTANCE

using namespace std;
using namespace parlay;

constexpr int NUM_SRC = 1;
constexpr int NUM_ROUND = 5;

bool essentiallyEqual(float a, float b) {
  constexpr float eps = 1e-3;
  return fabs(a - b) <= ((fabs(a) > fabs(b) ? fabs(b) : fabs(a)) * eps);
}

bool check_equal(const EdgeTy &d1, const EdgeTy &d2) {
  if constexpr (is_integral_v<EdgeTy>) {
    return d1 == d2;
  } else {
    return essentiallyEqual(d1, d2);
  }
}

sequence<double> read_coords(char const *filename, size_t n) {
  auto chars = chars_from_file(string(filename));
  auto tokens_seq = tokens(chars);
  ostringstream oss;
  oss << tokens_seq[0];
  auto header = oss.str();

  string default_header = "pbbs_sequencePoint";
  // pbbs_sequencePoint_${dim}_d
  assert(header.substr(0, default_header.size()) == default_header);
  size_t dim = stoull(header.substr(default_header.size(),
                                    header.size() - default_header.size() - 1));
  assert(tokens_seq.size() == dim * n + 1);

  sequence<double> coords(n * dim);
  parallel_for(0, n, [&](size_t i) {
    for (size_t j = 0; j < dim; ++j) {
      coords[i * dim + j] = chars_to_double(tokens_seq[dim * i + j + 1]);
    }
  });
  return coords;
}

template <class Func, class Algo>
void run(Func &&f, Algo &&algo) {
  double total_time = 0;
  for (int i = 0; i <= NUM_ROUND; i++) {
    algo.reset();
    internal::timer tm;
    f();
    tm.stop();
    if (i == 0) {
      printf("Warmup Round: %f\n", tm.total_time());
    } else {
      printf("Round %d: %f\n", i, tm.total_time());
      total_time += tm.total_time();
    }
  }
  double average_time = total_time / NUM_ROUND;
  printf("Average time: %f\n", average_time);

  ofstream ofs("ppsp.tsv", ios_base::app);
  ofs << fixed << setprecision(6);
  ofs << average_time << '\t';
  ofs.close();
}



#ifdef DISTANCE
sequence<pair<EdgeTy, NodeId>> select_dest(const Graph G, NodeId s) {
  sssp::Rho_Stepping solver(G);
  size_t n = G.n;
  solver.reset();
  auto dist = solver.sssp(s);
  auto indices = pack_index<NodeId>(
      delayed_seq<bool>(n, [&](size_t i) { return dist[i] != DIST_MAX; }));
  auto dist_node = tabulate(indices.size(), [&](size_t i) {
    return make_pair(dist[indices[i]], indices[i]);
  });
  sort_inplace(make_slice(dist_node));
  return dist_node;
}

void run(const Graph &G, [[maybe_unused]] const sequence<double> &coords,
         const sequence<NodeId> &largest_cc) {
  size_t m = largest_cc.size();
  for (int v = 0; v < NUM_SRC; v++) {
    NodeId s = largest_cc[hash32(v) % m];
    auto dist_node = select_dest(G, s);
    sequence<EdgeTy> exp_dist;
    {
      sssp::Rho_Stepping solver(G);
      run([&]() { exp_dist = solver.sssp(s); }, solver);
    }
    ofstream ofs("ppsp.tsv", ios_base::app);
    ofs << s << '\n';
    ofs.close();
    printf("\n");
    size_t k = dist_node.size();
    for (size_t kth = k - 1; kth < k; kth *= 2) {
      NodeId t = dist_node[kth].second;
      printf("source %d: %-10d\tdestination %d: %-10d\n", v, s, v, t);
      sequence<EdgeTy> dists;
      EdgeTy dist;

      //{
        //ppsp::Rho_Stepping verifier(G, false);
        //run([&]() { dist = verifier.ppsp(s, t); }, verifier);
        //assert(check_equal(dist, exp_dist[t]));
        //dists.push_back(dist);
      //}

      //{
        //ppsp::Rho_Stepping verifier(G, true);
        //run([&]() { dist = verifier.ppsp(s, t); }, verifier);
        //assert(check_equal(dist, exp_dist[t]));
        //dists.push_back(dist);
      //}

      if (!coords.empty()) {
        {
          astar::Rho_Stepping verifier(G, coords, false);
          run([&]() { dist = verifier.astar(s, t); }, verifier);
          assert(check_equal(dist, exp_dist[t]));
          dists.push_back(dist);
        }

        {
          astar::Rho_Stepping verifier(G, coords, true);
          run([&]() { dist = verifier.astar(s, t); }, verifier);
          assert(check_equal(dist, exp_dist[t]));
          dists.push_back(dist);
        }
      }

      for (auto d : dists) {
        assert(check_equal(d, dist));
      }

      ofstream ofs("ppsp.tsv", ios_base::app);
      ofs << fixed << setprecision(6);
      ofs << t << '\t' << kth << '\t' << dist << '\n';
      printf("\n");
    }
  }
}
#else
NodeId select_dest(const Graph G, NodeId s, int percentile) {
  sssp::Rho_Stepping solver(G);
  size_t n = G.n;
  solver.reset();
  auto dist = solver.sssp(s);
  auto indices = pack_index(
      delayed_seq<bool>(n, [&](size_t i) { return dist[i] != DIST_MAX; }));
  auto dist_node = tabulate(indices.size(), [&](size_t i) {
    return make_pair(dist[indices[i]], indices[i]);
  });
  sort_inplace(make_slice(dist_node));
  size_t id = 1.0 * percentile / 100 * indices.size();
  return dist_node[id].second;
}

void run(const Graph &G, [[maybe_unused]] const sequence<double> &coords,
         const sequence<NodeId> &largest_cc) {
  size_t m = largest_cc.size();
  auto percentiles = vector<int>{1, 50, 99};
  for (auto percentile : percentiles) {
    for (int v = 0; v < NUM_SRC; v++) {
      NodeId s = largest_cc[hash32(v) % m];
      NodeId t = select_dest(G, s, percentile);
      printf("source %d: %-10d\tdestination %d: %-10d\n", v, s, v, t);
      sequence<EdgeTy> dists;
      EdgeTy dist;

      {
        sssp::Rho_Stepping solver(G);
        run([&]() { dist = solver.sssp(s)[t]; }, solver);
        dists.push_back(dist);
      }

      {
        ppsp::Rho_Stepping verifier(G, false);
        run([&]() { dist = verifier.ppsp(s, t); }, verifier);
        dists.push_back(dist);
      }

      {
        ppsp::Rho_Stepping verifier(G, true);
        run([&]() { dist = verifier.ppsp(s, t); }, verifier);
        dists.push_back(dist);
      }

      if (!coords.empty()) {
        {
          astar::Rho_Stepping verifier(G, coords, false);
          run([&]() { dist = verifier.astar(s, t); }, verifier);
          dists.push_back(dist);
        }

        {
          astar::Rho_Stepping verifier(G, coords, true);
          run([&]() { dist = verifier.astar(s, t); }, verifier);
          dists.push_back(dist);
        }
      }

      for (auto d : dists) {
        assert(check_equal(d, dist));
      }

      ofstream ofs("ppsp.tsv", ios_base::app);
      ofs << s << '\t' << t << '\t' << dist << '\n';
      printf("\n");
    }
  }
}
#endif

int main(int argc, char *argv[]) {
  if (argc == 1) {
    fprintf(stderr,
            "Usage: %s [-i input_file] [-c coordinates] [-w] [-s]\n"
            "Options:\n"
            "\t-i,\tinput file path\n"
            "\t-c,\tcoordinate file path\n"
            "\t-w,\tweighted input graph\n"
            "\t-s,\tsymmetrized input graph\n",
            argv[0]);
    exit(EXIT_FAILURE);
  }
  char c;
  bool weighted = false;
  bool symmetrized = false;
  char const *FILEPATH = nullptr;
  char const *COORDPATH = nullptr;
  while ((c = getopt(argc, argv, "i:c:ws")) != -1) {
    switch (c) {
      case 'i':
        FILEPATH = optarg;
        break;
      case 'c':
        COORDPATH = optarg;
        break;
      case 'w':
        weighted = true;
        break;
      case 's':
        symmetrized = true;
        break;
      default:
        fprintf(stderr, "Error: Unknown option %c\n", optopt);
        exit(EXIT_FAILURE);
    }
  }
  Graph G(weighted, symmetrized);

  printf("Reading graph...\n");
  G.read_graph(FILEPATH);
  if (!weighted) {
    printf("Generating edge weights...\n");
    G.generate_weight();
  }

  sequence<double> coords;
  if (COORDPATH != nullptr && strcmp(COORDPATH, "")) {
    coords = read_coords(COORDPATH, G.n);
    // assert(coords.size() == G.n);
  }

  fprintf(stdout, "Running on %s: |V|=%zu, |E|=%zu, num_src=%d, num_round=%d\n",
          FILEPATH, G.n, G.m, NUM_SRC, NUM_ROUND);

  auto largest_cc = get_largest_cc(G);

  run(G, coords, largest_cc);
  return 0;
}