diff --git a/src/apiutils.jl b/src/apiutils.jl index 222e7bcb..f13c5d09 100644 --- a/src/apiutils.jl +++ b/src/apiutils.jl @@ -17,17 +17,19 @@ end ################################### # vector mode function evaluation # ################################### - -@generated function dualize(::Type{T}, x::SArray{S,V,D,N}) where {T,S,V,D,N} +@generated function dualize(::Type{T}, x::Union{FieldVector, SArray}) where {T} + N = length(x) + V = eltype(x) dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i}())) for i in 1:N]...) return quote - chunk = Chunk{N}() + chunk = Chunk{$N}() $(Expr(:meta, :inline)) - return SArray{S}($(dx)) + # return SArray{S}($(dx)) + similar_type($x,Dual{T,$V,$N})($(dx)) end end -@inline static_dual_eval(::Type{T}, f, x::SArray) where {T} = f(dualize(T, x)) +@inline static_dual_eval(::Type{T}, f, x::Union{FieldVector, SArray}) where {T} = f(dualize(T, x)) function vector_mode_dual_eval(f, x, cfg::Union{JacobianConfig,GradientConfig}) xdual = cfg.duals diff --git a/src/gradient.jl b/src/gradient.jl index 2cdb99bc..26cfe0b1 100644 --- a/src/gradient.jl +++ b/src/gradient.jl @@ -39,21 +39,23 @@ function gradient!(result::Union{AbstractArray,DiffResult}, f, x::AbstractArray, return result end -@inline gradient(f, x::SArray) = vector_mode_gradient(f, x) -@inline gradient(f, x::SArray, cfg::GradientConfig) = gradient(f, x) +@inline gradient(f, x::Union{FieldVector, SArray}) = vector_mode_gradient(f, x) +@inline gradient(f, x::Union{FieldVector, SArray}, cfg::GradientConfig) = gradient(f, x) -@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::SArray) = vector_mode_gradient!(result, f, x) -@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::SArray, cfg::GradientConfig) = gradient!(result, f, x) +@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::Union{FieldVector, SArray}) = vector_mode_gradient!(result, f, x) +@inline gradient!(result::Union{AbstractArray,DiffResult}, f, x::Union{FieldVector, SArray}, cfg::GradientConfig) = gradient!(result, f, x) ##################### # result extraction # ##################### -@generated function extract_gradient(::Type{T}, y::Real, ::SArray{S,X,D,N}) where {T,S,X,D,N} +@generated function extract_gradient(::Type{T}, y::Real, x::Union{FieldVector, SArray}) where {T} + N = length(x) + V = eltype(x) result = Expr(:tuple, [:(partials(T, y, $i)) for i in 1:N]...) return quote $(Expr(:meta, :inline)) - return SArray{S}($result) + return similar_type($x,$V)($result) end end @@ -102,16 +104,15 @@ function vector_mode_gradient!(result, f, x, cfg::GradientConfig{T}) where {T} return result end -@inline function vector_mode_gradient(f::F, x::SArray{S,V}) where {F,S,V} - T = typeof(Tag(f,V)) +@inline function vector_mode_gradient(f::F, x::Union{FieldVector, SArray}) where {F} + T = typeof(Tag(f,eltype(x))) return extract_gradient(T, static_dual_eval(T, f, x), x) end -@inline function vector_mode_gradient!(result, f::F, x::SArray{S,V}) where {F,S,V} - T = typeof(Tag(f,V)) +@inline function vector_mode_gradient!(result, f::F, x::Union{FieldVector, SArray}) where {F} + T = typeof(Tag(f,eltype(x))) return extract_gradient!(T, result, static_dual_eval(T, f, x)) end - ############## # chunk mode # ############## diff --git a/src/hessian.jl b/src/hessian.jl index 6bed53a7..01c0005e 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -55,17 +55,18 @@ function hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig{T} return result end -hessian(f, x::SArray) = jacobian(y -> gradient(f, y), x) +hessian(f, x::Union{FieldVector, SArray}) = jacobian(y -> gradient(f, y), x) -hessian(f, x::SArray, cfg::HessianConfig) = hessian(f, x) +hessian(f, x::Union{FieldVector, SArray}, cfg::HessianConfig) = hessian(f, x) -hessian!(result::AbstractArray, f, x::SArray) = jacobian!(result, y -> gradient(f, y), x) +hessian!(result::AbstractArray, f, x::Union{FieldVector, SArray}) = jacobian!(result, y -> gradient(f, y), x) -hessian!(result::MutableDiffResult, f, x::SArray) = hessian!(result, f, x, HessianConfig(f, result, x)) +hessian!(result::MutableDiffResult, f, x::Union{FieldVector, SArray}) = hessian!(result, f, x, HessianConfig(f, result, x)) -hessian!(result::ImmutableDiffResult, f, x::SArray, cfg::HessianConfig) = hessian!(result, f, x) +hessian!(result::ImmutableDiffResult, f, x::Union{FieldVector, SArray}, cfg::HessianConfig) = hessian!(result, f, x) -function hessian!(result::ImmutableDiffResult, f::F, x::SArray{S,V}) where {F,S,V} +function hessian!(result::ImmutableDiffResult, f::F, x::Union{FieldVector, SArray}) where {F} + V = eltype(x) T = typeof(Tag(f,V)) d1 = dualize(T, x) d2 = dualize(T, d1) diff --git a/src/jacobian.jl b/src/jacobian.jl index 7e8cf2fe..3f3be593 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -29,7 +29,7 @@ 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} - CHK && checktag(T, f!, x) + CHK && checktag(T, f!, x) if chunksize(cfg) == length(x) return vector_mode_jacobian(f!, y, x, cfg) else @@ -78,26 +78,29 @@ function jacobian!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray return result end -@inline jacobian(f, x::SArray) = vector_mode_jacobian(f, x) -@inline jacobian(f, x::SArray, cfg::JacobianConfig) = jacobian(f, x) +@inline jacobian(f, x::Union{FieldVector, SArray}) = vector_mode_jacobian(f, x) +@inline jacobian(f, x::Union{FieldVector, SArray}, cfg::JacobianConfig) = jacobian(f, x) -@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::SArray) = vector_mode_jacobian!(result, f, x) -@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::SArray, cfg::JacobianConfig) = jacobian!(result, f, x) +@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::Union{FieldVector, SArray}) = vector_mode_jacobian!(result, f, x) +@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::Union{FieldVector, SArray}, cfg::JacobianConfig) = jacobian!(result, f, x) ##################### # result extraction # ##################### -@generated function extract_jacobian(::Type{T}, ydual::SArray{SY,VY,DY,M}, - x::SArray{SX,VX,DX,N}) where {T,SY,VY,DY,M,SX,VX,DX,N} +@generated function extract_jacobian(::Type{T}, ydual::Union{FieldVector, SArray}, + x::Union{FieldVector, SArray}) where {T} + M = length(ydual) + N = length(x) result = Expr(:tuple, [:(partials(T, ydual[$i], $j)) for i in 1:M, j in 1:N]...) return quote $(Expr(:meta, :inline)) - return SArray{Tuple{M,N}}($result) + return SArray{Tuple{$M,$N}}($result) end end -function extract_jacobian(::Type{T}, ydual::AbstractArray, x::SArray{S,V,D,N}) where {T,S,V,D,N} +function extract_jacobian(::Type{T}, ydual::AbstractArray, x::Union{FieldVector, SArray}) where {T} + N = length(x) result = similar(ydual, valtype(eltype(ydual)), length(ydual), N) return extract_jacobian!(T, result, ydual, N) end @@ -165,12 +168,15 @@ function vector_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T,V,N}) return result end -@inline function vector_mode_jacobian(f::F, x::SArray{S,V}) where {F,S,V} +@inline function vector_mode_jacobian(f::F, x::Union{FieldVector, SArray}) where {F} + V = eltype(x) T = typeof(Tag(f,V)) return extract_jacobian(T, static_dual_eval(T, f, x), x) end -@inline function vector_mode_jacobian!(result, f::F, x::SArray{S,V,D,N}) where {F,S,V,D,N} +@inline function vector_mode_jacobian!(result, f::F, x::Union{FieldVector, SArray}) where {F} + N = length(x) + V = eltype(x) T = typeof(Tag(f,V)) ydual = static_dual_eval(T, f, x) result = extract_jacobian!(T, result, ydual, N) @@ -178,7 +184,8 @@ end return result end -@inline function vector_mode_jacobian!(result::ImmutableDiffResult, f::F, x::SArray{S,V,D,N}) where {F,S,V,D,N} +@inline function vector_mode_jacobian!(result::ImmutableDiffResult, f::F, x::Union{FieldVector, SArray}) where {F} + V = eltype(x) T = typeof(Tag(f,V)) ydual = static_dual_eval(T, f, x) result = DiffResults.jacobian!(result, extract_jacobian(T, ydual, x)) diff --git a/test/GradientTest.jl b/test/GradientTest.jl index c1f41d1a..c8e7a22f 100644 --- a/test/GradientTest.jl +++ b/test/GradientTest.jl @@ -136,5 +136,65 @@ sresult3 = ForwardDiff.gradient!(sresult3, prod, sx, scfg) @test DiffResults.gradient(sresult2) == DiffResults.gradient(result) @test DiffResults.gradient(sresult3) == DiffResults.gradient(result) +println(" ...testing specialized FieldVector codepaths") + +struct Point3D{R<:Real} <: FieldVector{3,R} + x::R + y::R + z::R +end +StaticArrays.similar_type(p::Type{P}, ::Type{R}, size::Size{(3,)}) where {P<:Point3D, R<:Real} = Point3D{R} + +x = rand(3, 1) +fx = Point3D(x) + +cfg = ForwardDiff.GradientConfig(nothing, x) +fcfg = ForwardDiff.GradientConfig(nothing, fx) + +actual = ForwardDiff.gradient(prod, x) +@test ForwardDiff.gradient(prod, fx) == actual[:] +@test ForwardDiff.gradient(prod, fx, cfg) == actual[:] +@test ForwardDiff.gradient(prod, fx, fcfg) == actual[:] + +out = similar(x) +ForwardDiff.gradient!(out, prod, fx) +@test out == actual + +out = similar(x) +ForwardDiff.gradient!(out, prod, fx, cfg) +@test out == actual + +out = similar(x) +ForwardDiff.gradient!(out, prod, fx, fcfg) +@test out == actual + +result = DiffResults.GradientResult(x) +result = ForwardDiff.gradient!(result, prod, x) + +result1 = DiffResults.GradientResult(x) +result2 = DiffResults.GradientResult(x) +result3 = DiffResults.GradientResult(x) +result1 = ForwardDiff.gradient!(result1, prod, fx) +result2 = ForwardDiff.gradient!(result2, prod, fx, cfg) +result3 = ForwardDiff.gradient!(result3, prod, fx, fcfg) +@test DiffResults.value(result1) == DiffResults.value(result) +@test DiffResults.value(result2) == DiffResults.value(result) +@test DiffResults.value(result3) == DiffResults.value(result) +@test DiffResults.gradient(result1) == DiffResults.gradient(result) +@test DiffResults.gradient(result2) == DiffResults.gradient(result) +@test DiffResults.gradient(result3) == DiffResults.gradient(result) + +fresult1 = DiffResults.GradientResult(fx) +fresult2 = DiffResults.GradientResult(fx) +fresult3 = DiffResults.GradientResult(fx) +fresult1 = ForwardDiff.gradient!(fresult1, prod, fx) +fresult2 = ForwardDiff.gradient!(fresult2, prod, fx, cfg) +fresult3 = ForwardDiff.gradient!(fresult3, prod, fx, fcfg) +@test DiffResults.value(fresult1) == DiffResults.value(result) +@test DiffResults.value(fresult2) == DiffResults.value(result) +@test DiffResults.value(fresult3) == DiffResults.value(result) +@test DiffResults.gradient(fresult1) == DiffResults.gradient(result)[:] +@test DiffResults.gradient(fresult2) == DiffResults.gradient(result)[:] +@test DiffResults.gradient(fresult3) == DiffResults.gradient(result)[:] end # module diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 256d3cab..bd178fba 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -153,4 +153,72 @@ sresult3 = ForwardDiff.hessian!(sresult3, prod, sx, ForwardDiff.HessianConfig(pr @test DiffResults.hessian(sresult2) == DiffResults.hessian(result) @test DiffResults.hessian(sresult3) == DiffResults.hessian(result) + +println(" ...testing specialized FieldVector codepaths") + +struct Point3D{R<:Real} <: FieldVector{3,R} + x::R + y::R + z::R +end +StaticArrays.similar_type(p::Type{P}, ::Type{R}, size::Size{(3,)}) where {P<:Point3D, R<:Real} = Point3D{R} + +x = rand(3, 1) +fx = Point3D(x) + +cfg = ForwardDiff.HessianConfig(nothing, x) +fcfg = ForwardDiff.HessianConfig(nothing, fx) + +actual = ForwardDiff.hessian(prod, x) +@test ForwardDiff.hessian(prod, fx) == actual +@test ForwardDiff.hessian(prod, fx, cfg) == actual +@test ForwardDiff.hessian(prod, fx, fcfg) == actual + +out = similar(x, length(x), length(x)) +ForwardDiff.hessian!(out, prod, fx) +@test out == actual + +out = similar(x, length(x), length(x)) +ForwardDiff.hessian!(out, prod, fx, cfg) +@test out == actual + +out = similar(x, length(x), length(x)) +ForwardDiff.hessian!(out, prod, fx, fcfg) +@test out == actual + +result = DiffResults.HessianResult(x) +result = ForwardDiff.hessian!(result, prod, x) + +result1 = DiffResults.HessianResult(x) +result2 = DiffResults.HessianResult(x) +result3 = DiffResults.HessianResult(x) +result1 = ForwardDiff.hessian!(result1, prod, fx) +result2 = ForwardDiff.hessian!(result2, prod, fx, ForwardDiff.HessianConfig(prod, result2, x, ForwardDiff.Chunk(x), nothing)) +result3 = ForwardDiff.hessian!(result3, prod, fx, ForwardDiff.HessianConfig(prod, result3, x, ForwardDiff.Chunk(x), nothing)) +@test DiffResults.value(result1) == DiffResults.value(result) +@test DiffResults.value(result2) == DiffResults.value(result) +@test DiffResults.value(result3) == DiffResults.value(result) +@test DiffResults.gradient(result1) == DiffResults.gradient(result) +@test DiffResults.gradient(result2) == DiffResults.gradient(result) +@test DiffResults.gradient(result3) == DiffResults.gradient(result) +@test DiffResults.hessian(result1) == DiffResults.hessian(result) +@test DiffResults.hessian(result2) == DiffResults.hessian(result) +@test DiffResults.hessian(result3) == DiffResults.hessian(result) + +fresult1 = DiffResults.HessianResult(fx) +fresult2 = DiffResults.HessianResult(fx) +fresult3 = DiffResults.HessianResult(fx) +fresult1 = ForwardDiff.hessian!(fresult1, prod, fx) +fresult2 = ForwardDiff.hessian!(fresult2, prod, fx, ForwardDiff.HessianConfig(prod, fresult2, x, ForwardDiff.Chunk(x), nothing)) +fresult3 = ForwardDiff.hessian!(fresult3, prod, fx, ForwardDiff.HessianConfig(prod, fresult3, x, ForwardDiff.Chunk(x), nothing)) +@test DiffResults.value(fresult1) == DiffResults.value(result) +@test DiffResults.value(fresult2) == DiffResults.value(result) +@test DiffResults.value(fresult3) == DiffResults.value(result) +@test DiffResults.gradient(fresult1) == DiffResults.gradient(result)[:] +@test DiffResults.gradient(fresult2) == DiffResults.gradient(result)[:] +@test DiffResults.gradient(fresult3) == DiffResults.gradient(result)[:] +@test DiffResults.hessian(fresult1) == DiffResults.hessian(result) +@test DiffResults.hessian(fresult2) == DiffResults.hessian(result) +@test DiffResults.hessian(fresult3) == DiffResults.hessian(result) + end # module diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index db204f25..69004fce 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -220,4 +220,66 @@ sresult3 = ForwardDiff.jacobian!(sresult3, diff, sx, scfg) @test DiffResults.jacobian(sresult2) == DiffResults.jacobian(result) @test DiffResults.jacobian(sresult3) == DiffResults.jacobian(result) + +println(" ...testing specialized FieldVector codepaths") + +struct Point3D{R<:Real} <: FieldVector{3,R} + x::R + y::R + z::R +end +StaticArrays.similar_type(p::Type{P}, ::Type{R}, size::Size{(3,)}) where {P<:Point3D, R<:Real} = Point3D{R} + +x = rand(3, 1) +fx = Point3D(x) + +cfg = ForwardDiff.JacobianConfig(nothing, x) +fcfg = ForwardDiff.JacobianConfig(nothing, fx) + +actual = ForwardDiff.jacobian(diff, x) +@test ForwardDiff.jacobian(diff, fx) == actual +@test ForwardDiff.jacobian(diff, fx, cfg) == actual +@test ForwardDiff.jacobian(diff, fx, fcfg) == actual + +out = similar(x, 2, 3) +ForwardDiff.jacobian!(out, diff, fx) +@test out == actual + +out = similar(x, 2, 3) +ForwardDiff.jacobian!(out, diff, fx, cfg) +@test out == actual + +out = similar(x, 2, 3) +ForwardDiff.jacobian!(out, diff, fx, fcfg) +@test out == actual + +result = DiffResults.JacobianResult(similar(x, 2), x) +result = ForwardDiff.jacobian!(result, diff, x) + +result1 = DiffResults.JacobianResult(similar(fx, 2), fx) +result2 = DiffResults.JacobianResult(similar(fx, 2), fx) +result3 = DiffResults.JacobianResult(similar(fx, 2), fx) +result1 = ForwardDiff.jacobian!(result1, diff, fx) +result2 = ForwardDiff.jacobian!(result2, diff, fx, cfg) +result3 = ForwardDiff.jacobian!(result3, diff, fx, fcfg) +@test DiffResults.value(result1) == DiffResults.value(result) +@test DiffResults.value(result2) == DiffResults.value(result) +@test DiffResults.value(result3) == DiffResults.value(result) +@test DiffResults.jacobian(result1) == DiffResults.jacobian(result) +@test DiffResults.jacobian(result2) == DiffResults.jacobian(result) +@test DiffResults.jacobian(result3) == DiffResults.jacobian(result) + +sy = @SVector fill(zero(eltype(fx)), 2) +fresult1 = DiffResults.JacobianResult(sy, fx) +fresult2 = DiffResults.JacobianResult(sy, fx) +fresult3 = DiffResults.JacobianResult(sy, fx) +fresult1 = ForwardDiff.jacobian!(fresult1, diff, fx) +fresult2 = ForwardDiff.jacobian!(fresult2, diff, fx, cfg) +fresult3 = ForwardDiff.jacobian!(fresult3, diff, fx, fcfg) +@test DiffResults.value(fresult1) == DiffResults.value(result) +@test DiffResults.value(fresult2) == DiffResults.value(result) +@test DiffResults.value(fresult3) == DiffResults.value(result) +@test DiffResults.jacobian(fresult1) == DiffResults.jacobian(result) +@test DiffResults.jacobian(fresult2) == DiffResults.jacobian(result) +@test DiffResults.jacobian(fresult3) == DiffResults.jacobian(result) end # module diff --git a/test/MiscTest.jl b/test/MiscTest.jl index 501f6b4e..13a94099 100644 --- a/test/MiscTest.jl +++ b/test/MiscTest.jl @@ -6,6 +6,7 @@ using Compat using Compat.Test using ForwardDiff using DiffTests +using StaticArrays include(joinpath(dirname(@__FILE__), "utils.jl")) @@ -137,6 +138,58 @@ for i in 1:3, j in 1:3 i != j && (@test h[i, j] ≈ 0.0) end +################ +# FieldVectors # +################ + +# issue 305 + +# Test mutable FieldVector + +mutable struct MPoint{R<:Real} <: FieldVector{2,R} + x::R + y::R +end + +StaticArrays.similar_type(p::Type{P}, ::Type{R}, size::Size{(2,)}) where {P<:MPoint, R<:Real} = MPoint{R} + +f305scalar(x::MPoint) = sum(x+x) +@test ForwardDiff.gradient(f305scalar, MPoint(1.0, 2.0)) == MPoint(2.0, 2.0) +@test ForwardDiff.hessian(f305scalar, MPoint(1.0, 2.0)) == zeros(2,2) + +f305vector(x::MPoint) = x+x +@test ForwardDiff.jacobian(f305vector, MPoint(1.0, 2.0)) == [2.0 0.0; 0.0 2.0] + +# Test correct dispatch in differentiation + +struct DifferentPoint{R<:Real} <: FieldVector{2,R} + x::R + y::R +end + +StaticArrays.similar_type(p::Type{P}, ::Type{R}, size::Size{(2,)}) where {P<:DifferentPoint, R<:Real} = DifferentPoint{R} + +f305dispatch(p) = sum(p.^2) +f305dispatch(p::MPoint) = p.x^2 +f305dispatch(p::DifferentPoint) = p.y^2 + +@test ForwardDiff.gradient(f305dispatch, [1.0, 2.0]) == [2.0, 4.0] +@test ForwardDiff.gradient(f305dispatch, MPoint(1.0, 2.0)) == [2.0, 0.0] +@test ForwardDiff.gradient(f305dispatch, DifferentPoint(1.0, 2.0)) == [0.0, 4.0] + +@test ForwardDiff.hessian(f305dispatch, [1.0, 2.0]) == 2.0*eye(2) +@test ForwardDiff.hessian(f305dispatch, MPoint(1.0, 2.0)) == diagm([2.0, 0.0]) # [2.0 0.0; 0.0 0.0] +@test ForwardDiff.hessian(f305dispatch, DifferentPoint(1.0, 2.0)) == diagm([0.0, 2.0]) # [0.0 0.0; 0.0 2.0] + +f305dispatchvec(p) = p.^2 +f305dispatchvec(p::MPoint) = MPoint(p.x^2, one(eltype(p))) +f305dispatchvec(p::DifferentPoint) = DifferentPoint(p.y^2, one(eltype(p))) + +@test ForwardDiff.jacobian(f305dispatchvec, [1.0, 2.0]) == diagm([2.0, 4.0]) +@test ForwardDiff.jacobian(f305dispatchvec, MPoint(1.0, 2.0)) == diagm([2.0, 0.0]) +@test ForwardDiff.jacobian(f305dispatchvec, DifferentPoint(1.0, 2.0)) == [0.0 4.0; 0.0 0.0] + + ######## # misc # ########