Repository: gonum/gonum
Pull request #1615 (2021-04-23 02:54:20)

graph/coloring: new package providing graph coloring functions
Please take a look.

/cc @mattConn 

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->


codecov-commenter: # [Codecov](https://codecov.io/gh/gonum/gonum/pull/1615?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) Report
> Merging [#1615](https://codecov.io/gh/gonum/gonum/pull/1615?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (aa50e37) into [master](https://codecov.io/gh/gonum/gonum/commit/ccc0da248a78266c0e8a38ac56e983d17c4ba15f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (ccc0da2) will **decrease** coverage by `0.12%`.
> The diff coverage is `83.92%`.

[![Impacted file tree graph](https://codecov.io/gh/gonum/gonum/pull/1615/graphs/tree.svg?width=650&height=150&src=pr&token=B7qHnqSPNL&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)](https://codecov.io/gh/gonum/gonum/pull/1615?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)

```diff
@@            Coverage Diff             @@
##           master    #1615      +/-   ##
==========================================
- Coverage   73.04%   72.91%   -0.13%     
==========================================
  Files         498      483      -15     
  Lines       58813    58444     -369     
==========================================
- Hits        42961    42617     -344     
+ Misses      13424    13417       -7     
+ Partials     2428     2410      -18     
```


| [Impacted Files](https://codecov.io/gh/gonum/gonum/pull/1615?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | Coverage Δ | |
|---|---|---|
| [graph/coloring/coloring.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvY29sb3JpbmcvY29sb3JpbmcuZ28=) | `83.92% <83.92%> (ø)` | |
| [internal/asm/c64/conj.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2M2NC9jb25qLmdv) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [optimize/listsearch.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-b3B0aW1pemUvbGlzdHNlYXJjaC5nbw==) | `72.22% <0.00%> (-1.86%)` | :arrow_down: |
| [graph/iterator/nodes\_map\_safe.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-Z3JhcGgvaXRlcmF0b3Ivbm9kZXNfbWFwX3NhZmUuZ28=) | | |
| [mat/offset\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L29mZnNldF9hcHBlbmdpbmUuZ28=) | | |
| [internal/asm/f32/ge\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2YzMi9nZV9ub2FzbS5nbw==) | | |
| [internal/asm/f64/l2norm\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9sMm5vcm1fbm9hc20uZ28=) | | |
| [internal/asm/f64/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9zdHVic19ub2FzbS5nbw==) | | |
| [internal/asm/f64/axpy.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9heHB5Lmdv) | | |
| [internal/asm/c128/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2MxMjgvc3R1YnNfbm9hc20uZ28=) | | |
| ... and [10 more](https://codecov.io/gh/gonum/gonum/pull/1615/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gonum/gonum/pull/1615?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gonum/gonum/pull/1615?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Last update [ccc0da2...aa50e37](https://codecov.io/gh/gonum/gonum/pull/1615?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).

mattConn: Looks interesting. I'll wrap up my review of `urna.go` in commit 3 of #1524 and then I will look through this.
kortschak: This is ready. There is another solver that is worth including, but it needs more investigation.

The golangci-lint failure can be ignored.
mattConn: I'll look at this after #1617 since that one is smaller and has a formal request.
mattConn: Just for an update, I am actively reviewing this now.
kortschak: Done
kortschak: The `DsaturExact` code is quite slow at the moment (unreasonably so). I've looked into this and it comes down to the default approach of using maps to hold edge and node lists, and colourings. After this is merged, I'll file an issue to make an internal representation based on bitmaps for edge and node lists within graphs and slices for colourings. The optimisation is likely to make the code considerably less readable, so I want to get the base code in first.

I also have a novel algorithmic optimisation for `DsaturExact` that I've been exploring. That is ready to go as soon as this is merged and gives quite impressive speed-ups for difficult graphs. The two approaches together should make this implementation competitive.
mattConn: For some kind of an update since it has been a while, I am almost done looking over `Dsatur` in `coloring.go` on the first commit. I am presently tracing this function and working on matching it with the steps outlined in Brélaz's text.
kortschak: Thanks for the review.

Patch:
diff --git a/graph/coloring/coloring.go b/graph/coloring/coloring.go
new file mode 100644
index 00000000..845bd2ba
--- /dev/null
+++ b/graph/coloring/coloring.go
@@ -0,0 +1,789 @@
+// Copyright ©2021 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.
+
+// See Lewis, A Guide to Graph Colouring: Algorithms and Applications
+// doi:10.1007/978-3-319-25730-3 for significant discussion of approaches.
+
+package coloring
+
+import (
+	"errors"
+	"sort"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/internal/set"
+	"gonum.org/v1/gonum/graph/iterator"
+	"gonum.org/v1/gonum/graph/topo"
+)
+
+// ErrInvalidPartialColoring is returned when a partial coloring
+// is provided for a graph with inadmissible color assignments.
+var ErrInvalidPartialColoring = errors.New("coloring: invalid partial coloring")
+
+// Sets returns the mapping from colors to sets of node IDs. Each set of
+// node IDs is sorted by ascending value.
+func Sets(colors map[int64]int) map[int][]int64 {
+	sets := make(map[int][]int64)
+	for id, c := range colors {
+		sets[c] = append(sets[c], id)
+	}
+	for _, s := range sets {
+		sort.Slice(s, func(i, j int) bool { return s[i] < s[j] })
+	}
+	return sets
+}
+
+// Dsatur returns an approximate minimal chromatic number of g and a
+// corresponding vertex coloring using the heuristic Dsatur coloring algorithm.
+// If a partial coloring is provided the coloring will be consistent with
+// that partial coloring if possible. Otherwise Dsatur will return
+// ErrInvalidPartialColoring.
+// See Brélaz doi:10.1145/359094.359101 for details of the algorithm.
+func Dsatur(g graph.Undirected, partial map[int64]int) (k int, colors map[int64]int, err error) {
+	nodes := g.Nodes()
+	n := nodes.Len()
+	if n == 0 {
+		return
+	}
+	partial, ok := newPartial(partial, g)
+	if !ok {
+		return -1, nil, ErrInvalidPartialColoring
+	}
+	order := bySaturationDegree(nodes, g, partial)
+	order.heuristic = order.dsatur
+	k, colors = greedyColoringOf(g, order, order.colors)
+	return k, colors, nil
+}
+
+// Terminator is a cancellation-only context type. A context.Context
+// may be used as a Terminator.
+type Terminator interface {
+	// Done returns a channel that is closed when work
+	// should be terminated. Done may return nil if this
+	// work can never be canceled.
+	// Successive calls to Done should return the same value.
+	Done() <-chan struct{}
+
+	// If Done is not yet closed, Err returns nil.
+	// If Done is closed, Err returns a non-nil error
+	// explaining why.
+	// After Err returns a non-nil error, successive
+	// calls to Err should return the same error.
+	Err() error
+}
+
+// DsaturExact returns the exact minimal chromatic number of g and a
+// corresponding vertex coloring using the branch-and-bound Dsatur coloring
+// algorithm of Brélaz. If the provided terminator is cancelled or times out
+// before completion, the terminator's reason for termination will be returned
+// along with a potentially sub-optimal chromatic number and coloring. If
+// term is nil, DsaturExact will run to completion.
+// See Brélaz doi:10.1145/359094.359101 for details of the algorithm.
+func DsaturExact(term Terminator, g graph.Undirected) (k int, colors map[int64]int, err error) {
+	// This is implemented essentially as described in algorithm 1 of
+	// doi:10.1002/net.21716 with the exception that we obtain a
+	// tighter upper bound by doing a single run of an approximate
+	// Brélaz Dsatur coloring, using the result if the recurrence is
+	// cancelled.
+
+	nodes := g.Nodes()
+	n := nodes.Len()
+	if n == 0 {
+		return
+	}
+
+	lb, maxClique := maximumClique(g)
+	if lb == n {
+		return lb, colorClique(maxClique), nil
+	}
+
+	order := bySaturationDegree(nodes, g, make(map[int64]int))
+	order.heuristic = order.dsatur
+	// Find initial coloring via Dsatur heuristic.
+	ub, initial := greedyColoringOf(g, order, order.colors)
+	if lb == ub {
+		return ub, initial, nil
+	}
+
+	selector := &order.saturationDegree
+	cand := newDsaturColoring(order.nodes, make(map[int64]int))
+	k, colors, err = dSaturExact(term, selector, cand, 0, ub, nil)
+	if colors == nil {
+		return ub, initial, err
+	}
+	if k == lb {
+		err = nil
+	}
+	return k, colors, err
+}
+
+// dSaturColoring is a partial graph coloring.
+type dSaturColoring struct {
+	colors    map[int64]int
+	uncolored set.Int64s
+}
+
+// newDsaturColoring returns a dSaturColoring representing a partial coloring
+// of a graph with the given nodes and colors.
+func newDsaturColoring(nodes []graph.Node, colors map[int64]int) dSaturColoring {
+	uncolored := make(set.Int64s)
+	for _, v := range nodes {
+		vid := v.ID()
+		if _, ok := colors[vid]; !ok {
+			uncolored.Add(vid)
+		}
+	}
+	return dSaturColoring{
+		colors:    colors,
+		uncolored: uncolored,
+	}
+}
+
+// color moves a node from the uncolored set to the colored set.
+func (c dSaturColoring) color(id int64) {
+	if !c.uncolored.Has(id) {
+		if _, ok := c.colors[id]; ok {
+			panic("coloring: coloring already colored node")
+		}
+		panic("coloring: coloring non-existent node")
+	}
+	// The node has its uncolored mark removed, but is
+	// not explicitly colored until the dSaturExact
+	// caller has completed its recursive exploration
+	// of the feasible colors.
+	c.uncolored.Remove(id)
+}
+
+// uncolor moves a node from the colored set to the uncolored set.
+func (c dSaturColoring) uncolor(id int64) {
+	if _, ok := c.colors[id]; !ok {
+		if c.uncolored.Has(id) {
+			panic("coloring: uncoloring already uncolored node")
+		}
+		panic("coloring: uncoloring non-existent node")
+	}
+	delete(c.colors, id)
+	c.uncolored.Add(id)
+}
+
+// dSaturExact recursively searches for an exact mimimum vertex coloring of the
+// full graph in cand. If no chromatic number lower than ub is found, colors is
+// returned as nil.
+func dSaturExact(term Terminator, selector *saturationDegree, cand dSaturColoring, k, ub int, best map[int64]int) (newK int, colors map[int64]int, err error) {
+	if len(cand.uncolored) == 0 {
+		// In the published algorithm, this is guarded by k < ub,
+		// but dSaturExact is never called with k >= ub; in the
+		// initial call we have excluded cases where k == ub and
+		// it cannot be greater, and in the recursive call, we
+		// have already checked that k < ub.
+		return k, clone(cand.colors), nil
+	}
+
+	if term != nil {
+		select {
+		case <-term.Done():
+			if best == nil {
+				return -1, nil, term.Err()
+			}
+			colors := make(set.Ints)
+			for _, c := range best {
+				colors.Add(c)
+			}
+			return colors.Count(), best, term.Err()
+		default:
+		}
+	}
+
+	// Select the next node.
+	selector.reset(cand.colors)
+	vid := selector.nodes[selector.dsatur()].ID()
+	cand.color(vid)
+	// If uncolor panics, we have failed to find a
+	// feasible color. This should never happen.
+	defer cand.uncolor(vid)
+
+	// Keep the adjacent colors set as it will be
+	// overwritten by child recurrences.
+	adjColors := selector.adjColors[selector.indexOf[vid]]
+
+	// Collect all feasible existing colors plus one, remembering it.
+	feasible := make(set.Ints)
+	for _, c := range cand.colors {
+		if adjColors.Has(c) {
+			continue
+		}
+		feasible.Add(c)
+	}
+	var newCol int
+	for c := 0; c < ub; c++ {
+		if feasible.Has(c) || adjColors.Has(c) {
+			continue
+		}
+		feasible.Add(c)
+		newCol = c
+		break
+	}
+
+	// Recur over every feasible color.
+	for c := range feasible {
+		cand.colors[vid] = c
+		effK := k
+		if c == newCol {
+			effK++
+		}
+		// In the published algorithm, the expression max(effK, lb) < ub is
+		// used, but lb < ub always since it is not updated and dSaturExact
+		// is not called if lb == ub, and it cannot be greater.
+		if effK < ub {
+			ub, best, err = dSaturExact(term, selector, cand, effK, ub, best)
+			if err != nil {
+				return ub, best, err
+			}
+		}
+	}
+
+	return ub, best, nil
+}
+
+// maximumClique returns a maximum clique in g and its order.
+func maximumClique(g graph.Undirected) (k int, maxClique []graph.Node) {
+	for _, c := range topo.BronKerbosch(g) {
+		if len(c) > len(maxClique) {
+			maxClique = c
+		}
+	}
+	return len(maxClique), maxClique
+}
+
+// colorClique returns a valid coloring for the given clique.
+func colorClique(clique []graph.Node) map[int64]int {
+	colors := make(map[int64]int, len(clique))
+	for c, u := range clique {
+		colors[u.ID()] = c
+	}
+	return colors
+}
+
+// Randomized returns an approximate minimal chromatic number of g and a
+// corresponding vertex coloring using a greedy coloring algorithm with a
+// random node ordering. If src is non-nil it will be used as the random
+// source, otherwise the global random source will be used. If a partial
+// coloring is provided the coloring will be consistent with that partial
+// coloring if possible. Otherwise Randomized will return
+// ErrInvalidPartialColoring.
+func Randomized(g graph.Undirected, partial map[int64]int, src rand.Source) (k int, colors map[int64]int, err error) {
+	nodes := g.Nodes()
+	n := nodes.Len()
+	if n == 0 {
+		return
+	}
+	partial, ok := newPartial(partial, g)
+	if !ok {
+		return -1, nil, ErrInvalidPartialColoring
+	}
+	k, colors = greedyColoringOf(g, randomize(nodes, src), partial)
+	return k, colors, nil
+}
+
+// randomize returns a graph.Node iterator that returns nodes in a random
+// order.
+func randomize(it graph.Nodes, src rand.Source) graph.Nodes {
+	nodes := graph.NodesOf(it)
+	var shuffle func(int, func(i, j int))
+	if src == nil {
+		shuffle = rand.Shuffle
+	} else {
+		shuffle = rand.New(src).Shuffle
+	}
+	shuffle(len(nodes), func(i, j int) {
+		nodes[i], nodes[j] = nodes[j], nodes[i]
+	})
+	return iterator.NewOrderedNodes(nodes)
+}
+
+// RecursiveLargestFirst returns an approximate minimal chromatic number
+// of g and a corresponding vertex coloring using the Recursive Largest
+// First coloring algorithm.
+// See Leighton doi:10.6028/jres.084.024 for details of the algorithm.
+func RecursiveLargestFirst(g graph.Undirected) (k int, colors map[int64]int) {
+	it := g.Nodes()
+	n := it.Len()
+	if n == 0 {
+		return
+	}
+	nodes := graph.NodesOf(it)
+	colors = make(map[int64]int)
+
+	// The names of variable here have been changed from the original PL-1
+	// for clarity, but the correspondence is as follows:
+	//  E   -> isolated
+	//  F   -> boundary
+	//  L   -> current
+	//  COL -> k
+	//  C   -> colors
+
+	// Initialize the boundary vector to the node degrees.
+	boundary := make([]int, len(nodes))
+	indexOf := make(map[int64]int)
+	for i, u := range nodes {
+		uid := u.ID()
+		indexOf[uid] = i
+		boundary[i] = g.From(uid).Len()
+	}
+	deleteFrom := func(vec []int, idx int) {
+		vec[idx] = -1
+		to := g.From(nodes[idx].ID())
+		for to.Next() {
+			vec[indexOf[to.Node().ID()]]--
+		}
+	}
+	isolated := make([]int, len(nodes))
+
+	// If there are any uncolored nodes, initiate the assignment of the next color.
+	// Incrementing color happens at the end of the loop in this implementation.
+	var current int
+	for j := 0; j < n; {
+		// Reinitialize the isolated vector.
+		copy(isolated, boundary)
+
+		// Select the node in U₁ with maximal degree in U₁.
+		for i := range nodes {
+			if boundary[i] > boundary[current] {
+				current = i
+			}
+		}
+
+		// Color the node just selected and continue to
+		// color nodes with color k until U₁ is empty.
+		for isolated[current] >= 0 {
+			// Color node and modify U₁ and U₂ accordingly.
+			deleteFrom(isolated, current)
+			deleteFrom(boundary, current)
+			colors[nodes[current].ID()] = k
+			j++
+			to := g.From(nodes[current].ID())
+			for to.Next() {
+				i := indexOf[to.Node().ID()]
+				if isolated[i] >= 0 {
+					deleteFrom(isolated, i)
+				}
+			}
+
+			// Find the first node in U₁, if any.
+			for i := range nodes {
+				if isolated[i] < 0 {
+					continue
+				}
+
+				// If U₁ is not empty, select the next node for coloring.
+				current = i
+				currAvail := boundary[current] - isolated[current]
+				for j := i; j < n; j++ {
+					if isolated[j] < 0 {
+						continue
+					}
+					nextAvail := boundary[j] - isolated[j]
+					switch {
+					case nextAvail > currAvail, nextAvail == currAvail && isolated[j] < isolated[current]:
+						current = j
+						currAvail = boundary[current] - isolated[current]
+					}
+				}
+				break
+			}
+		}
+
+		k++
+	}
+
+	return k, colors
+}
+
+// SanSegundo returns an approximate minimal chromatic number of g and a
+// corresponding vertex coloring using the PASS rule with a single run of a
+// greedy coloring algorithm. If a partial coloring is provided the coloring
+// will be consistent with that partial coloring if possible. Otherwise
+// SanSegundo will return ErrInvalidPartialColoring.
+// See San Segundo doi:10.1016/j.cor.2011.10.008 for details of the algorithm.
+func SanSegundo(g graph.Undirected, partial map[int64]int) (k int, colors map[int64]int, err error) {
+	nodes := g.Nodes()
+	n := nodes.Len()
+	if n == 0 {
+		return
+	}
+	partial, ok := newPartial(partial, g)
+	if !ok {
+		return -1, nil, ErrInvalidPartialColoring
+	}
+	order := bySaturationDegree(nodes, g, partial)
+	order.heuristic = order.pass
+	k, colors = greedyColoringOf(g, order, order.colors)
+	return k, colors, nil
+}
+
+// WelshPowell returns an approximate minimal chromatic number of g and a
+// corresponding vertex coloring using the Welsh and Powell coloring algorithm.
+// If a partial coloring is provided the coloring will be consistent with that
+// partial coloring if possible. Otherwise WelshPowell will return
+// ErrInvalidPartialColoring.
+// See Welsh and Powell doi:10.1093/comjnl/10.1.85 for details of the algorithm.
+func WelshPowell(g graph.Undirected, partial map[int64]int) (k int, colors map[int64]int, err error) {
+	nodes := g.Nodes()
+	n := nodes.Len()
+	if n == 0 {
+		return
+	}
+	partial, ok := newPartial(partial, g)
+	if !ok {
+		return -1, nil, ErrInvalidPartialColoring
+	}
+	k, colors = greedyColoringOf(g, byDescendingDegree(nodes, g), partial)
+	return k, colors, nil
+}
+
+// byDescendingDegree returns a graph.Node iterator that returns nodes
+// in order of descending degree.
+func byDescendingDegree(it graph.Nodes, g graph.Undirected) graph.Nodes {
+	nodes := graph.NodesOf(it)
+	n := byDescDegree{nodes: nodes, degrees: make([]int, len(nodes))}
+	for i, u := range nodes {
+		n.degrees[i] = g.From(u.ID()).Len()
+	}
+	sort.Sort(n)
+	return iterator.NewOrderedNodes(nodes)
+}
+
+// byDescDegree sorts a slice of graph.Node descending by the corresponding
+// value of the degrees slice.
+type byDescDegree struct {
+	nodes   []graph.Node
+	degrees []int
+}
+
+func (n byDescDegree) Len() int           { return len(n.nodes) }
+func (n byDescDegree) Less(i, j int) bool { return n.degrees[i] > n.degrees[j] }
+func (n byDescDegree) Swap(i, j int) {
+	n.nodes[i], n.nodes[j] = n.nodes[j], n.nodes[i]
+	n.degrees[i], n.degrees[j] = n.degrees[j], n.degrees[i]
+}
+
+// newPartial returns a new valid partial coloring is valid for g. An empty
+// partial coloring is valid. If the partial coloring is not valid, a nil map
+// is returned, otherwise a new non-nil map is returned. If the input partial
+// coloring is nil, a new map is created and returned.
+func newPartial(partial map[int64]int, g graph.Undirected) (map[int64]int, bool) {
+	if partial == nil {
+		return make(map[int64]int), true
+	}
+	for id, c := range partial {
+		if g.Node(id) == nil {
+			return nil, false
+		}
+		to := g.From(id)
+		for to.Next() {
+			if oc, ok := partial[to.Node().ID()]; ok && c == oc {
+				return nil, false
+			}
+		}
+	}
+	return clone(partial), true
+}
+
+func clone(colors map[int64]int) map[int64]int {
+	new := make(map[int64]int, len(colors))
+	for id, c := range colors {
+		new[id] = c
+	}
+	return new
+}
+
+// greedyColoringOf returns the chromatic number and a graph coloring of g
+// based on the sequential coloring of nodes given by order and starting from
+// the given partial coloring.
+func greedyColoringOf(g graph.Undirected, order graph.Nodes, partial map[int64]int) (k int, colors map[int64]int) {
+	colors = partial
+	constrained := false
+	for _, c := range colors {
+		if c > k {
+			k = c
+			constrained = true
+		}
+	}
+
+	// Next nodes are chosen by the specified heuristic in order.
+	for order.Next() {
+		uid := order.Node().ID()
+		used := colorsOf(g.From(uid), colors)
+		if c, ok := colors[uid]; ok {
+			if used.Has(c) {
+				return -1, nil
+			}
+			continue
+		}
+		// Color the chosen vertex with the least possible
+		// (lowest numbered) color.
+		for c := 0; c <= k+1; c++ {
+			if !used.Has(c) {
+				colors[uid] = c
+				if c > k {
+					k = c
+				}
+				break
+			}
+		}
+	}
+
+	if !constrained {
+		return k + 1, colors
+	}
+	seen := make(set.Ints)
+	for _, c := range colors {
+		seen.Add(c)
+	}
+	return seen.Count(), colors
+}
+
+// colorsOf returns all the colors in the coloring that are used by the
+// given nodes.
+func colorsOf(nodes graph.Nodes, coloring map[int64]int) set.Ints {
+	c := make(set.Ints, nodes.Len())
+	for nodes.Next() {
+		used, ok := coloring[nodes.Node().ID()]
+		if ok {
+			c.Add(used)
+		}
+	}
+	return c
+}
+
+// saturationDegreeIterator is a graph.Nodes iterator that returns nodes ordered
+// by decreasing saturation degree.
+type saturationDegreeIterator struct {
+	// cnt is the number of nodes that
+	// have been returned and curr is
+	// the current selection.
+	cnt, curr int
+
+	// heuristic determines the
+	// iterator's node selection
+	// heuristic. It can be either
+	// saturationDegree.dsatur or
+	// saturationDegree.pass.
+	heuristic func() int
+
+	saturationDegree
+}
+
+// bySaturationDegree returns a new saturationDegreeIterator that will
+// iterate over the node in it based on the given graph and partial coloring.
+// The saturationDegreeIterator holds a reference to colors allowing
+// greedyColoringOf to update its coloring.
+func bySaturationDegree(it graph.Nodes, g graph.Undirected, colors map[int64]int) *saturationDegreeIterator {
+	return &saturationDegreeIterator{
+		cnt: -1, curr: -1,
+		saturationDegree: newSaturationDegree(it, g, colors),
+	}
+}
+
+// Len returns the number of elements remaining in the iterator.
+// saturationDegreeIterator is an indeterminate iterator, so Len always
+// returns -1. This is required to satisfy the graph.Iterator interface.
+func (n *saturationDegreeIterator) Len() int { return -1 }
+
+// Next advances the iterator to the next node and returns whether
+// the next call to the Node method will return a valid Node.
+func (n *saturationDegreeIterator) Next() bool {
+	if uint(n.cnt)+1 < uint(len(n.nodes)) {
+		n.cnt++
+		switch n.cnt {
+		case 0:
+			max := -1
+			for i, d := range n.degrees {
+				if d > max {
+					max = d
+					n.curr = i
+				}
+			}
+		default:
+			prev := n.Node().ID()
+			c := n.colors[prev]
+			to := n.g.From(prev)
+			for to.Next() {
+				n.adjColors[n.indexOf[to.Node().ID()]].Add(c)
+			}
+
+			chosen := n.heuristic()
+			if chosen < 0 || chosen == n.curr {
+				return false
+			}
+			n.curr = chosen
+		}
+		return true
+	}
+	n.cnt = len(n.nodes)
+	return false
+}
+
+// Node returns the current node.
+func (n *saturationDegreeIterator) Node() graph.Node { return n.nodes[n.curr] }
+
+// Reset implements the graph.Iterator interface. It should not be called.
+func (n *saturationDegreeIterator) Reset() { panic("coloring: invalid call to Reset") }
+
+// saturationDegree is a saturation degree node choice heuristic.
+type saturationDegree struct {
+	// nodes is the set of nodes being
+	// iterated over.
+	nodes []graph.Node
+
+	// indexOf is a mapping between node
+	// IDs and elements of degree and
+	// adjColors. degrees holds the
+	// degree of each node and adjColors
+	// holds the current adjacent
+	// colors of each node.
+	indexOf   map[int64]int
+	degrees   []int
+	adjColors []set.Ints
+
+	// g and colors are the graph coloring.
+	// colors is held by both the iterator
+	// and greedyColoringOf.
+	g      graph.Undirected
+	colors map[int64]int
+
+	// work is a temporary workspace.
+	work []int
+}
+
+// newSaturationDegree returns a saturationDegree heuristic based on the
+// nodes in the given node iterator and graph, using the provided coloring.
+func newSaturationDegree(it graph.Nodes, g graph.Undirected, colors map[int64]int) saturationDegree {
+	nodes := graph.NodesOf(it)
+	sd := saturationDegree{
+		nodes:     nodes,
+		indexOf:   make(map[int64]int, len(nodes)),
+		degrees:   make([]int, len(nodes)),
+		adjColors: make([]set.Ints, len(nodes)),
+		g:         g,
+		colors:    colors,
+	}
+	for i, u := range nodes {
+		sd.degrees[i] = g.From(u.ID()).Len()
+		sd.adjColors[i] = make(set.Ints)
+		sd.indexOf[u.ID()] = i
+	}
+	for uid, c := range colors {
+		to := g.From(uid)
+		for to.Next() {
+			sd.adjColors[sd.indexOf[to.Node().ID()]].Add(c)
+		}
+	}
+	return sd
+}
+
+// reset re-initializes the saturation with the provided colors.
+func (sd *saturationDegree) reset(colors map[int64]int) {
+	sd.colors = colors
+	for i := range sd.nodes {
+		sd.adjColors[i] = make(set.Ints)
+	}
+	for uid, c := range colors {
+		to := sd.g.From(uid)
+		for to.Next() {
+			sd.adjColors[sd.indexOf[to.Node().ID()]].Add(c)
+		}
+	}
+}
+
+// dsatur implements the Dsatur heuristic from Brélaz doi:10.1145/359094.359101.
+func (sd *saturationDegree) dsatur() int {
+	maxSat, maxDeg, chosen := -1, -1, -1
+	for i, u := range sd.nodes {
+		uid := u.ID()
+		if _, ok := sd.colors[uid]; ok {
+			continue
+		}
+		s := saturationDegreeOf(uid, sd.g, sd.colors)
+		d := sd.degrees[i]
+		switch {
+		// Choose a vertex with a maximal saturation degree.
+		// If there is an equality, choose any vertex of maximal
+		// degree in the uncolored subgraph.
+		case s > maxSat, s == maxSat && d > maxDeg:
+			maxSat = s
+			maxDeg = d
+			chosen = i
+		}
+	}
+	return chosen
+}
+
+// pass implements the PASS heuristic from San Segundo doi:10.1016/j.cor.2011.10.008.
+func (sd *saturationDegree) pass() int {
+	maxSat, chosen := -1, -1
+	sd.work = sd.work[:0]
+	for i, u := range sd.nodes {
+		uid := u.ID()
+		if _, ok := sd.colors[uid]; ok {
+			continue
+		}
+		s := saturationDegreeOf(uid, sd.g, sd.colors)
+		switch {
+		case s > maxSat:
+			maxSat = s
+			sd.work = sd.work[:0]
+			fallthrough
+		case s == maxSat:
+			sd.work = append(sd.work, i)
+		}
+	}
+	maxAvail := -1
+	for _, vs := range sd.work {
+		var avail int
+		for _, v := range sd.work {
+			if v != vs {
+				avail += sd.same(sd.adjColors[vs], sd.adjColors[v])
+			}
+		}
+		switch {
+		case avail > maxAvail, avail == maxAvail && sd.nodes[chosen].ID() < sd.nodes[vs].ID():
+			maxAvail = avail
+			chosen = vs
+		}
+	}
+	return chosen
+}
+
+// same implements the same function from San Segundo doi:10.1016/j.cor.2011.10.008.
+func (sd *saturationDegree) same(vi, vj set.Ints) int {
+	valid := make(set.Ints)
+	for _, c := range sd.colors {
+		if !vi.Has(c) {
+			valid.Add(c)
+		}
+	}
+	for c := range vj {
+		valid.Remove(c)
+	}
+	return valid.Count()
+}
+
+// saturationDegreeOf returns the saturation degree of the node corresponding to
+// vid in g with the given coloring.
+func saturationDegreeOf(vid int64, g graph.Undirected, colors map[int64]int) int {
+	if _, ok := colors[vid]; ok {
+		panic("coloring: saturation degree not defined for colored node")
+	}
+	adjColors := make(set.Ints)
+	to := g.From(vid)
+	for to.Next() {
+		if c, ok := colors[to.Node().ID()]; ok {
+			adjColors.Add(c)
+		}
+	}
+	return adjColors.Count()
+}
diff --git a/graph/coloring/doc.go b/graph/coloring/doc.go
new file mode 100644
index 00000000..ff42379f
--- /dev/null
+++ b/graph/coloring/doc.go
@@ -0,0 +1,6 @@
+// Copyright ©2021 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 coloring provides graph coloring functions.
+package coloring // import "gonum.org/v1/gonum/graph/coloring"


Test patch:
diff --git a/graph/coloring/coloring_test.go b/graph/coloring/coloring_test.go
new file mode 100644
index 00000000..7a53f6e7
--- /dev/null
+++ b/graph/coloring/coloring_test.go
@@ -0,0 +1,773 @@
+// Copyright ©2021 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 coloring
+
+import (
+	"context"
+	"flag"
+	"testing"
+	"time"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/encoding/graph6"
+	"gonum.org/v1/gonum/graph/internal/set"
+	"gonum.org/v1/gonum/graph/simple"
+)
+
+var runLong = flag.Bool("color.long", false, "run long exact coloring tests")
+
+var coloringTests = []struct {
+	name   string
+	g      graph.Undirected
+	colors int
+
+	long bool
+
+	partial map[int64]int
+
+	dsatur      set.Ints
+	randomized  set.Ints
+	rlf         set.Ints
+	sanSegundo  set.Ints
+	welshPowell set.Ints
+}{
+	{
+		name:   "empty",
+		g:      simple.NewUndirectedGraph(),
+		colors: 0,
+	},
+	{
+		name:   "singleton", // https://hog.grinvin.org/ViewGraphInfo.action?id=1310
+		g:      graph6.Graph("@"),
+		colors: 1,
+	},
+	{
+		name:       "kite", // https://hog.grinvin.org/ViewGraphInfo.action?id=782
+		g:          graph6.Graph("DzC"),
+		colors:     3,
+		randomized: setOf(3, 4),
+
+		partial: map[int64]int{1: 2},
+	},
+	{
+		name:   "triangle+1",
+		g:      graph6.Graph("Cw"),
+		colors: 3,
+
+		partial: map[int64]int{1: 2},
+	},
+	{
+		name:       "bipartite halves",
+		g:          graph6.Graph("G?]uf?"),
+		colors:     2,
+		randomized: setOf(2, 3, 4),
+
+		partial: map[int64]int{1: 3},
+	},
+	{
+		name:        "bipartite alternating",
+		g:           graph6.Graph("GKUalO"),
+		colors:      2,
+		randomized:  setOf(2, 3, 4),
+		welshPowell: setOf(2, 3, 4),
+
+		partial: map[int64]int{1: 1},
+	},
+	{
+		name:   "3/4 bipartite", // https://hog.grinvin.org/ViewGraphInfo.action?id=466
+		g:      graph6.Graph("F?~v_"),
+		colors: 2,
+
+		partial: map[int64]int{1: 1},
+	},
+	{
+		name:        "cubical", // https://hog.grinvin.org/ViewGraphInfo.action?id=1022
+		g:           graph6.Graph("Gs@ipo"),
+		colors:      2,
+		randomized:  setOf(2, 3, 4),
+		welshPowell: setOf(2, 3),
+
+		partial: map[int64]int{1: 1},
+	},
+	{
+		name:       "HoG33106", // https://hog.grinvin.org/ViewGraphInfo.action?id=33106
+		g:          graph6.Graph("K???WWKKxf]C"),
+		colors:     2,
+		randomized: setOf(2, 3, 4),
+
+		partial: map[int64]int{1: 1},
+	},
+	{
+		name:        "HoG41237", // https://hog.grinvin.org/ViewGraphInfo.action?id=41237
+		g:           graph6.Graph("~?BCs?GO?@?O??@?A????G?A?????G??O????????????????????G????_???O??????????????????????OC????A???????????????@@?????A?????????????????????????C??????G??????O?????????????????????????G?@????????@???????D????????_??????CG???????@???????????????????????ACG???????EA???????O?G??????S????????C???_?????O??C?????@??????????W????????B??????????????O?????G??????????_????o????A????A_???@O_????????@?O????G???W???????????@G?????????A??????@?????S??????????????????????????????B????????????K????????????????????????@_????????????o??G????????????A?????????????_?????????_??G?????????G????@_???????????C?_????????????@????????A???A?????????A???O?????????????C???????????K??_??????????@O??_??????????C??@???????????G??????????E????????A??G??????????????????C????A???@???????????O??O???????????????_????????????@_?_????????????@O?C?????????????_?@?????????????G????C???O???????????????G?O????????????????O???????O???C????????????O?A?????????????????_??????????????@_C???????????????I?C??????????????@??G???????????????_???CA??????????????????_???????O????????O???????????????A??????????????????G?????????????????_C?????????????????OC?????????????????g????????D?????????????????@?G?????????????????GA??????????????????K??????????????????@O???????????????G????????G????????????????????P???????????????????OO???????????????????@_??????????????????????????????????????????????????????????????????@?????????????????????_???A?????????????????????????????????????????????????????????????O?????????????????????_??@????????????????????????????????????????G???????????@????????????????????????????????C?????????????????????????????????????????????????K??????????????????????I??????????????????????C_???????????????????_C????????????????????????K???????????????????????K???????????????????????`???????????????????????`??????????????????????I???????????????????????C?O?????????????????????A?O??????????????????????O?G??????????????????????@?C???????????????????????H????????????????????????OG???????????????????????GC?????????????????????????B?????????????????????????I?????????????????????????????????????????????????????????????????????????????????????????????????????????????c?????????????????????????CO???????????????????????A?A????????????????????????C?@????????????????????????C??_???????????????????????O??O?????????????????????????o???????????????????????????K???????????????????????????E??????????????????????????@_?????????????????????????D????????????????????????????`???????????????????????????????CG??????????????????????????@?O??????????????????????????@O???????????????????????????CG???????????????????????????C@???????????????????????????A?_?????????????????????C?????????????????????????????G?????????????????????????????@????????????????????????????????????AA?????????????????????????????B?????????????????????????????AG?????????????????????????????G_????????????????????????????CA?????????????????????????????CG??????????????????????????A???????????????????????????????O??????????????????????????????G????A?????????????????????????????AAG?????????????????????????????CCC??????????????????????????????cG?????????????????????????????@CA?????????????????????????????AG@??????????????????????????????P?"),
+		colors:      2,
+		randomized:  setOf(2, 3, 4),
+		welshPowell: setOf(2, 3),
+
+		partial: map[int64]int{1: 1},
+	},
+
+	// Test graphs from Leighton doi:10.6028/jres.084.024.
+	{
+		name:       "exams", // Figure 4 induced by the nodes in E.
+		g:          graph6.Graph("EDZw"),
+		colors:     3,
+		randomized: setOf(3, 4),
+
+		partial: map[int64]int{1: 2},
+	},
+	{
+		name:       "exam_scheduler_1", // Figure 4.
+		g:          graph6.Graph("JQ?HaWN{l~?"),
+		colors:     5,
+		randomized: setOf(5, 6),
+	},
+	{
+		name:       "exam_scheduler_2", // Figure 11.
+		g:          graph6.Graph("GTPHz{"),
+		colors:     4,
+		randomized: setOf(4, 5),
+
+		partial: map[int64]int{1: 3},
+	},
+
+	// Test graph from Brélaz doi:10.1145/359094.359101.
+	{
+		name:        "Brélaz", // Figure 1.
+		g:           graph6.Graph(`I??GG\rmg`),
+		colors:      3,
+		randomized:  setOf(3, 4),
+		rlf:         setOf(3, 4),
+		welshPowell: setOf(3, 4),
+
+		partial: map[int64]int{1: 2},
+	},
+
+	// Test graph from Lima and Carmo doi:10.22456/2175-2745.80721.
+	{
+		name:       "Lima Carmo", // Figure 2.
+		g:          graph6.Graph("Gh]S?G"),
+		colors:     3,
+		randomized: setOf(3, 4),
+		sanSegundo: setOf(3, 4),
+
+		partial: map[int64]int{1: 2},
+	},
+
+	// Test graph from San Segundo doi:10.1016/j.cor.2011.10.008.
+	{
+		name:       "San Segundo", // Figure 1 A.
+		g:          graph6.Graph(`HMn\r\v`),
+		colors:     5,
+		randomized: setOf(5, 6),
+
+		partial: map[int64]int{1: 4},
+	},
+
+	{
+		name:   "tetrahedron", // https://hog.grinvin.org/ViewGraphInfo.action?id=74
+		g:      graph6.Graph("C~"),
+		colors: 4,
+
+		partial: map[int64]int{1: 3},
+	},
+
+	{
+		name:   "triangle", // https://hog.grinvin.org/ViewGraphInfo.action?id=1374
+		g:      graph6.Graph("Bw"),
+		colors: 3,
+
+		partial: map[int64]int{1: 2},
+	},
+	{
+		name:   "square", // https://hog.grinvin.org/ViewGraphInfo.action?id=674
+		g:      graph6.Graph("Cl"),
+		colors: 2,
+
+		partial: map[int64]int{1: 1},
+	},
+	{
+		name:   "cycle-5", // https://hog.grinvin.org/ViewGraphInfo.action?id=340
+		g:      graph6.Graph("Dhc"),
+		colors: 3,
+
+		partial: map[int64]int{1: 2},
+	},
+	{
+		name:       "cycle-6", // https://hog.grinvin.org/ViewGraphInfo.action?id=670
+		g:          graph6.Graph("EhEG"),
+		colors:     2,
+		randomized: setOf(2, 3),
+
+		partial: map[int64]int{1: 1},
+	},
+
+	{
+		name:   "wheel-5", // https://hog.grinvin.org/ViewGraphInfo.action?id=442
+		g:      graph6.Graph("D|s"),
+		colors: 3,
+
+		partial: map[int64]int{1: 2},
+	},
+	{
+		name:   "wheel-6", // https://hog.grinvin.org/ViewGraphInfo.action?id=204
+		g:      graph6.Graph("E|fG"),
+		colors: 4,
+
+		partial: map[int64]int{1: 3},
+	},
+	{
+		name:       "wheel-7",
+		g:          graph6.Graph("F|eMG"),
+		colors:     3,
+		randomized: setOf(3, 4),
+
+		partial: map[int64]int{1: 2},
+	},
+
+	{
+		name:        "sudoku board", // The graph of the constraints of a sudoku puzzle board.
+		g:           graph6.Graph("~?@P~~~~~~wF?{BFFFFbbwF~?~{B~~?wF?wF_[BF?wwwFFb_[^?wF~?wF~_[B~{?_C?OA?OC_C?_W_C?_wOA?O{C?_C^?_C?fwA?OA~_C?_F~_C?_F?OA?OF?c?_CB_W_C?_FFA?OA?w{C?_CBbwC?_C?~wA?OA?~{?_C?_^~_C?_F?wA?OA?wF?c?_CB_[BC?_C?wFFA?OA?wFF__C?_[BbwC?_C?wF~?OA?OF?~{?_C?_[B~{?_C?_C?_A?OA?OA?OC_C?_C?_CBC?_C?_C?_wOA?OA?OAF__C?_C?_C^?_C?_C?_C~?OA?OA?OA~_C?_C?_C?~{?_C?_C?_F?OA?OA?OA?wC_C?_C?_CB_W_C?_C?_C?wwOA?OA?OA?w{C?_C?_C?_[^?_C?_C?_C?~wA?OA?OA?OF~_C?_C?_C?_^~_C?_C?_C?wF?OA?OA?OA?wF?c?_C?_C?_[B_W_C?_C?_C?wFFA?OA?OA?OF?w{C?_C?_C?_[BbwC?_C?_C?_F?~wA?OA?OA?OF?~{?_C?_C?_CB_^~"),
+		colors:      9,
+		randomized:  setOf(9, 10, 11, 12, 13),
+		sanSegundo:  setOf(9, 10),
+		welshPowell: setOf(9, 10, 11, 12, 13, 14),
+
+		partial: map[int64]int{1: 3},
+	},
+	{
+		name:        "sudoku problem", // The constraint graph for the problem in the sudoku example.
+		g:           graph6.Graph("~?@Y~~~~~~?F|_B?F?F_Bw?~?F{?^w?w??wC?[Bzwww?FF_?[^??F~??F~~kB~w?wF?v~?wFz{B_W?F?w|~F?w{?B_[^??F?~w??wF~_?B_^~?C?_C??A?OA?_?_C?_W?C?_CF??OA?O~|__C?bw??_C?fw??OA?V~}_C?_F~?C?_C?wD~OA?OF?_?_C?_[B??_C?_FF}wOA?OFF_?C?_CBbw??_C?_F~~oA?OA?~{??C?_CB~~^_C?_F?w??OA?OF?wC?C?_CB_[B|w_C?_F?ww?A?OA?wFF_?C?_CB_[^??C?_C?wF~??A?OA?wF~_??_C?_[B~w?_C?_C?_C??A?OA?OA?OC?C?_C?_C?_W?C?_C?_C?_~}A?OA?OA?O{??_C?_C?_C^vwC?_C?_C?f|~?OA?OA?OA~_??_C?_C?_F~|{?_C?_C?_F??A?OA?OA?OF?_?_C?_C?_CB_W?C?_C?_C?_FF??OA?OA?OA?w|~__C?_C?_CBbw??_C?_C?_C?~w??OA?OA?OA?~{??C?_C?_C?_^~?C?_C?_C?_F?w??OA?OA?OA?wF?_?_C?_C?_CB_[B??_C?_C?_C?wFF??OA?OA?OA?wFF_?C?_C?_C?_[Bbw??_C?_C?_C?wF~}wA?OA?OA?OF?~{??C?_C?_C?_[B~w"),
+		colors:      9,
+		dsatur:      setOf(9, 10, 11),
+		randomized:  setOf(9, 10, 11, 12, 13, 14, 15, 16, 17),
+		rlf:         setOf(9, 10, 11, 12),
+		sanSegundo:  setOf(9, 10, 11),
+		welshPowell: setOf(9, 10, 11, 12, 13),
+
+		partial: map[int64]int{1: 3},
+	},
+
+	// Test graphs from NetworkX.
+	{
+		name:       "cs_shc",
+		g:          graph6.Graph("Djs"),
+		colors:     3,
+		randomized: setOf(3, 4),
+	},
+	{
+		name:       "gis_hc",
+		g:          graph6.Graph("E?ow"),
+		colors:     2,
+		randomized: setOf(2, 3),
+	},
+	{
+		name:       "gis_shc|rs_shc", // https://hog.grinvin.org/ViewGraphInfo.action?id=594
+		g:          graph6.Graph("CR"),
+		colors:     2,
+		randomized: setOf(2, 3),
+	},
+	{
+		name:        "lf_hc",
+		g:           graph6.Graph(`F\^E?`),
+		colors:      3,
+		dsatur:      setOf(3, 4),
+		randomized:  setOf(3, 4),
+		rlf:         setOf(3, 4),
+		sanSegundo:  setOf(3, 4),
+		welshPowell: setOf(3, 4),
+	},
+	{
+		name:        "lf_shc",
+		g:           graph6.Graph("ELQ?"),
+		colors:      2,
+		randomized:  setOf(2, 3),
+		welshPowell: setOf(2, 3),
+	},
+	{
+		name:        "lfi_hc",
+		g:           graph6.Graph("Hhe[b@_"),
+		colors:      3,
+		dsatur:      setOf(3, 4),
+		randomized:  setOf(3, 4),
+		sanSegundo:  setOf(3, 4),
+		welshPowell: setOf(3, 4),
+	},
+	{
+		name:        "lfi_shc|slf_shc",
+		g:           graph6.Graph("FheZ?"),
+		colors:      3,
+		dsatur:      setOf(3, 4),
+		randomized:  setOf(3, 4),
+		welshPowell: setOf(3, 4),
+	},
+	{
+		name:       "no_solo", // https://hog.grinvin.org/ViewGraphInfo.action?id=264, https://hog.grinvin.org/ViewGraphInfo.action?id=498
+		g:          graph6.Graph("K????AccaQHG"),
+		colors:     2,
+		randomized: setOf(2, 3),
+	},
+	{
+		name:        "sl_hc",
+		g:           graph6.Graph("Gzg[Yk"),
+		colors:      4,
+		dsatur:      setOf(4, 5),
+		randomized:  setOf(4, 5),
+		sanSegundo:  setOf(4, 5),
+		welshPowell: setOf(4, 5),
+	},
+	{
+		name:        "sl_shc",
+		g:           graph6.Graph("E{Sw"),
+		colors:      3,
+		dsatur:      setOf(3, 4),
+		randomized:  setOf(3, 4),
+		welshPowell: setOf(3, 4),
+	},
+	{
+		name:        "slf_hc",
+		g:           graph6.Graph("G}`?W["),
+		colors:      3,
+		dsatur:      setOf(3, 4),
+		randomized:  setOf(3, 4),
+		rlf:         setOf(3, 4),
+		welshPowell: setOf(3, 4),
+	},
+	{
+		name:        "sli_hc",
+		g:           graph6.Graph("H{czYtt"),
+		colors:      4,
+		dsatur:      setOf(4, 5),
+		randomized:  setOf(4, 5),
+		welshPowell: setOf(4, 5),
+	},
+	{
+		name:        "sli_shc",
+		g:           graph6.Graph("FxdSW"),
+		colors:      3,
+		dsatur:      setOf(3, 4),
+		randomized:  setOf(3, 4),
+		sanSegundo:  setOf(3, 4),
+		welshPowell: setOf(3, 4),
+	},
+	{
+		name:       "rsi_shc",
+		g:          graph6.Graph("EheW"),
+		colors:     3,
+		randomized: setOf(3, 4),
+	},
+	{
+		name:       "V_plus_not_in_A_cal",
+		g:          graph6.Graph("HQQ?W__"),
+		colors:     2,
+		randomized: setOf(2, 3),
+	},
+
+	// DIMACS queens graphs
+	{
+		name:        "queen5_5",
+		g:           graph6.Graph("X~~FJk~F|KIxizS^dF{iWQjcdV[dFyQb}KiWOdVHAT\\acg~acg~"),
+		colors:      5,
+		dsatur:      setOf(5),
+		randomized:  setOf(5, 6, 7, 8, 9, 10),
+		rlf:         setOf(5),
+		sanSegundo:  setOf(5),
+		welshPowell: setOf(5, 6, 7),
+	},
+	{
+		name:        "queen6_6",
+		g:           graph6.Graph("c~~}FDrMw~`~goSwtMYhvIF{SN{dEAQfCehrcTMyPO~ca`~acgoPQSwcCtMWcahvaQIF|CcSN{KSdEAAIQfC__ehrCCcTMwSQPO~ogca`~"),
+		colors:      7,
+		dsatur:      setOf(7, 8, 9),
+		randomized:  setOf(7, 8, 9, 10, 11, 12),
+		rlf:         setOf(7, 8),
+		sanSegundo:  setOf(7, 8),
+		welshPowell: setOf(7, 8, 9),
+	},
+	{
+		name:        "queen7_7",
+		g:           graph6.Graph("p~~~}B`[XrbnB~@~sKAb`iMLS[xS\\wgN{IB~cSKAPPocibbcibfQDPvc`O^wcIB~aQIE@CcS[HCdS[WaQiM]GcIbnPC`O^xCQD@~ogca`_Ogcab`GCQTPpa@CdS[wQGcIbn`GaOgN|APC`O^{EDCcSKA@AaQIMC_OGcibbCA@CdS[wOHCQDPv_gQGcIB~_gQGcIB~"),
+		colors:      7,
+		dsatur:      setOf(7, 8, 9, 10, 11),
+		randomized:  setOf(7, 8, 9, 10, 11, 12, 13),
+		rlf:         setOf(7, 8, 9),
+		sanSegundo:  setOf(7, 8, 9),
+		welshPowell: setOf(7, 8, 9, 10, 11, 12),
+	},
+	{
+		name:        "queen8_8",
+		long:        true,
+		g:           graph6.Graph("~?@?~~~~~?wJ`fFFNBn_~wF~gK@OwLOwYg[[iFNDOzwSB~_gF~cIB?QDB_cdOw[ciFFQQg[{cDOzwcD?~wQA_^}GcIB?PC`OwHCQTB`aKciFFaCciFNPAOTBncOcD?~waC_gF~`GaOgK@APC`OwHAPCdOwW_GqQg[[GaCciFN`COcDOzyCPAOSB~cGaC_gF~_gQGcIB?OSHCQDB_c@APCdOwW_GAKciFFA?aGQQg[{C`COcDOz{C`COcD?~yAOaGQA_^}@_gQGcIB?OCDAPC`OwH?OCHCQTB`a?_GAKciFFA?_GaCciFN@?QCPAOTBn_SC`COcD?~{A_cGaC_gF~"),
+		colors:      9,
+		dsatur:      setOf(9, 10, 11, 12),
+		randomized:  setOf(9, 10, 11, 12, 13, 14, 15),
+		rlf:         setOf(9, 10),
+		sanSegundo:  setOf(9, 10, 11),
+		welshPowell: setOf(9, 10, 11, 12, 13, 14, 15),
+	},
+}
+
+func setOf(vals ...int) set.Ints {
+	s := make(set.Ints)
+	for _, v := range vals {
+		s.Add(v)
+	}
+	return s
+}
+
+func TestDsatur(t *testing.T) {
+	for _, test := range coloringTests {
+		for _, partial := range []map[int64]int{nil, test.partial} {
+			k, colors, err := Dsatur(test.g, partial)
+
+			if partial == nil && k != test.colors && !test.dsatur.Has(k) {
+				t.Errorf("unexpected chromatic number for %q: got:%d want:%d or in %v\ncolors:%v",
+					test.name, k, test.colors, test.dsatur, colors)
+			}
+			if s := Sets(colors); len(s) != k {
+				t.Errorf("mismatch between number of color sets and k: |sets|=%d k=%d", len(s), k)
+			}
+			if missing, ok := isCompleteColoring(colors, test.g); !ok {
+				t.Errorf("incomplete coloring for %q: missing %d\ngot:%v", test.name, missing, colors)
+			}
+			if xid, yid, ok := isValidColoring(colors, test.g); !ok {
+				t.Errorf("invalid coloring for %q: %d--%d match color\ncolors:%v",
+					test.name, xid, yid, colors)
+			}
+			if err != nil {
+				t.Errorf("unexpected error: %v", err)
+			}
+
+			for id, c := range partial {
+				if colors[id] != c {
+					t.Errorf("coloring not consistent with input partial for %q:\ngot:%v\nwant superset of:%v",
+						test.name, colors, partial)
+					break
+				}
+			}
+		}
+	}
+}
+
+func TestDsaturExact(t *testing.T) {
+	timeout := time.Microsecond
+	for _, test := range coloringTests {
+		for _, useTimeout := range []bool{false, true} {
+			if test.long && !*runLong && !useTimeout {
+				continue
+			}
+			var term Terminator
+			cancel := func() {}
+			if useTimeout {
+				term, cancel = context.WithTimeout(context.Background(), timeout)
+			}
+			k, colors, err := DsaturExact(term, test.g)
+			cancel()
+			if k != test.colors && (useTimeout && !test.dsatur.Has(k)) {
+				t.Errorf("unexpected chromatic number for %q timeout=%t: got:%d want:%d or in %v\ncolors:%v",
+					test.name, useTimeout, k, test.colors, test.dsatur, colors)
+			}
+			if s := Sets(colors); len(s) != k {
+				t.Errorf("mismatch between number of color sets and k: |sets|=%d k=%d", len(s), k)
+			}
+			if missing, ok := isCompleteColoring(colors, test.g); !ok {
+				t.Errorf("incomplete coloring for %q: missing %d\ngot:%v", test.name, missing, colors)
+			}
+			if xid, yid, ok := isValidColoring(colors, test.g); !ok {
+				t.Errorf("invalid coloring for %q: %d--%d match color\ncolors:%v",
+					test.name, xid, yid, colors)
+			}
+			if err != nil && !useTimeout {
+				t.Errorf("unexpected error: %v", err)
+			}
+		}
+	}
+}
+
+func TestRandomized(t *testing.T) {
+	for seed := uint64(1); seed <= 1000; seed++ {
+		for _, test := range coloringTests {
+			for _, partial := range []map[int64]int{nil, test.partial} {
+				k, colors, err := Randomized(test.g, partial, rand.NewSource(seed))
+
+				if partial == nil && k != test.colors && !test.randomized.Has(k) {
+					t.Errorf("unexpected chromatic number for %q with seed=%d: got:%d want:%d or in %v\ncolors:%v",
+						test.name, seed, k, test.colors, test.randomized, colors)
+				}
+				if s := Sets(colors); len(s) != k {
+					t.Errorf("mismatch between number of color sets and k: |sets|=%d k=%d", len(s), k)
+				}
+				if missing, ok := isCompleteColoring(colors, test.g); !ok {
+					t.Errorf("incomplete coloring for %q: missing %d\ngot:%v", test.name, missing, colors)
+				}
+				if xid, yid, ok := isValidColoring(colors, test.g); !ok {
+					t.Errorf("invalid coloring for %q: %d--%d match color\ncolors:%v",
+						test.name, xid, yid, colors)
+				}
+				if err != nil {
+					t.Errorf("unexpected error: %v", err)
+				}
+
+				for id, c := range partial {
+					if colors[id] != c {
+						t.Errorf("coloring not consistent with input partial for %q:\ngot:%v\nwant superset of:%v",
+							test.name, colors, partial)
+						break
+					}
+				}
+			}
+		}
+	}
+}
+
+func TestRecursiveLargestFirst(t *testing.T) {
+	for _, test := range coloringTests {
+		k, colors := RecursiveLargestFirst(test.g)
+		if k != test.colors && !test.rlf.Has(k) {
+			t.Errorf("unexpected chromatic number for %q: got:%d want:%d",
+				test.name, k, test.colors)
+		}
+		if s := Sets(colors); len(s) != k {
+			t.Errorf("mismatch between number of color sets and k: |sets|=%d k=%d", len(s), k)
+		}
+		if missing, ok := isCompleteColoring(colors, test.g); !ok {
+			t.Errorf("incomplete coloring for %q: missing %d\ngot:%v", test.name, missing, colors)
+		}
+		if xid, yid, ok := isValidColoring(colors, test.g); !ok {
+			t.Errorf("invalid coloring for %q: %d--%d match color\ncolors:%v",
+				test.name, xid, yid, colors)
+		}
+	}
+}
+
+func TestSanSegundo(t *testing.T) {
+	for _, test := range coloringTests {
+		for _, partial := range []map[int64]int{nil, test.partial} {
+			k, colors, err := SanSegundo(test.g, partial)
+
+			if partial == nil && k != test.colors && !test.sanSegundo.Has(k) {
+				t.Errorf("unexpected chromatic number for %q: got:%d want:%d\ncolors:%v",
+					test.name, k, test.colors, colors)
+			}
+			if s := Sets(colors); len(s) != k {
+				t.Errorf("mismatch between number of color sets and k: |sets|=%d k=%d", len(s), k)
+			}
+			if missing, ok := isCompleteColoring(colors, test.g); !ok {
+				t.Errorf("incomplete coloring for %q: missing %d\ngot:%v", test.name, missing, colors)
+			}
+			if xid, yid, ok := isValidColoring(colors, test.g); !ok {
+				t.Errorf("invalid coloring for %q: %d--%d match color\ncolors:%v",
+					test.name, xid, yid, colors)
+			}
+			if err != nil {
+				t.Errorf("unexpected error: %v", err)
+			}
+
+			for id, c := range partial {
+				if colors[id] != c {
+					t.Errorf("coloring not consistent with input partial for %q:\ngot:%v\nwant superset of:%v",
+						test.name, colors, partial)
+					break
+				}
+			}
+		}
+	}
+}
+
+func TestWelshPowell(t *testing.T) {
+	for _, test := range coloringTests {
+		for _, partial := range []map[int64]int{nil, test.partial} {
+			k, colors, err := WelshPowell(test.g, partial)
+
+			if partial == nil && k != test.colors && !test.welshPowell.Has(k) {
+				t.Errorf("unexpected chromatic number for %q: got:%d want:%d",
+					test.name, k, test.colors)
+			}
+			if s := Sets(colors); len(s) != k {
+				t.Errorf("mismatch between number of color sets and k: |sets|=%d k=%d", len(s), k)
+			}
+			if missing, ok := isCompleteColoring(colors, test.g); !ok {
+				t.Errorf("incomplete coloring for %q: missing %d\ngot:%v", test.name, missing, colors)
+			}
+			if xid, yid, ok := isValidColoring(colors, test.g); !ok {
+				t.Errorf("invalid coloring for %q: %d--%d match color\ncolors:%v",
+					test.name, xid, yid, colors)
+			}
+			if err != nil {
+				t.Errorf("unexpected error: %v", err)
+			}
+
+			for id, c := range partial {
+				if colors[id] != c {
+					t.Errorf("coloring not consistent with input partial for %q:\ngot:%v\nwant superset of:%v",
+						test.name, colors, partial)
+					break
+				}
+			}
+		}
+	}
+}
+
+var newPartialTests = []struct {
+	partial   map[int64]int
+	g         graph.Undirected
+	wantValid bool
+}{
+	{
+		partial: map[int64]int{0: 1, 1: 1},
+		g: undirectedGraphFrom([]intset{
+			0: linksTo(1, 2),
+			1: linksTo(2),
+			2: nil,
+		}),
+		wantValid: false,
+	},
+	{
+		partial: map[int64]int{0: 0, 1: 1},
+		g: undirectedGraphFrom([]intset{
+			0: linksTo(1, 2),
+			1: linksTo(2),
+			2: nil,
+		}),
+		wantValid: true,
+	},
+	{
+		partial: map[int64]int{0: 0, 1: 1, 3: 3},
+		g: undirectedGraphFrom([]intset{
+			0: linksTo(1, 2),
+			1: linksTo(2),
+			2: nil,
+		}),
+		wantValid: false,
+	},
+	{
+		partial: map[int64]int{0: 0, 1: 1, 2: 2},
+		g: undirectedGraphFrom([]intset{
+			0: linksTo(1, 2),
+			1: linksTo(2),
+			2: nil,
+		}),
+		wantValid: true,
+	},
+	{
+		partial: nil,
+		g: undirectedGraphFrom([]intset{
+			0: linksTo(1, 2),
+			1: linksTo(2),
+			2: nil,
+		}),
+		wantValid: true,
+	},
+}
+
+func TestNewPartial(t *testing.T) {
+	for i, test := range newPartialTests {
+		gotPartial, gotValid := newPartial(test.partial, test.g)
+		if gotValid != test.wantValid {
+			t.Errorf("unexpected validity for test %d: got:%t want:%t",
+				i, gotValid, test.wantValid)
+		}
+		xid, yid, ok := isValidColoring(gotPartial, test.g)
+		if !ok {
+			t.Errorf("invalid partial returned for test %d: %d--%d match color\ncolors:%v",
+				i, xid, yid, gotPartial)
+		}
+
+	}
+}
+
+func isCompleteColoring(colors map[int64]int, g graph.Undirected) (missing int64, ok bool) {
+	for _, n := range graph.NodesOf(g.Nodes()) {
+		if _, ok := colors[n.ID()]; !ok {
+			return n.ID(), false
+		}
+	}
+	return 0, true
+}
+
+func isValidColoring(colors map[int64]int, g graph.Undirected) (x, y int64, ok bool) {
+	for xid, c := range colors {
+		to := g.From(xid)
+		for to.Next() {
+			yid := to.Node().ID()
+			if oc, ok := colors[yid]; ok && c == oc {
+				return xid, yid, false
+			}
+		}
+	}
+	return 0, 0, true
+}
+
+// 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
+}
+
+func undirectedGraphFrom(g []intset) graph.Undirected {
+	dg := simple.NewUndirectedGraph()
+	for u, e := range g {
+		for v := range e {
+			dg.SetEdge(simple.Edge{F: simple.Node(u), T: simple.Node(v)})
+		}
+	}
+	return dg
+}
+
+func BenchmarkColoring(b *testing.B) {
+	for _, bench := range coloringTests {
+		b.Run(bench.name, func(b *testing.B) {
+			b.Run("Dsatur", func(b *testing.B) {
+				for i := 0; i < b.N; i++ {
+					_, _, err := Dsatur(bench.g, nil)
+					if err != nil {
+						b.Fatalf("coloring failed: %v", err)
+					}
+				}
+			})
+			if !bench.long || *runLong {
+				b.Run("DsaturExact", func(b *testing.B) {
+					for i := 0; i < b.N; i++ {
+						_, _, err := DsaturExact(nil, bench.g)
+						if err != nil {
+							b.Fatalf("coloring failed: %v", err)
+						}
+					}
+				})
+			}
+			b.Run("RecursiveLargestFirst", func(b *testing.B) {
+				for i := 0; i < b.N; i++ {
+					RecursiveLargestFirst(bench.g)
+				}
+			})
+			b.Run("SanSegundo", func(b *testing.B) {
+				for i := 0; i < b.N; i++ {
+					_, _, err := SanSegundo(bench.g, nil)
+					if err != nil {
+						b.Fatalf("coloring failed: %v", err)
+					}
+				}
+			})
+			b.Run("WelshPowell", func(b *testing.B) {
+				for i := 0; i < b.N; i++ {
+					_, _, err := WelshPowell(bench.g, nil)
+					if err != nil {
+						b.Fatalf("coloring failed: %v", err)
+					}
+				}
+			})
+		})
+	}
+}
diff --git a/graph/coloring/sudoku_example_test.go b/graph/coloring/sudoku_example_test.go
new file mode 100644
index 00000000..a9e96690
--- /dev/null
+++ b/graph/coloring/sudoku_example_test.go
@@ -0,0 +1,134 @@
+// Copyright ©2021 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 coloring_test
+
+import (
+	"fmt"
+	"log"
+
+	"gonum.org/v1/gonum/graph/coloring"
+	"gonum.org/v1/gonum/graph/graphs/gen"
+	"gonum.org/v1/gonum/graph/simple"
+)
+
+// A hard sudoku problem graded at a level of difficulty, "not fun".
+// https://dingo.sbs.arizona.edu/~sandiway/sudoku/examples.html
+var grid = [9][9]int{
+	{0, 2, 0 /**/, 0, 0, 0 /**/, 0, 0, 0},
+	{0, 0, 0 /**/, 6, 0, 0 /**/, 0, 0, 3},
+	{0, 7, 4 /**/, 0, 8, 0 /**/, 0, 0, 0},
+
+	{0, 0, 0 /**/, 0, 0, 3 /**/, 0, 0, 2},
+	{0, 8, 0 /**/, 0, 4, 0 /**/, 0, 1, 0},
+	{6, 0, 0 /**/, 5, 0, 0 /**/, 0, 0, 0},
+
+	{0, 0, 0 /**/, 0, 1, 0 /**/, 7, 8, 0},
+	{5, 0, 0 /**/, 0, 0, 9 /**/, 0, 0, 0},
+	{0, 0, 0 /**/, 0, 0, 0 /**/, 0, 4, 0},
+}
+
+func Example_sudoku() {
+	g := simple.NewUndirectedGraph()
+
+	// Build the sudoku board constraints.
+	for i := 0; i < 9; i++ {
+		gen.Complete(g, row(i))
+		gen.Complete(g, col(i))
+	}
+	for r := 0; r < 3; r++ {
+		for c := 0; c < 3; c++ {
+			gen.Complete(g, block{r, c})
+		}
+	}
+
+	// Add constraints for the digits.
+	gen.Complete(g, gen.IDRange{First: -9, Last: -1})
+
+	// Mark constraints onto the graph.
+	for r, row := range &grid {
+		for c, val := range &row {
+			if val == 0 {
+				continue
+			}
+			for i := 1; i <= 9; i++ {
+				if i != val {
+					g.SetEdge(simple.Edge{F: simple.Node(-i), T: simple.Node(id(r, c))})
+				}
+			}
+		}
+	}
+
+	k, colors, err := coloring.DsaturExact(nil, g)
+	if err != nil {
+		log.Fatal(err)
+	}
+	if k != 9 {
+		log.Fatalln("could not solve problem", k)
+	}
+	sets := coloring.Sets(colors)
+	for r := 0; r < 9; r++ {
+		if r != 0 && r%3 == 0 {
+			fmt.Println()
+		}
+		for c := 0; c < 9; c++ {
+			if c != 0 {
+				fmt.Print(" ")
+				if c%3 == 0 {
+					fmt.Print(" ")
+				}
+			}
+			got := -int(sets[colors[id(r, c)]][0])
+			if want := grid[r][c]; want != 0 && got != want {
+				log.Fatalf("mismatch at row=%d col=%d: %d != %d", r, c, got, want)
+			}
+			fmt.Print(got)
+		}
+		fmt.Println()
+	}
+
+	// Output:
+	//
+	// 1 2 6  4 3 7  9 5 8
+	// 8 9 5  6 2 1  4 7 3
+	// 3 7 4  9 8 5  1 2 6
+	//
+	// 4 5 7  1 9 3  8 6 2
+	// 9 8 3  2 4 6  5 1 7
+	// 6 1 2  5 7 8  3 9 4
+	//
+	// 2 6 9  3 1 4  7 8 5
+	// 5 4 8  7 6 9  2 3 1
+	// 7 3 1  8 5 2  6 4 9
+}
+
+// row is a gen.IDer that enumerates the IDs of graph
+// nodes representing a row of cells of a sudoku board.
+type row int
+
+func (r row) Len() int       { return 9 }
+func (r row) ID(i int) int64 { return id(int(r), i) }
+
+// col is a gen.IDer that enumerates the IDs of graph
+// nodes representing a column of cells of a sudoku board.
+type col int
+
+func (c col) Len() int       { return 9 }
+func (c col) ID(i int) int64 { return id(i, int(c)) }
+
+// block is a gen.IDer that enumerates the IDs of graph
+// nodes representing a 3×3 block of cells of a sudoku board.
+type block struct {
+	r, c int
+}
+
+func (b block) Len() int { return 9 }
+func (b block) ID(i int) int64 {
+	return id(b.r*3, b.c*3) + int64(i%3) + int64(i/3)*9
+}
+
+// id returns the graph node ID of a cell in a sudoku board.
+func id(row, col int) int64 {
+	return int64(row*9 + col)
+}