When we got to the end of my earlier post on the topic, we’d built a nice little graph data structure together with a core set of query operations to calculate vertex degrees and vertices’ neighbours. But all of this assumed we’d already got the graph data in place — when of course we haven’t! So this time round, we’ll examine the functions needed to add vertices and edges to the network.

Interruption!If you’re interested in following developer blogs on the D programming language, do check out John Colvin’s`foreach(hour; life) {/* do something here */}`

. John’s a lovely guy who contributes much helpful advice and insight on the D forums.

Let’s start by reminding ourselves of the basic data structure we’re dealing with:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
final class Graph(bool dir) { private: size_t[] _head; size_t[] _tail; size_t[] _indexHead; size_t[] _indexTail; size_t[] _sumHead = [0]; size_t[] _sumTail = [0]; public: enum bool directed = dir; // ... Yes, there was more than this, but this is what we need to recall for now ... } |

`_head`

and `_tail`

contain respectively the head and tail vertices of each edge; `_indexHead`

and `_indexTail`

contain indices of the edges sorted respectively according to head and tail IDs; and `_sumHead`

and `_sumTail`

contain cumulative sums of the number of edges with head/tail less than the given head/tail value. This is a transcription into D of the `igraph_t`

datatype from igraph; Gábor Csárdi, one of the authors of igraph, refers to it as an *indexed edge list*.

Let’s start with adding vertices. We could get a bit complicated and make the `vertexCount`

property writeable, but for now let’s keep it simple and just match igraph with a function to increase the number of vertices in the graph:

1 2 3 4 5 6 7 8 |
void addVertices(immutable size_t n) { immutable size_t l = _sumHead.length; _sumHead.length += n; _sumTail.length += n; _sumHead[l .. $] = _sumHead[l - 1]; _sumTail[l .. $] = _sumTail[l - 1]; } |

All we have to do is make sure that the cumulative sums of edges get extended, and since the new vertices don’t have any edges, this is a simple matter of copying the final value from the original sum array. We’re doing this using D’s slicing syntax: `arr[i .. j]`

simply means the sub-array going from the `i`

’th index up to the `j - 1`

’st index; the dollar sign `$`

is shorthand for the array’s length. If you’re interested, you can read in much more detail about slices in Ali Çehreli’s excellent online book.

Anyway, adding vertices is boring — it’s adding edges where the exciting stuff starts to happen. In principle there are two ways we can do this: we can add individual edges, one at a time, or we can add *many* edges, all in one go; igraph defines two different functions, `igraph_add_edge`

and `igraph_add_edges`

, to do this.

We’ll go with adding one at a time to start with. The task is simple — we need to add the head and tail values, place the new edge ID in its correct location in the head/tail-sorted indices, and increment the sum values:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
public: void addEdge(size_t head, size_t tail) { static if (!directed) { if (tail < head) { swap(head, tail); } } _head ~= head; _tail ~= tail; indexEdges(); ++_sumHead[head + 1 .. $]; ++_sumTail[tail + 1 .. $]; } private: void indexEdges() { _indexHead ~= iota(_indexHead.length, _head.length).array; _indexTail ~= iota(_indexTail.length, _tail.length).array; _indexHead.sort!((a, b) => _head[a] < _head[b]); _indexTail.sort!((a, b) => _tail[a] < _tail[b]); } |

I’m actually cheating a bit here, because when I originally wrote the `indexEdges()`

function, I used `schwartzSort!(a => _head[a], "a < b")`

and not the `sort`

formulation given here. `schwartzSort`

caches the results of the transformation `a => f(a)`

, making it cheaper where that transformation is costly; but as `a => _head[a]`

is so cheap, we’re actually better off using regular `sort`

, and the `(a, b)`

transformation here provides the rule for determining whether or not `a`

should come before `b`

in the sort order.

Note that we insist that in an undirected graph, the head of an edge has a vertex ID less than or equal to that of the tail. This may seem like a fussy convention but it’s actually useful — we’ll see why later.

Anyway, this all seems nice and in order, so why not test it? Something like this:

1 2 3 4 5 6 7 8 9 10 11 |
// Let's suppose we have a list of edges stored as an array ... size_t[] edges = [head1, tail1, head2, tail2 /* ... and so on */]; // We create an undirected graph (but it could also be directed if we wanted:-) auto g = new Graph!false; g.addVertices(/* however many we need */); foreach(i; 0 .. edges.length / 2) { g.addEdge(edges[2 * i], edges[2 * i + 1]); } |

Benchmarking this with a 50-vertex graph with 200 edges, the D code comes out more than 10 times faster than similar code written using igraph. A win, you think? You speak too soon — if we try 10,000 vertices and 20,000 edges, the D code is over 7 times slower! You might think that something must be seriously wrong here, and you’d be right — in fact, we’ve run into an unfortunate bug in D’s sort algorithm. In the corner case where an array is almost entirely sorted except for the last few entries, the sort displays worst-case behaviour and scales with *n*^{2} rather than the expected *n* log *n*. Fortunately for us there’s an alternative sorting algorithm available, which uses Timsort rather than the standard Quicksort:

1 2 |
_indexHead.sort!((a, b) => _head[a] < _head[b], SwapStrategy.semistable); _indexTail.sort!((a, b) => _tail[a] < _tail[b], SwapStrategy.semistable); |

This actually works out slightly slower on the smaller graph, but gives a major speedup for the larger one — suddely what was taking over a minute to run on my system is running in less than a second! Most likely the scaling bug wasn’t being triggered at all for the smaller network, but kicks in as the length of the arrays gets above a certain size, and *O*(*n*^{2}) as *n* → several thousand is going to leave a mark …

… but wait a minute. Aren’t we missing something here? We’re talking about adding one (or very few) new values to an already sorted list — we shouldn’t be using QuickSort or TimSort at all, we should be using an insertion sort. And it turns out that we can write this very simply indeed — even simpler than I’d anticipated, thanks to a technique pointed out to me by Dmitry Olshansky:

1 2 3 4 5 6 7 8 9 10 11 |
void indexEdgesInsertion() { immutable size_t l = _indexHead.length; foreach(e; iota(l, _head.length)) { size_t i = _indexHead.map!(a => _head[a]).assumeSorted.lowerBound(_head[e]).length; insertInPlace(_indexHead, i, e); i = _indexTail.map!(a => _tail[a]).assumeSorted.lowerBound(_tail[e]).length; insertInPlace(_indexTail, i, e); } } |

`lowerBound`

uses a binary search to extract from the already-sorted index list the largest left-hand part of the array such that these edges’ head/tail values are all less than the head/tail value of our new edge. `insertInPlace`

then drops the new index value into its correct location. This is going to be expensive if instead we’re adding multiple edges in one go, so let’s make it nicer by doing all the potential memory allocation in one go:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
void indexEdgesInsertion() { immutable size_t l = _indexHead.length; _indexHead.length = _head.length; // all the allocation _indexTail.length = _tail.length; // at once ... foreach(e; iota(l, _head.length)) { size_t i, j; i = _indexHead[0 .. e].map!(a => _head[a]).assumeSorted.lowerBound(_head[e]).length; for(j = e; j > i; --j) _indexHead[j] = _indexHead[j - 1]; _indexHead[i] = e; i = _indexTail[0 .. e].map!(a => _tail[a]).assumeSorted.lowerBound(_tail[e]).length; for(j = e; j > i; --j) _indexTail[j] = _indexTail[j - 1]; _indexTail[i] = e; } } |

… and suddenly our adding of new nodes is super-fast all round — whereas it takes igraph about 13 seconds on my machine to build a 10,000-vertex, 20,000-edge graph adding edges one at a time, the D code here delivers it in less than 0.1 seconds. Woop!

* * *

Let’s introduce a note of cold sobriety, before we go rushing off and opening the champagne. To compare with igraph like this may give a frisson of pleasure, but we’re not comparing like with like. igraph’s `igraph_add_edge`

function simply wraps its function for adding multiple edges:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
int igraph_add_edge(igraph_t *graph, igraph_integer_t from, igraph_integer_t to) { igraph_vector_t edges; int ret; IGRAPH_VECTOR_INIT_FINALLY(&edges, 2); VECTOR(edges)[0]=from; VECTOR(edges)[1]=to; IGRAPH_CHECK(ret=igraph_add_edges(graph, &edges, 0)); igraph_vector_destroy(&edges); IGRAPH_FINALLY_CLEAN(1); return ret; } |

This is bound to take a performance hit, because the function for adding multiple edges does a complete rearrangement of the sorted index and sum arrays — whereas our `addEdge`

function just tweaks the bits it needs to. The same could be implemented in igraph, and it would almost certainly beat the present D code in speed. The take-home point here isn’t that D is able to be faster — it’s that because we can write so quickly and simply and effectively in D, we had the time and mental space to implement a tailored solution like this from the get-go. (It also helps that D’s community is very friendly and helpful with suggestions!)

The second sobriety check comes when we try using igraph how it’s meant to be used — passing the graph all the edges in one big go as a vector. Suddenly it’s generating that big 10,000-vertex, 20,000-edge network in about 0.001 seconds! We clearly have some catching up to do.

Let’s assume, as igraph does (and as we already did in an example above) that our list of many edges is going to come in the form of an array `edges = [head1, tail1, head2, tail2, head3, tail3, ...]`

(maybe we’ll change this in future, but for now it serves). So now let’s write a variant of addEdge that will accept this:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 |
void addEdge(T : size_t)(T[] edgeList) { immutable size_t l = _head.length; _head.length += edgeList.length / 2; _tail.length += edgeList.length / 2; foreach(i; 0 .. edgeList.length / 2) { size_t head = edgeList[2 * i]; size_t tail = edgeList[2 * i + 1]; static if (!directed) { if (tail < head) { swap(head, tail); } } _head[l + i] = head; _tail[l + i] = tail; } indexEdgesInsertion(); sumEdges(_sumHead, _head, _indexHead); sumEdges(_sumTail, _tail, _indexTail); } void sumEdges(ref size_t[] sum, ref size_t[] vertex, ref size_t[] index) { size_t v = vertex[index[0]]; sum[0 .. v + 1] = 0; for(size_t i = 1; i < index.length; ++i) { size_t n = vertex[index[i]] - vertex[index[sum[v]]]; sum[v + 1 .. v + n + 1] = i; v += n; } sum[v + 1 .. $] = vertex.length; } |

So, we add the head/tail values to their respective arrays, we sort the indices, and we recalculate from scratch the `_sumHead`

and `_sumTail`

arrays. If a large number of edges has been added, the latter should be faster than just incrementing the sums one by one as we did when adding a single edge. For now we’ll stick with the insertion sort — it’s interesting to examine how it performs.

We’re obviously doing something right, because on my machine the time to create the smaller, 50-vertex graph is cut to less than a quarter of what it was before, placing it on a par with igraph. For the 10,000-vertex graph it’s almost 5 times faster than it was adding edges one at a time — but that’s still an order of magnitude slower than igraph. Presumably this is down to the sorting technique — if we’re adding many edges all at once, then insertion sort is probably no longer optimal. Let’s try out a regular sort:

1 2 3 4 5 6 7 8 |
void indexEdgesSort() { _indexHead ~= iota(_indexHead.length, _head.length).array; _indexTail ~= iota(_indexTail.length, _tail.length).array; assert(_indexHead.length == _indexTail.length); _indexHead.sort!((a, b) => _head[a] < _head[b]); _indexTail.sort!((a, b) => _tail[a] < _tail[b]); } |

This cuts the creation time for the 50-vertex graph very slightly, and for the 10,000-vertex graph very significantly — the latter is about 8 times faster than it was before! This is probably still very slightly slower than igraph, but only slightly — it’s about 0.002 seconds to generate the larger graph with D compared to about 0.0015 with igraph. (In case you’re wondering, this is an average over 1000 realizations, which takes about 1.5s for igraph, just over 2 for the D code.)

Clearly we’re standing on the shoulders of giants — we’re able to develop these solutions quickly and easily because we have the example of igraph already there in front of us to reverse engineer. The `sumEdges()`

function, for example, is a fairly straight D-ification of igraph’s `igraph_i_create_start()`

function. But as with the previous post, the key thing to note is how cleanly we can implement these ideas; there’s no unnecessary cruft or boilerplate needed to ensure the code operates safely, and many of the solutions are taken straight out of the D standard library, Phobos. Others, such as sorting by head or tail, can be implemented very simply because of the elegant syntax D offers us for delegates and lambdas.

* * *

Actually we’re not quite finished here. We sorted `_indexHead`

and `_indexTail`

according to `_head`

and `_tail`

values respectively, but in fact we want a double sort criterion as there is in igraph: `_indexHead`

should be sorted primarily with respect to `_head`

, but secondarily with respect to `_tail`

; and vice versa for `_indexTail`

. Thanks to advice from suggestions by John Colvin and bearophile on the D forums, we know there is a ready solution for this in the form of the multiSort function:

1 2 3 4 5 6 7 |
void indexEdgesSort() { _indexHead ~= iota(_indexHead.length, _head.length).array; _indexTail ~= iota(_indexTail.length, _tail.length).array; _indexHead.multiSort!((a, b) => _head[a] < _head[b], (a, b) => _tail[a] < _tail[b]); _indexTail.multiSort!((a, b) => _tail[a] < _tail[b], (a, b) => _head[a] < _head[b]); } |

… while a small tweak to the insertion sort in D ensures correct ordering there as well:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
void indexEdgesInsertion() { immutable size_t l = _indexHead.length; _indexHead.length = _head.length; _indexTail.length = _tail.length; foreach(e; l .. _head.length) { size_t i, j, lower, upper; upper = _indexHead[0 .. e].map!(a => _head[a]).assumeSorted.lowerBound(_head[e] + 1).length; lower = _indexHead[0 .. upper].map!(a => _head[a]).assumeSorted.lowerBound(_head[e]).length; i = lower + _indexHead[lower .. upper].map!(a => _tail[a]).assumeSorted.lowerBound(_tail[e]).length; for(j = e; j > i; --j) _indexHead[j] = _indexHead[j - 1]; _indexHead[i] = e; upper = _indexTail[0 .. e].map!(a => _tail[a]).assumeSorted.lowerBound(_tail[e] + 1).length; lower = _indexTail[0 .. upper].map!(a => _tail[a]).assumeSorted.lowerBound(_tail[e]).length; i = lower + _indexTail[lower .. upper].map!(a => _head[a]).assumeSorted.lowerBound(_head[e]).length; for(j = e; j > i; --j) _indexTail[j] = _indexTail[j - 1]; _indexTail[i] = e; } } |

We do take a performance hit here — we’re no longer level pegging with igraph, although we might get back that speed by using instead the radix sort that is used by igraph. But the hit is worth it: among other things, this extra level of sorting means we are guaranteed to get a sorted list out of the `neighbours`

function, but its main application is in making it easier to search for edges and their place in the edge list. We’ll examine this next time when we look in more detail at the range of queries we can make of a graph object, and how they shape up in performance.

**Edit 19-07-2013:** I’ve tweaked the article to make clear that the double sort criterion described here is also used by igraph and — unsurprisingly — was the inspiration for doing it here as well.