Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add Lapliacian and CrossLaplacian difference functions
- Loading branch information
Showing
7 changed files
with
564 additions
and
55 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,190 @@ | ||
// 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, LaplacianTwo 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("laplaciantwo: x has zero length") | ||
} | ||
if len(x) != len(y) { | ||
panic("laplaciantwo: 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) | ||
for _, pt := range formula.Stencil { | ||
if pt.Loc == 0 { | ||
evals -= n | ||
} | ||
} | ||
|
||
nWorkers := computeWorkers(concurrent, evals) | ||
if nWorkers == 1 { | ||
return crossLaplacianSerial(f, x, y, formula.Stencil, step, originKnown, originValue) | ||
} else { | ||
return crossLaplacianConcurrent(nWorkers, evals, f, x, y, formula.Stencil, step, originKnown, originValue) | ||
} | ||
panic("fix") | ||
} | ||
|
||
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 | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,111 @@ | ||
// 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 ( | ||
"testing" | ||
|
||
"gonum.org/v1/gonum/floats" | ||
"gonum.org/v1/gonum/mat" | ||
) | ||
|
||
type CrossLaplacianTester interface { | ||
Func(x, y []float64) float64 | ||
CrossLaplacian(x, y []float64) float64 | ||
} | ||
|
||
type WrapperCL struct { | ||
Tester HessianTester | ||
} | ||
|
||
func (WrapperCL) constructZ(x, y []float64) []float64 { | ||
z := make([]float64, len(x)+len(y)) | ||
copy(z, x) | ||
copy(z[len(x):], y) | ||
return z | ||
} | ||
|
||
func (w WrapperCL) Func(x, y []float64) float64 { | ||
z := w.constructZ(x, y) | ||
return w.Tester.Func(z) | ||
} | ||
|
||
func (w WrapperCL) CrossLaplacian(x, y []float64) float64 { | ||
z := w.constructZ(x, y) | ||
hess := mat.NewSymDense(len(z), nil) | ||
w.Tester.Hess(hess, z) | ||
// The LaplacianTwo is the trace of the off-diagonal block of the Hessian. | ||
var l float64 | ||
for i := 0; i < len(x); i++ { | ||
l += hess.At(i, i+len(x)) | ||
} | ||
return l | ||
} | ||
|
||
func TestCrossLaplacian(t *testing.T) { | ||
for cas, test := range []struct { | ||
l CrossLaplacianTester | ||
x, y []float64 | ||
settings *Settings | ||
tol float64 | ||
}{ | ||
{ | ||
l: WrapperCL{Watson{}}, | ||
x: []float64{0.2, 0.3}, | ||
y: []float64{0.1, 0.4}, | ||
tol: 1e-3, | ||
}, | ||
{ | ||
l: WrapperCL{Watson{}}, | ||
x: []float64{2, 3, 1}, | ||
y: []float64{1, 4, 1}, | ||
tol: 1e-3, | ||
}, | ||
{ | ||
l: WrapperCL{ConstFunc(6)}, | ||
x: []float64{2, -3, 1}, | ||
y: []float64{1, 4, -5}, | ||
tol: 1e-6, | ||
}, | ||
{ | ||
l: WrapperCL{LinearFunc{w: []float64{10, 6, -1, 5}, c: 5}}, | ||
x: []float64{3, 1}, | ||
y: []float64{8, 6}, | ||
tol: 1e-6, | ||
}, | ||
{ | ||
l: WrapperCL{QuadFunc{ | ||
a: mat.NewSymDense(4, []float64{ | ||
10, 2, 1, 9, | ||
2, 5, -3, 4, | ||
1, -3, 6, 2, | ||
9, 4, 2, -14, | ||
}), | ||
b: mat.NewVector(4, []float64{3, -2, -1, 4}), | ||
c: 5, | ||
}}, | ||
x: []float64{-1.6, -3}, | ||
y: []float64{1.8, 3.4}, | ||
tol: 1e-6, | ||
}, | ||
} { | ||
got := CrossLaplacian(test.l.Func, test.x, test.y, test.settings) | ||
want := test.l.CrossLaplacian(test.x, test.y) | ||
if !floats.EqualWithinAbsOrRel(got, want, test.tol, test.tol) { | ||
t.Errorf("Cas %d: LaplacianTwo mismatch serial. got %v, want %v", cas, got, want) | ||
} | ||
|
||
// Test that concurrency works. | ||
settings := test.settings | ||
if settings == nil { | ||
settings = &Settings{} | ||
} | ||
settings.Concurrent = true | ||
got2 := CrossLaplacian(test.l.Func, test.x, test.y, settings) | ||
if !floats.EqualWithinAbsOrRel(got, got2, 1e-6, 1e-6) { | ||
t.Errorf("Cas %d: Laplacian mismatch. got %v, want %v", cas, got2, got) | ||
} | ||
} | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.