Merge pull request #15 from YingboMa/myb/hessian

Hessian

Author: YingboMa

src/aosoa.jl renamed to experiment/aosoa.jl

@@ -8,6 +8,6 @@ include("dual_context.jl")
 include("api.jl")
 
 # Experimental
-#include("aosoa.jl")
+#include("experiment/aosoa.jl")
 
 end # module

src/ForwardDiff2.jl

@@ -40,12 +40,13 @@ function D(f)
         return diffres
     end
     # scalar
-    function deriv(arg)
-        dualrun() do
-            dualized = map(x->Dual(x, one(x)), arg)
-            res = f(dualized)
-            return map(partials, res)
+    function deriv(x)
+        dx = one(x)
+        res = dualrun() do
+            dualized = Dual(x, dx)
+            f(dualized)
         end
+        return map(partials, res)
     end
     return deriv
 end

src/api.jl

@@ -1,8 +1,6 @@
 using StaticArrays: SVector
 
 partial_type(::Dual{T,V,P}) where {T,V,P} = P
-# TODO: Tagging?
-# TODO: Integrate better with SVector. Maybe even use SIMD.jl?
 
 struct DualArray{T,E,M,V<:AbstractArray,D<:AbstractArray} <: AbstractArray{E,M}
     data::V

src/dualarray.jl

@@ -4,10 +4,15 @@ using StaticArrays
 
 @testset begin
     @test D(sin)(1.0) == cos(1.0)
+    # Gradient
     @test D(x->[x, x^2])(3) == [1, 6]
     @test D(sum)([1,2,3]) == ones(3)'
+    # Jacobian
     @test D(x->@SVector([x[1]^x[2], x[3]^3, x[3]*x[2]*x[1]]))(@SVector[1,2,3.]) === @SMatrix [2.0 0 0; 0 0 27; 6 3 2]
     @test D(cumsum)(@SVector([1,2,3])) == @SMatrix [1 0 0; 1 1 0; 1 1 1]
     @test D(cumsum)([1,2,3]) == [1 0 0; 1 1 0; 1 1 1]
     @test D(x->@SVector([x[1], x[2]]))(@SVector([1,2,3])) === @SMatrix [1 0 0; 0 1 0]
+    # Hessian
+    @test D(D(x->x[1]^x[2] + x[3]^3 + x[3]*x[2]*x[1]))(@SVector[1,2,3]) === @SMatrix [2 4 2; 4 0 1; 2 1 18.]
+    @test D(D(x->x[1]^x[2] + x[3]^3 + x[3]*x[2]*x[1]))([1,2,3]) == [2 4 2; 4 0 1; 2 1 18.]
 end