diff --git a/Project.toml b/Project.toml index 649895d..2b7d77e 100644 --- a/Project.toml +++ b/Project.toml @@ -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" diff --git a/src/derivatives.jl b/src/derivatives.jl index b624d21..3454a6e 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -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) @@ -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? @@ -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)) @@ -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; @@ -87,7 +89,7 @@ function finite_difference_derivative!( df, f, x, - fdtype = Val{:central}, + fdtype = Val(:central), returntype = eltype(x), fx = nothing, epsilon = nothing; @@ -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 @@ -142,10 +144,10 @@ 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 @@ -153,14 +155,14 @@ function finite_difference_derivative!( 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 diff --git a/src/epsilons.jl b/src/epsilons.jl index 1daef79..6be8c8c 100644 --- a/src/epsilons.jl +++ b/src/epsilons.jl @@ -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.") diff --git a/src/gradients.jl b/src/gradients.jl index c13f100..0ddf3c2 100644 --- a/src/gradients.jl +++ b/src/gradients.jl @@ -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 @@ -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 @@ -55,9 +57,9 @@ 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; @@ -65,10 +67,11 @@ function finite_difference_gradient( 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 @@ -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; @@ -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 @@ -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 @@ -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 @@ -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] @@ -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] @@ -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 @@ -320,13 +323,13 @@ 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) @@ -334,16 +337,16 @@ function finite_difference_gradient!( 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 @@ -351,9 +354,9 @@ function finite_difference_gradient!( _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) diff --git a/src/hessians.jl b/src/hessians.jl index 53410da..ada8fef 100644 --- a/src/hessians.jl +++ b/src/hessians.jl @@ -6,20 +6,24 @@ struct HessianCache{T,fdtype,inplace} end function HessianCache(xpp,xpm,xmp,xmm, - fdtype=Val{:hcentral}, - inplace = x isa StaticArray ? Val{false} : Val{true}) + fdtype=Val(:hcentral), + inplace = x isa StaticArray ? Val(false) : Val(true)) + fdtype isa Type && (fdtype = fdtype()) + inplace isa Type && (inplace = inplace()) HessianCache{typeof(xpp),fdtype,inplace}(xpp,xpm,xmp,xmm) end -function HessianCache(x, fdtype=Val{:hcentral}, - inplace = x isa StaticArray ? Val{false} : Val{true}) +function HessianCache(x, fdtype=Val(:hcentral), + inplace = x isa StaticArray ? Val(false) : Val(true)) cx = copy(x) + fdtype isa Type && (fdtype = fdtype()) + inplace isa Type && (inplace = inplace()) HessianCache{typeof(cx),fdtype,inplace}(cx, copy(x), copy(x), copy(x)) end function finite_difference_hessian(f, x, - fdtype = Val{:hcentral}, - inplace = x isa StaticArray ? Val{false} : Val{true}; + fdtype = Val(:hcentral), + inplace = x isa StaticArray ? Val(false) : Val(true); relstep = default_relstep(fdtype, eltype(x)), absstep = relstep) @@ -40,8 +44,8 @@ end function finite_difference_hessian!(H,f, x, - fdtype = Val{:hcentral}, - inplace = x isa StaticArray ? Val{false} : Val{true}; + fdtype = Val(:hcentral), + inplace = x isa StaticArray ? Val(false) : Val(true); relstep=default_relstep(fdtype, eltype(x)), absstep=relstep) @@ -54,20 +58,20 @@ function finite_difference_hessian!(H,f,x, relstep = default_relstep(fdtype, eltype(x)), absstep = relstep) where {T,fdtype,inplace} - @assert fdtype == Val{:hcentral} + @assert fdtype == Val(:hcentral) n = length(x) xpp, xpm, xmp, xmm = cache.xpp, cache.xpm, cache.xmp, cache.xmm fx = f(x) - if inplace === Val{true} + if inplace === Val(true) _xpp, _xpm, _xmp, _xmm = xpp, xpm, xmp, xmm end for i = 1:n xi = ArrayInterface.allowed_getindex(x,i) - epsilon = compute_epsilon(Val{:hcentral}, xi, relstep, absstep) + epsilon = compute_epsilon(Val(:hcentral), xi, relstep, absstep) - if inplace === Val{true} + if inplace === Val(true) ArrayInterface.allowed_setindex!(xpp,xi + epsilon,i) ArrayInterface.allowed_setindex!(xmm,xi - epsilon,i) else @@ -76,11 +80,11 @@ function finite_difference_hessian!(H,f,x, end ArrayInterface.allowed_setindex!(H,(f(_xpp) - 2*fx + f(_xmm)) / epsilon^2,i,i) - epsiloni = compute_epsilon(Val{:central}, xi, relstep, absstep) + epsiloni = compute_epsilon(Val(:central), xi, relstep, absstep) xp = xi + epsiloni xm = xi - epsiloni - if inplace === Val{true} + if inplace === Val(true) ArrayInterface.allowed_setindex!(xpp,xp,i) ArrayInterface.allowed_setindex!(xpm,xp,i) ArrayInterface.allowed_setindex!(xmp,xm,i) @@ -94,11 +98,11 @@ function finite_difference_hessian!(H,f,x, for j = i+1:n xj = ArrayInterface.allowed_getindex(x,j) - epsilonj = compute_epsilon(Val{:central}, xj, relstep, absstep) + epsilonj = compute_epsilon(Val(:central), xj, relstep, absstep) xp = xj + epsilonj xm = xj - epsilonj - if inplace === Val{true} + if inplace === Val(true) ArrayInterface.allowed_setindex!(xpp,xp,j) ArrayInterface.allowed_setindex!(xpm,xm,j) ArrayInterface.allowed_setindex!(xmp,xp,j) @@ -112,7 +116,7 @@ function finite_difference_hessian!(H,f,x, ArrayInterface.allowed_setindex!(H,(f(_xpp) - f(_xpm) - f(_xmp) + f(_xmm))/(4*epsiloni*epsilonj),i,j) - if inplace === Val{true} + if inplace === Val(true) ArrayInterface.allowed_setindex!(xpp,xj,j) ArrayInterface.allowed_setindex!(xpm,xj,j) ArrayInterface.allowed_setindex!(xmp,xj,j) @@ -125,7 +129,7 @@ function finite_difference_hessian!(H,f,x, end end - if inplace === Val{true} + if inplace === Val(true) ArrayInterface.allowed_setindex!(xpp,xi,i) ArrayInterface.allowed_setindex!(xpm,xi,i) ArrayInterface.allowed_setindex!(xmp,xi,i) diff --git a/src/jacobians.jl b/src/jacobians.jl index ce85685..ebcc6c5 100644 --- a/src/jacobians.jl +++ b/src/jacobians.jl @@ -8,13 +8,15 @@ end function JacobianCache( x, - fdtype :: Type{T1} = Val{:forward}, + fdtype :: Union{Val{T1},Type{T1}} = Val(:forward), returntype :: Type{T2} = eltype(x); - inplace :: Type{Val{T3}} = Val{true}, + inplace :: Union{Val{T3},Type{T3}} = Val(true), colorvec = 1:length(x), sparsity = nothing) where {T1,T2,T3} - if eltype(x) <: Real && fdtype==Val{:complex} + fdtype isa Type && (fdtype = fdtype()) + inplace isa Type && (inplace = inplace()) + if eltype(x) <: Real && fdtype==Val(:complex) x1 = false .* im .* x _fx = false .* im .* x else @@ -22,7 +24,7 @@ function JacobianCache( _fx = copy(x) end - if fdtype==Val{:complex} + if fdtype==Val(:complex) _fx1 = nothing else _fx1 = copy(x) @@ -34,25 +36,27 @@ end function JacobianCache( x , fx, - fdtype :: Type{T1} = Val{:forward}, + fdtype :: Union{Val{T1},Type{T1}} = Val(:forward), returntype :: Type{T2} = eltype(x); - inplace :: Type{Val{T3}} = Val{true}, + inplace :: Union{Val{T3},Type{T3}} = Val(true), colorvec = 1:length(x), sparsity = nothing) where {T1,T2,T3} - if eltype(x) <: Real && fdtype==Val{:complex} + fdtype isa Type && (fdtype = fdtype()) + inplace isa Type && (inplace = inplace()) + if eltype(x) <: Real && fdtype==Val(:complex) x1 = false .* im .* x else x1 = copy(x) end - if eltype(fx) <: Real && fdtype==Val{:complex} + if eltype(fx) <: Real && fdtype==Val(:complex) _fx = false .* im .* fx else _fx = copy(fx) end - if fdtype==Val{:complex} + if fdtype==Val(:complex) _fx1 = nothing else _fx1 = copy(fx) @@ -65,13 +69,15 @@ function JacobianCache( x1 , fx , fx1, - fdtype :: Type{T1} = Val{:forward}, + fdtype :: Union{Val{T1},Type{T1}} = Val(:forward), returntype :: Type{T2} = eltype(fx); - inplace :: Type{Val{T3}} = Val{true}, + inplace :: Union{Val{T3},Type{T3}} = Val(true), colorvec = 1:length(x1), sparsity = nothing) where {T1,T2,T3} - if fdtype==Val{:complex} + fdtype isa Type && (fdtype = fdtype()) + inplace isa Type && (inplace = inplace()) + if fdtype==Val(:complex) !(returntype<:Real) && fdtype_error(returntype) if eltype(fx) <: Real @@ -125,7 +131,7 @@ function _make_Ji(::AbstractArray, xtype, dx, color_i, nrows, ncols) end function finite_difference_jacobian(f, x, - fdtype = Val{:forward}, + fdtype = Val(:forward), returntype = eltype(x), f_in = nothing; relstep=default_relstep(fdtype, eltype(x)), @@ -162,9 +168,9 @@ function finite_difference_jacobian( if !(f_in isa Nothing) vecfx = _vec(f_in) - elseif fdtype == Val{:forward} + elseif fdtype == Val(:forward) vecfx = _vec(f(x)) - elseif fdtype == Val{:complex} && returntype <: Real + elseif fdtype == Val(:complex) && returntype <: Real vecfx = real(fx) else vecfx = _vec(fx) @@ -180,11 +186,11 @@ function finite_difference_jacobian( cols_index = [cols_index[i] for i in 1:length(cols_index)] end - if fdtype == Val{:forward} + if fdtype == Val(:forward) function calculate_Ji_forward(i) x_save = ArrayInterface.allowed_getindex(vecx, i) - epsilon = compute_epsilon(Val{:forward}, x_save, relstep, absstep, dir) + epsilon = compute_epsilon(Val(:forward), x_save, relstep, absstep, dir) _vecx1 = Base.setindex(vecx, x_save+epsilon, i) _x1 = reshape(_vecx1, axes(x)) vecfx1 = _vec(f(_x1)) @@ -202,7 +208,7 @@ function finite_difference_jacobian( J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) else tmp = norm(vecx .* (colorvec .== color_i)) - epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) + epsilon = compute_epsilon(Val(:forward), sqrt(tmp), relstep, absstep, dir) _vecx = @. vecx + epsilon * (colorvec == color_i) _x = reshape(_vecx, axes(x)) vecfx1 = _vec(f(_x)) @@ -212,12 +218,12 @@ function finite_difference_jacobian( end end end - elseif fdtype == Val{:central} + elseif fdtype == Val(:central) function calculate_Ji_central(i) x1_save = ArrayInterface.allowed_getindex(vecx1,i) x_save = ArrayInterface.allowed_getindex(vecx,i) - epsilon = compute_epsilon(Val{:forward}, x1_save, relstep, absstep, dir) + epsilon = compute_epsilon(Val(:forward), x1_save, relstep, absstep, dir) _vecx1 = Base.setindex(vecx1,x1_save+epsilon,i) _vecx = Base.setindex(vecx,x_save-epsilon,i) _x1 = reshape(_vecx1, axes(x)) @@ -238,7 +244,7 @@ function finite_difference_jacobian( J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) else tmp = norm(vecx1 .* (colorvec .== color_i)) - epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) + epsilon = compute_epsilon(Val(:forward), sqrt(tmp), relstep, absstep, dir) _vecx1 = @. vecx1 + epsilon * (colorvec == color_i) _vecx = @. vecx - epsilon * (colorvec == color_i) _x1 = reshape(_vecx1, axes(x)) @@ -251,7 +257,7 @@ function finite_difference_jacobian( end end end - elseif fdtype == Val{:complex} && returntype <: Real + elseif fdtype == Val(:complex) && returntype <: Real epsilon = eps(eltype(x)) function calculate_Ji_complex(i) @@ -290,14 +296,14 @@ end function finite_difference_jacobian!(J, f, x, - fdtype = Val{:forward}, + fdtype = Val(:forward), returntype = eltype(x), f_in = nothing; relstep=default_relstep(fdtype, eltype(x)), absstep=relstep, colorvec = 1:length(x), sparsity = ArrayInterface.has_sparsestruct(J) ? J : nothing) - if f_in isa Nothing && fdtype == Val{:forward} + if f_in isa Nothing && fdtype == Val(:forward) if size(J,1) == length(x) fx = zero(x) else @@ -342,7 +348,7 @@ function finite_difference_jacobian!( fill!(J,false) end - if fdtype == Val{:forward} + if fdtype == Val(:forward) vfx1 = _vec(fx1) if f_in isa Nothing @@ -355,7 +361,7 @@ function finite_difference_jacobian!( @inbounds for color_i ∈ 1:maximum(colorvec) if sparsity isa Nothing x1_save = ArrayInterface.allowed_getindex(x1,color_i) - epsilon = compute_epsilon(Val{:forward}, x1_save, relstep, absstep, dir) + epsilon = compute_epsilon(Val(:forward), x1_save, relstep, absstep, dir) ArrayInterface.allowed_setindex!(x1, x1_save + epsilon, color_i) f(fx1, x1) # J is dense, so either it is truly dense or this is the @@ -366,7 +372,7 @@ function finite_difference_jacobian!( else # Perturb along the colorvec vector @. fx1 = x1 * (_color == color_i) tmp = norm(fx1) - epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) + epsilon = compute_epsilon(Val(:forward), sqrt(tmp), relstep, absstep, dir) @. x1 = x1 + epsilon * (_color == color_i) f(fx1, x1) # J is a sparse matrix, so decompress on the fly @@ -390,12 +396,12 @@ function finite_difference_jacobian!( @. x1 = x1 - epsilon * (_color == color_i) end end #for ends here - elseif fdtype == Val{:central} + elseif fdtype == Val(:central) vfx1 = _vec(fx1) @inbounds for color_i ∈ 1:maximum(colorvec) if sparsity isa Nothing x_save = ArrayInterface.allowed_getindex(x, color_i) - epsilon = compute_epsilon(Val{:central}, x_save, relstep, absstep, dir) + epsilon = compute_epsilon(Val(:central), x_save, relstep, absstep, dir) ArrayInterface.allowed_setindex!(x1, x_save + epsilon, color_i) f(fx1, x1) ArrayInterface.allowed_setindex!(x1, x_save - epsilon, color_i) @@ -405,7 +411,7 @@ function finite_difference_jacobian!( else # Perturb along the colorvec vector @. fx1 = x1 * (_color == color_i) tmp = norm(fx1) - epsilon = compute_epsilon(Val{:central}, sqrt(tmp), relstep, absstep, dir) + epsilon = compute_epsilon(Val(:central), sqrt(tmp), relstep, absstep, dir) @. x1 = x1 + epsilon * (_color == color_i) @. x = x - epsilon * (_color == color_i) f(fx1, x1) @@ -424,7 +430,7 @@ function finite_difference_jacobian!( @. x = x + epsilon * (_color == color_i) end end - elseif fdtype==Val{:complex} && returntype<:Real + elseif fdtype==Val(:complex) && returntype<:Real epsilon = eps(eltype(x)) @inbounds for color_i ∈ 1:maximum(colorvec) if sparsity isa Nothing diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml index b28c55b..4905387 100644 --- a/test/downstream/Project.toml +++ b/test/downstream/Project.toml @@ -1,2 +1,3 @@ [deps] +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"