Repository: gonum/gonum
Pull request #2000 (2024-11-02 13:55:48)

graph: add algorithm for finding intervals in flow graphs
Hi,

This is the code for issue #796 allowing finding intervals in flow graphs. I've added some basic testing using an example graph given in C. Cifuentes, "Reverse Compilation Techniques", 1994. This is my first time contributing to an open source project and I look forward to your feedback!

Thanks,
Jack


harpend: Added the changes suggested

Also switched IntervalGraph to use to instead of from so the nodes to an interval header node can be accessed.
This is accessed through the To function.
kortschak: I'd like to go over the tests and make them more table driven, but I have not yet had time.
kortschak: I've gone through the tests to make them accept new test vectors when the need arises, and to more fully test the vectors. This raised issues for me with the semantics here.

The `Interval` type implements `graph.Graph`, but it looks to me like it is formally a directed graph. It would be nice to include this in the type if this is not too much work/weight.

In order to get the more rigorous testing to work and the `Interval` type to work with the rest of the graph packages, I needed to make the `IntervalGraph` hold and return `*Interval`. I'm not particularly fussed by this.

```diff
diff --git a/graph/flow/interval.go b/graph/flow/interval.go
index 5dcda8c7..73ceb1ed 100644
--- a/graph/flow/interval.go
+++ b/graph/flow/interval.go
@@ -48,7 +48,7 @@ terminates
 // of the directed graph in g from starting from the entry node identified
 // by eid.
 func Intervals(g graph.Directed, eid int64) IntervalGraph {
-       var intervals []Interval
+       var intervals []*Interval
 
        var ns linear.NodeStack
        visited := make(map[int64]bool)
@@ -62,7 +62,7 @@ func Intervals(g graph.Directed, eid int64) IntervalGraph {
        for worklist.Len() != 0 {
                var interval Interval
                maps.Copy(inInterval, interval.findInterval(worklist.Dequeue(), g))
-               intervals = append(intervals, interval)
+               intervals = append(intervals, &interval)
 
                for _, node := range ns {
                        if inInterval[node.ID()] != nil {
@@ -90,7 +90,7 @@ func Intervals(g graph.Directed, eid int64) IntervalGraph {
 
 // IntervalGraph contains the intervals and the edges between the intervals.
 type IntervalGraph struct {
-       Intervals []Interval
+       Intervals []*Interval
        headers   map[int64]graph.Node
        to        map[int64]map[int64]graph.Edge
 }
@@ -110,7 +110,7 @@ func (ig *IntervalGraph) To(id int64) graph.Nodes {
 type Interval struct {
        head  graph.Node
        nodes map[int64]graph.Node
-       from  map[int64]map[int64]graph.Edge
+       edges map[int64]map[int64]graph.Edge
 }
 
 // Head returns header node for an interval.
@@ -129,7 +129,7 @@ func (i *Interval) Nodes() graph.Nodes {
 // Edge returns the edge given 2 node id's if the edge exists.
 // Else it returns null.
 func (i *Interval) Edge(uid, vid int64) graph.Edge {
-       edge, ok := i.from[uid][vid]
+       edge, ok := i.edges[uid][vid]
        if !ok {
                return nil
        }
@@ -140,19 +140,19 @@ func (i *Interval) Edge(uid, vid int64) graph.Edge {
 // The returned graph.Nodes is only valid until the next mutation of
 // the receiver.
 func (i *Interval) From(id int64) graph.Nodes {
-       if len(i.from[id]) == 0 {
+       if len(i.edges[id]) == 0 {
                return graph.Empty
        }
-       return iterator.NewNodesByEdge(i.nodes, i.from[id])
+       return iterator.NewNodesByEdge(i.nodes, i.edges[id])
 }
 
 // HasEdgeBetween returns whether an edge exists between nodes x and y without
 // considering direction.
 func (i *Interval) HasEdgeBetween(xid int64, yid int64) bool {
-       if _, ok := i.from[xid][yid]; ok {
+       if _, ok := i.edges[xid][yid]; ok {
                return true
        }
-       _, ok := i.from[yid][xid]
+       _, ok := i.edges[yid][xid]
        return ok
 }
 
@@ -219,7 +219,7 @@ func dfsPostorder(g graph.Directed, eid int64, ns *linear.NodeStack, visited map
 }
 
 // linkIntervals computes the internal and external edges for the intervals.
-func linkIntervals(intervals []Interval, g graph.Directed) IntervalGraph {
+func linkIntervals(intervals []*Interval, g graph.Directed) IntervalGraph {
        ig := IntervalGraph{
                Intervals: intervals,
                headers:   make(map[int64]graph.Node),
@@ -228,17 +228,17 @@ func linkIntervals(intervals []Interval, g graph.Directed) IntervalGraph {
 
        for i, interval := range intervals {
                ig.headers[int64(i)] = interval.head
-               interval.from = make(map[int64]map[int64]graph.Edge)
+               interval.edges = make(map[int64]map[int64]graph.Edge)
                for _, node := range interval.nodes {
                        succs := g.From(node.ID())
                        for succs.Next() {
                                succNode := succs.Node()
                                if interval.nodes[succNode.ID()] != nil {
-                                       if interval.from[node.ID()] == nil {
-                                               interval.from[node.ID()] = make(map[int64]graph.Edge)
+                                       if interval.edges[node.ID()] == nil {
+                                               interval.edges[node.ID()] = make(map[int64]graph.Edge)
                                        }
 
-                                       interval.from[node.ID()][succNode.ID()] = g.Edge(node.ID(), succNode.ID())
+                                       interval.edges[node.ID()][succNode.ID()] = g.Edge(node.ID(), succNode.ID())
                                } else {
                                        if ig.to[succNode.ID()] == nil {
                                                ig.to[succNode.ID()] = make(map[int64]graph.Edge)
diff --git a/graph/flow/interval_test.go b/graph/flow/interval_test.go
index acfc96c5..692ec2b7 100644
--- a/graph/flow/interval_test.go
+++ b/graph/flow/interval_test.go
@@ -7,77 +7,127 @@ package flow
 import (
        "testing"
 
+       "gonum.org/v1/gonum/graph"
+       "gonum.org/v1/gonum/graph/encoding/dot"
        "gonum.org/v1/gonum/graph/simple"
+       "gonum.org/v1/gonum/graph/topo"
 )
 
-func TestInterval(t *testing.T) {
-       // Graph from C. Cifuentes, "Reverse Compilation Techniques", 1994 (figure 6-23)
-       g := simple.NewDirectedGraph()
-       edges := []struct{ from, to int }{
-               {1, 2},
-               {1, 5},
-               {2, 3},
-               {2, 4},
-               {3, 5},
-               {4, 5},
-               {5, 6},
-               {6, 7},
-               {6, 12},
-               {7, 8},
-               {7, 9},
-               {8, 9},
-               {8, 10},
-               {9, 10},
-               {10, 11},
-               {12, 13},
-               {14, 13},
-               {13, 14},
-               {14, 15},
-               {15, 6},
-       }
-
-       nodes := [][]int64{
-               {1, 2, 3, 4, 5},
-               {6, 7, 8, 9, 10, 11, 12},
-               {13, 14, 15},
-       }
+var intervalTests = []struct {
+       name  string
+       edges []struct{ from, to simple.Node }
+       root  int64
+       want  []map[int64]intset
+}{
+       {
+               // Graph from C. Cifuentes, "Reverse Compilation Techniques", 1994 (figure 6-23)
+               // Available via https://eprints.qut.edu.au/36820/
+               name: "cifuentes",
+               edges: []struct{ from, to simple.Node }{
+                       {1, 2},
+                       {1, 5},
+                       {2, 3},
+                       {2, 4},
+                       {3, 5},
+                       {4, 5},
+                       {5, 6},
+                       {6, 7},
+                       {6, 12},
+                       {7, 8},
+                       {7, 9},
+                       {8, 9},
+                       {8, 10},
+                       {9, 10},
+                       {10, 11},
+                       {12, 13},
+                       {14, 13},
+                       {13, 14},
+                       {14, 15},
+                       {15, 6},
+               },
+               root: 1,
+               want: []map[int64]intset{
+                       {
+                               1: linksTo(2, 5),
+                               2: linksTo(3, 4),
+                               3: linksTo(5),
+                               4: linksTo(5),
+                               5: nil,
+                       },
+                       {
+                               6:  linksTo(7, 12),
+                               7:  linksTo(8, 9),
+                               8:  linksTo(9, 10),
+                               9:  linksTo(10),
+                               10: linksTo(11),
+                               11: nil,
+                               12: nil,
+                       },
+                       {
+                               13: linksTo(14),
+                               14: linksTo(13, 15),
+                               15: nil,
+                       },
+               },
+       },
+}
 
-       intervalSize := []int64{
-               5, 7, 3,
-       }
+func TestInterval(t *testing.T) {
+       for _, test := range intervalTests {
+               t.Run(test.name, func(t *testing.T) {
+                       g := simple.NewDirectedGraph()
+                       for _, e := range test.edges {
+                               g.SetEdge(simple.Edge{F: e.from, T: e.to})
+                       }
 
-       for _, e := range edges {
-               g.SetEdge(g.NewEdge(simple.Node(e.from), simple.Node(e.to)))
-       }
+                       ig := Intervals(g, test.root)
+                       if len(ig.Intervals) != len(test.want) {
+                               t.Fatalf("unexpected interval count: got:%d want:%d", len(ig.Intervals), len(test.want))
+                       }
 
-       // Test number of intervals
-       ig := Intervals(g, 1)
-       if len(ig.Intervals) != 3 {
-               t.Fatalf("Expected 3 intervals, got %d", len(ig.Intervals))
+                       for i, iv := range ig.Intervals {
+                               var got graph.Graph = iv
+                               want := gFromIntsets(test.want[i])
+                               if !topo.Equal(got, want) {
+                                       gotDot, err := dot.Marshal(iv, "", "", "\t")
+                                       if err != nil {
+                                               t.Fatalf("unexpected error marshalling got DOT: %v", err)
+                                       }
+                                       wantDot, err := dot.Marshal(want, "", "", "\t")
+                                       if err != nil {
+                                               t.Fatalf("unexpected error marshalling want DOT: %v", err)
+                                       }
+                                       t.Errorf("unexpected topology of interval %d:\ngot:\n%s\nwant:\n%s", i, gotDot, wantDot)
+                               }
+                       }
+               })
        }
+}
 
-       for i, interval := range ig.Intervals {
-               // Test correct number of nodes are found in interval
-               if len(interval.nodes) != int(intervalSize[i]) {
-                       t.Errorf("Expected %d nodes in interval %d, got %d", intervalSize[i], i, len(interval.nodes))
+func gFromIntsets(s map[int64]intset) graph.Directed {
+       g := simple.NewDirectedGraph()
+       for u, e := range s {
+               // Add nodes that are not defined by an edge.
+               if g.Node(int64(u)) == nil {
+                       g.AddNode(simple.Node(u))
                }
-
-               // Test all expected in interval are present
-               for _, node := range nodes[i] {
-                       if interval.nodes[node] == nil {
-                               t.Errorf("Node %d not found in interval %d", node, i)
-                       }
+               for v := range e {
+                       g.SetEdge(simple.Edge{F: simple.Node(u), T: simple.Node(v)})
                }
        }
+       return g
+}
 
-       // test interval edges
-       toEdges := ig.To(6)
-       if toEdges.Len() != 2 {
-               t.Errorf("Expected 2 edges to node 6 in interval graph, got %d", len(ig.to))
-       }
+// intset is an integer set.
+type intset map[int64]struct{}
 
-       toEdges2 := ig.To(13)
-       if toEdges2.Len() != 1 {
-               t.Errorf("Expected 1 edge to node 13 in interval graph, got %d", len(ig.to))
+func linksTo(i ...int64) intset {
+       if len(i) == 0 {
+               return nil
+       }
+       s := make(intset)
+       for _, v := range i {
+               s[v] = struct{}{}
        }
+       return s
 }
```
harpend: Interval is updated to be a directed graph, is there any desire for interval graph to be a directed graph as well?
kortschak: I think that would make sense.
harpend: added tests for edges between intervals and IntervalGraph now implements directed graph
mewmew: With the updated tests and the change to have interval implement the directed graph interface, I think the PR looks good.

@kortschak has the last say before merging : )

Really happy to see this finally being added to Gonum. Awesome work @harpend!

Cheerful regards,
Robin
harpend: > I'll take a look on the weekend. Are you able to send the A+C before then? If so, I can merge then as well.

A+C PR is already opened.
harpend: The email in the header I only use for github, and do not wish to include a personal email.


kortschak: OK

Patch:
diff --git a/graph/flow/interval.go b/graph/flow/interval.go
new file mode 100644
index 00000000..930a9b62
--- /dev/null
+++ b/graph/flow/interval.go
@@ -0,0 +1,385 @@
+// Copyright ©2024 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package flow
+
+import (
+	"maps"
+	"slices"
+
+	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/internal/linear"
+	"gonum.org/v1/gonum/graph/iterator"
+	"gonum.org/v1/gonum/graph/simple"
+)
+
+/*
+Intervals are defined as the maximal, single entry subgraph for which h (head)
+is the entry node and in which all closed paths contain h.
+
+This is the algorithm given for computing them:
+
+1. Establish a set H for header nodes and initialize it with no, the
+unique entry node for the graph.
+
+2. For h in H find l(h) as follows:
+	2.1. Put h in l(h) as the first element of l(h).
+
+	2.2. Add to l(h) any node all of whose immediate predecessors
+	are already in l(h).
+
+	2.3. Repeat 2.2 until no more nodes can be added to l(h).
+
+3. Add to H all nodes in G which are not already in H and which
+are not in l(h) but which have immediate predecessors in l(h).
+Therefore a node is added to H the first time any (but not all)
+of its immediate predecessors become members of an interval.
+
+4. Add l(h) to a set S of intervals being developed.
+
+5. Select the next unprocessed node in H and repeat steps 2, 3, 4, 5.
+When there are no more unprocessed nodes in H, the procedure
+terminates
+
+(https://dl.acm.org/doi/pdf/10.1145/360018.360025)
+*/
+
+// Intervals returns the IntervalGraph containing the individual intervals
+// of the directed graph in g when starting from the entry node identified
+// by eid.
+func Intervals(g graph.Directed, eid int64) IntervalGraph {
+	var intervals []*Interval
+
+	var ns linear.NodeStack
+	visited := make(map[int64]bool)
+
+	dfsPostorder(g, eid, &ns, visited)
+	slices.Reverse(ns)
+	var worklist linear.NodeQueue
+	worklist.Enqueue(g.Node(eid))
+	inInterval := make(map[int64]graph.Node)
+	node2interval := make(map[int64]*Interval)
+	id := int64(0)
+
+	for worklist.Len() != 0 {
+		var interval Interval
+		map1, map2 := interval.findInterval(worklist.Dequeue(), g)
+		maps.Copy(inInterval, map1)
+		maps.Copy(node2interval, map2)
+
+		interval.id = id
+		id++
+		intervals = append(intervals, &interval)
+
+		for _, node := range ns {
+			if inInterval[node.ID()] != nil {
+				continue
+			}
+
+			preds := g.To(node.ID())
+			predsLength := preds.Len()
+			x := 0
+			for preds.Next() {
+				if interval.nodes[preds.Node().ID()] != nil {
+					x++
+				}
+			}
+
+			if 0 < x && x < predsLength {
+				worklist.Enqueue(node)
+				break
+			}
+		}
+	}
+
+	return linkIntervals(intervals, g, node2interval)
+}
+
+// IntervalGraph contains the intervals and the edges between the intervals.
+type IntervalGraph struct {
+	Intervals map[int64]*Interval
+	nodes     map[int64]graph.Node
+	head      graph.Node
+	from      map[int64]map[int64]graph.Edge
+	to        map[int64]map[int64]graph.Edge
+}
+
+// Head returns header node for an IntervalGraph.
+func (ig *IntervalGraph) Head() graph.Node {
+	return ig.head
+}
+
+// Nodes returns a node iterator for an IntervalGraph.
+func (ig *IntervalGraph) Nodes() graph.Nodes {
+	if len(ig.nodes) == 0 {
+		return graph.Empty
+	}
+	return iterator.NewNodes(ig.nodes)
+}
+
+// Edge returns the edge given 2 node id's if the edge exists.
+// Else it returns null.
+func (ig *IntervalGraph) Edge(uid, vid int64) graph.Edge {
+	edge, ok := ig.from[uid][vid]
+	if !ok {
+		return nil
+	}
+	return edge
+}
+
+// From returns all nodes in g that can be reached directly from n.
+//
+// The returned graph.Nodes is only valid until the next mutation of
+// the receiver.
+func (ig *IntervalGraph) From(id int64) graph.Nodes {
+	if len(ig.from[id]) == 0 {
+		return graph.Empty
+	}
+	return iterator.NewNodesByEdge(ig.nodes, ig.from[id])
+}
+
+// To returns all nodes in g that can reach directly to n.
+//
+// The returned graph.Nodes is only valid until the next mutation of
+// the receiver.
+func (ig *IntervalGraph) To(id int64) graph.Nodes {
+	if len(ig.to[id]) == 0 {
+		return graph.Empty
+	}
+	return iterator.NewNodesByEdge(ig.nodes, ig.to[id])
+}
+
+// HasEdgeBetween returns whether an edge exists between nodes x and y without
+// considering direction.
+func (ig *IntervalGraph) HasEdgeBetween(xid int64, yid int64) bool {
+	if _, ok := ig.from[xid][yid]; ok {
+		return true
+	}
+	_, ok := ig.from[yid][xid]
+	return ok
+}
+
+// Node returns the node with the given ID if it exists in the IntervalGraph,
+// and nil otherwise.
+func (ig *IntervalGraph) Node(id int64) graph.Node {
+	return ig.nodes[id]
+}
+
+// HasEdgeFromTo returns whether an edge exists in the graph from u to v.
+func (ig *IntervalGraph) HasEdgeFromTo(uid, vid int64) bool {
+	if _, ok := ig.from[uid][vid]; !ok {
+		return false
+	}
+	return true
+}
+
+// setEdge adds e, an edge from one node to another.
+func (ig *IntervalGraph) setEdge(e graph.Edge) {
+	var (
+		from = e.From()
+		fid  = from.ID()
+		to   = e.To()
+		tid  = to.ID()
+	)
+
+	ig.nodes[fid] = from
+	ig.nodes[tid] = to
+	ig.from[fid] = map[int64]graph.Edge{tid: e}
+	ig.to[tid] = map[int64]graph.Edge{fid: e}
+}
+
+// Interval I(h) is the maximal, single entry subgraph for which h (head)
+// is the entry node and in which all closed paths contain h.
+type Interval struct {
+	head  graph.Node
+	id    int64
+	nodes map[int64]graph.Node
+	from  map[int64]map[int64]graph.Edge
+	to    map[int64]map[int64]graph.Edge
+}
+
+// ID returns the ID number of the node.
+func (i *Interval) ID() int64 {
+	return i.id
+}
+
+// Head returns header node for an interval.
+func (i *Interval) Head() graph.Node {
+	return i.head
+}
+
+// Nodes returns a node iterator for an interval.
+func (i *Interval) Nodes() graph.Nodes {
+	if len(i.nodes) == 0 {
+		return graph.Empty
+	}
+	return iterator.NewNodes(i.nodes)
+}
+
+// Edge returns the edge given 2 node id's if the edge exists.
+// Else it returns null.
+func (i *Interval) Edge(uid, vid int64) graph.Edge {
+	edge, ok := i.from[uid][vid]
+	if !ok {
+		return nil
+	}
+	return edge
+}
+
+// From returns all nodes in g that can be reached directly from n.
+//
+// The returned graph.Nodes is only valid until the next mutation of
+// the receiver.
+func (i *Interval) From(id int64) graph.Nodes {
+	if len(i.from[id]) == 0 {
+		return graph.Empty
+	}
+	return iterator.NewNodesByEdge(i.nodes, i.from[id])
+}
+
+// To returns all nodes in g that can reach directly to n.
+//
+// The returned graph.Nodes is only valid until the next mutation of
+// the receiver.
+func (i *Interval) To(id int64) graph.Nodes {
+	if len(i.to[id]) == 0 {
+		return graph.Empty
+	}
+	return iterator.NewNodesByEdge(i.nodes, i.to[id])
+}
+
+// HasEdgeBetween returns whether an edge exists between nodes x and y without
+// considering direction.
+func (i *Interval) HasEdgeBetween(xid int64, yid int64) bool {
+	if _, ok := i.from[xid][yid]; ok {
+		return true
+	}
+	_, ok := i.from[yid][xid]
+	return ok
+}
+
+// Node returns the node with the given ID if it exists in the interval,
+// and nil otherwise.
+func (i *Interval) Node(id int64) graph.Node {
+	return i.nodes[id]
+}
+
+// HasEdgeFromTo returns whether an edge exists in the graph from u to v.
+func (i *Interval) HasEdgeFromTo(uid, vid int64) bool {
+	if _, ok := i.from[uid][vid]; !ok {
+		return false
+	}
+	return true
+}
+
+// findInterval finds all interval nodes.
+// Nodes are added to the interval if all their predecessors are in
+// the interval or they are the header node.
+func (i *Interval) findInterval(h graph.Node, g graph.Directed) (map[int64]graph.Node, map[int64]*Interval) {
+	i.head = h
+	var nq linear.NodeQueue
+	nq.Enqueue(h)
+	node2interval := make(map[int64]*Interval)
+	node2interval[h.ID()] = i
+	i.nodes = make(map[int64]graph.Node)
+	i.nodes[h.ID()] = h
+	var node graph.Node
+	for nq.Len() > 0 {
+		node = nq.Dequeue()
+		succs := g.From(node.ID())
+
+		for succs.Next() {
+			if i.nodes[succs.Node().ID()] != nil {
+				continue
+			}
+
+			preds := g.To(succs.Node().ID())
+			predsLength := preds.Len()
+			x := 0
+			for preds.Next() {
+				if i.nodes[preds.Node().ID()] != nil {
+					x++
+				}
+			}
+
+			if x == predsLength {
+				nq.Enqueue(succs.Node())
+				i.nodes[succs.Node().ID()] = succs.Node()
+				node2interval[succs.Node().ID()] = i
+			}
+		}
+	}
+
+	return i.nodes, node2interval
+}
+
+// dfsPosterorder puts nodes into a stack in postorder.
+func dfsPostorder(g graph.Directed, eid int64, ns *linear.NodeStack, visited map[int64]bool) {
+	succs := g.From(eid)
+	visited[eid] = true
+	for {
+		if !succs.Next() {
+			break
+		}
+
+		succ := succs.Node()
+		if !visited[succ.ID()] {
+			dfsPostorder(g, succ.ID(), ns, visited)
+		}
+	}
+
+	ns.Push(g.Node(eid))
+}
+
+// linkIntervals computes the internal and external edges for the intervals.
+func linkIntervals(intervals []*Interval, g graph.Directed, node2interval map[int64]*Interval) IntervalGraph {
+	ig := IntervalGraph{
+		Intervals: make(map[int64]*Interval),
+		nodes:     make(map[int64]graph.Node),
+		head:      intervals[0].head,
+		from:      make(map[int64]map[int64]graph.Edge),
+		to:        make(map[int64]map[int64]graph.Edge),
+	}
+
+	for i, interval := range intervals {
+		ig.nodes[int64(i)] = interval
+		ig.Intervals[int64(i)] = interval
+
+		interval.from = make(map[int64]map[int64]graph.Edge)
+		interval.to = make(map[int64]map[int64]graph.Edge)
+		for _, node := range interval.nodes {
+			succs := g.From(node.ID())
+			for succs.Next() {
+				succNode := succs.Node()
+				if interval.nodes[succNode.ID()] != nil {
+					// edge contained within interval
+					if interval.from[node.ID()] == nil {
+						interval.from[node.ID()] = make(map[int64]graph.Edge)
+					}
+
+					if interval.to[succNode.ID()] == nil {
+						interval.to[succNode.ID()] = make(map[int64]graph.Edge)
+					}
+
+					interval.from[node.ID()][succNode.ID()] = g.Edge(node.ID(), succNode.ID())
+					interval.to[succNode.ID()][node.ID()] = g.Edge(node.ID(), succNode.ID())
+				} else {
+					// edge from one interval to another interval
+					succInterval := node2interval[succNode.ID()]
+					if ig.to[succInterval.ID()] == nil {
+						ig.to[succInterval.ID()] = make(map[int64]graph.Edge)
+					}
+
+					if ig.from[interval.ID()] == nil {
+						ig.from[interval.ID()] = make(map[int64]graph.Edge)
+					}
+
+					ig.setEdge(simple.Edge{F: interval, T: succInterval})
+				}
+			}
+		}
+	}
+
+	return ig
+}


Test patch:
diff --git a/graph/flow/interval_test.go b/graph/flow/interval_test.go
new file mode 100644
index 00000000..6800d73a
--- /dev/null
+++ b/graph/flow/interval_test.go
@@ -0,0 +1,160 @@
+// Copyright ©2024 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package flow
+
+import (
+	"testing"
+
+	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/encoding/dot"
+	"gonum.org/v1/gonum/graph/simple"
+	"gonum.org/v1/gonum/graph/topo"
+)
+
+var intervalTests = []struct {
+	name          string
+	internalEdges []struct{ from, to simple.Node }
+	externalEdges []struct{ from, to simple.Node }
+	root          int64
+	internalWant  []map[int64]intset
+	externalWant  map[int64]intset
+}{
+	{
+		// Graph from C. Cifuentes, "Reverse Compilation Techniques", 1994 (figure 6-23)
+		// Available via https://eprints.qut.edu.au/36820/
+		name: "cifuentes",
+		internalEdges: []struct{ from, to simple.Node }{
+			{1, 2},
+			{1, 5},
+			{2, 3},
+			{2, 4},
+			{3, 5},
+			{4, 5},
+			{5, 6},
+			{6, 7},
+			{6, 12},
+			{7, 8},
+			{7, 9},
+			{8, 9},
+			{8, 10},
+			{9, 10},
+			{10, 11},
+			{12, 13},
+			{14, 13},
+			{13, 14},
+			{14, 15},
+			{15, 6},
+		},
+		externalEdges: []struct{ from, to simple.Node }{
+			{0, 1},
+			{1, 2},
+			{2, 1},
+		},
+		root: 1,
+		internalWant: []map[int64]intset{
+			{
+				1: linksTo(2, 5),
+				2: linksTo(3, 4),
+				3: linksTo(5),
+				4: linksTo(5),
+				5: nil,
+			},
+			{
+				6:  linksTo(7, 12),
+				7:  linksTo(8, 9),
+				8:  linksTo(9, 10),
+				9:  linksTo(10),
+				10: linksTo(11),
+				11: nil,
+				12: nil,
+			},
+			{
+				13: linksTo(14),
+				14: linksTo(13, 15),
+				15: nil,
+			},
+		},
+		externalWant: map[int64]intset{
+			0: linksTo(1),
+			1: linksTo(2),
+			2: linksTo(1),
+		},
+	},
+}
+
+func TestInterval(t *testing.T) {
+	for _, test := range intervalTests {
+		t.Run(test.name, func(t *testing.T) {
+			g := simple.NewDirectedGraph()
+			for _, e := range test.internalEdges {
+				g.SetEdge(simple.Edge{F: e.from, T: e.to})
+			}
+
+			ig := Intervals(g, test.root)
+			if len(ig.Intervals) != len(test.internalWant) {
+				t.Fatalf("unexpected interval count: got:%d want:%d", len(ig.Intervals), len(test.internalWant))
+			}
+
+			var igGot graph.Directed = &ig
+			igWant := gFromIntsets(test.externalWant)
+			if !topo.Equal(igGot, igWant) {
+				igGotDot, err := dot.Marshal(&ig, "", "", "\t")
+				if err != nil {
+					t.Fatalf("unexpected error marshalling got DOT: %v", err)
+				}
+				igWantDot, err := dot.Marshal(igWant, "", "", "\t")
+				if err != nil {
+					t.Fatalf("unexpected error marshalling want DOT: %v", err)
+				}
+				t.Errorf("unexpected topology of interval graph:\ngot:\n%s\nwant:\n%s", igGotDot, igWantDot)
+			}
+
+			for i, iv := range ig.Intervals {
+				var got graph.Directed = iv
+				iWant := gFromIntsets(test.internalWant[i])
+				if !topo.Equal(got, iWant) {
+					iGotDot, err := dot.Marshal(iv, "", "", "\t")
+					if err != nil {
+						t.Fatalf("unexpected error marshalling got DOT: %v", err)
+					}
+					iWantDot, err := dot.Marshal(iWant, "", "", "\t")
+					if err != nil {
+						t.Fatalf("unexpected error marshalling want DOT: %v", err)
+					}
+					t.Errorf("unexpected topology of interval %d:\ngot:\n%s\nwant:\n%s", i, iGotDot, iWantDot)
+				}
+			}
+		})
+	}
+}
+
+func gFromIntsets(s map[int64]intset) graph.Directed {
+	g := simple.NewDirectedGraph()
+	for u, e := range s {
+		// Add nodes that are not defined by an edge.
+		if g.Node(int64(u)) == nil {
+			g.AddNode(simple.Node(u))
+		}
+
+		for v := range e {
+			g.SetEdge(simple.Edge{F: simple.Node(u), T: simple.Node(v)})
+		}
+	}
+	return g
+}
+
+// intset is an integer set.
+type intset map[int64]struct{}
+
+func linksTo(i ...int64) intset {
+	if len(i) == 0 {
+		return nil
+	}
+	s := make(intset)
+	for _, v := range i {
+		s[v] = struct{}{}
+	}
+	return s
+}