From 2a0e58f891980f2892568d5f3511143c0179de58 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Sun, 7 Sep 2014 21:40:12 -0700
Subject: [PATCH 01/98] Added README

---
 README.md | 13 +++++++++++++
 1 file changed, 13 insertions(+)
 create mode 100644 README.md

diff --git a/README.md b/README.md
new file mode 100644
index 00000000..f0402485
--- /dev/null
+++ b/README.md
@@ -0,0 +1,13 @@
+# Gonum diff
+
+This is a package for computing derivatives of functions for the Go language.
+
+## Issues
+
+If you find any bugs, feel free to file an issue on the github issue tracker. Discussions on API changes, added features, code review, or similar requests are preferred on the gonum-dev Google Group.
+
+https://groups.google.com/forum/#!forum/gonum-dev
+
+## License
+
+Please see github.com/gonum/license for general license information, contributors, authors, etc on the Gonum suite of packages.

From b121b5dd1fdd3b73e59b55ef14255c7efac56e4f Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Sun, 7 Sep 2014 21:48:43 -0700
Subject: [PATCH 02/98] Initial commit for diff package. Includes
 FiniteDifference functions and several methods.

---
 diff.go      | 104 +++++++++++++++++++++++++++++++++++++++++++++++++++
 diff_test.go |  92 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 196 insertions(+)
 create mode 100644 diff.go
 create mode 100644 diff_test.go

diff --git a/diff.go b/diff.go
new file mode 100644
index 00000000..1b1f878e
--- /dev/null
+++ b/diff.go
@@ -0,0 +1,104 @@
+package diff
+
+import (
+	"math"
+	"sync"
+)
+
+// A Point is a stencil location in a difference method
+type Point struct {
+	Loc   float64
+	Coeff float64
+}
+
+// FDMethod is a specific finite difference method. FDMethod specifies the stencil,
+// that is, the function locations (relative to x) which will be used to estimate
+// the derivative. It also specifies the order of derivative it estimates. Order = 1
+// represents the derivative, Order = 2 represents the curvature, etc.
+type FDMethod struct {
+	Stencil []Point
+	Order   int // The order of the difference method (first derivative, second derivative, etc.)
+}
+
+// FDSettings is the settings structure for the FiniteDifference function
+type FDSettings struct {
+	OriginKnown bool     // Flag that the value at the origin x is known
+	OriginValue float64  // Value at the origin (only used if OriginKnown is true)
+	Step        float64  // step size
+	Concurrent  bool     // Should the function calls be executed concurrently
+	Method      FDMethod // Finite difference method to use
+}
+
+// DefaultFDSettings is a basic set of settings for the FiniteDifference
+// function. Computes a central difference approximation for the first derivative
+// of the function.
+func DefaultFDSettings() *FDSettings {
+	return &FDSettings{
+		Step:   1e-6,
+		Method: Central,
+	}
+}
+
+// FiniteDifference estimates a derivative of the function f at the given location.
+// The order of derivative, sample locations, and other options are specified
+// by settings.
+func FiniteDiffernce(f func(float64) float64, x float64, settings *FDSettings) float64 {
+	var deriv float64
+	method := settings.Method
+	if !settings.Concurrent {
+		for _, pt := range method.Stencil {
+			if settings.OriginKnown && pt.Loc == 0 {
+				deriv += pt.Coeff * settings.OriginValue
+				continue
+			}
+			deriv += pt.Coeff * f(x+settings.Step*pt.Loc)
+		}
+		return deriv / math.Pow(settings.Step, float64(method.Order))
+	}
+
+	wg := &sync.WaitGroup{}
+	mux := &sync.Mutex{}
+	for _, pt := range method.Stencil {
+		if settings.OriginKnown && pt.Loc == 0 {
+			mux.Lock()
+			deriv += pt.Coeff * settings.OriginValue
+			mux.Unlock()
+			continue
+		}
+		wg.Add(1)
+		go func(pt Point) {
+			defer wg.Done()
+			fofx := f(x + settings.Step*pt.Loc)
+			mux.Lock()
+			defer mux.Unlock()
+			deriv += pt.Coeff * fofx
+
+		}(pt)
+	}
+	wg.Wait()
+	return deriv / math.Pow(settings.Step, float64(method.Order))
+}
+
+// Forward represents a first-order forward difference.
+var Forward = FDMethod{
+	Stencil: []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
+	Order:   1,
+}
+
+// Backward represents a first-order backward difference
+var Backward = FDMethod{
+	Stencil: []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
+	Order:   1,
+}
+
+// Central represents a first-order central difference.
+var Central = FDMethod{
+	Stencil: []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
+	Order:   1,
+}
+
+// Central2nd represents a secord-order central difference.
+var Central2nd = FDMethod{
+	Stencil: []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
+	Order:   2,
+}
diff --git a/diff_test.go b/diff_test.go
new file mode 100644
index 00000000..ab6ff6d5
--- /dev/null
+++ b/diff_test.go
@@ -0,0 +1,92 @@
+package diff
+
+import (
+	"math"
+	"testing"
+)
+
+var xSquared = func(x float64) float64 { return x * x }
+
+type testPoint struct {
+	f   func(float64) float64
+	loc float64
+	ans float64
+}
+
+var testsFirst = []testPoint{
+	{
+		f:   xSquared,
+		loc: 0,
+		ans: 0,
+	},
+	{
+		f:   xSquared,
+		loc: 5,
+		ans: 10,
+	},
+	{
+		f:   xSquared,
+		loc: 2,
+		ans: 4,
+	},
+	{
+		f:   xSquared,
+		loc: -5,
+		ans: -10,
+	},
+}
+
+var testsSecond = []testPoint{
+	{
+		f:   xSquared,
+		loc: 0,
+		ans: 2,
+	},
+	{
+		f:   xSquared,
+		loc: 5,
+		ans: 2,
+	},
+	{
+		f:   xSquared,
+		loc: 2,
+		ans: 2,
+	},
+	{
+		f:   xSquared,
+		loc: -5,
+		ans: 2,
+	},
+}
+
+func testFirstOrder(t *testing.T, method FDMethod, tol float64, tests []testPoint) {
+	for _, test := range tests {
+		settings := DefaultFDSettings()
+		settings.Method = method
+		ans := FiniteDiffernce(test.f, test.loc, settings)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
+		}
+		settings.Concurrent = true
+		ans = FiniteDiffernce(test.f, test.loc, settings)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
+		}
+	}
+}
+
+func TestForward(t *testing.T) {
+	testFirstOrder(t, Forward, 1e-4, testsFirst)
+}
+
+func TestBackward(t *testing.T) {
+	testFirstOrder(t, Backward, 1e-4, testsFirst)
+}
+
+func TestCentral(t *testing.T) {
+	testFirstOrder(t, Central, 1e-6, testsFirst)
+}
+
+func TestCentralSecond(t *testing.T) {
+	testFirstOrder(t, Central2nd, 2e-3, testsSecond)
+}

From 65c7b50a89a18d0d51fa012fd9686545c14a62d0 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Sun, 7 Sep 2014 22:00:19 -0700
Subject: [PATCH 03/98] Fixed spelling of Difference

---
 diff.go      | 2 +-
 diff_test.go | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/diff.go b/diff.go
index 1b1f878e..64bda0ad 100644
--- a/diff.go
+++ b/diff.go
@@ -42,7 +42,7 @@ func DefaultFDSettings() *FDSettings {
 // FiniteDifference estimates a derivative of the function f at the given location.
 // The order of derivative, sample locations, and other options are specified
 // by settings.
-func FiniteDiffernce(f func(float64) float64, x float64, settings *FDSettings) float64 {
+func FiniteDifference(f func(float64) float64, x float64, settings *FDSettings) float64 {
 	var deriv float64
 	method := settings.Method
 	if !settings.Concurrent {
diff --git a/diff_test.go b/diff_test.go
index ab6ff6d5..81f90482 100644
--- a/diff_test.go
+++ b/diff_test.go
@@ -63,12 +63,12 @@ func testFirstOrder(t *testing.T, method FDMethod, tol float64, tests []testPoin
 	for _, test := range tests {
 		settings := DefaultFDSettings()
 		settings.Method = method
-		ans := FiniteDiffernce(test.f, test.loc, settings)
+		ans := FiniteDifference(test.f, test.loc, settings)
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
 		}
 		settings.Concurrent = true
-		ans = FiniteDiffernce(test.f, test.loc, settings)
+		ans = FiniteDifference(test.f, test.loc, settings)
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
 		}

From 97585958427a2e45dec1de55e9f8f7e71879f944 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Mon, 8 Sep 2014 21:57:11 -0700
Subject: [PATCH 04/98] Moved finite difference routines to a subpackage and
 removed stuttering names

---
 fd/diff.go      | 104 ++++++++++++++++++++++++++++++++++++++++++++++++
 fd/diff_test.go |  92 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 196 insertions(+)
 create mode 100644 fd/diff.go
 create mode 100644 fd/diff_test.go

diff --git a/fd/diff.go b/fd/diff.go
new file mode 100644
index 00000000..4d749692
--- /dev/null
+++ b/fd/diff.go
@@ -0,0 +1,104 @@
+package diff
+
+import (
+	"math"
+	"sync"
+)
+
+// A Point is a stencil location in a difference method
+type Point struct {
+	Loc   float64
+	Coeff float64
+}
+
+// Method is a specific finite difference method. Method specifies the stencil,
+// that is, the function locations (relative to x) which will be used to estimate
+// the derivative. It also specifies the order of derivative it estimates. Order = 1
+// represents the derivative, Order = 2 represents the curvature, etc.
+type Method struct {
+	Stencil []Point
+	Order   int // The order of the difference method (first derivative, second derivative, etc.)
+}
+
+// Settings is the settings structure for computing finite differences.
+type Settings struct {
+	OriginKnown bool    // Flag that the value at the origin x is known
+	OriginValue float64 // Value at the origin (only used if OriginKnown is true)
+	Step        float64 // step size
+	Concurrent  bool    // Should the function calls be executed concurrently
+	Method      Method  // Finite difference method to use
+}
+
+// DefaultSettings is a basic set of settings for computing finite differences.
+// Computes a central difference approximation for the first derivative
+// of the function.
+func DefaultSettings() *Settings {
+	return &Settings{
+		Step:   1e-6,
+		Method: Central,
+	}
+}
+
+// Derivative estimates the derivative of the function f at the given location.
+// The order of derivative, sample locations, and other options are specified
+// by settings.
+func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
+	var deriv float64
+	method := settings.Method
+	if !settings.Concurrent {
+		for _, pt := range method.Stencil {
+			if settings.OriginKnown && pt.Loc == 0 {
+				deriv += pt.Coeff * settings.OriginValue
+				continue
+			}
+			deriv += pt.Coeff * f(x+settings.Step*pt.Loc)
+		}
+		return deriv / math.Pow(settings.Step, float64(method.Order))
+	}
+
+	wg := &sync.WaitGroup{}
+	mux := &sync.Mutex{}
+	for _, pt := range method.Stencil {
+		if settings.OriginKnown && pt.Loc == 0 {
+			mux.Lock()
+			deriv += pt.Coeff * settings.OriginValue
+			mux.Unlock()
+			continue
+		}
+		wg.Add(1)
+		go func(pt Point) {
+			defer wg.Done()
+			fofx := f(x + settings.Step*pt.Loc)
+			mux.Lock()
+			defer mux.Unlock()
+			deriv += pt.Coeff * fofx
+
+		}(pt)
+	}
+	wg.Wait()
+	return deriv / math.Pow(settings.Step, float64(method.Order))
+}
+
+// Forward represents a first-order forward difference.
+var Forward = Method{
+	Stencil: []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
+	Order:   1,
+}
+
+// Backward represents a first-order backward difference
+var Backward = Method{
+	Stencil: []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
+	Order:   1,
+}
+
+// Central represents a first-order central difference.
+var Central = Method{
+	Stencil: []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
+	Order:   1,
+}
+
+// Central2nd represents a secord-order central difference.
+var Central2nd = Method{
+	Stencil: []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
+	Order:   2,
+}
diff --git a/fd/diff_test.go b/fd/diff_test.go
new file mode 100644
index 00000000..2cd98ee1
--- /dev/null
+++ b/fd/diff_test.go
@@ -0,0 +1,92 @@
+package diff
+
+import (
+	"math"
+	"testing"
+)
+
+var xSquared = func(x float64) float64 { return x * x }
+
+type testPoint struct {
+	f   func(float64) float64
+	loc float64
+	ans float64
+}
+
+var testsFirst = []testPoint{
+	{
+		f:   xSquared,
+		loc: 0,
+		ans: 0,
+	},
+	{
+		f:   xSquared,
+		loc: 5,
+		ans: 10,
+	},
+	{
+		f:   xSquared,
+		loc: 2,
+		ans: 4,
+	},
+	{
+		f:   xSquared,
+		loc: -5,
+		ans: -10,
+	},
+}
+
+var testsSecond = []testPoint{
+	{
+		f:   xSquared,
+		loc: 0,
+		ans: 2,
+	},
+	{
+		f:   xSquared,
+		loc: 5,
+		ans: 2,
+	},
+	{
+		f:   xSquared,
+		loc: 2,
+		ans: 2,
+	},
+	{
+		f:   xSquared,
+		loc: -5,
+		ans: 2,
+	},
+}
+
+func testFirstOrder(t *testing.T, method Method, tol float64, tests []testPoint) {
+	for _, test := range tests {
+		settings := DefaultSettings()
+		settings.Method = method
+		ans := Derivative(test.f, test.loc, settings)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
+		}
+		settings.Concurrent = true
+		ans = Derivative(test.f, test.loc, settings)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
+		}
+	}
+}
+
+func TestForward(t *testing.T) {
+	testFirstOrder(t, Forward, 1e-4, testsFirst)
+}
+
+func TestBackward(t *testing.T) {
+	testFirstOrder(t, Backward, 1e-4, testsFirst)
+}
+
+func TestCentral(t *testing.T) {
+	testFirstOrder(t, Central, 1e-6, testsFirst)
+}
+
+func TestCentralSecond(t *testing.T) {
+	testFirstOrder(t, Central2nd, 2e-3, testsSecond)
+}

From 19aaf536f1a525d40909ea9de0698ed8dd5f5e9b Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Tue, 9 Sep 2014 07:22:07 -0700
Subject: [PATCH 05/98] Removed extra file copies

---
 diff.go      | 104 ---------------------------------------------------
 diff_test.go |  92 ---------------------------------------------
 2 files changed, 196 deletions(-)
 delete mode 100644 diff.go
 delete mode 100644 diff_test.go

diff --git a/diff.go b/diff.go
deleted file mode 100644
index 64bda0ad..00000000
--- a/diff.go
+++ /dev/null
@@ -1,104 +0,0 @@
-package diff
-
-import (
-	"math"
-	"sync"
-)
-
-// A Point is a stencil location in a difference method
-type Point struct {
-	Loc   float64
-	Coeff float64
-}
-
-// FDMethod is a specific finite difference method. FDMethod specifies the stencil,
-// that is, the function locations (relative to x) which will be used to estimate
-// the derivative. It also specifies the order of derivative it estimates. Order = 1
-// represents the derivative, Order = 2 represents the curvature, etc.
-type FDMethod struct {
-	Stencil []Point
-	Order   int // The order of the difference method (first derivative, second derivative, etc.)
-}
-
-// FDSettings is the settings structure for the FiniteDifference function
-type FDSettings struct {
-	OriginKnown bool     // Flag that the value at the origin x is known
-	OriginValue float64  // Value at the origin (only used if OriginKnown is true)
-	Step        float64  // step size
-	Concurrent  bool     // Should the function calls be executed concurrently
-	Method      FDMethod // Finite difference method to use
-}
-
-// DefaultFDSettings is a basic set of settings for the FiniteDifference
-// function. Computes a central difference approximation for the first derivative
-// of the function.
-func DefaultFDSettings() *FDSettings {
-	return &FDSettings{
-		Step:   1e-6,
-		Method: Central,
-	}
-}
-
-// FiniteDifference estimates a derivative of the function f at the given location.
-// The order of derivative, sample locations, and other options are specified
-// by settings.
-func FiniteDifference(f func(float64) float64, x float64, settings *FDSettings) float64 {
-	var deriv float64
-	method := settings.Method
-	if !settings.Concurrent {
-		for _, pt := range method.Stencil {
-			if settings.OriginKnown && pt.Loc == 0 {
-				deriv += pt.Coeff * settings.OriginValue
-				continue
-			}
-			deriv += pt.Coeff * f(x+settings.Step*pt.Loc)
-		}
-		return deriv / math.Pow(settings.Step, float64(method.Order))
-	}
-
-	wg := &sync.WaitGroup{}
-	mux := &sync.Mutex{}
-	for _, pt := range method.Stencil {
-		if settings.OriginKnown && pt.Loc == 0 {
-			mux.Lock()
-			deriv += pt.Coeff * settings.OriginValue
-			mux.Unlock()
-			continue
-		}
-		wg.Add(1)
-		go func(pt Point) {
-			defer wg.Done()
-			fofx := f(x + settings.Step*pt.Loc)
-			mux.Lock()
-			defer mux.Unlock()
-			deriv += pt.Coeff * fofx
-
-		}(pt)
-	}
-	wg.Wait()
-	return deriv / math.Pow(settings.Step, float64(method.Order))
-}
-
-// Forward represents a first-order forward difference.
-var Forward = FDMethod{
-	Stencil: []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
-	Order:   1,
-}
-
-// Backward represents a first-order backward difference
-var Backward = FDMethod{
-	Stencil: []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
-	Order:   1,
-}
-
-// Central represents a first-order central difference.
-var Central = FDMethod{
-	Stencil: []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
-	Order:   1,
-}
-
-// Central2nd represents a secord-order central difference.
-var Central2nd = FDMethod{
-	Stencil: []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
-	Order:   2,
-}
diff --git a/diff_test.go b/diff_test.go
deleted file mode 100644
index 81f90482..00000000
--- a/diff_test.go
+++ /dev/null
@@ -1,92 +0,0 @@
-package diff
-
-import (
-	"math"
-	"testing"
-)
-
-var xSquared = func(x float64) float64 { return x * x }
-
-type testPoint struct {
-	f   func(float64) float64
-	loc float64
-	ans float64
-}
-
-var testsFirst = []testPoint{
-	{
-		f:   xSquared,
-		loc: 0,
-		ans: 0,
-	},
-	{
-		f:   xSquared,
-		loc: 5,
-		ans: 10,
-	},
-	{
-		f:   xSquared,
-		loc: 2,
-		ans: 4,
-	},
-	{
-		f:   xSquared,
-		loc: -5,
-		ans: -10,
-	},
-}
-
-var testsSecond = []testPoint{
-	{
-		f:   xSquared,
-		loc: 0,
-		ans: 2,
-	},
-	{
-		f:   xSquared,
-		loc: 5,
-		ans: 2,
-	},
-	{
-		f:   xSquared,
-		loc: 2,
-		ans: 2,
-	},
-	{
-		f:   xSquared,
-		loc: -5,
-		ans: 2,
-	},
-}
-
-func testFirstOrder(t *testing.T, method FDMethod, tol float64, tests []testPoint) {
-	for _, test := range tests {
-		settings := DefaultFDSettings()
-		settings.Method = method
-		ans := FiniteDifference(test.f, test.loc, settings)
-		if math.Abs(test.ans-ans) > tol {
-			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
-		}
-		settings.Concurrent = true
-		ans = FiniteDifference(test.f, test.loc, settings)
-		if math.Abs(test.ans-ans) > tol {
-			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
-		}
-	}
-}
-
-func TestForward(t *testing.T) {
-	testFirstOrder(t, Forward, 1e-4, testsFirst)
-}
-
-func TestBackward(t *testing.T) {
-	testFirstOrder(t, Backward, 1e-4, testsFirst)
-}
-
-func TestCentral(t *testing.T) {
-	testFirstOrder(t, Central, 1e-6, testsFirst)
-}
-
-func TestCentralSecond(t *testing.T) {
-	testFirstOrder(t, Central2nd, 2e-3, testsSecond)
-}

From d55939055a476d2c4e7924e6050ba11b793f86f9 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Tue, 14 Oct 2014 16:46:15 -0700
Subject: [PATCH 06/98] Added Gradient method and tests.

Gradient computes finite differences of vector-valued function using finite difference
---
 fd/diff.go          | 106 +++++++++++++++++++++++++++++++++++-
 fd/gradient_test.go | 130 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 234 insertions(+), 2 deletions(-)
 create mode 100644 fd/gradient_test.go

diff --git a/fd/diff.go b/fd/diff.go
index 4d749692..d4d216f8 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -2,7 +2,10 @@ package diff
 
 import (
 	"math"
+	"runtime"
 	"sync"
+
+	"github.com/gonum/floats"
 )
 
 // A Point is a stencil location in a difference method
@@ -26,6 +29,7 @@ type Settings struct {
 	OriginValue float64 // Value at the origin (only used if OriginKnown is true)
 	Step        float64 // step size
 	Concurrent  bool    // Should the function calls be executed concurrently
+	Workers     int     // Maximum number of concurrent executions when evaluating concurrently
 	Method      Method  // Finite difference method to use
 }
 
@@ -34,8 +38,9 @@ type Settings struct {
 // of the function.
 func DefaultSettings() *Settings {
 	return &Settings{
-		Step:   1e-6,
-		Method: Central,
+		Step:    1e-6,
+		Method:  Central,
+		Workers: runtime.GOMAXPROCS(0),
 	}
 }
 
@@ -79,6 +84,103 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 	return deriv / math.Pow(settings.Step, float64(method.Order))
 }
 
+// Gradient estimates the derivative of a vector-valued function f at the location
+// x. The resulting estimate is stored in-place into deriv. The order of derivative,
+// sample locations, and other options are specified by settings. Gradient panics
+// if len(deriv) != len(x).
+func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradient []float64) {
+	if len(gradient) != len(x) {
+		panic("fd: location and gradient length mismatch")
+	}
+	if settings == nil {
+		settings = DefaultSettings()
+	}
+	if !settings.Concurrent {
+		xcopy := make([]float64, len(x)) // So that x is not modified during the call
+		copy(xcopy, x)
+		for i := range xcopy {
+			var deriv float64
+			for _, pt := range settings.Method.Stencil {
+				if settings.OriginKnown && pt.Loc == 0 {
+					deriv += pt.Coeff * settings.OriginValue
+					continue
+				}
+				xcopy[i] += pt.Loc * settings.Step
+				deriv += pt.Coeff * f(xcopy)
+				xcopy[i] -= pt.Loc * settings.Step
+			}
+			gradient[i] = deriv / math.Pow(settings.Step, float64(settings.Method.Order))
+		}
+		return
+	}
+	quit := make(chan struct{})
+	sendChan := make(chan fdrun)
+	ansChan := make(chan fdrun)
+
+	nWorkers := settings.Workers
+	expect := len(settings.Method.Stencil) * len(x)
+	if nWorkers > expect {
+		nWorkers = expect
+	}
+
+	// Launch workers. Workers receive an index and a step, and compute the answer
+	for i := 0; i < settings.Workers; i++ {
+		go func() {
+			xcopy := make([]float64, len(x))
+			copy(xcopy, x)
+			for {
+				select {
+				case <-quit:
+					return
+				case run := <-sendChan:
+					xcopy[run.idx] += run.pt.Loc * settings.Step
+					run.result = f(xcopy)
+					xcopy[run.idx] -= run.pt.Loc * settings.Step
+					ansChan <- run
+				}
+			}
+		}()
+	}
+
+	// Launch the distributor. Distributor sends the cases to be computed
+	go func() {
+		for i := range x {
+			for _, pt := range settings.Method.Stencil {
+				if settings.OriginKnown && pt.Loc == 0 {
+					// Answer already known. Send the answer on the answer channel
+					ansChan <- fdrun{
+						idx:    i,
+						pt:     pt,
+						result: settings.OriginValue,
+					}
+					continue
+				}
+				// Answer not known, send the answer to be computed
+				sendChan <- fdrun{
+					idx: i,
+					pt:  pt,
+				}
+			}
+		}
+	}()
+	for i := range gradient {
+		gradient[i] = 0
+	}
+	// Read in all of the results
+	for i := 0; i < expect; i++ {
+		run := <-ansChan
+		gradient[run.idx] += run.pt.Coeff * run.result
+	}
+	floats.Scale(1/math.Pow(settings.Step, float64(settings.Method.Order)), gradient)
+	close(quit)
+}
+
+type fdrun struct {
+	idx    int
+	pt     Point
+	result float64
+}
+
 // Forward represents a first-order forward difference.
 var Forward = Method{
 	Stencil: []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
new file mode 100644
index 00000000..81afdbde
--- /dev/null
+++ b/fd/gradient_test.go
@@ -0,0 +1,130 @@
+package diff
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/floats"
+)
+
+type Rosenbrock struct {
+	nDim int
+}
+
+func (r Rosenbrock) F(x []float64) (sum float64) {
+	deriv := make([]float64, len(x))
+	return r.FDf(x, deriv)
+}
+
+func (r Rosenbrock) FDf(x []float64, deriv []float64) (sum float64) {
+	sum = 0
+
+	for i := range deriv {
+		deriv[i] = 0
+	}
+
+	for i := 0; i < len(x)-1; i++ {
+		sum += math.Pow(1-x[i], 2) + 100*math.Pow(x[i+1]-math.Pow(x[i], 2), 2)
+	}
+	for i := 0; i < len(x)-1; i++ {
+		deriv[i] += -1 * 2 * (1 - x[i])
+		deriv[i] += 2 * 100 * (x[i+1] - math.Pow(x[i], 2)) * (-2 * x[i])
+	}
+	for i := 1; i < len(x); i++ {
+		deriv[i] += 2 * 100 * (x[i] - math.Pow(x[i-1], 2))
+	}
+
+	return sum
+}
+
+func TestGradient(t *testing.T) {
+	for i, test := range []struct {
+		nDim   int
+		tol    float64
+		method Method
+	}{
+		{
+			nDim:   2,
+			tol:    1e-4,
+			method: Forward,
+		},
+		{
+			nDim:   2,
+			tol:    1e-6,
+			method: Central,
+		},
+		{
+			nDim:   40,
+			tol:    1e-4,
+			method: Forward,
+		},
+		{
+			nDim:   40,
+			tol:    1e-6,
+			method: Central,
+		},
+	} {
+		x := make([]float64, test.nDim)
+		for i := range x {
+			x[i] = rand.Float64()
+		}
+		xcopy := make([]float64, len(x))
+		copy(xcopy, x)
+
+		r := Rosenbrock{len(x)}
+		trueGradient := make([]float64, len(x))
+		r.FDf(x, trueGradient)
+
+		settings := DefaultSettings()
+		settings.Method = test.method
+		gradient := make([]float64, len(x))
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+
+		Gradient(r.F, x, settings, gradient)
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %v: x modified during call to gradient in serial")
+		}
+
+		// Try with known value
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		settings.OriginKnown = true
+		settings.OriginValue = r.F(x)
+		Gradient(r.F, x, settings, gradient)
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with known origin in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+
+		// Concurrently
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		settings.Concurrent = true
+		settings.OriginKnown = false
+		Gradient(r.F, x, settings, gradient)
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %v: x modified during call to gradient in parallel")
+		}
+
+		// Concurrently with origin known
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		settings.OriginKnown = true
+		Gradient(r.F, x, settings, gradient)
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with known origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+
+	}
+}

From fb46ec0bf4336507cb360d331ef5b7e90a213664 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Tue, 14 Oct 2014 16:54:16 -0700
Subject: [PATCH 07/98] Made the gradient channels buffered to avoid contention

---
 fd/diff.go | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index d4d216f8..9b62d17d 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -113,9 +113,6 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 		}
 		return
 	}
-	quit := make(chan struct{})
-	sendChan := make(chan fdrun)
-	ansChan := make(chan fdrun)
 
 	nWorkers := settings.Workers
 	expect := len(settings.Method.Stencil) * len(x)
@@ -123,6 +120,10 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 		nWorkers = expect
 	}
 
+	quit := make(chan struct{})
+	sendChan := make(chan fdrun, expect)
+	ansChan := make(chan fdrun, expect)
+
 	// Launch workers. Workers receive an index and a step, and compute the answer
 	for i := 0; i < settings.Workers; i++ {
 		go func() {

From e2ba9707219382fdd975101d5023f4e9384eac2b Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Wed, 15 Oct 2014 00:09:37 -0700
Subject: [PATCH 08/98] Fixed package name and added file headers

---
 fd/diff.go          | 6 +++++-
 fd/diff_test.go     | 6 +++++-
 fd/gradient_test.go | 6 +++++-
 3 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 9b62d17d..12c00fb5 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -1,4 +1,8 @@
-package diff
+// Copyright ©2014 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 fd
 
 import (
 	"math"
diff --git a/fd/diff_test.go b/fd/diff_test.go
index 2cd98ee1..04b13a8b 100644
--- a/fd/diff_test.go
+++ b/fd/diff_test.go
@@ -1,4 +1,8 @@
-package diff
+// Copyright ©2014 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 fd
 
 import (
 	"math"
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 81afdbde..08e87266 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -1,4 +1,8 @@
-package diff
+// Copyright ©2014 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 fd
 
 import (
 	"math"

From 98a1d1ac6b2c00916664264b70c981a00f7772b3 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Wed, 15 Oct 2014 11:11:22 -0700
Subject: [PATCH 09/98] Added channels as receiver arguments and used
 assignment rather than subtraction in resetting xcopy

---
 fd/diff.go | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 12c00fb5..f05c1b88 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -111,7 +111,7 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 				}
 				xcopy[i] += pt.Loc * settings.Step
 				deriv += pt.Coeff * f(xcopy)
-				xcopy[i] -= pt.Loc * settings.Step
+				xcopy[i] = x[i]
 			}
 			gradient[i] = deriv / math.Pow(settings.Step, float64(settings.Method.Order))
 		}
@@ -125,12 +125,13 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 	}
 
 	quit := make(chan struct{})
+	defer close(quit)
 	sendChan := make(chan fdrun, expect)
 	ansChan := make(chan fdrun, expect)
 
 	// Launch workers. Workers receive an index and a step, and compute the answer
 	for i := 0; i < settings.Workers; i++ {
-		go func() {
+		go func(sendChan <-chan fdrun, ansChan chan<- fdrun, quit <-chan struct{}) {
 			xcopy := make([]float64, len(x))
 			copy(xcopy, x)
 			for {
@@ -140,15 +141,15 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 				case run := <-sendChan:
 					xcopy[run.idx] += run.pt.Loc * settings.Step
 					run.result = f(xcopy)
-					xcopy[run.idx] -= run.pt.Loc * settings.Step
+					xcopy[run.idx] = x[run.idx]
 					ansChan <- run
 				}
 			}
-		}()
+		}(sendChan, ansChan, quit)
 	}
 
 	// Launch the distributor. Distributor sends the cases to be computed
-	go func() {
+	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
 		for i := range x {
 			for _, pt := range settings.Method.Stencil {
 				if settings.OriginKnown && pt.Loc == 0 {
@@ -167,7 +168,8 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 				}
 			}
 		}
-	}()
+	}(sendChan, ansChan)
+
 	for i := range gradient {
 		gradient[i] = 0
 	}
@@ -177,7 +179,6 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 		gradient[run.idx] += run.pt.Coeff * run.result
 	}
 	floats.Scale(1/math.Pow(settings.Step, float64(settings.Method.Order)), gradient)
-	close(quit)
 }
 
 type fdrun struct {

From 4fed8745268f77a7217cd57090c6c51b6edebcdb Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Wed, 15 Oct 2014 14:25:45 -0700
Subject: [PATCH 10/98] Made tests more robust to concurrent environment.

---
 fd/diff.go          |   8 ++--
 fd/diff_test.go     | 105 +++++++++++++++++++++++++++++---------------
 fd/gradient_test.go |   6 +--
 3 files changed, 77 insertions(+), 42 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index f05c1b88..097940cb 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -12,7 +12,7 @@ import (
 	"github.com/gonum/floats"
 )
 
-// A Point is a stencil location in a difference method
+// A Point is a stencil location in a difference method.
 type Point struct {
 	Loc   float64
 	Coeff float64
@@ -52,6 +52,9 @@ func DefaultSettings() *Settings {
 // The order of derivative, sample locations, and other options are specified
 // by settings.
 func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
+	if settings == nil {
+		settings = DefaultSettings()
+	}
 	var deriv float64
 	method := settings.Method
 	if !settings.Concurrent {
@@ -81,7 +84,6 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 			mux.Lock()
 			defer mux.Unlock()
 			deriv += pt.Coeff * fofx
-
 		}(pt)
 	}
 	wg.Wait()
@@ -89,7 +91,7 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 }
 
 // Gradient estimates the derivative of a vector-valued function f at the location
-// x. The resulting estimate is stored in-place into deriv. The order of derivative,
+// x. The resulting estimate is stored in-place into gradient. The order of derivative,
 // sample locations, and other options are specified by settings. Gradient panics
 // if len(deriv) != len(x).
 func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradient []float64) {
diff --git a/fd/diff_test.go b/fd/diff_test.go
index 04b13a8b..00848394 100644
--- a/fd/diff_test.go
+++ b/fd/diff_test.go
@@ -12,85 +12,120 @@ import (
 var xSquared = func(x float64) float64 { return x * x }
 
 type testPoint struct {
-	f   func(float64) float64
-	loc float64
-	ans float64
+	f    func(float64) float64
+	loc  float64
+	fofx float64
+	ans  float64
+	step float64
 }
 
 var testsFirst = []testPoint{
 	{
-		f:   xSquared,
-		loc: 0,
-		ans: 0,
+		f:    xSquared,
+		loc:  0,
+		fofx: 0,
+		ans:  0,
+		step: 1e-6,
 	},
 	{
-		f:   xSquared,
-		loc: 5,
-		ans: 10,
+		f:    xSquared,
+		loc:  5,
+		fofx: 25,
+		ans:  10,
+		step: 1e-6,
 	},
 	{
-		f:   xSquared,
-		loc: 2,
-		ans: 4,
+		f:    xSquared,
+		loc:  2,
+		fofx: 4,
+		ans:  4,
+		step: 1e-6,
 	},
 	{
-		f:   xSquared,
-		loc: -5,
-		ans: -10,
+		f:    xSquared,
+		loc:  -5,
+		fofx: 25,
+		ans:  -10,
+		step: 1e-6,
 	},
 }
 
 var testsSecond = []testPoint{
 	{
-		f:   xSquared,
-		loc: 0,
-		ans: 2,
+		f:    xSquared,
+		loc:  0,
+		fofx: 0,
+		ans:  2,
+		step: 1e-3,
 	},
 	{
-		f:   xSquared,
-		loc: 5,
-		ans: 2,
+		f:    xSquared,
+		loc:  5,
+		fofx: 25,
+		ans:  2,
+		step: 1e-3,
 	},
 	{
-		f:   xSquared,
-		loc: 2,
-		ans: 2,
+		f:    xSquared,
+		loc:  2,
+		fofx: 4,
+		ans:  2,
+		step: 1e-3,
 	},
 	{
-		f:   xSquared,
-		loc: -5,
-		ans: 2,
+		f:    xSquared,
+		loc:  -5,
+		fofx: 25,
+		ans:  2,
+		step: 1e-3,
 	},
 }
 
-func testFirstOrder(t *testing.T, method Method, tol float64, tests []testPoint) {
-	for _, test := range tests {
+func testDerivative(t *testing.T, method Method, tol float64, tests []testPoint) {
+	for i, test := range tests {
+
 		settings := DefaultSettings()
 		settings.Method = method
+		settings.Step = 1e-4
 		ans := Derivative(test.f, test.loc, settings)
 		if math.Abs(test.ans-ans) > tol {
-			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
+			t.Errorf("Case %v: ans mismatch serial: expected %v, found %v", i, test.ans, ans)
 		}
+
+		settings.OriginKnown = true
+		settings.OriginValue = test.fofx
+		ans = Derivative(test.f, test.loc, settings)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch serial origin known: expected %v, found %v", i, test.ans, ans)
+		}
+
+		settings.OriginKnown = false
 		settings.Concurrent = true
 		ans = Derivative(test.f, test.loc, settings)
 		if math.Abs(test.ans-ans) > tol {
-			t.Errorf("ans mismatch: expected %v, found %v", test.ans, ans)
+			t.Errorf("Case %v: ans mismatch concurrent: expected %v, found %v", i, test.ans, ans)
+		}
+
+		settings.OriginKnown = true
+		ans = Derivative(test.f, test.loc, settings)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch concurrent: expected %v, found %v", i, test.ans, ans)
 		}
 	}
 }
 
 func TestForward(t *testing.T) {
-	testFirstOrder(t, Forward, 1e-4, testsFirst)
+	testDerivative(t, Forward, 2e-4, testsFirst)
 }
 
 func TestBackward(t *testing.T) {
-	testFirstOrder(t, Backward, 1e-4, testsFirst)
+	testDerivative(t, Backward, 2e-4, testsFirst)
 }
 
 func TestCentral(t *testing.T) {
-	testFirstOrder(t, Central, 1e-6, testsFirst)
+	testDerivative(t, Central, 1e-6, testsFirst)
 }
 
 func TestCentralSecond(t *testing.T) {
-	testFirstOrder(t, Central2nd, 2e-3, testsSecond)
+	testDerivative(t, Central2nd, 1e-3, testsSecond)
 }
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 08e87266..17b39945 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -22,8 +22,6 @@ func (r Rosenbrock) F(x []float64) (sum float64) {
 }
 
 func (r Rosenbrock) FDf(x []float64, deriv []float64) (sum float64) {
-	sum = 0
-
 	for i := range deriv {
 		deriv[i] = 0
 	}
@@ -50,7 +48,7 @@ func TestGradient(t *testing.T) {
 	}{
 		{
 			nDim:   2,
-			tol:    1e-4,
+			tol:    2e-4,
 			method: Forward,
 		},
 		{
@@ -60,7 +58,7 @@ func TestGradient(t *testing.T) {
 		},
 		{
 			nDim:   40,
-			tol:    1e-4,
+			tol:    2e-4,
 			method: Forward,
 		},
 		{

From cda7abffbc2d17e7ecf1112fcc3f942a706979d0 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Wed, 15 Oct 2014 16:43:05 -0700
Subject: [PATCH 11/98] Fixed fmt statements with missing arguments

---
 fd/gradient_test.go | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 17b39945..8b432b67 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -90,7 +90,7 @@ func TestGradient(t *testing.T) {
 			t.Errorf("Case %v: gradient mismatch in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %v: x modified during call to gradient in serial")
+			t.Errorf("Case %v: x modified during call to gradient in serial", i)
 		}
 
 		// Try with known value
@@ -115,7 +115,7 @@ func TestGradient(t *testing.T) {
 			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %v: x modified during call to gradient in parallel")
+			t.Errorf("Case %v: x modified during call to gradient in parallel", i)
 		}
 
 		// Concurrently with origin known

From 2df99d0b8b1087af5f8032721f2ad0bbadce0b23 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Wed, 15 Oct 2014 16:56:41 -0700
Subject: [PATCH 12/98] fixed Gradient documentation

---
 fd/diff.go | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fd/diff.go b/fd/diff.go
index 097940cb..714705a1 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -90,7 +90,7 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 	return deriv / math.Pow(settings.Step, float64(method.Order))
 }
 
-// Gradient estimates the derivative of a vector-valued function f at the location
+// Gradient estimates the derivative of a multivariate function f at the location
 // x. The resulting estimate is stored in-place into gradient. The order of derivative,
 // sample locations, and other options are specified by settings. Gradient panics
 // if len(deriv) != len(x).

From f136cf3705acfce9f1ebe98fb5202dee5115b7fc Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Wed, 15 Oct 2014 17:43:42 -0700
Subject: [PATCH 13/98] Moved default step to be a function of the method

---
 fd/diff.go      | 41 ++++++++++++++++++++++++++++-------------
 fd/diff_test.go | 10 ----------
 2 files changed, 28 insertions(+), 23 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 714705a1..ee0e3ad8 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -24,7 +24,8 @@ type Point struct {
 // represents the derivative, Order = 2 represents the curvature, etc.
 type Method struct {
 	Stencil []Point
-	Order   int // The order of the difference method (first derivative, second derivative, etc.)
+	Order   int     // The order of the difference method (first derivative, second derivative, etc.)
+	Step    float64 // Default step size for the method.
 }
 
 // Settings is the settings structure for computing finite differences.
@@ -42,7 +43,6 @@ type Settings struct {
 // of the function.
 func DefaultSettings() *Settings {
 	return &Settings{
-		Step:    1e-6,
 		Method:  Central,
 		Workers: runtime.GOMAXPROCS(0),
 	}
@@ -50,11 +50,16 @@ func DefaultSettings() *Settings {
 
 // Derivative estimates the derivative of the function f at the given location.
 // The order of derivative, sample locations, and other options are specified
-// by settings.
+// by settings. If the step is zero, then the step size of the method will
+// be used.
 func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
 	if settings == nil {
 		settings = DefaultSettings()
 	}
+	step := settings.Step
+	if settings.Step == 0 {
+		step = settings.Method.Step
+	}
 	var deriv float64
 	method := settings.Method
 	if !settings.Concurrent {
@@ -63,9 +68,9 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 				deriv += pt.Coeff * settings.OriginValue
 				continue
 			}
-			deriv += pt.Coeff * f(x+settings.Step*pt.Loc)
+			deriv += pt.Coeff * f(x+step*pt.Loc)
 		}
-		return deriv / math.Pow(settings.Step, float64(method.Order))
+		return deriv / math.Pow(step, float64(method.Order))
 	}
 
 	wg := &sync.WaitGroup{}
@@ -80,20 +85,22 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 		wg.Add(1)
 		go func(pt Point) {
 			defer wg.Done()
-			fofx := f(x + settings.Step*pt.Loc)
+			fofx := f(x + step*pt.Loc)
 			mux.Lock()
 			defer mux.Unlock()
 			deriv += pt.Coeff * fofx
 		}(pt)
 	}
 	wg.Wait()
-	return deriv / math.Pow(settings.Step, float64(method.Order))
+	return deriv / math.Pow(step, float64(method.Order))
 }
 
 // Gradient estimates the derivative of a multivariate function f at the location
 // x. The resulting estimate is stored in-place into gradient. The order of derivative,
-// sample locations, and other options are specified by settings. Gradient panics
-// if len(deriv) != len(x).
+// sample locations, and other options are specified by settings.
+// If the step size is zero, then the step size of the method will
+// be used.
+// Gradient panics if len(deriv) != len(x).
 func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradient []float64) {
 	if len(gradient) != len(x) {
 		panic("fd: location and gradient length mismatch")
@@ -101,6 +108,10 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 	if settings == nil {
 		settings = DefaultSettings()
 	}
+	step := settings.Step
+	if settings.Step == 0 {
+		step = settings.Method.Step
+	}
 	if !settings.Concurrent {
 		xcopy := make([]float64, len(x)) // So that x is not modified during the call
 		copy(xcopy, x)
@@ -111,11 +122,11 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 					deriv += pt.Coeff * settings.OriginValue
 					continue
 				}
-				xcopy[i] += pt.Loc * settings.Step
+				xcopy[i] += pt.Loc * step
 				deriv += pt.Coeff * f(xcopy)
 				xcopy[i] = x[i]
 			}
-			gradient[i] = deriv / math.Pow(settings.Step, float64(settings.Method.Order))
+			gradient[i] = deriv / math.Pow(step, float64(settings.Method.Order))
 		}
 		return
 	}
@@ -141,7 +152,7 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 				case <-quit:
 					return
 				case run := <-sendChan:
-					xcopy[run.idx] += run.pt.Loc * settings.Step
+					xcopy[run.idx] += run.pt.Loc * step
 					run.result = f(xcopy)
 					xcopy[run.idx] = x[run.idx]
 					ansChan <- run
@@ -180,7 +191,7 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 		run := <-ansChan
 		gradient[run.idx] += run.pt.Coeff * run.result
 	}
-	floats.Scale(1/math.Pow(settings.Step, float64(settings.Method.Order)), gradient)
+	floats.Scale(1/math.Pow(step, float64(settings.Method.Order)), gradient)
 }
 
 type fdrun struct {
@@ -193,22 +204,26 @@ type fdrun struct {
 var Forward = Method{
 	Stencil: []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
 	Order:   1,
+	Step:    1e-6,
 }
 
 // Backward represents a first-order backward difference
 var Backward = Method{
 	Stencil: []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
 	Order:   1,
+	Step:    1e-6,
 }
 
 // Central represents a first-order central difference.
 var Central = Method{
 	Stencil: []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
 	Order:   1,
+	Step:    1e-6,
 }
 
 // Central2nd represents a secord-order central difference.
 var Central2nd = Method{
 	Stencil: []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
 	Order:   2,
+	Step:    1e-3,
 }
diff --git a/fd/diff_test.go b/fd/diff_test.go
index 00848394..39684fa0 100644
--- a/fd/diff_test.go
+++ b/fd/diff_test.go
@@ -16,7 +16,6 @@ type testPoint struct {
 	loc  float64
 	fofx float64
 	ans  float64
-	step float64
 }
 
 var testsFirst = []testPoint{
@@ -25,28 +24,24 @@ var testsFirst = []testPoint{
 		loc:  0,
 		fofx: 0,
 		ans:  0,
-		step: 1e-6,
 	},
 	{
 		f:    xSquared,
 		loc:  5,
 		fofx: 25,
 		ans:  10,
-		step: 1e-6,
 	},
 	{
 		f:    xSquared,
 		loc:  2,
 		fofx: 4,
 		ans:  4,
-		step: 1e-6,
 	},
 	{
 		f:    xSquared,
 		loc:  -5,
 		fofx: 25,
 		ans:  -10,
-		step: 1e-6,
 	},
 }
 
@@ -56,28 +51,24 @@ var testsSecond = []testPoint{
 		loc:  0,
 		fofx: 0,
 		ans:  2,
-		step: 1e-3,
 	},
 	{
 		f:    xSquared,
 		loc:  5,
 		fofx: 25,
 		ans:  2,
-		step: 1e-3,
 	},
 	{
 		f:    xSquared,
 		loc:  2,
 		fofx: 4,
 		ans:  2,
-		step: 1e-3,
 	},
 	{
 		f:    xSquared,
 		loc:  -5,
 		fofx: 25,
 		ans:  2,
-		step: 1e-3,
 	},
 }
 
@@ -86,7 +77,6 @@ func testDerivative(t *testing.T, method Method, tol float64, tests []testPoint)
 
 		settings := DefaultSettings()
 		settings.Method = method
-		settings.Step = 1e-4
 		ans := Derivative(test.f, test.loc, settings)
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("Case %v: ans mismatch serial: expected %v, found %v", i, test.ans, ans)

From b672c9a73e10fab290a22572ad997de866a8a7e0 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Wed, 15 Oct 2014 18:10:32 -0700
Subject: [PATCH 14/98] Changed step to only exist in Method

---
 fd/diff.go | 16 ++++------------
 1 file changed, 4 insertions(+), 12 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index ee0e3ad8..aa2d7088 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -32,7 +32,6 @@ type Method struct {
 type Settings struct {
 	OriginKnown bool    // Flag that the value at the origin x is known
 	OriginValue float64 // Value at the origin (only used if OriginKnown is true)
-	Step        float64 // step size
 	Concurrent  bool    // Should the function calls be executed concurrently
 	Workers     int     // Maximum number of concurrent executions when evaluating concurrently
 	Method      Method  // Finite difference method to use
@@ -50,16 +49,12 @@ func DefaultSettings() *Settings {
 
 // Derivative estimates the derivative of the function f at the given location.
 // The order of derivative, sample locations, and other options are specified
-// by settings. If the step is zero, then the step size of the method will
-// be used.
+// by settings.
 func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
 	if settings == nil {
 		settings = DefaultSettings()
 	}
-	step := settings.Step
-	if settings.Step == 0 {
-		step = settings.Method.Step
-	}
+	step := settings.Method.Step
 	var deriv float64
 	method := settings.Method
 	if !settings.Concurrent {
@@ -108,10 +103,7 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 	if settings == nil {
 		settings = DefaultSettings()
 	}
-	step := settings.Step
-	if settings.Step == 0 {
-		step = settings.Method.Step
-	}
+	step := settings.Method.Step
 	if !settings.Concurrent {
 		xcopy := make([]float64, len(x)) // So that x is not modified during the call
 		copy(xcopy, x)
@@ -207,7 +199,7 @@ var Forward = Method{
 	Step:    1e-6,
 }
 
-// Backward represents a first-order backward difference
+// Backward represents a first-order backward difference.
 var Backward = Method{
 	Stencil: []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
 	Order:   1,

From fb4a3a27cac951800fd36af592850a75aa995dd0 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 16:05:32 -0400
Subject: [PATCH 15/98] Create test-coverage.sh

---
 test-coverage.sh | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)
 create mode 100644 test-coverage.sh

diff --git a/test-coverage.sh b/test-coverage.sh
new file mode 100644
index 00000000..7d2be035
--- /dev/null
+++ b/test-coverage.sh
@@ -0,0 +1,30 @@
+#!/bin/bash
+
+# based on http://stackoverflow.com/questions/21126011/is-it-possible-to-post-coverage-for-multiple-packages-to-coveralls
+# with script found at https://github.com/gopns/gopns/blob/master/test-coverage.sh
+
+echo "mode: set" > acc.out
+for Dir in $(find ./* -maxdepth 10 -type d ); 
+do
+	if ls $Dir/*.go &> /dev/null;
+	then
+		returnval=`go test -v -coverprofile=profile.out $Dir`
+		echo ${returnval}
+		if [[ ${returnval} != *FAIL* ]]
+		then
+    		if [ -f profile.out ]
+    		then
+        		cat profile.out | grep -v "mode: set" >> acc.out 
+    		fi
+    	else
+    		exit 1
+    	fi	
+    fi
+done
+if [ -n "$COVERALLS_TOKEN" ]
+then
+	goveralls -coverprofile=acc.out $COVERALLS_TOKEN
+fi	
+
+rm -rf ./profile.out
+rm -rf ./acc.out

From 245f4dcc58a27db9cb6d2c4b548f1e6615198b08 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 16:07:19 -0400
Subject: [PATCH 16/98] Update test-coverage.sh

---
 test-coverage.sh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test-coverage.sh b/test-coverage.sh
index 7d2be035..9e7e50c0 100644
--- a/test-coverage.sh
+++ b/test-coverage.sh
@@ -23,7 +23,7 @@ do
 done
 if [ -n "$COVERALLS_TOKEN" ]
 then
-	goveralls -coverprofile=acc.out $COVERALLS_TOKEN
+	$HOME/gopath/bin/goveralls -coverprofile=acc.out -service=travis-ci -repotoken $COVERALLS_TOKEN
 fi	
 
 rm -rf ./profile.out

From 7ea8280f0283d44883e74e44da7cbaa637989d8e Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 16:07:42 -0400
Subject: [PATCH 17/98] Create .travis.yml

---
 .travis.yml | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)
 create mode 100644 .travis.yml

diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 00000000..9f6302c2
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,30 @@
+language: go
+
+go:
+ - release
+ - tip
+
+#env:
+#  global:
+#  - secure: "aq8rblr2+Ga10psfUJAKBCVtSIv/Uca7Ga1eTzRnRuGPsCHhC691jPmNAd8P7UhmhZ6ZHNN6yC+yYsS+lzCm62gisyIMA3+NmOM4Y2qf3wTSVEvnjAdVonge8aB3zWk5XimhdwiXWoBtW4+cVrNYWv+uNw7gNzGhND0pWFFU//c="
+
+
+before_install:
+ - go get code.google.com/p/go.tools/cmd/cover
+ - go get github.com/mattn/goveralls
+
+script:
+ - go get -d -v ./...
+ - go build -v ./...
+ - go test -v ./...
+ - diff <(gofmt -d .) <("")
+ - ./test-coverage.sh
+
+after_failure: failure
+
+notifications:
+  email:
+    recipients:
+      - jonathan.lawlor@gmail.com
+    on_success: change
+    on_failure: always

From 8fc45237f90c97108577066b9526e08b331d1d53 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 16:09:08 -0400
Subject: [PATCH 18/98] Update README.md

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index f0402485..780af4b7 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Gonum diff
+# Gonum diff  [![Build Status](https://travis-ci.org/jonlawlor/diff.svg)](https://travis-ci.org/jonlawlor/diff)
 
 This is a package for computing derivatives of functions for the Go language.
 

From f2c83425360dfcadaf11df3332c11c91f74f73a3 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 16:12:41 -0400
Subject: [PATCH 19/98] Update .travis.yml

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 9f6302c2..dc4d335b 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -18,7 +18,7 @@ script:
  - go build -v ./...
  - go test -v ./...
  - diff <(gofmt -d .) <("")
- - ./test-coverage.sh
+ - sudo ./test-coverage.sh
 
 after_failure: failure
 

From eb7d19bfec6e8a2d3222dc4af8fe046fcb6d23be Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 18:12:06 -0400
Subject: [PATCH 20/98] Update .travis.yml

---
 .travis.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.travis.yml b/.travis.yml
index dc4d335b..6dcfa12c 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -18,6 +18,8 @@ script:
  - go build -v ./...
  - go test -v ./...
  - diff <(gofmt -d .) <("")
+ - sudo ls -la
+ - sudo pwd
  - sudo ./test-coverage.sh
 
 after_failure: failure

From 3494d6c75a687102d05c927bd5cd6621fb3b60df Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 18:20:35 -0400
Subject: [PATCH 21/98] Update .travis.yml

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 6dcfa12c..bee1dab0 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -20,7 +20,7 @@ script:
  - diff <(gofmt -d .) <("")
  - sudo ls -la
  - sudo pwd
- - sudo ./test-coverage.sh
+ - sudo bash test-coverage.sh
 
 after_failure: failure
 

From eabf39fdc1ef7898f20cd3af38b793a6effd05b8 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 18:27:02 -0400
Subject: [PATCH 22/98] add coveralls.io encrypted key

---
 .travis.yml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index bee1dab0..78c455e2 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,10 +4,10 @@ go:
  - release
  - tip
 
-#env:
-#  global:
-#  - secure: "aq8rblr2+Ga10psfUJAKBCVtSIv/Uca7Ga1eTzRnRuGPsCHhC691jPmNAd8P7UhmhZ6ZHNN6yC+yYsS+lzCm62gisyIMA3+NmOM4Y2qf3wTSVEvnjAdVonge8aB3zWk5XimhdwiXWoBtW4+cVrNYWv+uNw7gNzGhND0pWFFU//c="
-
+env:
+  global:
+  - secure: "bbvnev0lw5pf9wY/Fv9CgDCynS0q8wvAiRXC/wPwaTR/X7bWwuQF7YmBhOlS9H/x1IxHW1HyqJMjOipIF2/MgPoTB4xVkk1jLHigd6MCt91o8aZvkaW4ITDKhP29CZNCyBKvPoEs/yLGEfaGYop3ivfJPWnDnMAd04m7+Yr1+3I="
+  
 
 before_install:
  - go get code.google.com/p/go.tools/cmd/cover

From 8b9cf0f6ed56877e45dc077d9e251a2d692353d7 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 18:27:42 -0400
Subject: [PATCH 23/98] Add coverage flair to readme.md

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 780af4b7..1ad829e3 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Gonum diff  [![Build Status](https://travis-ci.org/jonlawlor/diff.svg)](https://travis-ci.org/jonlawlor/diff)
+# Gonum diff  [![Build Status](https://travis-ci.org/jonlawlor/diff.svg)](https://travis-ci.org/jonlawlor/diff)  [![Coverage Status](https://img.shields.io/coveralls/gonum/diff.svg)](https://coveralls.io/r/gonum/diff)
 
 This is a package for computing derivatives of functions for the Go language.
 

From f0db9d07448aeafcd81bec16abf66d414db50f6f Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 18:33:44 -0400
Subject: [PATCH 24/98] Fix incorrect link in readme.me (s/jonlawlor/gonum)

Accidentally used the fork's link.
---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 1ad829e3..959e7322 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Gonum diff  [![Build Status](https://travis-ci.org/jonlawlor/diff.svg)](https://travis-ci.org/jonlawlor/diff)  [![Coverage Status](https://img.shields.io/coveralls/gonum/diff.svg)](https://coveralls.io/r/gonum/diff)
+# Gonum diff  [![Build Status](https://travis-ci.org/gonum/diff.svg)](https://travis-ci.org/gonum/diff)  [![Coverage Status](https://img.shields.io/coveralls/gonum/diff.svg)](https://coveralls.io/r/gonum/diff)
 
 This is a package for computing derivatives of functions for the Go language.
 

From 3ce5e2d8c0c86b701807286e9c37c4b23c387cb9 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 20:05:27 -0400
Subject: [PATCH 25/98] remove sudo from test-cover.sh call

No need to use sudo, and it resets the path so it can't call the go tool.
---
 .travis.yml | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index 78c455e2..e650c393 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -18,9 +18,7 @@ script:
  - go build -v ./...
  - go test -v ./...
  - diff <(gofmt -d .) <("")
- - sudo ls -la
- - sudo pwd
- - sudo bash test-coverage.sh
+ - bash test-coverage.sh
 
 after_failure: failure
 

From 3c09a87e65dd3af2a5ed4a90cf9e056e16d333d4 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Oct 2014 23:09:45 -0400
Subject: [PATCH 26/98] Update test-coverage.sh

added coverage for root package (currently unused)
---
 test-coverage.sh | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/test-coverage.sh b/test-coverage.sh
index 9e7e50c0..3df03f52 100644
--- a/test-coverage.sh
+++ b/test-coverage.sh
@@ -4,10 +4,23 @@
 # with script found at https://github.com/gopns/gopns/blob/master/test-coverage.sh
 
 echo "mode: set" > acc.out
+returnval=`go test -v -coverprofile=profile.out`
+echo ${returnval}
+if [[ ${returnval} != *FAIL* ]]
+then
+	if [ -f profile.out ]
+	then
+		cat profile.out | grep -v "mode: set" >> acc.out 
+	fi
+else
+	exit 1
+fi	
+
 for Dir in $(find ./* -maxdepth 10 -type d ); 
 do
 	if ls $Dir/*.go &> /dev/null;
 	then
+		echo $Dir
 		returnval=`go test -v -coverprofile=profile.out $Dir`
 		echo ${returnval}
 		if [[ ${returnval} != *FAIL* ]]

From bc1f097e7c6553c33126dd7ebdd3a4a1ca7dfd38 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Sat, 18 Oct 2014 22:49:28 -0400
Subject: [PATCH 27/98] add tests to improve test coverage

Added tests which cover cases with default settings, the panic in
Gradient, and one where more workers are asked for than can be used.
Test coverage is 100%.
---
 fd/diff_test.go     | 11 +++++++++++
 fd/gradient_test.go | 31 +++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+)

diff --git a/fd/diff_test.go b/fd/diff_test.go
index 39684fa0..d9f51c6b 100644
--- a/fd/diff_test.go
+++ b/fd/diff_test.go
@@ -119,3 +119,14 @@ func TestCentral(t *testing.T) {
 func TestCentralSecond(t *testing.T) {
 	testDerivative(t, Central2nd, 1e-3, testsSecond)
 }
+
+// TestDerivativeDefault checks that the derivative works when settings is set to nil.
+func TestDerivativeDefault(t *testing.T) {
+	tol := 1e-6
+	for i, test := range testsFirst {
+		ans := Derivative(test.f, test.loc, nil)
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch default: expected %v, found %v", i, test.ans, ans)
+		}
+	}
+}
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 8b432b67..4ee52db8 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -110,6 +110,7 @@ func TestGradient(t *testing.T) {
 		}
 		settings.Concurrent = true
 		settings.OriginKnown = false
+		settings.Workers = 1000
 		Gradient(r.F, x, settings, gradient)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
@@ -128,5 +129,35 @@ func TestGradient(t *testing.T) {
 			t.Errorf("Case %v: gradient mismatch with known origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 
+		// With default settings
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		settings = nil
+		Gradient(r.F, x, settings, gradient)
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with default settings. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+
+	}
+}
+
+func Panics(fun func()) (b bool) {
+	defer func() {
+		err := recover()
+		if err != nil {
+			b = true
+		}
+	}()
+	fun()
+	return
+}
+
+func TestGradientPanics(t *testing.T) {
+	// Test that it panics
+	if !Panics(func() {
+		Gradient(func(x []float64) float64 { return x[0] * x[0] }, []float64{0.0, 0.0}, nil, []float64{0.0})
+	}) {
+		t.Errorf("Gradient did not panic with length mismatch")
 	}
 }

From 8dce656a703931df344b6f12c6b9fc41ea82a572 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Tue, 4 Nov 2014 21:51:18 -0500
Subject: [PATCH 28/98] add conditional for test coverage

travis-ci only has secure variables available when a pull is from a
branch of gonum.  This conditional makes it so that coverage is only
calculated when the coveralls api key is available.
---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index e650c393..f041e041 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -18,7 +18,7 @@ script:
  - go build -v ./...
  - go test -v ./...
  - diff <(gofmt -d .) <("")
- - bash test-coverage.sh
+ - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash test-coverage.sh; fi
 
 after_failure: failure
 

From 4a2097a8f197fd93f9bdfcb11ab6c02fa0067bf5 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 13 Nov 2014 01:17:36 -0500
Subject: [PATCH 29/98] use new location for go cover tool

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index f041e041..391191a4 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -10,7 +10,7 @@ env:
   
 
 before_install:
- - go get code.google.com/p/go.tools/cmd/cover
+ - if ! go get code.google.com/p/go.tools/cmd/cover; then go get golang.org/x/tools/cmd/cover; fi
  - go get github.com/mattn/goveralls
 
 script:

From c4a060a0d4277c0b29aee130ce21d262f92e602f Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Tue, 9 Dec 2014 21:05:52 -0500
Subject: [PATCH 30/98] get packages only used in testing

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 391191a4..a4052d56 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -14,7 +14,7 @@ before_install:
  - go get github.com/mattn/goveralls
 
 script:
- - go get -d -v ./...
+ - go get -d -t -v ./...
  - go build -v ./...
  - go test -v ./...
  - diff <(gofmt -d .) <("")

From 793314b469710397dea45c208bdd9ab1075b5100 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 18 Dec 2014 21:42:57 -0500
Subject: [PATCH 31/98] disable sudo to use new travis servers

This change should speed up builds by ~1m and should reduce queueing.  See http://blog.travis-ci.com/2014-12-17-faster-builds-with-container-based-infrastructure/ for more details.
---
 .travis.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.travis.yml b/.travis.yml
index a4052d56..ad7dec80 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,3 +1,5 @@
+sudo: false
+
 language: go
 
 go:

From ba3cee67d1602215d85f9a151429f476ef32c2dc Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 1 Jan 2015 18:40:35 -0500
Subject: [PATCH 32/98] Add explicit go version 1.3.3 to build matrix

The build matrix represents the versions of go we are officially supporting.
---
 .travis.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.travis.yml b/.travis.yml
index ad7dec80..3225fb03 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,6 +5,7 @@ language: go
 go:
  - release
  - tip
+ - 1.3.3
 
 env:
   global:

From 21b8a5494f431dec9c61e6e7a659d0d5124e478b Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Fri, 9 Jan 2015 21:04:19 -0800
Subject: [PATCH 33/98] Changed gradient sig to move gradient first.

The convention is to put receivers first, and allocate if nil. Previously, Gradient had not followed this convention. This moves the []float64 receiver to be the first argument, allocates it if nil, and returns the value from the function. This brings it into convention with the other packages. This PR also adds a test for the nil case.
---
 fd/diff.go          |  8 ++++++--
 fd/gradient_test.go | 21 +++++++++++++++------
 2 files changed, 21 insertions(+), 8 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index aa2d7088..3639f05e 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -96,7 +96,10 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 // If the step size is zero, then the step size of the method will
 // be used.
 // Gradient panics if len(deriv) != len(x).
-func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradient []float64) {
+func Gradient(gradient []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
+	if gradient == nil {
+		gradient = make([]float64, len(x))
+	}
 	if len(gradient) != len(x) {
 		panic("fd: location and gradient length mismatch")
 	}
@@ -120,7 +123,7 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 			}
 			gradient[i] = deriv / math.Pow(step, float64(settings.Method.Order))
 		}
-		return
+		return gradient
 	}
 
 	nWorkers := settings.Workers
@@ -184,6 +187,7 @@ func Gradient(f func([]float64) float64, x []float64, settings *Settings, gradie
 		gradient[run.idx] += run.pt.Coeff * run.result
 	}
 	floats.Scale(1/math.Pow(step, float64(settings.Method.Order)), gradient)
+	return gradient
 }
 
 type fdrun struct {
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 8b432b67..e2fbb06d 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -80,17 +80,26 @@ func TestGradient(t *testing.T) {
 
 		settings := DefaultSettings()
 		settings.Method = test.method
-		gradient := make([]float64, len(x))
+
+		// try with gradient nil
+		gradient := Gradient(nil, r.F, x, settings)
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch in serial with nil. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %v: x modified during call to gradient in serial with nil.", i)
+		}
+
 		for i := range gradient {
 			gradient[i] = rand.Float64()
 		}
 
-		Gradient(r.F, x, settings, gradient)
+		Gradient(gradient, r.F, x, settings)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %v: x modified during call to gradient in serial", i)
+			t.Errorf("Case %v: x modified during call to gradient in serial with non-nil.", i)
 		}
 
 		// Try with known value
@@ -99,7 +108,7 @@ func TestGradient(t *testing.T) {
 		}
 		settings.OriginKnown = true
 		settings.OriginValue = r.F(x)
-		Gradient(r.F, x, settings, gradient)
+		Gradient(gradient, r.F, x, settings)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with known origin in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
@@ -110,7 +119,7 @@ func TestGradient(t *testing.T) {
 		}
 		settings.Concurrent = true
 		settings.OriginKnown = false
-		Gradient(r.F, x, settings, gradient)
+		Gradient(gradient, r.F, x, settings)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
@@ -123,7 +132,7 @@ func TestGradient(t *testing.T) {
 			gradient[i] = rand.Float64()
 		}
 		settings.OriginKnown = true
-		Gradient(r.F, x, settings, gradient)
+		Gradient(gradient, r.F, x, settings)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with known origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}

From b8526fb15c2bdc4f950a751dd23dc40018bfc250 Mon Sep 17 00:00:00 2001
From: btracey <bdtinternet@gmail.com>
Date: Fri, 9 Jan 2015 21:24:11 -0800
Subject: [PATCH 34/98] fixed code broken in pull from master

---
 fd/gradient_test.go      |  10 +--
 fd/gradient_test.go.orig | 176 ---------------------------------------
 2 files changed, 3 insertions(+), 183 deletions(-)
 delete mode 100644 fd/gradient_test.go.orig

diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 821bdba2..61d18e44 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -119,12 +119,8 @@ func TestGradient(t *testing.T) {
 		}
 		settings.Concurrent = true
 		settings.OriginKnown = false
-<<<<<<< HEAD
-		Gradient(gradient, r.F, x, settings)
-=======
 		settings.Workers = 1000
-		Gradient(r.F, x, settings, gradient)
->>>>>>> ba3cee67d1602215d85f9a151429f476ef32c2dc
+		Gradient(gradient, r.F, x, settings)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
@@ -147,7 +143,7 @@ func TestGradient(t *testing.T) {
 			gradient[i] = rand.Float64()
 		}
 		settings = nil
-		Gradient(r.F, x, settings, gradient)
+		Gradient(gradient, r.F, x, settings)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with default settings. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
@@ -169,7 +165,7 @@ func Panics(fun func()) (b bool) {
 func TestGradientPanics(t *testing.T) {
 	// Test that it panics
 	if !Panics(func() {
-		Gradient(func(x []float64) float64 { return x[0] * x[0] }, []float64{0.0, 0.0}, nil, []float64{0.0})
+		Gradient([]float64{0.0}, func(x []float64) float64 { return x[0] * x[0] }, []float64{0.0, 0.0}, nil)
 	}) {
 		t.Errorf("Gradient did not panic with length mismatch")
 	}
diff --git a/fd/gradient_test.go.orig b/fd/gradient_test.go.orig
deleted file mode 100644
index 821bdba2..00000000
--- a/fd/gradient_test.go.orig
+++ /dev/null
@@ -1,176 +0,0 @@
-// Copyright ©2014 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 fd
-
-import (
-	"math"
-	"math/rand"
-	"testing"
-
-	"github.com/gonum/floats"
-)
-
-type Rosenbrock struct {
-	nDim int
-}
-
-func (r Rosenbrock) F(x []float64) (sum float64) {
-	deriv := make([]float64, len(x))
-	return r.FDf(x, deriv)
-}
-
-func (r Rosenbrock) FDf(x []float64, deriv []float64) (sum float64) {
-	for i := range deriv {
-		deriv[i] = 0
-	}
-
-	for i := 0; i < len(x)-1; i++ {
-		sum += math.Pow(1-x[i], 2) + 100*math.Pow(x[i+1]-math.Pow(x[i], 2), 2)
-	}
-	for i := 0; i < len(x)-1; i++ {
-		deriv[i] += -1 * 2 * (1 - x[i])
-		deriv[i] += 2 * 100 * (x[i+1] - math.Pow(x[i], 2)) * (-2 * x[i])
-	}
-	for i := 1; i < len(x); i++ {
-		deriv[i] += 2 * 100 * (x[i] - math.Pow(x[i-1], 2))
-	}
-
-	return sum
-}
-
-func TestGradient(t *testing.T) {
-	for i, test := range []struct {
-		nDim   int
-		tol    float64
-		method Method
-	}{
-		{
-			nDim:   2,
-			tol:    2e-4,
-			method: Forward,
-		},
-		{
-			nDim:   2,
-			tol:    1e-6,
-			method: Central,
-		},
-		{
-			nDim:   40,
-			tol:    2e-4,
-			method: Forward,
-		},
-		{
-			nDim:   40,
-			tol:    1e-6,
-			method: Central,
-		},
-	} {
-		x := make([]float64, test.nDim)
-		for i := range x {
-			x[i] = rand.Float64()
-		}
-		xcopy := make([]float64, len(x))
-		copy(xcopy, x)
-
-		r := Rosenbrock{len(x)}
-		trueGradient := make([]float64, len(x))
-		r.FDf(x, trueGradient)
-
-		settings := DefaultSettings()
-		settings.Method = test.method
-
-		// try with gradient nil
-		gradient := Gradient(nil, r.F, x, settings)
-		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
-			t.Errorf("Case %v: gradient mismatch in serial with nil. Want: %v, Got: %v.", i, trueGradient, gradient)
-		}
-		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %v: x modified during call to gradient in serial with nil.", i)
-		}
-
-		for i := range gradient {
-			gradient[i] = rand.Float64()
-		}
-
-		Gradient(gradient, r.F, x, settings)
-		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
-			t.Errorf("Case %v: gradient mismatch in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
-		}
-		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %v: x modified during call to gradient in serial with non-nil.", i)
-		}
-
-		// Try with known value
-		for i := range gradient {
-			gradient[i] = rand.Float64()
-		}
-		settings.OriginKnown = true
-		settings.OriginValue = r.F(x)
-		Gradient(gradient, r.F, x, settings)
-		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
-			t.Errorf("Case %v: gradient mismatch with known origin in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
-		}
-
-		// Concurrently
-		for i := range gradient {
-			gradient[i] = rand.Float64()
-		}
-		settings.Concurrent = true
-		settings.OriginKnown = false
-<<<<<<< HEAD
-		Gradient(gradient, r.F, x, settings)
-=======
-		settings.Workers = 1000
-		Gradient(r.F, x, settings, gradient)
->>>>>>> ba3cee67d1602215d85f9a151429f476ef32c2dc
-		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
-			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
-		}
-		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %v: x modified during call to gradient in parallel", i)
-		}
-
-		// Concurrently with origin known
-		for i := range gradient {
-			gradient[i] = rand.Float64()
-		}
-		settings.OriginKnown = true
-		Gradient(gradient, r.F, x, settings)
-		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
-			t.Errorf("Case %v: gradient mismatch with known origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
-		}
-
-		// With default settings
-		for i := range gradient {
-			gradient[i] = rand.Float64()
-		}
-		settings = nil
-		Gradient(r.F, x, settings, gradient)
-		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
-			t.Errorf("Case %v: gradient mismatch with default settings. Want: %v, Got: %v.", i, trueGradient, gradient)
-		}
-
-	}
-}
-
-func Panics(fun func()) (b bool) {
-	defer func() {
-		err := recover()
-		if err != nil {
-			b = true
-		}
-	}()
-	fun()
-	return
-}
-
-func TestGradientPanics(t *testing.T) {
-	// Test that it panics
-	if !Panics(func() {
-		Gradient(func(x []float64) float64 { return x[0] * x[0] }, []float64{0.0, 0.0}, nil, []float64{0.0})
-	}) {
-		t.Errorf("Gradient did not panic with length mismatch")
-	}
-}

From 06e3f735faca04bee3885f4dff073db30f80609b Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Mon, 2 Feb 2015 17:40:41 -0500
Subject: [PATCH 35/98] Add explicit go version 1.4.1 to build matrix

---
 .travis.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index 3225fb03..1252883f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -3,9 +3,9 @@ sudo: false
 language: go
 
 go:
- - release
- - tip
  - 1.3.3
+ - 1.4.1
+ - tip
 
 env:
   global:

From d959ccddba00ca61d5347399146f5616d1d9ce40 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Wed, 18 Feb 2015 07:52:52 -0500
Subject: [PATCH 36/98] Change go version 1.4.1 -> 1.4.2 in build matrix

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 1252883f..994254a3 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 
 go:
  - 1.3.3
- - 1.4.1
+ - 1.4.2
  - tip
 
 env:

From 99e7f3656ffff3c2dcaf55bdc2ee4bd34deaba3b Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 4 Jun 2015 22:43:41 -0400
Subject: [PATCH 37/98] cleanup .travis.yml and associated files

---
 .travis.yml                                  | 25 +++++---------------
 test-coverage.sh => .travis/test-coverage.sh | 12 +++++-----
 2 files changed, 12 insertions(+), 25 deletions(-)
 rename test-coverage.sh => .travis/test-coverage.sh (82%)
 mode change 100644 => 100755

diff --git a/.travis.yml b/.travis.yml
index 994254a3..bcecab1d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,33 +1,20 @@
-sudo: false
-
 language: go
 
+# Versions of go that are explicitly supported by gonum.
 go:
  - 1.3.3
  - 1.4.2
- - tip
-
-env:
-  global:
-  - secure: "bbvnev0lw5pf9wY/Fv9CgDCynS0q8wvAiRXC/wPwaTR/X7bWwuQF7YmBhOlS9H/x1IxHW1HyqJMjOipIF2/MgPoTB4xVkk1jLHigd6MCt91o8aZvkaW4ITDKhP29CZNCyBKvPoEs/yLGEfaGYop3ivfJPWnDnMAd04m7+Yr1+3I="
-  
 
+# After go 1.3, go tools moved from code.google.com to golang.org
 before_install:
- - if ! go get code.google.com/p/go.tools/cmd/cover; then go get golang.org/x/tools/cmd/cover; fi
+ - if ! go get golang.org/x/tools/cmd/cover; then go get code.google.com/p/go.tools/cmd/cover; fi
  - go get github.com/mattn/goveralls
 
+# Get deps, build, test, and ensure the code is gofmt'ed.
+# If we are building as gonum, then we have access to the coveralls api key, so we can run coverage as well.
 script:
  - go get -d -t -v ./...
  - go build -v ./...
  - go test -v ./...
  - diff <(gofmt -d .) <("")
- - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash test-coverage.sh; fi
-
-after_failure: failure
-
-notifications:
-  email:
-    recipients:
-      - jonathan.lawlor@gmail.com
-    on_success: change
-    on_failure: always
+ - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash ./.travis/test-coverage.sh; fi
\ No newline at end of file
diff --git a/test-coverage.sh b/.travis/test-coverage.sh
old mode 100644
new mode 100755
similarity index 82%
rename from test-coverage.sh
rename to .travis/test-coverage.sh
index 3df03f52..854f8c3a
--- a/test-coverage.sh
+++ b/.travis/test-coverage.sh
@@ -10,13 +10,13 @@ if [[ ${returnval} != *FAIL* ]]
 then
 	if [ -f profile.out ]
 	then
-		cat profile.out | grep -v "mode: set" >> acc.out 
+		cat profile.out | grep -v "mode: set" >> acc.out
 	fi
 else
 	exit 1
-fi	
+fi
 
-for Dir in $(find ./* -maxdepth 10 -type d ); 
+for Dir in $(find ./* -maxdepth 10 -type d );
 do
 	if ls $Dir/*.go &> /dev/null;
 	then
@@ -27,17 +27,17 @@ do
 		then
     		if [ -f profile.out ]
     		then
-        		cat profile.out | grep -v "mode: set" >> acc.out 
+        		cat profile.out | grep -v "mode: set" >> acc.out
     		fi
     	else
     		exit 1
-    	fi	
+    	fi
     fi
 done
 if [ -n "$COVERALLS_TOKEN" ]
 then
 	$HOME/gopath/bin/goveralls -coverprofile=acc.out -service=travis-ci -repotoken $COVERALLS_TOKEN
-fi	
+fi
 
 rm -rf ./profile.out
 rm -rf ./acc.out

From 0d8bf9d3231381e12affb40a513e894966474dd9 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 4 Jun 2015 23:18:45 -0400
Subject: [PATCH 38/98] test without code.google.com

---
 .travis.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index bcecab1d..b93b9db7 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -7,7 +7,7 @@ go:
 
 # After go 1.3, go tools moved from code.google.com to golang.org
 before_install:
- - if ! go get golang.org/x/tools/cmd/cover; then go get code.google.com/p/go.tools/cmd/cover; fi
+ - go get golang.org/x/tools/cmd/cover
  - go get github.com/mattn/goveralls
 
 # Get deps, build, test, and ensure the code is gofmt'ed.
@@ -17,4 +17,4 @@ script:
  - go build -v ./...
  - go test -v ./...
  - diff <(gofmt -d .) <("")
- - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash ./.travis/test-coverage.sh; fi
\ No newline at end of file
+ - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash ./.travis/test-coverage.sh; fi

From 39807a8e6a3f4d285fdb581e50a22cb19671f23c Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 4 Jun 2015 23:23:03 -0400
Subject: [PATCH 39/98] Update .travis.yml

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index b93b9db7..43b937fc 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,7 @@ go:
  - 1.3.3
  - 1.4.2
 
-# After go 1.3, go tools moved from code.google.com to golang.org
+# Required for coverage.
 before_install:
  - go get golang.org/x/tools/cmd/cover
  - go get github.com/mattn/goveralls

From 9b7a802312846a6035104e6775eb93a3f9001a26 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Wed, 8 Jul 2015 22:57:59 -0400
Subject: [PATCH 40/98] add go 1.5beta1 to the build matrix

as per https://groups.google.com/forum/#!topic/gonum-dev/GKJjAtDxR9I
---
 .travis.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.travis.yml b/.travis.yml
index 43b937fc..d806dd0a 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,6 +4,7 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
+ - 1.5beta1
 
 # Required for coverage.
 before_install:

From 89c978481feb52601efe668379c6d7cc4732022f Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Jul 2015 18:25:53 -0400
Subject: [PATCH 41/98] update go 1.5beta1 -> go 1.5beta2

See https://groups.google.com/forum/#!topic/gonum-dev/GKJjAtDxR9I for discussion
---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index d806dd0a..54346c8e 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
- - 1.5beta1
+ - 1.5beta2
 
 # Required for coverage.
 before_install:

From e3ff293d10937069d866b021ca0d9a896e49b984 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 30 Jul 2015 11:12:56 -0400
Subject: [PATCH 42/98] update go 1.5beta2 -> go 1.5beta3

As per https://groups.google.com/forum/#!topic/gonum-dev/QS-1AaDyq4w
---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 54346c8e..59f02c47 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
- - 1.5beta2
+ - 1.5beta3
 
 # Required for coverage.
 before_install:

From 313fd6ce74b094f9e358bb26db8cadc04f38b2b9 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 6 Aug 2015 07:34:33 -0400
Subject: [PATCH 43/98] go 1.5beta3 -> go 1.5rc1

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 59f02c47..15e38b88 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
- - 1.5beta3
+ - 1.5rc1
 
 # Required for coverage.
 before_install:

From ee11909bd09862c0f6f696387dfcfcedf17bc86f Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Wed, 19 Aug 2015 13:30:52 -0400
Subject: [PATCH 44/98] go 1.5rc1 -> go 1.5

As per https://groups.google.com/forum/#!topic/gonum-dev/QS-1AaDyq4w
---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 15e38b88..9ec7a679 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
- - 1.5rc1
+ - 1.5
 
 # Required for coverage.
 before_install:

From d181d5b7839ae77031b1fbccca4798c1072903c5 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Wed, 9 Sep 2015 08:35:44 -0400
Subject: [PATCH 45/98] go 1.5 -> go 1.5.1

As per https://groups.google.com/forum/#!topic/gonum-dev/QS-1AaDyq4w
---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 9ec7a679..0d4cc2d5 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
- - 1.5
+ - 1.5.1
 
 # Required for coverage.
 before_install:

From faa2c47587d2e3d3542ba488b2c6a06fdd88413c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vladim=C3=ADr=20Chalupeck=C3=BD?=
 <vladimir-ch@users.noreply.github.com>
Date: Wed, 30 Sep 2015 11:35:04 +0900
Subject: [PATCH 46/98] Update Coveralls badge link

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 959e7322..56af6971 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Gonum diff  [![Build Status](https://travis-ci.org/gonum/diff.svg)](https://travis-ci.org/gonum/diff)  [![Coverage Status](https://img.shields.io/coveralls/gonum/diff.svg)](https://coveralls.io/r/gonum/diff)
+# Gonum diff  [![Build Status](https://travis-ci.org/gonum/diff.svg)](https://travis-ci.org/gonum/diff)  [![Coverage Status](https://coveralls.io/repos/gonum/diff/badge.svg?branch=master&service=github)](https://coveralls.io/github/gonum/diff?branch=master)
 
 This is a package for computing derivatives of functions for the Go language.
 

From ebb3682542fee745181c192917420846323d06d4 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 18 Dec 2015 16:11:19 -0500
Subject: [PATCH 47/98] add go 1.6beta1 to build matrix, go 1.5.1->go 1.5.2

---
 .travis.yml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 0d4cc2d5..04eb3abf 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,8 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
- - 1.5.1
+ - 1.5.2
+ - 1.6beta1
 
 # Required for coverage.
 before_install:

From b1683a7ec23ddd907046a1b617ab51b1df8edc4b Mon Sep 17 00:00:00 2001
From: kortschak <dan.kortschak@adelaide.edu.au>
Date: Fri, 15 Jan 2016 11:10:47 +1030
Subject: [PATCH 48/98] travis: go1.6beta1 -> go1.6beta2

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 04eb3abf..b27c03f4 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,7 @@ go:
  - 1.3.3
  - 1.4.2
  - 1.5.2
- - 1.6beta1
+ - 1.6beta2
 
 # Required for coverage.
 before_install:

From 1ef68c35efbb79e464aec965cee08bc6500470df Mon Sep 17 00:00:00 2001
From: Sebastien Binet <binet@cern.ch>
Date: Fri, 15 Jan 2016 22:08:13 +0100
Subject: [PATCH 49/98] doc: add godoc badge

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 56af6971..03ba1b8c 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Gonum diff  [![Build Status](https://travis-ci.org/gonum/diff.svg)](https://travis-ci.org/gonum/diff)  [![Coverage Status](https://coveralls.io/repos/gonum/diff/badge.svg?branch=master&service=github)](https://coveralls.io/github/gonum/diff?branch=master)
+# Gonum diff  [![Build Status](https://travis-ci.org/gonum/diff.svg)](https://travis-ci.org/gonum/diff)  [![Coverage Status](https://coveralls.io/repos/gonum/diff/badge.svg?branch=master&service=github)](https://coveralls.io/github/gonum/diff?branch=master) [![GoDoc](https://godoc.org/github.com/gonum/diff?status.svg)](https://godoc.org/github.com/gonum/diff)
 
 This is a package for computing derivatives of functions for the Go language.
 

From b53afd237b73c12b52e8738ddf24198a0f5efb64 Mon Sep 17 00:00:00 2001
From: kortschak <dan.kortschak@adelaide.edu.au>
Date: Fri, 29 Jan 2016 13:08:31 +1030
Subject: [PATCH 50/98] travis: update test matrix

go1.5.2 -> go1.5.3
go1.6beta2 -> go1.6rc1
---
 .travis.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index b27c03f4..c958670b 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,8 +4,8 @@ language: go
 go:
  - 1.3.3
  - 1.4.2
- - 1.5.2
- - 1.6beta2
+ - 1.5.3
+ - 1.6rc1
 
 # Required for coverage.
 before_install:

From cb4205294dd322e049073df43f6c74e8414804d0 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Sat, 6 Feb 2016 17:44:23 -0500
Subject: [PATCH 51/98] travis: go1.6rc1 -> go1.6rc2

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index c958670b..d20a6d2f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,7 @@ go:
  - 1.3.3
  - 1.4.2
  - 1.5.3
- - 1.6rc1
+ - 1.6rc2
 
 # Required for coverage.
 before_install:

From 9b0f67582ffabe794d3f536c253623b7dd68cd4f Mon Sep 17 00:00:00 2001
From: kortschak <dan.kortschak@adelaide.edu.au>
Date: Thu, 18 Feb 2016 09:07:38 +1030
Subject: [PATCH 52/98] travis: go1.6

Drop go1.3 support.
---
 .travis.yml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index d20a6d2f..137c7f25 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,10 +2,9 @@ language: go
 
 # Versions of go that are explicitly supported by gonum.
 go:
- - 1.3.3
  - 1.4.2
  - 1.5.3
- - 1.6rc2
+ - 1.6
 
 # Required for coverage.
 before_install:

From 92a6cbfa9483762391ff5316e2a2e6306f8892ee Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 11:08:07 +0900
Subject: [PATCH 53/98] fd: add package-level doc comment

---
 fd/diff.go | 1 +
 1 file changed, 1 insertion(+)

diff --git a/fd/diff.go b/fd/diff.go
index 3639f05e..817f6dd5 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -2,6 +2,7 @@
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
+// Package fd provides functions to approximate derivatives using finite differences.
 package fd
 
 import (

From 825987e5e14430df402673bfb41479f6bedaf3ef Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 11:47:13 +0900
Subject: [PATCH 54/98] fd: fix and update doc comment for Gradient

---
 fd/diff.go | 36 ++++++++++++++++++------------------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 817f6dd5..750cef6d 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -91,18 +91,18 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 	return deriv / math.Pow(step, float64(method.Order))
 }
 
-// Gradient estimates the derivative of a multivariate function f at the location
-// x. The resulting estimate is stored in-place into gradient. The order of derivative,
-// sample locations, and other options are specified by settings.
-// If the step size is zero, then the step size of the method will
-// be used.
-// Gradient panics if len(deriv) != len(x).
-func Gradient(gradient []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
-	if gradient == nil {
-		gradient = make([]float64, len(x))
+// Gradient estimates the gradient of the multivariate function f at the
+// location x. The result is stored in-place into dst if dst is not nil,
+// otherwise a new slice will be allocated and returned. Finite difference
+// kernel and other options are specified by settings. If settings is nil,
+// default settings will be used.
+// Gradient panics if the length of dst and x is not equal.
+func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
+	if dst == nil {
+		dst = make([]float64, len(x))
 	}
-	if len(gradient) != len(x) {
-		panic("fd: location and gradient length mismatch")
+	if len(dst) != len(x) {
+		panic("fd: slice length mismatch")
 	}
 	if settings == nil {
 		settings = DefaultSettings()
@@ -122,9 +122,9 @@ func Gradient(gradient []float64, f func([]float64) float64, x []float64, settin
 				deriv += pt.Coeff * f(xcopy)
 				xcopy[i] = x[i]
 			}
-			gradient[i] = deriv / math.Pow(step, float64(settings.Method.Order))
+			dst[i] = deriv / math.Pow(step, float64(settings.Method.Order))
 		}
-		return gradient
+		return dst
 	}
 
 	nWorkers := settings.Workers
@@ -179,16 +179,16 @@ func Gradient(gradient []float64, f func([]float64) float64, x []float64, settin
 		}
 	}(sendChan, ansChan)
 
-	for i := range gradient {
-		gradient[i] = 0
+	for i := range dst {
+		dst[i] = 0
 	}
 	// Read in all of the results
 	for i := 0; i < expect; i++ {
 		run := <-ansChan
-		gradient[run.idx] += run.pt.Coeff * run.result
+		dst[run.idx] += run.pt.Coeff * run.result
 	}
-	floats.Scale(1/math.Pow(step, float64(settings.Method.Order)), gradient)
-	return gradient
+	floats.Scale(1/math.Pow(step, float64(settings.Method.Order)), dst)
+	return dst
 }
 
 type fdrun struct {

From 7ea8f1c931a98f0d8b26136912f9a99b50711bbb Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 11:56:45 +0900
Subject: [PATCH 55/98] fd: update doc comment for Derivative

---
 fd/diff.go | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 750cef6d..7c6b0202 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -49,8 +49,8 @@ func DefaultSettings() *Settings {
 }
 
 // Derivative estimates the derivative of the function f at the given location.
-// The order of derivative, sample locations, and other options are specified
-// by settings.
+// The order of the derivative, sample locations, and other options are
+// specified by settings. If settings is nil, default settings will be used.
 func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
 	if settings == nil {
 		settings = DefaultSettings()

From 6192d3c92e6840a6f5b5e2f6e6dd25495bec120d Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 14:56:16 +0900
Subject: [PATCH 56/98] fd: clean up comment for Method

---
 fd/diff.go | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 7c6b0202..15c90fe9 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -13,19 +13,21 @@ import (
 	"github.com/gonum/floats"
 )
 
-// A Point is a stencil location in a difference method.
+// A Point is a stencil location in a finite difference formula.
 type Point struct {
 	Loc   float64
 	Coeff float64
 }
 
-// Method is a specific finite difference method. Method specifies the stencil,
-// that is, the function locations (relative to x) which will be used to estimate
-// the derivative. It also specifies the order of derivative it estimates. Order = 1
-// represents the derivative, Order = 2 represents the curvature, etc.
+// Method represents a finite difference formula that approximates
+// the derivative of order k of a function f at x as
+//  d^k f(x) ≅ (1 / h^k) * \sum_i Coeff_i * f(x + h * Loc_i),
+// where h is a small positive step.
 type Method struct {
+	// Stencil is the set of sampling Points which are used to estimate the
+	// derivative. The locations will be scaled by Step and are relative to x.
 	Stencil []Point
-	Order   int     // The order of the difference method (first derivative, second derivative, etc.)
+	Order   int     // The order of the approximated derivative.
 	Step    float64 // Default step size for the method.
 }
 

From 8e626f9827eafad8191512140a2b20b11e5b7b62 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 15:34:58 +0900
Subject: [PATCH 57/98] fd: rename Method to Formula

---
 fd/diff.go          | 42 +++++++++++++++++++++---------------------
 fd/diff_test.go     |  4 ++--
 fd/gradient_test.go | 32 ++++++++++++++++----------------
 3 files changed, 39 insertions(+), 39 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 15c90fe9..ac2bf5fc 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -19,16 +19,16 @@ type Point struct {
 	Coeff float64
 }
 
-// Method represents a finite difference formula that approximates
+// Formula represents a finite difference formula that approximates
 // the derivative of order k of a function f at x as
 //  d^k f(x) ≅ (1 / h^k) * \sum_i Coeff_i * f(x + h * Loc_i),
 // where h is a small positive step.
-type Method struct {
+type Formula struct {
 	// Stencil is the set of sampling Points which are used to estimate the
 	// derivative. The locations will be scaled by Step and are relative to x.
 	Stencil []Point
 	Order   int     // The order of the approximated derivative.
-	Step    float64 // Default step size for the method.
+	Step    float64 // Default step size for the formula.
 }
 
 // Settings is the settings structure for computing finite differences.
@@ -37,7 +37,7 @@ type Settings struct {
 	OriginValue float64 // Value at the origin (only used if OriginKnown is true)
 	Concurrent  bool    // Should the function calls be executed concurrently
 	Workers     int     // Maximum number of concurrent executions when evaluating concurrently
-	Method      Method  // Finite difference method to use
+	Formula     Formula // Finite difference formula to use
 }
 
 // DefaultSettings is a basic set of settings for computing finite differences.
@@ -45,7 +45,7 @@ type Settings struct {
 // of the function.
 func DefaultSettings() *Settings {
 	return &Settings{
-		Method:  Central,
+		Formula: Central,
 		Workers: runtime.GOMAXPROCS(0),
 	}
 }
@@ -57,23 +57,23 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 	if settings == nil {
 		settings = DefaultSettings()
 	}
-	step := settings.Method.Step
+	step := settings.Formula.Step
 	var deriv float64
-	method := settings.Method
+	formula := settings.Formula
 	if !settings.Concurrent {
-		for _, pt := range method.Stencil {
+		for _, pt := range formula.Stencil {
 			if settings.OriginKnown && pt.Loc == 0 {
 				deriv += pt.Coeff * settings.OriginValue
 				continue
 			}
 			deriv += pt.Coeff * f(x+step*pt.Loc)
 		}
-		return deriv / math.Pow(step, float64(method.Order))
+		return deriv / math.Pow(step, float64(formula.Order))
 	}
 
 	wg := &sync.WaitGroup{}
 	mux := &sync.Mutex{}
-	for _, pt := range method.Stencil {
+	for _, pt := range formula.Stencil {
 		if settings.OriginKnown && pt.Loc == 0 {
 			mux.Lock()
 			deriv += pt.Coeff * settings.OriginValue
@@ -90,7 +90,7 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 		}(pt)
 	}
 	wg.Wait()
-	return deriv / math.Pow(step, float64(method.Order))
+	return deriv / math.Pow(step, float64(formula.Order))
 }
 
 // Gradient estimates the gradient of the multivariate function f at the
@@ -109,13 +109,13 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	if settings == nil {
 		settings = DefaultSettings()
 	}
-	step := settings.Method.Step
+	step := settings.Formula.Step
 	if !settings.Concurrent {
 		xcopy := make([]float64, len(x)) // So that x is not modified during the call
 		copy(xcopy, x)
 		for i := range xcopy {
 			var deriv float64
-			for _, pt := range settings.Method.Stencil {
+			for _, pt := range settings.Formula.Stencil {
 				if settings.OriginKnown && pt.Loc == 0 {
 					deriv += pt.Coeff * settings.OriginValue
 					continue
@@ -124,13 +124,13 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 				deriv += pt.Coeff * f(xcopy)
 				xcopy[i] = x[i]
 			}
-			dst[i] = deriv / math.Pow(step, float64(settings.Method.Order))
+			dst[i] = deriv / math.Pow(step, float64(settings.Formula.Order))
 		}
 		return dst
 	}
 
 	nWorkers := settings.Workers
-	expect := len(settings.Method.Stencil) * len(x)
+	expect := len(settings.Formula.Stencil) * len(x)
 	if nWorkers > expect {
 		nWorkers = expect
 	}
@@ -162,7 +162,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	// Launch the distributor. Distributor sends the cases to be computed
 	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
 		for i := range x {
-			for _, pt := range settings.Method.Stencil {
+			for _, pt := range settings.Formula.Stencil {
 				if settings.OriginKnown && pt.Loc == 0 {
 					// Answer already known. Send the answer on the answer channel
 					ansChan <- fdrun{
@@ -189,7 +189,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		run := <-ansChan
 		dst[run.idx] += run.pt.Coeff * run.result
 	}
-	floats.Scale(1/math.Pow(step, float64(settings.Method.Order)), dst)
+	floats.Scale(1/math.Pow(step, float64(settings.Formula.Order)), dst)
 	return dst
 }
 
@@ -200,28 +200,28 @@ type fdrun struct {
 }
 
 // Forward represents a first-order forward difference.
-var Forward = Method{
+var Forward = Formula{
 	Stencil: []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
 	Order:   1,
 	Step:    1e-6,
 }
 
 // Backward represents a first-order backward difference.
-var Backward = Method{
+var Backward = Formula{
 	Stencil: []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
 	Order:   1,
 	Step:    1e-6,
 }
 
 // Central represents a first-order central difference.
-var Central = Method{
+var Central = Formula{
 	Stencil: []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
 	Order:   1,
 	Step:    1e-6,
 }
 
 // Central2nd represents a secord-order central difference.
-var Central2nd = Method{
+var Central2nd = Formula{
 	Stencil: []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
 	Order:   2,
 	Step:    1e-3,
diff --git a/fd/diff_test.go b/fd/diff_test.go
index d9f51c6b..580b03da 100644
--- a/fd/diff_test.go
+++ b/fd/diff_test.go
@@ -72,11 +72,11 @@ var testsSecond = []testPoint{
 	},
 }
 
-func testDerivative(t *testing.T, method Method, tol float64, tests []testPoint) {
+func testDerivative(t *testing.T, formula Formula, tol float64, tests []testPoint) {
 	for i, test := range tests {
 
 		settings := DefaultSettings()
-		settings.Method = method
+		settings.Formula = formula
 		ans := Derivative(test.f, test.loc, settings)
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("Case %v: ans mismatch serial: expected %v, found %v", i, test.ans, ans)
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 61d18e44..e1107b75 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -42,29 +42,29 @@ func (r Rosenbrock) FDf(x []float64, deriv []float64) (sum float64) {
 
 func TestGradient(t *testing.T) {
 	for i, test := range []struct {
-		nDim   int
-		tol    float64
-		method Method
+		nDim    int
+		tol     float64
+		formula Formula
 	}{
 		{
-			nDim:   2,
-			tol:    2e-4,
-			method: Forward,
+			nDim:    2,
+			tol:     2e-4,
+			formula: Forward,
 		},
 		{
-			nDim:   2,
-			tol:    1e-6,
-			method: Central,
+			nDim:    2,
+			tol:     1e-6,
+			formula: Central,
 		},
 		{
-			nDim:   40,
-			tol:    2e-4,
-			method: Forward,
+			nDim:    40,
+			tol:     2e-4,
+			formula: Forward,
 		},
 		{
-			nDim:   40,
-			tol:    1e-6,
-			method: Central,
+			nDim:    40,
+			tol:     1e-6,
+			formula: Central,
 		},
 	} {
 		x := make([]float64, test.nDim)
@@ -79,7 +79,7 @@ func TestGradient(t *testing.T) {
 		r.FDf(x, trueGradient)
 
 		settings := DefaultSettings()
-		settings.Method = test.method
+		settings.Formula = test.formula
 
 		// try with gradient nil
 		gradient := Gradient(nil, r.F, x, settings)

From 6709b8f0a9d4da3e9ff263b9f08de38be3af3d19 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 15:39:03 +0900
Subject: [PATCH 58/98] fd: remove Settings.Workers

Use runtime.NumCPU() instead
---
 fd/diff.go          | 18 ++++++++----------
 fd/gradient_test.go |  1 -
 2 files changed, 8 insertions(+), 11 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index ac2bf5fc..cf427a4b 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -35,8 +35,7 @@ type Formula struct {
 type Settings struct {
 	OriginKnown bool    // Flag that the value at the origin x is known
 	OriginValue float64 // Value at the origin (only used if OriginKnown is true)
-	Concurrent  bool    // Should the function calls be executed concurrently
-	Workers     int     // Maximum number of concurrent executions when evaluating concurrently
+	Concurrent  bool    // Should the function calls be executed concurrently.
 	Formula     Formula // Finite difference formula to use
 }
 
@@ -46,7 +45,6 @@ type Settings struct {
 func DefaultSettings() *Settings {
 	return &Settings{
 		Formula: Central,
-		Workers: runtime.GOMAXPROCS(0),
 	}
 }
 
@@ -129,19 +127,19 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		return dst
 	}
 
-	nWorkers := settings.Workers
-	expect := len(settings.Formula.Stencil) * len(x)
-	if nWorkers > expect {
-		nWorkers = expect
-	}
-
 	quit := make(chan struct{})
 	defer close(quit)
+
+	expect := len(settings.Formula.Stencil) * len(x)
 	sendChan := make(chan fdrun, expect)
 	ansChan := make(chan fdrun, expect)
 
 	// Launch workers. Workers receive an index and a step, and compute the answer
-	for i := 0; i < settings.Workers; i++ {
+	nWorkers := runtime.NumCPU()
+	if nWorkers > expect {
+		nWorkers = expect
+	}
+	for i := 0; i < nWorkers; i++ {
 		go func(sendChan <-chan fdrun, ansChan chan<- fdrun, quit <-chan struct{}) {
 			xcopy := make([]float64, len(x))
 			copy(xcopy, x)
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index e1107b75..3796ba4c 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -119,7 +119,6 @@ func TestGradient(t *testing.T) {
 		}
 		settings.Concurrent = true
 		settings.OriginKnown = false
-		settings.Workers = 1000
 		Gradient(gradient, r.F, x, settings)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)

From 989bdb3ddc1adb72a55d1586a0d2aea1b7dcb6e9 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 16:13:01 +0900
Subject: [PATCH 59/98] fd: add user-settable step size to Settings

---
 fd/diff.go | 34 ++++++++++++++++++++++++----------
 1 file changed, 24 insertions(+), 10 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index cf427a4b..8884a719 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -21,8 +21,7 @@ type Point struct {
 
 // Formula represents a finite difference formula that approximates
 // the derivative of order k of a function f at x as
-//  d^k f(x) ≅ (1 / h^k) * \sum_i Coeff_i * f(x + h * Loc_i),
-// where h is a small positive step.
+//  d^k f(x) ≅ (1 / Step^k) * \sum_i Coeff_i * f(x + Step * Loc_i).
 type Formula struct {
 	// Stencil is the set of sampling Points which are used to estimate the
 	// derivative. The locations will be scaled by Step and are relative to x.
@@ -33,10 +32,16 @@ type Formula struct {
 
 // Settings is the settings structure for computing finite differences.
 type Settings struct {
+	// Formula is the finite difference formula used
+	// for approximating the derivative.
+	Formula Formula
+	// Step is the distance between points of the stencil.
+	// If equal to 0, formula's default step will be used.
+	Step float64
+
 	OriginKnown bool    // Flag that the value at the origin x is known
 	OriginValue float64 // Value at the origin (only used if OriginKnown is true)
 	Concurrent  bool    // Should the function calls be executed concurrently.
-	Formula     Formula // Finite difference formula to use
 }
 
 // DefaultSettings is a basic set of settings for computing finite differences.
@@ -55,9 +60,13 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 	if settings == nil {
 		settings = DefaultSettings()
 	}
-	step := settings.Formula.Step
-	var deriv float64
 	formula := settings.Formula
+	step := settings.Step
+	if step == 0 {
+		step = formula.Step
+	}
+
+	var deriv float64
 	if !settings.Concurrent {
 		for _, pt := range formula.Stencil {
 			if settings.OriginKnown && pt.Loc == 0 {
@@ -107,13 +116,18 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	if settings == nil {
 		settings = DefaultSettings()
 	}
-	step := settings.Formula.Step
+	formula := settings.Formula
+	step := settings.Step
+	if step == 0 {
+		step = formula.Step
+	}
+
 	if !settings.Concurrent {
 		xcopy := make([]float64, len(x)) // So that x is not modified during the call
 		copy(xcopy, x)
 		for i := range xcopy {
 			var deriv float64
-			for _, pt := range settings.Formula.Stencil {
+			for _, pt := range formula.Stencil {
 				if settings.OriginKnown && pt.Loc == 0 {
 					deriv += pt.Coeff * settings.OriginValue
 					continue
@@ -122,7 +136,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 				deriv += pt.Coeff * f(xcopy)
 				xcopy[i] = x[i]
 			}
-			dst[i] = deriv / math.Pow(step, float64(settings.Formula.Order))
+			dst[i] = deriv / math.Pow(step, float64(formula.Order))
 		}
 		return dst
 	}
@@ -160,7 +174,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	// Launch the distributor. Distributor sends the cases to be computed
 	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
 		for i := range x {
-			for _, pt := range settings.Formula.Stencil {
+			for _, pt := range formula.Stencil {
 				if settings.OriginKnown && pt.Loc == 0 {
 					// Answer already known. Send the answer on the answer channel
 					ansChan <- fdrun{
@@ -187,7 +201,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		run := <-ansChan
 		dst[run.idx] += run.pt.Coeff * run.result
 	}
-	floats.Scale(1/math.Pow(step, float64(settings.Formula.Order)), dst)
+	floats.Scale(1/math.Pow(step, float64(formula.Order)), dst)
 	return dst
 }
 

From e56f8dc55f9ae3f7f631945a2e705815810ef833 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 16:17:03 +0900
Subject: [PATCH 60/98] fd: panic in Gradient if the derivative order is not 1

---
 fd/diff.go | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 8884a719..36da442f 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -105,7 +105,8 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 // otherwise a new slice will be allocated and returned. Finite difference
 // kernel and other options are specified by settings. If settings is nil,
 // default settings will be used.
-// Gradient panics if the length of dst and x is not equal.
+// Gradient panics if the length of dst and x is not equal, or if the
+// derivative order of the formula is not 1.
 func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
 	if dst == nil {
 		dst = make([]float64, len(x))
@@ -116,6 +117,9 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	if settings == nil {
 		settings = DefaultSettings()
 	}
+	if settings.Formula.Order != 1 {
+		panic("fd: invalid derivative order")
+	}
 	formula := settings.Formula
 	step := settings.Step
 	if step == 0 {
@@ -136,7 +140,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 				deriv += pt.Coeff * f(xcopy)
 				xcopy[i] = x[i]
 			}
-			dst[i] = deriv / math.Pow(step, float64(formula.Order))
+			dst[i] = deriv / step
 		}
 		return dst
 	}
@@ -201,7 +205,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		run := <-ansChan
 		dst[run.idx] += run.pt.Coeff * run.result
 	}
-	floats.Scale(1/math.Pow(step, float64(formula.Order)), dst)
+	floats.Scale(1/step, dst)
 	return dst
 }
 

From bd9fca8b99f2fd4b32082e265e6d1274ad2d53e2 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 16:52:04 +0900
Subject: [PATCH 61/98] fd: use Central formula if Stencil is nil

---
 fd/diff.go | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/fd/diff.go b/fd/diff.go
index 36da442f..39cd2af3 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -61,6 +61,9 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 		settings = DefaultSettings()
 	}
 	formula := settings.Formula
+	if formula.Stencil == nil {
+		formula = Central
+	}
 	step := settings.Step
 	if step == 0 {
 		step = formula.Step
@@ -121,6 +124,9 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		panic("fd: invalid derivative order")
 	}
 	formula := settings.Formula
+	if formula.Stencil == nil {
+		formula = Central
+	}
 	step := settings.Step
 	if step == 0 {
 		step = formula.Step

From ecece58210b25d7a4c0730e614f1cef4f2148508 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 16 Mar 2016 17:09:09 +0900
Subject: [PATCH 62/98] fd: add example for Derivative

---
 fd/example_test.go | 47 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 fd/example_test.go

diff --git a/fd/example_test.go b/fd/example_test.go
new file mode 100644
index 00000000..75043b40
--- /dev/null
+++ b/fd/example_test.go
@@ -0,0 +1,47 @@
+// Copyright ©2016 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 fd_test
+
+import (
+	"fmt"
+	"math"
+
+	"github.com/gonum/diff/fd"
+)
+
+func ExampleDerivative() {
+	f := func(x float64) float64 {
+		return math.Sin(x)
+	}
+	// Compute the first derivative of f at 0 using the default settings.
+	fmt.Println("f'(0) ≅", fd.Derivative(f, 0, nil))
+	// Compute the first derivative of f at 0 using the forward approximation
+	// with a custom step size.
+	df := fd.Derivative(f, 0, &fd.Settings{
+		Formula: fd.Forward,
+		Step:    1e-8,
+	})
+	fmt.Println("f'(0) ≅", df)
+
+	f = func(x float64) float64 {
+		return math.Pow(math.Cos(x), 3)
+	}
+	// Compute the second derivative of f at 0 using the centered
+	// approximation, a custom step size, concurrent evaluation, and a known
+	// function value at x.
+	df = fd.Derivative(f, 0, &fd.Settings{
+		Formula:     fd.Central2nd,
+		Step:        1e-4,
+		Concurrent:  2,
+		OriginKnown: true,
+		OriginValue: f(0),
+	})
+	fmt.Println("f''(0) ≅", df)
+
+	// Output:
+	// f'(0) ≅ 0.9999999999998334
+	// f'(0) ≅ 1
+	// f''(0) ≅ -2.999999981767587
+}

From 4f9ba9d95f16c12fcf1ca8aaa10d71c815d39ff9 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Thu, 17 Mar 2016 09:49:14 +0900
Subject: [PATCH 63/98] fd: update comment for Formula

---
 fd/diff.go | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 39cd2af3..6c245a1f 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -19,8 +19,8 @@ type Point struct {
 	Coeff float64
 }
 
-// Formula represents a finite difference formula that approximates
-// the derivative of order k of a function f at x as
+// Formula represents a finite difference formula on a regularly spaced grid
+// that approximates the derivative of order k of a function f at x as
 //  d^k f(x) ≅ (1 / Step^k) * \sum_i Coeff_i * f(x + Step * Loc_i).
 type Formula struct {
 	// Stencil is the set of sampling Points which are used to estimate the

From 27a027e9c96bb1857305a9b652593b7e115c17a9 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Thu, 17 Mar 2016 09:51:44 +0900
Subject: [PATCH 64/98] fd: finish comments with dot

---
 fd/diff.go | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 6c245a1f..d517c5b7 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -39,9 +39,10 @@ type Settings struct {
 	// If equal to 0, formula's default step will be used.
 	Step float64
 
-	OriginKnown bool    // Flag that the value at the origin x is known
-	OriginValue float64 // Value at the origin (only used if OriginKnown is true)
-	Concurrent  bool    // Should the function calls be executed concurrently.
+	OriginKnown bool    // Flag that the value at the origin x is known.
+	OriginValue float64 // Value at the origin (only used if OriginKnown is true).
+
+	Concurrent bool // Should the function calls be executed concurrently.
 }
 
 // DefaultSettings is a basic set of settings for computing finite differences.
@@ -133,7 +134,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	}
 
 	if !settings.Concurrent {
-		xcopy := make([]float64, len(x)) // So that x is not modified during the call
+		xcopy := make([]float64, len(x)) // So that x is not modified during the call.
 		copy(xcopy, x)
 		for i := range xcopy {
 			var deriv float64
@@ -158,7 +159,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	sendChan := make(chan fdrun, expect)
 	ansChan := make(chan fdrun, expect)
 
-	// Launch workers. Workers receive an index and a step, and compute the answer
+	// Launch workers. Workers receive an index and a step, and compute the answer.
 	nWorkers := runtime.NumCPU()
 	if nWorkers > expect {
 		nWorkers = expect
@@ -181,12 +182,12 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		}(sendChan, ansChan, quit)
 	}
 
-	// Launch the distributor. Distributor sends the cases to be computed
+	// Launch the distributor. Distributor sends the cases to be computed.
 	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
 		for i := range x {
 			for _, pt := range formula.Stencil {
 				if settings.OriginKnown && pt.Loc == 0 {
-					// Answer already known. Send the answer on the answer channel
+					// Answer already known. Send the answer on the answer channel.
 					ansChan <- fdrun{
 						idx:    i,
 						pt:     pt,
@@ -194,7 +195,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 					}
 					continue
 				}
-				// Answer not known, send the answer to be computed
+				// Answer not known, send the answer to be computed.
 				sendChan <- fdrun{
 					idx: i,
 					pt:  pt,
@@ -206,7 +207,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	for i := range dst {
 		dst[i] = 0
 	}
-	// Read in all of the results
+	// Read in all of the results.
 	for i := 0; i < expect; i++ {
 		run := <-ansChan
 		dst[run.idx] += run.pt.Coeff * run.result

From af9cab08caebfb9266be242b35d14d88e21bf579 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Thu, 17 Mar 2016 11:31:51 +0900
Subject: [PATCH 65/98] fd: add tests for zero value of Settings

---
 fd/diff_test.go     | 8 +++++++-
 fd/gradient_test.go | 8 ++++++++
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/fd/diff_test.go b/fd/diff_test.go
index 580b03da..58029fac 100644
--- a/fd/diff_test.go
+++ b/fd/diff_test.go
@@ -120,7 +120,8 @@ func TestCentralSecond(t *testing.T) {
 	testDerivative(t, Central2nd, 1e-3, testsSecond)
 }
 
-// TestDerivativeDefault checks that the derivative works when settings is set to nil.
+// TestDerivativeDefault checks that the derivative works when settings is nil
+// or zero value.
 func TestDerivativeDefault(t *testing.T) {
 	tol := 1e-6
 	for i, test := range testsFirst {
@@ -128,5 +129,10 @@ func TestDerivativeDefault(t *testing.T) {
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("Case %v: ans mismatch default: expected %v, found %v", i, test.ans, ans)
 		}
+
+		ans = Derivative(test.f, test.loc, &Settings{})
+		if math.Abs(test.ans-ans) > tol {
+			t.Errorf("Case %v: ans mismatch zero value: expected %v, found %v", i, test.ans, ans)
+		}
 	}
 }
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 3796ba4c..3f5a020a 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -147,6 +147,14 @@ func TestGradient(t *testing.T) {
 			t.Errorf("Case %v: gradient mismatch with default settings. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 
+		// With zero settings
+		for i := range gradient {
+			gradient[i] = rand.Float64()
+		}
+		Gradient(gradient, r.F, x, &Settings{})
+		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
+			t.Errorf("Case %v: gradient mismatch with zero settings. Want: %v, Got: %v.", i, trueGradient, gradient)
+		}
 	}
 }
 

From 95f9eb6b08812d21613eedc9ce5a5128cf63b902 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Thu, 17 Mar 2016 11:39:25 +0900
Subject: [PATCH 66/98] fd: test formula's Order after setting defaults in
 Gradient

---
 fd/diff.go | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index d517c5b7..5ec10791 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -121,13 +121,13 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	if settings == nil {
 		settings = DefaultSettings()
 	}
-	if settings.Formula.Order != 1 {
-		panic("fd: invalid derivative order")
-	}
 	formula := settings.Formula
 	if formula.Stencil == nil {
 		formula = Central
 	}
+	if formula.Order != 1 {
+		panic("fd: invalid derivative order")
+	}
 	step := settings.Step
 	if step == 0 {
 		step = formula.Step

From 9376339a029e32cc49624acd1521ae2b87bddd6c Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Thu, 17 Mar 2016 14:47:03 +0900
Subject: [PATCH 67/98] fd: remove DefaultSettings

---
 fd/diff.go          | 19 ++++++-----------
 fd/diff_test.go     | 29 ++++++++++++++++----------
 fd/example_test.go  |  2 +-
 fd/gradient_test.go | 51 ++++++++++++++++++++++++++-------------------
 4 files changed, 54 insertions(+), 47 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 5ec10791..89ec8471 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -34,6 +34,7 @@ type Formula struct {
 type Settings struct {
 	// Formula is the finite difference formula used
 	// for approximating the derivative.
+	// Zero value will default to the Central formula.
 	Formula Formula
 	// Step is the distance between points of the stencil.
 	// If equal to 0, formula's default step will be used.
@@ -45,21 +46,13 @@ type Settings struct {
 	Concurrent bool // Should the function calls be executed concurrently.
 }
 
-// DefaultSettings is a basic set of settings for computing finite differences.
-// Computes a central difference approximation for the first derivative
-// of the function.
-func DefaultSettings() *Settings {
-	return &Settings{
-		Formula: Central,
-	}
-}
-
 // Derivative estimates the derivative of the function f at the given location.
-// The order of the derivative, sample locations, and other options are
-// specified by settings. If settings is nil, default settings will be used.
+// The finite difference formula, the step size, and other options are
+// specified by settings. If settings is nil, the first derivative will be
+// estimated using the Central formula and a suitable step size.
 func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
 	if settings == nil {
-		settings = DefaultSettings()
+		settings = &Settings{}
 	}
 	formula := settings.Formula
 	if formula.Stencil == nil {
@@ -119,7 +112,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		panic("fd: slice length mismatch")
 	}
 	if settings == nil {
-		settings = DefaultSettings()
+		settings = &Settings{}
 	}
 	formula := settings.Formula
 	if formula.Stencil == nil {
diff --git a/fd/diff_test.go b/fd/diff_test.go
index 58029fac..35d4f0f3 100644
--- a/fd/diff_test.go
+++ b/fd/diff_test.go
@@ -75,29 +75,36 @@ var testsSecond = []testPoint{
 func testDerivative(t *testing.T, formula Formula, tol float64, tests []testPoint) {
 	for i, test := range tests {
 
-		settings := DefaultSettings()
-		settings.Formula = formula
-		ans := Derivative(test.f, test.loc, settings)
+		ans := Derivative(test.f, test.loc, &Settings{
+			Formula: formula,
+		})
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("Case %v: ans mismatch serial: expected %v, found %v", i, test.ans, ans)
 		}
 
-		settings.OriginKnown = true
-		settings.OriginValue = test.fofx
-		ans = Derivative(test.f, test.loc, settings)
+		ans = Derivative(test.f, test.loc, &Settings{
+			Formula:     formula,
+			OriginKnown: true,
+			OriginValue: test.fofx,
+		})
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("Case %v: ans mismatch serial origin known: expected %v, found %v", i, test.ans, ans)
 		}
 
-		settings.OriginKnown = false
-		settings.Concurrent = true
-		ans = Derivative(test.f, test.loc, settings)
+		ans = Derivative(test.f, test.loc, &Settings{
+			Formula:    formula,
+			Concurrent: true,
+		})
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("Case %v: ans mismatch concurrent: expected %v, found %v", i, test.ans, ans)
 		}
 
-		settings.OriginKnown = true
-		ans = Derivative(test.f, test.loc, settings)
+		ans = Derivative(test.f, test.loc, &Settings{
+			Formula:     formula,
+			OriginKnown: true,
+			OriginValue: test.fofx,
+			Concurrent:  true,
+		})
 		if math.Abs(test.ans-ans) > tol {
 			t.Errorf("Case %v: ans mismatch concurrent: expected %v, found %v", i, test.ans, ans)
 		}
diff --git a/fd/example_test.go b/fd/example_test.go
index 75043b40..e46538ba 100644
--- a/fd/example_test.go
+++ b/fd/example_test.go
@@ -34,7 +34,7 @@ func ExampleDerivative() {
 	df = fd.Derivative(f, 0, &fd.Settings{
 		Formula:     fd.Central2nd,
 		Step:        1e-4,
-		Concurrent:  2,
+		Concurrent:  true,
 		OriginKnown: true,
 		OriginValue: f(0),
 	})
diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index 3f5a020a..d065cb3f 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -78,11 +78,10 @@ func TestGradient(t *testing.T) {
 		trueGradient := make([]float64, len(x))
 		r.FDf(x, trueGradient)
 
-		settings := DefaultSettings()
-		settings.Formula = test.formula
-
-		// try with gradient nil
-		gradient := Gradient(nil, r.F, x, settings)
+		// Try with gradient nil.
+		gradient := Gradient(nil, r.F, x, &Settings{
+			Formula: test.formula,
+		})
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch in serial with nil. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
@@ -90,11 +89,13 @@ func TestGradient(t *testing.T) {
 			t.Errorf("Case %v: x modified during call to gradient in serial with nil.", i)
 		}
 
+		// Try with provided gradient.
 		for i := range gradient {
 			gradient[i] = rand.Float64()
 		}
-
-		Gradient(gradient, r.F, x, settings)
+		Gradient(gradient, r.F, x, &Settings{
+			Formula: test.formula,
+		})
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
@@ -102,24 +103,27 @@ func TestGradient(t *testing.T) {
 			t.Errorf("Case %v: x modified during call to gradient in serial with non-nil.", i)
 		}
 
-		// Try with known value
+		// Try with known value.
 		for i := range gradient {
 			gradient[i] = rand.Float64()
 		}
-		settings.OriginKnown = true
-		settings.OriginValue = r.F(x)
-		Gradient(gradient, r.F, x, settings)
+		Gradient(gradient, r.F, x, &Settings{
+			Formula:     test.formula,
+			OriginKnown: true,
+			OriginValue: r.F(x),
+		})
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with known origin in serial. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 
-		// Concurrently
+		// Try with concurrent evaluation.
 		for i := range gradient {
 			gradient[i] = rand.Float64()
 		}
-		settings.Concurrent = true
-		settings.OriginKnown = false
-		Gradient(gradient, r.F, x, settings)
+		Gradient(gradient, r.F, x, &Settings{
+			Formula:    test.formula,
+			Concurrent: true,
+		})
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with unknown origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
@@ -127,27 +131,30 @@ func TestGradient(t *testing.T) {
 			t.Errorf("Case %v: x modified during call to gradient in parallel", i)
 		}
 
-		// Concurrently with origin known
+		// Try with concurrent evaluation with origin known.
 		for i := range gradient {
 			gradient[i] = rand.Float64()
 		}
-		settings.OriginKnown = true
-		Gradient(gradient, r.F, x, settings)
+		Gradient(gradient, r.F, x, &Settings{
+			Formula:     test.formula,
+			Concurrent:  true,
+			OriginKnown: true,
+			OriginValue: r.F(x),
+		})
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with known origin in parallel. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 
-		// With default settings
+		// Try with nil settings.
 		for i := range gradient {
 			gradient[i] = rand.Float64()
 		}
-		settings = nil
-		Gradient(gradient, r.F, x, settings)
+		Gradient(gradient, r.F, x, nil)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
 			t.Errorf("Case %v: gradient mismatch with default settings. Want: %v, Got: %v.", i, trueGradient, gradient)
 		}
 
-		// With zero settings
+		// Try with zero-valued settings.
 		for i := range gradient {
 			gradient[i] = rand.Float64()
 		}

From 236271b55bbae6527b1edd1df31424d41ce9cf67 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Thu, 17 Mar 2016 14:49:37 +0900
Subject: [PATCH 68/98] =?UTF-8?q?fd:=20replace=20=E2=89=85=20with=20?=
 =?UTF-8?q?=E2=89=88?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 fd/diff.go         |  2 +-
 fd/example_test.go | 12 ++++++------
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 89ec8471..69acc9cd 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -21,7 +21,7 @@ type Point struct {
 
 // Formula represents a finite difference formula on a regularly spaced grid
 // that approximates the derivative of order k of a function f at x as
-//  d^k f(x) ≅ (1 / Step^k) * \sum_i Coeff_i * f(x + Step * Loc_i).
+//  d^k f(x) ≈ (1 / Step^k) * \sum_i Coeff_i * f(x + Step * Loc_i).
 type Formula struct {
 	// Stencil is the set of sampling Points which are used to estimate the
 	// derivative. The locations will be scaled by Step and are relative to x.
diff --git a/fd/example_test.go b/fd/example_test.go
index e46538ba..b0d974ab 100644
--- a/fd/example_test.go
+++ b/fd/example_test.go
@@ -16,14 +16,14 @@ func ExampleDerivative() {
 		return math.Sin(x)
 	}
 	// Compute the first derivative of f at 0 using the default settings.
-	fmt.Println("f'(0) ≅", fd.Derivative(f, 0, nil))
+	fmt.Println("f'(0) ≈", fd.Derivative(f, 0, nil))
 	// Compute the first derivative of f at 0 using the forward approximation
 	// with a custom step size.
 	df := fd.Derivative(f, 0, &fd.Settings{
 		Formula: fd.Forward,
 		Step:    1e-8,
 	})
-	fmt.Println("f'(0) ≅", df)
+	fmt.Println("f'(0) ≈", df)
 
 	f = func(x float64) float64 {
 		return math.Pow(math.Cos(x), 3)
@@ -38,10 +38,10 @@ func ExampleDerivative() {
 		OriginKnown: true,
 		OriginValue: f(0),
 	})
-	fmt.Println("f''(0) ≅", df)
+	fmt.Println("f''(0) ≈", df)
 
 	// Output:
-	// f'(0) ≅ 0.9999999999998334
-	// f'(0) ≅ 1
-	// f''(0) ≅ -2.999999981767587
+	// f'(0) ≈ 0.9999999999998334
+	// f'(0) ≈ 1
+	// f''(0) ≈ -2.999999981767587
 }

From 063ade24cbf17dd393b2afa437d86b4607e9284f Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Tue, 22 Mar 2016 17:25:41 +0900
Subject: [PATCH 69/98] fd: change wording in Derivative comment

---
 fd/diff.go | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fd/diff.go b/fd/diff.go
index 69acc9cd..f62f3054 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -49,7 +49,7 @@ type Settings struct {
 // Derivative estimates the derivative of the function f at the given location.
 // The finite difference formula, the step size, and other options are
 // specified by settings. If settings is nil, the first derivative will be
-// estimated using the Central formula and a suitable step size.
+// estimated using the Central formula and a default step size.
 func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
 	if settings == nil {
 		settings = &Settings{}

From c25ba45005e7e1954b10579110f2b89cf9f7a333 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Tue, 22 Mar 2016 17:32:13 +0900
Subject: [PATCH 70/98] fd: check properly whether Formula is valid

---
 fd/diff.go | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index f62f3054..d22b36b5 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -30,6 +30,10 @@ type Formula struct {
 	Step    float64 // Default step size for the formula.
 }
 
+func (f Formula) isZero() bool {
+	return f.Stencil == nil && f.Order == 0 && f.Step == 0
+}
+
 // Settings is the settings structure for computing finite differences.
 type Settings struct {
 	// Formula is the finite difference formula used
@@ -55,9 +59,12 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 		settings = &Settings{}
 	}
 	formula := settings.Formula
-	if formula.Stencil == nil {
+	if formula.isZero() {
 		formula = Central
 	}
+	if formula.Order == 0 || formula.Stencil == nil || formula.Step == 0 {
+		panic("fd: bad formula")
+	}
 	step := settings.Step
 	if step == 0 {
 		step = formula.Step
@@ -115,9 +122,12 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		settings = &Settings{}
 	}
 	formula := settings.Formula
-	if formula.Stencil == nil {
+	if formula.isZero() {
 		formula = Central
 	}
+	if formula.Order == 0 || formula.Stencil == nil || formula.Step == 0 {
+		panic("fd: bad formula")
+	}
 	if formula.Order != 1 {
 		panic("fd: invalid derivative order")
 	}

From 8f1a0ce5965239eae0b1050e2f6966ee7322b084 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 23 Mar 2016 11:36:03 +0900
Subject: [PATCH 71/98] fd: change condition for serial evaluation in Gradient

---
 fd/diff.go | 23 ++++++++++++++---------
 1 file changed, 14 insertions(+), 9 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index d22b36b5..f560f14f 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -121,6 +121,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	if settings == nil {
 		settings = &Settings{}
 	}
+
 	formula := settings.Formula
 	if formula.isZero() {
 		formula = Central
@@ -131,12 +132,22 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	if formula.Order != 1 {
 		panic("fd: invalid derivative order")
 	}
+
 	step := settings.Step
 	if step == 0 {
 		step = formula.Step
 	}
 
-	if !settings.Concurrent {
+	expect := len(formula.Stencil) * len(x)
+	nWorkers := 1
+	if settings.Concurrent {
+		nWorkers = runtime.GOMAXPROCS(0)
+		if nWorkers > expect {
+			nWorkers = expect
+		}
+	}
+
+	if nWorkers == 1 {
 		xcopy := make([]float64, len(x)) // So that x is not modified during the call.
 		copy(xcopy, x)
 		for i := range xcopy {
@@ -155,18 +166,12 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		return dst
 	}
 
+	sendChan := make(chan fdrun, expect)
+	ansChan := make(chan fdrun, expect)
 	quit := make(chan struct{})
 	defer close(quit)
 
-	expect := len(settings.Formula.Stencil) * len(x)
-	sendChan := make(chan fdrun, expect)
-	ansChan := make(chan fdrun, expect)
-
 	// Launch workers. Workers receive an index and a step, and compute the answer.
-	nWorkers := runtime.NumCPU()
-	if nWorkers > expect {
-		nWorkers = expect
-	}
 	for i := 0; i < nWorkers; i++ {
 		go func(sendChan <-chan fdrun, ansChan chan<- fdrun, quit <-chan struct{}) {
 			xcopy := make([]float64, len(x))

From 55ff30754425e33f24aeddfe4956e56999e9b5cd Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 23 Mar 2016 11:49:55 +0900
Subject: [PATCH 72/98] fd: rename Formula.Order to Formula.Derivative

---
 fd/diff.go | 42 +++++++++++++++++++++---------------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index f560f14f..80c2544e 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -25,13 +25,13 @@ type Point struct {
 type Formula struct {
 	// Stencil is the set of sampling Points which are used to estimate the
 	// derivative. The locations will be scaled by Step and are relative to x.
-	Stencil []Point
-	Order   int     // The order of the approximated derivative.
-	Step    float64 // Default step size for the formula.
+	Stencil    []Point
+	Derivative int     // The order of the approximated derivative.
+	Step       float64 // Default step size for the formula.
 }
 
 func (f Formula) isZero() bool {
-	return f.Stencil == nil && f.Order == 0 && f.Step == 0
+	return f.Stencil == nil && f.Derivative == 0 && f.Step == 0
 }
 
 // Settings is the settings structure for computing finite differences.
@@ -62,7 +62,7 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 	if formula.isZero() {
 		formula = Central
 	}
-	if formula.Order == 0 || formula.Stencil == nil || formula.Step == 0 {
+	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
 		panic("fd: bad formula")
 	}
 	step := settings.Step
@@ -79,7 +79,7 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 			}
 			deriv += pt.Coeff * f(x+step*pt.Loc)
 		}
-		return deriv / math.Pow(step, float64(formula.Order))
+		return deriv / math.Pow(step, float64(formula.Derivative))
 	}
 
 	wg := &sync.WaitGroup{}
@@ -101,7 +101,7 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 		}(pt)
 	}
 	wg.Wait()
-	return deriv / math.Pow(step, float64(formula.Order))
+	return deriv / math.Pow(step, float64(formula.Derivative))
 }
 
 // Gradient estimates the gradient of the multivariate function f at the
@@ -126,10 +126,10 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	if formula.isZero() {
 		formula = Central
 	}
-	if formula.Order == 0 || formula.Stencil == nil || formula.Step == 0 {
+	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
 		panic("fd: bad formula")
 	}
-	if formula.Order != 1 {
+	if formula.Derivative != 1 {
 		panic("fd: invalid derivative order")
 	}
 
@@ -232,28 +232,28 @@ type fdrun struct {
 
 // Forward represents a first-order forward difference.
 var Forward = Formula{
-	Stencil: []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
-	Order:   1,
-	Step:    1e-6,
+	Stencil:    []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
+	Derivative: 1,
+	Step:       1e-6,
 }
 
 // Backward represents a first-order backward difference.
 var Backward = Formula{
-	Stencil: []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
-	Order:   1,
-	Step:    1e-6,
+	Stencil:    []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
+	Derivative: 1,
+	Step:       1e-6,
 }
 
 // Central represents a first-order central difference.
 var Central = Formula{
-	Stencil: []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
-	Order:   1,
-	Step:    1e-6,
+	Stencil:    []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
+	Derivative: 1,
+	Step:       1e-6,
 }
 
 // Central2nd represents a secord-order central difference.
 var Central2nd = Formula{
-	Stencil: []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
-	Order:   2,
-	Step:    1e-3,
+	Stencil:    []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
+	Derivative: 2,
+	Step:       1e-3,
 }

From 53278edd298dd11dcdae73d14d1765dd50bd02d1 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 23 Mar 2016 13:24:14 +0900
Subject: [PATCH 73/98] fd: change condition for serial evaluation in
 Derivative

---
 fd/diff.go | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fd/diff.go b/fd/diff.go
index 80c2544e..d8c7c4e8 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -71,7 +71,7 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 	}
 
 	var deriv float64
-	if !settings.Concurrent {
+	if !settings.Concurrent || runtime.GOMAXPROCS(0) == 1 {
 		for _, pt := range formula.Stencil {
 			if settings.OriginKnown && pt.Loc == 0 {
 				deriv += pt.Coeff * settings.OriginValue

From 85a23c46f700184ad068b9ad9db3890789a4efe6 Mon Sep 17 00:00:00 2001
From: btracey <tracey.brendan@gmail.com>
Date: Wed, 23 Mar 2016 09:00:00 -0600
Subject: [PATCH 74/98] Protect against user function modifying input slice.

Fixes #24.
---
 fd/diff.go | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index d8c7c4e8..34d12151 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -149,7 +149,6 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 
 	if nWorkers == 1 {
 		xcopy := make([]float64, len(x)) // So that x is not modified during the call.
-		copy(xcopy, x)
 		for i := range xcopy {
 			var deriv float64
 			for _, pt := range formula.Stencil {
@@ -157,9 +156,9 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 					deriv += pt.Coeff * settings.OriginValue
 					continue
 				}
+				copy(xcopy, x)
 				xcopy[i] += pt.Loc * step
 				deriv += pt.Coeff * f(xcopy)
-				xcopy[i] = x[i]
 			}
 			dst[i] = deriv / step
 		}
@@ -175,15 +174,14 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	for i := 0; i < nWorkers; i++ {
 		go func(sendChan <-chan fdrun, ansChan chan<- fdrun, quit <-chan struct{}) {
 			xcopy := make([]float64, len(x))
-			copy(xcopy, x)
 			for {
 				select {
 				case <-quit:
 					return
 				case run := <-sendChan:
+					copy(xcopy, x)
 					xcopy[run.idx] += run.pt.Loc * step
 					run.result = f(xcopy)
-					xcopy[run.idx] = x[run.idx]
 					ansChan <- run
 				}
 			}

From 8a35fbb281254d58d7c00548e357728c4b1d4b62 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Fri, 25 Mar 2016 12:12:14 +0900
Subject: [PATCH 75/98] fd: evaluate f only once at x in Gradient

---
 fd/diff.go | 23 ++++++++++++++++++-----
 1 file changed, 18 insertions(+), 5 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 34d12151..2980d96e 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -147,13 +147,26 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 		}
 	}
 
+	var hasOrigin bool
+	for _, pt := range formula.Stencil {
+		if pt.Loc == 0 {
+			hasOrigin = true
+			break
+		}
+	}
+	xcopy := make([]float64, len(x)) // So that x is not modified during the call.
+	originValue := settings.OriginValue
+	if hasOrigin && !settings.OriginKnown {
+		copy(xcopy, x)
+		originValue = f(xcopy)
+	}
+
 	if nWorkers == 1 {
-		xcopy := make([]float64, len(x)) // So that x is not modified during the call.
 		for i := range xcopy {
 			var deriv float64
 			for _, pt := range formula.Stencil {
-				if settings.OriginKnown && pt.Loc == 0 {
-					deriv += pt.Coeff * settings.OriginValue
+				if pt.Loc == 0 {
+					deriv += pt.Coeff * originValue
 					continue
 				}
 				copy(xcopy, x)
@@ -192,12 +205,12 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
 		for i := range x {
 			for _, pt := range formula.Stencil {
-				if settings.OriginKnown && pt.Loc == 0 {
+				if pt.Loc == 0 {
 					// Answer already known. Send the answer on the answer channel.
 					ansChan <- fdrun{
 						idx:    i,
 						pt:     pt,
-						result: settings.OriginValue,
+						result: originValue,
 					}
 					continue
 				}

From 401825754083cd3f49e239cec3b3e4354a0cc5f0 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Fri, 25 Mar 2016 15:00:08 +0900
Subject: [PATCH 76/98] fd: change default step in predefined Formulas

---
 fd/diff.go | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 2980d96e..83cfcfa8 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -245,26 +245,26 @@ type fdrun struct {
 var Forward = Formula{
 	Stencil:    []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
 	Derivative: 1,
-	Step:       1e-6,
+	Step:       2e-8,
 }
 
 // Backward represents a first-order backward difference.
 var Backward = Formula{
 	Stencil:    []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
 	Derivative: 1,
-	Step:       1e-6,
+	Step:       2e-8,
 }
 
 // Central represents a first-order central difference.
 var Central = Formula{
 	Stencil:    []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
 	Derivative: 1,
-	Step:       1e-6,
+	Step:       6e-6,
 }
 
 // Central2nd represents a secord-order central difference.
 var Central2nd = Formula{
 	Stencil:    []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
 	Derivative: 2,
-	Step:       1e-3,
+	Step:       1e-4,
 }

From 3824ecf96b80efbaddaf4cbdcfb1f067f87faee6 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Fri, 25 Mar 2016 15:04:13 +0900
Subject: [PATCH 77/98] fd: adapt example output due to changed default step

---
 fd/example_test.go | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fd/example_test.go b/fd/example_test.go
index b0d974ab..f5a01d6a 100644
--- a/fd/example_test.go
+++ b/fd/example_test.go
@@ -41,7 +41,7 @@ func ExampleDerivative() {
 	fmt.Println("f''(0) ≈", df)
 
 	// Output:
-	// f'(0) ≈ 0.9999999999998334
+	// f'(0) ≈ 0.999999999994
 	// f'(0) ≈ 1
 	// f''(0) ≈ -2.999999981767587
 }

From 2a2be36611d7133f3bd9e3bcb543dc3bea6888d8 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Fri, 25 Mar 2016 15:04:52 +0900
Subject: [PATCH 78/98] fd: update docs for predefined Formulas

---
 fd/diff.go | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 83cfcfa8..64aadb38 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -241,28 +241,32 @@ type fdrun struct {
 	result float64
 }
 
-// Forward represents a first-order forward difference.
+// Forward represents a first-order accurate forward approximation
+// to the first derivative.
 var Forward = Formula{
 	Stencil:    []Point{{Loc: 0, Coeff: -1}, {Loc: 1, Coeff: 1}},
 	Derivative: 1,
 	Step:       2e-8,
 }
 
-// Backward represents a first-order backward difference.
+// Backward represents a first-order accurate backward approximation
+// to the first derivative.
 var Backward = Formula{
 	Stencil:    []Point{{Loc: -1, Coeff: -1}, {Loc: 0, Coeff: 1}},
 	Derivative: 1,
 	Step:       2e-8,
 }
 
-// Central represents a first-order central difference.
+// Central represents a second-order accurate centered approximation
+// to the first derivative.
 var Central = Formula{
 	Stencil:    []Point{{Loc: -1, Coeff: -0.5}, {Loc: 1, Coeff: 0.5}},
 	Derivative: 1,
 	Step:       6e-6,
 }
 
-// Central2nd represents a secord-order central difference.
+// Central2nd represents a secord-order accurate centered approximation
+// to the second derivative.
 var Central2nd = Formula{
 	Stencil:    []Point{{Loc: -1, Coeff: 1}, {Loc: 0, Coeff: -2}, {Loc: 1, Coeff: 1}},
 	Derivative: 2,

From 3177c4548585234e429da55b1b92ec1ed8409bfc Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Mon, 28 Mar 2016 09:52:27 +0900
Subject: [PATCH 79/98] fd: change default formula to Forward in Derivative and
 Gradient

---
 fd/diff.go | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 64aadb38..078d8c76 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -38,7 +38,7 @@ func (f Formula) isZero() bool {
 type Settings struct {
 	// Formula is the finite difference formula used
 	// for approximating the derivative.
-	// Zero value will default to the Central formula.
+	// Zero value indicates a default formula.
 	Formula Formula
 	// Step is the distance between points of the stencil.
 	// If equal to 0, formula's default step will be used.
@@ -53,14 +53,14 @@ type Settings struct {
 // Derivative estimates the derivative of the function f at the given location.
 // The finite difference formula, the step size, and other options are
 // specified by settings. If settings is nil, the first derivative will be
-// estimated using the Central formula and a default step size.
+// estimated using the Forward formula and a default step size.
 func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
 	if settings == nil {
 		settings = &Settings{}
 	}
 	formula := settings.Formula
 	if formula.isZero() {
-		formula = Central
+		formula = Forward
 	}
 	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
 		panic("fd: bad formula")
@@ -124,7 +124,7 @@ func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *S
 
 	formula := settings.Formula
 	if formula.isZero() {
-		formula = Central
+		formula = Forward
 	}
 	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
 		panic("fd: bad formula")

From 4a5cb38255a496ba7cceb942f7cc21d6065b6db8 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Mon, 28 Mar 2016 09:53:16 +0900
Subject: [PATCH 80/98] fd: update comment for Gradient

---
 fd/diff.go | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/fd/diff.go b/fd/diff.go
index 078d8c76..3fb45b0a 100644
--- a/fd/diff.go
+++ b/fd/diff.go
@@ -105,12 +105,14 @@ func Derivative(f func(float64) float64, x float64, settings *Settings) float64
 }
 
 // Gradient estimates the gradient of the multivariate function f at the
-// location x. The result is stored in-place into dst if dst is not nil,
-// otherwise a new slice will be allocated and returned. Finite difference
-// kernel and other options are specified by settings. If settings is nil,
-// default settings will be used.
-// Gradient panics if the length of dst and x is not equal, or if the
-// derivative order of the formula is not 1.
+// location x. If dst is not nil, the result will be stored in-place into dst
+// and returned, otherwise a new slice will be allocated first. Finite
+// difference kernel and other options are specified by settings. If settings is
+// nil, the gradient will be estimated using the Forward formula and a default
+// step size.
+//
+// Gradient panics if the length of dst and x is not equal, or if the derivative
+// order of the formula is not 1.
 func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
 	if dst == nil {
 		dst = make([]float64, len(x))

From 47012cd15b339d0b3cdddbbb5929289847b3e1de Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Mon, 28 Mar 2016 09:53:30 +0900
Subject: [PATCH 81/98] fd: update example for Derivative

---
 fd/example_test.go | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/fd/example_test.go b/fd/example_test.go
index f5a01d6a..0b659991 100644
--- a/fd/example_test.go
+++ b/fd/example_test.go
@@ -21,19 +21,18 @@ func ExampleDerivative() {
 	// with a custom step size.
 	df := fd.Derivative(f, 0, &fd.Settings{
 		Formula: fd.Forward,
-		Step:    1e-8,
+		Step:    1e-3,
 	})
 	fmt.Println("f'(0) ≈", df)
 
 	f = func(x float64) float64 {
 		return math.Pow(math.Cos(x), 3)
 	}
-	// Compute the second derivative of f at 0 using the centered
-	// approximation, a custom step size, concurrent evaluation, and a known
-	// function value at x.
+	// Compute the second derivative of f at 0 using
+	// the centered approximation, concurrent evaluation,
+	// and a known function value at x.
 	df = fd.Derivative(f, 0, &fd.Settings{
 		Formula:     fd.Central2nd,
-		Step:        1e-4,
 		Concurrent:  true,
 		OriginKnown: true,
 		OriginValue: f(0),
@@ -41,7 +40,7 @@ func ExampleDerivative() {
 	fmt.Println("f''(0) ≈", df)
 
 	// Output:
-	// f'(0) ≈ 0.999999999994
 	// f'(0) ≈ 1
+	// f'(0) ≈ 0.9999998333333416
 	// f''(0) ≈ -2.999999981767587
 }

From 79be9fe65fc0acd50cc591927cc703564add93aa Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Tue, 29 Mar 2016 16:28:37 +0900
Subject: [PATCH 82/98] fd: make Gradient tests predictable

---
 fd/gradient_test.go | 1 +
 1 file changed, 1 insertion(+)

diff --git a/fd/gradient_test.go b/fd/gradient_test.go
index d065cb3f..09a335e8 100644
--- a/fd/gradient_test.go
+++ b/fd/gradient_test.go
@@ -41,6 +41,7 @@ func (r Rosenbrock) FDf(x []float64, deriv []float64) (sum float64) {
 }
 
 func TestGradient(t *testing.T) {
+	rand.Seed(1)
 	for i, test := range []struct {
 		nDim    int
 		tol     float64

From 3a47b20eaa22da80cd932da539a95407e59e7321 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Mon, 18 Apr 2016 12:52:49 -0400
Subject: [PATCH 83/98] go 1.5.3 -> 1.5.4, 1.6 -> 1.6.1

---
 .travis.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index 137c7f25..bae08859 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -3,8 +3,8 @@ language: go
 # Versions of go that are explicitly supported by gonum.
 go:
  - 1.4.2
- - 1.5.3
- - 1.6
+ - 1.5.4
+ - 1.6.1
 
 # Required for coverage.
 before_install:

From b0943569635c0bba5c38e2d8515146fe902dc3e7 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Sun, 3 Apr 2016 01:53:23 +0900
Subject: [PATCH 84/98] fd: add Jacobian

---
 fd/jacobian.go | 187 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 187 insertions(+)
 create mode 100644 fd/jacobian.go

diff --git a/fd/jacobian.go b/fd/jacobian.go
new file mode 100644
index 00000000..aa2f9d32
--- /dev/null
+++ b/fd/jacobian.go
@@ -0,0 +1,187 @@
+// Copyright ©2016 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 fd
+
+import (
+	"runtime"
+	"sync"
+
+	"github.com/gonum/floats"
+	"github.com/gonum/matrix/mat64"
+)
+
+type JacobianSettings struct {
+	Formula     Formula
+	OriginValue []float64
+	Step        float64
+	Concurrent  bool
+}
+
+// Jacobian approximates the Jacobian matrix of a vector-valued function f at
+// the location x.
+//
+// The Jacobian matrix J is the matrix of all first-order partial derivatives of f.
+// If f maps an n-dimensional vector x to an m-dimensional vector y = f(x), J is
+// an m×n matrix whose elements are given as
+//  J_{i,j} = ∂f_i/∂x_j,
+// or expanded out
+//      [ ∂f_1/∂x_1 ... ∂f_1/∂x_n ]
+//      [     .  .          .     ]
+//  J = [     .      .      .     ]
+//      [     .          .  .     ]
+//      [ ∂f_m/∂x_1 ... ∂f_m/∂x_n ]
+//
+// If dst is not nil, the result will be stored in-place into dst and returned,
+// otherwise a new matrix will be allocated first. Finite difference formula and
+// other options are specified by settings. If settings is nil, the Jacobian
+// will be estimated using the Forward formula and a default step size.
+//
+// Jacobian panics if dst is not nil and its size is not m × len(x), or if the
+// derivative order of the formula is not 1.
+func Jacobian(dst *mat64.Dense, f func(y, x []float64), m int, x []float64, settings *JacobianSettings) *mat64.Dense {
+	n := len(x)
+	if dst == nil {
+		dst = mat64.NewDense(m, n, nil)
+	}
+	r, c := dst.Dims()
+	if r != m || c != n {
+		panic("jacobian: mismatched matrix size")
+	}
+
+	if settings == nil {
+		settings = &JacobianSettings{}
+	}
+	if settings.OriginValue != nil && len(settings.OriginValue) != m {
+		panic("jacobian: mismatched OriginValue slice length")
+	}
+
+	formula := settings.Formula
+	if formula.isZero() {
+		formula = Forward
+	}
+	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
+		panic("jacobian: bad formula")
+	}
+	if formula.Derivative != 1 {
+		panic("jacobian: invalid derivative order")
+	}
+
+	step := settings.Step
+	if step == 0 {
+		step = formula.Step
+	}
+
+	var hasOrigin bool
+	for _, pt := range formula.Stencil {
+		if pt.Loc == 0 {
+			hasOrigin = true
+			break
+		}
+	}
+
+	xcopy := make([]float64, n)
+	origin := settings.OriginValue
+	if hasOrigin && origin == nil {
+		origin = make([]float64, m)
+		copy(xcopy, x)
+		f(origin, xcopy)
+	}
+
+	evals := n * len(formula.Stencil)
+	if hasOrigin {
+		evals -= n
+	}
+	nWorkers := 1
+	if settings.Concurrent {
+		nWorkers = runtime.GOMAXPROCS(0)
+		if nWorkers > evals {
+			nWorkers = evals
+		}
+	}
+
+	if nWorkers == 1 {
+		jacobianSerial(dst, f, x, xcopy, origin, formula, step)
+	} else {
+		jacobianConcurrent(dst, f, x, origin, formula, step, nWorkers)
+	}
+	return dst
+}
+
+func jacobianSerial(dst *mat64.Dense, f func([]float64, []float64), x, xcopy, origin []float64, formula Formula, step float64) {
+	r, c := dst.Dims()
+	y := make([]float64, r)
+	col := make([]float64, r)
+	for j := 0; j < c; j++ {
+		for i := range col {
+			col[i] = 0
+		}
+		for _, pt := range formula.Stencil {
+			if pt.Loc == 0 {
+				floats.AddScaled(col, pt.Coeff, origin)
+			} else {
+				copy(xcopy, x)
+				xcopy[j] += pt.Loc * step
+				f(y, xcopy)
+				floats.AddScaled(col, pt.Coeff, y)
+			}
+		}
+		dst.SetCol(j, col)
+	}
+	dst.Scale(1/step, dst)
+}
+
+func jacobianConcurrent(dst *mat64.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64, nWorkers int) {
+	r, c := dst.Dims()
+	originVec := mat64.NewVector(r, origin)
+	for i := 0; i < r; i++ {
+		for j := 0; j < c; j++ {
+			dst.Set(i, j, 0)
+		}
+	}
+
+	var wg sync.WaitGroup
+	jobs := make(chan jacJob, nWorkers)
+	mus := make([]sync.Mutex, c) // Guard access to individual columns.
+	for i := 0; i < nWorkers; i++ {
+		wg.Add(1)
+		go func() {
+			defer wg.Done()
+			xcopy := make([]float64, len(x))
+			y := make([]float64, r)
+			yVec := mat64.NewVector(r, y)
+			for job := range jobs {
+				if job.pt.Loc != 0 {
+					// Not an origin, we need to evaluate.
+					copy(xcopy, x)
+					xcopy[job.j] += job.pt.Loc * step
+					f(y, xcopy)
+				}
+
+				col := dst.ColView(job.j)
+				mus[job.j].Lock()
+				if job.pt.Loc == 0 {
+					col.AddScaledVec(col, job.pt.Coeff, originVec)
+				} else {
+					col.AddScaledVec(col, job.pt.Coeff, yVec)
+				}
+				mus[job.j].Unlock()
+			}
+		}()
+	}
+	for _, pt := range formula.Stencil {
+		for j := 0; j < c; j++ {
+			jobs <- jacJob{j, pt}
+		}
+	}
+	close(jobs)
+	wg.Wait()
+
+	dst.Scale(1/step, dst)
+}
+
+type jacJob struct {
+	j  int
+	pt Point
+}

From f325327cb37ddf53a8b005ea552c1a449930106d Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Sun, 3 Apr 2016 01:53:45 +0900
Subject: [PATCH 85/98] fd: add tests for Jacobian

---
 fd/jacobian_test.go | 286 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 286 insertions(+)
 create mode 100644 fd/jacobian_test.go

diff --git a/fd/jacobian_test.go b/fd/jacobian_test.go
new file mode 100644
index 00000000..7266f84b
--- /dev/null
+++ b/fd/jacobian_test.go
@@ -0,0 +1,286 @@
+// Copyright ©2016 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 fd
+
+import (
+	"math"
+	"math/rand"
+	"testing"
+
+	"github.com/gonum/floats"
+	"github.com/gonum/matrix/mat64"
+)
+
+func vecFunc13(y, x []float64) {
+	y[0] = 5*x[0] + x[2]*math.Sin(x[1]) + 1
+}
+func vecFunc13Jac(jac *mat64.Dense, x []float64) {
+	jac.Set(0, 0, 5)
+	jac.Set(0, 1, x[2]*math.Cos(x[1]))
+	jac.Set(0, 2, math.Sin(x[1]))
+}
+
+func vecFunc22(y, x []float64) {
+	y[0] = x[0]*x[0]*x[1] + 1
+	y[1] = 5*x[0] + math.Sin(x[1]) + 1
+}
+func vecFunc22Jac(jac *mat64.Dense, x []float64) {
+	jac.Set(0, 0, 2*x[0]*x[1])
+	jac.Set(0, 1, x[0]*x[0])
+	jac.Set(1, 0, 5)
+	jac.Set(1, 1, math.Cos(x[1]))
+}
+
+func vecFunc43(y, x []float64) {
+	y[0] = x[0] + 1
+	y[1] = 5*x[2] + 1
+	y[2] = 4*x[1]*x[1] - 2*x[2] + 1
+	y[3] = x[2]*math.Sin(x[0]) + 1
+}
+func vecFunc43Jac(jac *mat64.Dense, x []float64) {
+	jac.Set(0, 0, 1)
+	jac.Set(0, 1, 0)
+	jac.Set(0, 2, 0)
+	jac.Set(1, 0, 0)
+	jac.Set(1, 1, 0)
+	jac.Set(1, 2, 5)
+	jac.Set(2, 0, 0)
+	jac.Set(2, 1, 8*x[1])
+	jac.Set(2, 2, -2)
+	jac.Set(3, 0, x[2]*math.Cos(x[0]))
+	jac.Set(3, 1, 0)
+	jac.Set(3, 2, math.Sin(x[0]))
+}
+
+func TestJacobian(t *testing.T) {
+	rand.Seed(1)
+
+	// Test with default settings.
+	for tc, test := range []struct {
+		m, n int
+		f    func([]float64, []float64)
+		jac  func(*mat64.Dense, []float64)
+	}{
+		{
+			m:   1,
+			n:   3,
+			f:   vecFunc13,
+			jac: vecFunc13Jac,
+		},
+		{
+			m:   2,
+			n:   2,
+			f:   vecFunc22,
+			jac: vecFunc22Jac,
+		},
+		{
+			m:   4,
+			n:   3,
+			f:   vecFunc43,
+			jac: vecFunc43Jac,
+		},
+	} {
+		const tol = 1e-6
+
+		x := randomSlice(test.n, 10)
+		xcopy := make([]float64, test.n)
+		copy(xcopy, x)
+
+		want := mat64.NewDense(test.m, test.n, nil)
+		test.jac(want, x)
+
+		got := Jacobian(nil, test.f, test.m, x, nil)
+		if !mat64.EqualApprox(want, got, tol) {
+			t.Errorf("Case %d (nil dst, default settings): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (nil dst, default settings): x modified", tc)
+		}
+
+		fillNaNDense(got)
+		Jacobian(got, test.f, test.m, x, nil)
+		if !mat64.EqualApprox(want, got, tol) {
+			t.Errorf("Case %d (default settings): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (default settings): x modified", tc)
+		}
+	}
+
+	// Test with non-default settings.
+	for tc, test := range []struct {
+		m, n    int
+		f       func([]float64, []float64)
+		jac     func(*mat64.Dense, []float64)
+		tol     float64
+		formula Formula
+	}{
+		{
+			m:       1,
+			n:       3,
+			f:       vecFunc13,
+			jac:     vecFunc13Jac,
+			tol:     1e-6,
+			formula: Forward,
+		},
+		{
+			m:       1,
+			n:       3,
+			f:       vecFunc13,
+			jac:     vecFunc13Jac,
+			tol:     1e-6,
+			formula: Backward,
+		},
+		{
+			m:       1,
+			n:       3,
+			f:       vecFunc13,
+			jac:     vecFunc13Jac,
+			tol:     1e-9,
+			formula: Central,
+		},
+		{
+			m:       2,
+			n:       2,
+			f:       vecFunc22,
+			jac:     vecFunc22Jac,
+			tol:     1e-6,
+			formula: Forward,
+		},
+		{
+			m:       2,
+			n:       2,
+			f:       vecFunc22,
+			jac:     vecFunc22Jac,
+			tol:     1e-6,
+			formula: Backward,
+		},
+		{
+			m:       2,
+			n:       2,
+			f:       vecFunc22,
+			jac:     vecFunc22Jac,
+			tol:     1e-9,
+			formula: Central,
+		},
+		{
+			m:       4,
+			n:       3,
+			f:       vecFunc43,
+			jac:     vecFunc43Jac,
+			tol:     1e-6,
+			formula: Forward,
+		},
+		{
+			m:       4,
+			n:       3,
+			f:       vecFunc43,
+			jac:     vecFunc43Jac,
+			tol:     1e-6,
+			formula: Backward,
+		},
+		{
+			m:       4,
+			n:       3,
+			f:       vecFunc43,
+			jac:     vecFunc43Jac,
+			tol:     1e-9,
+			formula: Central,
+		},
+	} {
+		x := randomSlice(test.n, 10)
+		xcopy := make([]float64, test.n)
+		copy(xcopy, x)
+
+		want := mat64.NewDense(test.m, test.n, nil)
+		test.jac(want, x)
+
+		got := Jacobian(nil, test.f, test.m, x, &JacobianSettings{
+			Formula: test.formula,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d (nil dst): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (nil dst): x modified", tc)
+		}
+
+		fillNaNDense(got)
+		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+			Formula: test.formula,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d: unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d: x modified", tc)
+		}
+
+		fillNaNDense(got)
+		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+			Formula:    test.formula,
+			Concurrent: true,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d (concurrent): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (concurrent): x modified", tc)
+		}
+
+		fillNaNDense(got)
+		origin := make([]float64, test.m)
+		test.f(origin, x)
+		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+			Formula:     test.formula,
+			OriginValue: origin,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d (origin): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (origin): x modified", tc)
+		}
+
+		fillNaNDense(got)
+		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+			Formula:     test.formula,
+			OriginValue: origin,
+			Concurrent:  true,
+		})
+		if !mat64.EqualApprox(want, got, test.tol) {
+			t.Errorf("Case %d (concurrent, origin): unexpected Jacobian.\nwant: %v\ngot:  %v",
+				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
+		}
+		if !floats.Equal(x, xcopy) {
+			t.Errorf("Case %d (concurrent, origin): x modified", tc)
+		}
+	}
+}
+
+// randomSlice returns a slice of n elements from the interval [-bound,bound).
+func randomSlice(n int, bound float64) []float64 {
+	x := make([]float64, n)
+	for i := range x {
+		x[i] = 2*bound*rand.Float64() - bound
+	}
+	return x
+}
+
+// fillNaNDense fills the matrix m with NaN values.
+func fillNaNDense(m *mat64.Dense) {
+	r, c := m.Dims()
+	for i := 0; i < r; i++ {
+		for j := 0; j < c; j++ {
+			m.Set(i, j, math.NaN())
+		}
+	}
+}

From 9638d94d1acbb87bb7621c27b12b22e6618116ea Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Sun, 3 Apr 2016 02:08:36 +0900
Subject: [PATCH 86/98] fd: add example for Jacobian

---
 fd/example_test.go | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/fd/example_test.go b/fd/example_test.go
index 0b659991..ea4dbd5a 100644
--- a/fd/example_test.go
+++ b/fd/example_test.go
@@ -9,6 +9,7 @@ import (
 	"math"
 
 	"github.com/gonum/diff/fd"
+	"github.com/gonum/matrix/mat64"
 )
 
 func ExampleDerivative() {
@@ -44,3 +45,24 @@ func ExampleDerivative() {
 	// f'(0) ≈ 0.9999998333333416
 	// f''(0) ≈ -2.999999981767587
 }
+
+func ExampleJacobian() {
+	f := func(dst, x []float64) {
+		dst[0] = x[0] + 1
+		dst[1] = 5 * x[2]
+		dst[2] = 4*x[1]*x[1] - 2*x[2]
+		dst[3] = x[2] * math.Sin(x[0])
+	}
+	x := []float64{1, 2, 3}
+	jac := fd.Jacobian(nil, f, 4, x, &fd.JacobianSettings{
+		Formula:    fd.Central,
+		Concurrent: true,
+	})
+	fmt.Printf("J ≈ %.6v\n", mat64.Formatted(jac, mat64.Prefix("    ")))
+
+	// Output:
+	// J ≈ ⎡       1         0         0⎤
+	//     ⎢       0         0         5⎥
+	//     ⎢       0        16        -2⎥
+	//     ⎣ 1.62091         0  0.841471⎦
+}

From f988c6263a01d7cd643244d44b2e002e5edde7e8 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Wed, 13 Apr 2016 16:50:15 +0900
Subject: [PATCH 87/98] fd: make evaluations in Jacobian completely concurrent

Previously, for simplicity f was evaluated serially even when
JacobianSettings.Concurrent == true.
---
 fd/jacobian.go | 122 ++++++++++++++++++++++++++++---------------------
 1 file changed, 69 insertions(+), 53 deletions(-)

diff --git a/fd/jacobian.go b/fd/jacobian.go
index aa2f9d32..cb39584b 100644
--- a/fd/jacobian.go
+++ b/fd/jacobian.go
@@ -73,26 +73,13 @@ func Jacobian(dst *mat64.Dense, f func(y, x []float64), m int, x []float64, sett
 		step = formula.Step
 	}
 
-	var hasOrigin bool
+	evals := n * len(formula.Stencil)
 	for _, pt := range formula.Stencil {
 		if pt.Loc == 0 {
-			hasOrigin = true
+			evals -= n - 1
 			break
 		}
 	}
-
-	xcopy := make([]float64, n)
-	origin := settings.OriginValue
-	if hasOrigin && origin == nil {
-		origin = make([]float64, m)
-		copy(xcopy, x)
-		f(origin, xcopy)
-	}
-
-	evals := n * len(formula.Stencil)
-	if hasOrigin {
-		evals -= n
-	}
 	nWorkers := 1
 	if settings.Concurrent {
 		nWorkers = runtime.GOMAXPROCS(0)
@@ -100,25 +87,31 @@ func Jacobian(dst *mat64.Dense, f func(y, x []float64), m int, x []float64, sett
 			nWorkers = evals
 		}
 	}
-
 	if nWorkers == 1 {
-		jacobianSerial(dst, f, x, xcopy, origin, formula, step)
+		jacobianSerial(dst, f, x, settings.OriginValue, formula, step)
 	} else {
-		jacobianConcurrent(dst, f, x, origin, formula, step, nWorkers)
+		jacobianConcurrent(dst, f, x, settings.OriginValue, formula, step, nWorkers)
 	}
+
 	return dst
 }
 
-func jacobianSerial(dst *mat64.Dense, f func([]float64, []float64), x, xcopy, origin []float64, formula Formula, step float64) {
-	r, c := dst.Dims()
-	y := make([]float64, r)
-	col := make([]float64, r)
-	for j := 0; j < c; j++ {
+func jacobianSerial(dst *mat64.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64) {
+	m, n := dst.Dims()
+	xcopy := make([]float64, n)
+	y := make([]float64, m)
+	col := make([]float64, m)
+	for j := 0; j < n; j++ {
 		for i := range col {
 			col[i] = 0
 		}
 		for _, pt := range formula.Stencil {
 			if pt.Loc == 0 {
+				if origin == nil {
+					origin = make([]float64, m)
+					copy(xcopy, x)
+					f(origin, xcopy)
+				}
 				floats.AddScaled(col, pt.Coeff, origin)
 			} else {
 				copy(xcopy, x)
@@ -133,51 +126,74 @@ func jacobianSerial(dst *mat64.Dense, f func([]float64, []float64), x, xcopy, or
 }
 
 func jacobianConcurrent(dst *mat64.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64, nWorkers int) {
-	r, c := dst.Dims()
-	originVec := mat64.NewVector(r, origin)
-	for i := 0; i < r; i++ {
-		for j := 0; j < c; j++ {
+	m, n := dst.Dims()
+	for i := 0; i < m; i++ {
+		for j := 0; j < n; j++ {
 			dst.Set(i, j, 0)
 		}
 	}
 
-	var wg sync.WaitGroup
+	var (
+		wg sync.WaitGroup
+		mu = make([]sync.Mutex, n) // Guard access to individual columns.
+	)
+	worker := func(jobs <-chan jacJob) {
+		defer wg.Done()
+		xcopy := make([]float64, n)
+		y := make([]float64, m)
+		yVec := mat64.NewVector(m, y)
+		for job := range jobs {
+			copy(xcopy, x)
+			xcopy[job.j] += job.pt.Loc * step
+			f(y, xcopy)
+			col := dst.ColView(job.j)
+			mu[job.j].Lock()
+			col.AddScaledVec(col, job.pt.Coeff, yVec)
+			mu[job.j].Unlock()
+		}
+	}
 	jobs := make(chan jacJob, nWorkers)
-	mus := make([]sync.Mutex, c) // Guard access to individual columns.
 	for i := 0; i < nWorkers; i++ {
 		wg.Add(1)
-		go func() {
-			defer wg.Done()
-			xcopy := make([]float64, len(x))
-			y := make([]float64, r)
-			yVec := mat64.NewVector(r, y)
-			for job := range jobs {
-				if job.pt.Loc != 0 {
-					// Not an origin, we need to evaluate.
-					copy(xcopy, x)
-					xcopy[job.j] += job.pt.Loc * step
-					f(y, xcopy)
-				}
-
-				col := dst.ColView(job.j)
-				mus[job.j].Lock()
-				if job.pt.Loc == 0 {
-					col.AddScaledVec(col, job.pt.Coeff, originVec)
-				} else {
-					col.AddScaledVec(col, job.pt.Coeff, yVec)
-				}
-				mus[job.j].Unlock()
-			}
-		}()
+		go worker(jobs)
 	}
+	var (
+		hasOrigin   bool
+		originCoeff float64
+	)
 	for _, pt := range formula.Stencil {
-		for j := 0; j < c; j++ {
+		if pt.Loc == 0 {
+			hasOrigin = true
+			originCoeff = pt.Coeff
+			continue
+		}
+		for j := 0; j < n; j++ {
 			jobs <- jacJob{j, pt}
 		}
 	}
+	if hasOrigin && origin == nil {
+		wg.Add(1)
+		go func() {
+			defer wg.Done()
+			origin = make([]float64, m)
+			xcopy := make([]float64, n)
+			copy(xcopy, x)
+			f(origin, xcopy)
+		}()
+	}
 	close(jobs)
 	wg.Wait()
 
+	if hasOrigin {
+		// The formula evaluated at x, we need to add scaled origin to
+		// all columns of dst.
+		originVec := mat64.NewVector(m, origin)
+		for j := 0; j < n; j++ {
+			col := dst.ColView(j)
+			col.AddScaledVec(col, originCoeff, originVec)
+		}
+	}
+
 	dst.Scale(1/step, dst)
 }
 

From 43de7bc3d5acaad7b17ee58dc72d0b61f6300c1e Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Tue, 19 Apr 2016 09:49:22 +0900
Subject: [PATCH 88/98] fd: change signature of Jacobian, dst must be non-nil

---
 fd/example_test.go  |  4 ++--
 fd/jacobian.go      | 27 ++++++++++++---------------
 fd/jacobian_test.go | 32 +++++++-------------------------
 3 files changed, 21 insertions(+), 42 deletions(-)

diff --git a/fd/example_test.go b/fd/example_test.go
index ea4dbd5a..1856ffa6 100644
--- a/fd/example_test.go
+++ b/fd/example_test.go
@@ -53,8 +53,8 @@ func ExampleJacobian() {
 		dst[2] = 4*x[1]*x[1] - 2*x[2]
 		dst[3] = x[2] * math.Sin(x[0])
 	}
-	x := []float64{1, 2, 3}
-	jac := fd.Jacobian(nil, f, 4, x, &fd.JacobianSettings{
+	jac := mat64.NewDense(4, 3, nil)
+	fd.Jacobian(jac, f, []float64{1, 2, 3}, &fd.JacobianSettings{
 		Formula:    fd.Central,
 		Concurrent: true,
 	})
diff --git a/fd/jacobian.go b/fd/jacobian.go
index cb39584b..89d52e15 100644
--- a/fd/jacobian.go
+++ b/fd/jacobian.go
@@ -20,7 +20,11 @@ type JacobianSettings struct {
 }
 
 // Jacobian approximates the Jacobian matrix of a vector-valued function f at
-// the location x.
+// the location x and stores the result in-place into dst.
+//
+// Finite difference formula and other options are specified by settings. If
+// settings is nil, the Jacobian will be estimated using the Forward formula and
+// a default step size.
 //
 // The Jacobian matrix J is the matrix of all first-order partial derivatives of f.
 // If f maps an n-dimensional vector x to an m-dimensional vector y = f(x), J is
@@ -33,20 +37,15 @@ type JacobianSettings struct {
 //      [     .          .  .     ]
 //      [ ∂f_m/∂x_1 ... ∂f_m/∂x_n ]
 //
-// If dst is not nil, the result will be stored in-place into dst and returned,
-// otherwise a new matrix will be allocated first. Finite difference formula and
-// other options are specified by settings. If settings is nil, the Jacobian
-// will be estimated using the Forward formula and a default step size.
-//
-// Jacobian panics if dst is not nil and its size is not m × len(x), or if the
-// derivative order of the formula is not 1.
-func Jacobian(dst *mat64.Dense, f func(y, x []float64), m int, x []float64, settings *JacobianSettings) *mat64.Dense {
+// dst must be non-nil, the number of its columns must equal the length of x, and
+// the derivative order of the formula must be 1, otherwise Jacobian will panic.
+func Jacobian(dst *mat64.Dense, f func(y, x []float64), x []float64, settings *JacobianSettings) {
 	n := len(x)
-	if dst == nil {
-		dst = mat64.NewDense(m, n, nil)
+	if n == 0 {
+		panic("jacobian: x has zero length")
 	}
-	r, c := dst.Dims()
-	if r != m || c != n {
+	m, c := dst.Dims()
+	if c != n {
 		panic("jacobian: mismatched matrix size")
 	}
 
@@ -92,8 +91,6 @@ func Jacobian(dst *mat64.Dense, f func(y, x []float64), m int, x []float64, sett
 	} else {
 		jacobianConcurrent(dst, f, x, settings.OriginValue, formula, step, nWorkers)
 	}
-
-	return dst
 }
 
 func jacobianSerial(dst *mat64.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64) {
diff --git a/fd/jacobian_test.go b/fd/jacobian_test.go
index 7266f84b..89f1af11 100644
--- a/fd/jacobian_test.go
+++ b/fd/jacobian_test.go
@@ -91,17 +91,9 @@ func TestJacobian(t *testing.T) {
 		want := mat64.NewDense(test.m, test.n, nil)
 		test.jac(want, x)
 
-		got := Jacobian(nil, test.f, test.m, x, nil)
-		if !mat64.EqualApprox(want, got, tol) {
-			t.Errorf("Case %d (nil dst, default settings): unexpected Jacobian.\nwant: %v\ngot:  %v",
-				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
-		}
-		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %d (nil dst, default settings): x modified", tc)
-		}
-
+		got := mat64.NewDense(test.m, test.n, nil)
 		fillNaNDense(got)
-		Jacobian(got, test.f, test.m, x, nil)
+		Jacobian(got, test.f, x, nil)
 		if !mat64.EqualApprox(want, got, tol) {
 			t.Errorf("Case %d (default settings): unexpected Jacobian.\nwant: %v\ngot:  %v",
 				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
@@ -199,19 +191,9 @@ func TestJacobian(t *testing.T) {
 		want := mat64.NewDense(test.m, test.n, nil)
 		test.jac(want, x)
 
-		got := Jacobian(nil, test.f, test.m, x, &JacobianSettings{
-			Formula: test.formula,
-		})
-		if !mat64.EqualApprox(want, got, test.tol) {
-			t.Errorf("Case %d (nil dst): unexpected Jacobian.\nwant: %v\ngot:  %v",
-				tc, mat64.Formatted(want, mat64.Prefix("      ")), mat64.Formatted(got, mat64.Prefix("      ")))
-		}
-		if !floats.Equal(x, xcopy) {
-			t.Errorf("Case %d (nil dst): x modified", tc)
-		}
-
+		got := mat64.NewDense(test.m, test.n, nil)
 		fillNaNDense(got)
-		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+		Jacobian(got, test.f, x, &JacobianSettings{
 			Formula: test.formula,
 		})
 		if !mat64.EqualApprox(want, got, test.tol) {
@@ -223,7 +205,7 @@ func TestJacobian(t *testing.T) {
 		}
 
 		fillNaNDense(got)
-		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+		Jacobian(got, test.f, x, &JacobianSettings{
 			Formula:    test.formula,
 			Concurrent: true,
 		})
@@ -238,7 +220,7 @@ func TestJacobian(t *testing.T) {
 		fillNaNDense(got)
 		origin := make([]float64, test.m)
 		test.f(origin, x)
-		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+		Jacobian(got, test.f, x, &JacobianSettings{
 			Formula:     test.formula,
 			OriginValue: origin,
 		})
@@ -251,7 +233,7 @@ func TestJacobian(t *testing.T) {
 		}
 
 		fillNaNDense(got)
-		Jacobian(got, test.f, test.m, x, &JacobianSettings{
+		Jacobian(got, test.f, x, &JacobianSettings{
 			Formula:     test.formula,
 			OriginValue: origin,
 			Concurrent:  true,

From 0efc47f30990bc2948809406e9d7a5d58282d068 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir.chalupecky@gmail.com>
Date: Tue, 19 Apr 2016 09:57:04 +0900
Subject: [PATCH 89/98] fd: handle correctly repeated origins in Jacobian

---
 fd/jacobian.go | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/fd/jacobian.go b/fd/jacobian.go
index 89d52e15..1aab264c 100644
--- a/fd/jacobian.go
+++ b/fd/jacobian.go
@@ -154,20 +154,17 @@ func jacobianConcurrent(dst *mat64.Dense, f func([]float64, []float64), x, origi
 		wg.Add(1)
 		go worker(jobs)
 	}
-	var (
-		hasOrigin   bool
-		originCoeff float64
-	)
+	var hasOrigin bool
 	for _, pt := range formula.Stencil {
 		if pt.Loc == 0 {
 			hasOrigin = true
-			originCoeff = pt.Coeff
 			continue
 		}
 		for j := 0; j < n; j++ {
 			jobs <- jacJob{j, pt}
 		}
 	}
+	close(jobs)
 	if hasOrigin && origin == nil {
 		wg.Add(1)
 		go func() {
@@ -178,16 +175,22 @@ func jacobianConcurrent(dst *mat64.Dense, f func([]float64, []float64), x, origi
 			f(origin, xcopy)
 		}()
 	}
-	close(jobs)
 	wg.Wait()
 
 	if hasOrigin {
 		// The formula evaluated at x, we need to add scaled origin to
-		// all columns of dst.
+		// all columns of dst. Iterate again over all Formula points
+		// because we don't forbid repeated locations.
+
 		originVec := mat64.NewVector(m, origin)
-		for j := 0; j < n; j++ {
-			col := dst.ColView(j)
-			col.AddScaledVec(col, originCoeff, originVec)
+		for _, pt := range formula.Stencil {
+			if pt.Loc != 0 {
+				continue
+			}
+			for j := 0; j < n; j++ {
+				col := dst.ColView(j)
+				col.AddScaledVec(col, pt.Coeff, originVec)
+			}
 		}
 	}
 

From 88667fadb3e02688fcd98d5a281e5040ea20c3b1 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 21 Apr 2016 10:11:12 -0400
Subject: [PATCH 90/98] travis: go 1.6.1 -> 1.6.2

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index bae08859..769ef94a 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 go:
  - 1.4.2
  - 1.5.4
- - 1.6.1
+ - 1.6.2
 
 # Required for coverage.
 before_install:

From ab821a9b706ca871177efe652dae36e49646218f Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Thu, 2 Jun 2016 10:04:59 -0400
Subject: [PATCH 91/98] travis: go1.7beta1

Add go1.7beta1 to build matrix, PTAL.
---
 .travis.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.travis.yml b/.travis.yml
index 769ef94a..3c849815 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,6 +5,7 @@ go:
  - 1.4.2
  - 1.5.4
  - 1.6.2
+ - 1.7beta1
 
 # Required for coverage.
 before_install:

From 40819595ff6164231c6ed349ffc6d277fdf26f15 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 17 Jun 2016 14:44:00 -0400
Subject: [PATCH 92/98] go 1.7beta1 -> go 1.7beta2

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 3c849815..aa349807 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,7 @@ go:
  - 1.4.2
  - 1.5.4
  - 1.6.2
- - 1.7beta1
+ - 1.7beta2
 
 # Required for coverage.
 before_install:

From 5ca3bd81dea1fb567f18a48f533b55c8482651e5 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Tue, 19 Jul 2016 16:37:03 -0400
Subject: [PATCH 93/98] go 1.7beta2 -> go 1.7rc2, go 1.6.2 -> go 1.6.3

---
 .travis.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index aa349807..d773cdd9 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,8 +4,8 @@ language: go
 go:
  - 1.4.2
  - 1.5.4
- - 1.6.2
- - 1.7beta2
+ - 1.6.3
+ - 1.7rc2
 
 # Required for coverage.
 before_install:

From 8ac35ff5ad00b0ecf6183dafc57dcf36c8de2f13 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Fri, 22 Jul 2016 09:21:22 -0400
Subject: [PATCH 94/98] go 1.7rc2 -> go 1.7rc3

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index d773cdd9..d41d68aa 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,7 +5,7 @@ go:
  - 1.4.2
  - 1.5.4
  - 1.6.3
- - 1.7rc2
+ - 1.7rc3
 
 # Required for coverage.
 before_install:

From 2f501cd00c9d73638d7fd0fdd925eb9e60fab85c Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Tue, 16 Aug 2016 08:44:50 -0400
Subject: [PATCH 95/98] go 1.7rc3 -> go 1.7, drop go 1.4.2

---
 .travis.yml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index d41d68aa..c5d2fa7d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,10 +2,9 @@ language: go
 
 # Versions of go that are explicitly supported by gonum.
 go:
- - 1.4.2
  - 1.5.4
  - 1.6.3
- - 1.7rc3
+ - 1.7
 
 # Required for coverage.
 before_install:

From 44c60c262ed183adbb3a0b9ebfa6bda466113110 Mon Sep 17 00:00:00 2001
From: Jonathan J Lawlor <jonathan.lawlor@gmail.com>
Date: Tue, 25 Oct 2016 08:45:43 -0400
Subject: [PATCH 96/98] go 1.7 -> go 1.7.3

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index c5d2fa7d..869de376 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,7 @@ language: go
 go:
  - 1.5.4
  - 1.6.3
- - 1.7
+ - 1.7.3
 
 # Required for coverage.
 before_install:

From 069b4cb67446c7e03c842bbffef2b9024460c40c Mon Sep 17 00:00:00 2001
From: Evert Lammerts <evert@lucipher.net>
Date: Thu, 17 Nov 2016 21:36:48 +0100
Subject: [PATCH 97/98] travis and coverage script changes

---
 .travis.yml              |  2 +-
 .travis/test-coverage.sh | 68 ++++++++++++++++++----------------------
 2 files changed, 31 insertions(+), 39 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index 869de376..621c6e19 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -17,5 +17,5 @@ script:
  - go get -d -t -v ./...
  - go build -v ./...
  - go test -v ./...
- - diff <(gofmt -d .) <("")
+ - test -z "`gofmt -d .`"
  - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash ./.travis/test-coverage.sh; fi
diff --git a/.travis/test-coverage.sh b/.travis/test-coverage.sh
index 854f8c3a..c143abfc 100755
--- a/.travis/test-coverage.sh
+++ b/.travis/test-coverage.sh
@@ -1,43 +1,35 @@
 #!/bin/bash
 
-# based on http://stackoverflow.com/questions/21126011/is-it-possible-to-post-coverage-for-multiple-packages-to-coveralls
-# with script found at https://github.com/gopns/gopns/blob/master/test-coverage.sh
+PROFILE_OUT=`pwd`/profile.out
+ACC_OUT=`pwd`/acc.out
 
-echo "mode: set" > acc.out
-returnval=`go test -v -coverprofile=profile.out`
-echo ${returnval}
-if [[ ${returnval} != *FAIL* ]]
-then
-	if [ -f profile.out ]
-	then
-		cat profile.out | grep -v "mode: set" >> acc.out
-	fi
-else
-	exit 1
-fi
+testCover() {
+	# set the return value to 0 (succesful)
+	retval=0
+	# get the directory to check from the parameter. Default to '.'
+	d=${1:-.}
+	# skip if there are no Go files here
+	ls $d/*.go &> /dev/null || return $retval
+	# switch to the directory to check
+	pushd $d > /dev/null
+	# create the coverage profile
+	coverageresult=`go test -v -coverprofile=$PROFILE_OUT`
+	# output the result so we can check the shell output
+	echo ${coverageresult}
+	# append the results to acc.out if coverage didn't fail, else set the retval to 1 (failed)
+	( [[ ${coverageresult} == *FAIL* ]] && retval=1 ) || ( [ -f $PROFILE_OUT ] && grep -v "mode: set" $PROFILE_OUT >> $ACC_OUT )
+	# return to our working dir
+	popd > /dev/null
+	# return our return value
+	return $retval
+}
 
-for Dir in $(find ./* -maxdepth 10 -type d );
-do
-	if ls $Dir/*.go &> /dev/null;
-	then
-		echo $Dir
-		returnval=`go test -v -coverprofile=profile.out $Dir`
-		echo ${returnval}
-		if [[ ${returnval} != *FAIL* ]]
-		then
-    		if [ -f profile.out ]
-    		then
-        		cat profile.out | grep -v "mode: set" >> acc.out
-    		fi
-    	else
-    		exit 1
-    	fi
-    fi
-done
-if [ -n "$COVERALLS_TOKEN" ]
-then
-	$HOME/gopath/bin/goveralls -coverprofile=acc.out -service=travis-ci -repotoken $COVERALLS_TOKEN
-fi
+# Init acc.out
+echo "mode: set" > $ACC_OUT
+
+# Run test coverage on all directories containing go files
+find . -maxdepth 10 -type d | while read d; do testCover $d || exit; done
+
+# Upload the coverage profile to coveralls.io
+[ -n "$COVERALLS_TOKEN" ] && goveralls -coverprofile=$ACC_OUT -service=travis-ci -repotoken $COVERALLS_TOKEN
 
-rm -rf ./profile.out
-rm -rf ./acc.out

From 2d52176cf20386e1c628cb9843cdb29e8b678747 Mon Sep 17 00:00:00 2001
From: Evert Lammerts <evert@lucipher.net>
Date: Fri, 18 Nov 2016 21:28:26 +0100
Subject: [PATCH 98/98] removed backticks and replaced pwd

---
 .travis.yml              | 2 +-
 .travis/test-coverage.sh | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/.travis.yml b/.travis.yml
index 621c6e19..e0b47f93 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -17,5 +17,5 @@ script:
  - go get -d -t -v ./...
  - go build -v ./...
  - go test -v ./...
- - test -z "`gofmt -d .`"
+ - test -z "$(gofmt -d .)"
  - if [[ $TRAVIS_SECURE_ENV_VARS = "true" ]]; then bash ./.travis/test-coverage.sh; fi
diff --git a/.travis/test-coverage.sh b/.travis/test-coverage.sh
index c143abfc..7df8aa6a 100755
--- a/.travis/test-coverage.sh
+++ b/.travis/test-coverage.sh
@@ -1,7 +1,7 @@
 #!/bin/bash
 
-PROFILE_OUT=`pwd`/profile.out
-ACC_OUT=`pwd`/acc.out
+PROFILE_OUT=$PWD/profile.out
+ACC_OUT=$PWD/acc.out
 
 testCover() {
 	# set the return value to 0 (succesful)