compute_f = x -> [x[1]^2; x[2]^2 + x[3]; -x[3]^2 + x[4] + x[5]^3]

compute_Jf = x -> [
  2x[1] 0 0 0 0;
  0 2x[2] 1 0 0;
  0 0 -2x[3] 1 3*x[5]^2
]

function compute_Hf(x)
  Hf = zeros(3, 5, 5)

  Hf[1, 1, 1] = 2
  Hf[2, 2, 2] = 2
  Hf[3, 3, 3] = -2
  Hf[3, 5, 5] = 6 * x[5]

  return Hf
end

@testset "Tensors" begin
  m, n = 3, 5

  x = 2rand(n) .- 0.5

  @test TaylorTest.check(compute_f, compute_Jf, x)
  @test TaylorTest.check(compute_Jf, compute_Hf, x)

end

TP = TO.tensorproduct

function prod(A, B)
  return TP([1, 2, 3], A, [1, 2], B, [3])
end

function revprod(A, B)
  return TP([1, 3, 2], B, [3], A, [1, 2])
end

function symprod(A, B)
  return prod(A, B) + revprod(A, B)
end

@testset "Product of scalar and vector functions" begin
  v = x -> [exp(2 * x[1]), cos(x[1] * x[2]), 1.0]
  Jv = x -> [
    2*exp(2 * x[1]) 0;
    -x[2]*sin(x[1] * x[2]) -x[1]*sin(x[1] * x[2]);
    0 0
  ]
  Hv = x -> begin
    h = zeros(3, 2, 2)
    h[1, :, :] = [
      4*exp(2 * x[1]) 0;
      0 0
    ]

    h[2, :, :] = [
      -x[2]^2*cos(x[1] * x[2]) -sin(x[1] * x[2])-x[1]*x[2]*cos(x[1] * x[2]);
      -sin(x[1] * x[2])-x[1]*x[2]*cos(x[1] * x[2]) -x[1]^2*cos(x[1] * x[2])
    ]
    h
  end

  φ = x -> x[1] * x[2]^2
  ∇φ = x -> [x[2]^2; 2x[1] * x[2]]
  Hφ = x -> [
    0 2x[2];
    2x[2] 2*x[1]
  ]

  f = x -> v(x) * φ(x)
  Jf = x -> φ(x) * Jv(x) + v(x) * ∇φ(x)'
  Hf = x -> φ(x) * Hv(x) + symprod(Jv(x), ∇φ(x)) + TP(v(x), [1], Hφ(x), [2, 3])

  x = 2rand(2) .- 0.5

  @test TaylorTest.check(v, Jv, x)
  @test TaylorTest.check(Jv, Hv, x)

  @test TaylorTest.check(φ, ∇φ, x)
  @test TaylorTest.check(∇φ, Hφ, x)

  @test TaylorTest.check(f, Jf, x)
  @test TaylorTest.check(Jf, Hf, x)
end