Repository: gonum/gonum
Pull request #1662 (2021-06-11 20:07:41)

interp: add ClampedCubic interpolator
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/graph/formats/rdf/iso_canonical.go b/graph/formats/rdf/iso_canonical.go
index 7f66ce1b..2f31e6a8 100644
--- a/graph/formats/rdf/iso_canonical.go
+++ b/graph/formats/rdf/iso_canonical.go
@@ -10,7 +10,6 @@ import (
 	"fmt"
 	"hash"
 	"sort"
-	"strings"
 )
 
 // See "Canonical Forms for Isomorphic and Equivalent RDF Graphs: Algorithms
@@ -160,7 +159,7 @@ func C14n(dst, src []*Statement, terms map[string]map[string]bool) ([]*Statement
 			s.Object.Value,
 			s.Label.Value,
 		} {
-			if !strings.HasPrefix(t, "_:") {
+			if !isBlank(t) {
 				continue
 			}
 			need[t] = true
@@ -172,7 +171,7 @@ func C14n(dst, src []*Statement, terms map[string]map[string]bool) ([]*Statement
 	for h, m := range terms {
 		var ok bool
 		for t := range m {
-			if strings.HasPrefix(t, "_:") {
+			if isBlank(t) {
 				ok = true
 				break
 			}
@@ -194,7 +193,7 @@ func C14n(dst, src []*Statement, terms map[string]map[string]bool) ([]*Statement
 			return nil, fmt.Errorf("rdf: no term for blank with hash %x", b)
 		}
 		for t := range terms[b] {
-			if !strings.HasPrefix(t, "_:") {
+			if !isBlank(t) {
 				continue
 			}
 			if _, exists := c14n[t]; exists {
@@ -275,13 +274,13 @@ func hashBNodes(statements []*Statement, h hash.Hash, zero []byte, hash0 map[str
 			switch {
 			case i == 3 && t == "":
 				continue
-			case strings.HasPrefix(t, "_:"):
+			case isBlank(t):
 				if hash0 == nil {
 					curr.set(t, zero)
 				} else {
 					curr.set(t, hash0[t])
 				}
-			case strings.HasPrefix(t, "<") && strings.HasSuffix(t, ">"):
+			case isIRI(t):
 				h.Reset()
 				h.Write([]byte(t[1 : len(t)-1])) //nolint:errcheck
 				curr.set(t, h.Sum(nil))
@@ -298,7 +297,7 @@ func hashBNodes(statements []*Statement, h hash.Hash, zero []byte, hash0 map[str
 	for {
 		curr, last = last, curr
 		for _, s := range statements {
-			if strings.HasPrefix(s.Subject.Value, "_:") {
+			if isBlank(s.Subject.Value) {
 				var lab []byte
 				if s.Label.Value != "" {
 					lab = last.hashOf[s.Label.Value]
@@ -307,7 +306,7 @@ func hashBNodes(statements []*Statement, h hash.Hash, zero []byte, hash0 map[str
 				bag.add(s.Subject.Value, c)
 			}
 
-			if strings.HasPrefix(s.Object.Value, "_:") {
+			if isBlank(s.Object.Value) {
 				var lab []byte
 				if s.Label.Value != "" {
 					lab = last.hashOf[s.Label.Value]
@@ -320,7 +319,7 @@ func hashBNodes(statements []*Statement, h hash.Hash, zero []byte, hash0 map[str
 			// required for RDF dataset hashing as described in
 			// https://doi.org/10.5281/zenodo.3154322 v1.0
 			// Readme.md#adaptation-of-the-algorithms-to-handle-datasets.
-			if strings.HasPrefix(s.Label.Value, "_:") {
+			if isBlank(s.Label.Value) {
 				c := hashTuple(h, last.hashOf[s.Subject.Value], last.hashOf[s.Predicate.Value], last.hashOf[s.Object.Value], []byte{'.'})
 				bag.add(s.Label.Value, c)
 			}
@@ -430,7 +429,7 @@ func (t *table) set(term string, hash []byte) {
 		t.termsFor[string(hash)] = map[string]bool{term: true}
 	}
 
-	if !t.wasCloned() && strings.HasPrefix(term, "_:") {
+	if !t.wasCloned() && isBlank(term) {
 		// We are in the original table, so note
 		// any blank node label that we see.
 		t.isBlank[term] = true
@@ -571,7 +570,7 @@ func appendOrdered(parts byLengthHash, partSets map[string]map[string]bool) byLe
 	for h, s := range partSets {
 		var p []string
 		for e := range s {
-			if strings.HasPrefix(e, "_:") {
+			if isBlank(e) {
 				p = append(p, e)
 			}
 		}
@@ -631,7 +630,7 @@ func split(statements []*Statement) [][]*Statement {
 	for _, s := range statements {
 		ds.add(s.Subject.Value)
 		ds.add(s.Object.Value)
-		if strings.HasPrefix(s.Subject.Value, "_:") && strings.HasPrefix(s.Object.Value, "_:") {
+		if isBlank(s.Subject.Value) && isBlank(s.Object.Value) {
 			ds.union(ds.find(s.Subject.Value), ds.find(s.Object.Value))
 		}
 	}
@@ -644,9 +643,9 @@ func split(statements []*Statement) [][]*Statement {
 	for _, s := range statements {
 		var t string
 		switch {
-		case strings.HasPrefix(s.Subject.Value, "_:"):
+		case isBlank(s.Subject.Value):
 			t = s.Subject.Value
-		case strings.HasPrefix(s.Object.Value, "_:"):
+		case isBlank(s.Object.Value):
 			t = s.Object.Value
 		default:
 			ground = append(ground, s)
@@ -929,15 +928,15 @@ type relabeledTerm struct {
 }
 
 func (a relabeledTerm) less(b relabeledTerm) bool {
-	aBlank := strings.HasPrefix(a.term.Value, "_:")
-	bBlank := strings.HasPrefix(b.term.Value, "_:")
+	aIsBlank := isBlank(a.term.Value)
+	bIsBlank := isBlank(b.term.Value)
 	switch {
-	case aBlank && bBlank:
+	case aIsBlank && bIsBlank:
 		return bytes.Compare(a.labels[a.term.Value], b.labels[b.term.Value]) < 0
-	case aBlank:
-		return "_:" < unquoteIRI(b.term.Value)
-	case bBlank:
-		return unquoteIRI(a.term.Value) < "_:"
+	case aIsBlank:
+		return blankPrefix < unquoteIRI(b.term.Value)
+	case bIsBlank:
+		return unquoteIRI(a.term.Value) < blankPrefix
 	default:
 		return unquoteIRI(a.term.Value) < unquoteIRI(b.term.Value)
 	}
@@ -951,7 +950,7 @@ func unquoteIRI(s string) string {
 }
 
 func (t relabeledTerm) String() string {
-	if !strings.HasPrefix(t.term.Value, "_:") {
+	if !isBlank(t.term.Value) {
 		return t.term.Value
 	}
 	h, ok := t.labels[t.term.Value]
diff --git a/graph/formats/rdf/rdf.go b/graph/formats/rdf/rdf.go
index 63e4ddc3..e82f767c 100644
--- a/graph/formats/rdf/rdf.go
+++ b/graph/formats/rdf/rdf.go
@@ -75,7 +75,13 @@ func NewBlankTerm(label string) (Term, error) {
 	if err != nil {
 		return Term{}, err
 	}
-	return Term{Value: "_:" + label}, nil
+	return Term{Value: blankPrefix + label}, nil
+}
+
+const blankPrefix = "_:"
+
+func isBlank(s string) bool {
+	return strings.HasPrefix(s, blankPrefix)
 }
 
 // NewIRITerm returns a Term based on the provided IRI which must
@@ -89,6 +95,10 @@ func NewIRITerm(iri string) (Term, error) {
 	return Term{Value: escape("<", iri, ">")}, nil
 }
 
+func isIRI(s string) bool {
+	return strings.HasPrefix(s, "<") && strings.HasSuffix(s, ">")
+}
+
 // NewLiteralTerm returns a Term based on the literal text and an
 // optional qualifier which may either be a "@"-prefixed language
 // tag or a valid IRI. The text will be escaped if necessary and quoted,
diff --git a/graph/formats/rdf/urna.go b/graph/formats/rdf/urna.go
index c51f17d9..a7b6dc5a 100644
--- a/graph/formats/rdf/urna.go
+++ b/graph/formats/rdf/urna.go
@@ -12,7 +12,6 @@ import (
 	"fmt"
 	"hash"
 	"sort"
-	"strings"
 
 	"gonum.org/v1/gonum/stat/combin"
 )
@@ -112,7 +111,7 @@ func (u *urna) relabel(dst, src []*Statement) ([]*Statement, error) {
 			s.Object.Value,
 			s.Label.Value,
 		} {
-			if !strings.HasPrefix(t, "_:") {
+			if !isBlank(t) {
 				continue
 			}
 			for _, e := range u.statementsFor[t] {
@@ -258,7 +257,7 @@ func (u *urna) hashFirstDegreeQuads(b string) string {
 // replaceBlank implements 3.1 of the algorithm described at
 // https://json-ld.github.io/rdf-dataset-canonicalization/spec/index.html#algorithm-1.
 func replaceBlank(b, matching, label string) string {
-	if !strings.HasPrefix(b, "_:") { // 3,1
+	if !isBlank(b) { // 3.1
 		return b
 	}
 	if label != "" { // URGNA2012 modification.
@@ -372,7 +371,7 @@ func (u *urna) hashToRelatedURDNA2015(b string, names *issuer) map[string][]stri
 			s.Object.Value,
 			s.Label.Value,
 		} {
-			if !strings.HasPrefix(term, "_:") || term == b {
+			if !isBlank(term) || term == b {
 				continue
 			}
 
@@ -404,10 +403,10 @@ func (u *urna) hashToRelatedURGNA2012(b string, names *issuer) map[string][]stri
 			pos  byte
 		)
 		switch {
-		case strings.HasPrefix(s.Subject.Value, "_:") && s.Subject.Value != b: // 1.1
+		case isBlank(s.Subject.Value) && s.Subject.Value != b: // 1.1
 			term = s.Subject.Value
 			pos = 'p'
-		case strings.HasPrefix(s.Object.Value, "_:") && s.Object.Value != b: // 1.2
+		case isBlank(s.Object.Value) && s.Object.Value != b: // 1.2
 			term = s.Object.Value
 			pos = 'r'
 		default:
diff --git a/interp/cubic.go b/interp/cubic.go
index 4afb5e59..e6f0ed0a 100644
--- a/interp/cubic.go
+++ b/interp/cubic.go
@@ -377,8 +377,9 @@ func makeCubicSplineSecondDerivativeEquations(xs, ys []float64) (*mat.Tridiag, m
 
 // NaturalCubic is a piecewise cubic 1-dimensional interpolator with
 // continuous value, first and second derivatives, which can be fitted to (X, Y)
-// value pairs without providing derivatives. See e.g. https://www.math.drexel.edu/~tolya/cubicspline.pdf
-// for details.
+// value pairs without providing derivatives. It uses the boundary conditions
+// Y''(left end ) = Y''(right end) = 0.
+// See e.g. https://www.math.drexel.edu/~tolya/cubicspline.pdf for details.
 type NaturalCubic struct {
 	cubic PiecewiseCubic
 }
@@ -412,3 +413,48 @@ func (nc *NaturalCubic) Fit(xs, ys []float64) error {
 	}
 	return err
 }
+
+// ClampedCubic is a piecewise cubic 1-dimensional interpolator with
+// continuous value, first and second derivatives, which can be fitted to (X, Y)
+// value pairs without providing derivatives. It uses the boundary conditions
+// Y'(left end ) = Y'(right end) = 0.
+type ClampedCubic struct {
+	cubic PiecewiseCubic
+}
+
+// Predict returns the interpolation value at x.
+func (cc *ClampedCubic) Predict(x float64) float64 {
+	return cc.cubic.Predict(x)
+}
+
+// PredictDerivative returns the predicted derivative at x.
+func (cc *ClampedCubic) PredictDerivative(x float64) float64 {
+	return cc.cubic.PredictDerivative(x)
+}
+
+// Fit fits a predictor to (X, Y) value pairs provided as two slices.
+// It panics if len(xs) < 2, elements of xs are not strictly increasing
+// or len(xs) != len(ys). It returns an error if solving the required system
+// of linear equations fails.
+func (cc *ClampedCubic) Fit(xs, ys []float64) error {
+	a, b := makeCubicSplineSecondDerivativeEquations(xs, ys)
+	// Add boundary conditions y''(left) = y''(right) = 0:
+	n := len(xs)
+	// Condition Y'(left end) = 0:
+	dxL := xs[1] - xs[0]
+	b.SetVec(0, (ys[1]-ys[0])/dxL)
+	a.SetBand(0, 0, dxL/3)
+	a.SetBand(0, 1, dxL/6)
+	// Condition Y'(right end) = 0:
+	m := n - 1
+	dxR := xs[m] - xs[m-1]
+	b.SetVec(m, (ys[m]-ys[m-1])/dxR)
+	a.SetBand(m, m, -dxR/3)
+	a.SetBand(m, m-1, -dxR/6)
+	x := mat.NewVecDense(n, nil)
+	err := a.SolveVecTo(x, false, b)
+	if err == nil {
+		cc.cubic.fitWithSecondDerivatives(xs, ys, x.RawVector().Data)
+	}
+	return err
+}
diff --git a/lapack/gonum/dlasy2.go b/lapack/gonum/dlasy2.go
index d670f40d..d95402ca 100644
--- a/lapack/gonum/dlasy2.go
+++ b/lapack/gonum/dlasy2.go
@@ -261,7 +261,7 @@ func (impl Implementation) Dlasy2(tranl, tranr bool, isgn, n1, n2 int, tl []floa
 
 		maxbtmp := math.Max(math.Abs(btmp[0]), math.Abs(btmp[1]))
 		maxbtmp = math.Max(maxbtmp, math.Max(math.Abs(btmp[2]), math.Abs(btmp[3])))
-		scale = 1 / 8 / maxbtmp
+		scale = (1.0 / 8.0) / maxbtmp
 		btmp[0] *= scale
 		btmp[1] *= scale
 		btmp[2] *= scale
diff --git a/lapack/gonum/dormbr.go b/lapack/gonum/dormbr.go
index e1b37373..20f3ca57 100644
--- a/lapack/gonum/dormbr.go
+++ b/lapack/gonum/dormbr.go
@@ -145,7 +145,7 @@ func (impl Implementation) Dormbr(vect lapack.ApplyOrtho, side blas.Side, trans
 				i1 = 1
 				i2 = 0
 			}
-			impl.Dormqr(side, trans, mi, ni, nq-1, a[1*lda:], lda, tau[:nq-1], c[i1*ldc+i2:], ldc, work, lwork)
+			impl.Dormqr(side, trans, mi, ni, nq-1, a[lda:], lda, tau[:nq-1], c[i1*ldc+i2:], ldc, work, lwork)
 		}
 		work[0] = float64(lworkopt)
 		return


Test patch:
diff --git a/graph/formats/rdf/iso_canonical_test.go b/graph/formats/rdf/iso_canonical_test.go
index 580b1e2c..74b94b13 100644
--- a/graph/formats/rdf/iso_canonical_test.go
+++ b/graph/formats/rdf/iso_canonical_test.go
@@ -14,7 +14,6 @@ import (
 	"path/filepath"
 	"reflect"
 	"sort"
-	"strings"
 	"testing"
 	"text/tabwriter"
 	"time"
@@ -287,7 +286,7 @@ func permuteBlanks(s []*Statement, src rand.Source) ([]*Statement, map[string]st
 
 	var blanks []string
 	for t := range terms {
-		if strings.HasPrefix(t, "_:") {
+		if isBlank(t) {
 			blanks = append(blanks, t)
 		}
 	}
@@ -309,7 +308,7 @@ func hashBlanks(s []*Statement, h hash.Hash) ([]*Statement, map[string]string) {
 			e.Object.Value,
 			e.Label.Value,
 		} {
-			if !strings.HasPrefix(t, "_:") {
+			if !isBlank(t) {
 				continue
 			}
 			h.Reset()
@@ -331,7 +330,7 @@ func mangleFirstIL(s []*Statement, h hash.Hash) ([]*Statement, map[string]string
 			e.Object.Value,
 			e.Label.Value,
 		} {
-			if strings.HasPrefix(t, "_:") {
+			if isBlank(t) {
 				continue
 			}
 			h.Reset()
@@ -354,7 +353,7 @@ func mergeFirst2B(s []*Statement) ([]*Statement, map[string]string) {
 			e.Object.Value,
 			e.Label.Value,
 		} {
-			if !strings.HasPrefix(t, "_:") {
+			if !isBlank(t) {
 				continue
 			}
 			terms[t] = t
@@ -444,7 +443,7 @@ func TestLexicalStatements(t *testing.T) {
 func termsFor(hashes map[string][]byte, hash hash.Hash) map[string]string {
 	terms := make(map[string]string)
 	for t, h := range hashes {
-		if strings.HasPrefix(t, "_:") {
+		if isBlank(t) {
 			terms[t] = fmt.Sprintf("_:%0*x", 2*hash.Size(), h)
 		}
 	}
diff --git a/interp/cubic_test.go b/interp/cubic_test.go
index b59df3e8..63b24cb1 100644
--- a/interp/cubic_test.go
+++ b/interp/cubic_test.go
@@ -799,6 +799,48 @@ func TestNaturalCubicFit(t *testing.T) {
 		if math.Abs(got-want[i]) > 1e-14 {
 			t.Errorf("Mismatch in predicted Y value for x = %g: got %v, want %g", testXs[i], got, want[i])
 		}
+		got = nc.PredictDerivative(testXs[i])
+		want := discrDerivPredict(&nc, xs[0], xs[len(xs)-1], testXs[i], 1e-9)
+		if math.Abs(got-want) > 1e-6 {
+			t.Errorf("Mismatch in predicted dY/dX value for x = %g: got %v, want %g", testXs[i], got, want)
+		}
 	}
+}
 
+func TestClampedCubicFit(t *testing.T) {
+	t.Parallel()
+	xs := []float64{-1, 0, 2, 3.5}
+	ys := []float64{2, 0, 2, 1.5}
+	var cc ClampedCubic
+	err := cc.Fit(xs, ys)
+	if err != nil {
+		t.Errorf("Error when fitting ClampedCubic: %v", err)
+	}
+	testXs := []float64{-1, -0.99, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.49, 3.5}
+	// From scipy.interpolate.CubicSpline:
+	want := []float64{
+		2.0,
+		1.999564111111111,
+		1.2021604938271606,
+		0.0,
+		-0.20833333333333343,
+		0.41975308641975284,
+		1.337962962962962,
+		2.0,
+		2.026748971193416,
+		1.7078189300411522,
+		1.5001129711934151,
+		1.4999999999999996,
+	}
+	for i := 0; i < len(testXs); i++ {
+		got := cc.Predict(testXs[i])
+		if math.Abs(got-want[i]) > 1e-14 {
+			t.Errorf("Mismatch in predicted Y value for x = %g: got %v, want %g", testXs[i], got, want[i])
+		}
+		got = cc.PredictDerivative(testXs[i])
+		want := discrDerivPredict(&cc, xs[0], xs[len(xs)-1], testXs[i], 1e-9)
+		if math.Abs(got-want) > 1e-6 {
+			t.Errorf("Mismatch in predicted dY/dX value for x = %g: got %v, want %g", testXs[i], got, want)
+		}
+	}
 }