diff --git a/Project.toml b/Project.toml index 1effa854..11971be4 100644 --- a/Project.toml +++ b/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" diff --git a/README.md b/README.md index 75ba06ad..0e380f89 100644 --- a/README.md +++ b/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 +``` diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 88f82914..8a54eae9 100644 --- a/docs/Manifest.toml +++ b/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" diff --git a/docs/src/pages/api.md b/docs/src/pages/api.md index f4d051d0..cf70b88f 100644 --- a/docs/src/pages/api.md +++ b/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 ``` diff --git a/src/FiniteDifferences.jl b/src/FiniteDifferences.jl index 2cf2b3d3..bc67c23d 100644 --- a/src/FiniteDifferences.jl +++ b/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 diff --git a/src/methods.jl b/src/methods.jl index cae1c4be..2fb6a7f3 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -1,4 +1,13 @@ -export FiniteDifferenceMethod, fdm, backward_fdm, forward_fdm, central_fdm +export FiniteDifferenceMethod, fdm, backward_fdm, forward_fdm, central_fdm, extrapolate_fdm + +""" + add_tiny(x::Union{AbstractFloat, Integer}) + +Add a tiny number, 10^{-40}, to a real floating point number `x`, preserving the type. If +`x` is an `Integer`, it is promoted to a suitable floating point type. +""" +add_tiny(x::T) where T<:AbstractFloat = x + convert(T, 1e-40) +add_tiny(x::Integer) = add_tiny(float(x)) """ FiniteDifferences.DEFAULT_CONDITION @@ -10,160 +19,187 @@ derivatives. const DEFAULT_CONDITION = 100 """ - add_tiny(x::Real) + FiniteDifferenceMethod{G<:AbstractVector, C<:AbstractVector, E<:Function} -Add a tiny number, 10^{-20}, to `x`, preserving the type. If `x` is an `Integer`, it is -promoted to a suitable floating point type. -""" -add_tiny(x::T) where {T<:Real} = x + convert(T, 1e-20) -add_tiny(x::Integer) = add_tiny(float(x)) +A finite difference method. -forward_grid(p::Int) = 0:(p - 1) -backward_grid(p::Int) = (1 - p):0 -function central_grid(p::Int) - if isodd(p) - return div(1 - p, 2):div(p - 1, 2) - else - return vcat(div(-p, 2):-1, 1:div(p, 2)) - end +# Fields +- `grid::G`: Multiples of the step size that the function will be evaluated at. +- `q::Int`: Order of derivative to estimate. +- `coefs::C`: Coefficients corresponding to the grid functions that the function evaluations + will be weighted by. +- `bound_estimator::Function`: A function that takes in the function and the evaluation + point and returns a bound on the magnitude of the `length(grid)`th derivative. +""" +struct FiniteDifferenceMethod{G<:AbstractVector, C<:AbstractVector, E<:Function} + grid::G + q::Int + coefs::C + bound_estimator::E end """ - History + FiniteDifferenceMethod( + grid::AbstractVector, + q::Int; + condition::Real=DEFAULT_CONDITION + ) + +Construct a finite difference method. + +# Arguments +- `grid::Abstract`: The grid. See [`FiniteDifferenceMethod`](@ref). +- `q::Int`: Order of the derivative to estimate. +- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref). -A mutable type that tracks several values during adaptive bound computation. +# Returns +- `FiniteDifferenceMethod`: Specified finite difference method. """ -mutable struct History - adapt::Int - eps::Real - bound::Real - step::Real - accuracy::Real - - function History(; kwargs...) - h = new() - for (k, v) in kwargs - setfield!(h, k, v) - end - return h - end +function FiniteDifferenceMethod( + grid::AbstractVector, + q::Int; + condition::Real=DEFAULT_CONDITION +) + p = length(grid) + _check_p_q(p, q) + return FiniteDifferenceMethod( + grid, + q, + _coefs(grid, q), + _make_default_bound_estimator(condition=condition) + ) end """ - FiniteDifferenceMethod + (m::FiniteDifferenceMethod)( + f::Function, + x::T; + factor::Real=1, + max_step::Real=0.1 * max(abs(x), one(x)) + ) where T<:AbstractFloat -Abstract type for all finite differencing method types. -Subtypes of `FiniteDifferenceMethod` are callable with the signature +Estimate the derivative of `f` at `x` using the finite differencing method `m` and an +automatically determined step size. -``` -method(f, x; kwargs...) -``` +# Arguments +- `f::Function`: Function to estimate derivative of. +- `x::T`: Input to estimate derivative at. -where the keyword arguments can be any of +# Keywords +- `factor::Real=1`: Factor to amplify the estimated round-off error by. This can be used + to force a more conservative step size. +- `max_step::Real=0.1 * max(abs(x), one(x))`: Maximum step size. -* `adapt`: The number of adaptive steps to use improve the estimate of `bound`. -* `bound`: Bound on the value of the function and its derivatives at `x`. -* `condition`: The condition number. See [`DEFAULT_CONDITION`](@ref). -* `eps`: The assumed roundoff error. Defaults to `add_tiny(eps())`. -""" -abstract type FiniteDifferenceMethod end +# Returns +- Estimate of the derivative. -function Base.show(io::IO, x::FiniteDifferenceMethod) - @printf io "FiniteDifferenceMethod:\n" - @printf io " order of method: %d\n" x.p - @printf io " order of derivative: %d\n" x.q - @printf io " grid: %s\n" x.grid - @printf io " coefficients: %s\n" x.coefs - h = x.history - if all(p->isdefined(h, p), propertynames(h)) - @printf io " roundoff error: %.2e\n" h.eps - @printf io " bounds on derivatives: %.2e\n" h.bound - @printf io " step size: %.2e\n" h.step - @printf io " accuracy: %.2e\n" h.accuracy - end -end +# Examples -for D in (:Forward, :Backward, :Central, :Nonstandard) - @eval begin - struct $D{G<:AbstractVector, C<:AbstractVector} <: FiniteDifferenceMethod - p::Int - q::Int - grid::G - coefs::C - history::History - end - (d::$D)(f, x=0.0; kwargs...) = fdm(d, f, x; kwargs...) - end -end +```julia-repl +julia> fdm = central_fdm(5, 1) +FiniteDifferenceMethod: + order of method: 5 + order of derivative: 1 + grid: [-2, -1, 0, 1, 2] + coefficients: [0.08333333333333333, -0.6666666666666666, 0.0, 0.6666666666666666, -0.08333333333333333] -# The below does not apply to Nonstandard, as it has its own constructor -for D in (:Forward, :Backward, :Central) - lcname = lowercase(String(D)) - gridf = Symbol(lcname, "_grid") - fdmf = Symbol(lcname, "_fdm") - - @eval begin - # Compatibility layer over the "old" API - function $fdmf(p::Integer, q::Integer; adapt=1, kwargs...) - _dep_kwarg(kwargs) - return $D(p, q; adapt=adapt, kwargs...) - end +julia> fdm(sin, 1) +0.5403023058681155 - function $D(p::Integer, q::Integer; adapt=1, kwargs...) - _check_p_q(p, q) - grid = $gridf(p) - coefs = _coefs(grid, p, q) - hist = History(; adapt=adapt, kwargs...) - return $D{typeof(grid), typeof(coefs)}(Int(p), Int(q), grid, coefs, hist) - end +julia> fdm(sin, 1) - cos(1) # Check the error. +-2.4313884239290928e-14 - @doc """ - FiniteDifferences.$($(Meta.quot(D)))(p, q; kwargs...) - $($(Meta.quot(fdmf)))(p, q; kwargs...) - - Construct a $($lcname) finite difference method of order `p` to compute the `q`th - derivative. - See [`FiniteDifferenceMethod`](@ref) for more details. - """ - ($D, $fdmf) - end +julia> FiniteDifferences.estimate_step(fdm, sin, 1.0) # Computes step size and estimates the error. +(0.0010632902144695163, 1.9577610541734626e-13) +``` +""" +@inline function (m::FiniteDifferenceMethod)(f::Function, x::Real; kw_args...) + # Assume that converting to float is desired. + return _call_method(m, f, float(x); kw_args...) +end +@inline function _call_method( + m::FiniteDifferenceMethod, + f::Function, + x::T; + factor::Real=1, + max_step::Real=0.1 * max(abs(x), one(x)) +) where T<:AbstractFloat + # The automatic step size calculation fails if `m.q == 0`, so handle that edge case. + iszero(m.q) && return f(x) + h, _ = estimate_step(m, f, x, factor=factor, max_step=max_step) + return _eval_method(m, f, x, h) end """ - FiniteDifferences.Nonstandard(grid, q; kwargs...) + (m::FiniteDifferenceMethod)(f::Function, x::T, h::Real) where T<:AbstractFloat + +Estimate the derivative of `f` at `x` using the finite differencing method `m` and a given +step size. + +# Arguments +- `f::Function`: Function to estimate derivative of. +- `x::T`: Input to estimate derivative at. +- `h::Real`: Step size. + +# Returns +- Estimate of the derivative. -An finite differencing method which is constructed based on a user-defined grid. It is -nonstandard in the sense that it represents neither forward, backward, nor central -differencing. -See [`FiniteDifferenceMethod`](@ref) for further details. +# Examples + +```julia-repl +julia> fdm = central_fdm(5, 1) +FiniteDifferenceMethod: + order of method: 5 + order of derivative: 1 + grid: [-2, -1, 0, 1, 2] + coefficients: [0.08333333333333333, -0.6666666666666666, 0.0, 0.6666666666666666, -0.08333333333333333] + +julia> fdm(sin, 1, 1e-3) +0.5403023058679624 + +julia> fdm(sin, 1, 1e-3) - cos(1) # Check the error. +-1.7741363933510002e-13 +``` """ -function Nonstandard(grid::AbstractVector{<:Real}, q::Integer; adapt=0, kwargs...) - p = length(grid) - _check_p_q(p, q) - coefs = _coefs(grid, p, q) - hist = History(; adapt=adapt, kwargs...) - return Nonstandard{typeof(grid), typeof(coefs)}(Int(p), Int(q), grid, coefs, hist) +@inline function (m::FiniteDifferenceMethod)(f::Function, x::Real, h::Real) + # Assume that converting to float is desired. + return _eval_method(m, f, float(x), h) +end +@inline function _eval_method( + m::FiniteDifferenceMethod, + f::Function, + x::T, + h::Real +) where T<:AbstractFloat + return sum( + i -> convert(T, m.coefs[i]) * f(T(x + h * m.grid[i])), + eachindex(m.grid) + ) / h^m.q end -# Check the method and derivative orders for consistency +# Check the method and derivative orders for consistency. function _check_p_q(p::Integer, q::Integer) - q >= 0 || throw(ArgumentError("order of derivative must be nonnegative")) - q < p || throw(ArgumentError("order of the method must be strictly greater than that " * - "of the derivative")) - # Check whether the method can be computed. We require the factorial of the - # method order to be computable with regular `Int`s, but `factorial` will overflow - # after 20, so 20 is the largest we can allow. - p > 20 && throw(ArgumentError("order of the method is too large to be computed")) + q >= 0 || throw(DomainError(q, "order of derivative (`q`) must be non-negative")) + q < p || throw(DomainError( + (q, p), + "order of the method (q) must be strictly greater than that of the derivative (p)", + )) + # Check whether the method can be computed. We require the factorial of the method order + # to be computable with regular `Int`s, but `factorial` will after 20, so 20 is the + # largest we can allow. + p > 20 && throw(DomainError(p, "order of the method (`p`) is too large to be computed")) return end -# Compute coefficients for the method +const _COEFFS_CACHE = Dict{Tuple{AbstractVector{<:Real}, Integer}, Vector{Float64}}() -const _COEFFS_CACHE = Dict{Tuple{AbstractVector{<:Real}, Integer, Integer}, Vector{Float64}}() -function _coefs(grid::AbstractVector{<:Real}, p::Integer, q::Integer) - return get!(_COEFFS_CACHE, (grid, p, q)) do - # For high precision on the \ we use Rational, and to prevent overflows we use Int128 - # At the end we go to Float64 for fast floating point math (rather than rational math) +# Compute coefficients for the method and cache the result. +function _coefs(grid::AbstractVector{<:Real}, q::Integer) + return get!(_COEFFS_CACHE, (grid, q)) do + p = length(grid) + # For high precision on the `\`, we use `Rational`, and to prevent overflows we use + # `Int128`. At the end we go to `Float64` for fast floating point math, rather than + # rational math. C = [Rational{Int128}(g^i) for i in 0:(p - 1), g in grid] x = zeros(Rational{Int128}, p) x[q + 1] = factorial(q) @@ -171,157 +207,206 @@ function _coefs(grid::AbstractVector{<:Real}, p::Integer, q::Integer) end end +# Estimate the bound on the derivative by amplifying the ∞-norm. +function _make_default_bound_estimator(; condition::Real=DEFAULT_CONDITION) + default_bound_estimator(f, x) = condition * maximum(abs, f(x)) + return default_bound_estimator +end -# Estimate the bound on the function value and its derivatives at a point -_estimate_bound(x, cond) = add_tiny(cond * maximum(abs, x)) +function Base.show(io::IO, m::MIME"text/plain", x::FiniteDifferenceMethod) + @printf io "FiniteDifferenceMethod:\n" + @printf io " order of method: %d\n" length(x.grid) + @printf io " order of derivative: %d\n" x.q + @printf io " grid: %s\n" x.grid + @printf io " coefficients: %s\n" x.coefs +end """ - fdm(m::FiniteDifferenceMethod, f, x[, Val(false)]; kwargs...) -> Real - fdm(m::FiniteDifferenceMethod, f, x, Val(true); kwargs...) -> Tuple{FiniteDifferenceMethod, Real} - -Compute the derivative of `f` at `x` using the finite differencing method `m`. -The optional `Val` argument dictates whether the method should be returned alongside the -derivative value, which can be useful for examining the step size used and other such -parameters. - -The recognized keywords are: - -* `adapt`: The number of adaptive steps to use improve the estimate of `bound`. -* `bound`: Bound on the value of the function and its derivatives at `x`. -* `condition`: The condition number. See [`DEFAULT_CONDITION`](@ref). -* `eps`: The assumed roundoff error. Defaults to `add_tiny(eps())`. - -!!! warning - Bounds can't be adaptively computed over nonstandard grids; passing a value for - `adapt` greater than 0 when `m::Nonstandard` results in an error. - -!!! note - Calling [`FiniteDifferenceMethod`](@ref) objects is equivalent to passing them to `fdm`. - -# Examples - -```julia-repl -julia> fdm(central_fdm(5, 1), sin, 1; adapt=2) -0.5403023058681039 - -julia> fdm(central_fdm(2, 1), exp, 0, Val(true)) -(FiniteDifferenceMethod: - order of method: 2 - order of derivative: 1 - grid: [-1, 1] - coefficients: [-0.5, 0.5] - roundoff error: 1.42e-14 - bounds on derivatives: 1.00e+02 - step size: 1.69e-08 - accuracy: 1.69e-06 -, 1.0000000031817473) -``` + function estimate_step( + m::FiniteDifferenceMethod, + f::Function, + x::T; + factor::Real=1, + max_step::Real=0.1 * max(abs(x), one(x)) + ) where T<:AbstractFloat + +Estimate the step size for a finite difference method `m`. Also estimates the error of the +estimate of the derivative. + +# Arguments +- `m::FiniteDifferenceMethod`: Finite difference method to estimate the step size for. +- `f::Function`: Function to evaluate the derivative of. +- `x::T`: Point to estimate the derivative at. + +# Keywords +- `factor::Real=1`. Factor to amplify the estimated round-off error by. This can be used + to force a more conservative step size. +- `max_step::Real=0.1 * max(abs(x), one(x))`: Maximum step size. + +# Returns +- `Tuple{T, <:AbstractFloat}`: Estimated step size and an estimate of the error of the + finite difference estimate. """ -function fdm( - m::M, - f, - x::T, - ::Val{true}; - condition=DEFAULT_CONDITION, - bound=_estimate_bound(f(x), condition), - eps=add_tiny(Base.eps(bound)), - adapt=m.history.adapt, - max_step=convert(T, 0.1), -) where {T<:AbstractFloat, M<:FiniteDifferenceMethod} - if M <: Nonstandard && adapt > 0 - throw(ArgumentError("can't adaptively compute bounds over Nonstandard grids")) - end - eps > 0 || throw(ArgumentError("eps must be positive, got $eps")) - bound > 0 || throw(ArgumentError("bound must be positive, got $bound")) - 0 <= adapt < 20 - m.p || throw(ArgumentError("can't perform $adapt adaptation steps")) - - # The below calculation can fail for `m.q == 0`, because then `C₂` may be zero. We - # therefore handle this edge case here. - m.q == 0 && return m, f(x) - - p = m.p +function estimate_step( + m::FiniteDifferenceMethod, + f::Function, + x::T; + factor::Real=1, + max_step::Real=0.1 * max(abs(x), one(x)) +) where T<:AbstractFloat + p = length(m.coefs) q = m.q - grid = m.grid - coefs = m.coefs - - # Adaptively compute the bound on the function and derivative values, if applicable. - if adapt > 0 - newm = (M.name.wrapper)(p + 1, p) - dfdx = fdm( - newm, - f, - x; - condition=condition, - eps=eps, - bound=bound, - max_step=max_step, - adapt=(adapt - 1), - ) - bound = _estimate_bound(dfdx, condition) - end + f_x = float(f(x)) + + # Estimate the bound and round-off error. + ε = add_tiny(maximum(eps, f_x)) * factor + M = add_tiny(m.bound_estimator(f, x)) # Set the step size by minimising an upper bound on the error of the estimate. - C₁ = eps * sum(abs, coefs) - C₂ = bound * sum(n->abs(coefs[n] * grid[n]^p), eachindex(coefs)) / factorial(p) - ĥ = convert(T, min((q / (p - q) * C₁ / C₂)^(1 / p), max_step)) + C₁ = ε * sum(abs, m.coefs) + C₂ = M * sum(n -> abs(m.coefs[n] * m.grid[n]^p), eachindex(m.coefs)) / factorial(p) + # Type inference fails on this, so we annotate it, which gives big performance benefits. + h::T = convert(T, min((q / (p - q) * C₁ / C₂)^(1 / p), max_step)) # Estimate the accuracy of the method. - accuracy = ĥ^(-q) * C₁ + ĥ^(p - q) * C₂ + accuracy = h^(-q) * C₁ + h^(p - q) * C₂ - # Estimate the value of the derivative. - dfdx = sum(i -> convert(T, coefs[i]) * f(T(x + ĥ * grid[i])), eachindex(grid)) / ĥ^q + return h, accuracy +end - m.history.eps = eps - m.history.bound = bound - m.history.step = ĥ - m.history.accuracy = accuracy +for direction in [:forward, :central, :backward] + fdm_fun = Symbol(direction, "_fdm") + grid_fun = Symbol("_", direction, "_grid") + @eval begin function $fdm_fun( + p::Int, + q::Int; + adapt::Int=1, + condition::Real=DEFAULT_CONDITION, + geom::Bool=false + ) + _check_p_q(p, q) + grid = $grid_fun(p) + geom && (grid = _exponentiate_grid(grid)) + coefs = _coefs(grid, q) + return FiniteDifferenceMethod( + grid, + q, + coefs, + _make_adaptive_bound_estimator($fdm_fun, p, adapt, condition, geom=geom), + ) + end - return m, dfdx + @doc """ + $($(Meta.quot(fdm_fun)))( + p::Int, + q::Int; + adapt::Int=1, + condition::Real=DEFAULT_CONDITION, + geom::Bool=false + ) + +Contruct a finite difference method at a $($(Meta.quot(direction))) grid of `p` linearly +spaced points. + +# Arguments +- `p::Int`: Number of grid points. +- `q::Int`: Order of the derivative to estimate. + +# Keywords +- `adapt::Int=1`: Use another finite difference method to estimate the magnitude of the + `p`th order derivative, which is important for the step size computation. Recurse + this procedure `adapt` times. +- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref). +- `geom::Bool`: Use geometrically spaced points instead of linearly spaced points. + +# Returns +- `FiniteDifferenceMethod`: The specified finite difference method. + """ $fdm_fun + end end -# Handle inputs that aren't `AbstractFloat`s -- assume what you wanted was a `Float64`. This -# most common way that this method gets invoked is when `x` is an `Integer`. -function fdm(m::FiniteDifferenceMethod, f, x::Real, val; kwargs...) - return fdm(m, f, Float64(x), val; kwargs...) +function _make_adaptive_bound_estimator( + constructor::Function, + q::Int, + adapt::Int, + condition::Int; + kw_args... +) + if adapt >= 1 + estimate_derivative = constructor( + q + 1, q, adapt=adapt - 1, condition=condition; kw_args... + ) + return (f, x) -> maximum(abs, estimate_derivative(f, x)) + else + return _make_default_bound_estimator(condition=condition) + end end -function fdm(m::FiniteDifferenceMethod, f, x::Real, ::Val{false}=Val(false); kwargs...) - _, dfdx = fdm(m, f, x, Val(true); kwargs...) - return dfdx -end +_forward_grid(p::Int) = collect(0:(p - 1)) +_backward_grid(p::Int) = collect((1 - p):0) -## Deprecations +function _central_grid(p::Int) + if isodd(p) + return collect(div(1 - p, 2):div(p - 1, 2)) + else + return vcat(div(-p, 2):-1, 1:div(p, 2)) + end +end -# Used for keyword argument name deprecations -function _dep_kwarg(kwargs) - for (old, new) in [(:ε, :eps), (:M, :bound)] - haskey(kwargs, old) || continue - val = kwargs[old] - error( - "keyword argument `", old, "` should now be passed as `", new, "` upon ", - "application of the method. For example:\n ", - "central_fdm(5, 1)(f, x; $new=$val)\n", - "not\n ", - "central_fdm(5, 1; $old=$val)(f, x)" - ) +_exponentiate_grid(grid::Vector, base::Int=3) = sign.(grid) .* base .^ abs.(grid) ./ base + +function _is_symmetric(vec::Vector; centre_zero::Bool=false, negate_half::Bool=false) + half_sign = negate_half ? -1 : 1 + if isodd(length(vec)) + centre_zero && vec[end ÷ 2 + 1] != 0 && return false + return vec[1:end ÷ 2] == half_sign .* reverse(vec[(end ÷ 2 + 2):end]) + else + return vec[1:end ÷ 2] == half_sign .* reverse(vec[(end ÷ 2 + 1):end]) end end -function fdm( - grid::AbstractVector{<:Real}, - q::Int, - ::Union{Val{true}, Val{false}}=Val(false); - kwargs..., -) - error("to use a custom grid, use `Nonstandard(grid, q)` and pass the result to `fdm`") +function _is_symmetric(m::FiniteDifferenceMethod) + grid_symmetric = _is_symmetric(m.grid, centre_zero=true, negate_half=true) + coefs_symmetric =_is_symmetric(m.coefs, negate_half=true) + return grid_symmetric && coefs_symmetric end -for fdmf in (:central_fdm, :backward_fdm, :forward_fdm) - @eval function $fdmf(p::Int, q::Int, ::Union{Val{true}, Val{false}}; kwargs...) - error( - "the `Val` argument should now be passed directly to `fdm` after ", - "constructing the method, not to the method constructor itself" - ) - end +""" + extrapolate_fdm( + m::FiniteDifferenceMethod, + f::Function, + x::T, + h::Real=0.1 * max(abs(x), one(x)); + power=nothing, + breaktol=Inf, + kw_args... + ) where T<:AbstractFloat + +Use Richardson extrapolation to refine a finite difference method. + +Takes further in keyword arguments for `Richardson.extrapolate`. This method +automatically sets `power = 2` if `m` is symmetric and `power = 1`. Moreover, it defaults +`breaktol = Inf`. + +# Arguments +- `m::FiniteDifferenceMethod`: Finite difference method to estimate the step size for. +- `f::Function`: Function to evaluate the derivative of. +- `x::T`: Point to estimate the derivative at. +- `h::Real=0.1 * max(abs(x), one(x))`: Initial step size. + +# Returns +- `Tuple{<:AbstractFloat, <:AbstractFloat}`: Estimate of the derivative and error. +""" +function extrapolate_fdm( + m::FiniteDifferenceMethod, + f::Function, + x::T, + h::Real=0.1 * max(abs(x), one(x)); + power::Int=1, + breaktol::Real=Inf, + kw_args... +) where T<:AbstractFloat + (power == 1 && _is_symmetric(m)) && (power = 2) + return extrapolate(h -> m(f, x, h), h; power=power, breaktol=breaktol, kw_args...) end diff --git a/test/grad.jl b/test/grad.jl index 29beee2d..4374553e 100644 --- a/test/grad.jl +++ b/test/grad.jl @@ -57,7 +57,7 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec @test J_fdm ≈ J_exact @test J_fdm == jacobian(fdm, f, x)[1] - # Check that the estimated jvp and j′vp are consistent with their definitions. + # Check that the estimated jvp and j′vp are consistent with their definitions. @test _jvp(fdm, f, x, ẋ) ≈ J_exact * ẋ @test _j′vp(fdm, f, ȳ, x) ≈ transpose(J_exact) * ȳ @@ -85,7 +85,7 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec @testset "multi vars jacobian/grad" begin rng, fdm = MersenneTwister(123456), central_fdm(5, 1) - + f1(x, y) = x * y + x f2(x, y) = sum(x * y + x) f3(x::Tuple) = sum(x[1]) + x[2] @@ -126,7 +126,7 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec x, y = rand(rng, 3, 3), 2.0 dxs = grad(fdm, f3, (x, y))[1] @test dxs[1] ≈ grad(fdm, x->f3((x, y)), x)[1] - @test dxs[2] ≈ grad(fdm, y->f3((x, y)), y)[1] + @test dxs[2] ≈ grad(fdm, y->f3((x, y)), y)[1] end @testset "check dict" begin @@ -160,7 +160,7 @@ using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec rng = MersenneTwister(123456) x = randn(rng, T) y = randn(rng, T) - fdm = FiniteDifferences.Central(5, 1) + fdm = FiniteDifferences.central_fdm(5, 1) # Addition. dx, dy = FiniteDifferences.jacobian(fdm, +, x, y) diff --git a/test/methods.jl b/test/methods.jl index 317c21d1..c938e615 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -1,4 +1,4 @@ -using FiniteDifferences: Forward, Backward, Central, Nonstandard, add_tiny +using FiniteDifferences: add_tiny @testset "Methods" begin @@ -31,33 +31,32 @@ using FiniteDifferences: Forward, Backward, Central, Nonstandard, add_tiny # Test all combinations of the above settings, i.e. differentiate all functions using # all methods and data types. @testset "foo=$(foo.f), method=$m, type=$(T)" for foo in foos, m in methods, T in types - - @testset "method-order=$order" for order in [1, 2, 3] - @test m(order, 0; adapt=2)(foo.f, T(1)) isa T - @test m(order, 0; adapt=2)(foo.f, T(1)) == T(foo.f(1)) + @testset "method-order=$order" for order in [1, 2, 3, 4, 5] + @test m(order, 0, adapt=2)(foo.f, T(1)) isa T + @test m(order, 0, adapt=2)(foo.f, T(1)) == T(foo.f(1)) end @test m(10, 1)(foo.f, T(1)) isa T @test m(10, 1)(foo.f, T(1)) ≈ T(foo.d1) - @test m(10, 2; bound=1)(foo.f, T(1)) isa T + @test m(10, 2)(foo.f, T(1)) isa T if T == Float64 - @test m(10, 2; bound=1)(foo.f, T(1)) ≈ T(foo.d2) + @test m(10, 2)(foo.f, T(1)) ≈ T(foo.d2) end end # Integration test to ensure that Integer-output functions can be tested. - @testset "Integer Output" begin + @testset "Integer output" begin @test isapprox(central_fdm(5, 1)(x -> 5, 0), 0; rtol=1e-12, atol=1e-12) end @testset "Adaptation improves estimate" begin - @test forward_fdm(5, 1)(log, 0.001; adapt=0) ≈ 969.2571703 - @test forward_fdm(5, 1)(log, 0.001; adapt=1) ≈ 1000 + @test forward_fdm(5, 1, adapt=0)(log, 0.001) ≈ 997.077814 + @test forward_fdm(5, 1, adapt=1)(log, 0.001) ≈ 1000 end @testset "Limiting step size" begin - @test !isfinite(central_fdm(5, 1)(abs, 0.001; max_step=0)) + @test !isfinite(central_fdm(5, 1)(abs, 0.001, max_step=0)) @test central_fdm(5, 1)(abs, 0.001) ≈ 1.0 end @@ -65,62 +64,71 @@ using FiniteDifferences: Forward, Backward, Central, Nonstandard, add_tiny # Regression test against issues with precision during computation of _coeffs # see https://github.com/JuliaDiff/FiniteDifferences.jl/issues/64 - @test fdm(central_fdm(9, 5), exp, 1.0, adapt=4) ≈ exp(1) atol=1e-7 + @test central_fdm(9, 5, adapt=4, condition=1)(exp, 1.0) ≈ exp(1) atol=1e-7 poly(x) = 4x^3 + 3x^2 + 2x + 1 - @test fdm(central_fdm(9, 3), poly, 1.0, adapt=4) ≈ 24 atol=1e-11 + @test central_fdm(9, 3, adapt=4)(poly, 1.0) ≈ 24 atol=1e-11 end - @testset "Printing FiniteDifferenceMethods" begin - @test sprint(show, central_fdm(2, 1)) == """ + @test sprint(show, "text/plain", central_fdm(2, 1)) == """ FiniteDifferenceMethod: order of method: 2 order of derivative: 1 grid: [-1, 1] coefficients: [-0.5, 0.5] """ - m, _ = fdm(central_fdm(2, 1), sin, 1, Val(true)) - report = sprint(show, m) - regex_float = r"[\d\.\+-e]+" - regex_array = r"\[([\d.+-e]+(, )?)+\]" - @test occursin(Regex(join(map(x -> x.pattern, - [ - r"FiniteDifferenceMethod:", - r"order of method:", r"\d+", - r"order of derivative:", r"\d+", - r"grid:", regex_array, - r"coefficients:", regex_array, - r"roundoff error:", regex_float, - r"bounds on derivatives:", regex_float, - r"step size:", regex_float, - r"accuracy:", regex_float, - r"" - ] - ), r"\s*".pattern)), report) end - @testset "Breaking deprecations" begin - @test_throws ErrorException fdm([1,2,3], 4) # Custom grids need Nonstandard - for f in (forward_fdm, backward_fdm, central_fdm) - @test_throws ErrorException f(2, 1; M=1) # Old kwarg, now misplaced - @test_throws ErrorException f(2, 1, Val(true)) # Ask fdm for reports instead + @testset "_is_symmetric" begin + # Test odd grids: + @test FiniteDifferences._is_symmetric([2, 1, 0, 1, 2]) + @test !FiniteDifferences._is_symmetric([2, 1, 0, 3, 2]) + @test !FiniteDifferences._is_symmetric([4, 1, 0, 1, 2]) + + # Test even grids: + @test FiniteDifferences._is_symmetric([2, 1, 1, 2]) + @test !FiniteDifferences._is_symmetric([2, 1, 3, 2]) + @test !FiniteDifferences._is_symmetric([4, 1, 1, 2]) + + # Test zero at centre: + @test FiniteDifferences._is_symmetric([2, 1, 4, 1, 2]) + @test !FiniteDifferences._is_symmetric([2, 1, 4, 1, 2], centre_zero=true) + @test FiniteDifferences._is_symmetric([2, 1, 1, 2], centre_zero=true) + + # Test negation of a half: + @test !FiniteDifferences._is_symmetric([2, 1, -1, -2]) + @test FiniteDifferences._is_symmetric([2, 1, -1, -2], negate_half=true) + @test FiniteDifferences._is_symmetric([2, 1, 0, -1, -2], negate_half=true) + @test FiniteDifferences._is_symmetric([2, 1, 4, -1, -2], negate_half=true) + @test !FiniteDifferences._is_symmetric( + [2, 1, 4, -1, -2], + negate_half=true, + centre_zero=true + ) + + # Test symmetry of `central_fdm`. + for p in 2:10 + m = central_fdm(p, 1) + @test FiniteDifferences._is_symmetric(m) end - end - @testset "Types" begin - @testset "$T" for T in (Forward, Backward, Central) - @test T(5, 1)(sin, 1; adapt=4) ≈ cos(1) - @test_throws ArgumentError T(3, 3) - @test_throws ArgumentError T(3, 4) - @test_throws ArgumentError T(40, 5) - @test_throws ArgumentError T(5, 1)(sin, 1; adapt=200) - @test_throws ArgumentError T(5, 1)(sin, 1; eps=0.0) - @test_throws ArgumentError T(5, 1)(sin, 1; bound=0.0) + # Test asymmetry of `forward_fdm` and `backward_fdm`. + for p in 2:10 + for f in [forward_fdm, backward_fdm] + m = f(p, 1) + @test !FiniteDifferences._is_symmetric(m) + end end - @testset "Nonstandard" begin - @test Nonstandard([-2, -1, 1], 1)(sin, 1) ≈ cos(1) - @test_throws ArgumentError Nonstandard([-2, -1, 1], 1)(sin, 1; adapt=2) + end + + @testset "extrapolate_fdm" begin + # Also test an `Integer` argument as input. + for x in [1, 1.0] + for f in [forward_fdm, central_fdm, backward_fdm] + estimate, _ = extrapolate_fdm(f(4, 3), exp, x, contract=0.8) + @test estimate ≈ exp(1.0) atol=1e-7 + end end end end