text: // This code is part of the project "Theoretically Efficient Parallel Graph
  // Algorithms Can Be Fast and Scalable", presented at Symposium on Parallelism
  // in Algorithms and Architectures, 2018.
  // Copyright (c) 2018 Laxman Dhulipala, Guy Blelloch, and Julian Shun
  //
  // 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.
  
  #pragma once
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  #include <cmath>
  #include <fstream>
  #include <iostream>
  #include <tuple>
  
  #include "pbbslib/binary_search.h"
  #include "pbbslib/seq.h"
  #include "pbbslib/sequence_ops.h"
  #include "pbbslib/utilities.h"
  
  #include "ligra/macros.h"
  
  namespace bytepd {
  
  template <class W>
  inline std::tuple<uintE, W> get_ith_neighbor(uchar* edge_start, uintE source,
                                               uintE degree, size_t i) {}
  
  // Read default weight (expects pbbslib::empty)
  template <class W,
            typename std::enable_if<!std::is_same<W, intE>::value, int>::type = 0>
  inline W eatWeight(uchar*& start) {
    return (W)pbbslib::empty();
  }
  
  // Read integer weight
  template <class W,
            typename std::enable_if<std::is_same<W, intE>::value, int>::type = 0>
  inline W eatWeight(uchar*& start) {
    uchar fb = *start++;
    intE edgeRead = (fb & 0x3f);
    if (LAST_BIT_SET(fb)) {
      int shiftAmount = 6;
      while (1) {
        uchar b = *start;
        edgeRead |= ((b & 0x7f) << shiftAmount);
        start++;
        if (LAST_BIT_SET(b))
          shiftAmount += EDGE_SIZE_PER_BYTE;
        else
          break;
      }
    }
    return (fb & 0x40) ? -edgeRead : edgeRead;
  }
  
  inline intE eatFirstEdge(uchar*& start, uintE source) {
    uchar fb = *start++;
    intE edgeRead = (fb & 0x3f);
    if (LAST_BIT_SET(fb)) {
      int shiftAmount = 6;
      while (1) {
        uchar b = *start;
        edgeRead |= ((b & 0x7f) << shiftAmount);
        start++;
        if (LAST_BIT_SET(b))
          shiftAmount += EDGE_SIZE_PER_BYTE;
        else
          break;
      }
    }
    return (fb & 0x40) ? source - edgeRead : source + edgeRead;
  }
  
  /*
    Reads any edge of an out-edge list after the first edge.
  */
  inline uintE eatEdge(uchar*& start) {
    uintE edgeRead = 0;
    int shiftAmount = 0;
  
    while (1) {
      uchar b = *start;
      edgeRead += ((b & 0x7f) << shiftAmount);
      start++;
      if (LAST_BIT_SET(b))
        shiftAmount += EDGE_SIZE_PER_BYTE;
      else
        break;
    }
    return edgeRead;
  }
  
  template <class W>
  struct iter {
    uchar* base;
    uchar* finger;
    uintE src;
    uintT degree;
  
    uintE num_blocks;
    uintE cur_chunk;
  
    std::tuple<uintE, W> last_edge;
    uintE proc;
  
    iter(uchar* _base, uintT _degree, uintE _src)
        : base(_base), degree(_degree), src(_src) {
      num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      finger = base + (num_blocks - 1) * sizeof(uintE);
      cur_chunk = 0;
      if (degree == 0) {
        return;
      }
  
      std::get<0>(last_edge) = eatFirstEdge(finger, src);
      std::get<1>(last_edge) = eatWeight<W>(finger);
      proc = 1;
    }
  
    inline std::tuple<uintE, W> cur() { return last_edge; }
  
    inline std::tuple<uintE, W> next() {
      if (proc == PARALLEL_DEGREE) {
        cur_chunk++;
        std::get<0>(last_edge) = eatFirstEdge(finger, src);
        std::get<1>(last_edge) = eatWeight<W>(finger);
        proc = 1;
      } else {
        std::get<0>(last_edge) += eatEdge(finger);
        std::get<1>(last_edge) = eatWeight<W>(finger);
        proc++;
      }
      return last_edge;
    }
  
    inline bool has_next() {
      return (cur_chunk * PARALLEL_DEGREE + proc) < degree;
    }
  };
  
  template <class W>
  struct simple_iter {
    uchar* base;
    uchar* finger;
    uintE src;
    uintT degree;
  
    uintE num_blocks;
    uintE cur_chunk;
  
    std::tuple<uintE, W> last_edge;
    uintE proc;
  
    simple_iter(uchar* _base, uintT _degree, uintE _src)
        : base(_base), degree(_degree), src(_src) {
      num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      finger = base + (num_blocks - 1) * sizeof(uintE);
      cur_chunk = 0;
      if (degree == 0) {
        return;
      }
  
      std::get<0>(last_edge) = eatFirstEdge(finger, src);
      std::get<1>(last_edge) = eatWeight<W>(finger);
      proc = 1;
    }
  
    inline std::tuple<uintE, W> cur() { return last_edge; }
  
    inline std::tuple<uintE, W> next() {
      if (proc == PARALLEL_DEGREE) {
        cur_chunk++;
        std::get<0>(last_edge) = eatFirstEdge(finger, src);
        std::get<1>(last_edge) = eatWeight<W>(finger);
        proc = 1;
      } else {
        std::get<0>(last_edge) += eatEdge(finger);
        std::get<1>(last_edge) = eatWeight<W>(finger);
        proc++;
      }
      return last_edge;
    }
  
    inline bool has_next() {
      return (cur_chunk * PARALLEL_DEGREE + proc) < degree;
    }
  };
  
  // Decode unweighted edges
  template <
      class W, class T,
      typename std::enable_if<std::is_same<W, pbbslib::empty>::value, int>::type = 0>
  inline void decode(T t, uchar* edge_start, const uintE& source,
                     const uintT& degree, const bool par = true) {
    if (degree > 0) {
      size_t num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)
          edge_start;  // use beginning of edgeArray for offsets into edge list
      uchar* finger = edge_start + (num_blocks - 1) * sizeof(uintE);
      // do first chunk
      size_t end = std::min<long>(PARALLEL_DEGREE, degree);
  
      // Eat first edge, which is compressed specially
      uintE ngh = eatFirstEdge(finger, source);
      if (!t(source, ngh, pbbslib::empty(), 0)) return;
      for (uintE edgeID = 1; edgeID < end; edgeID++) {
        // Eat the next 'edge', which is a difference, and reconstruct edge.
        ngh += eatEdge(finger);
        if (!t(source, ngh, pbbslib::empty(), edgeID)) return;
      }
      // do remaining chunks in parallel
      par_for(1, num_blocks, 1, [&] (size_t i) {
        size_t o = i * PARALLEL_DEGREE;
        size_t end = std::min<long>(o + PARALLEL_DEGREE, degree);
        uchar* finger = edge_start + block_offsets[i - 1];
        // Eat first edge, which is compressed specially
        uintE ngh = eatFirstEdge(finger, source);
        if (!t(source, ngh, pbbslib::empty(), o)) end = 0;
        for (size_t edgeID = o + 1; edgeID < end; edgeID++) {
          // Eat the next 'edge', which is a difference, and reconstruct edge.
          ngh += eatEdge(finger);
          if (!t(source, ngh, pbbslib::empty(), edgeID)) break;
        }
      }, par);
    }
  }
  
  // Decode weighted edges
  template <class W, class T,
            typename std::enable_if<!std::is_same<W, pbbslib::empty>::value,
                                    int>::type = 0>
  inline void decode(T t, uchar* edge_start, const uintE& source,
                     const uintT& degree, const bool par = true) {
    if (degree > 0) {
      size_t num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)
          edge_start;  // use beginning of edgeArray for offsets into edge list
      uchar* finger = edge_start + (num_blocks - 1) * sizeof(uintE);
      // do first chunk
      size_t end = std::min<long>(PARALLEL_DEGREE, degree);
  
      // Eat first edge, which is compressed specially
      uintE ngh = eatFirstEdge(finger, source);
      W weight = eatWeight<W>(finger);
      if (!t(source, ngh, weight, 0)) return;
      for (size_t edgeID = 1; edgeID < end; edgeID++) {
        // Eat the next 'edge', which is a difference, and reconstruct edge.
        ngh += eatEdge(finger);
        W weight = eatWeight<W>(finger);
        if (!t(source, ngh, weight, edgeID)) return;
      }
      // do remaining chunks in parallel
      par_for(1, num_blocks, 1, [&] (size_t i) {
        size_t o = i * PARALLEL_DEGREE;
        size_t end = std::min<long>(o + PARALLEL_DEGREE, degree);
        uchar* finger = edge_start + block_offsets[i - 1];
        // Eat first edge, which is compressed specially
        uintE ngh = eatFirstEdge(finger, source);
        W weight = eatWeight<W>(finger);
        if (!t(source, ngh, weight, o)) end = 0;
        for (size_t edgeID = o + 1; edgeID < end; edgeID++) {
          // Eat the next 'edge', which is a difference, and reconstruct edge.
          ngh += eatEdge(finger);
          W weight = eatWeight<W>(finger);
          if (!t(source, ngh, weight, edgeID)) break;
        }
      }, par);
    }
  }
  
  // Map_reduce
  template <class W, class E, class M, class R>
  inline E map_reduce(uchar* edge_start, const uintE& source, const uintT& degree,
                      E id, M& m, R& r, const bool par = true) {
    if (degree > 0) {
      size_t num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)edge_start;
  
      E stk[1000];
      E* block_outputs;
      if (num_blocks > 1000) {
        block_outputs = pbbslib::new_array_no_init<E>(num_blocks);
      } else {
        block_outputs = (E*)stk;
      }
  
      par_for(0, num_blocks, 2, [&] (size_t i) {
        size_t start = i * PARALLEL_DEGREE;
        size_t end = std::min<long>(start + PARALLEL_DEGREE, degree);
        uchar* finger = edge_start + ((i == 0) ? (num_blocks - 1) * sizeof(uintE)
                                               : block_offsets[i - 1]);
  
        // Eat first edge, which is compressed specially
        uintE ngh = eatFirstEdge(finger, source);
        W wgh = eatWeight<W>(finger);
        E cur = m(source, ngh, wgh);
        for (size_t j = start + 1; j < end; j++) {
          ngh += eatEdge(finger);
          W wgh = eatWeight<W>(finger);
          cur = r(cur, m(source, ngh, wgh));
        }
        block_outputs[i] = cur;
      }, par);
  
      auto im = pbbslib::make_sequence(block_outputs, num_blocks);
      E res = pbbslib::reduce(im, r);
      if (num_blocks > 1000) {
        pbbslib::free_array(block_outputs);
      }
      return res;
    } else {
      return id;
    }
  }
  
  template <class W>
  uchar* decode_block(uchar* finger, std::tuple<uintE, W>* out, size_t start,
                      size_t end, const uintE& source) {
    uintE ngh = eatFirstEdge(finger, source);
    W wgh = eatWeight<W>(finger);
    out[start] = std::make_tuple(ngh, wgh);
    for (size_t i = start + 1; i < end; i++) {
      // Eat the next 'edge', which is a difference, and reconstruct edge.
      ngh += eatEdge(finger);
      W wgh = eatWeight<W>(finger);
      out[i] = std::make_tuple(ngh, wgh);
    }
    return finger;
  }
  
  template <class W, class T>
  inline void decode_block_seq(T t, uchar* edge_start, const uintE& source,
                               const uintT& degree, uintE block_size,
                               uintE block_num) {
    assert(false);
  }
  
  inline size_t get_virtual_degree(uintE d, uchar* ngh_arr) { return d; }
  
  #define SEQ_THRESH 10
  
  // Represents the sequence from [...[|cur_block|...|cur_block+num_blocks|]...]
  struct seq_info {
    uchar* edge_start;
    uintE degree;
    uintE source_id;
    uintE total_blocks;
  
    uintE start;
    uintE end;
  
    seq_info(uchar* es, uintE d, uintE sid, uintE tb, uintE s, uintE e)
        : edge_start(es),
          degree(d),
          source_id(sid),
          total_blocks(tb),
          start(s),
          end(e) {}
  
    uchar* get_start_of_block(const uintE& block_id) {
      if (total_blocks == 1) {
        return edge_start;
      }
      uintE* offs = (uintE*)edge_start;
      return edge_start + offs[block_id - 1];
    }
  
    uintE get_pivot() {
      uintE* offs = (uintE*)edge_start;
      uintE mid_block = (end + start) / 2;
      uchar* finger = edge_start + offs[mid_block];
      return eatFirstEdge(finger, source_id);
    }
  
    uintE pivot_block() { return (end + start) / 2; }
  
    uintE binary_search(uintE pivot) {
      uintE* offs = (uintE*)edge_start;
      auto start_f = [&](size_t i) {
        uchar* finger = edge_start + offs[start + i];
        return eatFirstEdge(finger, source_id);
      };
      auto start_im = pbbslib::make_sequence<uintE>(size(), start_f);
      uintE ind =
          pbbslib::binary_search(start_im, pivot, std::greater<uintE>());  // check
      // ind is the first block index (from start) <= our pivot.
      //    decode_block<
  
      return ind;
    }
  
    seq_info cut(uintE l, uintE r) {
      return seq_info(edge_start, degree, source_id, total_blocks, l, r);
    }
  
    uintE size() { return end - start; }
  };
  
  inline uintE seq_intersect_full(seq_info u, seq_info v) {
    assert(false);  // not implemented
    return 0;
  }
  
  inline uintE seq_intersect(seq_info u, seq_info v) {
    assert(false);  // not implemented
    return 0;
    //  decode_block<pbbslib::empty>(finger, (std::tuple<uintE, pbbslib::empty>*)ngh_u,
    //  0,
  }
  
  inline uintE intersect(seq_info u, seq_info v) {
    // Might need to swap here
    uintE nA = u.size();
    uintE nB = v.size();
    if (nA + nB < SEQ_THRESH) {
      return seq_intersect_full(u, v);
    } else if (nA == 1) {  // merge base
      return seq_intersect(u, v);
    } else if (nB == 0) {
      return 0;
    } else {  // (large, large)
      uintE pivot = u.get_pivot();
      uintE mU = u.pivot_block();
      uintE mV = v.binary_search(pivot);
      uintE lA = 0, rA = 0;
      par_do(true, [&]() { lA = intersect(u.cut(0, mU), v.cut(0, mV)); },
             [&]() { rA = intersect(u.cut(mU, mU), v.cut(0, mV)); });
      return lA + rA;
    }
  }
  
  /*
    Compresses the first edge, writing target-source and a sign bit.
  */
  long compressFirstEdge(uchar* start, long offset, uintE source, uintE target) {
    intE preCompress = (intE)target - source;
    int bytesUsed = 0;
    uchar firstByte = 0;
    intE toCompress = std::abs(preCompress);
    firstByte = toCompress & 0x3f;  // 0011|1111
    if (preCompress < 0) {
      firstByte |= 0x40;
    }
    toCompress = toCompress >> 6;
    if (toCompress > 0) {
      firstByte |= 0x80;
    }
    start[offset] = firstByte;
    offset++;
  
    uchar curByte = toCompress & 0x7f;
    while ((curByte > 0) || (toCompress > 0)) {
      bytesUsed++;
      uchar toWrite = curByte;
      toCompress = toCompress >> 7;
      // Check to see if there's any bits left to represent
      curByte = toCompress & 0x7f;
      if (toCompress > 0) {
        toWrite |= 0x80;
      }
      start[offset] = toWrite;
      offset++;
    }
    return offset;
  }
  
  template <class W,
            typename std::enable_if<!std::is_same<W, intE>::value, int>::type = 0>
  inline long compressWeight(uchar* start, long offset, W weight) {
    return offset;
  }
  
  template <class W,
            typename std::enable_if<std::is_same<W, intE>::value, int>::type = 0>
  inline long compressWeight(uchar* start, long offset, W weight) {
    return compressFirstEdge(start, offset, 0, weight);
  }
  
  /*
    Should provide the difference between this edge and the previous edge
  */
  
  long compressEdge(uchar* start, long curOffset, uintE e) {
    uchar curByte = e & 0x7f;
    int bytesUsed = 0;
    while ((curByte > 0) || (e > 0)) {
      bytesUsed++;
      uchar toWrite = curByte;
      e = e >> 7;
      // Check to see if there's any bits left to represent
      curByte = e & 0x7f;
      if (e > 0) {
        toWrite |= 0x80;
      }
      start[curOffset] = toWrite;
      curOffset++;
    }
    return curOffset;
  }
  
  template <class W>
  size_t compute_block_size(std::tuple<uintE, W>* edges, size_t start, size_t end,
                            const uintE& source) {
    uchar tmp[16];
    size_t total_bytes = 0;
    auto nw = edges[start];
    total_bytes += compressFirstEdge(tmp, 0, source, std::get<0>(nw));
    total_bytes += compressWeight<W>(tmp, 0, std::get<1>(nw));
    uintE last_ngh = std::get<0>(nw);
    for (size_t i = start + 1; i < end; i++) {
      nw = edges[i];
      uintE difference = std::get<0>(nw) - last_ngh;
      total_bytes += compressEdge(tmp, 0, difference);
      total_bytes += compressWeight<W>(tmp, 0, std::get<1>(nw));
      last_ngh = std::get<0>(nw);
    }
    return total_bytes;
  }
  
  template <class W>
  size_t compute_size_in_bytes(std::tuple<uintE, W>* edges, const uintE& source,
                               const uintE& d) {
    if (d > 0) {
      size_t num_blocks = 1 + (d - 1) / PARALLEL_DEGREE;
      uintE stk[100];
      uintE* block_bytes;
      if (num_blocks > 100) {
        block_bytes = pbbslib::new_array_no_init<uintE>(num_blocks);
      } else {
        block_bytes = (uintE*)stk;
      }
      par_for(0, num_blocks, 10, [&] (size_t i) {
        // calculate size in bytes of this block
        uintE start = i * PARALLEL_DEGREE;
        uintE end = start + std::min<uintE>(PARALLEL_DEGREE, d - start);
        block_bytes[i] = compute_block_size(edges, start, end, source);
      });
      auto bytes_imap = pbbslib::make_sequence(block_bytes, num_blocks);
      size_t total_space = pbbslib::scan_add_inplace(bytes_imap);
  
      // add in space for storing offsets to the start of each block
      total_space += sizeof(uintE) * (num_blocks - 1);
      if (num_blocks > 100) {
        pbbslib::free_array(block_bytes);
      }
      return total_space;
    } else {
      return 0;
    }
  }
  
  template <class W, class I>
  long sequentialCompressEdgeSet(uchar* edgeArray, size_t current_offset,
                                 uintT degree, uintE source, I& it) {
    if (degree > 0) {
      size_t start_offset = current_offset;
      size_t num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)edgeArray;
      current_offset += (num_blocks - 1) * sizeof(uintE);  // block_offs
      for (long i = 0; i < num_blocks; i++) {
        size_t o = i * PARALLEL_DEGREE;
        size_t end = std::min<size_t>(PARALLEL_DEGREE, degree - o);
  
        if (i > 0)
          block_offsets[i - 1] =
              current_offset -
              start_offset;  // store offset for all chunks but the first
  
        std::tuple<uintE, W> lst = (i == 0) ? it.cur() : it.next();
        uintE last_ngh = std::get<0>(lst);
        current_offset =
            compressFirstEdge(edgeArray, current_offset, source, last_ngh);
        current_offset =
            compressWeight<W>(edgeArray, current_offset, std::get<1>(lst));
        for (size_t edgeI = 1; edgeI < end; edgeI++) {
          std::tuple<uintE, W> nxt = it.next();
          uintE difference = std::get<0>(nxt) - last_ngh;
          current_offset = compressEdge(edgeArray, current_offset, difference);
          current_offset =
              compressWeight<W>(edgeArray, current_offset, std::get<1>(nxt));
          last_ngh = std::get<0>(nxt);
        }
      }
    }
    return current_offset;
  }
  
  /*
    Takes:
      1. The edge array of chars to write into
      2. The current offset into this array
      3. The vertices degree
      4. The vertices vertex number
      5. The array of saved out-edges we're compressing
    Returns:
      The new offset into the edge array
  */
  long sequentialCompressEdgeSet(uchar* edgeArray, size_t current_offset,
                                 uintT degree, uintE source, uintE* edges) {
    if (degree > 0) {
      long start_offset = current_offset;
      long num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)
          edgeArray;  // use beginning of edgeArray for offsets into edge list
      current_offset += (num_blocks - 1) * sizeof(uintE);
      for (long i = 0; i < num_blocks; i++) {
        long o = i * PARALLEL_DEGREE;
        long end = std::min<long>(PARALLEL_DEGREE, degree - o);
        uintE* edges_i = edges + o;
        if (i > 0)
          block_offsets[i - 1] =
              current_offset -
              start_offset;  // store offset for all chunks but the first
        // Compress the first edge whole, which is signed difference coded
        current_offset =
            compressFirstEdge(edgeArray, current_offset, source, edges_i[0]);
        for (uintT edgeI = 1; edgeI < end; edgeI++) {
          // Store difference between cur and prev edge.
          uintE difference = edges_i[edgeI] - edges_i[edgeI - 1];
          current_offset = compressEdge(edgeArray, current_offset, difference);
        }
      }
    }
    return current_offset;
  }
  
  // Packs out the block sequentially.
  template <class W, class P>
  size_t packBlock(uchar* edge_start, const uintE& source, const uintE& degree,
                   P& pred) {
    uchar* finger = edge_start;
  
    //  auto map_f = [&] (const uintE& src, const uintE& ngh, const W& w) {
    //    assert(src < 3072627);
    //    assert(ngh < 3072627);
    //    return 1;
    //  };
    //  auto red_f = [&] (size_t l, size_t r) { return l + r; };
    //  map_reduce<W>(edge_start, source, degree, 0, map_f, red_f, true);
    //
    size_t ct = 0;
  
    std::tuple<uintE, W> tmp[PARALLEL_DEGREE];
    uintE ngh = eatFirstEdge(finger, source);
    W wgh = eatWeight<W>(finger);
    if (pred(source, ngh, wgh)) {
      tmp[ct++] = std::make_tuple(ngh, wgh);
    }
    for (size_t i = 1; i < degree; i++) {
      ngh += eatEdge(finger);
      wgh = eatWeight<W>(finger);
      if (pred(source, ngh, wgh)) {
        tmp[ct++] = std::make_tuple(ngh, wgh);
      }
    }
  
    if (ct == 0) {
      return 0;
    }
  
    size_t offset = 0;  // the write offset
    offset = compressFirstEdge(edge_start, offset, source, std::get<0>(tmp[0]));
    offset = compressWeight<W>(edge_start, offset, std::get<1>(tmp[0]));
    uintE last_ngh = std::get<0>(tmp[0]);
    for (size_t i = 1; i < ct; i++) {
      const auto& e = tmp[i];
      uintE difference = std::get<0>(e) - last_ngh;
      offset = compressEdge(edge_start, offset, difference);
      offset = compressWeight<W>(edge_start, offset, std::get<1>(e));
      last_ngh = std::get<0>(e);
    }
    assert(offset <= (finger - edge_start));
  
    //  size_t ctt = map_reduce<W>(edge_start, source, ct, 0, map_f, red_f, true);
  
    return ct;
  }
  
  template <class W, class P, class O>
  inline size_t filterBlock(P pred, uchar* edge_start, const uintE& source,
                            const uintE& degree, O& out) {
    size_t num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
    size_t n_offsets = num_blocks - 1;
    uintE edges_offset = (n_offsets) * sizeof(uintE);
    uintE* block_offsets = (uintE*)edge_start;
    size_t k = 0;
    for (size_t i = 0; i < num_blocks; i++) {
      uchar* finger =
          edge_start + ((i == 0) ? edges_offset : block_offsets[i - 1]);
      uintE start = i * PARALLEL_DEGREE;
      uintE end = std::min<uintE>(start + PARALLEL_DEGREE, degree);
  
      uintE ngh = eatFirstEdge(finger, source);
      W wgh = eatWeight<W>(finger);
      if (pred(source, ngh, wgh)) {
        out(k++, std::make_tuple(ngh, wgh));
      }
      for (size_t j = start + 1; j < end; j++) {
        ngh += eatEdge(finger);
        W wgh = eatWeight<W>(finger);
        if (pred(source, ngh, wgh)) {
          out(k++, std::make_tuple(ngh, wgh));
        }
      }
    }
    return k;
  }
  
  template <class W>
  void compress_edges(uchar* edgeArray, const uintE& source, const uintE& d,
                      std::tuple<uintE, W>* edges, uchar* last_finger,
                      bool par = true) {
    size_t num_blocks = 1 + (d - 1) / PARALLEL_DEGREE;
    uintE stk[101];
    uintE* block_bytes;
    if (num_blocks > 100) {
      block_bytes = pbbslib::new_array_no_init<uintE>(num_blocks + 1);
    } else {
      block_bytes = (uintE*)stk;
    }
  
    par_for(0, num_blocks, 1, [&] (size_t i) {
      uintE start = i * PARALLEL_DEGREE;
      uintE end = std::min<uintE>(start + PARALLEL_DEGREE, d);
      block_bytes[i] = compute_block_size(edges, start, end, source);
    });
    block_bytes[num_blocks] = 0;
  
    uintE edges_offset = (num_blocks - 1) * sizeof(uintE);
    uintE* block_offsets = (uintE*)edgeArray;
  
    auto bytes_imap = pbbslib::make_sequence(block_bytes, num_blocks + 1);
    uintE total_space = pbbslib::scan_add_inplace(bytes_imap);
  
    if (total_space > (last_finger - edgeArray)) {
      std::cout << "Space error!"
                << "\n";
      exit(0);
    }
    if (total_space == (last_finger - edgeArray)) {
      std::cout << "d = " << d << " to exactly the same space"
                << "\n";
    }
  
    par_for(1, num_blocks, 1, [&] (size_t i) {
      // store offset to start of this block from edgeArray
      block_offsets[i - 1] = edges_offset + bytes_imap[i];
    });
  
    par_for(0, num_blocks, 1, [&] (size_t i) {
      uintE start = i * PARALLEL_DEGREE;
      uintE end = std::min<uintE>(start + PARALLEL_DEGREE, d);
      long write_offset = (i == 0) ? edges_offset : block_offsets[i - 1];
      auto nw = edges[start];
      write_offset =
          compressFirstEdge(edgeArray, write_offset, source, std::get<0>(nw));
      write_offset = compressWeight<W>(edgeArray, write_offset, std::get<1>(nw));
      uintE last_edge = std::get<0>(nw);
      for (size_t i = start + 1; i < end; i++) {
        nw = edges[i];
        write_offset =
            compressEdge(edgeArray, write_offset, std::get<0>(nw) - last_edge);
        write_offset =
            compressWeight<W>(edgeArray, write_offset, std::get<1>(nw));
        last_edge = std::get<0>(nw);
      }
    });
    if (num_blocks > 100) {
      pbbslib::free_array(block_bytes);
    }
  }
  
  template <class W, class P>
  inline size_t pack(P& pred, uchar* edge_start, const uintE& source,
                     const uintE& degree, std::tuple<uintE, W>* tmp) {
    if (degree <= PARALLEL_DEGREE) {
      // Pack small degree vertices sequentially
      // return packSequential<W, P>(edge_start, source, degree, pred);
      return packBlock<W, P>(edge_start, source, degree, pred);
    } else {
      return degree;
  
      // Pack large degree vertices in parallel
      // Decode nghs into temporary space, filter, then recompress.
      size_t num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)edge_start;
  
      std::tuple<uintE, W> tmp_A[200000];
      std::tuple<uintE, W>* our_tmp = (std::tuple<uintE, W>*)tmp_A;
      if (degree > PD_PACK_THRESHOLD) {
        our_tmp = tmp;
      }
  
      //    auto map_f = [&] (const uintE& src, const uintE& ngh, const W& w) {
      //      assert(src < 3072627);
      //      assert(ngh < 3072627);
      //      return 1;
      //    };
      //    auto red_f = [&] (size_t l, size_t r) { return l + r; };
      //    size_t ct = map_reduce<W>(edge_start, source, degree, 0, map_f, red_f,
      //    true);
  
      uchar* mf = 0;
      par_for(0, num_blocks, 1, [&] (size_t i) {
        size_t start = i * PARALLEL_DEGREE;
        size_t end = std::min<long>(start + PARALLEL_DEGREE, degree);
        uchar* finger = edge_start + ((i == 0) ? (num_blocks - 1) * sizeof(uintE)
                                               : block_offsets[i - 1]);
        uchar* last_finger = decode_block(finger, our_tmp, start, end, source);
        if (i == num_blocks - 1) {
          mf = last_finger;
        }
      });
  
      // filter edges into tmp2
      std::tuple<uintE, W>* tmp2 = our_tmp + degree;
      auto pd = [&](const std::tuple<uintE, W>& nw) {
        return pred(source, std::get<0>(nw), std::get<1>(nw));
      };
      uintE k = pbbslib::filterf(our_tmp, tmp2, degree, pd);
      if (k == degree || k == 0) {
        return k;
      }
  
      for (size_t i = 0; i < k; i++) {
        assert(std::get<0>(tmp2[i]) < 3072627);
      }
  
      // pack blocks out in parallel
      compress_edges<W>(edge_start, source, k, tmp2, mf);
  
      //    ct = map_reduce<W>(edge_start, source, k, 0, map_f, red_f, true);
      //    assert(ct == k);
  
      return k;
    }
  }
  
  template <class W, class P, class O>
  inline size_t filter(P pred, uchar* edge_start, const uintE& source,
                       const uintE& degree, std::tuple<uintE, W>* tmp, O& out) {
    return filterBlock<W>(pred, edge_start, source, degree, out);
    if (degree <= PD_PACK_THRESHOLD && degree > 0) {
      // Filter small degree vertices sequentially
      uchar* finger = edge_start;
      uintE ngh = eatFirstEdge(finger, source);
      W wgh = eatWeight<W>(finger);
      uintE k = 0;
      if (pred(source, ngh, wgh)) {
        out(k++, std::make_tuple(ngh, wgh));
      }
      for (size_t i = 1; i < degree; i++) {
        ngh += eatEdge(finger);
        W wgh = eatWeight<W>(finger);
        if (pred(source, ngh, wgh)) {
          out(k++, std::make_tuple(ngh, wgh));
        }
      }
      return k;
    } else {
      // Pack large degree vertices in parallel
      // Decode nghs into temporary space, then filter
      size_t num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)edge_start;
  
      par_for(0, num_blocks, 1, [&] (size_t i) {
        size_t start = i * PARALLEL_DEGREE;
        size_t end = std::min<long>(start + PARALLEL_DEGREE, degree);
        uchar* finger = edge_start + ((i == 0) ? (num_blocks - 1) * sizeof(uintE)
                                               : block_offsets[i - 1]);
        decode_block(finger, tmp, start, end, source);
      });
  
      // filter edges into tmp2
      std::tuple<uintE, W>* tmp2 = tmp + degree;
      auto pd = [&](const std::tuple<uintE, W>& nw) {
        return pred(source, std::get<0>(nw), std::get<1>(nw));
      };
      uintE k = pbbslib::filterf(tmp, tmp2, degree, pd);
      par_for(0, k, 2000, [&] (size_t i) { out(i, tmp2[i]); });
      return k;
    }
    return 0;
  }
  
  /*
    Compresses the edge set in parallel.
  */
  uintE* parallelCompressEdges(uintE* edges, uintT* offsets, long n, long m,
                               uintE* Degrees) {
    std::cout << "parallel compressing, (n,m) = (" << n << "," << m << ")"
              << "\n";
    uintE** edgePts = pbbslib::new_array_no_init<uintE*>(n);
    long* charsUsedArr = pbbslib::new_array_no_init<long>(n+1);
    auto charsUsed = pbbslib::make_sequence(charsUsedArr, n+1);
  //  long* compressionStarts = pbbslib::new_array_no_init<long>(n + 1);
    {
      par_for(0, n, pbbslib::kSequentialForThreshold, [&] (size_t i)
                      { charsUsed[i] = ceil((Degrees[i] * 9) / 8) + 4; });
    }
    charsUsed[n] = 0;
    long toAlloc = pbbslib::scan_inplace(charsUsed, pbbslib::addm<long>());
    uintE* iEdges = pbbslib::new_array_no_init<uintE>(toAlloc);
  
    {
      par_for(0, n, [&] (size_t i) {
        edgePts[i] = iEdges + charsUsed[i];
        long charsUsed =
            sequentialCompressEdgeSet((uchar*)(iEdges + charsUsed[i]), 0,
                                      Degrees[i], i, edges + offsets[i]);
        charsUsed[i] = charsUsed;
      });
    }
  
    // produce the total space needed for all compressed lists in chars.
    auto compressionStarts = pbbslib::scan(charsUsedArr, pbbslib::addm<long>());
    size_t totalSpace = compressionStarts[n];
    pbbslib::free_array(charsUsedArr);
  
    uchar* finalArr = pbbslib::new_array_no_init<uchar>(totalSpace);
    std::cout << "total space requested is : " << totalSpace << "\n";
    float avgBitsPerEdge = (float)totalSpace * 8 / (float)m;
    std::cout << "Average bits per edge: " << avgBitsPerEdge << "\n";
  
    {
      par_for(0, n, [&] (size_t i) {
        long o = compressionStarts[i];
        memcpy(finalArr + o, (uchar*)(edgePts[i]), compressionStarts[i + 1] - o);
        offsets[i] = o;
      });
    }
    offsets[n] = totalSpace;
    pbbslib::free_array(iEdges);
    pbbslib::free_array(edgePts);
    std::cout << "finished compressing, bytes used = " << totalSpace << "\n";
    std::cout << "would have been, " << (m * 4) << "\n";
    return ((uintE*)finalArr);
  }
  
  long sequentialCompressWeightedEdgeSet(uchar* edgeArray, long current_offset,
                                         const uintE& degree, const uintE& source,
                                         std::tuple<uintE, intE>* edges) {
    if (degree > 0) {
      uchar* our_start = edgeArray + current_offset;
      long start_offset = current_offset;
      long num_blocks = 1 + (degree - 1) / PARALLEL_DEGREE;
      uintE* block_offsets = (uintE*)our_start;
      current_offset += (num_blocks - 1) * sizeof(uintE);
  
      for (long i = 0; i < num_blocks; i++) {
        long start = i * PARALLEL_DEGREE;
        long end = start + std::min<long>(PARALLEL_DEGREE, degree - start);
        if (i > 0) block_offsets[i - 1] = current_offset - start_offset;
  
        auto nw = edges[start];
        current_offset =
            compressFirstEdge(edgeArray, current_offset, source, std::get<0>(nw));
        current_offset =
            compressWeight<intE>(edgeArray, current_offset, std::get<1>(nw));
  
        uintE last_ngh = std::get<0>(nw);
        for (size_t i = start + 1; i < end; i++) {
          auto nw = edges[i];
          uintE difference = std::get<0>(nw) - last_ngh;
          current_offset = compressEdge(edgeArray, current_offset, difference);
          current_offset =
              compressWeight<intE>(edgeArray, current_offset, std::get<1>(nw));
          last_ngh = std::get<0>(nw);
        }
      }
    }
    return current_offset;
  }
  
  uchar* parallelCompressWeightedEdges(std::tuple<uintE, intE>* edges,
                                       uintT* offsets, long n, long m,
                                       uintE* Degrees) {
    std::cout << "parallel compressing, (n,m) = (" << n << "," << m << ")"
              << "\n";
    auto bytes_used = pbbslib::new_array_no_init<size_t>(n + 1);
    auto bytes_seq = pbbslib::make_sequence<size_t>(bytes_used, n+1);
  
    par_for(0, n, [&] (size_t i) {
      bytes_used[i] = compute_size_in_bytes(edges + offsets[i], i, Degrees[i]);
    });
    bytes_used[n] = 0;
    size_t total_bytes = pbbslib::scan_inplace(bytes_used, pbbslib::addm<size_t>());
  
    uchar* edges_c = pbbslib::new_array_no_init<uchar>(total_bytes);
  
    par_for(0, n, [&] (size_t i) {
      sequentialCompressWeightedEdgeSet(edges_c, bytes_used[i], Degrees[i], i,
                                        edges + offsets[i]);
    });
  
    float avgBitsPerEdge = (float)total_bytes * 8 / (float)m;
    std::cout << "Average bits per edge: " << avgBitsPerEdge << "\n";
  
    size_t end = n+1;
    par_for(0, end, pbbslib::kSequentialForThreshold, [&] (size_t i)
                    { offsets[i] = bytes_used[i]; });
    std::cout << "finished compressing, bytes used = " << total_bytes << "\n";
    std::cout << "would have been, " << (m * 8) << "\n";
    return edges_c;
  }
  
  };  // namespace bytepd
  