Skip to content

Commit

Permalink
Merge pull request #64 from JuliaDiffEq/myb/dir
Browse files Browse the repository at this point in the history 
Add `dir` keyword argument
  • Loading branch information
@ChrisRackauckas
ChrisRackauckas committed Jul 1, 2019
2 parents e6fbdde + 4298303 commit c4be5ff
Show file tree
Hide file tree
Showing 5 changed files with 59 additions and 35 deletions.
17 changes: 10 additions & 7 deletions src/derivatives.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,10 @@ function finite_difference_derivative(
returntype::Type{T2}=eltype(x),
f_x::Union{Nothing,T}=nothing;
relstep=default_relstep(fdtype, T),
absstep=relstep) where {T<:Number,T1,T2}
absstep=relstep,
dir=true) where {T<:Number,T1,T2}

epsilon = compute_epsilon(fdtype, x, relstep, absstep)
epsilon = compute_epsilon(fdtype, x, relstep, absstep, dir)
if fdtype==Val{:forward}
return (f(x+epsilon) - f(x)) / epsilon
elseif fdtype==Val{:central}
Expand Down Expand Up @@ -103,11 +104,12 @@ function finite_difference_derivative!(
x::AbstractArray{<:Number},
cache::DerivativeCache{T1,T2,fdtype,returntype};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,fdtype,returntype}
absstep=relstep,
dir=true) where {T1,T2,fdtype,returntype}

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

if fdtype == Val{:forward}
fx = cache.fx
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(Val{:forward}, x[i], relstep, absstep)
epsilon = compute_epsilon(Val{:forward}, x[i], relstep, absstep, dir)
x_plus = x[i] + epsilon
if typeof(fx) == Nothing
df[i] = (f(x_plus) - f(x[i])) / epsilon
Expand All @@ -152,7 +155,7 @@ function finite_difference_derivative!(
end
elseif fdtype == Val{:central}
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(Val{:central}, x[i], relstep, absstep)
epsilon = compute_epsilon(Val{:central}, x[i], relstep, absstep, dir)
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
10 changes: 5 additions & 5 deletions src/finitediff.jl
View file
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, relstep::Real, absstep::Real) where T<:Number
return max(relstep*abs(x), absstep)
@inline function compute_epsilon(::Type{Val{:forward}}, x::T, relstep::Real, absstep::Real, dir::Real) where T<:Number
return max(relstep*abs(x), absstep)*dir
end

@inline function compute_epsilon(::Type{Val{:central}}, x::T, relstep::Real, absstep::Real) where T<:Number
@inline function compute_epsilon(::Type{Val{:central}}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
return max(relstep*abs(x), absstep)
end

@inline function compute_epsilon(::Type{Val{:hcentral}}, x::T, relstep::Real, absstep::Real) where T<:Number
@inline function compute_epsilon(::Type{Val{:hcentral}}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
return max(relstep*abs(x), absstep)
end

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

Expand Down
31 changes: 18 additions & 13 deletions src/gradients.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,8 @@ function finite_difference_gradient(
c1::Union{Nothing,AbstractArray{<:Number}}=nothing,
c2::Union{Nothing,AbstractArray{<:Number}}=nothing;
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3}
absstep=relstep,
dir=true) where {T1,T2,T3}

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

function finite_difference_gradient!(
Expand All @@ -159,14 +160,15 @@ function finite_difference_gradient(
x,
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}
absstep=relstep,
dir=true) 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, relstep=relstep, absstep=absstep)
finite_difference_gradient!(df, f, x, cache, relstep=relstep, absstep=absstep, dir=dir)
df
end

Expand All @@ -178,18 +180,19 @@ function finite_difference_gradient!(
x::AbstractArray{<:Number},
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}
absstep=relstep,
dir=true) 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, relstep, absstep)
@. c2 = compute_epsilon(fdtype, x, relstep, absstep, dir)
copyto!(c1,x)
end
if fdtype == Val{:forward}
@inbounds for i ∈ eachindex(x)
epsilon = c2[i]
epsilon = c2[i]*dir
c1_old = c1[i]
c1[i] += epsilon
if typeof(fx) != Nothing
Expand Down Expand Up @@ -251,7 +254,8 @@ function finite_difference_gradient!(
x::StridedVector{<:Number},
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}
absstep=relstep,
dir=true) 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 @@ -263,7 +267,7 @@ function finite_difference_gradient!(
end
if fdtype == Val{:forward}
for i ∈ eachindex(x)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep, dir)
x_old = x[i]
if typeof(fx) != Nothing
x[i] += epsilon
Expand Down Expand Up @@ -300,7 +304,7 @@ function finite_difference_gradient!(
end
elseif fdtype == Val{:central}
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep, dir)
x_old = x[i]
x[i] += epsilon
dfi = f(x)
Expand Down Expand Up @@ -349,14 +353,15 @@ function finite_difference_gradient!(
x::Number,
cache::GradientCache{T1,T2,T3,fdtype,returntype,inplace};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep) where {T1,T2,T3,fdtype,returntype,inplace}
absstep=relstep,
dir=true) 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, relstep, absstep)
epsilon = compute_epsilon(Val{:forward}, x, relstep, absstep, dir)
if inplace == Val{true}
f(c1, x+epsilon)
else
Expand All @@ -373,7 +378,7 @@ function finite_difference_gradient!(
@. df = (c1 - c2) / epsilon
end
elseif fdtype == Val{:central}
epsilon = compute_epsilon(Val{:central}, x, relstep, absstep)
epsilon = compute_epsilon(Val{:central}, x, relstep, absstep, dir)
if inplace == Val{true}
f(c1, x+epsilon)
f(c2, x-epsilon)
Expand Down
21 changes: 12 additions & 9 deletions src/jacobians.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,10 +132,11 @@ function finite_difference_jacobian(f, x::AbstractArray{<:Number},
f_in :: Union{T2,Nothing}=nothing;
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep,
color = eachindex(x)) where {T1,T2,T3}
color = eachindex(x),
dir=true) where {T1,T2,T3}

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

function finite_difference_jacobian(
Expand All @@ -145,10 +146,11 @@ function finite_difference_jacobian(
f_in=nothing;
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep,
color = cache.color) where {T1,T2,T3,cType,fdtype,returntype,inplace}
color = cache.color,
dir=true) where {T1,T2,T3,cType,fdtype,returntype,inplace}
_J = false .* x .* x'
_J isa SMatrix ? J = MArray(_J) : J = _J
finite_difference_jacobian!(J, f, x, cache, f_in; relstep=relstep, absstep=absstep, color=color)
finite_difference_jacobian!(J, f, x, cache, f_in; relstep=relstep, absstep=absstep, color=color, dir=dir)
_J isa SMatrix ? SArray(J) : J
end

Expand All @@ -160,7 +162,8 @@ function finite_difference_jacobian!(
f_in::Union{T2,Nothing}=nothing;
relstep = default_relstep(fdtype, eltype(x)),
absstep=relstep,
color = cache.color) where {T1,T2,T3,cType,fdtype,returntype,inplace}
color = cache.color,
dir = true) where {T1,T2,T3,cType,fdtype,returntype,inplace}

m, n = size(J)
x1, fx, fx1 = cache.x1, cache.fx, cache.fx1
Expand Down Expand Up @@ -190,7 +193,7 @@ function finite_difference_jacobian!(
@inbounds for color_i ∈ 1:maximum(color)

if color isa Base.OneTo || color isa StaticArrays.SOneTo # Dense matrix
epsilon = compute_epsilon(Val{:forward}, x1[color_i], relstep, absstep)
epsilon = compute_epsilon(Val{:forward}, x1[color_i], relstep, absstep, dir)
x1_save = x1[color_i]
if inplace == Val{true}
x1[color_i] += epsilon
Expand All @@ -204,7 +207,7 @@ function finite_difference_jacobian!(
tmp += abs2(x1[i])
end
end
epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep)
epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir)

if inplace != Val{true}
_x1 = copy(x1)
Expand Down Expand Up @@ -275,7 +278,7 @@ function finite_difference_jacobian!(
@inbounds for color_i ∈ 1:maximum(color)

if color isa Base.OneTo || color isa StaticArrays.SOneTo # Dense matrix
epsilon = compute_epsilon(Val{:central}, x[color_i], relstep, absstep)
epsilon = compute_epsilon(Val{:central}, x[color_i], relstep, absstep, dir)
x1_save = x1[color_i]
x_save = x[color_i]
if inplace == Val{true}
Expand All @@ -292,7 +295,7 @@ function finite_difference_jacobian!(
tmp += abs2(x1[i])
end
end
epsilon = compute_epsilon(Val{:central}, sqrt(tmp), relstep, absstep)
epsilon = compute_epsilon(Val{:central}, sqrt(tmp), relstep, absstep, dir)

if inplace != Val{true}
_x1 = copy(x1)
Expand Down
15 changes: 14 additions & 1 deletion test/finitedifftests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@ err_func(a,b) = maximum(abs.(a-b))

@time @testset "Derivative single point f : R -> R tests" begin
@test err_func(DiffEqDiffTools.finite_difference_derivative(sin, π/4, Val{:forward}), √2/2) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_derivative(x->x > π/4 ? error() : sin(x), π/4, Val{:forward}, dir=-1), √2/2) < 1e-4
@test_throws Any err_func(DiffEqDiffTools.finite_difference_derivative(x->x >= π/4 ? error() : sin(x), π/4, Val{:forward}), √2/2) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_derivative(sin, π/4, Val{:central}), √2/2) < 1e-8
@test err_func(DiffEqDiffTools.finite_difference_derivative(sin, π/4, Val{:complex}), √2/2) < 1e-15
end
Expand Down Expand Up @@ -167,6 +169,8 @@ end
# Gradient tests
f(x) = 2x[1] + x[2]^2
x = rand(2)
z = copy(x)
ff(k) = !all(k .<= z) ? error() : f(k)
fx = f(x)
df = fill(0.,2)
df_ref = [2., 2*x[2]]
Expand All @@ -176,6 +180,8 @@ complex_cache = DiffEqDiffTools.GradientCache(df,x,Val{:complex})

@time @testset "Gradient of f:vector->scalar real-valued tests" begin
@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:forward}), df_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_gradient(ff, x, Val{:forward}, dir=-1), df_ref) < 1e-4
@test_throws Any err_func(DiffEqDiffTools.finite_difference_gradient(ff, x, Val{:forward}), df_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:forward}, relstep=sqrt(eps())), df_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:central}), df_ref) < 1e-8
@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:complex}), df_ref) < 1e-15
Expand Down Expand Up @@ -252,9 +258,10 @@ central_cache = DiffEqDiffTools.GradientCache(df,x,Val{:central},eltype(df))
end

f(df,x) = (df[1]=sin(x); df[2]=cos(x); df)
x = 2π * rand()
z = x = 2π * rand()
fx = fill(0.,2)
f(fx,x)
ff(df, x) = !all(x .<= z) ? error() : f(df, x)
df = fill(0.,2)
df_ref = [cos(x), -sin(x)]
forward_cache = DiffEqDiffTools.GradientCache(df,x,Val{:forward})
Expand All @@ -265,6 +272,8 @@ complex_cache = DiffEqDiffTools.GradientCache(df,x,Val{:complex})
@time @testset "Gradient of f:scalar->vector real-valued tests" begin
@test_broken err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:forward}), df_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:forward}, eltype(x), Val{true}, fx), df_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_gradient(ff, x, Val{:forward}, eltype(x), Val{true}, fx, dir=-1), df_ref) < 1e-4
@test_throws Any err_func(DiffEqDiffTools.finite_difference_gradient(ff, x, Val{:forward}), df_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:central}, eltype(x), Val{true}, fx), df_ref) < 1e-8
@test err_func(DiffEqDiffTools.finite_difference_gradient(f, x, Val{:complex}, eltype(x), Val{true}, fx), df_ref) < 1e-15

Expand Down Expand Up @@ -303,6 +312,8 @@ function f(fvec,x)
fvec[2] = sin(x[2]*exp(x[1])-1)
end
x = rand(2); y = rand(2)
z = copy(x)
ff(df, x) = !all(x .<= z) ? error() : f(df, x)
f(y,x)
J_ref = [[-7+x[2]^3 3*(3+x[1])*x[2]^2]; [exp(x[1])*x[2]*cos(1-exp(x[1])*x[2]) exp(x[1])*cos(1-exp(x[1])*x[2])]]
J = zero(J_ref)
Expand All @@ -316,6 +327,8 @@ f_in = copy(y)

@time @testset "Jacobian StridedArray real-valued tests" begin
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, forward_cache), J_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_jacobian(ff, x, forward_cache, dir=-1), J_ref) < 1e-4
@test_throws Any err_func(DiffEqDiffTools.finite_difference_jacobian(ff, x, forward_cache), J_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, forward_cache, relstep=sqrt(eps())), J_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, forward_cache, f_in), J_ref) < 1e-4
@test err_func(DiffEqDiffTools.finite_difference_jacobian(f, x, central_cache), J_ref) < 1e-8
Expand Down

0 comments on commit c4be5ff

Please sign in to comment.