diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 56350777..a7536455 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -21,7 +21,7 @@ jobs: - '1.0' - '1.6' - '1' - - 'nightly' + # - 'nightly' os: - ubuntu-latest arch: diff --git a/.gitignore b/.gitignore index 78fb62c9..69feffb1 100644 --- a/.gitignore +++ b/.gitignore @@ -4,5 +4,6 @@ *.DS_Store /docs/build/ /docs/site/ +/docs/Manifest.toml /benchmark_data/ /Manifest.toml diff --git a/Project.toml b/Project.toml index f7c66517..2ee06921 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ CommonSubexpressions = "0.3" DiffResults = "0.0.1, 0.0.2, 0.0.3, 0.0.4, 1.0.1" DiffRules = "1.4.0" DiffTests = "0.0.1, 0.1" +IrrationalConstants = "0.1, 0.2" LogExpFunctions = "0.3" NaNMath = "0.2.2, 0.3, 1" Preferences = "1" @@ -35,12 +36,13 @@ ForwardDiffStaticArraysExt = "StaticArrays" Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" DiffTests = "de460e47-3fe3-5279-bb4a-814414816d5d" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Calculus", "DiffTests", "SparseArrays", "Test", "InteractiveUtils", "StaticArrays"] +test = ["Calculus", "DiffTests", "IrrationalConstants", "SparseArrays", "Test", "InteractiveUtils", "StaticArrays"] [weakdeps] -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" \ No newline at end of file +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/docs/Project.toml b/docs/Project.toml index 050425cf..9b9a6718 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + +[compat] +Documenter = "1" diff --git a/docs/make.jl b/docs/make.jl index e653afca..4bbc19af 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,7 +11,8 @@ makedocs(modules=[ForwardDiff], "Upgrading from Older Versions" => "user/upgrade.md"], "Developer Documentation" => [ "How ForwardDiff Works" => "dev/how_it_works.md", - "How to Contribute" => "dev/contributing.md"]]) + "How to Contribute" => "dev/contributing.md"]], + checkdocs=:exports) deploydocs( repo = "github.com/JuliaDiff/ForwardDiff.jl.git" diff --git a/ext/ForwardDiffStaticArraysExt.jl b/ext/ForwardDiffStaticArraysExt.jl index cb9222ff..25b73fb8 100644 --- a/ext/ForwardDiffStaticArraysExt.jl +++ b/ext/ForwardDiffStaticArraysExt.jl @@ -7,7 +7,7 @@ using ForwardDiff: Dual, partials, GradientConfig, JacobianConfig, HessianConfig gradient, hessian, jacobian, gradient!, hessian!, jacobian!, extract_gradient!, extract_jacobian!, extract_value!, vector_mode_gradient, vector_mode_gradient!, - vector_mode_jacobian, vector_mode_jacobian!, valtype, value, _lyap_div! + vector_mode_jacobian, vector_mode_jacobian!, valtype, value using DiffResults: DiffResult, ImmutableDiffResult, MutableDiffResult @generated function dualize(::Type{T}, x::StaticArray) where T @@ -23,27 +23,25 @@ end @inline static_dual_eval(::Type{T}, f, x::StaticArray) where T = f(dualize(T, x)) +# To fix method ambiguity issues: function LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix}) where {Tg,T<:Real,N} - λ,Q = eigen(Symmetric(value.(parent(A)))) - parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N) - Dual{Tg}.(λ, tuple.(parts...)) + return ForwardDiff._eigvals(A) end - function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}, <:StaticArrays.StaticMatrix}) where {Tg,T<:Real,N} - λ = eigvals(A) - _,Q = eigen(Symmetric(value.(parent(A)))) - parts = ntuple(j -> Q*ForwardDiff._lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N) - Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) + return ForwardDiff._eigen(A) end +# For `MMatrix` we can use the in-place method +ForwardDiff._lyap_div!!(A::StaticArrays.MMatrix, λ::AbstractVector) = ForwardDiff._lyap_div!(A, λ) + # Gradient -@inline ForwardDiff.gradient(f, x::StaticArray) = vector_mode_gradient(f, x) -@inline ForwardDiff.gradient(f, x::StaticArray, cfg::GradientConfig) = gradient(f, x) -@inline ForwardDiff.gradient(f, x::StaticArray, cfg::GradientConfig, ::Val) = gradient(f, x) +@inline ForwardDiff.gradient(f::F, x::StaticArray) where F = vector_mode_gradient(f, x) +@inline ForwardDiff.gradient(f::F, x::StaticArray, cfg::GradientConfig) where F = gradient(f, x) +@inline ForwardDiff.gradient(f::F, x::StaticArray, cfg::GradientConfig, ::Val) where F = gradient(f, x) -@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray) = vector_mode_gradient!(result, f, x) -@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::GradientConfig) = gradient!(result, f, x) -@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::GradientConfig, ::Val) = gradient!(result, f, x) +@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f::F, x::StaticArray) where F = vector_mode_gradient!(result, f, x) +@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f::F, x::StaticArray, cfg::GradientConfig) where F = gradient!(result, f, x) +@inline ForwardDiff.gradient!(result::Union{AbstractArray,DiffResult}, f::F, x::StaticArray, cfg::GradientConfig, ::Val) where F = gradient!(result, f, x) @generated function extract_gradient(::Type{T}, y::Real, x::S) where {T,S<:StaticArray} result = Expr(:tuple, [:(partials(T, y, $i)) for i in 1:length(x)]...) @@ -65,13 +63,13 @@ end end # Jacobian -@inline ForwardDiff.jacobian(f, x::StaticArray) = vector_mode_jacobian(f, x) -@inline ForwardDiff.jacobian(f, x::StaticArray, cfg::JacobianConfig) = jacobian(f, x) -@inline ForwardDiff.jacobian(f, x::StaticArray, cfg::JacobianConfig, ::Val) = jacobian(f, x) +@inline ForwardDiff.jacobian(f::F, x::StaticArray) where F = vector_mode_jacobian(f, x) +@inline ForwardDiff.jacobian(f::F, x::StaticArray, cfg::JacobianConfig) where F = jacobian(f, x) +@inline ForwardDiff.jacobian(f::F, x::StaticArray, cfg::JacobianConfig, ::Val) where F = jacobian(f, x) -@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray) = vector_mode_jacobian!(result, f, x) -@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::JacobianConfig) = jacobian!(result, f, x) -@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f, x::StaticArray, cfg::JacobianConfig, ::Val) = jacobian!(result, f, x) +@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f::F, x::StaticArray) where F = vector_mode_jacobian!(result, f, x) +@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f::F, x::StaticArray, cfg::JacobianConfig) where F = jacobian!(result, f, x) +@inline ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f::F, x::StaticArray, cfg::JacobianConfig, ::Val) where F = jacobian!(result, f, x) @generated function extract_jacobian(::Type{T}, ydual::StaticArray, x::S) where {T,S<:StaticArray} M, N = length(ydual), length(x) @@ -110,18 +108,18 @@ end end # Hessian -ForwardDiff.hessian(f, x::StaticArray) = jacobian(y -> gradient(f, y), x) -ForwardDiff.hessian(f, x::StaticArray, cfg::HessianConfig) = hessian(f, x) -ForwardDiff.hessian(f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian(f, x) +ForwardDiff.hessian(f::F, x::StaticArray) where F = jacobian(y -> gradient(f, y), x) +ForwardDiff.hessian(f::F, x::StaticArray, cfg::HessianConfig) where F = hessian(f, x) +ForwardDiff.hessian(f::F, x::StaticArray, cfg::HessianConfig, ::Val) where F = hessian(f, x) -ForwardDiff.hessian!(result::AbstractArray, f, x::StaticArray) = jacobian!(result, y -> gradient(f, y), x) +ForwardDiff.hessian!(result::AbstractArray, f::F, x::StaticArray) where F = jacobian!(result, y -> gradient(f, y), x) -ForwardDiff.hessian!(result::MutableDiffResult, f, x::StaticArray) = hessian!(result, f, x, HessianConfig(f, result, x)) +ForwardDiff.hessian!(result::MutableDiffResult, f::F, x::StaticArray) where F = hessian!(result, f, x, HessianConfig(f, result, x)) -ForwardDiff.hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig) = hessian!(result, f, x) -ForwardDiff.hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian!(result, f, x) +ForwardDiff.hessian!(result::ImmutableDiffResult, f::F, x::StaticArray, cfg::HessianConfig) where F = hessian!(result, f, x) +ForwardDiff.hessian!(result::ImmutableDiffResult, f::F, x::StaticArray, cfg::HessianConfig, ::Val) where F = hessian!(result, f, x) -function ForwardDiff.hessian!(result::ImmutableDiffResult, f, x::StaticArray) +function ForwardDiff.hessian!(result::ImmutableDiffResult, f::F, x::StaticArray) where F T = typeof(Tag(f, eltype(x))) d1 = dualize(T, x) d2 = dualize(T, d1) diff --git a/src/ForwardDiff.jl b/src/ForwardDiff.jl index a27d6dba..b2670451 100644 --- a/src/ForwardDiff.jl +++ b/src/ForwardDiff.jl @@ -7,7 +7,12 @@ if VERSION >= v"1.6" end using Random using LinearAlgebra - +if VERSION < v"1.2.0-DEV.125" # 1da48c2e4028c1514ed45688be727efbef1db884 + require_one_based_indexing(A...) = !Base.has_offset_axes(A...) || throw(ArgumentError( + "offset arrays are not supported but got an array with index other than 1")) +else + using Base: require_one_based_indexing +end import Printf import NaNMath import SpecialFunctions diff --git a/src/derivative.jl b/src/derivative.jl index 5d7e5dab..b39e2a48 100644 --- a/src/derivative.jl +++ b/src/derivative.jl @@ -22,8 +22,9 @@ stored in `y`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -@inline function derivative(f!, y::AbstractArray, x::Real, - cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK} +@inline function derivative(f!::F, y::AbstractArray, x::Real, + cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F, T, CHK} + require_one_based_indexing(y) CHK && checktag(T, f!, x) ydual = cfg.duals seed!(ydual, y) @@ -42,6 +43,7 @@ This method assumes that `isa(f(x), Union{Real,AbstractArray})`. """ @inline function derivative!(result::Union{AbstractArray,DiffResult}, f::F, x::R) where {F,R<:Real} + result isa DiffResult || require_one_based_indexing(result) T = typeof(Tag(f, R)) ydual = f(Dual{T}(x, one(x))) result = extract_value!(T, result, ydual) @@ -58,8 +60,9 @@ called as `f!(y, x)` where the result is stored in `y`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ @inline function derivative!(result::Union{AbstractArray,DiffResult}, - f!, y::AbstractArray, x::Real, - cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK} + f!::F, y::AbstractArray, x::Real, + cfg::DerivativeConfig{T} = DerivativeConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F, T, CHK} + result isa DiffResult ? require_one_based_indexing(y) : require_one_based_indexing(result, y) CHK && checktag(T, f!, x) ydual = cfg.duals seed!(ydual, y) diff --git a/src/dual.jl b/src/dual.jl index 5afb2144..cc78b59b 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -78,6 +78,10 @@ end @inline Dual{T,V,N}(x::Number) where {T,V,N} = convert(Dual{T,V,N}, x) @inline Dual{T,V}(x) where {T,V} = convert(Dual{T,V}, x) +# Fix method ambiguity issue by adapting the definition in Base to `Dual`s +Dual{T,V,N}(x::Base.TwicePrecision) where {T,V,N} = + (Dual{T,V,N}(x.hi) + Dual{T,V,N}(x.lo))::Dual{T,V,N} + ############################## # Utility/Accessor Functions # ############################## @@ -340,7 +344,6 @@ else Base.div(x::Dual, y::Dual) = div(value(x), value(y)) end -Base.hash(d::Dual) = hash(value(d)) Base.hash(d::Dual, hsh::UInt) = hash(value(d), hsh) function Base.read(io::IO, ::Type{Dual{T,V,N}}) where {T,V,N} @@ -416,7 +419,7 @@ function Base.promote_rule(::Type{Dual{T,A,N}}, return Dual{T,promote_type(A, B),N} end -for R in (Irrational, Real, BigFloat, Bool) +for R in (AbstractIrrational, Real, BigFloat, Bool) if isconcretetype(R) # issue #322 @eval begin Base.promote_rule(::Type{$R}, ::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T,promote_type($R, V),N} @@ -703,7 +706,11 @@ end # Symmetric eigvals # #-------------------# -function LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} +# To be able to reuse this default definition in the StaticArrays extension +# (has to be re-defined to avoid method ambiguity issues) +# we forward the call to an internal method that can be shared and reused +LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigvals(A) +function _eigvals(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} λ,Q = eigen(Symmetric(value.(parent(A)))) parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N) Dual{Tg}.(λ, tuple.(parts...)) @@ -721,8 +728,19 @@ function LinearAlgebra.eigvals(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:R Dual{Tg}.(λ, tuple.(parts...)) end -# A ./ (λ - λ') but with diag special cased -function _lyap_div!(A, λ) +# A ./ (λ' .- λ) but with diag special cased +# Default out-of-place method +function _lyap_div!!(A::AbstractMatrix, λ::AbstractVector) + return map( + (a, b, idx) -> a / (idx[1] == idx[2] ? oneunit(b) : b), + A, + λ' .- λ, + CartesianIndices(A), + ) +end +# For `Matrix` (and e.g. `StaticArrays.MMatrix`) we can use an in-place method +_lyap_div!!(A::Matrix, λ::AbstractVector) = _lyap_div!(A, λ) +function _lyap_div!(A::AbstractMatrix, λ::AbstractVector) for (j,μ) in enumerate(λ), (k,λ) in enumerate(λ) if k ≠ j A[k,j] /= μ - λ @@ -731,17 +749,21 @@ function _lyap_div!(A, λ) A end -function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} +# To be able to reuse this default definition in the StaticArrays extension +# (has to be re-defined to avoid method ambiguity issues) +# we forward the call to an internal method that can be shared and reused +LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} = _eigen(A) +function _eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} λ = eigvals(A) _,Q = eigen(Symmetric(value.(parent(A)))) - parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N) + parts = ntuple(j -> Q*_lyap_div!!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N) Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) end function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N} λ = eigvals(A) _,Q = eigen(SymTridiagonal(value.(parent(A)))) - parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N) + parts = ntuple(j -> Q*_lyap_div!!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N) Eigen(λ,Dual{Tg}.(Q, tuple.(parts...))) end diff --git a/src/gradient.jl b/src/gradient.jl index 614bea67..6d4cc791 100644 --- a/src/gradient.jl +++ b/src/gradient.jl @@ -13,7 +13,8 @@ This method assumes that `isa(f(x), Real)`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function gradient(f, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK} +function gradient(f::F, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {F, T, CHK} + require_one_based_indexing(x) CHK && checktag(T, f, x) if chunksize(cfg) == length(x) return vector_mode_gradient(f, x, cfg) @@ -32,6 +33,7 @@ This method assumes that `isa(f(x), Real)`. """ function gradient!(result::Union{AbstractArray,DiffResult}, f::F, x::AbstractArray, cfg::GradientConfig{T} = GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK, F} + result isa DiffResult ? require_one_based_indexing(x) : require_one_based_indexing(result, x) CHK && checktag(T, f, x) if chunksize(cfg) == length(x) vector_mode_gradient!(result, f, x, cfg) diff --git a/src/hessian.jl b/src/hessian.jl index 9a9c4f77..9c755c9a 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -11,7 +11,8 @@ This method assumes that `isa(f(x), Real)`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function hessian(f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK} +function hessian(f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F, T,CHK} + require_one_based_indexing(x) CHK && checktag(T, f, x) ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}()) return jacobian(∇f, x, cfg.jacobian_config, Val{false}()) @@ -27,7 +28,8 @@ This method assumes that `isa(f(x), Real)`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK} +function hessian!(result::AbstractArray, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} + require_one_based_indexing(result, x) CHK && checktag(T, f, x) ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}()) jacobian!(result, ∇f, x, cfg.jacobian_config, Val{false}()) @@ -61,7 +63,7 @@ because `isa(result, DiffResult)`, `cfg` is constructed as `HessianConfig(f, res Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {T,CHK} +function hessian!(result::DiffResult, f::F, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} CHK && checktag(T, f, x) ∇f! = InnerGradientForHess(result, cfg, f) jacobian!(DiffResults.hessian(result), ∇f!, DiffResults.gradient(result), x, cfg.jacobian_config, Val{false}()) diff --git a/src/jacobian.jl b/src/jacobian.jl index cd417a79..e7ca96f5 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -15,7 +15,8 @@ This method assumes that `isa(f(x), AbstractArray)`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function jacobian(f, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK} +function jacobian(f::F, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} + require_one_based_indexing(x) CHK && checktag(T, f, x) if chunksize(cfg) == length(x) return vector_mode_jacobian(f, x, cfg) @@ -32,7 +33,8 @@ stored in `y`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function jacobian(f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T, CHK} +function jacobian(f!::F, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F,T, CHK} + require_one_based_indexing(y, x) CHK && checktag(T, f!, x) if chunksize(cfg) == length(x) return vector_mode_jacobian(f!, y, x, cfg) @@ -52,7 +54,8 @@ This method assumes that `isa(f(x), AbstractArray)`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function jacobian!(result::Union{AbstractArray,DiffResult}, f, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK} +function jacobian!(result::Union{AbstractArray,DiffResult}, f::F, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x), ::Val{CHK}=Val{true}()) where {F,T, CHK} + result isa DiffResult ? require_one_based_indexing(x) : require_one_based_indexing(result, x) CHK && checktag(T, f, x) if chunksize(cfg) == length(x) vector_mode_jacobian!(result, f, x, cfg) @@ -72,7 +75,8 @@ This method assumes that `isa(f(x), AbstractArray)`. Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care. """ -function jacobian!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T,CHK} +function jacobian!(result::Union{AbstractArray,DiffResult}, f!::F, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {F,T,CHK} + result isa DiffResult ? require_one_based_indexing(y, x) : require_one_based_indexing(result, y, x) CHK && checktag(T, f!, x) if chunksize(cfg) == length(x) vector_mode_jacobian!(result, f!, y, x, cfg) diff --git a/src/partials.jl b/src/partials.jl index fce67b0a..a5316e3e 100644 --- a/src/partials.jl +++ b/src/partials.jl @@ -52,8 +52,7 @@ Base.:(==)(a::Partials{N}, b::Partials{N}) where {N} = a.values == b.values const PARTIALS_HASH = hash(Partials) -Base.hash(partials::Partials) = hash(partials.values, PARTIALS_HASH) -Base.hash(partials::Partials, hsh::UInt64) = hash(hash(partials), hsh) +Base.hash(partials::Partials, hsh::UInt64) = hash(hash(partials.values, PARTIALS_HASH), hsh) @inline Base.copy(partials::Partials) = partials diff --git a/src/prelude.jl b/src/prelude.jl index 1a5d0cfa..5fe0a9c1 100644 --- a/src/prelude.jl +++ b/src/prelude.jl @@ -14,12 +14,9 @@ const BINARY_PREDICATES = Symbol[:isequal, :isless, :<, :>, :(==), :(!=), :(<=), struct Chunk{N} end -const CHUNKS = [Chunk{i}() for i in 1:DEFAULT_CHUNK_THRESHOLD] - function Chunk(input_length::Integer, threshold::Integer = DEFAULT_CHUNK_THRESHOLD) N = pickchunksize(input_length, threshold) - 0 < N <= DEFAULT_CHUNK_THRESHOLD && return CHUNKS[N] - return Chunk{N}() + Base.@nif 12 d->(N == d) d->(Chunk{d}()) d->(Chunk{N}()) end function Chunk(x::AbstractArray, threshold::Integer = DEFAULT_CHUNK_THRESHOLD) diff --git a/test/AllocationsTest.jl b/test/AllocationsTest.jl index 2a0075a7..ad1832d0 100644 --- a/test/AllocationsTest.jl +++ b/test/AllocationsTest.jl @@ -24,17 +24,25 @@ convert_test_574() = convert(ForwardDiff.Dual{Nothing,ForwardDiff.Dual{Nothing,F index = 1 alloc = @allocated ForwardDiff.seed!(duals, x, index, seeds) alloc = @allocated ForwardDiff.seed!(duals, x, index, seeds) - @test alloc == 0 + if VERSION < v"1.9" || VERSION >= v"1.11" + @test alloc == 0 + else + @test_broken alloc == 0 + end index = 1 alloc = @allocated ForwardDiff.seed!(duals, x, index, seed) alloc = @allocated ForwardDiff.seed!(duals, x, index, seed) - @test alloc == 0 - + if VERSION < v"1.9" || VERSION >= v"1.11" + @test alloc == 0 + else + @test_broken alloc == 0 + end + alloc = @allocated convert_test_574() alloc = @allocated convert_test_574() @test alloc == 0 - + end end diff --git a/test/DerivativeTest.jl b/test/DerivativeTest.jl index dfdd8ed2..4b7463c8 100644 --- a/test/DerivativeTest.jl +++ b/test/DerivativeTest.jl @@ -17,8 +17,7 @@ Random.seed!(1) const x = 1 -for f in DiffTests.NUMBER_TO_NUMBER_FUNCS - println(" ...testing $f") +@testset "$f" for f in DiffTests.NUMBER_TO_NUMBER_FUNCS v = f(x) d = ForwardDiff.derivative(f, x) @test isapprox(d, Calculus.derivative(f, x), atol=FINITEDIFF_ERROR) @@ -29,8 +28,7 @@ for f in DiffTests.NUMBER_TO_NUMBER_FUNCS @test isapprox(DiffResults.derivative(out), d) end -for f in DiffTests.NUMBER_TO_ARRAY_FUNCS - println(" ...testing $f") +@testset "$f" for f in DiffTests.NUMBER_TO_ARRAY_FUNCS v = f(x) d = ForwardDiff.derivative(f, x) @@ -47,8 +45,7 @@ for f in DiffTests.NUMBER_TO_ARRAY_FUNCS @test isapprox(DiffResults.derivative(out), d) end -for f! in DiffTests.INPLACE_NUMBER_TO_ARRAY_FUNCS - println(" ...testing $f!") +@testset "$(f!)" for f! in DiffTests.INPLACE_NUMBER_TO_ARRAY_FUNCS m, n = 3, 2 y = fill(0.0, m, n) f = x -> (tmp = similar(y, promote_type(eltype(y), typeof(x)), m, n); f!(tmp, x); tmp) diff --git a/test/DualTest.jl b/test/DualTest.jl index 8d87f90b..786d3bda 100644 --- a/test/DualTest.jl +++ b/test/DualTest.jl @@ -30,8 +30,7 @@ ForwardDiff.:≺(::Int, ::Type{TestTag()}) = false ForwardDiff.:≺(::Type{TestTag}, ::Type{OuterTestTag}) = true ForwardDiff.:≺(::Type{OuterTestTag}, ::Type{TestTag}) = false -for N in (0,3), M in (0,4), V in (Int, Float32) - println(" ...testing Dual{TestTag(),$V,$N} and Dual{TestTag(),Dual{TestTag(),$V,$M},$N}") +@testset "Dual{Z,$V,$N} and Dual{Z,Dual{Z,$V,$M},$N}" for N in (0,3), M in (0,4), V in (Int, Float32) PARTIALS = Partials{N,V}(ntuple(n -> intrand(V), N)) PRIMAL = intrand(V) @@ -446,13 +445,12 @@ for N in (0,3), M in (0,4), V in (Int, Float32) @test abs(NESTED_FDNUM) === NESTED_FDNUM if V != Int - for (M, f, arity) in DiffRules.diffrules(filter_modules = nothing) + @testset "$(M).$(f) with $arity arguments" for (M, f, arity) in DiffRules.diffrules(filter_modules = nothing) if f in (:/, :rem2pi) continue # Skip these rules elseif !(isdefined(@__MODULE__, M) && isdefined(getfield(@__MODULE__, M), f)) continue # Skip rules for methods not defined in the current scope end - println(" ...auto-testing $(M).$(f) with $arity arguments") if arity == 1 deriv = DiffRules.diffrule(M, f, :x) modifier = if in(f, (:asec, :acsc, :asecd, :acscd, :acosh, :acoth)) @@ -637,4 +635,8 @@ end @test ForwardDiff.derivative(float, 1)::Float64 === 1.0 end +@testset "TwicePrecision" begin + @test ForwardDiff.derivative(x -> sum(1 .+ x .* (0:0.1:1)), 1) == 5.5 +end + end # module diff --git a/test/GradientTest.jl b/test/GradientTest.jl index 707fcd68..d806a967 100644 --- a/test/GradientTest.jl +++ b/test/GradientTest.jl @@ -19,8 +19,7 @@ x = [0.1, 0.2, 0.3] v = f(x) g = [-9.4, 15.6, 52.0] -for c in (1, 2, 3), tag in (nothing, Tag(f, eltype(x))) - println(" ...running hardcoded test with chunk size = $c and tag = $(repr(tag))") +@testset "Rosenbrock, chunk size = $c and tag = $(repr(tag))" for c in (1, 2, 3), tag in (nothing, Tag(f, eltype(x))) cfg = ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{c}(), tag) @test eltype(cfg) == Dual{typeof(tag), eltype(x), c} @@ -59,8 +58,7 @@ for f in DiffTests.VECTOR_TO_NUMBER_FUNCS v = f(X) g = ForwardDiff.gradient(f, X) @test isapprox(g, Calculus.gradient(f, X), atol=FINITEDIFF_ERROR) - for c in CHUNK_SIZES, tag in (nothing, Tag(f, eltype(x))) - println(" ...testing $f with chunk size = $c and tag = $(repr(tag))") + @testset "... with chunk size = $c and tag = $(repr(tag))" for c in CHUNK_SIZES, tag in (nothing, Tag(f, eltype(x))) cfg = ForwardDiff.GradientConfig(f, X, ForwardDiff.Chunk{c}(), tag) out = ForwardDiff.gradient(f, X, cfg) @@ -81,11 +79,9 @@ end # test specialized StaticArray codepaths # ########################################## -println(" ...testing specialized StaticArray codepaths") +@testset "Specialized StaticArray codepaths: $T" for T in (StaticArrays.SArray, StaticArrays.MArray) + x = rand(3, 3) -x = rand(3, 3) - -for T in (StaticArrays.SArray, StaticArrays.MArray) sx = T{Tuple{3,3}}(x) cfg = ForwardDiff.GradientConfig(nothing, x) @@ -139,6 +135,9 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.gradient(sresult1) == DiffResults.gradient(result) @test DiffResults.gradient(sresult2) == DiffResults.gradient(result) @test DiffResults.gradient(sresult3) == DiffResults.gradient(result) + + # make sure this is not a source of type instability + @inferred ForwardDiff.GradientConfig(f, sx) end @testset "exponential function at base zero" begin @@ -166,7 +165,7 @@ end function f(p) sum(collect(0.0:p[1]:p[2])) end - @test ForwardDiff.gradient(f, [0.2,25.0]) == [7875.0, 0.0] + @test ForwardDiff.gradient(f, [0.3, 25.0]) == [3486.0, 0.0] end end # module diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 1df4232b..64c289f2 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -22,8 +22,7 @@ h = [-66.0 -40.0 0.0; -40.0 130.0 -80.0; 0.0 -80.0 200.0] -for c in HESSIAN_CHUNK_SIZES, tag in (nothing, Tag((f,ForwardDiff.gradient), eltype(x))) - println(" ...running hardcoded test with chunk size = $c and tag = $(repr(tag))") +@testset "running hardcoded test with chunk size = $c and tag = $(repr(tag))" for c in HESSIAN_CHUNK_SIZES, tag in (nothing, Tag((f,ForwardDiff.gradient), eltype(x))) cfg = ForwardDiff.HessianConfig(f, x, ForwardDiff.Chunk{c}(), tag) resultcfg = ForwardDiff.HessianConfig(f, DiffResults.HessianResult(x), x, ForwardDiff.Chunk{c}(), tag) @@ -68,8 +67,7 @@ for f in DiffTests.VECTOR_TO_NUMBER_FUNCS h = ForwardDiff.hessian(f, X) # finite difference approximation error is really bad for Hessians... @test isapprox(h, Calculus.hessian(f, X), atol=0.02) - for c in HESSIAN_CHUNK_SIZES, tag in (nothing, Tag((f,ForwardDiff.gradient), eltype(x))) - println(" ...testing $f with chunk size = $c and tag = $(repr(tag))") + @testset "$f with chunk size = $c and tag = $(repr(tag))" for c in HESSIAN_CHUNK_SIZES, tag in (nothing, Tag((f,ForwardDiff.gradient), eltype(x))) cfg = ForwardDiff.HessianConfig(f, X, ForwardDiff.Chunk{c}(), tag) resultcfg = ForwardDiff.HessianConfig(f, DiffResults.HessianResult(X), X, ForwardDiff.Chunk{c}(), tag) @@ -92,7 +90,7 @@ end # test specialized StaticArray codepaths # ########################################## -println(" ...testing specialized StaticArray codepaths") +@info "testing specialized StaticArray codepaths" x = rand(3, 3) for T in (StaticArrays.SArray, StaticArrays.MArray) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 69f9409c..0c24b018 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -103,11 +103,10 @@ for f in DiffTests.ARRAY_TO_ARRAY_FUNCS v = f(X) j = ForwardDiff.jacobian(f, X) @test isapprox(j, Calculus.jacobian(x -> vec(f(x)), X, :forward), atol=1.3FINITEDIFF_ERROR) - for c in CHUNK_SIZES, tag in (nothing, Tag) + @testset "$f with chunk size = $c and tag = $(repr(tag))" for c in CHUNK_SIZES, tag in (nothing, Tag) if tag == Tag tag = Tag(f, eltype(X)) end - println(" ...testing $f with chunk size = $c and tag = $(repr(tag))") cfg = JacobianConfig(f, X, ForwardDiff.Chunk{c}(), tag) out = ForwardDiff.jacobian(f, X, cfg) @@ -129,8 +128,7 @@ for f! in DiffTests.INPLACE_ARRAY_TO_ARRAY_FUNCS f!(v, X) j = ForwardDiff.jacobian(f!, fill!(similar(Y), 0.0), X) @test isapprox(j, Calculus.jacobian(x -> (y = fill!(similar(Y), 0.0); f!(y, x); vec(y)), X, :forward), atol=FINITEDIFF_ERROR) - for c in CHUNK_SIZES, tag in (nothing, Tag(f!, eltype(X))) - println(" ...testing $(f!) with chunk size = $c and tag = $(repr(tag))") + @testset "$(f!) with chunk size = $c and tag = $(repr(tag))" for c in CHUNK_SIZES, tag in (nothing, Tag(f!, eltype(X))) ycfg = JacobianConfig(f!, fill!(similar(Y), 0.0), X, ForwardDiff.Chunk{c}(), tag) y = fill!(similar(Y), 0.0) @@ -164,7 +162,7 @@ end # test specialized StaticArray codepaths # ########################################## -println(" ...testing specialized StaticArray codepaths") +@info "testing specialized StaticArray codepaths" x = rand(3, 3) for T in (StaticArrays.SArray, StaticArrays.MArray) @@ -224,6 +222,9 @@ for T in (StaticArrays.SArray, StaticArrays.MArray) @test DiffResults.jacobian(sresult1) == DiffResults.jacobian(result) @test DiffResults.jacobian(sresult2) == DiffResults.jacobian(result) @test DiffResults.jacobian(sresult3) == DiffResults.jacobian(result) + + # make sure this is not a source of type instability + @inferred ForwardDiff.JacobianConfig(f, sx) end @testset "dimension errors for jacobian" begin @@ -235,14 +236,18 @@ end @testset "eigen" begin @test ForwardDiff.jacobian(x -> eigvals(SymTridiagonal(x, x[1:end-1])), [1.,2.]) ≈ [(1 - 3/sqrt(5))/2 (1 - 1/sqrt(5))/2 ; (1 + 3/sqrt(5))/2 (1 + 1/sqrt(5))/2] @test ForwardDiff.jacobian(x -> eigvals(Symmetric(x*x')), [1.,2.]) ≈ [0 0; 2 4] - + x0 = [1.0, 2.0]; ev1(x) = eigen(Symmetric(x*x')).vectors[:,1] @test ForwardDiff.jacobian(ev1, x0) ≈ Calculus.finite_difference_jacobian(ev1, x0) ev2(x) = eigen(SymTridiagonal(x, x[1:end-1])).vectors[:,1] @test ForwardDiff.jacobian(ev2, x0) ≈ Calculus.finite_difference_jacobian(ev2, x0) - x0_static = SVector{2}(x0) - @test ForwardDiff.jacobian(ev1, x0_static) ≈ Calculus.finite_difference_jacobian(ev1, x0) + + x0_svector = SVector{2}(x0) + @test ForwardDiff.jacobian(ev1, x0_svector) ≈ Calculus.finite_difference_jacobian(ev1, x0) + + x0_mvector = MVector{2}(x0) + @test ForwardDiff.jacobian(ev1, x0_mvector) ≈ Calculus.finite_difference_jacobian(ev1, x0) end @testset "type stability" begin diff --git a/test/MiscTest.jl b/test/MiscTest.jl index 2a22cb2f..96971560 100644 --- a/test/MiscTest.jl +++ b/test/MiscTest.jl @@ -6,6 +6,9 @@ using Test using ForwardDiff using DiffTests using SparseArrays: sparse +using StaticArrays +using IrrationalConstants +using LinearAlgebra include(joinpath(dirname(@__FILE__), "utils.jl")) @@ -141,6 +144,24 @@ let i, j end end +# AbstractIrrational numbers # +#----------------------------# +struct TestTag end +intrand(V) = V == Int ? rand(2:10) : rand(V) + +for N in (0,3), V in (Int, Float32), I in (Irrational, AbstractIrrational) + PARTIALS = ForwardDiff.Partials{N,V}(ntuple(n -> intrand(V), N)) + PRIMAL = intrand(V) + FDNUM = ForwardDiff.Dual{TestTag()}(PRIMAL, PARTIALS) + + @test promote_rule(typeof(FDNUM), I) == promote_rule(I, typeof(FDNUM)) + # π::Irrational, twoπ::AbstractIrrational + for IRR in (π, twoπ) + val_dual, val_irr = promote(FDNUM, IRR) + @test (val_irr, val_dual) == promote(IRR, FDNUM) + end +end + ######## # misc # ######## @@ -158,4 +179,41 @@ end @test ForwardDiff.derivative(x -> rem2pi(x, RoundUp), rand()) == 1 @test ForwardDiff.derivative(x -> rem2pi(x, RoundDown), rand()) == 1 +@testset "_lyap_div!!" begin + # In-place version for `Matrix` + A = rand(3, 3) + Acopy = copy(A) + λ = rand(3) + B = @inferred(ForwardDiff._lyap_div!!(A, λ)) + @test B === A + @test B[diagind(B)] == Acopy[diagind(Acopy)] + no_diag(X) = [X[i] for i in eachindex(X) if !(i in diagind(X))] + @test no_diag(B) == no_diag(Acopy ./ (λ' .- λ)) + + # Immutable static arrays + A_smatrix = SMatrix{3,3}(Acopy) + λ_svector = SVector{3}(λ) + B_smatrix = @inferred(ForwardDiff._lyap_div!!(A_smatrix, λ_svector)) + @test B_smatrix !== A_smatrix + @test B_smatrix isa SMatrix{3,3} + @test B_smatrix == B + λ_mvector = MVector{3}(λ) + B_smatrix = @inferred(ForwardDiff._lyap_div!!(A_smatrix, λ_mvector)) + @test B_smatrix !== A_smatrix + @test B_smatrix isa SMatrix{3,3} + @test B_smatrix == B + + # Mutable static arrays + A_mmatrix = MMatrix{3,3}(Acopy) + λ_svector = SVector{3}(λ) + B_mmatrix = @inferred(ForwardDiff._lyap_div!!(A_mmatrix, λ_svector)) + @test B_mmatrix === A_mmatrix + @test B_mmatrix == B + A_mmatrix = MMatrix{3,3}(Acopy) + λ_mvector = MVector{3}(λ) + B_mmatrix = @inferred(ForwardDiff._lyap_div!!(A_mmatrix, λ_mvector)) + @test B_mmatrix === A_mmatrix + @test B_mmatrix == B +end + end # module diff --git a/test/PartialsTest.jl b/test/PartialsTest.jl index 39fb05d7..8372a53e 100644 --- a/test/PartialsTest.jl +++ b/test/PartialsTest.jl @@ -7,8 +7,7 @@ using ForwardDiff: Partials samerng() = MersenneTwister(1) -for N in (0, 3), T in (Int, Float32, Float64) - println(" ...testing Partials{$N,$T}") +@testset "Partials{$N,$T}" for N in (0, 3), T in (Int, Float32, Float64) VALUES = (rand(T,N)...,) PARTIALS = Partials{N,T}(VALUES) @@ -62,7 +61,6 @@ for N in (0, 3), T in (Int, Float32, Float64) @test isequal(PARTIALS, copy(PARTIALS)) @test isequal(PARTIALS, PARTIALS2) == (N == 0) - @test hash(PARTIALS) == hash(VALUES, ForwardDiff.PARTIALS_HASH) @test hash(PARTIALS) == hash(copy(PARTIALS)) @test hash(PARTIALS, hash(1)) == hash(copy(PARTIALS), hash(1)) @test hash(PARTIALS, hash(1)) == hash(copy(PARTIALS), hash(1))