WIP: Rework and Simplify Internals (#109)

* Simplify the internals

* Fix test and remove redundant print statement

* Add extrapolate_fdm

* Add Richardson

* Add method of extrapolate_fdm with given h

* Exploit symmetry of FDM in extrapolate

* Update README and add section about multivariate derivatives

* Update Manifest.toml for docs

* Fix docs build

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Fix style of short-form functions

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Update src/methods.jl

* Update src/methods.jl

* Fix test for Base.show

* Improve typing

* Update src/methods.jl

Co-authored-by: Lyndon White <[email protected]>

* Fix typo in comment

* Apply Lyndon's style suggestion

* Inline simple methods

* Add types for power and breaktol

* Increase version number and add Richardson to compat

* Let grid functions always return Vectors

* Restrict docstring of add_tiny

Co-authored-by: Lyndon White <[email protected]>

Project.toml

@@ -1,16 +1,18 @@
 name = "FiniteDifferences"
 uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
-version = "0.10.9"
+version = "0.11.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Richardson = "708f8203-808e-40c0-ba2d-98a6953ed40d"
 
 [compat]
 ChainRulesCore = "0.9"
 julia = "1"
+Richardson = "1.2"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

README.md

@@ -24,38 +24,42 @@ Right now here are the differences:
  - FiniteDiff.jl is carefully optimized to minimize allocations
  - FiniteDiff.jl supports coloring vectors for efficient calculation of sparse Jacobians
 
-## Examples
+
+#### FDM.jl
+This package was formerly called FDM.jl. We recommend users of FDM.jl [update to FiniteDifferences.jl](https://github.com/JuliaDiff/FiniteDifferences.jl/issues/37).
+
+
+## Example: Scalar Derivatives
 
 Compute the first derivative of `sin` with a 5th order central method:
 
 ```julia
 julia> central_fdm(5, 1)(sin, 1) - cos(1)
--1.247890679678676e-13
+-2.4313884239290928e-14
 ```
+
 Compute the second derivative of `sin` with a 5th order central method:
 
 ```julia
 julia> central_fdm(5, 2)(sin, 1) + sin(1)
-9.747314066999024e-12
+-4.367495254342657e-11
 ```
 
-Construct a FiniteDifferences on a custom grid:
+The functions `forward_fdm` and `backward_fdm` can be used to construct
+forward differences and backward differences respectively.
+
+Alternatively, you can construct a finite difference method on a custom grid:
 
 ```julia
-julia> method, report = fdm([-2, 0, 5], 1, report=true)
-(FiniteDifferences.method, FiniteDifferencesReport:
+julia> method = FiniteDifferenceMethod([-2, 0, 5], 1)
+FiniteDifferenceMethod:
   order of method:       3
   order of derivative:   1
   grid:                  [-2, 0, 5]
-  coefficients:          [-0.357143, 0.3, 0.0571429]
-  roundoff error:        2.22e-16
-  bounds on derivatives: 1.00e+00
-  step size:             3.62e-06
-  accuracy:              6.57e-11
-)
+  coefficients:          [-0.35714285714285715, 0.3, 0.05714285714285714]
 
 julia> method(sin, 1) - cos(1)
--2.05648831297367e-11
+-8.423706177040913e-11
 ```
 
 Compute a directional derivative:
@@ -64,10 +68,78 @@ Compute a directional derivative:
 julia> f(x) = sum(x)
 f (generic function with 1 method)
 
-julia> central_fdm(5, 1)(ε -> f([1, 1, 1] + ε * [1, 2, 3]), 0) - 6
--2.922107000813412e-13
+julia> central_fdm(5, 1)(ε -> f([1, 1, 1] .+ ε .* [1, 2, 3]), 0) - 6
+-6.217248937900877e-15
 ```
 
-## FDM.jl
+## Example: Multivariate Derivatives
 
-This package was formerly called FDM.jl. We recommend users of FDM.jl [update to FiniteDifferences.jl](https://github.com/JuliaDiff/FiniteDifferences.jl/issues/37).
+Consider a quadratic function:
+
+```julia
+julia> a = randn(3, 3); a = a * a'
+3×3 Array{Float64,2}:
+  8.0663   -1.12965   1.68556
+ -1.12965   3.55005  -3.10405
+  1.68556  -3.10405   3.77251
+
+julia> f(x) = 0.5 * x' * a * x
+```
+
+Compute the gradient:
+
+```julia
+julia> grad(central_fdm(5, 1), f, x)[1] - a * x
+3-element Array{Float64,1}:
+ -1.2612133559741778e-12
+ -3.526068326209497e-13
+ -2.3305801732931286e-12
+```
+
+Compute the Jacobian:
+
+```julia
+julia> jacobian(central_fdm(5, 1), f, x)[1] - (a * x)'
+1×3 Array{Float64,2}:
+ -1.26121e-12  -3.52607e-13  -2.33058e-12
+```
+
+The Jacobian can also be computed for non-scalar functions:
+
+```julia
+julia> a = randn(3, 3)
+3×3 Array{Float64,2}:
+ -0.343783   1.5708     0.723289
+ -0.425706  -0.478881  -0.306881
+  1.27326   -0.171606   2.23671
+
+julia> f(x) = a * x
+
+julia> jacobian(central_fdm(5, 1), f, x)[1] - a
+3×3 Array{Float64,2}:
+ -2.81331e-13   2.77556e-13  1.28342e-13
+ -3.34732e-14  -6.05072e-15  6.05072e-15
+ -2.24709e-13   1.88821e-13  1.06581e-14
+```
+
+To compute Jacobian--vector products, use `jvp` and `j′vp`:
+
+```julia
+julia> v = randn(3)
+3-element Array{Float64,1}:
+ -1.290782164377614
+ -0.37701592844250903
+ -1.4288108966380777
+
+julia> jvp(central_fdm(5, 1), f, (x, v)) - a * v
+3-element Array{Float64,1}:
+ -1.3233858453531866e-13
+  9.547918011776346e-15
+  3.632649736573512e-13
+
+julia> j′vp(central_fdm(5, 1), f, x, v)[1] - a'x
+ 3-element Array{Float64,1}:
+  3.5704772471945034e-13
+  4.04121180963557e-13
+  1.2807532812075806e-12
+```

docs/Manifest.toml

@@ -30,10 +30,10 @@ uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 version = "0.25.1"
 
 [[FiniteDifferences]]
-deps = ["ChainRulesCore", "LinearAlgebra", "Printf", "Random"]
+deps = ["ChainRulesCore", "LinearAlgebra", "Printf", "Random", "Richardson"]
 path = ".."
 uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
-version = "0.10.7"
+version = "0.10.8"
 
 [[InteractiveUtils]]
 deps = ["Markdown"]
@@ -93,6 +93,12 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[[Richardson]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "74d2cf4de9eda38175c3f94bd94c755a023d5623"
+uuid = "708f8203-808e-40c0-ba2d-98a6953ed40d"
+version = "1.1.0"
+
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 

docs/src/pages/api.md

@@ -2,10 +2,12 @@
 
 ```@docs
 FiniteDifferenceMethod
-fdm
-backward_fdm
-central_fdm
+FiniteDifferenceMethod(::AbstractVector, ::Int; ::Int)
+FiniteDifferences.estimate_step
 forward_fdm
+central_fdm
+backward_fdm
+extrapolate_fdm
 assert_approx_equal
 FiniteDifferences.DEFAULT_CONDITION
 ```

src/FiniteDifferences.jl

@@ -4,6 +4,7 @@ module FiniteDifferences
     using LinearAlgebra
     using Printf
     using Random
+    using Richardson
 
     export to_vec, grad, jacobian, jvp, j′vp