Skip to content

Precompute grids and coeffs #61

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 3 commits into from
Closed

Precompute grids and coeffs #61

Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7d0f13f
Precompute grids and coeffs
oxinabox Feb 2, 2020
041cb3c
sort out adapt
oxinabox Feb 2, 2020
0983c77
Update src/methods.jl
oxinabox Mar 28, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 71 additions & 43 deletions src/methods.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -97,20 +122,35 @@ 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...)
_dep_kwarg(kwargs)
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
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be defined in a constant?
Or should it be declared outside the eval, and then inlined in as a local literal
and thus should constant fold

p == 1 && return nothing # no grid define for p==1
grid = $gridf(p)
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)
grid = $gridf(p)
coefs = _coefs(grid, p, q)
if p <= $num_precomputed
grid, all_coefs = $precomputes[p]
coefs = all_coefs[q]
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
Expand Down Expand Up @@ -143,29 +183,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}
Expand All @@ -181,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
Expand Down Expand Up @@ -211,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
Expand All @@ -237,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,
Expand All @@ -256,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

Expand Down