From 3106983e27a7dd5e7049c861d16eb36772558dbd Mon Sep 17 00:00:00 2001 From: Jarrett Revels Date: Thu, 13 Apr 2017 14:28:19 -0400 Subject: [PATCH 1/4] add specialized StaticArrays-based methods --- REQUIRE | 1 + src/ForwardDiff.jl | 24 ++++++++++++++++++++++++ src/config.jl | 27 --------------------------- src/derivative.jl | 2 +- src/dual.jl | 2 +- src/gradient.jl | 16 ++++++++++++++++ src/jacobian.jl | 26 ++++++++++++++++++++++++++ src/partials.jl | 2 +- src/utils.jl | 15 ++++++++++++++- 9 files changed, 84 insertions(+), 31 deletions(-) diff --git a/REQUIRE b/REQUIRE index 7f2f5e18..5260c389 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,4 +1,5 @@ julia 0.6- +StaticArrays 0.4 DiffBase 0.0.3 Calculus 0.2.0 NaNMath 0.2.2 diff --git a/src/ForwardDiff.jl b/src/ForwardDiff.jl index 429c6291..fc86116f 100644 --- a/src/ForwardDiff.jl +++ b/src/ForwardDiff.jl @@ -4,6 +4,7 @@ module ForwardDiff using DiffBase using DiffBase: DiffResult +using StaticArrays import Calculus import NaNMath @@ -41,6 +42,29 @@ const REAL_TYPES = (AbstractFloat, Irrational, Integer, Rational, Real) const DEFAULT_CHUNK_THRESHOLD = 10 +struct Chunk{N} end + +function Chunk(input_length::Integer, threshold::Integer = DEFAULT_CHUNK_THRESHOLD) + N = pickchunksize(input_length, threshold) + return Chunk{N}() +end + +function Chunk(x::AbstractArray, threshold::Integer = DEFAULT_CHUNK_THRESHOLD) + return Chunk(length(x), threshold) +end + +# Constrained to `N <= threshold`, minimize (in order of priority): +# 1. the number of chunks that need to be computed +# 2. the number of "left over" perturbations in the final chunk +function pickchunksize(input_length, threshold = DEFAULT_CHUNK_THRESHOLD) + if input_length <= threshold + return input_length + else + nchunks = round(Int, input_length / DEFAULT_CHUNK_THRESHOLD, RoundUp) + return round(Int, input_length / nchunks, RoundUp) + end +end + ############ # includes # ############ diff --git a/src/config.jl b/src/config.jl index c43077aa..eedb8b66 100644 --- a/src/config.jl +++ b/src/config.jl @@ -18,33 +18,6 @@ struct Tag{F,H} end end end -######### -# Chunk # -######### - -struct Chunk{N} end - -function Chunk(input_length::Integer, threshold::Integer = DEFAULT_CHUNK_THRESHOLD) - N = pickchunksize(input_length, threshold) - return Chunk{N}() -end - -function Chunk(x::AbstractArray, threshold::Integer = DEFAULT_CHUNK_THRESHOLD) - return Chunk(length(x), threshold) -end - -# Constrained to `N <= threshold`, minimize (in order of priority): -# 1. the number of chunks that need to be computed -# 2. the number of "left over" perturbations in the final chunk -function pickchunksize(input_length, threshold = DEFAULT_CHUNK_THRESHOLD) - if input_length <= threshold - return input_length - else - nchunks = round(Int, input_length / DEFAULT_CHUNK_THRESHOLD, RoundUp) - return round(Int, input_length / nchunks, RoundUp) - end -end - ################## # AbstractConfig # ################## diff --git a/src/derivative.jl b/src/derivative.jl index 4d8449a5..9c04171a 100644 --- a/src/derivative.jl +++ b/src/derivative.jl @@ -22,7 +22,7 @@ end @inline extract_derivative(y::Dual{T,V,1}) where {T,V} = partials(y, 1) @inline extract_derivative(y::Real) = zero(y) -@inline extract_derivative(y::AbstractArray) = extract_derivative!(similar(y, valtype(eltype(y))), y) +@inline extract_derivative(y::AbstractArray) = map(extract_derivative, y) # mutating # #----------# diff --git a/src/dual.jl b/src/dual.jl index 0c49347b..af01c9d3 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -21,7 +21,7 @@ end @inline (::Type{Dual{T}})(value::Real, partials::Tuple) where {T} = Dual{T}(value, Partials(partials)) @inline (::Type{Dual{T}})(value::Real, partials::Tuple{}) where {T} = Dual{T}(value, Partials{0,typeof(value)}(partials)) @inline (::Type{Dual{T}})(value::Real, partials::Real...) where {T} = Dual{T}(value, partials) -@inline (::Type{Dual{T}})(value::V, ::Type{Val{N}}, ::Type{Val{i}}) where {T,V<:Real,N,i} = Dual{T}(value, single_seed(Partials{N,V}, Val{i})) +@inline (::Type{Dual{T}})(value::V, ::Chunk{N}, p::Val{i}) where {T,V<:Real,N,i} = Dual{T}(value, single_seed(Partials{N,V}, p)) @inline Dual(args...) = Dual{Void}(args...) diff --git a/src/gradient.jl b/src/gradient.jl index c611d22a..6368fa08 100644 --- a/src/gradient.jl +++ b/src/gradient.jl @@ -24,10 +24,26 @@ function gradient!(out, f::F, x, cfg::AllowedGradientConfig{F,H} = GradientConfi return out end +@inline function gradient(f::F, x::SArray) where F + return extract_gradient(vector_mode_dual_eval(f, x), x) +end + +@inline function gradient!(out, f::F, x::SArray) where F + return extract_gradient!(out, vector_mode_dual_eval(f, x)) +end + ##################### # result extraction # ##################### +@generated function extract_gradient(y::Real, ::SArray{S,X,D,N}) where {S,X,D,N} + result = Expr(:tuple, [:(partials(y, $i)) for i in 1:N]...) + return quote + $(Expr(:meta, :inline)) + return SArray{S}($result) + end +end + function extract_gradient!(out::DiffResult, y::Real) DiffBase.value!(out, y) grad = DiffBase.gradient(out) diff --git a/src/jacobian.jl b/src/jacobian.jl index 96f73de4..b62ab669 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -43,10 +43,36 @@ function jacobian!(out, f!::F, y, x, cfg::AllowedJacobianConfig{F,H} = JacobianC return out end +@inline function jacobian(f::F, x::SArray) where F + return extract_jacobian(vector_mode_dual_eval(f, x), x) +end + +@inline function jacobian!(out, f::F, x::SArray) where F + return extract_jacobian!(out, vector_mode_dual_eval(f, x)) +end + ##################### # result extraction # ##################### +@generated function extract_jacobian(ydual::SArray{SY,VY,DY,M}, + x::SArray{SX,VX,DX,N}) where {SY,VY,DY,M,SX,VX,DX,N} + elements = Vector{Any}(0) + for i in 1:M, j in 1:N + push!(elements, :(partials(ydual[$i], $j))) + end + result = Expr(:tuple, elements...) + return quote + $(Expr(:meta, :inline)) + return SArray{Tuple{M,N}}($result) + end +end + +function extract_jacobian(ydual::AbstractArray, x::SArray{S,V,D,N}) where {S,V,D,N} + out = similar(ydual, valtype(eltype(ydual)), length(ydual), N) + return extract_jacobian!(out, ydual, N) +end + function extract_jacobian!(out::AbstractArray, ydual::AbstractArray, n) out_reshaped = reshape(out, length(ydual), n) for col in 1:size(out_reshaped, 2), row in 1:size(out_reshaped, 1) diff --git a/src/partials.jl b/src/partials.jl index 1485364a..4c0deef0 100644 --- a/src/partials.jl +++ b/src/partials.jl @@ -6,7 +6,7 @@ end # Utility/Accessor Functions # ############################## -@generated function single_seed(::Type{Partials{N,V}}, ::Type{Val{i}}) where {N,V,i} +@generated function single_seed(::Type{Partials{N,V}}, ::Val{i}) where {N,V,i} ex = Expr(:tuple, [ifelse(i === j, :(one(V)), :(zero(V))) for j in 1:N]...) return :(Partials($(ex))) end diff --git a/src/utils.jl b/src/utils.jl index a9d3dbd4..62124e7b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -17,6 +17,19 @@ end # vector mode function evaluation # ################################### +@generated function dualize(::F, x::SArray{S,V,D,N}) where {F,S,V,D,N} + tag = Tag(F, x) + dx = Expr(:tuple, [:(Dual{T}(x[$i], chunk, Val{$i}())) for i in 1:N]...) + return quote + chunk = Chunk{N}() + T = typeof($tag) + $(Expr(:meta, :inline)) + return SArray{S}($(dx)) + end +end + +@inline vector_mode_dual_eval(f::F, x::SArray) where {F} = f(dualize(f, x)) + function vector_mode_dual_eval(f::F, x, cfg::Union{JacobianConfig,GradientConfig}) where F xdual = cfg.duals seed!(xdual, x, cfg.seeds) @@ -36,7 +49,7 @@ end ################################## @generated function construct_seeds(::Type{Partials{N,V}}) where {N,V} - return Expr(:tuple, [:(single_seed(Partials{N,V}, Val{$i})) for i in 1:N]...) + return Expr(:tuple, [:(single_seed(Partials{N,V}, Val{$i}())) for i in 1:N]...) end function seed!(duals::AbstractArray{Dual{T,V,N}}, x, From 679604ed5334e1988840a54eb687d8472ff96706 Mon Sep 17 00:00:00 2001 From: Jarrett Revels Date: Tue, 18 Apr 2017 11:09:45 -0400 Subject: [PATCH 2/4] add StaticArray specializations for Hessians --- src/hessian.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/hessian.jl b/src/hessian.jl index c74cf3a8..77e6b204 100644 --- a/src/hessian.jl +++ b/src/hessian.jl @@ -29,3 +29,9 @@ function hessian!(out::DiffResult, f::F, x, cfg::AllowedHessianConfig{F,H} = Hes jacobian!(DiffBase.hessian(out), ∇f!, DiffBase.gradient(out), x, cfg.jacobian_config) return out end + +hessian(f::F, x::SArray) where {F} = jacobian(y -> gradient(f, y), x) + +hessian!(out, f::F, x::SArray) where {F} = jacobian!(out, y -> gradient(f, y), x) + +hessian!(out::DiffResult, f::F, x::SArray) where {F} = hessian!(out, f, x, HessianConfig(f, out, x)) From 0454e02197ed92a1ac26b9fe99850c13207053ca Mon Sep 17 00:00:00 2001 From: Jarrett Revels Date: Tue, 18 Apr 2017 13:11:14 -0400 Subject: [PATCH 3/4] add tests/fix bugs --- REQUIRE | 2 +- src/gradient.jl | 16 ++++++++++------ src/jacobian.jl | 25 ++++++++++++++----------- test/GradientTest.jl | 25 +++++++++++++++++++++++++ test/HessianTest.jl | 26 ++++++++++++++++++++++++++ test/JacobianTest.jl | 25 +++++++++++++++++++++++++ 6 files changed, 101 insertions(+), 18 deletions(-) diff --git a/REQUIRE b/REQUIRE index 5260c389..db757995 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,5 +1,5 @@ julia 0.6- -StaticArrays 0.4 +StaticArrays 0.5.0 DiffBase 0.0.3 Calculus 0.2.0 NaNMath 0.2.2 diff --git a/src/gradient.jl b/src/gradient.jl index 6368fa08..2e708c09 100644 --- a/src/gradient.jl +++ b/src/gradient.jl @@ -24,13 +24,9 @@ function gradient!(out, f::F, x, cfg::AllowedGradientConfig{F,H} = GradientConfi return out end -@inline function gradient(f::F, x::SArray) where F - return extract_gradient(vector_mode_dual_eval(f, x), x) -end +@inline gradient(f::F, x::SArray) where {F} = vector_mode_gradient(f, x) -@inline function gradient!(out, f::F, x::SArray) where F - return extract_gradient!(out, vector_mode_dual_eval(f, x)) -end +@inline gradient!(out, f::F, x::SArray) where {F} = vector_mode_gradient!(out, f, x) ##################### # result extraction # @@ -89,6 +85,14 @@ function vector_mode_gradient!(out, f::F, x, cfg) where {F} return out end +@inline function vector_mode_gradient(f::F, x::SArray) where F + return extract_gradient(vector_mode_dual_eval(f, x), x) +end + +@inline function vector_mode_gradient!(out, f::F, x::SArray) where F + return extract_gradient!(out, vector_mode_dual_eval(f, x)) +end + ############## # chunk mode # ############## diff --git a/src/jacobian.jl b/src/jacobian.jl index b62ab669..9218554d 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -43,13 +43,9 @@ function jacobian!(out, f!::F, y, x, cfg::AllowedJacobianConfig{F,H} = JacobianC return out end -@inline function jacobian(f::F, x::SArray) where F - return extract_jacobian(vector_mode_dual_eval(f, x), x) -end +@inline jacobian(f::F, x::SArray) where {F} = vector_mode_jacobian(f, x) -@inline function jacobian!(out, f::F, x::SArray) where F - return extract_jacobian!(out, vector_mode_dual_eval(f, x)) -end +@inline jacobian!(out, f::F, x::SArray) where {F} = vector_mode_jacobian!(out, f, x) ##################### # result extraction # @@ -57,11 +53,7 @@ end @generated function extract_jacobian(ydual::SArray{SY,VY,DY,M}, x::SArray{SX,VX,DX,N}) where {SY,VY,DY,M,SX,VX,DX,N} - elements = Vector{Any}(0) - for i in 1:M, j in 1:N - push!(elements, :(partials(ydual[$i], $j))) - end - result = Expr(:tuple, elements...) + result = Expr(:tuple, [:(partials(ydual[$i], $j)) for i in 1:M, j in 1:N]...) return quote $(Expr(:meta, :inline)) return SArray{Tuple{M,N}}($result) @@ -136,6 +128,17 @@ function vector_mode_jacobian!(out, f!::F, y, x, cfg::JacobianConfig{T,V,N}) whe return out end +@inline function vector_mode_jacobian(f::F, x::SArray) where F + return extract_jacobian(vector_mode_dual_eval(f, x), x) +end + +@inline function vector_mode_jacobian!(out, f::F, x::SArray{S,V,D,N}) where {F,S,V,D,N} + ydual = vector_mode_dual_eval(f, x) + extract_jacobian!(out, ydual, N) + extract_value!(out, ydual) + return out +end + # chunk mode # #------------# diff --git a/test/GradientTest.jl b/test/GradientTest.jl index 487ad6b2..6e4ae6c3 100644 --- a/test/GradientTest.jl +++ b/test/GradientTest.jl @@ -4,6 +4,7 @@ import Calculus using Base.Test using ForwardDiff +using StaticArrays include(joinpath(dirname(@__FILE__), "utils.jl")) @@ -67,4 +68,28 @@ for f in DiffBase.VECTOR_TO_NUMBER_FUNCS end end +########################################## +# test specialized StaticArray codepaths # +########################################## + +x = rand(3, 3) +sx = StaticArrays.SArray{Tuple{3,3}}(x) +out = similar(x) +actual = ForwardDiff.gradient(prod, x) + +@test ForwardDiff.gradient(prod, sx) == actual + +ForwardDiff.gradient!(out, prod, sx) + +@test out == actual + +result = DiffBase.GradientResult(x) +sresult = DiffBase.GradientResult(sx) + +ForwardDiff.gradient!(result, prod, x) +ForwardDiff.gradient!(sresult, prod, sx) + +@test DiffBase.value(sresult) == DiffBase.value(result) +@test DiffBase.gradient(sresult) == DiffBase.gradient(result) + end # module diff --git a/test/HessianTest.jl b/test/HessianTest.jl index 5d697b16..88fdac67 100644 --- a/test/HessianTest.jl +++ b/test/HessianTest.jl @@ -4,6 +4,7 @@ import Calculus using Base.Test using ForwardDiff +using StaticArrays include(joinpath(dirname(@__FILE__), "utils.jl")) @@ -78,4 +79,29 @@ for f in DiffBase.VECTOR_TO_NUMBER_FUNCS end end +########################################## +# test specialized StaticArray codepaths # +########################################## + +x = rand(3, 3) +sx = StaticArrays.SArray{Tuple{3,3}}(x) +out = similar(x, 9, 9) +actual = ForwardDiff.hessian(prod, x) + +@test ForwardDiff.hessian(prod, sx) == actual + +ForwardDiff.hessian!(out, prod, sx) + +@test out == actual + +result = DiffBase.HessianResult(x) +sresult = DiffBase.HessianResult(sx) + +ForwardDiff.hessian!(result, prod, x) +ForwardDiff.hessian!(sresult, prod, sx) + +@test DiffBase.value(sresult) == DiffBase.value(result) +@test DiffBase.gradient(sresult) == DiffBase.gradient(result) +@test DiffBase.hessian(sresult) == DiffBase.hessian(result) + end # module diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 83ee0a47..e83547c8 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -5,6 +5,7 @@ import Calculus using Base.Test using ForwardDiff using ForwardDiff: JacobianConfig +using StaticArrays include(joinpath(dirname(@__FILE__), "utils.jl")) @@ -146,4 +147,28 @@ for f! in DiffBase.INPLACE_ARRAY_TO_ARRAY_FUNCS end end +########################################## +# test specialized StaticArray codepaths # +########################################## + +x = rand(3, 3) +sx = StaticArrays.SArray{Tuple{3,3}}(x) +out = similar(x, 6, 9) +actual = ForwardDiff.jacobian(diff, x) + +@test ForwardDiff.jacobian(diff, sx) == actual + +ForwardDiff.jacobian!(out, diff, sx) + +@test out == actual + +result = DiffBase.JacobianResult(similar(x, 6), x) +sresult = DiffBase.JacobianResult(similar(sx, 6), sx) + +ForwardDiff.jacobian!(result, diff, x) +ForwardDiff.jacobian!(sresult, diff, sx) + +@test DiffBase.value(sresult) == DiffBase.value(result) +@test DiffBase.jacobian(sresult) == DiffBase.jacobian(result) + end # module From 5947f2a4010b386070359573df57e83efbc0319d Mon Sep 17 00:00:00 2001 From: Jarrett Revels Date: Thu, 27 Apr 2017 10:15:57 -0400 Subject: [PATCH 4/4] drop 32-bit windows support because it's annoying --- appveyor.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/appveyor.yml b/appveyor.yml index 0a6a1962..21bdcad2 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,6 +1,5 @@ environment: matrix: - - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe" - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe" branches: