Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
203 lines (186 sloc) 4.77 KB
// 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 (
"sync"
"gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/mat"
)
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 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
// 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 ]
//
// 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 *mat.Dense, f func(y, x []float64), x []float64, settings *JacobianSettings) {
n := len(x)
if n == 0 {
panic("jacobian: x has zero length")
}
m, c := dst.Dims()
if c != n {
panic("jacobian: mismatched matrix size")
}
// Default settings.
formula := Forward
step := formula.Step
var originValue []float64
var concurrent bool
// Use user settings if provided.
if settings != nil {
if !settings.Formula.isZero() {
formula = settings.Formula
step = formula.Step
checkFormula(formula)
if formula.Derivative != 1 {
panic(badDerivOrder)
}
}
if settings.Step != 0 {
step = settings.Step
}
originValue = settings.OriginValue
if originValue != nil && len(originValue) != m {
panic("jacobian: mismatched OriginValue slice length")
}
concurrent = settings.Concurrent
}
evals := n * len(formula.Stencil)
for _, pt := range formula.Stencil {
if pt.Loc == 0 {
evals -= n - 1
break
}
}
nWorkers := computeWorkers(concurrent, evals)
if nWorkers == 1 {
jacobianSerial(dst, f, x, originValue, formula, step)
return
}
jacobianConcurrent(dst, f, x, originValue, formula, step, nWorkers)
}
func jacobianSerial(dst *mat.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)
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 *mat.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64, nWorkers int) {
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
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 := mat.NewVecDense(m, y)
var col mat.VecDense
for job := range jobs {
copy(xcopy, x)
xcopy[job.j] += job.pt.Loc * step
f(y, xcopy)
col.ColViewOf(dst, job.j)
mu[job.j].Lock()
col.AddScaledVec(&col, job.pt.Coeff, yVec)
mu[job.j].Unlock()
}
}
jobs := make(chan jacJob, nWorkers)
for i := 0; i < nWorkers; i++ {
wg.Add(1)
go worker(jobs)
}
var hasOrigin bool
for _, pt := range formula.Stencil {
if pt.Loc == 0 {
hasOrigin = true
continue
}
for j := 0; j < n; j++ {
jobs <- jacJob{j, pt}
}
}
close(jobs)
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)
}()
}
wg.Wait()
if hasOrigin {
// The formula evaluated at x, we need to add scaled origin to
// all columns of dst. Iterate again over all Formula points
// because we don't forbid repeated locations.
originVec := mat.NewVecDense(m, origin)
for _, pt := range formula.Stencil {
if pt.Loc != 0 {
continue
}
var col mat.VecDense
for j := 0; j < n; j++ {
col.ColViewOf(dst, j)
col.AddScaledVec(&col, pt.Coeff, originVec)
}
}
}
dst.Scale(1/step, dst)
}
type jacJob struct {
j int
pt Point
}
You can’t perform that action at this time.