add Hessian for composed Gauss function

test/gauss.jl

@@ -38,10 +38,25 @@ end
   @test TaylorTest.check(J_gauss, H_gauss, p; x=x)
 end
 
-@testset "Composed Gauss function" begin
-  p = 1 .+ 5 * rand(6)
+function check_symmetry(a)
+  R = true
+  if ndims(a) == 2
+    R = norm(a - a') < 1e-10
+  else
+    for i in axes(a, 1)
+      r = a[i, :, :]' == a[i, :, :]
+      R &= r
+      if !r
+        n = maximum(abs, a[i, :, :]' - a[i, :, :])
+        @warn "Tensor not symmetric at i=$i: $n"
+      end
+    end
+  end
+  return R
+end
 
-  x = p[3] + p[4] + 2p[5] + sqrt(p[6]) * (2rand() - 1)
+@testset "Composed Gauss function" begin
+  p = 1 .+ 5rand(6)
 
   p_aux = _p -> [_p[1] / _p[2], _p[3] + _p[4] + 2_p[5], _p[6]]
 
@@ -49,9 +64,28 @@ end
     0 0 1 1 2 0;
     0 0 0 0 0 1]
 
-  g = (_p; x=x) -> gauss(p_aux(_p); x=x)
-  Jg = (_p; x=x) -> J_gauss(p_aux(_p); x=x) * J_p_aux(p)
-  Hg = (_p; x=x) -> J_gauss(p_aux(_p); x=x) * J_p_aux(p)
+  H_p_aux = _p -> begin
+    h = zeros(3, 6, 6)
+    h[1, 1, 2] = h[1, 2, 1] = -1 / _p[2]^2
+    h[1, 2, 2] = 2_p[1] / _p[2]^3
+    h
+  end
+
+  g = (_p; x = 0.0) -> gauss(p_aux(_p); x=x)
+  Jg = (_p; x = 0.0) -> J_gauss(p_aux(_p); x=x) * J_p_aux(_p)
+  Hg = (_p; x = 0.0) -> begin
+    p = p_aux(_p)
+    _J_gauss = vec(J_gauss(p; x=x))
+    _H_p = H_p_aux(_p)
+    _J_p = J_p_aux(_p)
+    r = _J_p' * H_gauss(p; x=x) * _J_p
+    r .+= TO.tensorcontract([2, 3], _J_gauss, [1], _H_p, [1, 2, 3])
+    r
+  end
+
+  _p = p_aux(p)
+  x = _p[2] + _p[2] * 2 * (rand() - 0.5)
 
   @test TaylorTest.check(g, Jg, p; x=x)
+  @test TaylorTest.check(Jg, Hg, p; x=x)
 end