From 0603acadf297f005541e22232f660492b1c9ba98 Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Thu, 16 May 2019 16:08:58 -0700 Subject: [PATCH 1/6] Add adaptive function bound computation This introduces a type-based approach that works slightly differently than the current API. Upon construction of a `FDMethod` subtype, the order of the method and the order of the derivative are checked for agreement, and the grid and coefficients are computed and stored. Instances of the types are callable, and doing so produces the value of the derivative. It works recursively when adaptively computing the bound, which results in a better step size and thus a (hopefully) more robust estimate of the derivative. --- src/methods.jl | 83 +++++++++++++++++++++++++++++++++++++++++++++++++ test/methods.jl | 10 ++++++ 2 files changed, 93 insertions(+) diff --git a/src/methods.jl b/src/methods.jl index 451aa925..447a8d2a 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -1,5 +1,88 @@ export FDMReport, fdm, backward_fdm, forward_fdm, central_fdm +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 +end + +abstract type FDMethod end + +for D in (:Forward, :Backward, :Central) + gridf = Symbol(lowercase(String(D)), "_grid") + @eval begin + struct $D{G<:AbstractVector, C<:AbstractVector} <: FDMethod + p::Int + q::Int + grid::G + coefs::C + end + function $D(p::Integer, q::Integer) + 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")) + + grid = $gridf(p) + + C = [g^i for i in 0:(p - 1), g in grid] + x = zeros(Int, p) + x[q + 1] = factorial(q) + coefs = C \ x + + return $D{typeof(grid), typeof(coefs)}(p, q, grid, coefs) + end + (d::$D)(f, x; kwargs...) = fdm(d, f, x; kwargs...) + end +end + +maxabs(x) = maximum(abs, x) + +function fdm( + m::M, + f, + x; + bound=maxabs(f(x)), + eps=Base.eps(bound), + adapt=1, + report=false, +) where M<:FDMethod + p = m.p + q = m.q + grid = m.grid + coefs = m.coefs + + # Adaptively compute the bound on the function and derivative values, if applicable. + if 0 < adapt < p + newm = (M.name.wrapper)(p + 1, q) + dfdx, = fdm(newm, f, x; eps=eps, bound=bound, adapt=(adapt - 1)) + bound = maxabs(dfdx) + end + + # 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) + ĥ = (q / (p - q) * C₁ / C₂) ^ (1 / p) + + # Estimate the accuracy of the method. + accuracy = ĥ^(-q) * C₁ + ĥ^(p - q) * C₂ + + # Estimate the value of the derivative. + dfdx = sum(i->coefs[i] * f(x + ĥ * grid[i]), eachindex(grid)) / ĥ^q + + if report + return dfdx, FDMReport(p, q, grid, coefs, eps, bound, ĥ, accuracy) + else + return dfdx + end +end + """ FDMReport diff --git a/test/methods.jl b/test/methods.jl index 22ff57e1..6fb3f2f3 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -1,3 +1,5 @@ +using FDM: Forward, Backward, Central + @testset "Methods" begin for f in [:forward_fdm, :backward_fdm, :central_fdm] @eval @test $f(10, 1; M=1)(sin, 1) ≈ cos(1) @@ -41,4 +43,12 @@ @test_throws ArgumentError fdm([1,2,3], 4) @test_throws ArgumentError fdm(zeros(40), 5) 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, 4) + @test_throws ArgumentError T(40, 5) + end + end end From a2511fee3fabe2abff3d2842965b702f9f71873c Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Fri, 17 May 2019 11:13:54 -0700 Subject: [PATCH 2/6] WIP: Add more things --- src/methods.jl | 29 +++++++++++++++++++++++++---- test/methods.jl | 3 +++ 2 files changed, 28 insertions(+), 4 deletions(-) diff --git a/src/methods.jl b/src/methods.jl index 447a8d2a..d8ad76ef 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -1,5 +1,21 @@ export FDMReport, fdm, backward_fdm, forward_fdm, central_fdm +""" + FDM.DEFAULT_CONDITION + +The default [condition number](https://en.wikipedia.org/wiki/Condition_number) used when +computing bounds. It provides amplification of the ∞-norm when passed to the function's +derivatives. +""" +const DEFAULT_CONDITION = 100 + +""" + FDM.TINY + +A tiny number added to some quantities to ensure that division by 0 does not occur. +""" +const TINY = 1e-20 + forward_grid(p::Int) = 0:(p - 1) backward_grid(p::Int) = (1 - p):0 function central_grid(p::Int) @@ -42,17 +58,22 @@ for D in (:Forward, :Backward, :Central) end end -maxabs(x) = maximum(abs, x) +estimate_bound(x, cond) = cond * maximum(abs, x) + TINY function fdm( m::M, f, x; - bound=maxabs(f(x)), + condition=DEFAULT_CONDITION, + bound=estimate_bound(f(x), condition), eps=Base.eps(bound), adapt=1, report=false, ) where M<:FDMethod + 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")) + p = m.p q = m.q grid = m.grid @@ -60,9 +81,9 @@ function fdm( # Adaptively compute the bound on the function and derivative values, if applicable. if 0 < adapt < p - newm = (M.name.wrapper)(p + 1, q) + newm = (M.name.wrapper)(p + 1, p) dfdx, = fdm(newm, f, x; eps=eps, bound=bound, adapt=(adapt - 1)) - bound = maxabs(dfdx) + bound = estimate_bound(dfdx, condition) end # Set the step size by minimising an upper bound on the error of the estimate. diff --git a/test/methods.jl b/test/methods.jl index 6fb3f2f3..eb45b841 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -49,6 +49,9 @@ using FDM: Forward, Backward, Central @test T(5, 1)(sin, 1; adapt=4) ≈ cos(1) @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) end end end From dbdfe2af0958d30feeced51ebb07791714047cfd Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Fri, 17 May 2019 18:55:44 -0700 Subject: [PATCH 3/6] WIP: Store computational history in a History object --- src/methods.jl | 43 +++++++++++++++++++++++++++++++------------ 1 file changed, 31 insertions(+), 12 deletions(-) diff --git a/src/methods.jl b/src/methods.jl index d8ad76ef..89f68b64 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -26,6 +26,22 @@ function central_grid(p::Int) end end +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 +end + abstract type FDMethod end for D in (:Forward, :Backward, :Central) @@ -36,8 +52,9 @@ for D in (:Forward, :Backward, :Central) q::Int grid::G coefs::C + history::History end - function $D(p::Integer, q::Integer) + function $D(p::Integer, q::Integer; adapt=1, kwargs...) 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 @@ -52,7 +69,9 @@ for D in (:Forward, :Backward, :Central) x[q + 1] = factorial(q) coefs = C \ x - return $D{typeof(grid), typeof(coefs)}(p, q, grid, coefs) + hist = History(; adapt=adapt, kwargs...) + + return $D{typeof(grid), typeof(coefs)}(p, q, grid, coefs, hist) end (d::$D)(f, x; kwargs...) = fdm(d, f, x; kwargs...) end @@ -66,9 +85,8 @@ function fdm( x; condition=DEFAULT_CONDITION, bound=estimate_bound(f(x), condition), - eps=Base.eps(bound), - adapt=1, - report=false, + eps=(Base.eps(bound) + TINY), + adapt=m.history.adapt, ) where M<:FDMethod eps > 0 || throw(ArgumentError("eps must be positive, got $eps")) bound > 0 || throw(ArgumentError("bound must be positive, got $bound")) @@ -80,9 +98,9 @@ function fdm( coefs = m.coefs # Adaptively compute the bound on the function and derivative values, if applicable. - if 0 < adapt < p + if adapt > 0 newm = (M.name.wrapper)(p + 1, p) - dfdx, = fdm(newm, f, x; eps=eps, bound=bound, adapt=(adapt - 1)) + dfdx = fdm(newm, f, x; condition=condition, eps=eps, bound=bound, adapt=(adapt - 1)) bound = estimate_bound(dfdx, condition) end @@ -97,11 +115,12 @@ function fdm( # Estimate the value of the derivative. dfdx = sum(i->coefs[i] * f(x + ĥ * grid[i]), eachindex(grid)) / ĥ^q - if report - return dfdx, FDMReport(p, q, grid, coefs, eps, bound, ĥ, accuracy) - else - return dfdx - end + m.history.eps = eps + m.history.bound = bound + m.history.step = ĥ + m.history.accuracy = accuracy + + return dfdx end """ From 2638f4766cc7f54ff0b88ba369a71069e63dcbaf Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Mon, 27 May 2019 14:18:27 -0700 Subject: [PATCH 4/6] Restructure to use the new API, add deprecations This change consistutes a near rewrite of method handling. --- src/methods.jl | 353 ++++++++++++++++++++++++++---------------------- test/methods.jl | 88 +++++++----- 2 files changed, 240 insertions(+), 201 deletions(-) diff --git a/src/methods.jl b/src/methods.jl index 89f68b64..1eb33b35 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -1,4 +1,4 @@ -export FDMReport, fdm, backward_fdm, forward_fdm, central_fdm +export fdm, backward_fdm, forward_fdm, central_fdm """ FDM.DEFAULT_CONDITION @@ -26,6 +26,11 @@ function central_grid(p::Int) end end +""" + History + +A mutable type that tracks several values during adaptive bound computation. +""" mutable struct History adapt::Int eps::Real @@ -42,10 +47,41 @@ mutable struct History end end +""" + FDMethod + +Abstract type for all finite differencing method types. +Subtypes of `FDMethod` are callable with the signature + +``` +method(f, x; kwargs...) +``` + +where the keyword arguments can be any of + +* `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 `eps()` plus [`TINY`](@ref). +""" abstract type FDMethod end -for D in (:Forward, :Backward, :Central) - gridf = Symbol(lowercase(String(D)), "_grid") +function Base.show(io::IO, x::FDMethod) + @printf io "FDMethod:\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 + +for D in (:Forward, :Backward, :Central, :Nonstandard) @eval begin struct $D{G<:AbstractVector, C<:AbstractVector} <: FDMethod p::Int @@ -54,40 +90,136 @@ for D in (:Forward, :Backward, :Central) coefs::C history::History end - function $D(p::Integer, q::Integer; adapt=1, kwargs...) - 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")) + (d::$D)(f, x; kwargs...) = fdm(d, f, x; kwargs...) + end +end - grid = $gridf(p) +# 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") - C = [g^i for i in 0:(p - 1), g in grid] - x = zeros(Int, p) - x[q + 1] = factorial(q) - coefs = C \ x + @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 + 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)}(p, q, grid, coefs, hist) + return $D{typeof(grid), typeof(coefs)}(Int(p), Int(q), grid, coefs, hist) end - (d::$D)(f, x; kwargs...) = fdm(d, f, x; kwargs...) + + @doc """ + FDM.$($(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 [`FDMethod`](@ref) for more details. + """ + ($D, $fdmf) end end -estimate_bound(x, cond) = cond * maximum(abs, x) + TINY +""" + FDM.Nonstandard(grid, q; kwargs...) + +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 [`FDMethod`](@ref) for further details. +""" +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) +end + +# Check the method and derivative orders for consistency +function _check_p_q(p::Integer, q::Integer) + 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")) + return +end +# Compute coefficients for the method +function _coefs(grid::AbstractVector{<:Real}, p::Integer, q::Integer) + C = [g^i for i in 0:(p - 1), g in grid] + x = zeros(Int, p) + x[q + 1] = factorial(q) + return C \ x +end + +# Estimate the bound on the function value and its derivatives at a point +_estimate_bound(x, cond) = cond * maximum(abs, x) + TINY + +""" + fdm(m::FDMethod, f, x[, Val(false)]; kwargs...) -> Real + fdm(m::FDMethod, f, x, Val(true); kwargs...) -> Tuple{FDMethod, 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 `eps()` plus [`TINY`](@ref). + +!!! 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 [`FDMethod`](@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)) +(FDMethod: + 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 fdm( m::M, f, - x; + x, + ::Val{true}; condition=DEFAULT_CONDITION, - bound=estimate_bound(f(x), condition), - eps=(Base.eps(bound) + TINY), + bound=_estimate_bound(f(x), condition), + eps=(Base.eps(float(bound)) + TINY), adapt=m.history.adapt, ) where M<:FDMethod + 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")) @@ -101,7 +233,7 @@ function fdm( if adapt > 0 newm = (M.name.wrapper)(p + 1, p) dfdx = fdm(newm, f, x; condition=condition, eps=eps, bound=bound, adapt=(adapt - 1)) - bound = estimate_bound(dfdx, condition) + bound = _estimate_bound(dfdx, condition) end # Set the step size by minimising an upper bound on the error of the estimate. @@ -120,155 +252,46 @@ function fdm( m.history.step = ĥ m.history.accuracy = accuracy - return dfdx + return m, dfdx end -""" - FDMReport - -Details of a finite difference method to estimate a derivative. Instances of `FDMReport` -`Base.show` nicely. - -# Fields -- `p::Int`: Order of the method. -- `q::Int`: Order of the derivative that is estimated. -- `grid::Vector{<:Real}`: Relative spacing of samples of `f` that are used by the method. -- `coefs::Vector{<:Real}`: Weights of the samples of `f`. -- `ε::Real`: Absolute roundoff error of the function evaluations. -- `M::Real`: Assumed upper bound of `f` and all its derivatives at `x`. -- `ĥ::Real`: Step size. -- `err::Real`: Estimated absolute accuracy. -""" -struct FDMReport{Tε, TM, Tĥ, Tacc} - p::Int - q::Int - grid::Vector{<:Real} - coefs::Vector{<:Real} - ε::Tε - M::TM - ĥ::Tĥ - acc::Tacc -end -function Base.show(io::IO, x::FDMReport) - @printf io "FDMReport:\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 - @printf io " roundoff error: %.2e\n" x.ε - @printf io " bounds on derivatives: %.2e\n" x.M - @printf io " step size: %.2e\n" x.ĥ - @printf io " accuracy: %.2e\n" x.acc -end - -""" - function fdm( - grid::Vector{<:Real}, - q::Int; - ε::Real=eps(), - M::Real=1, - report::Bool=false - ) - -Construct a function `method(f, x::Real, h::Real=ĥ)` that takes in a -function `f`, a point `x` in the domain of `f`, and optionally a step size `h`, and -estimates the `q`'th order derivative of `f` at `x` with a `length(grid)`'th order -finite difference method. - -# Arguments -- `grid::Vector{<:Real}`: Relative spacing of samples of `f` that are used by the method. - The length of `grid` determines the order of the method. -- `q::Int`: Order of the derivative to estimate. `q` must be strictly less than the order - of the method. - -# Keywords -- `ε::Real=eps()`: Absolute roundoff error on the function evaluations. -- `M::Real=1`: Upper bound on `f` and all its derivatives. -- `::Bool=false`: Also return an instance of `FDMReport` containing information - about the method constructed. -""" -function fdm(grid::AbstractVector{T}, q::Int, ::Val{false}; ε::Real=eps(), M::Real=one(T), - ) where {T<:Real} - method, p, q, grid, coefs, ε, M, ĥ, acc = _fdm(grid, q, ε, M) - return method -end -function fdm(grid::AbstractVector{T}, q::Int, ::Val{true}; ε::Real=eps(), M::Real=one(T), - ) where {T<:Real} - method, p, q, grid, coefs, ε, M, ĥ, acc = _fdm(grid, q, ε, M) - return method, FDMReport(p, q, grid, coefs, ε, M, ĥ, acc) -end -function fdm(grid::AbstractVector{T}, q::Int; ε::Real=eps(), M::Real=one(T)) where T<:Real - return fdm(grid, q, Val(false); ε=ε, M=M) +function fdm(m::FDMethod, f, x, ::Val{false}=Val(false); kwargs...) + _, dfdx = fdm(m, f, x, Val(true); kwargs...) + return dfdx end -function _fdm(grid::AbstractVector{T}, q::Int, ε::Real, M::Real) where T - p = length(grid) # Order of the method. - 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")) - factp = factorial(p) - # Compute the coefficients of the FDM. - C = [grid[n]^i for i in 0:p - 1, n in eachindex(grid)] - x = zeros(Int, p) - x[q+1] = factorial(q) - coefs = C \ x - - # Set the step size by minimising an upper bound on the error of the estimate. - C₁ = ε * sum(abs, coefs) - C₂ = M * sum(n->abs(coefs[n] * grid[n]^p), eachindex(coefs)) / factp - ĥ = (q / (p - q) * C₁ / C₂) ^ (1 / p) - - # Estimate the accuracy of the method. - acc = ĥ^(-q) * C₁ + ĥ^(p - q) * C₂ - - # Construct the FDM. - function method(f, x::Real=zero(T), h::Real=ĥ) - return sum(n->coefs[n] * f(x + h * grid[n]), eachindex(grid)) / h^q +## Deprecations + +# 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)" + ) end - return method, p, q, grid, coefs, ε, M, ĥ, acc 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`") +end -""" - backward_fdm(p::Int, ...) - -Construct a backward finite difference method of order `p`. See `fdm` for further details. - -# Arguments -- `p::Int`: Order of the method. - -Further takes, in the following order, the arguments `q`, `ε`, `M`, and `report` from `fdm`. -""" -backward_fdm(p::Int, args...; kws...) = fdm(1 - p:0, args...; kws...) - -""" - forward_fdm(p::Int, ...) - -Construct a forward finite difference method of order `p`. See `fdm` for further details. - -# Arguments -- `p::Int`: Order of the method. - -Further takes, in the following order, the arguments `q`, `ε`, `M`, and `report` from `fdm`. -""" -forward_fdm(p::Int, args...; kws...) = fdm(0:p - 1, args...; kws...) - -""" - central_fdm(p::Int, ...) - -Construct a central finite difference method of order `p`. See `fdm` for further details. - -# Arguments -- `p::Int`: Order of the method. - -Further takes, in the following order, the arguments `q`, `ε`, `M`, and `report` from `fdm`. -""" -function central_fdm(p::Int, args...; kws...) - return isodd(p) ? fdm(Int(-(p - 1) / 2):Int((p - 1) / 2), args...; kws...) : - fdm([Int(-p/2):-1; 1:Int(p/2)], args...; kws...) +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 end diff --git a/test/methods.jl b/test/methods.jl index eb45b841..2529af45 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -1,47 +1,59 @@ -using FDM: Forward, Backward, Central +using FDM: Forward, Backward, Central, Nonstandard @testset "Methods" begin for f in [:forward_fdm, :backward_fdm, :central_fdm] - @eval @test $f(10, 1; M=1)(sin, 1) ≈ cos(1) - @eval @test $f(10, 2; M=1)(sin, 1) ≈ -sin(1) + @eval @test $f(10, 1; bound=1)(sin, 1) ≈ cos(1) + @eval @test $f(10, 2; bound=1)(sin, 1) ≈ -sin(1) - @eval @test $f(10, 1; M=1)(exp, 1) ≈ exp(1) - @eval @test $f(10, 2; M=1)(exp, 1) ≈ exp(1) + @eval @test $f(10, 1; bound=1)(exp, 1) ≈ exp(1) + @eval @test $f(10, 2; bound=1)(exp, 1) ≈ exp(1) - @eval @test $f(10, 1; M=1)(abs2, 1) ≈ 2 - @eval @test $f(10, 2; M=1)(abs2, 1) ≈ 2 + @eval @test $f(10, 1; bound=1)(abs2, 1) ≈ 2 + @eval @test $f(10, 2; bound=1)(abs2, 1) ≈ 2 - @eval @test $f(10, 1; M=1)(sqrt, 1) ≈ .5 - @eval @test $f(10, 2; M=1)(sqrt, 1) ≈ -.25 + @eval @test $f(10, 1; bound=1)(sqrt, 1) ≈ .5 + @eval @test $f(10, 2; bound=1)(sqrt, 1) ≈ -.25 end - @test_throws ArgumentError central_fdm(100, 1) - - # Test that printing an instance of `FDMReport` contains the information that it should - # contain. - buffer = IOBuffer() - show(buffer, central_fdm(2, 1, Val(true))[2]) - report = String(take!(copy(buffer))) - regex_float = r"[\d\.\+-e]+" - regex_array = r"\[([\d.+-e]+(, )?)+\]" - @test occursin(Regex(join(map(x -> x.pattern, - [ - r"FDMReport:", - 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) - - @testset "Error conditions" begin - @test_throws ArgumentError fdm([1,2,3], 4) - @test_throws ArgumentError fdm(zeros(40), 5) + @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 + end + + @testset "Printing FDMethods" begin + @test sprint(show, central_fdm(2, 1)) == """ + FDMethod: + 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"FDMethod:", + 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 + end end @testset "Types" begin @@ -53,5 +65,9 @@ using FDM: Forward, Backward, Central @test_throws ArgumentError T(5, 1)(sin, 1; eps=0.0) @test_throws ArgumentError T(5, 1)(sin, 1; bound=0.0) 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 end end From 142d50e0ec9d5b8ef6ef55a933163cbc4e4f84c1 Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Tue, 28 May 2019 09:55:09 -0700 Subject: [PATCH 5/6] Add a max_step parameter to cap the step size If it's too big, Bad Things can happen. See e.g. the added test. --- src/methods.jl | 5 +++-- test/methods.jl | 5 +++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/methods.jl b/src/methods.jl index 1eb33b35..6eb8fe8f 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -90,7 +90,7 @@ for D in (:Forward, :Backward, :Central, :Nonstandard) coefs::C history::History end - (d::$D)(f, x; kwargs...) = fdm(d, f, x; kwargs...) + (d::$D)(f, x=0.0; kwargs...) = fdm(d, f, x; kwargs...) end end @@ -216,6 +216,7 @@ function fdm( bound=_estimate_bound(f(x), condition), eps=(Base.eps(float(bound)) + TINY), adapt=m.history.adapt, + max_step=0.1, ) where M<:FDMethod if M <: Nonstandard && adapt > 0 throw(ArgumentError("can't adaptively compute bounds over Nonstandard grids")) @@ -239,7 +240,7 @@ function fdm( # 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) - ĥ = (q / (p - q) * C₁ / C₂) ^ (1 / p) + ĥ = min((q / (p - q) * C₁ / C₂)^(1 / p), max_step) # Estimate the accuracy of the method. accuracy = ĥ^(-q) * C₁ + ĥ^(p - q) * C₂ diff --git a/test/methods.jl b/test/methods.jl index 2529af45..2fad844e 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -20,6 +20,11 @@ using FDM: Forward, Backward, Central, Nonstandard @test forward_fdm(5, 1)(log, 0.001; adapt=1) ≈ 1000 end + @testset "Limiting step size" begin + @test !isfinite(central_fdm(5, 1)(abs, 0.001; max_step=0)) + @test central_fdm(5, 1)(abs, 0.001) ≈ 1.0 + end + @testset "Printing FDMethods" begin @test sprint(show, central_fdm(2, 1)) == """ FDMethod: From cc9bf6325501075e0a1595c5a33f0723c265a4a5 Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Tue, 28 May 2019 10:00:57 -0700 Subject: [PATCH 6/6] Pass max_step during adaptive steps as well --- src/methods.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/methods.jl b/src/methods.jl index 6eb8fe8f..65331201 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -233,7 +233,16 @@ function fdm( # 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, adapt=(adapt - 1)) + dfdx = fdm( + newm, + f, + x; + condition=condition, + eps=eps, + bound=bound, + max_step=max_step, + adapt=(adapt - 1), + ) bound = _estimate_bound(dfdx, condition) end