Add Laplacian and CrossLaplacian difference functions (#154)

* Add Laplacian and CrossLaplacian difference functions

* use usesOrigin

diff/fd/crosslaplacian.go

@@ -0,0 +1,186 @@
+// Copyright ©2017 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"
+	"sync"
+)
+
+// CrossLaplacian computes a Laplacian-like quantity for a function of two vectors
+// at the locations x and y.
+// It computes
+//  ∇_y · ∇_x f(x,y) = \sum_i ∂^2 f(x,y)/∂x_i ∂y_i
+// The two input vector lengths must be the same.
+//
+// Finite difference formula and other options are specified by settings. If
+// settings is nil, CrossLaplacian will be estimated using the Forward formula and
+// a default step size.
+//
+// CrossLaplacian panics if the two input vectors are not the same length, or if
+// the derivative order of the formula is not 1.
+func CrossLaplacian(f func(x, y []float64) float64, x, y []float64, settings *Settings) float64 {
+	n := len(x)
+	if n == 0 {
+		panic("crosslaplacian: x has zero length")
+	}
+	if len(x) != len(y) {
+		panic("crosslaplacian: input vector length mismatch")
+	}
+
+	// Default settings.
+	formula := Forward
+	step := math.Sqrt(formula.Step) // Use the sqrt because taking derivatives of derivatives.
+	var originValue float64
+	var originKnown, concurrent bool
+
+	// Use user settings if provided.
+	if settings != nil {
+		if !settings.Formula.isZero() {
+			formula = settings.Formula
+			step = math.Sqrt(formula.Step)
+			checkFormula(formula)
+			if formula.Derivative != 1 {
+				panic(badDerivOrder)
+			}
+		}
+		if settings.Step != 0 {
+			if settings.Step < 0 {
+				panic(negativeStep)
+			}
+			step = settings.Step
+		}
+		originKnown = settings.OriginKnown
+		originValue = settings.OriginValue
+		concurrent = settings.Concurrent
+	}
+
+	evals := n * len(formula.Stencil) * len(formula.Stencil)
+	if usesOrigin(formula.Stencil) {
+		evals -= n
+	}
+
+	nWorkers := computeWorkers(concurrent, evals)
+	if nWorkers == 1 {
+		return crossLaplacianSerial(f, x, y, formula.Stencil, step, originKnown, originValue)
+	}
+	return crossLaplacianConcurrent(nWorkers, evals, f, x, y, formula.Stencil, step, originKnown, originValue)
+}
+
+func crossLaplacianSerial(f func(x, y []float64) float64, x, y []float64, stencil []Point, step float64, originKnown bool, originValue float64) float64 {
+	n := len(x)
+	xCopy := make([]float64, len(x))
+	yCopy := make([]float64, len(y))
+	fo := func() float64 {
+		// Copy x and y in case they are modified during the call.
+		copy(xCopy, x)
+		copy(yCopy, y)
+		return f(x, y)
+	}
+	origin := getOrigin(originKnown, originValue, fo, stencil)
+
+	is2 := 1 / (step * step)
+	var laplacian float64
+	for i := 0; i < n; i++ {
+		for _, pty := range stencil {
+			for _, ptx := range stencil {
+				var v float64
+				if ptx.Loc == 0 && pty.Loc == 0 {
+					v = origin
+				} else {
+					// Copying the data anew has two benefits. First, it
+					// avoids floating point issues where adding and then
+					// subtracting the step don't return to the exact same
+					// location. Secondly, it protects against the function
+					// modifying the input data.
+					copy(yCopy, y)
+					copy(xCopy, x)
+					yCopy[i] += pty.Loc * step
+					xCopy[i] += ptx.Loc * step
+					v = f(xCopy, yCopy)
+				}
+				laplacian += v * ptx.Coeff * pty.Coeff * is2
+			}
+		}
+	}
+	return laplacian
+}
+
+func crossLaplacianConcurrent(nWorkers, evals int, f func(x, y []float64) float64, x, y []float64, stencil []Point, step float64, originKnown bool, originValue float64) float64 {
+	n := len(x)
+	type run struct {
+		i          int
+		xIdx, yIdx int
+		result     float64
+	}
+
+	send := make(chan run, evals)
+	ans := make(chan run, evals)
+
+	var originWG sync.WaitGroup
+	hasOrigin := usesOrigin(stencil)
+	if hasOrigin {
+		originWG.Add(1)
+		// Launch worker to compute the origin.
+		go func() {
+			defer originWG.Done()
+			xCopy := make([]float64, len(x))
+			yCopy := make([]float64, len(y))
+			copy(xCopy, x)
+			copy(yCopy, y)
+			originValue = f(xCopy, yCopy)
+		}()
+	}
+
+	var workerWG sync.WaitGroup
+	// Launch workers.
+	for i := 0; i < nWorkers; i++ {
+		workerWG.Add(1)
+		go func(send <-chan run, ans chan<- run) {
+			defer workerWG.Done()
+			xCopy := make([]float64, len(x))
+			yCopy := make([]float64, len(y))
+			for r := range send {
+				if stencil[r.xIdx].Loc == 0 && stencil[r.yIdx].Loc == 0 {
+					originWG.Wait()
+					r.result = originValue
+				} else {
+					// See crossLaplacianSerial for comment on the copy.
+					copy(xCopy, x)
+					copy(yCopy, y)
+					xCopy[r.i] += stencil[r.xIdx].Loc * step
+					yCopy[r.i] += stencil[r.yIdx].Loc * step
+					r.result = f(xCopy, yCopy)
+				}
+				ans <- r
+			}
+		}(send, ans)
+	}
+
+	// Launch the distributor, which sends all of runs.
+	go func(send chan<- run) {
+		for i := 0; i < n; i++ {
+			for xIdx := range stencil {
+				for yIdx := range stencil {
+					send <- run{
+						i: i, xIdx: xIdx, yIdx: yIdx,
+					}
+				}
+			}
+		}
+		close(send)
+		// Wait for all the workers to quit, then close the ans channel.
+		workerWG.Wait()
+		close(ans)
+	}(send)
+
+	// Read in the results.
+	is2 := 1 / (step * step)
+	var laplacian float64
+	for r := range ans {
+		laplacian += r.result * stencil[r.xIdx].Coeff * stencil[r.yIdx].Coeff * is2
+	}
+	return laplacian
+}