Repository: gonum/gonum
Pull request #1727 (2021-09-27 21:04:49)

README: fix line link syntax
Please take a look.

<!--
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 _____."
-->



Patch:
diff --git a/README.md b/README.md
index f7ab18ea..2808820d 100644
--- a/README.md
+++ b/README.md
@@ -18,7 +18,7 @@ go get -u gonum.org/v1/gonum/...
 
 ## Supported Go versions
 
-Gonum supports and tests using the gc compiler on the [two most recent Go releases](https://github.com/gonum/gonum/blob/master/.github/workflows/ci.yml#L14-15) on Linux (386, amd64 and arm64), macOS and Windows (both on amd64).
+Gonum supports and tests using the gc compiler on the [two most recent Go releases](https://github.com/gonum/gonum/blob/master/.github/workflows/ci.yml#L14-L15) on Linux (386, amd64 and arm64), macOS and Windows (both on amd64).
 
 ## Release schedule
 
diff --git a/graph/formats/rdf/equi_canonical.go b/graph/formats/rdf/equi_canonical.go
new file mode 100644
index 00000000..9d16033e
--- /dev/null
+++ b/graph/formats/rdf/equi_canonical.go
@@ -0,0 +1,576 @@
+// 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 rdf
+
+import (
+	"errors"
+	"sort"
+)
+
+// Throughout, the comments refer to doi:10.1145/3068333 which should be
+// understood as a synonym for http://aidanhogan.com/docs/rdf-canonicalisation.pdf
+// although there are differences between the two, see http://aidanhogan.com/#errataH17.
+// Where there are differences, the document at http://aidanhogan.com/ is the
+// canonical truth. The DOI reference is referred to for persistence.
+
+// Lean returns an RDF core of g that entails g. If g contains any non-zero
+// labels, Lean will return a non-nil error and a core of g assuming no graph
+// labels exist.
+//
+// See http://aidanhogan.com/docs/rdf-canonicalisation.pdf for details of
+// the algorithm.
+func Lean(g []*Statement) ([]*Statement, error) {
+	// BUG(kortschak): Graph leaning does not take into account graph label terms
+	// since the formal semantics for a multiple graph data model have not been
+	// defined. See https://www.w3.org/TR/rdf11-datasets/#declaring.
+
+	var (
+		hasBlanks bool
+		err       error
+	)
+	for _, s := range g {
+		if isBlank(s.Subject.Value) || isBlank(s.Object.Value) {
+			hasBlanks = true
+			if err != nil {
+				break
+			}
+		}
+		if s.Label.Value != "" && err == nil {
+			err = errors.New("rdf: data-set contains graph names")
+			if hasBlanks {
+				break
+			}
+		}
+	}
+	if hasBlanks {
+		g = lean(&dfs{}, g)
+	}
+	return g, err
+}
+
+// removeRedundantBnodes removes blank nodes whose edges are a subset of
+// another term in the RDF graph.
+//
+// This is algorithm 4 in doi:10.1145/3068333.
+func removeRedundantBnodes(g []*Statement) []*Statement {
+	g = append(g[:0:0], g...)
+	for {
+		edges := make(map[string]map[triple]bool)
+		for _, s := range g {
+			for i, t := range []string{
+				s.Subject.Value,
+				s.Object.Value,
+			} {
+				e, ok := edges[t]
+				if !ok {
+					e = make(map[triple]bool)
+					edges[t] = e
+				}
+				switch i {
+				case 0:
+					e[triple{s.Predicate.Value, s.Object.Value, "+"}] = true
+				case 1:
+					e[triple{s.Predicate.Value, s.Subject.Value, "-"}] = true
+				}
+			}
+		}
+
+		seen := make(map[string]bool)
+		bNodes := make(map[string]bool)
+		terms := make(map[string]bool)
+		for _, s := range g {
+			for _, t := range []string{
+				s.Subject.Value,
+				s.Predicate.Value,
+				s.Object.Value,
+			} {
+				terms[t] = true
+				if isBlank(t) {
+					bNodes[t] = true
+				} else {
+					seen[t] = true
+				}
+			}
+		}
+
+		redundant := make(map[string]bool)
+		for x := range bNodes {
+			for xp := range terms {
+				if isProperSubset(edges[x], edges[xp]) || (seen[xp] && isEqualEdges(edges[x], edges[xp])) {
+					redundant[x] = true
+					break
+				}
+			}
+			seen[x] = true
+		}
+
+		n := len(g)
+		for i := 0; i < len(g); {
+			if !redundant[g[i].Subject.Value] && !redundant[g[i].Object.Value] {
+				i++
+				continue
+			}
+			g[i], g = g[len(g)-1], g[:len(g)-1]
+		}
+		if n == len(g) {
+			return g
+		}
+	}
+}
+
+type triple [3]string
+
+func isProperSubset(a, b map[triple]bool) bool {
+	for k := range a {
+		if !b[k] {
+			return false
+		}
+	}
+	return len(a) < len(b)
+}
+
+func isEqualEdges(a, b map[triple]bool) bool {
+	if len(a) != len(b) {
+		return false
+	}
+	for k := range a {
+		if !b[k] {
+			return false
+		}
+	}
+	return true
+}
+
+// findCandidates finds candidates for blank nodes and blank nodes that are fixed.
+//
+// This is algorithm 5 in doi:10.1145/3068333.
+func findCandidates(g []*Statement) ([]*Statement, map[string]bool, map[string]map[string]bool, bool) {
+	g = removeRedundantBnodes(g)
+
+	edges := make(map[triple]bool)
+	f := make(map[string]bool)
+	for _, s := range g {
+		sub := s.Subject.Value
+		prd := s.Predicate.Value
+		obj := s.Object.Value
+
+		edges[triple{sub, prd, obj}] = true
+		edges[triple{sub, prd, "*"}] = true
+		edges[triple{"*", prd, obj}] = true
+		switch {
+		case isBlank(sub) && isBlank(obj):
+			f[sub] = false
+			f[obj] = false
+		case isBlank(sub):
+			if _, ok := f[sub]; !ok {
+				f[sub] = true
+			}
+		case isBlank(obj):
+			if _, ok := f[obj]; !ok {
+				f[obj] = true
+			}
+		}
+	}
+	for k, v := range f {
+		if !v {
+			delete(f, k)
+		}
+	}
+	if len(f) == 0 {
+		f = nil
+	}
+
+	cands := make(map[string]map[string]bool)
+	bnodes := make(map[string]bool)
+	for _, s := range g {
+		for _, b := range []string{
+			s.Subject.Value,
+			s.Object.Value,
+		} {
+			if !isBlank(b) {
+				continue
+			}
+			bnodes[b] = true
+			if f[b] {
+				cands[b] = map[string]bool{b: true}
+			} else {
+				terms := make(map[string]bool)
+				for _, s := range g {
+					for _, t := range []string{
+						s.Subject.Value,
+						s.Predicate.Value,
+						s.Object.Value,
+					} {
+						terms[t] = true
+					}
+				}
+				cands[b] = terms
+			}
+		}
+	}
+	if isEqualTerms(f, bnodes) {
+		return g, f, cands, true
+	}
+
+	for {
+		bb := make(map[string]bool)
+		for b := range bnodes {
+			if !f[b] {
+				bb[b] = true
+			}
+		}
+		for b := range bb {
+			for x := range cands[b] {
+				if x == b {
+					continue
+				}
+				for _, s := range g {
+					if s.Subject.Value != b {
+						continue
+					}
+					prd := s.Predicate.Value
+					obj := s.Object.Value
+					if (inILF(obj, f) && !edges[triple{x, prd, obj}]) || (bb[obj] && !edges[triple{x, prd, "*"}]) {
+						delete(cands[b], x)
+						break
+					}
+				}
+				if !cands[b][x] {
+					continue
+				}
+				for _, s := range g {
+					if s.Object.Value != b {
+						continue
+					}
+					sub := s.Subject.Value
+					prd := s.Predicate.Value
+					if (inIF(sub, f) && !edges[triple{sub, prd, x}]) || (bb[sub] && !edges[triple{"*", prd, x}]) {
+						delete(cands[b], x)
+						break
+					}
+				}
+			}
+		}
+
+		fp := f
+		f = make(map[string]bool)
+		for b := range fp {
+			f[b] = true
+		}
+		for b := range bb { // Mark newly fixed blank nodes.
+			if len(cands[b]) == 1 && cands[b][b] {
+				f[b] = true
+			}
+		}
+		allFixed := isEqualTerms(f, bnodes)
+		if isEqualTerms(fp, f) || allFixed {
+			if len(f) == 0 {
+				f = nil
+			}
+			return g, f, cands, allFixed
+		}
+	}
+}
+
+// inILF returns whether t is in IL or F.
+func inILF(t string, f map[string]bool) bool {
+	return isIRI(t) || isLiteral(t) || f[t]
+}
+
+// inIF returns whether t is in I or F.
+func inIF(t string, f map[string]bool) bool {
+	return isIRI(t) || f[t]
+}
+
+// dfs is a depth-first search strategy.
+type dfs struct{}
+
+// lean returns a core of the RDF graph g using the given strategy.
+//
+// This is lines 1-9 of algorithm 6 in doi:10.1145/3068333.
+func lean(strategy *dfs, g []*Statement) []*Statement {
+	foundBnode := false
+search:
+	for _, s := range g {
+		for _, t := range []string{
+			s.Subject.Value,
+			s.Object.Value,
+		} {
+			if isBlank(t) {
+				foundBnode = true
+				break search
+			}
+		}
+	}
+	if !foundBnode {
+		return g
+	}
+	g, fixed, cands, allFixed := findCandidates(g)
+	if allFixed {
+		return g
+	}
+	for _, s := range g {
+		if isBlank(s.Subject.Value) && isBlank(s.Object.Value) {
+			mu := make(map[string]string, len(fixed))
+			for b := range fixed {
+				mu[b] = b
+			}
+			mu = findCoreEndomorphism(strategy, g, cands, mu)
+			return applyMu(g, mu)
+		}
+	}
+	return g
+}
+
+// findCoreEndomorphism returns a core solution using the given strategy.
+//
+// This is lines 10-14 of algorithm 6 in doi:10.1145/3068333.
+func findCoreEndomorphism(strategy *dfs, g []*Statement, cands map[string]map[string]bool, mu map[string]string) map[string]string {
+	var q []*Statement
+	preds := make(map[string]int)
+	seen := make(map[triple]bool)
+	for _, s := range g {
+		preds[s.Predicate.Value]++
+		if isBlank(s.Subject.Value) && isBlank(s.Object.Value) {
+			if seen[triple{s.Subject.Value, s.Predicate.Value, s.Object.Value}] {
+				continue
+			}
+			seen[triple{s.Subject.Value, s.Predicate.Value, s.Object.Value}] = true
+			q = append(q, s)
+		}
+	}
+	sort.Slice(q, func(i, j int) bool {
+		return selectivity(q[i], cands, preds) < selectivity(q[j], cands, preds)
+	})
+	return strategy.evaluate(g, q, cands, mu)
+}
+
+// selectivity returns the selectivity heuristic score for s. Lower scores
+// are more selective.
+func selectivity(s *Statement, cands map[string]map[string]bool, preds map[string]int) int {
+	return min(len(cands[s.Subject.Value])*len(cands[s.Object.Value]), preds[s.Predicate.Value])
+}
+
+// evaluate returns an endomorphism using a DFS strategy.
+//
+// This is lines 25-32 of algorithm 6 in doi:10.1145/3068333.
+func (st *dfs) evaluate(g, q []*Statement, cands map[string]map[string]bool, mu map[string]string) map[string]string {
+	mu = st.search(g, q, cands, mu)
+	for len(mu) != len(codom(mu)) {
+		mupp := fixedFrom(cands)
+		mup := findCoreEndomorphism(st, applyMu(g, mu), cands, mupp)
+		if isAutomorphism(mup) {
+			return mu
+		}
+		for b, x := range mu {
+			if _, ok := mup[b]; !ok {
+				mup[b] = x
+			}
+		}
+		mu = mup
+	}
+	return mu
+}
+
+func fixedFrom(cands map[string]map[string]bool) map[string]string {
+	fixed := make(map[string]string)
+	for b, m := range cands {
+		if len(m) == 1 && m[b] {
+			fixed[b] = b
+		}
+	}
+	return fixed
+}
+
+// applyMu applies mu to g returning the result.
+func applyMu(g []*Statement, mu map[string]string) []*Statement {
+	back := make([]Statement, 0, len(g))
+	dst := make([]*Statement, 0, len(g))
+	seen := make(map[Statement]bool)
+	for _, s := range g {
+		n := Statement{
+			Subject:   Term{Value: translate(s.Subject.Value, mu)},
+			Predicate: Term{Value: s.Predicate.Value},
+			Object:    Term{Value: translate(s.Object.Value, mu)},
+			Label:     Term{Value: s.Label.Value},
+		}
+		if seen[n] {
+			continue
+		}
+		seen[n] = true
+		back = append(back, n)
+		dst = append(dst, &back[len(back)-1])
+	}
+	return dst
+}
+
+// search returns a minimum endomorphism using a DFS strategy.
+//
+// This is lines 33-46 of algorithm 6 in doi:10.1145/3068333.
+func (st *dfs) search(g, q []*Statement, cands map[string]map[string]bool, mu map[string]string) map[string]string {
+	qMin := q[0]
+	m := st.join(qMin, g, cands, mu)
+	if len(m) == 0 {
+		// Early exit if no mapping found.
+		return nil
+	}
+	sortByCodom(m)
+	mMin := m[0]
+	qp := q[1:]
+	if len(qp) != 0 {
+		for len(m) != 0 {
+			mMin = m[0]
+			mup := st.search(g, qp, cands, mMin)
+			if !isAutomorphism(mup) {
+				return mup
+			}
+			m = m[1:]
+		}
+	}
+	return mMin
+}
+
+// isAutomorphism returns whether mu is an automorphism, this is equivalent to
+// dom(mu) == codom(mu).
+func isAutomorphism(mu map[string]string) bool {
+	return isEqualTerms(dom(mu), codom(mu))
+}
+
+// dom returns the domain of mu.
+func dom(mu map[string]string) map[string]bool {
+	d := make(map[string]bool, len(mu))
+	for v := range mu {
+		d[v] = true
+	}
+	return d
+}
+
+// codom returns the codomain of mu.
+func codom(mu map[string]string) map[string]bool {
+	cd := make(map[string]bool, len(mu))
+	for _, v := range mu {
+		cd[v] = true
+	}
+	return cd
+}
+
+// isEqualTerms returns whether a and b are identical.
+func isEqualTerms(a, b map[string]bool) bool {
+	if len(a) != len(b) {
+		return false
+	}
+	for k := range a {
+		if !b[k] {
+			return false
+		}
+	}
+	return true
+}
+
+// sortByCodom performs a sort of maps ordered by fewest blank nodes in
+// codomain, then fewest self mappings.
+func sortByCodom(maps []map[string]string) {
+	m := orderedByCodom{
+		maps:  maps,
+		attrs: make([]attrs, len(maps)),
+	}
+	for i, mu := range maps {
+		m.attrs[i].blanks = make(map[string]bool)
+		for x, y := range mu {
+			if isBlank(y) {
+				m.attrs[i].blanks[y] = true
+			}
+			if x == y {
+				m.attrs[i].selfs++
+			}
+		}
+	}
+	sort.Sort(m)
+}
+
+type orderedByCodom struct {
+	maps  []map[string]string
+	attrs []attrs
+}
+
+type attrs struct {
+	blanks map[string]bool
+	selfs  int
+}
+
+func (m orderedByCodom) Len() int { return len(m.maps) }
+func (m orderedByCodom) Less(i, j int) bool {
+	attrI := m.attrs[i]
+	attrJ := m.attrs[j]
+	switch {
+	case len(attrI.blanks) < len(attrJ.blanks):
+		return true
+	case len(attrI.blanks) > len(attrJ.blanks):
+		return false
+	default:
+		return attrI.selfs < attrJ.selfs
+	}
+}
+func (m orderedByCodom) Swap(i, j int) {
+	m.maps[i], m.maps[j] = m.maps[j], m.maps[i]
+	m.attrs[i], m.attrs[j] = m.attrs[j], m.attrs[i]
+}
+
+// join evaluates the given pattern, q, joining with solutions in m.
+// This takes only a single mapping and so only works for the DFS strategy.
+//
+// This is lines 47-51 of algorithm 6 in doi:10.1145/3068333.
+func (st *dfs) join(q *Statement, g []*Statement, cands map[string]map[string]bool, m map[string]string) []map[string]string {
+	var mp []map[string]string
+	isLoop := q.Subject.Value == q.Object.Value
+	for _, s := range g {
+		// Line 45: M_q ← {µ | µ(q) ∈ G}
+		//  | µ(q) ∈ G
+		//
+		//    µ(q) ∈ G ↔ (µ(q_s),q_p,µ(q_o)) ∈ G
+		if q.Predicate.Value != s.Predicate.Value {
+			continue
+		}
+		//    q_s = q_o ↔ µ(q_s) =_µ(q_o)
+		if isLoop && s.Subject.Value != s.Object.Value {
+			continue
+		}
+
+		// Line 46: M_q' ← {µ ∈ M_q | for all b ∈ bnodes({q}), µ(b) ∈ cands[b]}
+		//  | for all b ∈ bnodes({q}), µ(b) ∈ cands[b]
+		if !cands[q.Subject.Value][s.Subject.Value] || !cands[q.Object.Value][s.Object.Value] {
+			continue
+		}
+
+		// Line 47: M' ← M_q' ⋈ M
+		// M₁ ⋈ M₂ = {μ₁ ∪ μ₂ | μ₁ ∈ M₁, μ₂ ∈ M₂ and μ₁, μ₂ are compatible mappings}
+		//  | μ₁ ∈ M₁, μ₂ ∈ M₂ and μ₁, μ₂ are compatible mappings
+		if mq, ok := m[q.Subject.Value]; ok && mq != s.Subject.Value {
+			continue
+		}
+		if !isLoop {
+			if mq, ok := m[q.Object.Value]; ok && mq != s.Object.Value {
+				continue
+			}
+		}
+		// Line 47: μ₁ ∪ μ₂
+		var mu map[string]string
+		if isLoop {
+			mu = map[string]string{
+				q.Subject.Value: s.Subject.Value,
+			}
+		} else {
+			mu = map[string]string{
+				q.Subject.Value: s.Subject.Value,
+				q.Object.Value:  s.Object.Value,
+			}
+		}
+		for b, mb := range m {
+			mu[b] = mb
+		}
+		mp = append(mp, mu)
+	}
+	return mp
+}
diff --git a/graph/formats/rdf/rdf.go b/graph/formats/rdf/rdf.go
index e82f767c..d8af33a1 100644
--- a/graph/formats/rdf/rdf.go
+++ b/graph/formats/rdf/rdf.go
@@ -134,6 +134,10 @@ func checkIRIText(iri string) error {
 	}
 }
 
+func isLiteral(s string) bool {
+	return strings.HasPrefix(s, `"`) && strings.HasSuffix(s, `"`)
+}
+
 // Parts returns the pars of the term and the kind of the term.
 // IRI node text is returned as a valid IRI with the quoting angle
 // brackets removed and escape sequences interpreted, and blank
diff --git a/lapack/gonum/dorgr2.go b/lapack/gonum/dorgr2.go
index 7011565c..1cdb6e80 100644
--- a/lapack/gonum/dorgr2.go
+++ b/lapack/gonum/dorgr2.go
@@ -9,39 +9,34 @@ import (
 	"gonum.org/v1/gonum/blas/blas64"
 )
 
-// Dorgr2 generates an m×n real matrix Q with orthonormal rows,
-// which is defined as the last m rows of a product of k elementary
-// reflectors of order n
+// Dorgr2 generates an m×n real matrix Q with orthonormal rows, which is defined
+// as the last m rows of a product of k elementary reflectors of order n
 //  Q = H_0 * H_1 * ... * H_{k-1}
 // as returned by Dgerqf.
 //
-// Each entry of tau contains the scalar factor of the elementary reflector H_i
-// and the length of tau must be at least k. The length of work must be at least m.
-// a is a matrix of dimensions (n,lda). On entry the [m-k+i-1]-th row must contain (counting from zero)
-// the vector which defines the elementary reflector H_i, for i = 0,1,2,...,k-1, as
-// returned by Dgerqf in the last k rows of its array argument A.
-// On exit, the m×n matrix Q.
-//  n >= m >= k >= 0
+// On entry, the (m-k+i)-th row of A must contain the vector which defines the
+// elementary reflector H_i, for i = 0,1,...,k, as returned by Dgerqf. On
+// return, A will contain the m×n matrix Q.
 //
-// Dorgr2 will panic if the conditions on input values are not met.
+// The i-th element of tau must contain the scalar factor of the elementary
+// reflector H_i, as returned by Dgerqf.
+//
+// It must hold that
+//  n >= m >= k >= 0,
+// the length of tau must be k and the length of work must be m, otherwise
+// Dorgr2 will panic.
 //
 // Dorgr2 is an internal routine. It is exported for testing purposes.
 func (impl Implementation) Dorgr2(m, n, k int, a []float64, lda int, tau, work []float64) {
 	switch {
-	case m < 0:
-		panic(mLT0)
-	case n < 0:
-		panic(nLT0)
-	case m > n:
-		panic(mGTN)
 	case k < 0:
 		panic(kLT0)
-	case k > m:
+	case m < k:
 		panic(kGTM)
+	case n < m:
+		panic(mGTN)
 	case lda < max(1, n):
 		panic(badLdA)
-	case len(work) < m:
-		panic(shortWork)
 	}
 
 	// Quick return if possible.
@@ -54,16 +49,17 @@ func (impl Implementation) Dorgr2(m, n, k int, a []float64, lda int, tau, work [
 		panic(badLenTau)
 	case len(a) < (m-1)*lda+n:
 		panic(shortA)
+	case len(work) < m:
+		panic(shortWork)
 	}
 
-	if k < m {
-		// Initialise rows 0:m-k to rows of the unit matrix.
-		for l := 0; l < m-k; l++ {
-			for j := 0; j < n; j++ {
-				a[l*lda+j] = 0
-			}
-			a[l*lda+n-m+l] = 1
+	// Initialise rows 0:m-k to rows of the unit matrix.
+	for l := 0; l < m-k; l++ {
+		row := a[l*lda : l*lda+n]
+		for j := range row {
+			row[j] = 0
 		}
+		a[l*lda+n-m+l] = 1
 	}
 	bi := blas64.Implementation()
 	for i := 0; i < k; i++ {


Test patch:
diff --git a/graph/formats/rdf/equi_canonical_example_test.go b/graph/formats/rdf/equi_canonical_example_test.go
new file mode 100644
index 00000000..202219ec
--- /dev/null
+++ b/graph/formats/rdf/equi_canonical_example_test.go
@@ -0,0 +1,155 @@
+// 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 rdf_test
+
+import (
+	"fmt"
+	"strings"
+
+	"gonum.org/v1/gonum/graph/formats/rdf"
+)
+
+func ExampleLean() {
+	for i, statements := range []string{
+		0: `
+_:author1 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:gonum .
+_:author3 <ex:contributesTo> _:gonum .
+`,
+		1: `
+_:author1 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:gonum .
+_:author3 <ex:contributesTo> _:gonum .
+_:gonum <ex:dependsOn> _:go .
+`,
+		2: `
+_:author1 <ex:contributesTo> _:go .
+_:author1 <ex:notContributesTo> _:gonum .
+_:author2 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:gonum .
+_:author3 <ex:contributesTo> _:gonum .
+_:gonum <ex:dependsOn> _:go .
+`,
+		3: `
+_:author1 <ex:is> "Alice" .
+_:author1 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:gonum .
+_:author3 <ex:contributesTo> _:gonum .
+_:gonum <ex:dependsOn> _:go .
+`,
+		4: `
+_:author1 <ex:contributesTo> _:go .
+_:author1 <ex:notContributesTo> _:gonum .
+_:author2 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:gonum .
+_:author3 <ex:contributesTo> _:gonum .
+_:author3 <ex:notContributesTo> _:go .
+_:gonum <ex:dependsOn> _:go .
+`,
+		5: `
+_:author1 <ex:is> "Alice" .
+_:author2 <ex:is> "Bob" .
+_:author1 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:gonum .
+_:author3 <ex:contributesTo> _:gonum .
+_:gonum <ex:dependsOn> _:go .
+`,
+		6: `
+_:author1 <ex:is> "Alice" .
+_:author2 <ex:is> "Bob" .
+_:author3 <ex:is> "Charlie" .
+_:author1 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:go .
+_:author2 <ex:contributesTo> _:gonum .
+_:author3 <ex:contributesTo> _:gonum .
+_:gonum <ex:dependsOn> _:go .
+`,
+	} {
+		// Decode the statement stream.
+		dec := rdf.NewDecoder(strings.NewReader(statements))
+		var s []*rdf.Statement
+		for {
+			l, err := dec.Unmarshal()
+			if err != nil {
+				break
+			}
+			s = append(s, l)
+		}
+
+		// Lean the graph to remove redundant statements.
+		lean, err := rdf.Lean(s)
+		if err != nil {
+			fmt.Println(err)
+		}
+
+		// Canonicalize the blank nodes in-place.
+		_, err = rdf.URDNA2015(lean, lean)
+		if err != nil {
+			fmt.Println(err)
+			continue
+		}
+
+		fmt.Printf("%d:\n", i)
+		for _, s := range lean {
+			fmt.Println(s)
+		}
+		fmt.Println()
+	}
+
+	// Output:
+	//
+	// 0:
+	// _:c14n0 <ex:contributesTo> _:c14n1 .
+	//
+	// 1:
+	// _:c14n0 <ex:contributesTo> _:c14n1 .
+	// _:c14n0 <ex:contributesTo> _:c14n2 .
+	// _:c14n2 <ex:dependsOn> _:c14n1 .
+	//
+	// 2:
+	// _:c14n0 <ex:contributesTo> _:c14n1 .
+	// _:c14n0 <ex:contributesTo> _:c14n3 .
+	// _:c14n2 <ex:contributesTo> _:c14n1 .
+	// _:c14n2 <ex:notContributesTo> _:c14n3 .
+	// _:c14n3 <ex:dependsOn> _:c14n1 .
+	//
+	// 3:
+	// _:c14n0 <ex:contributesTo> _:c14n1 .
+	// _:c14n0 <ex:contributesTo> _:c14n3 .
+	// _:c14n2 <ex:contributesTo> _:c14n1 .
+	// _:c14n2 <ex:is> "Alice" .
+	// _:c14n3 <ex:dependsOn> _:c14n1 .
+	//
+	// 4:
+	// _:c14n0 <ex:contributesTo> _:c14n1 .
+	// _:c14n0 <ex:contributesTo> _:c14n2 .
+	// _:c14n2 <ex:dependsOn> _:c14n1 .
+	// _:c14n3 <ex:contributesTo> _:c14n1 .
+	// _:c14n3 <ex:notContributesTo> _:c14n2 .
+	// _:c14n4 <ex:contributesTo> _:c14n2 .
+	// _:c14n4 <ex:notContributesTo> _:c14n1 .
+	//
+	// 5:
+	// _:c14n1 <ex:contributesTo> _:c14n0 .
+	// _:c14n1 <ex:contributesTo> _:c14n3 .
+	// _:c14n1 <ex:is> "Bob" .
+	// _:c14n2 <ex:contributesTo> _:c14n0 .
+	// _:c14n2 <ex:is> "Alice" .
+	// _:c14n3 <ex:dependsOn> _:c14n0 .
+	//
+	// 6:
+	// _:c14n0 <ex:dependsOn> _:c14n1 .
+	// _:c14n2 <ex:contributesTo> _:c14n0 .
+	// _:c14n2 <ex:contributesTo> _:c14n1 .
+	// _:c14n2 <ex:is> "Bob" .
+	// _:c14n3 <ex:contributesTo> _:c14n1 .
+	// _:c14n3 <ex:is> "Alice" .
+	// _:c14n4 <ex:contributesTo> _:c14n0 .
+	// _:c14n4 <ex:is> "Charlie" .
+}
diff --git a/graph/formats/rdf/equi_canonical_test.go b/graph/formats/rdf/equi_canonical_test.go
new file mode 100644
index 00000000..c7b7b59e
--- /dev/null
+++ b/graph/formats/rdf/equi_canonical_test.go
@@ -0,0 +1,877 @@
+// 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 rdf
+
+import (
+	"fmt"
+	"io"
+	"reflect"
+	"strings"
+	"testing"
+)
+
+func TestRemoveRedundantNodes(t *testing.T) {
+	tests := []struct {
+		name       string
+		statements string
+		want       string
+	}{
+		{
+			name: "Example 5.1",
+			statements: `
+<ex:Chile> <ex:cabinet> _:b1 .
+<ex:Chile> <ex:cabinet> _:b2 .
+<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:cabinet> _:b4 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b2 <ex:members> "23" .
+_:b3 <ex:members> "23" .
+`,
+			want: `<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b3 <ex:members> "23" .
+`,
+		},
+		{
+			name: "Example 5.2",
+			statements: `
+_:a <ex:p> _:b .
+_:a <ex:p> _:d .
+_:c <ex:p> _:b .
+_:c <ex:p> _:f .
+_:e <ex:p> _:b .
+_:e <ex:p> _:d .
+_:e <ex:p> _:f .
+_:e <ex:p> _:h .
+_:g <ex:p> _:d .
+_:g <ex:p> _:h .
+_:i <ex:p> _:f .
+_:i <ex:p> _:h .
+`,
+			want: `_:e <ex:p> _:b .
+`,
+		},
+	}
+
+	for _, test := range tests {
+		g := parseStatements(strings.NewReader(test.statements))
+		gWant := parseStatements(strings.NewReader(test.want))
+
+		g = removeRedundantBnodes(g)
+
+		got := canonicalStatements(g)
+		want := canonicalStatements(gWant)
+		if got != want {
+			got = formatStatements(g)
+			t.Errorf("unexpected result for %s:\ngot: \n%s\nwant:\n%s",
+				test.name, got, test.want)
+		}
+
+	}
+}
+func TestFindCandidates(t *testing.T) {
+	tests := []struct {
+		name         string
+		statements   string
+		want         string
+		wantFixed    []map[string]bool
+		wantAllFixed bool
+		wantCands    []map[string]map[string]bool
+	}{
+		{
+			name: "Example 5.1",
+			statements: `
+<ex:Chile> <ex:cabinet> _:b1 .
+<ex:Chile> <ex:cabinet> _:b2 .
+<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:cabinet> _:b4 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b2 <ex:members> "23" .
+_:b3 <ex:members> "23" .
+`,
+			want: `<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b3 <ex:members> "23" .
+`,
+			// Hence, in this particular case, we have managed to fix all blank
+			// nodes, and the graph is thus lean, and we need to go no further.
+			// In other cases we will look at later, however, some blank nodes
+			// may maintain multiple candidates.
+			//
+			// Note that there are two valid labellings of the graph since _:b2
+			// and _:b3 are not distinguishable.
+			wantFixed: []map[string]bool{
+				{
+					"_:a1": true, "_:a2": true, "_:a3": true, "_:a4": true,
+					"_:b2": true, "_:c": true,
+				},
+				{
+					"_:a1": true, "_:a2": true, "_:a3": true, "_:a4": true,
+					"_:b3": true, "_:c": true,
+				},
+			},
+			wantAllFixed: true,
+			wantCands: []map[string]map[string]bool{
+				{
+					"_:a1": {"_:a1": true},
+					"_:a2": {"_:a2": true},
+					"_:a3": {"_:a3": true},
+					"_:a4": {"_:a4": true},
+					"_:b2": {"_:b2": true},
+					"_:c":  {"_:c": true},
+				},
+				{
+					"_:a1": {"_:a1": true},
+					"_:a2": {"_:a2": true},
+					"_:a3": {"_:a3": true},
+					"_:a4": {"_:a4": true},
+					"_:b3": {"_:b3": true},
+					"_:c":  {"_:c": true},
+				},
+			},
+		},
+		{
+			name: "Example 5.6", // This is 5.1, but simplified.
+			statements: `
+<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b3 <ex:members> "23" .
+`,
+			want: `<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b3 <ex:members> "23" .
+`,
+			// Hence, in this particular case, we have managed to fix all blank
+			// nodes, and the graph is thus lean, and we need to go no further.
+			// In other cases we will look at later, however, some blank nodes
+			// may maintain multiple candidates.
+			wantFixed: []map[string]bool{{
+				"_:a1": true, "_:a2": true, "_:a3": true, "_:a4": true,
+				"_:b3": true, "_:c": true,
+			}},
+			wantAllFixed: true,
+			wantCands: []map[string]map[string]bool{{
+				"_:a1": {"_:a1": true},
+				"_:a2": {"_:a2": true},
+				"_:a3": {"_:a3": true},
+				"_:a4": {"_:a4": true},
+				"_:b3": {"_:b3": true},
+				"_:c":  {"_:c": true},
+			}},
+		},
+		{
+			name: "Example 5.9",
+			statements: `
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+		`,
+			want: `<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+`,
+			wantFixed:    []map[string]bool{nil},
+			wantAllFixed: true,
+			wantCands: []map[string]map[string]bool{{
+				"_:a1": {"_:a1": true, "_:a2": true, "_:a3": true},
+				"_:a2": {"_:a2": true, "_:a3": true},
+				"_:a3": {"_:a2": true, "_:a3": true},
+				"_:a4": {"_:a2": true, "_:a3": true, "_:a4": true},
+			}},
+		},
+		{
+			name: "Example 5.10",
+			statements: `
+_:a <ex:p> _:b .
+_:a <ex:p> _:d .
+_:b <ex:q> _:e .
+_:c <ex:p> _:b .
+_:c <ex:p> _:f .
+_:d <ex:q> _:e .
+_:f <ex:q> _:e .
+_:g <ex:p> _:d .
+_:g <ex:p> _:h .
+_:h <ex:q> _:e .
+_:i <ex:p> _:f .
+_:i <ex:p> _:h .
+`,
+			want: `_:a <ex:p> _:b .
+_:a <ex:p> _:d .
+_:b <ex:q> _:e .
+_:c <ex:p> _:b .
+_:c <ex:p> _:f .
+_:d <ex:q> _:e .
+_:f <ex:q> _:e .
+_:g <ex:p> _:d .
+_:g <ex:p> _:h .
+_:h <ex:q> _:e .
+_:i <ex:p> _:f .
+_:i <ex:p> _:h .
+`,
+			wantFixed:    []map[string]bool{{"_:e": true}},
+			wantAllFixed: false,
+			wantCands: []map[string]map[string]bool{{
+				"_:a": {"_:a": true, "_:c": true, "_:g": true, "_:i": true},
+				"_:b": {"_:b": true, "_:d": true, "_:f": true, "_:h": true},
+				"_:c": {"_:a": true, "_:c": true, "_:g": true, "_:i": true},
+				"_:d": {"_:b": true, "_:d": true, "_:f": true, "_:h": true},
+				"_:e": {"_:e": true},
+				"_:f": {"_:b": true, "_:d": true, "_:f": true, "_:h": true},
+				"_:g": {"_:a": true, "_:c": true, "_:g": true, "_:i": true},
+				"_:h": {"_:b": true, "_:d": true, "_:f": true, "_:h": true},
+				"_:i": {"_:a": true, "_:c": true, "_:g": true, "_:i": true},
+			}},
+		},
+	}
+
+	for _, test := range tests[:1] {
+		g := parseStatements(strings.NewReader(test.statements))
+		gWant := parseStatements(strings.NewReader(test.want))
+
+		g, fixed, cands, allFixed := findCandidates(g)
+
+		got := canonicalStatements(g)
+		want := canonicalStatements(gWant)
+		if got != want {
+			got = formatStatements(g)
+			t.Errorf("unexpected result for %s:\ngot: \n%s\nwant:\n%s",
+				test.name, got, test.want)
+		}
+
+		matchedFixed := false
+		for _, wantFixed := range test.wantFixed {
+			if reflect.DeepEqual(fixed, wantFixed) {
+				matchedFixed = true
+				break
+			}
+		}
+		if !matchedFixed {
+			t.Errorf("unexpected fixed result for %s:\ngot: \n%v\nwant:\n%v",
+				test.name, fixed, test.wantFixed)
+		}
+
+		if allFixed != test.wantAllFixed {
+			t.Errorf("unexpected all-fixed result for %s:\ngot:%t\nwant:%t",
+				test.name, allFixed, test.wantAllFixed)
+		}
+
+		matchedCands := false
+		for _, wantCands := range test.wantCands {
+			if reflect.DeepEqual(cands, wantCands) {
+				matchedCands = true
+				break
+			}
+		}
+		if !matchedCands {
+			t.Errorf("unexpected candidates result for %s:\ngot: \n%v\nwant:\n%v",
+				test.name, cands, test.wantCands)
+		}
+	}
+}
+
+func TestLean(t *testing.T) {
+	var tests = []struct {
+		name       string
+		statements string
+		want       string
+		wantErr    error
+	}{
+		{
+			name: "Example 5.1",
+			statements: `
+<ex:Chile> <ex:cabinet> _:b1 .
+<ex:Chile> <ex:cabinet> _:b2 .
+<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:cabinet> _:b4 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b2 <ex:members> "23" .
+_:b3 <ex:members> "23" .
+`,
+			want: `<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b3 <ex:members> "23" .
+`,
+		},
+		{
+			name: "Example 5.6",
+			statements: `
+<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b3 <ex:members> "23" .
+`,
+			want: `<ex:Chile> <ex:cabinet> _:b3 .
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+<ex:MBachelet> <ex:spouse> _:c .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a2 <ex:president> <ex:MBachelet> .
+_:a3 <ex:next> _:a4 .
+_:a4 <ex:president> <ex:MBachelet> .
+_:b3 <ex:members> "23" .
+`,
+		},
+		{
+			name: "Example 5.9",
+			statements: `
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+		`,
+			want: `<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+`,
+		},
+		{
+			name: "Example 5.10",
+			statements: `
+_:a <ex:p> _:b .
+_:a <ex:p> _:d .
+_:b <ex:q> _:e .
+_:c <ex:p> _:b .
+_:c <ex:p> _:f .
+_:d <ex:q> _:e .
+_:f <ex:q> _:e .
+_:g <ex:p> _:d .
+_:g <ex:p> _:h .
+_:h <ex:q> _:e .
+_:i <ex:p> _:f .
+_:i <ex:p> _:h .
+`,
+			want: `_:a <ex:p> _:b .
+_:b <ex:q> _:e .
+`,
+		},
+		{
+			name: "Example 5.10 halved",
+			statements: `
+_:a <ex:p> _:b .
+_:a <ex:p> _:d .
+_:b <ex:q> _:e .
+_:c <ex:p> _:b .
+_:c <ex:p> _:f .
+_:d <ex:q> _:e .
+_:f <ex:q> _:e .
+`,
+			want: `_:a <ex:p> _:b .
+_:b <ex:q> _:e .
+`,
+		},
+		{
+			name: "Example 5.10 quartered",
+			statements: `
+_:a <ex:p> _:b .
+_:a <ex:p> _:d .
+_:b <ex:q> _:e .
+_:d <ex:q> _:e .
+`,
+			want: `_:a <ex:p> _:b .
+_:b <ex:q> _:e .
+`,
+		},
+	}
+
+	for _, test := range tests {
+		g := parseStatements(strings.NewReader(test.statements))
+		gWant := parseStatements(strings.NewReader(test.want))
+
+		lean, err := Lean(g)
+		if err != test.wantErr {
+			t.Errorf("unexpected error for %v: got:%v want:%v",
+				test.name, err, test.wantErr)
+		}
+
+		got := canonicalStatements(lean)
+		want := canonicalStatements(gWant)
+
+		if got != want {
+			got = formatStatements(g)
+			t.Errorf("unexpected result for %s:\ngot: \n%s\nwant:\n%s",
+				test.name, got, test.want)
+		}
+	}
+}
+
+func TestJoin(t *testing.T) {
+	var tests = []struct {
+		name       string
+		q          string
+		statements string
+		cands      map[string]map[string]bool
+		mu         map[string]string
+		want       []map[string]string
+	}{
+		{
+			name: "Indentity",
+			q:    `_:a <ex:p> _:b .`,
+			statements: `
+_:a <ex:p> _:b .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true},
+				"_:b": {"_:b": true},
+			},
+			mu: nil,
+			want: []map[string]string{
+				{"_:a": "_:a", "_:b": "_:b"},
+			},
+		},
+		{
+			name: "Cross identity",
+			q:    `_:a <ex:p> _:b .`,
+			statements: `
+_:a <ex:p> _:b .
+_:b <ex:p> _:a .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu: nil,
+			want: []map[string]string{
+				{"_:a": "_:a", "_:b": "_:b"},
+				{"_:a": "_:b", "_:b": "_:a"},
+			},
+		},
+		{
+			name: "Cross identity with restriction",
+			q:    `_:a <ex:p> _:b .`,
+			statements: `
+_:a <ex:p> _:b .
+_:b <ex:p> _:a .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu: map[string]string{"_:a": "_:a"},
+			want: []map[string]string{
+				{"_:a": "_:a", "_:b": "_:b"},
+			},
+		},
+		{
+			name: "Cross identity with complete restriction",
+			q:    `_:a <ex:p> _:b .`,
+			statements: `
+_:a <ex:p> _:b .
+_:b <ex:p> _:a .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu:   map[string]string{"_:a": "_:a", "_:b": "_:a"},
+			want: nil,
+		},
+		{
+			name: "Cross identity with extension",
+			q:    `_:a <ex:p> _:b .`,
+			statements: `
+_:a <ex:p> _:b .
+_:b <ex:p> _:a .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu: map[string]string{"_:c": "_:a"},
+			want: []map[string]string{
+				{"_:a": "_:a", "_:b": "_:b", "_:c": "_:a"},
+				{"_:a": "_:b", "_:b": "_:a", "_:c": "_:a"},
+			},
+		},
+
+		{
+			name: "Loop",
+			q:    `_:a <ex:p> _:a .`,
+			statements: `
+_:a <ex:p> _:a .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true},
+			},
+			mu: nil,
+			want: []map[string]string{
+				{"_:a": "_:a"},
+			},
+		},
+		{
+			name: "Cross identity loop",
+			q:    `_:a <ex:p> _:a .`,
+			statements: `
+_:a <ex:p> _:a .
+_:b <ex:p> _:b .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu: nil,
+			want: []map[string]string{
+				{"_:a": "_:a"},
+				{"_:a": "_:b"},
+			},
+		},
+		{
+			name: "Cross identity loop with restriction",
+			q:    `_:a <ex:p> _:a .`,
+			statements: `
+_:a <ex:p> _:a .
+_:b <ex:p> _:b .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu: map[string]string{"_:a": "_:a"},
+			want: []map[string]string{
+				{"_:a": "_:a"},
+			},
+		},
+		{
+			name: "Cross identity loop with complete restriction",
+			q:    `_:a <ex:p> _:a .`,
+			statements: `
+_:a <ex:p> _:a .
+_:b <ex:p> _:b .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu: map[string]string{"_:a": "_:b", "_:b": "_:a"},
+			want: []map[string]string{
+				{"_:a": "_:b", "_:b": "_:a"},
+			},
+		},
+		{
+			name: "Cross identity loop with extension",
+			q:    `_:a <ex:p> _:a .`,
+			statements: `
+_:a <ex:p> _:a .
+_:b <ex:p> _:b .
+`,
+			cands: map[string]map[string]bool{
+				"_:a": {"_:a": true, "_:b": true},
+				"_:b": {"_:b": true, "_:a": true},
+			},
+			mu: map[string]string{"_:c": "_:a"},
+			want: []map[string]string{
+				{"_:a": "_:a", "_:c": "_:a"},
+				{"_:a": "_:b", "_:c": "_:a"},
+			},
+		},
+
+		{
+			name: "Example 5.9 step 1",
+			q:    `_:a1 <ex:next> _:a2 .`,
+			statements: `
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+`,
+			cands: map[string]map[string]bool{
+				"_:a1": {"_:a1": true, "_:a2": true, "_:a3": true},
+				"_:a2": {"_:a2": true, "_:a3": true},
+				"_:a3": {"_:a2": true, "_:a3": true},
+				"_:a4": {"_:a2": true, "_:a3": true, "_:a4": true},
+			},
+			mu: map[string]string{},
+			want: []map[string]string{
+				{"_:a1": "_:a1", "_:a2": "_:a2"},
+				{"_:a1": "_:a2", "_:a2": "_:a3"},
+			},
+		},
+		{
+			name: "Example 5.9 step 2",
+			q:    `_:a2 <ex:next> _:a3 .`,
+			statements: `
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+`,
+			cands: map[string]map[string]bool{
+				"_:a1": {"_:a1": true, "_:a2": true, "_:a3": true},
+				"_:a2": {"_:a2": true, "_:a3": true},
+				"_:a3": {"_:a2": true, "_:a3": true},
+				"_:a4": {"_:a2": true, "_:a3": true, "_:a4": true},
+			},
+			mu:   map[string]string{"_:a1": "_:a2", "_:a2": "_:a3"},
+			want: nil,
+		},
+		{
+			name: "Example 5.9 step 3",
+			q:    `_:a2 <ex:next> _:a3 .`,
+			statements: `
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+`,
+			cands: map[string]map[string]bool{
+				"_:a1": {"_:a1": true, "_:a2": true, "_:a3": true},
+				"_:a2": {"_:a2": true, "_:a3": true},
+				"_:a3": {"_:a2": true, "_:a3": true},
+				"_:a4": {"_:a2": true, "_:a3": true, "_:a4": true},
+			},
+			mu: map[string]string{"_:a1": "_:a1", "_:a2": "_:a2"},
+			want: []map[string]string{
+				{"_:a1": "_:a1", "_:a2": "_:a2", "_:a3": "_:a3"},
+			},
+		},
+		{
+			name: "Example 5.9 step 4",
+			q:    `_:a3 <ex:next> _:a4 .`,
+			statements: `
+<ex:Chile> <ex:presidency> _:a1 .
+<ex:Chile> <ex:presidency> _:a2 .
+<ex:Chile> <ex:presidency> _:a3 .
+<ex:Chile> <ex:presidency> _:a4 .
+_:a1 <ex:next> _:a2 .
+_:a2 <ex:next> _:a3 .
+_:a3 <ex:next> _:a4 .
+`,
+			cands: map[string]map[string]bool{
+				"_:a1": {"_:a1": true, "_:a2": true, "_:a3": true},
+				"_:a2": {"_:a2": true, "_:a3": true},
+				"_:a3": {"_:a2": true, "_:a3": true},
+				"_:a4": {"_:a2": true, "_:a3": true, "_:a4": true},
+			},
+			mu: map[string]string{"_:a1": "_:a1", "_:a2": "_:a2", "_:a3": "_:a3"},
+			want: []map[string]string{
+				{"_:a1": "_:a1", "_:a2": "_:a2", "_:a3": "_:a3", "_:a4": "_:a4"},
+			},
+		},
+	}
+
+	for _, test := range tests {
+		q := parseStatement(strings.NewReader(test.q))
+		g := parseStatements(strings.NewReader(test.statements))
+
+		st := dfs{}
+		got := st.join(q, g, test.cands, test.mu)
+		if !reflect.DeepEqual(got, test.want) {
+			t.Errorf("unexpected result for %s:\ngot:\n%#v\nwant:\n%#v",
+				test.name, got, test.want)
+		}
+
+		naive := joinNaive(q, g, test.cands, []map[string]string{test.mu})
+		if !reflect.DeepEqual(naive, test.want) {
+			t.Errorf("unexpected naive result for %s:\ngot:\n%#v\nwant:\n%#v",
+				test.name, naive, test.want)
+		}
+
+	}
+}
+
+// joinNaive is a direct translation of lines 47-51 of algorithm 6 in doi:10.1145/3068333.
+func joinNaive(q *Statement, G []*Statement, cands map[string]map[string]bool, M []map[string]string) []map[string]string {
+	isLoop := q.Subject.Value == q.Object.Value
+	// Line 48: M_q ← {µ | µ(q) ∈ G}
+	var M_q []map[string]string
+	for _, s := range G {
+		// µ(q) ∈ G ↔ (µ(q_s),q_p,µ(q_o)) ∈ G
+		if q.Predicate.Value != s.Predicate.Value {
+			continue
+		}
+		// q_s = q_o ↔ µ(q_s) =_µ(q_o)
+		if isLoop && s.Subject.Value != s.Object.Value {
+			continue
+		}
+
+		var µ map[string]string
+		if isLoop {
+			µ = map[string]string{
+				q.Subject.Value: s.Subject.Value,
+			}
+		} else {
+			µ = map[string]string{
+				q.Subject.Value: s.Subject.Value,
+				q.Object.Value:  s.Object.Value,
+			}
+		}
+		M_q = append(M_q, µ)
+	}
+
+	// Line 49: M_q' ← {µ ∈ M_q | for all b ∈ bnodes({q}), µ(b) ∈ cands[b]}
+	var M_qPrime []map[string]string
+	for _, µ := range M_q {
+		if !cands[q.Subject.Value][µ[q.Subject.Value]] {
+			continue
+		}
+		if !cands[q.Object.Value][µ[q.Object.Value]] {
+			continue
+		}
+		M_qPrime = append(M_qPrime, µ)
+	}
+
+	// Line 50: M' ← M_q' ⋈ M
+	// M₁ ⋈ M₂ = {μ₁ ∪ μ₂ | μ₁ ∈ M₁, μ₂ ∈ M₂ and μ₁, μ₂ are compatible mappings}
+	var MPrime []map[string]string
+	for _, µ := range M {
+	join:
+		for _, µ_qPrime := range M_qPrime {
+			for b, x_qPrime := range µ_qPrime {
+				if x, ok := µ[b]; ok && x != x_qPrime {
+					continue join
+				}
+			}
+			// Line 50: μ₁ ∪ μ₂
+			for b, x := range µ {
+				µ_qPrime[b] = x
+			}
+			MPrime = append(MPrime, µ_qPrime)
+		}
+	}
+	return MPrime
+}
+
+func parseStatement(r io.Reader) *Statement {
+	g := parseStatements(r)
+	if len(g) != 1 {
+		panic(fmt.Sprintf("invalid statement stream length %d != 1", len(g)))
+	}
+	return g[0]
+}
+
+func parseStatements(r io.Reader) []*Statement {
+	var g []*Statement
+	dec := NewDecoder(r)
+	for {
+		s, err := dec.Unmarshal()
+		if err != nil {
+			if err == io.EOF {
+				break
+			}
+			panic(err)
+		}
+		g = append(g, s)
+	}
+	return g
+}
+
+func canonicalStatements(g []*Statement) string {
+	g, _ = URDNA2015(nil, g)
+	return formatStatements(g)
+}
+
+func formatStatements(g []*Statement) string {
+	var buf strings.Builder
+	for _, s := range g {
+		fmt.Fprintln(&buf, s)
+	}
+	return buf.String()
+}
diff --git a/lapack/testlapack/dorgr2.go b/lapack/testlapack/dorgr2.go
index 37def8d8..4f778d73 100644
--- a/lapack/testlapack/dorgr2.go
+++ b/lapack/testlapack/dorgr2.go
@@ -6,12 +6,14 @@ package testlapack
 
 import (
 	"fmt"
+	"math"
 	"testing"
 
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/floats"
 	"gonum.org/v1/gonum/lapack"
 )
 
@@ -26,7 +28,7 @@ func Dorgr2Test(t *testing.T, impl Dorgr2er) {
 	for _, k := range []int{0, 1, 2, 5} {
 		for _, m := range []int{k, k + 1, k + 2, k + 4} {
 			for _, n := range []int{m, m + 1, m + 2, m + 4, m + 7} {
-				for _, lda := range []int{n, n + 5} {
+				for _, lda := range []int{max(1, n), n + 5} {
 					dorgr2Test(t, impl, rnd, m, n, k, lda)
 				}
 			}
@@ -35,60 +37,72 @@ func Dorgr2Test(t *testing.T, impl Dorgr2er) {
 }
 
 func dorgr2Test(t *testing.T, impl Dorgr2er, rnd *rand.Rand, m, n, k, lda int) {
-	const tol = 1e-12
-	name := fmt.Sprintf("m=%v,n=%v,k=%v,lda=%v", m, n, k, lda)
+	const tol = 1e-14
 
-	if lda == 0 {
-		lda = n
-	}
+	name := fmt.Sprintf("m=%v,n=%v,k=%v,lda=%v", m, n, k, lda)
 
+	// Generate a random m×n matrix A.
 	a := randomGeneral(m, n, lda, rnd)
-	aCopy := cloneGeneral(a)
+
 	// Compute the RQ decomposition of A.
+	rq := cloneGeneral(a)
 	tau := make([]float64, m)
 	work := make([]float64, 1)
-	impl.Dgerqf(m, n, a.Data, a.Stride, tau, work, -1)
+	impl.Dgerqf(m, n, rq.Data, rq.Stride, tau, work, -1)
 	work = make([]float64, int(work[0]))
-	impl.Dgerqf(m, n, a.Data, a.Stride, tau, work, len(work))
+	impl.Dgerqf(m, n, rq.Data, rq.Stride, tau, work, len(work))
+
+	tauCopy := make([]float64, len(tau))
+	copy(tauCopy, tau)
 
-	// Generate the upper triangular matrix R from subarray A[m-k:m, n-m:n]
-	r := zeros(k, m, m)
-	for i := 0; i < k; i++ {
-		for j := 0; j < k; j++ {
-			ia := i + a.Rows - k
-			ja := j + a.Cols - k
-			jr := j + r.Cols - k
-			if i <= j {
-				r.Data[i*r.Stride+jr] = a.Data[ia*a.Stride+ja]
-			}
-		}
-	}
 	// Compute the matrix Q using Dorg2r.
-	impl.Dorgr2(m, n, k, a.Data, a.Stride, tau[m-k:m], work)
+	q := cloneGeneral(rq)
+	impl.Dorgr2(m, n, k, q.Data, q.Stride, tau[m-k:m], work)
+
 	if m == 0 {
 		return
 	}
-	q := a
-	// Test Q orthogonality.
+
+	// Check that tau hasn't been modified.
+	if !floats.Equal(tau, tauCopy) {
+		t.Errorf("%v: unexpected modification in tau", name)
+	}
+
+	// Check that Q has orthonormal rows.
 	res := residualOrthogonal(q, true)
-	if res > tol {
-		t.Errorf("%v: |I - Q * Qᵀ| residual too large (%g)", name, res)
+	if res > tol || math.IsNaN(res) {
+		t.Errorf("%v: residual |I - Q*Qᵀ| too large, got %v, want <= %v", name, res, tol)
 	}
+
 	if k == 0 {
 		return
 	}
 
-	// Reconstruct last rows of A.
+	// Extract the k×m upper triangular matrix R from RQ[m-k:m,n-k:n].
+	r := zeros(k, m, m)
+	for i := 0; i < k; i++ {
+		for j := 0; j < k; j++ {
+			ii := rq.Rows - k + i
+			jj := rq.Cols - k + j
+			jr := r.Cols - k + j
+			if i <= j {
+				r.Data[i*r.Stride+jr] = rq.Data[ii*rq.Stride+jj]
+			}
+		}
+	}
+
+	// Construct a view A[m-k:m,0:n] of the last k rows of A.
 	aRec := blas64.General{
 		Rows:   k,
-		Cols:   aCopy.Cols,
-		Data:   aCopy.Data[(a.Rows-k)*lda:],
-		Stride: lda,
+		Cols:   n,
+		Data:   a.Data[(m-k)*a.Stride:],
+		Stride: a.Stride,
 	}
-	// Test |A[m-k:m,0:n] - R[m-k:m,0:m] * Q| is small
-	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, r, q, -1, aRec)
+	// Compute A - R*Q.
+	blas64.Gemm(blas.NoTrans, blas.NoTrans, -1, r, q, 1, aRec)
+	// Check that |A - R*Q| is small.
 	res = dlange(lapack.MaxColumnSum, aRec.Rows, aRec.Cols, aRec.Data, aRec.Stride)
-	if res > tol {
-		t.Errorf("%v: |A[m-k:m,0:n] - R[m-k:m,0:m] * Q| residual too large (%g)", name, res)
+	if res > tol || math.IsNaN(res) {
+		t.Errorf("%v: residual |A - R*Q| too large, got %v, want <= %v", name, res, tol)
 	}
 }