Merge branch 'master' into ck/aqua

.github/workflows/Invalidations.yml

@@ -16,15 +16,15 @@ jobs:
     if: github.base_ref == github.event.repository.default_branch
     runs-on: ubuntu-latest
     steps:
-    - uses: julia-actions/setup-julia@v1
+    - uses: julia-actions/setup-julia@v2
       with:
         version: '1'
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - uses: julia-actions/julia-buildpkg@v1
     - uses: julia-actions/julia-invalidations@v1
       id: invs_pr
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         ref: ${{ github.event.repository.default_branch }}
     - uses: julia-actions/julia-buildpkg@v1

.github/workflows/ci.yml

@@ -19,18 +19,18 @@ jobs:
         version:
           - '1.6'
           - '1'
-          - 'nightly'
+          # - 'nightly'
         os:
           - ubuntu-latest
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v3
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:
@@ -43,15 +43,15 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
   docs:
     name: Documentation
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
       - run: |

.gitignore

@@ -4,5 +4,6 @@
 *.DS_Store
 /docs/build/
 /docs/site/
+/docs/Manifest.toml
 /benchmark_data/
 /Manifest.toml

Project.toml

@@ -1,6 +1,6 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.11-DEV"
+version = "0.11.0-DEV"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"
@@ -30,6 +30,7 @@ DiffRules = "1.4"
 DiffTests = "0.1"
 InteractiveUtils = "1"
 LinearAlgebra = "1"
+IrrationalConstants = "0.1, 0.2"
 LogExpFunctions = "0.3"
 NaNMath = "1"
 Preferences = "1"
@@ -46,9 +47,10 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 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 = ["Aqua", "Calculus", "DiffTests", "SparseArrays", "StaticArrays", "Test", "InteractiveUtils"]
+test = ["Aqua", "Calculus", "DiffTests", "IrrationalConstants", "SparseArrays", "StaticArrays", "Test", "InteractiveUtils"]

docs/Project.toml

@@ -1,3 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+
+[compat]
+Documenter = "1"

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"

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,19 +23,17 @@ 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*_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)

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 #
 ##############################
@@ -384,17 +388,31 @@ end
 # Before PR#481 this loop ran over this list:
 # BINARY_PREDICATES = Symbol[:isequal, :isless, :<, :>, :(==), :(!=), :(<=), :(>=)]
 # Not a minimal set, as Base defines some in terms of others.
-for pred in [:isless, :<, :>, :(<=), :(>=)]
+for pred in [:<, :>]
+    predeq = Symbol(pred, :(=))
     @eval begin
         @define_binary_dual_op(
             Base.$(pred),
-            $(pred)(value(x), value(y)),
-            $(pred)(value(x), y),
-            $(pred)(x, value(y)),
+            $(pred)(value(x), value(y)) || (value(x) == value(y) && $(pred)(partials(x), partials(y))),
+            $(pred)(value(x), y) || (value(x) == y && $(pred)(partials(x), zero(partials(x)))),
+            $(pred)(x, value(y)) || (x == value(y) && $(pred)(zero(partials(y)), partials(y))),
+        )
+        @define_binary_dual_op(
+            Base.$(predeq),
+            $(pred)(value(x), value(y)) || (value(x) == value(y) && $(predeq)(partials(x), partials(y))),
+            $(pred)(value(x), y) || (value(x) == y && $(predeq)(partials(x), zero(partials(x)))),
+            $(pred)(x, value(y)) || (x == value(y) && $(predeq)(zero(partials(y)), partials(y))),
         )
     end
 end
 
+@define_binary_dual_op(
+    Base.isless,
+    isless(value(x), value(y)) || (isequal(value(x), value(y)) && isless(partials(x), partials(y))),
+    isless(value(x), y)        || (isequal(value(x), y) && isless(partials(x), zero(partials(x)))),
+    isless(x, value(y))        || (isequal(x, value(y)) && isless(zero(partials(y)), partials(y))),
+)
+
 Base.iszero(x::Dual) = iszero(value(x)) && iszero(partials(x))  # shortcut, equivalent to x == zero(x)
 
 for pred in [:isequal, :(==)]
@@ -434,7 +452,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}
@@ -449,6 +467,7 @@ for R in (Irrational, Real, BigFloat, Bool)
 end
 
 @inline Base.convert(::Type{Dual{T,V,N}}, d::Dual{T}) where {T,V,N} = Dual{T}(V(value(d)), convert(Partials{N,V}, partials(d)))
+@inline Base.convert(::Type{Dual{T,Dual{T,V,M},N}}, d::Dual{T,V,M}) where {T,V,N,M} = Dual{T}(d, Partials{N,Dual{T,V,M}}(zero_tuple(NTuple{N,Dual{T,V,M}})))
 @inline Base.convert(::Type{Dual{T,V,N}}, x) where {T,V,N} = Dual{T}(V(x), zero(Partials{N,V}))
 @inline Base.convert(::Type{Dual{T,V,N}}, x::Number) where {T,V,N} = Dual{T}(V(x), zero(Partials{N,V}))
 Base.convert(::Type{D}, d::D) where {D<:Dual} = d
@@ -719,7 +738,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...))
@@ -737,8 +760,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] /= μ - λ
@@ -747,17 +781,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
 

src/prelude.jl

@@ -5,16 +5,11 @@ const AMBIGUOUS_TYPES = (AbstractFloat, Irrational, Integer, Rational, Real, Rou
 
 const UNARY_PREDICATES = Symbol[:isinf, :isnan, :isfinite, :iseven, :isodd, :isreal, :isinteger]
 
-const DEFAULT_CHUNK_THRESHOLD = 12
-
 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)

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

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)