Skip to content

Introduce absstep to specify the stepsize to use close to zero #55

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

Merged
merged 1 commit into from
May 26, 2019
Merged

Introduce absstep to specify the stepsize to use close to zero #55

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
34 changes: 22 additions & 12 deletions src/derivatives.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,16 @@
#=
Single-point derivatives of scalar->scalar maps.
=#
function finite_difference_derivative(f, x::T, fdtype::Type{T1}=Val{:central},
returntype::Type{T2}=eltype(x), f_x::Union{Nothing,T}=nothing) where {T<:Number,T1,T2}
function finite_difference_derivative(
f,
x::T,
fdtype::Type{T1}=Val{:central},
returntype::Type{T2}=eltype(x),
f_x::Union{Nothing,T}=nothing;
relstep=default_relstep(fdtype, T),
absstep=relstep) where {T<:Number,T1,T2}
Copy link
Member

@ChrisRackauckas ChrisRackauckas May 26, 2019

Choose a reason for hiding this comment

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

why does this default make sense? Does it recover the previous behavior?

Copy link
Member Author

@andreasnoack andreasnoack May 26, 2019

Choose a reason for hiding this comment

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

Does it recover the previous behavior?

Yes, it's to match the current behavior.


epsilon = compute_epsilon(fdtype, x)
epsilon = compute_epsilon(fdtype, x, relstep, absstep)
if fdtype==Val{:forward}
return (f(x+epsilon) - f(x)) / epsilon
elseif fdtype==Val{:central}
Expand Down Expand Up @@ -69,10 +75,11 @@ function finite_difference_derivative(
returntype :: Type{T2} = eltype(x), # return type of f
fx :: Union{Nothing,AbstractArray{<:Number}} = nothing,
epsilon :: Union{Nothing,AbstractArray{<:Real}} = nothing;
epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2}

df = fill(zero(returntype), size(x))
finite_difference_derivative!(df, f, x, fdtype, returntype, fx, epsilon; epsilon_factor=epsilon_factor)
finite_difference_derivative!(df, f, x, fdtype, returntype, fx, epsilon; relstep=relstep, absstep=absstep)
end

function finite_difference_derivative!(
Expand All @@ -83,22 +90,24 @@ function finite_difference_derivative!(
returntype :: Type{T2} = eltype(x),
fx :: Union{Nothing,AbstractArray{<:Number}} = nothing,
epsilon :: Union{Nothing,AbstractArray{<:Real}} = nothing;
epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2}

cache = DerivativeCache(x, fx, epsilon, fdtype, returntype)
finite_difference_derivative!(df, f, x, cache; epsilon_factor=epsilon_factor)
finite_difference_derivative!(df, f, x, cache; relstep=relstep, absstep=absstep)
end

function finite_difference_derivative!(
df::AbstractArray{<:Number},
f,
x::AbstractArray{<:Number},
cache::DerivativeCache{T1,T2,fdtype,returntype};
epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,fdtype,returntype}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,fdtype,returntype}

fx, epsilon = cache.fx, cache.epsilon
if typeof(epsilon) != Nothing
@. epsilon = compute_epsilon(fdtype, x, epsilon_factor)
@. epsilon = compute_epsilon(fdtype, x, relstep, absstep)
end
if fdtype == Val{:forward}
if typeof(fx) == Nothing
Expand Down Expand Up @@ -127,12 +136,13 @@ function finite_difference_derivative!(
f,
x::StridedArray,
cache::DerivativeCache{T1,T2,fdtype,returntype};
epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,fdtype,returntype}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,fdtype,returntype}

if fdtype == Val{:forward}
fx = cache.fx
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
epsilon = compute_epsilon(Val{:forward}, x[i], relstep, absstep)
x_plus = x[i] + epsilon
if typeof(fx) == Nothing
df[i] = (f(x_plus) - f(x[i])) / epsilon
Expand All @@ -142,7 +152,7 @@ function finite_difference_derivative!(
end
elseif fdtype == Val{:central}
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
epsilon = compute_epsilon(Val{:central}, x[i], relstep, absstep)
epsilon_double_inv = one(typeof(epsilon)) / (2*epsilon)
x_plus, x_minus = x[i]+epsilon, x[i]-epsilon
df[i] = (f(x_plus) - f(x_minus)) * epsilon_double_inv
Expand Down
14 changes: 7 additions & 7 deletions src/finitediff.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,19 +6,19 @@ Very heavily inspired by Calculus.jl, but with an emphasis on performance and Di
Compute the finite difference interval epsilon.
Reference: Numerical Recipes, chapter 5.7.
=#
@inline function compute_epsilon(::Type{Val{:forward}}, x::T, eps_sqrt=sqrt(eps(real(T)))) where T<:Number
eps_sqrt * max(one(real(T)), abs(x))
@inline function compute_epsilon(::Type{Val{:forward}}, x::T, relstep::Real, absstep::Real) where T<:Number
return relstep*abs(x) + absstep
end

@inline function compute_epsilon(::Type{Val{:central}}, x::T, eps_cbrt=cbrt(eps(real(T)))) where T<:Number
eps_cbrt * max(one(real(T)), abs(x))
@inline function compute_epsilon(::Type{Val{:central}}, x::T, relstep::Real, absstep::Real) where T<:Number
return relstep*abs(x) + absstep
end

@inline function compute_epsilon(::Type{Val{:complex}}, x::T, ::Union{Nothing,T}=nothing) where T<:Real
eps(T)
@inline function compute_epsilon(::Type{Val{:complex}}, x::T, ::Union{Nothing,T}=nothing, ::Union{Nothing,T}=nothing) where T<:Real
return eps(T)
end

@inline function compute_epsilon_factor(fdtype::DataType, ::Type{T}) where T<:Number
@inline function default_relstep(fdtype::DataType, ::Type{T}) where T<:Number
if fdtype==Val{:forward}
return sqrt(eps(real(T)))
elseif fdtype==Val{:central}
Expand Down
36 changes: 20 additions & 16 deletions src/gradients.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,8 @@ function finite_difference_gradient(
fx::Union{Nothing,AbstractArray{<:Number}}=nothing,
c1::Union{Nothing,AbstractArray{<:Number}}=nothing,
c2::Union{Nothing,AbstractArray{<:Number}}=nothing;
epsilon_factor=compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,T3}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3}

if typeof(x) <: AbstractArray
df = fill(zero(returntype), size(x))
Expand All @@ -133,7 +134,7 @@ function finite_difference_gradient(
end
end
cache = GradientCache(df, x, fdtype, returntype, inplace)
finite_difference_gradient!(df, f, x, cache, epsilon_factor=epsilon_factor)
finite_difference_gradient!(df, f, x, cache, relstep=relstep, absstep=absstep)
end

function finite_difference_gradient!(
Expand All @@ -146,24 +147,26 @@ function finite_difference_gradient!(
fx::Union{Nothing,AbstractArray{<:Number}}=nothing,
c1::Union{Nothing,AbstractArray{<:Number}}=nothing,
c2::Union{Nothing,AbstractArray{<:Number}}=nothing;
epsilon_factor=compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,T3}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3}

cache = GradientCache(df, x, fdtype, returntype, inplace)
finite_difference_gradient!(df, f, x, cache, epsilon_factor=epsilon_factor)
finite_difference_gradient!(df, f, x, cache, relstep=relstep, absstep=absstep)
end

function finite_difference_gradient(
f,
x,
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
epsilon_factor=compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,T3,fdtype,returntype,inplace}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}

if typeof(x) <: AbstractArray
df = fill(zero(returntype), size(x))
else
df = zero(cache.c1)
end
finite_difference_gradient!(df, f, x, cache, epsilon_factor=epsilon_factor)
finite_difference_gradient!(df, f, x, cache, relstep=relstep, absstep=absstep)
df
end

Expand All @@ -174,13 +177,14 @@ function finite_difference_gradient!(
f,
x::AbstractArray{<:Number},
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
epsilon_factor=compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,T3,fdtype,returntype,inplace}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}

# NOTE: in this case epsilon is a vector, we need two arrays for epsilon and x1
# c1 denotes x1, c2 is epsilon
fx, c1, c2 = cache.fx, cache.c1, cache.c2
if fdtype != Val{:complex}
@. c2 = compute_epsilon(fdtype, x, epsilon_factor)
@. c2 = compute_epsilon(fdtype, x, relstep, absstep)
copyto!(c1,x)
end
if fdtype == Val{:forward}
Expand Down Expand Up @@ -246,8 +250,8 @@ function finite_difference_gradient!(
f,
x::StridedVector{<:Number},
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
epsilon_factor=compute_epsilon_factor(fdtype, eltype(x)),
) where {T1,T2,T3,fdtype,returntype,inplace}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}

# c1 is x1 if we need a complex copy of x, otherwise Nothing
# c2 is Nothing
Expand All @@ -259,7 +263,7 @@ function finite_difference_gradient!(
end
if fdtype == Val{:forward}
for i ∈ eachindex(x)
epsilon = compute_epsilon(fdtype, x[i], epsilon_factor)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep)
x_old = x[i]
if typeof(fx) != Nothing
x[i] += epsilon
Expand Down Expand Up @@ -296,7 +300,7 @@ function finite_difference_gradient!(
end
elseif fdtype == Val{:central}
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(fdtype, x[i], epsilon_factor)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep)
x_old = x[i]
x[i] += epsilon
dfi = f(x)
Expand Down Expand Up @@ -344,15 +348,15 @@ function finite_difference_gradient!(
f,
x::Number,
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
epsilon_factor=compute_epsilon_factor(fdtype, eltype(x))
) where {T1,T2,T3,fdtype,returntype,inplace}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}

# NOTE: in this case epsilon is a scalar, we need two arrays for fx1 and fx2
# c1 denotes fx1, c2 is fx2, sizes guaranteed by the cache constructor
fx, c1, c2 = cache.fx, cache.c1, cache.c2

if fdtype == Val{:forward}
epsilon = compute_epsilon(Val{:forward}, x, epsilon_factor)
epsilon = compute_epsilon(Val{:forward}, x, relstep, absstep)
if inplace == Val{true}
f(c1, x+epsilon)
else
Expand All @@ -369,7 +373,7 @@ function finite_difference_gradient!(
@. df = (c1 - c2) / epsilon
end
elseif fdtype == Val{:central}
epsilon = compute_epsilon(Val{:central}, x, epsilon_factor)
epsilon = compute_epsilon(Val{:central}, x, relstep, absstep)
if inplace == Val{true}
f(c1, x+epsilon)
f(c2, x-epsilon)
Expand Down
17 changes: 10 additions & 7 deletions src/jacobians.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,21 +90,23 @@ function finite_difference_jacobian(f, x::AbstractArray{<:Number},
returntype :: Type{T2}=eltype(x),
inplace :: Type{Val{T3}}=Val{true},
f_in :: Union{T2,Nothing}=nothing;
epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,T3}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3}

cache = JacobianCache(x, fdtype, returntype, inplace)
finite_difference_jacobian(f, x, cache, f_in; epsilon_factor=epsilon_factor)
finite_difference_jacobian(f, x, cache, f_in; relstep=relstep, absstep=absstep)
end

function finite_difference_jacobian(
f,
x,
cache::JacobianCache{T1,T2,T3,fdtype,returntype,inplace},
f_in=nothing;
epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,T3,fdtype,returntype,inplace}
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}

J = fill(zero(eltype(x)), length(x), length(x))
finite_difference_jacobian!(J, f, x, cache, f_in; epsilon_factor=epsilon_factor)
finite_difference_jacobian!(J, f, x, cache, f_in; relstep=relstep, absstep=absstep)
J
end

Expand All @@ -114,7 +116,8 @@ function finite_difference_jacobian!(
x::AbstractArray{<:Number},
cache::JacobianCache{T1,T2,T3,fdtype,returntype,inplace},
f_in::Union{T2,Nothing}=nothing;
epsilon_factor = compute_epsilon_factor(fdtype, eltype(x))) where {T1,T2,T3,fdtype,returntype,inplace}
relstep = default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}

m, n = size(J)
x1, fx, fx1 = cache.x1, cache.fx, cache.fx1
Expand All @@ -123,7 +126,7 @@ function finite_difference_jacobian!(
if fdtype == Val{:forward}
vfx1 = vec(fx1)
@inbounds for i ∈ 1:n
epsilon = compute_epsilon(Val{:forward}, x[i], epsilon_factor)
epsilon = compute_epsilon(Val{:forward}, x[i], relstep, absstep)
x1_save = x1[i]
x1[i] += epsilon
if inplace == Val{true}
Expand All @@ -148,7 +151,7 @@ function finite_difference_jacobian!(
elseif fdtype == Val{:central}
vfx1 = vec(fx1)
@inbounds for i ∈ 1:n
epsilon = compute_epsilon(Val{:central}, x[i], epsilon_factor)
epsilon = compute_epsilon(Val{:central}, x[i], relstep, absstep)
x1_save = x1[i]
x_save = x[i]
x1[i] += epsilon
Expand Down
Loading