Finding existing circles in a set of edges in an efficient manner

Question

Summary

Given a set of edges (anywhere from 3 to 1,000,000), assemble closed loops (I call them circles for ease) efficently. I am running this as a part of a larger Python project so I expect that the best solution will be written in C++/CUDA with python bindings (which is what I've tried).

---

The problem

Assemble rings, given a set of edges (defined by two vertex indices) that satisfy the following:

• there are no stray edges (all edges are used to create circles),
• all rings close (ex. if there is an edge 1-2 then there is an edge n-1 that will close it),
• edges have unordered indices (ex. an edge can be 1-2 or 2-1).

---

My general approach

• pick the edge at the back of a list and set one vertex to the start of the new ring (pStart) and the other vertex to the next point in the chain (pEnd) adding both to a new ring list,
• search for pEnd in the edges list,
• once edge is found, update pEnd to the vertex that is not equal to pEnd and add this to the ring list,
• repeat the above two steps until pStart==pEnd,
• if no more edges stop, else repeat above.

---

My Implementation

In C++, I implemented a serial and parallel version. I tested it using a set of 45,000 edges and obtain the following results:

• serial (105 seconds)
• parallel - CUDA Thrust (28 seconds)

Serial:

#include <algorithm>
#include <vector>
#include <string>
#include <boost/algorithm/string.hpp>
#include <fstream>
#include <iostream>
#include <chrono>

std::vector< std::vector<int> > rings_from_edges(std::vector<std::vector<int>> edges)
{
    int pStart, pEnd;

    std::vector<int> temp;
    std::vector<std::vector<int>> rings;

    temp = edges.back();
    edges.pop_back();

    pStart = temp[0];
    pEnd = temp[1];

    int p1, p2;

    while(not edges.empty())
        // Scan edges list until pEnd is found.
        for(auto const& pts: edges)
        {
            p1 = pts[0];
            p2 = pts[1];

            // Check if the start of the edge corresponds with the end of the ring.
            if(p1 == pEnd)
            {
                temp.push_back(p2);
                pEnd = p2;
                edges.erase(std::remove(edges.begin(), edges.end(), pts), edges.end());

                // Check if the beginning of the ring is the same as the end of the newly appended edge.
                if (p2 == pStart)
                {
                    // Add the newly created ring to the rings list.
                    rings.push_back(temp);
                    temp.clear();

                    // If the edges list contains more edges, reset the starting and end points to search for a new ring.
                    if(not edges.empty())
                    {
                        temp = edges.back();
                        edges.pop_back();

                        pStart = temp[0];
                        pEnd = temp[1];
                    }
                }
                break;
            }
                // Check if the end of the edge corresponds with the end of the ring.
            else if(p2 == pEnd)
            {
                temp.push_back(p1);
                pEnd = p1;
                edges.erase(std::remove(edges.begin(), edges.end(), pts), edges.end());

                // Check if the beginning of the ring is the same as the end of the newly appended edge.
                if (p1 == pStart)
                {
                    // Add the newly created ring to the rings list.
                    rings.push_back(temp);
                    temp.clear();

                    // If the edges list contains more edges, reset the starting and end points to search for a new ring.
                    if(not edges.empty())
                    {
                        temp = edges.back();
                        edges.pop_back();

                        pStart = temp[0];
                        pEnd = temp[1];
                    }
                }
                break;
            }
        }
    return rings;
}

int main() {

    std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

    std::vector< std::vector<int> > vectIN, vectOUT;


    std::string fileName = "PATH TO CSV FILE";

    std::string delimeter = ",";

    std::ifstream file(fileName);

    std::string line = "";

    while (getline(file, line))
    {
        std::vector<std::string> vec;
        boost::algorithm::split(vec, line, boost::is_any_of(delimeter));
        std::vector<int> vec2;
        vec2.emplace_back(std::stoi(vec.data()[0]));
        vec2.emplace_back(std::stoi(vec.data()[1]));

        vectIN.push_back(vec2);
    }

    file.close();

    std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();

    vectOUT = rings_from_edges(vectIN);

    std::chrono::steady_clock::time_point t3 = std::chrono::steady_clock::now();

    for (auto const& ring:vectOUT)
    {
        for(auto const& pt:ring)
        {
            if(pt>=0)
                std::cout << pt << " ";
        }
        std::cout << std::endl;
    }

    std::chrono::steady_clock::time_point t4 = std::chrono::steady_clock::now();

    long t1_t2 = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
    long t2_t3 = std::chrono::duration_cast<std::chrono::milliseconds>(t3 - t2).count();
    long t3_t4 = std::chrono::duration_cast<std::chrono::milliseconds>(t4 - t3).count();

    std::cout << "Load csv:      " << t1_t2 << std::endl;
    std::cout << "Ring assembly: " << t2_t3 << std::endl;
    std::cout << "Output:        " << t3_t4 << std::endl;

    std::cout << "----------------- THAT'S ALL FOLKS!!! -----------------" << std::endl;

    return 0;
}

CUDA Thrust version of above:

#include <thrust/remove.h>
#include <thrust/device_vector.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/tuple.h>
#include <thrust/copy.h>
#include <thrust/sort.h>
#include <thrust/find.h>
#include <thrust/functional.h>

#include <iostream>
#include <string>
#include <vector>
#include <fstream>

#include <iterator>
#include <algorithm>
#include <boost/algorithm/string.hpp>
#include <chrono>
#include <thrust/binary_search.h>

#include <thrust/uninitialized_copy.h>
#include <thrust/device_malloc.h>
#include <thrust/device_vector.h>


int main(){
    std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

    std::string fileName = "PATH TO CSV HERE";

    std::string delimeter = ",";

    std::ifstream file(fileName);

    std::vector<std::vector<int>> vectIN;

    std::string line = "";

    while (getline(file, line))
    {
        std::vector<std::string> vec;
        boost::algorithm::split(vec, line, boost::is_any_of(delimeter));
        std::vector<int> vec2;
        vec2.emplace_back(std::stoi(vec.data()[0]));
        vec2.emplace_back(std::stoi(vec.data()[1]));

        vectIN.push_back(vec2);
    }

    file.close();
    std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();


    std::vector<int> h_edge1, h_edge2;

    h_edge1.reserve(vectIN.size());
    h_edge2.reserve(vectIN.size());

    for(auto const& pts: vectIN)
    {
        h_edge1.emplace_back(pts[0]);
        h_edge2.emplace_back(pts[1]);
    }
    std::chrono::steady_clock::time_point t3 = std::chrono::steady_clock::now();

    thrust::device_vector<int> d_pStart(1);
    thrust::device_vector<int> d_pEnd(1);

    thrust::host_vector<int> h_rings;
    thrust::device_vector<int> d_rings;

    // Initialize edge vectors / pStart / pEnd /  while minimizing copying with CudaMalloc

    thrust::device_vector<int> d_edge1(vectIN.size());
    thrust::device_vector<int> d_edge2(vectIN.size());

    thrust::copy(thrust::make_zip_iterator(thrust::make_tuple(h_edge1.begin(), h_edge2.begin())),
                 thrust::make_zip_iterator(thrust::make_tuple(h_edge1.end(),   h_edge2.end())),
                 thrust::make_zip_iterator(thrust::make_tuple(d_edge1.begin(), d_edge2.begin())));

    // Arrange edges with edge1 as key and edge2 as value
    thrust::sort_by_key(d_edge1.begin(), d_edge1.end(), d_edge2.begin());

    d_rings.push_back(d_edge1.back());
    d_rings.push_back(d_edge2.back());

    d_edge1.pop_back();
    d_edge2.pop_back();

    d_pStart[0] = d_rings[0];
    d_pEnd[0] = d_rings[1];

    thrust::device_vector<int> element(1), p1(1), p2(1);

    while(not d_edge1.empty())
    {
        element.clear();

        int temp = d_pEnd[0];

        auto iter1 = thrust::equal_range(thrust::device, d_edge1.begin(), d_edge1.end(), temp);


        if(iter1.first != iter1.second)
        {
            element[0] = thrust::distance(d_edge1.begin(), iter1.first);
        }
        else
        {
            auto iter2 = thrust::find(thrust::device, d_edge2.begin(), d_edge2.end(), d_pEnd[0]);
            element[0] = thrust::distance(d_edge2.begin(), iter2);
        }

        // EDGE START INDEX (P1) AND END INDEX (P2)
        p1[0] = d_edge1[element[0]];
        p2[0] = d_edge2[element[0]];

        // ERASE THE EDGE FROM DEVICE LIST
        d_edge1.erase(d_edge1.begin()+element[0]);
        d_edge2.erase(d_edge2.begin()+element[0]);

        if(p1[0] == d_pEnd[0])
        {
            d_pEnd[0] = p2[0];

            if( d_pStart[0] == d_pEnd[0])
            {
                d_rings.push_back(-p2[0]);

                if(not d_edge1.empty())
                {
                    d_pStart[0] = d_edge1.back();
                    d_pEnd[0]   = d_edge2.back();

                    d_rings.push_back(d_pStart[0]);
                    d_rings.push_back(d_pEnd[0]);

                    d_edge1.pop_back();
                    d_edge2.pop_back();

                }
            }
            else
            {
                d_rings.push_back(p2[0]);
            }
        }
        else if(p2[0] == d_pEnd[0])
        {
            d_pEnd[0] = p1[0];

            if(d_pStart[0] == d_pEnd[0])
            {
                d_rings.push_back(-p1[0]);

                if(not d_edge1.empty())
                {
                    d_pStart[0] = d_edge1.back();
                    d_pEnd[0]   = d_edge2.back();

                    d_rings.push_back(d_pStart[0]);
                    d_rings.push_back(d_pEnd[0]);

                    d_edge1.pop_back();
                    d_edge2.pop_back();
                }
            }
            else
            {
                d_rings.push_back(p1[0]);
            }
        }
    }
    std::chrono::steady_clock::time_point t4 = std::chrono::steady_clock::now();

    // Copy rings to host and print them.
    h_rings = d_rings;

    for(auto const& pt:h_rings)
    {
        if(pt>=0)
            std::cout << pt << " ";
        else
            std::cout << -pt << std::endl;
    }
    std::cout << std::endl;

    long t1_t2 = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
    long t2_t3 = std::chrono::duration_cast<std::chrono::milliseconds>(t3 - t2).count();
    long t3_t4 = std::chrono::duration_cast<std::chrono::milliseconds>(t4 - t3).count();

    std::cout << "Load csv:      " << t1_t2 << std::endl;
    std::cout << "Create vector: " << t2_t3 << std::endl;
    std::cout << "Ring assembly: " << t3_t4 << std::endl;

    std::cout << "----------------- THAT'S ALL FOLKS!!! -----------------" << std::endl;

    return 0;
}

Other

I have implemented something similar to the above CUDA code but organized data into buckets so that searches would only have to be done on a limited amount of data. Unfortunately, I have not gotten it to work fully yet.

Recently I've been looking into graphing libraries to see if I could maybe do it that way, but I have not been successful in getting this way to work yet either. I know the CUDA toolkit has one as well as boost.

---

Final remarks

I would like to run this in at least less then 10 seconds, but ideally I would like it under a second for upto a million edges. I don't know if this is realistic but I am hoping that accelerating it with Cuda could achieve this or finding a different algorithm all together. I am reaching out to see if someone could please help me out to achieve this.

Answer 1

I bet a serial C++ implementation of Hierholzer's algorithm to find an Euler cycle would run in less than a second on 10^6 edges no problem, since the asymptotic running time is O(|E|). Given the tour, we still need to decompose it into simple cycles, which we can do with Python code like this (warning: untested).

def simple_cycles(tour_vertices):
    stack = []
    index = {}
    for v in tour_vertices:
        stack.append(v)
        i = index.get(v)
        if i is None:
            index[v] = len(stack) - 1
            continue
        yield stack[i:]
        for w in stack[i+1:]:
            del index[w]
        del stack[i+1:]

Here's complete C++ code for what I had in mind. Compiles but is otherwise completely untested. Absolutely no guarantees of any kind.

#include <list>
#include <stack>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

struct HalfEdge {
  int head;
  std::list<HalfEdge>::iterator opposite;
};

std::vector<std::vector<int>> SimpleCyclesFromEdges(
    const std::vector<std::pair<int, int>>& edges) {
  std::unordered_map<int, std::list<HalfEdge>> graph;
  for (const auto& edge : edges) {
    auto& first_neighbors = graph[edge.first];
    auto& second_neighbors = graph[edge.second];
    auto forward = first_neighbors.emplace(first_neighbors.begin());
    auto backward = second_neighbors.emplace(second_neighbors.begin());
    *forward = {edge.second, backward};
    *backward = {edge.first, forward};
  }
  std::vector<std::vector<int>> simple_cycles;
  for (auto& item : graph) {
    int v = item.first;
    std::unordered_set<int> on_stack;
    std::stack<int> stack;
    while (true) {
      if (on_stack.insert(v).second) {
        stack.push(v);
      } else {
        std::vector<int> cycle = {v};
        while (stack.top() != v) {
          cycle.push_back(stack.top());
          on_stack.erase(stack.top());
          stack.pop();
        }
        simple_cycles.push_back(std::move(cycle));
      }
      auto& neighbors = graph[v];
      if (neighbors.empty()) {
        break;
      }
      auto forward = neighbors.begin();
      v = forward->head;
      graph[v].erase(forward->opposite);
      neighbors.erase(forward);
    }
  }
  return simple_cycles;
}