From 7d0f13fece64bd4877c03e0f2a5136d4b38cde41 Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Sun, 2 Feb 2020 22:33:49 +0000 Subject: [PATCH 1/3] Precompute grids and coeffs --- src/methods.jl | 70 +++++++++++++++++++++++++++++++------------------- 1 file changed, 43 insertions(+), 27 deletions(-) diff --git a/src/methods.jl b/src/methods.jl index 20b56ea9..2f6ca727 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -26,6 +26,31 @@ function central_grid(p::Int) end end + +# 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")) + 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 + + """ History @@ -97,9 +122,11 @@ end # 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") + gridf = Symbol(lcname, :_grid) + fdmf = Symbol(lcname, :_fdm) + precomputes = Symbol(lcname, :_precomputed) + num_precomputed = 8 @eval begin # Compatibility layer over the "old" API function $fdmf(p::Integer, q::Integer; adapt=1, kwargs...) @@ -107,10 +134,22 @@ for D in (:Forward, :Backward, :Central) return $D(p, q; adapt=adapt, kwargs...) end + # precompute a bunch of grids and coefs for the q=1 case + const $precomputes = ntuple($num_precomputed) do p + p == 1 && return nothing # no grid define for p==1 + grid = $gridf(p) + coefs = _coefs(grid, p, 1) + (grid, coefs) + end + function $D(p::Integer, q::Integer; adapt=1, kwargs...) _check_p_q(p, q) - grid = $gridf(p) - coefs = _coefs(grid, p, q) + if p <= $num_precomputed && q == 1 + grid, coefs = $precomputes[p] + else + grid = $gridf(p) + coefs = _coefs(grid, p, q) + end hist = History(; adapt=adapt, kwargs...) return $D{typeof(grid), typeof(coefs)}(Int(p), Int(q), grid, coefs, hist) end @@ -143,29 +182,6 @@ function Nonstandard(grid::AbstractVector{<:Real}, q::Integer; adapt=0, 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 >= 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")) - 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::FiniteDifferenceMethod, f, x[, Val(false)]; kwargs...) -> Real fdm(m::FiniteDifferenceMethod, f, x, Val(true); kwargs...) -> Tuple{FiniteDifferenceMethod, Real} From 041cb3ce7e84fae8f9789d4d6ce4038cb486f802 Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Sun, 2 Feb 2020 23:08:30 +0000 Subject: [PATCH 2/3] sort out adapt --- src/methods.jl | 52 +++++++++++++++++++++++++++++++------------------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/src/methods.jl b/src/methods.jl index 2f6ca727..8d9a099c 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -138,14 +138,15 @@ for D in (:Forward, :Backward, :Central) const $precomputes = ntuple($num_precomputed) do p p == 1 && return nothing # no grid define for p==1 grid = $gridf(p) - coefs = _coefs(grid, p, 1) - (grid, coefs) + all_coefs = ntuple(q->_coefs(grid, p, q), p - 1) + (grid, all_coefs) end function $D(p::Integer, q::Integer; adapt=1, kwargs...) _check_p_q(p, q) - if p <= $num_precomputed && q == 1 - grid, coefs = $precomputes[p] + if p <= $num_precomputed + grid, all_coefs = $precomputes[p] + coefs = all_coefs[q] else grid = $gridf(p) coefs = _coefs(grid, p, q) @@ -197,6 +198,7 @@ The recognized keywords are: * `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). +* `track_history`: wether to update the history of the method `m` with e.g. accuracy stats. !!! warning Bounds can't be adaptively computed over nonstandard grids; passing a value for @@ -227,14 +229,15 @@ julia> fdm(central_fdm(2, 1), exp, 0, Val(true)) function fdm( m::M, f, - x, + x::T, ::Val{true}; condition=DEFAULT_CONDITION, - bound=_estimate_bound(f(x), condition), - eps=(Base.eps(float(bound)) + TINY), + bound::T=_estimate_bound(f(x)::T, condition)::T, + eps::T=(Base.eps(float(bound)) + TINY)::T, adapt=m.history.adapt, max_step=0.1, -) where M<:FiniteDifferenceMethod + track_history=true # will be set to false if `Val{false}()` used +)::Tuple{M, T} where {M<:FiniteDifferenceMethod, T} if M <: Nonstandard && adapt > 0 throw(ArgumentError("can't adaptively compute bounds over Nonstandard grids")) end @@ -253,7 +256,7 @@ function fdm( # Adaptively compute the bound on the function and derivative values, if applicable. if adapt > 0 - newm = (M.name.wrapper)(p + 1, p) + newm = adapted_FDM(m) dfdx = fdm( newm, f, @@ -272,25 +275,34 @@ function fdm( C₂ = bound * sum(n->abs(coefs[n] * grid[n]^p), eachindex(coefs)) / factorial(p) ĥ = min((q / (p - q) * C₁ / C₂)^(1 / p), max_step) - # 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 - - m.history.eps = eps - m.history.bound = bound - m.history.step = ĥ - m.history.accuracy = accuracy - + dfdx_s = sum(eachindex(grid)) do i + (@inbounds coefs[i] * f(x + ĥ * grid[i]))::T + end + dfdx = dfdx_s / ĥ^q + + if track_history + # Estimate the accuracy of the method. + accuracy = ĥ^(-q) * C₁ + ĥ^(p - q) * C₂ + m.history.eps = eps + m.history.bound = bound + m.history.step = ĥ + m.history.accuracy = accuracy + end return m, dfdx end function fdm(m::FiniteDifferenceMethod, f, x, ::Val{false}=Val(false); kwargs...) - _, dfdx = fdm(m, f, x, Val(true); kwargs...) + _, dfdx = fdm(m, f, x, Val{true}(); track_history=false, kwargs...) return dfdx end +_adapted_FDM(m::M) where M = (M.name.wrapper)(m.p + 1, m.p) +adapted_FDM(m::M) where M = _adapted_FDM(m) +# Deal with type-instability in Central: +adapted_FDM(m::Central{Vector{Int}}) where M = _adapted_FDM(m)::Central{UnitRange{Int}} +adapted_FDM(m::Central{UnitRange{Int}}) where M = _adapted_FDM(m)::Central{Vector{Int}} + ## Deprecations From 0983c77efac9d5a04721fb1c698d8d427eebfc04 Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Sat, 28 Mar 2020 16:25:18 +0000 Subject: [PATCH 3/3] Update src/methods.jl --- src/methods.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods.jl b/src/methods.jl index 8d9a099c..511ec95e 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -134,7 +134,7 @@ for D in (:Forward, :Backward, :Central) return $D(p, q; adapt=adapt, kwargs...) end - # precompute a bunch of grids and coefs for the q=1 case + # precompute a bunch of grids and coefs const $precomputes = ntuple($num_precomputed) do p p == 1 && return nothing # no grid define for p==1 grid = $gridf(p)