How can we produce a graph layout with hierarchical edge bundling in Mathematica?

I'm looking to reproduce the types of graph layout that are possible in D3 (http://bl.ocks.org/mbostock/7607999) in Mathematica. Since Mathematica doesn't explicitly support force directed edge bundling, but does seem to support HierarchicalEdgeBundling (at least it's referred to once in the documentation) I wondered if anyone has had any success in doing something similar to the D3 visualisation linked above.

Update: I wrapped all this up into a small package for those who don't want to go through all the steps but would like to try this out anyway. Warning: There's not a lot of error checking and it may be very slow for more than a couple of hundred nodes.

The other answer I wrote shows how to use the builtin functionality. In this answer I am going to show how to implement such a graph layout from scratch. I hope that people will find this useful both from an educational point of view and to be able to customize the layout to their taste.

On the way we are going to get a little help from the IGraph/M package, the igraph interface for Mathematica. IGraph/M was in turn made possible by the LTemplate package.

How does the layout work?

This type of layout is useful because it makes the community structure in the graph evident. It is based on hierarchical community detection. A detailed description can be found in Y Jia, M Garland, JC Hart: Hierarchial edge bundles for general graphs. I will admit that I didn't actually read this paper. I only looked at the figures for inspiration. After all I am doing this for fun, not for a perfect result. However, if you want to get deeper into the topic, it is probably a good idea to read it.

The first reason why we need IGraph/M is that we are going to need to use a dendrogram output by some hierarchical community detection algorithm. IGraph/M has several such functions:

<<IGraphM


All these functions return an IGClusterData object. We can then query several "Properties" of such an object.

As an example, let us analyse the following network:

g = ExampleData[{"NetworkGraph", "DolphinSocialNetwork"}]


This is a social network between dolphins. Like most social networks, it has a relatively clear community structure.

cl = IGCommunitiesEdgeBetweenness[g]


Running the edge betweenness based (Girvan-Newman) community detection algorithm on it yields 5 communities. Note that we could also get this using the builtin FindGraphCommunities[..., Method -> "Centrality"], however, this function doesn't give us the full dendrogram.

Notice in the above screenshot of the IGClusterData object that it says: "Hierarchical: True". Not all algorithms included in igraph will produce a hierarchical structure, but this one does. We can visualize the dendrogram like this:

<< HierarchicalClustering
DendrogramPlot[cl["HierarchicalClusters"], LeafLabels -> (Rotate[#, Pi/2] &)]


Note: In IGraph/M 0.3.0 (not released as of this writing), it will be possible to plot the dendrogram using Dendrogram[cl["Tree"]]. Dendrogram is new in M10.4.

Other clustering algorithms in IGraph/M that can produce a dendrogram are Walktrap and Greedy, both being faster than EdgeBetweenness.

We can also obtain this dendrogram as a clustering tree, i.e. a Graph object.

tree = cl["Tree"]


So how does hierarchical edge bundling work? Instead of layout out the graph g, it constructs a version of the clustering tree and lays it out radially. The leaves correspond to the nodes of the original graph g. Then it uses this tree as a skeleton to route the edges of g between its nodes. The following figure from the paper gives an reasonably clear idea:

Let's go ahead then an implement this in Mathematica.

Implementing the layout

As a first step, we must simplify the clustering tree by identifying the subtrees corresponding to each cluster, and collapsing them into a single node, with all members of the community connecting directly to it.

Using a radial visualization of the tree, the graph on the left is what the algorithm gave us, and the one on the right is what we can use as a practical skeleton:

First of all, we must identify the root of the tree. Since this is a binary tree, we can just take the single node with degree 2. Leaves will have degree 1 and intermediate nodes will have degree 3.

root = First@Pick[VertexList[tree], VertexDegree[tree], 2]
(* 123 *)


Let us also get rid of the vertex labels and put back the nodes into the visualization, so we can see better what we are doing:

tree = SetProperty[RemoveProperty[tree, VertexLabels],
VertexShapeFunction -> "Circle"]


The clustering tree has integers as its nodes. We are going to need to forward and reverse mapping between the leaves of this tree and the graph g.

map = AssociationThread[Range@VertexCount[g], VertexList[g]];
revmap = Association[Reverse /@ Normal[map]];


Note: In IGraph/M 0.3.0 (not released as of this writing) it will be possible to simply use map = PropertyValue[tree, "LeafLabels"];

Now we group the leaves of the clustering tree based on which community they belong to. Note the use of Lookup and Map instead of ReplaceAll to prevent unpredictable replacements, especially in cases when some graph nodes are lists themselves.

communities = Lookup[revmap, #] & /@ cl["Communities"]
(* {{1, 3, 11, 29, 31, 43, 48},
{2, 6, 7, 8, 10, 14, 18, 20, 23, 26, 27, 28, 32, 33, 40, 42, 49, 55, 57, 58, 61},
{4, 9, 13, 15, 17, 21, 34, 35, 37, 38, 39, 41, 44, 45, 47, 50, 51, 53, 59, 60}, {5, 12, 16, 19, 22, 24, 25, 30, 36, 46, 52, 56},
{54, 62}} *)


To extract the subtree of a community, we make use of the fact that in an (undirected) tree there is precisely one path between any two nodes. Mapping enough paths to connect any two leaves will give us the subtree.

Clear[subtree]
subtree[tree_, {el_}] := {el}
subtree[tree_, community_] :=
Union @@ (First@FindPath[UndirectedGraph[tree], ##] &) @@@ Partition[community, 2, 1]

HighlightGraph[tree, subtree[tree, #] & /@ communities,


The root of a subtree is the node that appears first in the breadth-first ordering. Alternatively we could look for the only degree-2 node in the subtree, but let's use BFS here.  We can create is using [BreadthFirtScan](http://reference.wolfram.com/language/ref/BreadthFirtScan.html).

ord = First@Last@Reap[
];


For those not familiar with what breadth-first ordering is, the following animation will be educational. Nodes are visited in the order of their distance from the starting point—in this case from the tree root.

Animate[
HighlightGraph[
SetProperty[RemoveProperty[tree, VertexLabels], VertexShapeFunction -> "Circle"],
Take[ord, k]
],
{k, 1, VertexCount[tree], 1}
]


So let's get the roots of the subtrees:

posIndex = First /@ PositionIndex[ord];
subtreeRoots = First@MinimalBy[subtree[tree, #], posIndex] & /@ communities
(* {101, 119, 118, 112, 63} *)

HighlightGraph[tree, subtreeRoots, GraphHighlightStyle -> "Thick"]


New let's extract the tree which has these vertices as its leaves:

rootTree =
Subgraph[tree, VertexInComponent[tree, subtreeRoots],
GraphLayout -> "LayeredDigraphEmbedding"]


Notice that this is not going to make a useful skeleton because the different leaves have different distances from the root. Let us augment it with intermediate nodes to force all leaves to the same distance. A better method would not modify the tree—it would modify how the tree is embedded in space instead. But I am lazy so I will just do the augmentation here. I will use the symbolic wrapper a to generate names for the new nodes.

(* renaming the community roots will make things easier for cases where we have single-node communities *)
rootTree = VertexReplace[rootTree, Thread[subtreeRoots -> (a[#][0]&) /@ subtreeRoots]];
subtreeRoots = (a[#][0]&) /@ subtreeRoots;

dist = AssociationMap[GraphDistance[rootTree, root, #] &, subtreeRoots];
maxd = Max[dist];

augmentedRootTree =
Graph[
GraphUnion[
rootTree,
GraphUnion @@ Table[
PathGraph[
Prepend[Head[node] /@ Range[maxd - dist[node]], node],
DirectedEdges -> True
],
{node, subtreeRoots}
]
],
GraphLayout -> "LayeredDigraphEmbedding"
]


This is what we need. Now we just need to "hang" all nodes of g on the leaves of this tree. Let's extract the leaves in the same order as the communities:

rootTreeLeaves =

augmentedTree = SetProperty[
GraphUnion[
augmentedRootTree,
],
GraphLayout -> {"RadialDrawing", "RootVertex" -> root}
]


This is the tree we can use as a skeleton. But if we use the built-in radial tree visualization, then the leaves are not equispaced on a circle. Thus we employ the help of IGraph/M again and use its implementation of the Reingold-Tilford algorithm:

skeleton = IGLayoutReingoldTilfordCircular[augmentedTree,
"RootVertices" -> {root}]


We want to route the edges of g guided by the paths between them in the clustering tree, e.g.

spf = FindShortestPath[UndirectedGraph@augmentedTree, All, All];
HighlightGraph[
IGLayoutReingoldTilfordCircular[UndirectedGraph[augmentedTree],
"RootVertices" -> {root}],
Part[PathGraph /@ spf @@@ Map[revmap, EdgeList[g], {2}], 9],
GraphHighlightStyle -> "Thick"]


One simple way to do this is to use a BSplineCurve, with the intermediate points in the clustering tree as control points. But this will result in a very tight bundling of edges. To counter this, we can create new control points by interpolating between the original ones and a straight line going through the end vertices. The following function will do this:

Clear[smoothen]
smoothen[curve_, v_] :=
Module[{s = First[curve], t = Last[curve], line},
line = Table[s + (t - s) u, {u, 0, 1, 1/(Length[curve] - 1)}];
v line + (1 - v) curve
]


This function constructs the B splines:

Clear[plotGraph]
plotGraph[augmentedTree_, g_, v_: 0, sz_: 0.02] :=
Module[{pts, paths, spf},
spf = FindShortestPath[UndirectedGraph@augmentedTree, All, All];
pts = GraphEmbedding[augmentedTree];
paths = spf @@@ Map[revmap, EdgeList[g], {2}];
Graphics[
{
{
Opacity[0.5], ColorData["Legacy"]["RoyalBlue"],
Table[
With[
{curve =
PropertyValue[{augmentedTree, #}, VertexCoordinates] & /@
path},
BSplineCurve[smoothen[curve, v],
SplineDegree -> Length[curve] - 1]
],
{path, paths}
]
},
{
PointSize[sz], Black,
Point[PropertyValue[{augmentedTree, #}, VertexCoordinates]] & /@
Range@VertexCount[g]
}
}
]
]


Now let's take our skeleton and apply the visualization:

plotGraph[skeleton, g, 0.1]


It's interesting to use a Manipulate to control the tightness of the edge bundling:

Manipulate[plotGraph[skeleton, g, v], {v, 0, 1}]


Now with everything in place, we can try to use other skeletons as well. Using a standard top-to-bottom tree visualization the vertices are placed on a line:

skeleton = SetProperty[augmentedTree, GraphLayout -> "LayeredDigraphEmbedding"];
Show[
plotGraph[skeleton, g, 0.2, 0.01],
AspectRatio -> 1/3
]


skeleton = IGLayoutKamadaKawai[augmentedTree];
plotGraph[skeleton, g, 0.1]


With a bit more work we can also make this dynamic with clickable vertex labels, like in the example at http://bl.ocks.org/mbostock/7607999 The code is after the break.

I hope you enjoyed this little demo. I apologize for the messy code and for not packaging this up into a single function. I might do that on a better day when I have more time. Any feedback is welcome.

Code for Dynamic version

v = 0.15;
paths = spf @@@ Map[revmap, EdgeList[g], {2}];
With[{augmentedTree =
IGLayoutReingoldTilfordCircular[augmentedTree,
"RootVertices" -> {root}]},
DynamicModule[{flags =
ConstantArray[False, VertexCount[g]]]},
Graphics[
{
{
Opacity[0.5], ColorData["Legacy"]["RoyalBlue"],
Table[
With[{curve =
PropertyValue[{augmentedTree, #}, VertexCoordinates] & /@
path},
BSplineCurve[smoothen[curve, v],
SplineDegree -> Length[curve] - 1]
],
{path, paths}
],
Table[
With[{curve =
PropertyValue[{augmentedTree, #}, VertexCoordinates] & /@
path, f = First[path], l = Last[path]},
Dynamic@
Style[#,
If[flags[f] || flags[l],
Directive[Thickness[0.007], Red, Opacity[1]],
Opacity[0]]] &@
BSplineCurve[smoothen[curve, v],
SplineDegree -> Length[curve] - 1]
],
{path, paths}
]
},
{
PointSize[0.02],
With[{pt = PropertyValue[{augmentedTree, #}, VertexCoordinates],
offset = 1.05},
EventHandler[
{Dynamic[If[flags[#], Red, Black]],

Rotate[Text[VertexList[g][[#]],
offset pt, {-Sign@First[pt], 0}],
If[First[pt] < 0, Pi, 0] + ArcTan @@ pt, offset pt],
Point[pt]
},
{"MouseClicked" :> (flags[#] = ! flags[#])}
]
] & /@ Range@VertexCount[g]
}
},
ImageSize -> Large
]
]
]


Package

This is a small package containing the above functionality. To get started, simply apply HEBPlot to a connected graph which is not too large. Check Options[HEBPlot] for more control over the output. Warning: This is a rudimentary package with limited options and very little error checking.

Example:

HEBPlot[ExampleData[{"NetworkGraph", "JazzMusicians"}],
"CommunityDetectionFunction" -> IGCommunitiesGreedy,
EdgeStyle -> Directive[Opacity[0.2], RGBColor[0.254906, 0.411802, 0.882397]],
VertexSize -> 0.01]


BeginPackage["HierarchicalEdgeBundling", {"IGraphM"}];

HEBSkeleton::usage = "HEBSkeleton[clusterData] constructs a skeleton usable for hierarhical edge bundling. The argument must be an IGClusterData object.";
HEBEmbedSkeleton::usage = "HEBEmbedSkeleton[tree, layout] lays out a skeleton tree using \"Circular\" or \"Linear\" layouts.";
HEBLayout::usage = "HEBLayout[graph, tree] visualizes graph based on a precomputed skeleton tree.";
HEBPlot::usage = "HEBPlot[graph] visualizes graph using hierarchical edge bundling.";

Begin["Private"];

a; (* symbol for naming vertices augmenting the skeleton tree *)

subtree[tree_, {el_}] := {el}
subtree[tree_, community_?ListQ] :=
Union @@ (First@FindPath[UndirectedGraph[tree], ##]&) @@@ Partition[community, 2, 1]

HEBSkeleton::nohr = "The cluster data does not contain hierarchical information.";

HEBSkeleton[cl_IGClusterData] :=
Module[{tree, root, ord, revmap, communities, posIndex, subtreeRoots, rootTree, dist, maxd, augmentedRootTree, rootTreeLeaves},
If[Not@MemberQ[cl["Properties"], "Merges"],
Return@Failure["NotHierarchical", <|"MessageTemplate" -> HEBSkeleton::nohr|>]
];
tree = cl["Tree"];
root = First@Pick[VertexList[tree], VertexDegree[tree], 2];
communities = Lookup[revmap, #]& /@ cl["Communities"];
posIndex = First /@ PositionIndex[ord];
subtreeRoots = First@MinimalBy[subtree[tree,#], posIndex]& /@ communities;
rootTree = Subgraph[tree, VertexInComponent[tree, subtreeRoots]];
rootTree = VertexReplace[rootTree, Thread[subtreeRoots -> (a[#][0]&) /@ subtreeRoots]];
subtreeRoots = (a[#][0]&) /@ subtreeRoots;
dist = AssociationMap[GraphDistance[rootTree, root, #]&, subtreeRoots];
maxd = Max[dist];
augmentedRootTree = Graph[
GraphUnion[
rootTree,
GraphUnion @@ Table[
PathGraph[
DirectedEdges->True
],
{node, subtreeRoots}
]
]
];
GraphUnion[
augmentedRootTree,
]
]

HEBEmbedSkeleton[tree_?TreeGraphQ, layout_String : "Circular"] :=
Switch[layout,
"Circular", IGLayoutReingoldTilfordCircular[tree, "RootVertices" -> Pick[VertexList[tree], VertexInDegree[tree], 0]],
"Linear", SetProperty[tree, GraphLayout -> "LayeredDigraphEmbedding"]
]

smoothen[curve_,v_]:=
Module[{s = First[curve], t = Last[curve], line},
line = Table[s+(t-s)u,{u,0,1,1/(Length[curve]-1)}];
v line+(1-v) curve
]

Options[HEBLayout] = {
"BundleTightness" -> 0.1,
VertexSize -> 0.02,
VertexStyle -> Black,
EdgeStyle -> Directive[Opacity[0.5], ColorData["Legacy"]["RoyalBlue"]]
};

HEBLayout[g_?GraphQ, tree_?TreeGraphQ, opt : OptionsPattern[]] :=
Module[{paths, spf, v = OptionValue["BundleTightness"], revmap},
spf = FindShortestPath[UndirectedGraph[tree], All, All];
paths = spf@@@Map[revmap, EdgeList[g], {2}];
Graphics[
{
{
OptionValue[EdgeStyle],
Table[
With[
{curve=PropertyValue[{tree, #}, VertexCoordinates]&/@path},
BSplineCurve[smoothen[curve,v], SplineDegree -> Length[curve]-1]
],
{path,paths}
]
},
{
PointSize@OptionValue[VertexSize], OptionValue[VertexStyle],
Point[PropertyValue[{tree, #},VertexCoordinates]]&/@Range@VertexCount[g]
}
}
]
]

Options[HEBPlot] = Options[HEBLayout] ~Join~ {"Layout" -> "Circular", "CommunityDetectionFunction" -> IGCommunitiesEdgeBetweenness};

HEBPlot[g_?ConnectedGraphQ, opt : OptionsPattern[]] :=
Module[{cl, skeleton},
cl = OptionValue["CommunityDetectionFunction"][g];
skeleton = HEBSkeleton[cl];
skeleton = HEBEmbedSkeleton[skeleton, OptionValue["Layout"]];
HEBLayout[g, skeleton, FilterRules[{opt}, Keys@Options[HEBLayout]]]
]

End[];

EndPackage[];

• Fantastic work! – DPF Oct 5 '16 at 16:09
• Re-posted on Wolfram Community here. Any updates I may make in the future will first appear here on SE due to convenience. – Szabolcs Oct 6 '16 at 19:59

Using the "HierarchicalEdgeBundling" layout is a bit tricky. It is not very well documented, so it takes a bit of a guesswork to figure out the syntax. It is described here:

You can do it like this:

g = ExampleData[{"NetworkGraph", "DolphinSocialNetwork"}]

SetProperty[g, GraphLayout -> {"EdgeLayout" -> "HierarchicalEdgeBundling"}]


The key is to set the "EdgeLayout" suboption of GraphLayout and not to set a vertex layout. I was confused by this at first as edge and vertex layouts are separate for the other options, and we must set a vertex layout too. Doing so here will prevent this from working.

According to the comments in the linked thread, there is one working suboption, "LCARemove", which is short for "lowest common ancestor removal". More info can be found in this paper. Valid values are All and None.

SetProperty[g,
GraphLayout -> {"EdgeLayout" -> {"HierarchicalEdgeBundling",
"LCARemove" -> #}}] & /@ {All, None}


This type of visualization is rather spectacular so it is a pity that it is not better documented.

• I presume this could be used for this question as well? – J. M.'s discontentment Oct 5 '16 at 15:42
• @J.M. Maybe ... not sure how well it will do with that many nodes and edges. There's also the thing that the vertex ordering around the circle is determined automatically, based on the communities in the graph. I don't think it would work well, but it's worth a try. – Szabolcs Oct 7 '16 at 8:57