Skip to content

Make FiniteDiff take both instances and types #116

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 3 commits into from
Sep 16, 2020
Merged

Make FiniteDiff take both instances and types #116

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
39da4b5
Use Val(fdtype)
YingboMa Sep 16, 2020
de47945
Inplace: Val{true} -> Val(true)
YingboMa Sep 16, 2020
eae102d
New minor release
YingboMa Sep 16, 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
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "FiniteDiff"
uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
version = "2.6.0"
version = "2.7.0"

[deps]
ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
Expand Down
42 changes: 22 additions & 20 deletions src/derivatives.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,19 +4,20 @@ Single-point derivatives of scalar->scalar maps.
function finite_difference_derivative(
f,
x::T,
fdtype=Val{:central},
fdtype=Val(:central),
returntype=eltype(x),
f_x=nothing;
relstep=default_relstep(fdtype, T),
absstep=relstep,
dir=true) where {T<:Number}

fdtype isa Type && (fdtype = fdtype())
epsilon = compute_epsilon(fdtype, x, relstep, absstep, dir)
if fdtype==Val{:forward}
if fdtype==Val(:forward)
return (f(x+epsilon) - f(x)) / epsilon
elseif fdtype==Val{:central}
elseif fdtype==Val(:central)
return (f(x+epsilon) - f(x-epsilon)) / (2*epsilon)
elseif fdtype==Val{:complex} && returntype<:Real
elseif fdtype==Val(:complex) && returntype<:Real
return imag(f(x+im*epsilon)) / epsilon
end
fdtype_error(returntype)
Expand All @@ -36,15 +37,16 @@ function DerivativeCache(
x :: AbstractArray{<:Number},
fx :: Union{Nothing,AbstractArray{<:Number}} = nothing,
epsilon :: Union{Nothing,AbstractArray{<:Real}} = nothing,
fdtype :: Type{T1} = Val{:central},
fdtype :: Type{T1} = Val(:central),
returntype :: Type{T2} = eltype(x)) where {T1,T2}

if fdtype==Val{:complex} && !(eltype(returntype)<:Real)
fdtype isa Type && (fdtype = fdtype())
if fdtype==Val(:complex) && !(eltype(returntype)<:Real)
fdtype_error(returntype)
end

if fdtype!=Val{:forward} && typeof(fx)!=Nothing
@warn("Pre-computed function values are only useful for fdtype==Val{:forward}.")
if fdtype!=Val(:forward) && typeof(fx)!=Nothing
@warn("Pre-computed function values are only useful for fdtype==Val(:forward).")
_fx = nothing
else
# more runtime sanity checks?
Expand All @@ -54,8 +56,8 @@ function DerivativeCache(
if typeof(epsilon)!=Nothing && typeof(x)<:StridedArray && typeof(fx)<:Union{Nothing,StridedArray} && 1==2
@warn("StridedArrays don't benefit from pre-allocating epsilon.")
_epsilon = nothing
elseif typeof(epsilon)!=Nothing && fdtype==Val{:complex}
@warn("Val{:complex} makes the epsilon array redundant.")
elseif typeof(epsilon)!=Nothing && fdtype==Val(:complex)
@warn("Val(:complex) makes the epsilon array redundant.")
_epsilon = nothing
else
if typeof(epsilon)==Nothing || eltype(epsilon)!=real(eltype(x))
Expand All @@ -72,7 +74,7 @@ Compute the derivative df of a scalar-valued map f at a collection of points x.
function finite_difference_derivative(
f,
x,
fdtype = Val{:central},
fdtype = Val(:central),
returntype = eltype(x), # return type of f
fx = nothing,
epsilon = nothing;
Expand All @@ -87,7 +89,7 @@ function finite_difference_derivative!(
df,
f,
x,
fdtype = Val{:central},
fdtype = Val(:central),
returntype = eltype(x),
fx = nothing,
epsilon = nothing;
Expand All @@ -111,15 +113,15 @@ function finite_difference_derivative!(
if typeof(epsilon) != Nothing
@. epsilon = compute_epsilon(fdtype, x, relstep, absstep, dir)
end
if fdtype == Val{:forward}
if fdtype == Val(:forward)
if typeof(fx) == Nothing
@. df = (f(x+epsilon) - f(x)) / epsilon
else
@. df = (f(x+epsilon) - fx) / epsilon
end
elseif fdtype == Val{:central}
elseif fdtype == Val(:central)
@. df = (f(x+epsilon) - f(x-epsilon)) / (2 * epsilon)
elseif fdtype == Val{:complex} && returntype<:Real
elseif fdtype == Val(:complex) && returntype<:Real
epsilon_complex = eps(eltype(x))
@. df = imag(f(x+im*epsilon_complex)) / epsilon_complex
else
Expand All @@ -142,25 +144,25 @@ function finite_difference_derivative!(
absstep=relstep,
dir=true) where {T1,T2,fdtype,returntype}

if fdtype == Val{:forward}
if fdtype == Val(:forward)
fx = cache.fx
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(Val{:forward}, x[i], relstep, absstep, dir)
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
else
df[i] = (f(x_plus) - fx[i]) / epsilon
end
end
elseif fdtype == Val{:central}
elseif fdtype == Val(:central)
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(Val{:central}, x[i], relstep, absstep, dir)
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
end
elseif fdtype == Val{:complex}
elseif fdtype == Val(:complex)
epsilon_complex = eps(eltype(x))
@inbounds for i ∈ eachindex(x)
df[i] = imag(f(x[i]+im*epsilon_complex)) / epsilon_complex
Expand Down
23 changes: 12 additions & 11 deletions src/epsilons.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,38 +6,39 @@ 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, dir::Real) where T<:Number
@inline function compute_epsilon(::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, dir=nothing) where T<:Number
@inline function compute_epsilon(::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, dir=nothing) where T<:Number
@inline function compute_epsilon(::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, dir=nothing) where T<:Real
@inline function compute_epsilon(::Val{:complex}, x::T, ::Union{Nothing,T}=nothing, ::Union{Nothing,T}=nothing, dir=nothing) where T<:Real
return eps(T)
end

@inline function default_relstep(fdtype::DataType, ::Type{T}) where T<:Number
if fdtype==Val{:forward}
default_relstep(v::Type, T) = default_relstep(v(), T)
@inline function default_relstep(::Val{fdtype}, ::Type{T}) where {fdtype,T<:Number}
if fdtype==:forward
return sqrt(eps(real(T)))
elseif fdtype==Val{:central}
elseif fdtype==:central
return cbrt(eps(real(T)))
elseif fdtype==Val{:hcentral}
elseif fdtype==:hcentral
eps(T)^(1/4)
else
return one(real(T))
end
end

function fdtype_error(funtype::Type{T}=Float64) where T
if funtype<:Real
function fdtype_error(::Type{T}=Float64) where T
if T<:Real
error("Unrecognized fdtype: valid values are Val{:forward}, Val{:central} and Val{:complex}.")
elseif funtype<:Complex
elseif T<:Complex
error("Unrecognized fdtype: valid values are Val{:forward} or Val{:central}.")
else
error("Unrecognized returntype: should be a subtype of Real or Complex.")
Expand Down
57 changes: 30 additions & 27 deletions src/gradients.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,14 @@ end
function GradientCache(
df,
x,
fdtype = Val{:central},
fdtype = Val(:central),
returntype = eltype(df),
inplace = Val{true})
inplace = Val(true))

fdtype isa Type && (fdtype = fdtype())
inplace isa Type && (inplace = inplace())
if typeof(x)<:AbstractArray # the vector->scalar case
if fdtype!=Val{:complex} # complex-mode FD only needs one cache, for x+eps*im
if fdtype!=Val(:complex) # complex-mode FD only needs one cache, for x+eps*im
if typeof(x)<:StridedVector
if eltype(df)<:Complex && !(eltype(x)<:Complex)
_c1 = zero(Complex{eltype(x)}) .* x
Expand All @@ -37,7 +39,7 @@ function GradientCache(
_c3 = similar(x)
else # the scalar->vector case
# need cache arrays for fx1 and fx2, except in complex mode, which needs one complex array
if fdtype != Val{:complex}
if fdtype != Val(:complex)
_c1 = similar(df)
_c2 = similar(df)
else
Expand All @@ -55,20 +57,21 @@ end
function finite_difference_gradient(
f,
x,
fdtype = Val{:central},
fdtype = Val(:central),
returntype = eltype(x),
inplace = Val{true},
inplace = Val(true),
fx = nothing,
c1 = nothing,
c2 = nothing;
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep,
dir=true)

inplace isa Type && (inplace = inplace())
if typeof(x) <: AbstractArray
df = zero(returntype) .* x
else
if inplace == Val{true}
if inplace == Val(true)
if typeof(fx)==Nothing && typeof(c1)==Nothing && typeof(c2)==Nothing
error("In the scalar->vector in-place map case, at least one of fx, c1 or c2 must be provided, otherwise we cannot infer the return size.")
else
Expand All @@ -89,9 +92,9 @@ function finite_difference_gradient!(
df,
f,
x,
fdtype=Val{:central},
fdtype=Val(:central),
returntype=eltype(df),
inplace=Val{true},
inplace=Val(true),
fx=nothing,
c1=nothing,
c2=nothing;
Expand Down Expand Up @@ -133,12 +136,12 @@ function finite_difference_gradient!(
# 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, c3 = cache.fx, cache.c1, cache.c2, cache.c3
if fdtype != Val{:complex} && ArrayInterface.fast_scalar_indexing(c2)
if fdtype != Val(:complex) && ArrayInterface.fast_scalar_indexing(c2)
@. c2 = compute_epsilon(fdtype, x, relstep, absstep, dir)
copyto!(c1,x)
end
copyto!(c3,x)
if fdtype == Val{:forward}
if fdtype == Val(:forward)
@inbounds for i ∈ eachindex(x)
if ArrayInterface.fast_scalar_indexing(c2)
epsilon = ArrayInterface.allowed_getindex(c2,i)*dir
Expand Down Expand Up @@ -168,7 +171,7 @@ function finite_difference_gradient!(
ArrayInterface.allowed_setindex!(c1,c1_old,i)
end
end
elseif fdtype == Val{:central}
elseif fdtype == Val(:central)
@inbounds for i ∈ eachindex(x)
if ArrayInterface.fast_scalar_indexing(c2)
epsilon = ArrayInterface.allowed_getindex(c2,i)*dir
Expand All @@ -191,7 +194,7 @@ function finite_difference_gradient!(
ArrayInterface.allowed_setindex!(c1,c1_old, i)
ArrayInterface.allowed_setindex!(c3,x_old,i)
end
elseif fdtype == Val{:complex} && returntype <: Real
elseif fdtype == Val(:complex) && returntype <: Real
copyto!(c1,x)
epsilon_complex = eps(real(eltype(x)))
# we use c1 here to avoid typing issues with x
Expand Down Expand Up @@ -219,13 +222,13 @@ function finite_difference_gradient!(
# c1 is x1 if we need a complex copy of x, otherwise Nothing
# c2 is Nothing
fx, c1, c2, c3 = cache.fx, cache.c1, cache.c2, cache.c3
if fdtype != Val{:complex}
if fdtype != Val(:complex)
if eltype(df)<:Complex && !(eltype(x)<:Complex)
copyto!(c1,x)
end
end
copyto!(c3,x)
if fdtype == Val{:forward}
if fdtype == Val(:forward)
for i ∈ eachindex(x)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep, dir)
x_old = x[i]
Expand Down Expand Up @@ -262,7 +265,7 @@ function finite_difference_gradient!(
df[i] -= im * imag(dfi)
end
end
elseif fdtype == Val{:central}
elseif fdtype == Val(:central)
@inbounds for i ∈ eachindex(x)
epsilon = compute_epsilon(fdtype, x[i], relstep, absstep, dir)
x_old = x[i]
Expand All @@ -289,7 +292,7 @@ function finite_difference_gradient!(
df[i] -= im*imag(dfi / (2*im*epsilon))
end
end
elseif fdtype==Val{:complex} && returntype<:Real && eltype(df)<:Real && eltype(x)<:Real
elseif fdtype==Val(:complex) && returntype<:Real && eltype(df)<:Real && eltype(x)<:Real
copyto!(c1,x)
epsilon_complex = eps(real(eltype(x)))
# we use c1 here to avoid typing issues with x
Expand Down Expand Up @@ -320,40 +323,40 @@ function finite_difference_gradient!(
# c1 denotes fx1, c2 is fx2, sizes guaranteed by the cache constructor
fx, c1, c2 = cache.fx, cache.c1, cache.c2

if inplace == Val{true}
if inplace == Val(true)
_c1, _c2 = c1, c2
end

if fdtype == Val{:forward}
epsilon = compute_epsilon(Val{:forward}, x, relstep, absstep, dir)
if inplace == Val{true}
if fdtype == Val(:forward)
epsilon = compute_epsilon(Val(:forward), x, relstep, absstep, dir)
if inplace == Val(true)
f(c1, x+epsilon)
else
_c1 = f(x+epsilon)
end
if typeof(fx) != Nothing
@. df = (_c1 - fx) / epsilon
else
if inplace == Val{true}
if inplace == Val(true)
f(c2, x)
else
_c2 = f(x)
end
@. df = (_c1 - _c2) / epsilon
end
elseif fdtype == Val{:central}
epsilon = compute_epsilon(Val{:central}, x, relstep, absstep, dir)
if inplace == Val{true}
elseif fdtype == Val(:central)
epsilon = compute_epsilon(Val(:central), x, relstep, absstep, dir)
if inplace == Val(true)
f(c1, x+epsilon)
f(c2, x-epsilon)
else
_c1 = f(x+epsilon)
_c2 = f(x-epsilon)
end
@. df = (_c1 - _c2) / (2*epsilon)
elseif fdtype == Val{:complex} && returntype <: Real
elseif fdtype == Val(:complex) && returntype <: Real
epsilon_complex = eps(real(eltype(x)))
if inplace == Val{true}
if inplace == Val(true)
f(c1, x+im*epsilon_complex)
else
_c1 = f(x+im*epsilon_complex)
Expand Down
Loading