Skip to content

Commit

Permalink
Merge pull request #712 from JuliaDiff/dw/backport
Browse files Browse the repository at this point in the history
Backport bug and test fixes to release-0.10 to prepare for new 0.10 release
  • Loading branch information
devmotion authored Nov 1, 2024
2 parents 2ff6808 + c7f62ca commit 1674ee9
Show file tree
Hide file tree
Showing 22 changed files with 203 additions and 99 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ jobs:
- '1.0'
- '1.6'
- '1'
- 'nightly'
# - 'nightly'
os:
- ubuntu-latest
arch:
Expand Down
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -4,5 +4,6 @@
*.DS_Store
/docs/build/
/docs/site/
/docs/Manifest.toml
/benchmark_data/
/Manifest.toml
6 changes: 4 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
3 changes: 3 additions & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
[deps]
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

[compat]
Documenter = "1"
3 changes: 2 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
56 changes: 27 additions & 29 deletions ext/ForwardDiffStaticArraysExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)]...)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
7 changes: 6 additions & 1 deletion src/ForwardDiff.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
11 changes: 7 additions & 4 deletions src/derivative.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
38 changes: 30 additions & 8 deletions src/dual.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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 #
##############################
Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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...))
Expand All @@ -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] /= μ - λ
Expand All @@ -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

Expand Down
4 changes: 3 additions & 1 deletion src/gradient.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down
8 changes: 5 additions & 3 deletions src/hessian.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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}())
Expand All @@ -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}())
Expand Down Expand Up @@ -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}())
Expand Down
Loading

0 comments on commit 1674ee9

Please sign in to comment.