diff/fd/crosslaplacian.go

// 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
}