From f6bd88423f095c4cd5d284cb48338253e5b61b93 Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Sat, 1 Feb 2020 21:55:26 +0000 Subject: [PATCH 1/6] Do gradient via mutating and unmutating cell --- src/grad.jl | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/grad.jl b/src/grad.jl index 815a10a5..447cbd55 100644 --- a/src/grad.jl +++ b/src/grad.jl @@ -8,19 +8,25 @@ Approximate the gradient of `f` at `xs...` using `fdm`. Assumes that `f(xs...)` """ function grad end -function grad(fdm, f, x::AbstractArray{T}) where T <: Number +function _grad(fdm, f, x::AbstractArray{T}) where T <: Number + # x must be mutable, we will mutate it and then mutate it back. dx = similar(x) - tmp = similar(x) for k in eachindex(x) dx[k] = fdm(zero(T)) do ϵ - tmp .= x - tmp[k] += ϵ - return f(tmp) + xk = x[k] + x[k] = xk + ϵ + ret = f(x) + x[k] = xk # Can't do `x[k] -= ϵ` as floating-point math is not associative + return ret end end return (dx, ) end +grad(fdm, f, x::Array{<:Number}) = _grad(fdm, f, x) +# Fallback for when we don't know `x` will be mutable: +grad(fdm, f, x::AbstractArray{<:Number}) = _grad(fdm, f, similar(x).=x) + grad(fdm, f, x::Real) = (fdm(f, x), ) grad(fdm, f, x::Tuple) = (grad(fdm, (xs...)->f(xs), x...), ) From 37f63af4d978302d0120eebb798e0baebaea8020 Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Sat, 1 Feb 2020 22:19:44 +0000 Subject: [PATCH 2/6] Remove copy from Dict grad also --- src/grad.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/grad.jl b/src/grad.jl index 447cbd55..d3443cb9 100644 --- a/src/grad.jl +++ b/src/grad.jl @@ -31,14 +31,15 @@ grad(fdm, f, x::Real) = (fdm(f, x), ) grad(fdm, f, x::Tuple) = (grad(fdm, (xs...)->f(xs), x...), ) function grad(fdm, f, d::Dict{K, V}) where {K, V} - dd = Dict{K, V}() + ∇d = Dict{K, V}() for (k, v) in d + dk = d[k] function f′(x) - tmp = copy(d) - tmp[k] = x - return f(tmp) + d[k] = x + return f(d) end - dd[k] = grad(fdm, f′, v)[1] + ∇d[k] = grad(fdm, f′, v)[1] + d[k] = dk end return (dd, ) end From f677183c5e1f8723952f28dff1fe9d997c039a12 Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Sat, 1 Feb 2020 23:23:01 +0000 Subject: [PATCH 3/6] correct dict return name --- src/grad.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grad.jl b/src/grad.jl index d3443cb9..75242b95 100644 --- a/src/grad.jl +++ b/src/grad.jl @@ -41,7 +41,7 @@ function grad(fdm, f, d::Dict{K, V}) where {K, V} ∇d[k] = grad(fdm, f′, v)[1] d[k] = dk end - return (dd, ) + return (∇d, ) end function grad(fdm, f, x) From 6b6b2c991e54e10c88461f54522dd872632e40ce Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Sun, 2 Feb 2020 22:14:35 +0000 Subject: [PATCH 4/6] make tracking history not happen and specify some types. --- src/grad.jl | 2 +- src/methods.jl | 35 ++++++++++++++++++++--------------- 2 files changed, 21 insertions(+), 16 deletions(-) diff --git a/src/grad.jl b/src/grad.jl index 75242b95..aef1811c 100644 --- a/src/grad.jl +++ b/src/grad.jl @@ -17,7 +17,7 @@ function _grad(fdm, f, x::AbstractArray{T}) where T <: Number x[k] = xk + ϵ ret = f(x) x[k] = xk # Can't do `x[k] -= ϵ` as floating-point math is not associative - return ret + return ret::T end end return (dx, ) diff --git a/src/methods.jl b/src/methods.jl index 20b56ea9..5a3b8c86 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -181,6 +181,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 @@ -211,14 +212,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 @@ -256,22 +258,25 @@ 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 From ff18572248bd6f441f71c510625c96ba7095b25e Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Fri, 7 Feb 2020 14:38:07 +0000 Subject: [PATCH 5/6] Update src/methods.jl Co-Authored-By: willtebbutt --- src/methods.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods.jl b/src/methods.jl index 5a3b8c86..45a99122 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -181,7 +181,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. +* `track_history`: whether 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 From bb43e4262e0c64f99d577839c741938f55a12c06 Mon Sep 17 00:00:00 2001 From: Lyndon White Date: Fri, 7 Feb 2020 14:51:51 +0000 Subject: [PATCH 6/6] only check type on eps --- src/methods.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/methods.jl b/src/methods.jl index 45a99122..e7911e32 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -212,11 +212,11 @@ julia> fdm(central_fdm(2, 1), exp, 0, Val(true)) function fdm( m::M, f, - x::T, + x, ::Val{true}; condition=DEFAULT_CONDITION, - bound::T=_estimate_bound(f(x)::T, condition)::T, - eps::T=(Base.eps(float(bound)) + TINY)::T, + bound=_estimate_bound(f(x), condition), + eps::T=(Base.eps(float(bound)) + TINY), adapt=m.history.adapt, max_step=0.1, track_history=true # will be set to false if `Val{false}()` used