= Graph Visualization =

I have a soft spot for graph layout algorithms because as an undergraduate, I
signed up for a project that led me to implement what
http://www.math.ntua.gr/~symvonis/[my professor] called "the spring algorithm"
(Eades 1984). These days, it appears to be called
https://en.wikipedia.org/wiki/Force-directed_graph_drawing[force-directed graph
drawing].

By _graph_, I mean what Sylvester meant when he coined the term in 1878: a set
of _nodes_ (or _vertices_ or _points_) connected by _edges_ (or _lines_ or
_links_). I wish he had called it a "web" instead, because "graph" already has
a meaning in mathematics, and anybody who's surfed the worldwide *web* already
understands nodes connected with links.

Hidden in the physics building, at the end of the hallway on the second floor,
there was a computer lab that housed a bevy of Silicon Graphics workstations,
each capable of real-time ray-tracing. Behind a glass window in a darkened room
loomed a Cray supercomputer, which looked like an altar. Perhaps it really was
one, as only the high priest who oversaw the lab could perform rites on it, and
ask it for blessings on our behalf.

Like a secret agent, I was issued a nondescript security card that I'd swipe at
a nondescript door to gain entry, while the rest of the computer science
students were condemned to crummy terminals in the stuffy, crowded, squalid
basement of the Madsen building, sharing a single Sun server. With the aid of a
futuristic API named OpenGL, I was able to render graphs as beautifully shaded
spheres and cylinders in three dimensions. What a time to be alive!

Today, those miracles I witnessed in that lab are commonplace. Mobile phones
are more powerful than the vaunted Cray, whose manufacturer is long gone, as
are Sun Microsystems, Inc. and Silicon Graphics, Inc. (I later spent years in
the former SGI headquarters during my stint at Google; one of my many
co-workers was one of SGI's earliest employees.)

My brain also seems to have gone through similar upgrades somewhere along the
way. Back then I painstakingly edited page after page of code over a semester
to lay out a graph; today, I can get it done by busting out a few lines.

== Spring time ==

Working in a room in the physics building all those years ago was appropriate,
because I was in fact writing a simple physics engine. Given a connected graph,
we place the vertices randomly. We imagine all nodes to be positive charges,
and all edges to be springs of the same length, then apply the laws of physics.

We copy the pseudo-random number generator from link:../play/netwalk.html[our NetWalk page]:

data PCG = PCG Word64 Word64
pcg a b = PCG (a + b') b' where b' = 2*b + 1

next :: PCG -> (Word, PCG)
next (PCG x inc) = (r, PCG x' inc) where
  x' = 6364136223846793005*x + inc
  r = (fromIntegral (x `xor` (x `shiftR` 18) `shiftR` 27) :: Word) `rotateR` (fromIntegral $ x `shiftR` 59)

fromPCG p = map (fromIntegral . fst) $ tail $ iterate (next . snd) (undefined, p)

split :: PCG -> (PCG, PCG)
split p = (PCG (lohi a b) (lohi c d), PCG (lohi e f) (lohi g h)) where
  [a,b,c,d,e,f,g,h] = take 8 $ fromPCG p
  lohi = Word64

chunksOf i ls = map (take i) (go ls) where
  go [] = []
  go l  = l : go (drop i l)

We generate coordinates using a seed that happens to work nicely with our
examples:

points42 = chunksOf 2 $ (\n -> doubleFromWord n / 2^25) <$> fromPCG (pcg 42 42)

In our examples we use integers to refer to the nodes, though our code is more
general.

We store the nodes and their coordinates in an association list. The edges
are represented by a function that returns whether two nodes are connected.
For example, the following is \(K_5\), the complete 5-node graph, with an
initial random layout:

graphK5 = (\_ _ -> True, zip [1..5] points42)

We display the graphs with SVG elements which we insert into the HTML of
this page:

allPairs as = [(x, y)
  | x:t <- takeWhile ((> 1) . length) $ iterate tail as, y <- t]

svgline [x1, y1] [x2, y2] = ("<line x1='"<>)
  . shows x1 . ("' y1='"<>) . shows y1 . ("' x2='"<>)
  . shows x2 . ("' y2='"<>) . shows y2 . ("' stroke='red' />"<>)

svgcircle [x, y] = ("<circle cx='"<>)
  . shows x . ("' cy='"<>)
  . shows y . ("' r='2' />"<>)

svgWrap (isEdge, vs) = ("<svg viewBox='"<>)
  . (foldr (.) id $ intersperse (' ':) $ shows <$> [x0 - mx, y0 - my, dx + 2*mx, dy + 2*my])
  . ("' width=300 height=200>"<>)
  . foldr (.) id [ svgline p q | ((v, p), (w, q)) <- allPairs vs, isEdge v w]
  . foldr (.) id (svgcircle . snd <$> vs)
  $ "</svg>"
  where
  xs = (!!0) . snd <$> vs
  ys = (!!1) . snd <$> vs
  x0 = foldr1 min xs
  x1 = foldr1 max xs
  dx = x1 - x0
  y0 = foldr1 min ys
  y1 = foldr1 max ys
  dy = y1 - y0
  mx = max 2 $ dx * 0.05
  my = max 2 $ dy * 0.05

draw g = jsEval $ "repl.outdiv.insertAdjacentHTML('beforeend',`" ++ svgWrap g ++ "`);"

draw graphK5

== Hooked on Coulomb ==

Let \(d\) be the distance between two nodes. If they're connected by an edge,
the force exerted on each node away from the other node is given by
Coulomb's Law and Hooke's Law:

\[
F = \frac{k_1}{d^2} + k_2(L - d)
\]

where \(L\) is the length of each spring (at equilibrium), and \(k_1, k_2\) are
constants. If the nodes have no edge between them, then we drop the second
term, as then we only wish to simulate the repulsion of the two charges.

Each frame, for each node, we sum the forces on it due to every other node,
multiply the total by a tunable tiny constant, and adjust its position by the
result.

This makes sense, because for constant acceleration, the distance traveled over
a given time is proportional to the acceleration; of course, the forces
continuously vary as the nodes move around, which is why we try to find a
duration short enough for a passable approximation, though not so short that
we lose our patience waiting for the simulation.

In link:re.html[our regular expressions demo], zero-length springs (\(L = 0\))
happen to work with my example, but it appears to perform poorly for the graphs
below. I settled on the constants in the code after some experimentation.

sim (isEdge, vs) = (isEdge, nudge <$> vs) where
  nudge (v, p0) = (v, foldr ($) p0 $ [zipWith (+) $ delta (isEdge v w) $ zipWith (-) p1 p0 | (w, p1) <- vs, v /= w])

delta sprung dp = (f*) <$> dp where
  d = sqrt $ sum $ map (^2) dp
  f = -0.2 * (9000/(d*d) + bool 0 (20 - d) sprung)/d

Let's look at the first few frames when drawing \(K_5\).

mapM draw $ take 8 $ iterate sim graphK5

The result is already pleasing, though complete graphs only exercise one code
path.

== Good Will Hunting II ==

The movie _Good Will Hunting_ features a couple of combinatorics problem sets.
link:goodwill.html[We explore the first set elsewhere]. Here, we address the
second set, as it's a good excuse to show off our code.

The first problem asks for the number of \(n\)-node labeled trees.

A _tree_ is a connected graph with no cycles, that is, given any two nodes,
there is a unique path between them. We label the nodes 1 to \(n\) and consider
two trees to be equivalent if two nodes are adjacent on one tree exactly when
the identically labeled nodes are adjacent on the other.

Counting the first 5 cases is easy, and provides enough data for a good
educated guess at the general formula. I urge you to try before reading on.

If \(n = 1, 2\) there is only one possible unlabeled tree, and only one way
to label it uniquely. If \(n = 3\), the only possible unlabeled tree is
3 nodes connected in series, but there are 3 different labelings, each
determined by the choice of the label for the node in the middle.

fromEdgeStr s = (edgeP, zip vs $ points42) where
  edgeP v w = [min v w, max v w] `elem`
    chunksOf 2 (fromIntegral . readInteger <$> words s)
  ns = fromIntegral . readInteger <$> words s :: [Int]
  vs = nub ns

nub = foldr go [] where
  go x xs
    | elem x xs = xs
    | otherwise = x:xs

mapM draw
  [ (undefined, [(1, [0, 0])])
  , fromEdgeStr "1 2"
  , (!!50) $ iterate sim $ fromEdgeStr "1 2 2 3"
  ]

If \(n = 4\) then there are two shapes: either all 4 nodes are connected in
series, or there is one node of degree 3 that connects to the others.
In the latter case, there are 4 labelings, where again each is determined by
the choice of the label of the central node. In the former case, there are \(4!
= 24\) permutations we can arrange 4 labels in a row, but this counts each
labeling twice: rotating the line of nodes by 180 degrees yields the reverse of
their permutation. Thus there are a total of \(4 + 24/2 = 16\) trees.

mapM (draw . (!!50) . iterate sim . fromEdgeStr)
  [ "1 2 2 3 3 4"
  , "1 2 1 3 1 4"
  ]

If \(n = 5\), then one of the following occurs:

  * we have 5 nodes in series, yielding \(5! / 2 = 60\) unique labelings.

  * we have one node of degree 4 connected to the others, yielding 5 unique
  labelings.

  * we have one node of degree 3, one of degree 2, and the other nodes have degree 1, yielding \(5\times 4\times 3 = 60\) unique labelings which are determined by the choices of labels of the degree 2 node and its neighbours.

mapM (draw . (!!50) . iterate sim . fromEdgeStr)
  [ "1 2 2 3 3 4 4 5"
  , "1 2 1 3 1 4 1 5"
  , "1 2 1 3 1 4 2 5"
  ]

This gives a total of 125 labeled trees. So we now have the sequence
\(1, 1, 3, 16 = 4^2, 125 = 5^3\), which leads to the guess:

\[
n^{n-2}
\]

This guess turns out to be correct, and is known as Cayley's formula, who
proved it in the 19th century
(see
https://rcin.org.pl/impan/dlibra/publication/176630/edition/143722/content[Cayley,
_A theorem on trees_]).
Will Hunting writes it on the board with no proof. However, since that boy is
wicked smart, he would likely be able to supply one.

== Good Arthur Cayley ==

Many famous mathematicians have tackled Cayley's formula. Aigner and Ziegler,
_Proofs from THE BOOK_, walk through four distinct proofs of Cayley's formula,
each of which is ingenious and elegant. Because of my background in
cryptography, my favourite is the following due to Joyal.

On a labeled tree with \(n\) nodes, choose any node to be the "start", and any
node to be the "end"; these may be the same node. Since it's a tree, there is a
unique path from the start node to the end node, which we'll call the "spine".

Let \(a b ... z\) be the labels of the nodes of the spine sorted in ascending
order. Then the path from the start node to the end node describes a
permutation of \(a b ... z\) which we denote \(\pi\).

Define a function \(f\) from the set of naturals \([1..n]\) to itself as
follows. If \(k \in [1..n] \) is on the spine, then \(f(k) = \pi(k)\).
Otherwise, there is a unique path from \(k\) to the start node (in fact, we can
use any node on the spine). Start walking along this path, and define \(f(k)\)
to be the first node you reach. Different labelings and choices of start and
end nodes lead to different functions.

There are \(n\) choices for the start node, and same for the end node. As there
are \(n^n\) different functions from \([1..n]\) to itself, we have shown there
are \(n^{n-2}\) different labeled trees...provided we fill in a gap.

To prove we have constructed a bijection, we need the flip side: we must show
that any function \(f\) from \([1..n]\) to itself corresponds to a labeled
tree. But consider repeatedly applying \(f\) on some \(k \in [1..n]\). Either
\(k\) is in a cycle, or we'll eventually reach a value that is in a cycle.

Thus we can reverse the above process. Together, the cycles describe a
permutation \(\pi\) of some of the elements of \([1..n]\), which we use to
build a spine. As for the other elements, applying \(f\) describes how to
connect them via subtrees rooted to the spine.

In cryptography,
https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm[Pollard's rho
algorithm] repeatedly applies a function from a finite set to itself,
eventually causing the values to cycle, which is exploited to factor an integer
or find a discrete logarithm.

== Robot Will Hunting ==

The last task is to draw all homeomorphically irreducible trees with 10 nodes.
This is easy, so we up the ante with a layer of indirection: our goal is to
write a program that draws all homeomorphically irreducible trees with 10
nodes.

Let's briefly review some definitions. A node is _external_ if it has degree 1,
and is _internal_ otherwise. A graph is _homeomorphically irreducible_ or just
_irreducible_ if no node has degree 2. This is the same as saying all internal
nodes have degree 3 or more.

We define two graphs to be equivalent (_isomorphic_) if both graphs have \(n\)
nodes, and we can label the nodes of each graph with \([1..n]\) such that two
nodes are connected in one graph if and only if in the other graph, the nodes
with the same labels are connected.

Now that we understand the problem, we simply exhaust all possibilities and
feed any trees we find into our code above to draw them.

We first notice that the internal nodes must form a tree. Suppose there are 4
internal nodes. As above, they are either they are connected in series, or one
node is connected to each of the other three. We find adding 6 external nodes
are barely enough to obtain irreducible trees, thus we have drawn 2 of the
desired trees, and have also discovered there can be no more than 4 internal
nodes (we're waving our hands here, but this argument can easily be made more
rigorous).

t4_1 = fromEdgeStr "1 2 2 3 3 4 1 5 1 6 2 7 3 8 4 9 4 10"
t4_2 = fromEdgeStr "1 2 1 3 1 4 2 5 2 6 3 7 3 8 4 9 4 10"
mapM (draw . (!!80) . iterate sim ) [t4_1, t4_2]

We continue bashing through cases. If there is 1 internal node, then there's
only one way to attach 9 external nodes to it. If there are 2 internal nodes,
they must be adjacent, and we find 3 different ways to attach external nodes to
reach 10 nodes total.

t1 = (\v w -> v == 1 || w == 1, zip [1..10] $ points42)
t2_1 = fromEdgeStr "1 2 1 3 1 4 1 5 1 6 1 7 1 8 2 9 2 10"
t2_2 = fromEdgeStr "1 2 1 3 1 4 1 5 1 6 1 7 2 8 2 9 2 10"
t2_3 = fromEdgeStr "1 2 1 3 1 4 1 5 1 6 2 7 2 8 2 9 2 10"
mapM (draw . (!!50) . iterate sim ) [t1, t2_1, t2_2, t2_3]

If there are 3 internal nodes, they must be connected in series, and we find 4
different ways to attach the external nodes: the degrees of the 3 nodes are
5-3-3, 4-4-3, 4-3-4, and 3-5-3.

t3_1 = fromEdgeStr "1 2 2 3 1 4 1 5 2 6 3 7 3 8 1 9 1 10"
t3_2 = fromEdgeStr "1 2 2 3 1 4 1 5 2 6 3 7 3 8 1 9 2 10"
t3_3 = fromEdgeStr "1 2 2 3 1 4 1 5 2 6 3 7 3 8 1 9 3 10"
t3_4 = fromEdgeStr "1 2 2 3 1 4 1 5 2 6 3 7 3 8 2 9 2 10"
mapM (draw . (!!50) . iterate sim ) [t3_1, t3_2, t3_3, t3_4]

In the story, Will Hunting draws only 8 trees before he is interrupted and
chased off. We see the board is missing the case with 4 internal nodes in
series and the case with 3 internal nodes with degrees 5-3-3. Speaking of
which, I applaud the choice of this question, because it looks good on screen.
The viewer enjoys intricate and intriguing drawings of trees, and is reminded
that mathematics is not numbers plus algebraic gobbledygook.

== Here be dragons ==

While finding the above 10 trees is easy, there are deep, difficult questions
lurking only a hair's breadth away. For starters, what if we generalize and ask
for the number of irreducible trees with \(n\) nodes? It seems there is no neat
formula, and we must settle for link:genfun.html[an arduous trek through
combinatorics to arrive at a generating function whose coefficients we can
compute]. When doing so, we discover this is also the case for counting
unlabled trees.

What about writing code that determines if two graphs are equivalent? This
results in an algorithm for the _graph isomorphism problem_, which lies in NP
because we can succinctly describe how to relabel vertices to prove two graphs
are the same. But is there an quick way to find such a relabeling? In other
words, is the graph isomorphism problem in P? Nobody knows.

Murkier still, nobody knows if graph isomorphism is NP-complete (roughly: no
easier than any other NP problem). Thus graph isomorphism is a rare problem
that lies somewhere between P and NP-complete, but no one has proved it's one
or the other. Or strictly in between, assuming P ≠ NP. In other words,
graph isomorphism might be a member of the exclusive club known as
https://en.wikipedia.org/wiki/NP-intermediate[NP-intermediate].

On the one hand, we can generalize slightly to obtain a known NP-complete
problem, the _subgraph isomorphism problem_, which asks if a given graph is
equivalent to some subgraph of another graph. On the other hand, Babai
published (then retracted then restored) a paper showing that graph isomorphism
can be solved in quasi-polynomial time; the best known algorithms for solving
NP-complete problems are exponential-time.

Cryptography mysteriously appears again. Other famous suspected members of
NP-intermediate are integer factorization and discrete log. Computer security
relies on these problems being difficult to solve. To the dismay of those who
build cryptosystems, while nobody can solve them in polynomial time (without a
quantum computer that is perpetually ten years away), there exist
subexponential-time algorithms.