text:
"""
Overview: collections of shape functions and their derivatives
"""

function select_quadrature_rule(etype::Symbol)
  #
  if (etype==:TE4)
    return get_quadrature_points(Val{:GLTET15})
  elseif (etype==:T10)
    return get_quadrature_points(Val{:GLTET15})
  elseif (etype==:HE8)
    return get_quadrature_points(Val{:GLHEX8})
  elseif (etype==:H20)
    return get_quadrature_points(Val{:GLHEX27})
  elseif (etype==:H27)
    return get_quadrature_points(Val{:GLHEX27})
  elseif (etype==:PE6)
    return get_quadrature_points(Val{:GLWED21})
  elseif (etype==:P15)
    return get_quadrature_points(Val{:GLWED21})
  elseif (etype==:P18)
    return get_quadrature_points(Val{:GLWED21})
  end
  #
  return nothing
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:T10}})
It fills N with the shape functions of a T10 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:T10}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:T10}})
  #
  gp4 = 1.0-gp[1]-gp[2]-gp[3]
  #

  N[1] = gp[1]*(2.0*gp[1]-1)
  N[2] = gp[2]*(2.0*gp[2]-1)
  N[3] = gp[3]*(2.0*gp[3]-1)
  N[4] = gp4*(2.0*gp4-1)
  #
  N[5]  =  4.0*gp[1]*gp[2]
  N[6]  =  4.0*gp[2]*gp[3]
  N[7]  =  4.0*gp[1]*gp[3]
  N[8]  =  4.0*gp[1]*gp4
  N[9]  =  4.0*gp[2]*gp4
  N[10] =  4.0*gp[3]*gp4
  #
  return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:T10}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:T10}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:T10}})
  #
  qa1 = 4.0*gp[1]
  qa2 = 4.0*gp[2]
  qa3 = 4.0*gp[3]
  qa4 = 4.0*(1-gp[1]-gp[2]-gp[3])
  #
  ∂N∂a[1,1] = -1.0+qa1
  ∂N∂a[1,2] = 0.0
  ∂N∂a[1,3] = 0.0
  #
  ∂N∂a[2,1] = 0.0
  ∂N∂a[2,2] = -1.0+qa2
  ∂N∂a[2,3] = 0.0
  #
  ∂N∂a[3,1] = 0.0
  ∂N∂a[3,2] = 0.0
  ∂N∂a[3,3] = -1.0+qa3
  #
  ∂N∂a[4,1] = 1.0-qa4
  ∂N∂a[4,2] = 1.0-qa4
  ∂N∂a[4,3] = 1.0-qa4
  #
  ∂N∂a[5,1] = qa2
  ∂N∂a[5,2] = qa1
  ∂N∂a[5,3] = 0.0
  #
  ∂N∂a[6,1] = 0.0
  ∂N∂a[6,2] = qa3
  ∂N∂a[6,3] = qa2
  #
  ∂N∂a[7,1] = qa3
  ∂N∂a[7,2] = 0.0
  ∂N∂a[7,3] = qa1
  #
  ∂N∂a[8,1] = qa4-qa1
  ∂N∂a[8,2] = -qa1
  ∂N∂a[8,3] = -qa1
  #
  ∂N∂a[9,1] = -qa2
  ∂N∂a[9,2] = qa4-qa2
  ∂N∂a[9,3] = -qa2
  #
  ∂N∂a[10,1] = -qa3
  ∂N∂a[10,2] = -qa3
  ∂N∂a[10,3] = qa4-qa3
  #
  return nothing
  #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:HE8}})
It fills N with the shape functions of a H8 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:HE8}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:HE8}})
  #
  x,y,z = gp
  # -1, +1, -1
  N[1] = 0.125*(1 - x)*(1 + y)*(1 - z)
  # +1, +1, -1
  N[2] = 0.125*(1 + x)*(1 + y)*(1 - z)
  # +1, -1, -1
  N[3] = 0.125*(1 + x)*(1 - y)*(1 - z)
  # -1, -1, -1
  N[4] = 0.125*(1 - x)*(1 - y)*(1 - z)
  # -1, +1, +1
  N[5] = 0.125*(1 - x)*(1 + y)*(1 + z)
  # +1, +1, +1
  N[6] = 0.125*(1 + x)*(1 + y)*(1 + z)
  # +1, -1, +1
  N[7] = 0.125*(1 + x)*(1 - y)*(1 + z)
  # -1, -1, +1
  N[8] = 0.125*(1 - x)*(1 - y)*(1 + z)
  #
  return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:HE8}})
It fills ∂N∂a with the derivatives of the shape functions of a H8 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:HE8}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:HE8}})
  #
  x,y,z = gp
  #
  ∂N∂a[1,1] = -0.125*(1 + y)*(1 - z)
  ∂N∂a[1,2] = +0.125*(1 - x)*(1 - z)
  ∂N∂a[1,3] = -0.125*(1 - x)*(1 + y)
  #
  ∂N∂a[2,1] = +0.125*(1 + y)*(1 - z)
  ∂N∂a[2,2] = +0.125*(1 + x)*(1 - z)
  ∂N∂a[2,3] = -0.125*(1 + x)*(1 + y)
  #
  ∂N∂a[3,1] = +0.125*(1 - y)*(1 - z)
  ∂N∂a[3,2] = -0.125*(1 + x)*(1 - z)
  ∂N∂a[3,3] = -0.125*(1 + x)*(1 - y)
  #
  ∂N∂a[4,1] = -0.125*(1 - y)*(1 - z)
  ∂N∂a[4,2] = -0.125*(1 - x)*(1 - z)
  ∂N∂a[4,3] = -0.125*(1 - x)*(1 - y)
  #
  ∂N∂a[5,1] = -0.125*(1 + y)*(1 + z)
  ∂N∂a[5,2] = +0.125*(1 - x)*(1 + z)
  ∂N∂a[5,3] = +0.125*(1 - x)*(1 + y)
  #
  ∂N∂a[6,1] = +0.125*(1 + y)*(1 + z)
  ∂N∂a[6,2] = +0.125*(1 + x)*(1 + z)
  ∂N∂a[6,3] = +0.125*(1 + x)*(1 + y)
  #
  ∂N∂a[7,1] = +0.125*(1 - y)*(1 + z)
  ∂N∂a[7,2] = -0.125*(1 + x)*(1 + z)
  ∂N∂a[7,3] = +0.125*(1 + x)*(1 - y)
  #
  ∂N∂a[8,1] = -0.125*(1 - y)*(1 + z)
  ∂N∂a[8,2] = -0.125*(1 - x)*(1 + z)
  ∂N∂a[8,3] = +0.125*(1 - x)*(1 - y)
  #
  return nothing
  #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:TE4}})
It fills N with the shape functions of a T4 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:TE4}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:TE4}})
  #
  N[1] = gp[1]
  N[2] = gp[2]
  N[3] = gp[3]
  N[4] = 1-gp[1]-gp[2]-gp[3]
  #
  return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:TE4}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:TE4}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:TE4}})
  #
  ∂N∂a[1,1] = 1.0
  ∂N∂a[1,2] = 0.0
  ∂N∂a[1,3] = 0.0
  #
  ∂N∂a[2,1] = 0.0
  ∂N∂a[2,2] = 1.0
  ∂N∂a[2,3] = 0.0
  #
  ∂N∂a[3,1] = 0.0
  ∂N∂a[3,2] = 0.0
  ∂N∂a[3,3] = 1.0
  #
  ∂N∂a[4,1] = -1.0
  ∂N∂a[4,2] = -1.0
  ∂N∂a[4,3] = -1.0
  #
  return nothing
  #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:PE6}})
It fills N with the shape functions of a PE6 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:PE6}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:PE6}})
  #
  x,y,z = gp
  #
  N[1] = x * 0.5 * (1.0 - z)
  N[2] = y * 0.5 * (1.0 - z)
  N[3] = (1.0-x-y) * 0.5 * (1.0 - z)
  N[4] = x * 0.5 * (1.0 + z)
  N[5] = y * 0.5 * (1.0 + z)
  N[6] = (1.0-x-y) * 0.5 * (1.0 + z)
  #
  return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:PE6}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:PE6}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:PE6}})
  #
  x,y,z = gp
  #
  ∂N∂a[1,1] = 0.5*(1 - z)
  ∂N∂a[1,2] = 0.0
  ∂N∂a[1,3] = -0.5*x
  #
  ∂N∂a[2,1] = 0.0
  ∂N∂a[2,2] = 0.5*(1 - z)
  ∂N∂a[2,3] = -0.5*y
  #
  ∂N∂a[3,1] = -0.5*(1 - z)
  ∂N∂a[3,2] = -0.5*(1 - z)
  ∂N∂a[3,3] = -0.5*(1 - x - y)
  #
  ∂N∂a[4,1] = +0.5*(1 + z)
  ∂N∂a[4,2] = +0.0
  ∂N∂a[4,3] = +0.5*x
  #
  ∂N∂a[5,1] = +0.0
  ∂N∂a[5,2] = +0.5*(1 + z)
  ∂N∂a[5,3] = +0.5*y
  #
  ∂N∂a[6,1] = -0.5*(1 + z)
  ∂N∂a[6,2] = -0.5*(1 + z)
  ∂N∂a[6,3] = +0.5*(1 - x - y)
  #
  return nothing
  #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:P15}})
It fills N with the shape functions of a PE6 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:P15}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:P15}})
    #
    @inbounds c1 = 1.0 - gp[1] - gp[2]
    @inbounds c2 = 1.0 - gp[3]
    @inbounds c3 = 1.0 + gp[3]
    #
    @inbounds da1 = 2.0*gp[1]
    @inbounds da2 = 2.0*gp[2]
    #
    @inbounds N[1] = -c1*c2*(da1+da2+gp[3])/2.0
    @inbounds N[2] = -gp[1]*c2*(2.0-da1+gp[3])/2.0
    @inbounds N[3] = -gp[2]*c2*(2.0-da2+gp[3])/2.0
    @inbounds N[4] = -c1*c3*(da1+da2-gp[3])/2.0
    @inbounds N[5] = -gp[1]*c3*(2.0-da1-gp[3])/2.0
    @inbounds N[6] = -gp[2]*c3*(2.0-da2-gp[3])/2.0
    @inbounds N[7] = da1*c1*c2
    @inbounds N[8] = da1*gp[2]*c2
    @inbounds N[9] = da2*c1*c2
    @inbounds N[10] = da1*c1*c3
    @inbounds N[11] = da1*gp[2]*c3
    @inbounds N[12] = da2*c1*c3
    @inbounds N[13] = c1*c2*c3
    @inbounds N[14] = gp[1]*c2*c3
    @inbounds N[15] = gp[2]*c2*c3
    return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:P15}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:P15}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:P15}})
  @inbounds  c1 = 1.0 - gp[1] - gp[2]
  @inbounds  c2 = 1.0 - gp[3]
  @inbounds  c3 = 1.0 + gp[3]
  #
  @inbounds  da1 = 2.0*gp[1]
  @inbounds  da2 = 2.0*gp[2]
  @inbounds  da3 = 2.0*gp[3]
  #
  @inbounds  qa1 = 4.0*gp[1]
  @inbounds  qa2 = 4.0*gp[2]
  #
  @inbounds  ∂N∂a[1,1] = c2*(qa1+qa2+gp[3]-2.0)/2.0
  @inbounds  ∂N∂a[1,2] = c2*(qa1+qa2+gp[3]-2.0)/2.0
  @inbounds  ∂N∂a[1,3] = c1*(da1+da2+da3-1.0)/2.0
  #
  @inbounds ∂N∂a[2,1] = c2*(qa1-gp[3]-2.0)/2.0
  @inbounds ∂N∂a[2,2] = 0.0
  @inbounds ∂N∂a[2,3] = gp[1]*(-da1+da3+1.0)/2.0
  #
  @inbounds ∂N∂a[3,1] = 0.0
  @inbounds ∂N∂a[3,2] = c2*(qa2-gp[3]-2.0)/2.0
  @inbounds ∂N∂a[3,3] = gp[2]*(-da2+da3+1.0)/2.0
  #
  @inbounds ∂N∂a[4,1] = c3*(qa1+qa2-gp[3]-2.0)/2.0
  @inbounds ∂N∂a[4,2] = c3*(qa1+qa2-gp[3]-2.0)/2.0
  @inbounds ∂N∂a[4,3] = -c1*(da1+da2-da3-1.0)/2.0
  #
  @inbounds ∂N∂a[5,1] = c3*(qa1+gp[3]-2.0)/2.0
  @inbounds ∂N∂a[5,2] = 0.0
  @inbounds ∂N∂a[5,3] = gp[1]*(da1+da3-1.0)/2.0
  #
  @inbounds ∂N∂a[6,1] = 0.0
  @inbounds ∂N∂a[6,2] = c3*(qa2+gp[3]-2.0)/2.0
  @inbounds ∂N∂a[6,3] = gp[2]*(da2+da3-1.0)/2.0
  #
  @inbounds ∂N∂a[7,1] = -c2*(da1+gp[2]-1.0)*2.0
  @inbounds ∂N∂a[7,2] = -gp[1]*c2*2.0
  @inbounds ∂N∂a[7,3] = -gp[1]*c1*2.0
  #
  @inbounds ∂N∂a[8,1] = gp[2]*c2*2.0
  @inbounds ∂N∂a[8,2] = gp[1]*c2*2.0
  @inbounds ∂N∂a[8,3] = -gp[1]*gp[2]*2.0
  #
  @inbounds ∂N∂a[9,1] = -gp[2]*c2*2.0
  @inbounds ∂N∂a[9,2] = -c2*(gp[1]+da2-1.0)*2.0
  @inbounds ∂N∂a[9,3] = -gp[2]*c1*2.0
  #
  @inbounds ∂N∂a[10,1] = -c3*(da1+gp[2]-1.0)*2.0
  @inbounds ∂N∂a[10,2] = -gp[1]*c3*2.0
  @inbounds ∂N∂a[10,3] = gp[1]*c1*2.0
  #
  @inbounds ∂N∂a[11,1] = gp[2]*c3*2.0
  @inbounds ∂N∂a[11,2] = gp[1]*c3*2.0
  @inbounds ∂N∂a[11,3] = gp[1]*gp[2]*2.0
  #
  @inbounds ∂N∂a[12,1] = -gp[2]*c3*2.0
  @inbounds ∂N∂a[12,2] = -c3*(gp[1]+da2-1.0)*2.0
  @inbounds ∂N∂a[12,3] = gp[2]*c1*2.0
  #
  @inbounds ∂N∂a[13,1] = -c2*c3
  @inbounds ∂N∂a[13,2] = -c2*c3
  @inbounds ∂N∂a[13,3] = -da3*c1
  #
  @inbounds ∂N∂a[14,1] = c2*c3
  @inbounds ∂N∂a[14,2] = 0.0
  @inbounds ∂N∂a[14,3] = -da1*gp[3]
  #
  @inbounds ∂N∂a[15,1] = 0.0
  @inbounds ∂N∂a[15,2] = c2*c3
  @inbounds ∂N∂a[15,3] = -da2*gp[3]
  #
  return nothing
  #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:H27}})
It fills N with the shape functions of a PE6 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:H27}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:H27}})
  #
  x1 = gp[1]*(gp[1]-1.0)/2.0
  x2 = (1-gp[1])*(1+gp[1])
  x3 = gp[1]*(gp[1]+1.0)/2.0
  #
  y1 = gp[2]*(gp[2]-1.0)/2.0
  y2 = (1-gp[2])*(1+gp[2])
  y3 = gp[2]*(gp[2]+1.0)/2.0
  #
  z1 = gp[3]*(gp[3]-1.0)/2.0
  z2 = (1-gp[3])*(1+gp[3])
  z3 = gp[3]*(gp[3]+1.0)/2.0
  #
  N[1] = x1*y1*z1
  N[2] = x3*y1*z1
  N[3] = x1*y3*z1
  N[4] = x3*y3*z1
  #
  N[5] = x1*y1*z3
  N[6] = x3*y1*z3
  N[7] = x1*y3*z3
  N[8] = x3*y3*z3
  #
  N[9] = x2*y1*z1
  N[10] = x1*y2*z1
  N[11] = x2*y2*z1
  N[12] = x3*y2*z1
  #
  N[13] = x2*y3*z1
  N[14] = x1*y1*z2
  N[15] = x2*y1*z2
  N[16] = x3*y1*z2
  #
  N[17] = x1*y2*z2
  N[18] = x2*y2*z2
  N[19] = x3*y2*z2
  N[20] = x1*y3*z2
  #
  N[21] = x2*y3*z2
  N[22] = x3*y3*z2
  N[23] = x2*y1*z3
  N[24] = x1*y2*z3
  N[25] = x2*y2*z3
  N[26] = x3*y2*z3
  N[27] = x2*y3*z3
  #
    return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:H27}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:H27}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:H27}})
  #
  x1 = gp[1]*(gp[1]-1.0)/2.0
  x2 = (1-gp[1])*(1+gp[1])
  x3 = gp[1]*(gp[1]+1.0)/2.0
  #
  y1 = gp[2]*(gp[2]-1.0)/2.0
  y2 = (1-gp[2])*(1+gp[2])
  y3 = gp[2]*(gp[2]+1.0)/2.0
  #
  z1 = gp[3]*(gp[3]-1.0)/2.0
  z2 = (1-gp[3])*(1+gp[3])
  z3 = gp[3]*(gp[3]+1.0)/2.0
  #
  dx1 = (2.0*gp[1]-1.0)/2.0
  dx2 = -2.0*gp[1]
  dx3 = (2.0*gp[1]+1.0)/2.0
  #
  dy1 = (2.0*gp[2]-1.0)/2.0
  dy2 = -2.0*gp[2]
  dy3 = (2.0*gp[2]+1.0)/2.0
  #
  dz1 = (2.0*gp[3]-1.0)/2.0
  dz2 = -2.0*gp[3]
  dz3 = (2.0*gp[3]+1.0)/2.0
  # 1
  ∂N∂a[1,1] = dx1*y1*z1
  ∂N∂a[1,2] = x1*dy1*z1
  ∂N∂a[1,3] = x1*y1*dz1
  # 2
  ∂N∂a[2,1] = dx3*y1*z1
  ∂N∂a[2,2] = x3*dy1*z1
  ∂N∂a[2,3] = x3*y1*dz1
  # 3
  ∂N∂a[3,1] = dx1*y3*z1
  ∂N∂a[3,2] = x1*dy3*z1
  ∂N∂a[3,3] = x1*y3*dz1
  # 4
  ∂N∂a[4,1] = dx3*y3*z1
  ∂N∂a[4,2] = x3*dy3*z1
  ∂N∂a[4,3] = x3*y3*dz1
  # 5
  ∂N∂a[5,1] = dx1*y1*z3
  ∂N∂a[5,2] = x1*dy1*z3
  ∂N∂a[5,3] = x1*y1*dz3
  # 6
  ∂N∂a[6,1] = dx3*y1*z3
  ∂N∂a[6,2] = x3*dy1*z3
  ∂N∂a[6,3] = x3*y1*dz3
  # 7
  ∂N∂a[7,1] = dx1*y3*z3
  ∂N∂a[7,2] = x1*dy3*z3
  ∂N∂a[7,3] = x1*y3*dz3
  # 8
  ∂N∂a[8,1] = dx3*y3*z3
  ∂N∂a[8,2] = x3*dy3*z3
  ∂N∂a[8,3] = x3*y3*dz3
  # 9
  ∂N∂a[9,1] = dx2*y1*z1
  ∂N∂a[9,2] = x2*dy1*z1
  ∂N∂a[9,3] = x2*y1*dz1
  # 10
  ∂N∂a[10,1] = dx1*y2*z1
  ∂N∂a[10,2] = x1*dy2*z1
  ∂N∂a[10,3] = x1*y2*dz1
  # 11
  ∂N∂a[11,1] = dx2*y2*z1
  ∂N∂a[11,2] = x2*dy2*z1
  ∂N∂a[11,3] = x2*y2*dz1
  # 12
  ∂N∂a[12,1] = dx3*y2*z1
  ∂N∂a[12,2] = x3*dy2*z1
  ∂N∂a[12,3] = x3*y2*dz1
  # 13
  ∂N∂a[13,1] = dx2*y3*z1
  ∂N∂a[13,2] = x2*dy3*z1
  ∂N∂a[13,3] = x2*y3*dz1
  # 14
  ∂N∂a[14,1] = dx1*y1*z2
  ∂N∂a[14,2] = x1*dy1*z2
  ∂N∂a[14,3] = x1*y1*dz2
  # 15
  ∂N∂a[15,1] = dx2*y1*z2
  ∂N∂a[15,2] = x2*dy1*z2
  ∂N∂a[15,3] = x2*y1*dz2
  # 16
  ∂N∂a[16,1] = dx3*y1*z2
  ∂N∂a[16,2] = x3*dy1*z2
  ∂N∂a[16,3] = x3*y1*dz2
  # 17
  ∂N∂a[17,1] = dx1*y2*z2
  ∂N∂a[17,2] = x1*dy2*z2
  ∂N∂a[17,3] = x1*y2*dz2
  # 18
  ∂N∂a[18,1] = dx2*y2*z2
  ∂N∂a[18,2] = x2*dy2*z2
  ∂N∂a[18,3] = x2*y2*dz2
  # 19
  ∂N∂a[19,1] = dx3*y2*z2
  ∂N∂a[19,2] = x3*dy2*z2
  ∂N∂a[19,3] = x3*y2*dz2
  # 20
  ∂N∂a[20,1] = dx1*y3*z2
  ∂N∂a[20,2] = x1*dy3*z2
  ∂N∂a[20,3] = x1*y3*dz2
  # 21
  ∂N∂a[21,1] = dx2*y3*z2
  ∂N∂a[21,2] = x2*dy3*z2
  ∂N∂a[21,3] = x2*y3*dz2
  # 22
  ∂N∂a[22,1] = dx3*y3*z2
  ∂N∂a[22,2] = x3*dy3*z2
  ∂N∂a[22,3] = x3*y3*dz2
  # 23
  ∂N∂a[23,1] = dx2*y1*z3
  ∂N∂a[23,2] = x2*dy1*z3
  ∂N∂a[23,3] = x2*y1*dz3
  # 24
  ∂N∂a[24,1] = dx1*y2*z3
  ∂N∂a[24,2] = x1*dy2*z3
  ∂N∂a[24,3] = x1*y2*dz3
  # 25
  ∂N∂a[25,1] = dx2*y2*z3
  ∂N∂a[25,2] = x2*dy2*z3
  ∂N∂a[25,3] = x2*y2*dz3
  # 26
  ∂N∂a[26,1] = dx3*y2*z3
  ∂N∂a[26,2] = x3*dy2*z3
  ∂N∂a[26,3] = x3*y2*dz3
  # 27
  ∂N∂a[27,1] = dx2*y3*z3
  ∂N∂a[27,2] = x2*dy3*z3
  ∂N∂a[27,3] = x2*y3*dz3
  #
  return nothing
    #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:H20}})
It fills N with the shape functions of a PE6 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:H20}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:H20}})
  #
  x,y,z = gp
  # [-1,-1,-1] [+1,-1,-1] [+1,+1,-1] [-1,+1,-1]
  N[1]  = (1.0-x) * (1.0-y) * (1.0-z) * (-x-y-z-2.0) / 8.0
  N[4]  = (1.0+x) * (1.0-y) * (1.0-z) * (+x-y-z-2.0) / 8.0
  N[3]  = (1.0+x) * (1.0+y) * (1.0-z) * (+x+y-z-2.0) / 8.0
  N[2]  = (1.0-x) * (1.0+y) * (1.0-z) * (-x+y-z-2.0) / 8.0
  # [-1,-1,+1] [+1,-1,+1] [+1,+1,+1] [-1,+1,+1]
  N[5]  = (1.0-x) * (1.0-y) * (1.0+z) * (-x-y+z-2.0) / 8.0
  N[8]  = (1.0+x) * (1.0-y) * (1.0+z) * (+x-y+z-2.0) / 8.0
  N[7]  = (1.0+x) * (1.0+y) * (1.0+z) * (+x+y+z-2.0) / 8.0
  N[6]  = (1.0-x) * (1.0+y) * (1.0+z) * (-x+y+z-2.0) / 8.0
  # [0,-1,-1] [0,+1,-1] [0,+1,+1] [0,-1,+1]
  N[12]  = (1-x^2) * (1.0-y) * (1.0-z) * 0.25
  N[10] = (1-x^2) * (1.0+y) * (1.0-z) * 0.25
  N[14] = (1-x^2) * (1.0+y) * (1.0+z) * 0.25
  N[16] = (1-x^2) * (1.0-y) * (1.0+z) * 0.25
  # [-1,0,-1] [+1,0,-1] [+1,0,+1] [-1,0,+1]
  N[9] = (1-y^2) * (1.0-x) * (1.0-z) * 0.25
  N[11] = (1-y^2) * (1.0+x) * (1.0-z) * 0.25
  N[15] = (1-y^2) * (1.0+x) * (1.0+z) * 0.25
  N[13] = (1-y^2) * (1.0-x) * (1.0+z) * 0.25
  # [-1,-1,0] [+1,-1,0] [+1,+1,0] [-1,+1,0]
  N[17] = (1-z^2) * (1.0-y) * (1.0-x) * 0.25
  N[20] = (1-z^2) * (1.0+y) * (1.0-x) * 0.25
  N[19] = (1-z^2) * (1.0+y) * (1.0+x) * 0.25
  N[18] = (1-z^2) * (1.0-y) * (1.0+x) * 0.25
  #
  return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:H20}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:H20}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:H20}})
  #
  x,y,z = gp
  #
  ∂N∂a[1,1] = 1/8*(-1 + y)*(-1 + z)*(1 + 2*x + y + z)
  ∂N∂a[1,2] = 1/8*(-1 + x)*(-1 + z)*(1 + x + 2*y + z)
  ∂N∂a[1,3] = 1/8*(-1 + x)*(-1 + y)*(1 + x + y + 2*z)
  # 2
  ∂N∂a[4,1] = -(1/8)*(-1 + y)*(-1 + z)*(1 - 2*x + y + z)
  ∂N∂a[4,2] = 1/8*(1 + x)*(-1 + x - 2*y - z)*(-1 + z)
  ∂N∂a[4,3] = 1/8*(1 + x)*(-1 + y)*(-1 + x - y - 2*z)
  # 3
  ∂N∂a[3,1] = -(1/8)*(1 + y)*(-1 + 2*x + y - z)*(-1 + z)
  ∂N∂a[3,2] = -(1/8)*(1 + x)*(-1 + x + 2*y - z)*(-1 + z)
  ∂N∂a[3,3] = -(1/8)*(1 + x)*(1 + y)*(-1 + x + y - 2*z)
  # 4
  ∂N∂a[2,1] = 1/8*(1 + y)*(-1 - 2*x + y - z)*(-1 + z)
  ∂N∂a[2,2] = -(1/8)*(-1 + x)*(-1 + z)*(1 + x - 2*y + z)
  ∂N∂a[2,3] = -(1/8)*(-1 + x)*(1 + y)*(1 + x - y + 2*z)
  # 5
  ∂N∂a[5,1] = -(1/8)*(-1 + y)*(1 + 2*x + y - z)*(1 + z)
  ∂N∂a[5,2] = -(1/8)*(-1 + x)*(1 + x + 2*y - z)*(1 + z)
  ∂N∂a[5,3] = -(1/8)*(-1 + x)*(-1 + y)*(1 + x + y - 2*z)
  # 6
  ∂N∂a[8,1] = 1/8*(-1 + y)*(1 - 2*x + y - z)*(1 + z)
  ∂N∂a[8,2] = -(1/8)*(1 + x)*(1 + z)*(-1 + x - 2*y + z)
  ∂N∂a[8,3] = -(1/8)*(1 + x)*(-1 + y)*(-1 + x - y + 2*z)
  # 7
  ∂N∂a[7,1] = 1/8*(1 + y)*(1 + z)*(-1 + 2*x + y + z)
  ∂N∂a[7,2] = 1/8*(1 + x)*(1 + z)*(-1 + x + 2*y + z)
  ∂N∂a[7,3] = 1/8*(1 + x)*(1 + y)*(-1 + x + y + 2*z)
  # 8
  ∂N∂a[6,1] = -(1/8)*(1 + y)*(1 + z)*(-1 - 2*x + y + z)
  ∂N∂a[6,2] = 1/8*(-1 + x)*(1 + x - 2*y - z)*(1 + z)
  ∂N∂a[6,3] = 1/8*(-1 + x)*(1 + y)*(1 + x - y - 2*z)
  # 9
  ∂N∂a[12,1] = -(1/2)*x*(-1 + y)*(-1 + z)
  ∂N∂a[12,2] = -(1/4)*(-1 + x^2)*(-1 + z)
  ∂N∂a[12,3] = -(1/4)*(-1 + x^2)*(-1 + y)
  # 10
  ∂N∂a[10,1] = 1/2*x*(1 + y)*(-1 + z)
  ∂N∂a[10,2] = 1/4*(-1 + x^2)*(-1 + z)
  ∂N∂a[10,3] = 1/4*(-1 + x^2)*(1 + y)
  # 11
  ∂N∂a[14,1] = -(1/2)*x*(1 + y)*(1 + z)
  ∂N∂a[14,2] = -(1/4)*(-1 + x^2)*(1 + z)
  ∂N∂a[14,3] = -(1/4)*(-1 + x^2)*(1 + y)
  # 12
  ∂N∂a[16,1] = 1/2*x*(-1 + y)*(1 + z)
  ∂N∂a[16,2] = 1/4*(-1 + x^2)*(1 + z)
  ∂N∂a[16,3] = 1/4*(-1 + x^2)*(-1 + y)
  # 13
  ∂N∂a[9,1] = -(1/4)*(-1 + y^2)*(-1 + z)
  ∂N∂a[9,2] = -(1/2)*(-1 + x)*y*(-1 + z)
  ∂N∂a[9,3] = -(1/4)*(-1 + x)*(-1 + y^2)
  # 14
  ∂N∂a[11,1] = 1/4*(-1 + y^2)*(-1 + z)
  ∂N∂a[11,2] = 1/2*(1 + x)*y*(-1 + z)
  ∂N∂a[11,3] = 1/4*(1 + x)*(-1 + y^2)
  # 15
  ∂N∂a[15,1] = -(1/4)*(-1 + y^2)*(1 + z)
  ∂N∂a[15,2] = -(1/2)*(1 + x)*y*(1 + z)
  ∂N∂a[15,3] = -(1/4)*(1 + x)*(-1 + y^2)
  # 16
  ∂N∂a[13,1] = 1/4*(-1 + y^2)*(1 + z)
  ∂N∂a[13,2] = 1/2*(-1 + x)*y*(1 + z)
  ∂N∂a[13,3] = 1/4*(-1 + x)*(-1 + y^2)
  # 17
  ∂N∂a[17,1] = -(1/4)*(-1 + y)*(-1 + z^2)
  ∂N∂a[17,2] = -(1/4)*(-1 + x)*(-1 + z^2)
  ∂N∂a[17,3] = -(1/2)*(-1 + x)*(-1 + y)*z
  # 18
  ∂N∂a[20,1] = 1/4*(-1 + y)*(-1 + z^2)
  ∂N∂a[20,2] = 1/4*(1 + x)*(-1 + z^2)
  ∂N∂a[20,3] = 1/2*(1 + x)*(-1 + y)*z
  # 19
  ∂N∂a[19,1] = -(1/4)*(1 + y)*(-1 + z^2)
  ∂N∂a[19,2] = -(1/4)*(1 + x)*(-1 + z^2)
  ∂N∂a[19,3] = -(1/2)*(1 + x)*(1 + y)*z
  # 20
  ∂N∂a[18,1] = 1/4*(1 + y)*(-1 + z^2)
  ∂N∂a[18,2] = 1/4*(-1 + x)*(-1 + z^2)
  ∂N∂a[18,3] = 1/2*(-1 + x)*(1 + y)*z
  #
  return nothing
    #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:P18}})
It fills N with the shape functions of a PE6 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:P18}} : element type specifier
"""
function N!(N::Array{Float64},gp::Tuple,::Type{Val{:P18}})
  #
  x,y,z = gp
  a1 = x
  a2 = y
  a3 = 1-x-y
  #
  N[1] = (2.0*a1-1.0)*a1 * z*(z-1)/2.0
  N[2] = (2.0*a2-1.0)*a2 * z*(z-1)/2.0
  N[3] = (2.0*a3-1.0)*a3 * z*(z-1)/2.0
  #
  N[4] = 4.0*a1*a2       * z*(z-1)/2.0
  N[5] = 4.0*a3*a2       * z*(z-1)/2.0
  N[6] = 4.0*a1*a3       * z*(z-1)/2.0
  #
  N[7] = (2.0*a1-1.0)*a1 * z*(z+1)/2.0
  N[8] = (2.0*a2-1.0)*a2 * z*(z+1)/2.0
  N[9] = (2.0*a3-1.0)*a3 * z*(z+1)/2.0
  #
  N[10] = 4.0*a1*a2       * z*(z+1)/2.0
  N[11] = 4.0*a3*a2       * z*(z+1)/2.0
  N[12] = 4.0*a1*a3       * z*(z+1)/2.0
  #
  N[13] = (2.0*a1-1.0)*a1 * (1+z)*(1-z)
  N[14] = (2.0*a2-1.0)*a2 * (1+z)*(1-z)
  N[15] = (2.0*a3-1.0)*a3 * (1+z)*(1-z)
  #
  N[16] = 4.0*a1*a2       * (1+z)*(1-z)
  N[17] = 4.0*a3*a2       * (1+z)*(1-z)
  N[18] = 4.0*a1*a3       * (1+z)*(1-z)
  #
  return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:P18}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:P18}} : element type specifier
"""
function ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:P18}})
  #
  x,y,z = gp
  #
  ∂N∂a[1,1] = 1/2.0*(-1.0 + 4.0* x)*(-1.0 + z)*z
  ∂N∂a[1,2] = 0.0
  ∂N∂a[1,3] = 1/2*x*(-1 + 2*x)*(-1 + 2*z)
  #
  ∂N∂a[2,1] = 0.0
  ∂N∂a[2,2] = 1/2*(-1 + 4*y)*(-1 + z)*z
  ∂N∂a[2,3] = 1/2*y*(-1 + 2*y)*(-1 + 2*z)
  #
  ∂N∂a[3,1] = 1/2*(-3 + 4*x + 4*y)*(-1 + z)*z
  ∂N∂a[3,2] = 1/2*(-3 + 4*x + 4*y)*(-1 + z)*z
  ∂N∂a[3,3] = 1/2*(1 + 2*x^2 - 3*y + 2*y^2 + x*(-3 + 4*y))*(-1 + 2*z)
  #
  ∂N∂a[4,1] = 2*y*(-1 + z)*z
  ∂N∂a[4,2] = 2*x*(-1 + z)*z
  ∂N∂a[4,3] = 2*x*y*(-1 + 2*z)
  #
  ∂N∂a[5,1] = -2*y*(-1 + z)*z
  ∂N∂a[5,2] = -2*(-1 + x + 2*y)*(-1 + z)*z
  ∂N∂a[5,3] = -2*y*(-1 + x + y)*(-1 + 2*z)
  #
  ∂N∂a[6,1] = -2*(-1 + 2*x + y)*(-1 + z)*z
  ∂N∂a[6,2] = -2*x*(-1 + z)*z
  ∂N∂a[6,3] = -2*x*(-1 + x + y)*(-1 + 2*z)
  #
  ∂N∂a[7,1] = 1/2*(-1 + 4*x)*z*(1 + z)
  ∂N∂a[7,2] = 0.0
  ∂N∂a[7,3] = 1/2*x*(-1 + 2*x)*(1 + 2*z)
  #
  ∂N∂a[8,1] = 0.0
  ∂N∂a[8,2] = 1/2*(-1 + 4*y)*z*(1 + z)
  ∂N∂a[8,3] = 1/2*y*(-1 + 2*y)*(1 + 2*z)
  #
  ∂N∂a[9,1] = 1/2*(-3 + 4*x + 4*y)*z*(1 + z)
  ∂N∂a[9,2] = 1/2*(-3 + 4*x + 4*y)*z*(1 + z)
  ∂N∂a[9,3] = 1/2*(1 + 2*x^2 - 3*y + 2*y^2 + x*(-3 + 4*y))*(1 + 2*z)
  #
  ∂N∂a[10,1] = 2*y*z*(1 + z)
  ∂N∂a[10,2] = 2*x*z*(1 + z)
  ∂N∂a[10,3] = 2*x*y*(1 + 2*z)
  #
  ∂N∂a[11,1] = -2*y*z*(1 + z)
  ∂N∂a[11,2] = -2*(-1 + x + 2*y)*z*(1 + z)
  ∂N∂a[11,3] = -2*y*(-1 + x + y)*(1 + 2*z)
  #
  ∂N∂a[12,1] = -2*(-1 + 2*x + y)*z*(1 + z)
  ∂N∂a[12,2] = -2*x*z*(1 + z)
  ∂N∂a[12,3] = -2*x*(-1 + x + y)*(1 + 2*z)
  #
  ∂N∂a[13,1] = -((-1 + 4*x)*(-1 + z^2))
  ∂N∂a[13,2] = 0.0
  ∂N∂a[13,3] = 2* (1 - 2*x)*x*z
  #
  ∂N∂a[14,1] = 0.0
  ∂N∂a[14,2] = -((-1 + 4*y)*(-1 + z^2))
  ∂N∂a[14,3] = 2*(1 - 2*y)*y*z
  #
  ∂N∂a[15,1] = -((-3 + 4*x + 4*y)*(-1 + z^2))
  ∂N∂a[15,2] = -((-3 + 4*x + 4*y)*(-1 + z^2))
  ∂N∂a[15,3] = -2*(1 + 2*x^2 - 3*y + 2*y^2 + x*(-3 + 4*y))*z
  #
  ∂N∂a[16,1] = -4*y*(-1 + z)*(1 + z)
  ∂N∂a[16,2] = -4*x*(-1 + z)*(1 + z)
  ∂N∂a[16,3] = -8*x*y*z
  #
  ∂N∂a[17,1] = 4*y*(-1 + z)*(1 + z)
  ∂N∂a[17,2] = 4*(-1 + x + 2*y)*(-1 + z^2)
  ∂N∂a[17,3] = 8*y*(-1 + x + y)*z
  #
  ∂N∂a[18,1] = 4*(-1 + 2*x + y)*(-1 + z^2)
  ∂N∂a[18,2] = 4*x*(-1 + z)*(1 + z)
  ∂N∂a[18,3] = 8*x*(-1 + x + y)*z
  #
  return nothing
    #
end

"""
> N!(N::Array{Float64},gp::Tuple,::Type{Val{:H27}})
It fills N with the shape functions of a PE6 element evaluated at gauss point gp.
- N : shape functions Array
- gp : gauss point coordinates
- :Type{Val{:H27}} : element type specifier
"""
function N27!(N::Array{Float64},gp::Tuple)
  #
  x1 = gp[1]*(gp[1]-1.0)/2.0
  x2 = (1-gp[1])*(1+gp[1])
  x3 = gp[1]*(gp[1]+1.0)/2.0
  #
  y1 = gp[2]*(gp[2]-1.0)/2.0
  y2 = (1-gp[2])*(1+gp[2])
  y3 = gp[2]*(gp[2]+1.0)/2.0
  #
  z1 = gp[3]*(gp[3]-1.0)/2.0
  z2 = (1-gp[3])*(1+gp[3])
  z3 = gp[3]*(gp[3]+1.0)/2.0
  #
  N[1] = x1*y1*z1
  N[2] = x3*y1*z1
  N[3] = x1*y3*z1
  N[4] = x3*y3*z1
  #
  N[5] = x1*y1*z3
  N[6] = x3*y1*z3
  N[7] = x1*y3*z3
  N[8] = x3*y3*z3
  #
  N[9] = x2*y1*z1
  N[10] = x1*y2*z1
  N[11] = x2*y2*z1
  N[12] = x3*y2*z1
  #
  N[13] = x2*y3*z1
  N[14] = x1*y1*z2
  N[15] = x2*y1*z2
  N[16] = x3*y1*z2
  #
  N[17] = x1*y2*z2
  N[18] = x2*y2*z2
  N[19] = x3*y2*z2
  N[20] = x1*y3*z2
  #
  N[21] = x2*y3*z2
  N[22] = x3*y3*z2
  N[23] = x2*y1*z3
  N[24] = x1*y2*z3
  N[25] = x2*y2*z3
  N[26] = x3*y2*z3
  N[27] = x2*y3*z3
  #
    return nothing
end

"""
> ∂N∂a!(∂N∂a::Array{Float64},gp::Tuple,::Type{Val{:H27}})
It fills ∂N∂a with the derivatives of the shape functions of a T4 element evaluated at gauss point gp.
- ∂N∂a : shape functions derivatives Array
- gp : gauss point coordinates
- :Type{Val{:H27}} : element type specifier
"""
function ∂N∂a27!(∂N∂a::Array{Float64},gp::Tuple)
  #
  x1 = gp[1]*(gp[1]-1.0)/2.0
  x2 = (1-gp[1])*(1+gp[1])
  x3 = gp[1]*(gp[1]+1.0)/2.0
  #
  y1 = gp[2]*(gp[2]-1.0)/2.0
  y2 = (1-gp[2])*(1+gp[2])
  y3 = gp[2]*(gp[2]+1.0)/2.0
  #
  z1 = gp[3]*(gp[3]-1.0)/2.0
  z2 = (1-gp[3])*(1+gp[3])
  z3 = gp[3]*(gp[3]+1.0)/2.0
  #
  dx1 = (2.0*gp[1]-1.0)/2.0
  dx2 = -2.0*gp[1]
  dx3 = (2.0*gp[1]+1.0)/2.0
  #
  dy1 = (2.0*gp[2]-1.0)/2.0
  dy2 = -2.0*gp[2]
  dy3 = (2.0*gp[2]+1.0)/2.0
  #
  dz1 = (2.0*gp[3]-1.0)/2.0
  dz2 = -2.0*gp[3]
  dz3 = (2.0*gp[3]+1.0)/2.0
  # 1
  ∂N∂a[1,1] = dx1*y1*z1
  ∂N∂a[1,2] = x1*dy1*z1
  ∂N∂a[1,3] = x1*y1*dz1
  # 2
  ∂N∂a[2,1] = dx3*y1*z1
  ∂N∂a[2,2] = x3*dy1*z1
  ∂N∂a[2,3] = x3*y1*dz1
  # 3
  ∂N∂a[3,1] = dx1*y3*z1
  ∂N∂a[3,2] = x1*dy3*z1
  ∂N∂a[3,3] = x1*y3*dz1
  # 4
  ∂N∂a[4,1] = dx3*y3*z1
  ∂N∂a[4,2] = x3*dy3*z1
  ∂N∂a[4,3] = x3*y3*dz1
  # 5
  ∂N∂a[5,1] = dx1*y1*z3
  ∂N∂a[5,2] = x1*dy1*z3
  ∂N∂a[5,3] = x1*y1*dz3
  # 6
  ∂N∂a[6,1] = dx3*y1*z3
  ∂N∂a[6,2] = x3*dy1*z3
  ∂N∂a[6,3] = x3*y1*dz3
  # 7
  ∂N∂a[7,1] = dx1*y3*z3
  ∂N∂a[7,2] = x1*dy3*z3
  ∂N∂a[7,3] = x1*y3*dz3
  # 8
  ∂N∂a[8,1] = dx3*y3*z3
  ∂N∂a[8,2] = x3*dy3*z3
  ∂N∂a[8,3] = x3*y3*dz3
  # 9
  ∂N∂a[9,1] = dx2*y1*z1
  ∂N∂a[9,2] = x2*dy1*z1
  ∂N∂a[9,3] = x2*y1*dz1
  # 10
  ∂N∂a[10,1] = dx1*y2*z1
  ∂N∂a[10,2] = x1*dy2*z1
  ∂N∂a[10,3] = x1*y2*dz1
  # 11
  ∂N∂a[11,1] = dx2*y2*z1
  ∂N∂a[11,2] = x2*dy2*z1
  ∂N∂a[11,3] = x2*y2*dz1
  # 12
  ∂N∂a[12,1] = dx3*y2*z1
  ∂N∂a[12,2] = x3*dy2*z1
  ∂N∂a[12,3] = x3*y2*dz1
  # 13
  ∂N∂a[13,1] = dx2*y3*z1
  ∂N∂a[13,2] = x2*dy3*z1
  ∂N∂a[13,3] = x2*y3*dz1
  # 14
  ∂N∂a[14,1] = dx1*y1*z2
  ∂N∂a[14,2] = x1*dy1*z2
  ∂N∂a[14,3] = x1*y1*dz2
  # 15
  ∂N∂a[15,1] = dx2*y1*z2
  ∂N∂a[15,2] = x2*dy1*z2
  ∂N∂a[15,3] = x2*y1*dz2
  # 16
  ∂N∂a[16,1] = dx3*y1*z2
  ∂N∂a[16,2] = x3*dy1*z2
  ∂N∂a[16,3] = x3*y1*dz2
  # 17
  ∂N∂a[17,1] = dx1*y2*z2
  ∂N∂a[17,2] = x1*dy2*z2
  ∂N∂a[17,3] = x1*y2*dz2
  # 18
  ∂N∂a[18,1] = dx2*y2*z2
  ∂N∂a[18,2] = x2*dy2*z2
  ∂N∂a[18,3] = x2*y2*dz2
  # 19
  ∂N∂a[19,1] = dx3*y2*z2
  ∂N∂a[19,2] = x3*dy2*z2
  ∂N∂a[19,3] = x3*y2*dz2
  # 20
  ∂N∂a[20,1] = dx1*y3*z2
  ∂N∂a[20,2] = x1*dy3*z2
  ∂N∂a[20,3] = x1*y3*dz2
  # 21
  ∂N∂a[21,1] = dx2*y3*z2
  ∂N∂a[21,2] = x2*dy3*z2
  ∂N∂a[21,3] = x2*y3*dz2
  # 22
  ∂N∂a[22,1] = dx3*y3*z2
  ∂N∂a[22,2] = x3*dy3*z2
  ∂N∂a[22,3] = x3*y3*dz2
  # 23
  ∂N∂a[23,1] = dx2*y1*z3
  ∂N∂a[23,2] = x2*dy1*z3
  ∂N∂a[23,3] = x2*y1*dz3
  # 24
  ∂N∂a[24,1] = dx1*y2*z3
  ∂N∂a[24,2] = x1*dy2*z3
  ∂N∂a[24,3] = x1*y2*dz3
  # 25
  ∂N∂a[25,1] = dx2*y2*z3
  ∂N∂a[25,2] = x2*dy2*z3
  ∂N∂a[25,3] = x2*y2*dz3
  # 26
  ∂N∂a[26,1] = dx3*y2*z3
  ∂N∂a[26,2] = x3*dy2*z3
  ∂N∂a[26,3] = x3*y2*dz3
  # 27
  ∂N∂a[27,1] = dx2*y3*z3
  ∂N∂a[27,2] = x2*dy3*z3
  ∂N∂a[27,3] = x2*y3*dz3
  #
  return nothing
    #
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ dcb19828-89b7-11eb-11f9-b354cb1e45c5
begin
	using PlutoUI: LocalResource
	LocalResource("../figures/universal_quantum_workflow.png", :alt => "Software stack on a gate-model quantum computer", :width => 500)
end

# ╔═╡ ab1f2d3a-4a37-48eb-8f05-e12263660b7a
begin
	using Yao, YaoPlots
	circuit = chain(2, put(1=>H), control(1,2=>X))
	plot(circuit)
end

# ╔═╡ bc3234fa-89b6-11eb-0b7e-f13e208505b3
include("./essentials.jl")

# ╔═╡ 457000db-5d94-47aa-8370-5d795c96e307
md"""
So far you mastered the notation of quantum mechanics and quantum computing, understood as much physics as needed to perform various operations on quantum states, and now you are ready to build quantum algorithms. In this notebook, we look at the basics of gate-model quantum computing, which is sometimes also referred to as universal quantum computing. Most academic and commercial efforts to build a quantum computer focus on this model: Alibaba, Baidu, Google, HP, IBM Q, Intel, IonQ, Microsoft, Rigetti Computing, and Tencent all aim at this, and the list keeps expanding. It remains unclear which implementation will prove scalable: superconducting chips, photonic systems, and ion traps are the most common types, each having its own advantages and disadvantages. We abstract away, we focus on the quantum algorithms irrespective of the physical implementation.

To get there, first we have to familiarize ourselves with some gates and what happens to those gates on quantum computers. The following diagram shows the software stack that bridges a problem we want to solve with the actual computational back-end [[1](#1)]:
"""

# ╔═╡ de0c2aa8-89b7-11eb-00ff-6d3a927ac4ab
md"""
First, we define the problem at a high-level and a suitable quantum algorithm is chosen. Then, we express the quantum algorithm as a quantum circuit composed of gates. This in turn has to be compiled to a specific quantum gate set available. The last step is to execute the final circuit either on a quantum processor or on a simulator.

The quantum algorithms we are interested in are about machine learning. In this notebook, we look at the levels below algorithms: the definition of circuits, their compilation, and the mapping to the hardware or a simulator.

# Defining circuits

Circuits are composed of qubit registers, gates acting on them, and measurements on the registers. To store the outcome of registers, many quantum computing libraries add classical registers to the circuits. Even by this language, you can tell that this is a very low level of programming a computer. It resembles the assembly language of digital computers, in which a program consists of machine code instructions.

Qubit registers are indexed from 0. We often just say qubit 0, qubit 1, and so on, to refer to the register containing a qubit. This is not to be confused with the actual state of the qubit, which can be $|0\rangle$, $|1\rangle$, or any superposition thereof. For instance, qubit 0 can be in the state $|1\rangle$.

Let's take a look at the gates. In digital computing, a processor transforms bit strings to bit strings with logical gates. Any bit string can be achieved with just two gates, which makes universal computations possible with simple operations composed only of these two types of gates. It is remarkable and surprising that the same is also true for quantum computers: any unitary operation can be decomposed into elementary gates, and three types of gates are sufficient. This is remarkable since we are talking about transforming continuous-valued probability amplitudes, not just discrete elements. Yet, this result is what provides the high-level theoretical foundation for being able to build a universal quantum computer at all.

Let's look at some common gates, some of which we have already seen. Naturally, all of these are unitary.

| Gate |Name                | $\;\;\;\;\;\;\;\;\;\;\;\;\;\;$ Matrix |
|------|--------------------|---------------------------------------------------------------------|
|   X  | Pauli-X or NOT gate|$\begin{bmatrix}0 & 1\\ 1& 0\end{bmatrix}$|
|   Z  | Pauli-Z  gate      |$\begin{bmatrix}1 & 0\\ 0& -1\end{bmatrix}$|
|  H   | Hadamard gate      |$\frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1\\ 1& -1\end{bmatrix}$|
| Rx($\theta$)| Rotation around X|$\begin{bmatrix}\cos(\theta/2) & -i \sin(\theta/2)\\ -i \sin(\theta / 2) & \cos(\theta / 2)\end{bmatrix}$|
| Ry($\theta$)| Rotation around Y|$\begin{bmatrix}\cos(\theta/2) & -\sin(\theta/2)\\ \sin(\theta / 2) & \cos(\theta / 2)\end{bmatrix}$|
| CNOT, CX | Controlled-NOT | $\begin{bmatrix}1 & 0 & 0 &0\\ 0 & 1 & 0 &0\\ 0 & 0 & 0 &1\\ 0 & 0 & 1 &0\end{bmatrix}$|

As we have seen before, the rotations correspond to axis defined in the Bloch sphere.

There should be one thing immediately apparent from the table: there are many, in fact, infinitely many single-qubit operations. The rotations, for instance, are parametrized by a continuous value. This is in stark contrast with digital circuits, where the only non-trivial single-bit gate is the NOT gate.

The CNOT gate is the only two-qubit gate in this list. It has a special role: we need two-qubit interactions to create entanglement. Let's repeat the circuit for creating the $|\phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$. We will have two qubit registers and two classical registers for measurement output. First, let's define the circuit and plot it:
"""

# ╔═╡ 0d4a96ff-ef37-41a3-bcdc-c67211a9ac53
md"""
Note that we can't just initialize the qubit registers in a state we fancy. All registers are initialized in $|0\rangle$ and creating a desired state is **part** of the circuit. In a sense, arbitrary state preparation is the same as universal quantum computation: the end of the calculation is a state that we desired to prepare. Some states are easier to prepare than others. The above circuit has only two gates to prepare our target state, so it is considered very easy.

Let us see what happens in this circuit. The Hadamard gate prepares an equal superposition $\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ in qubit 0. This qubit controls an X gate on qubit 1. Since qubit 0 is in the equal superposition after the Hadamard gate, it will not apply the X gate for the first part of the superposition ($|0\rangle$) and it will apply the X gate for the second part of the superposition ($|1\rangle$). Thus we create the final state $\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$, and we entangle the two qubit registers.

A digital computer's processing unit typically has 64-bit registers and it is able to perform universal calculations on bit strings. Any complex calculation is broken down into elementary 64-bit operations, either sequentially or in parallel execution. So you may wonder what is the deal with the thousands of qubits we expect from a quantum computer. Why can't a 64-qubit quantum computer be enough?

Entanglement is the easiest way to understand why we need so many qubits. Entanglement is a key resource in quantum computing and we want to make use of it. If we have 64-qubits and we want to entangle another one outside these 64 registers, we would have get rid of the qubit in one of the registers, potentially destroying a superposition and definitely destroying entanglement between that register and any other qubit on the chip. The only way to make use of superpositions and the strong correlations provided by entanglement is if the entire problem is on the quantum processing unit for the duration of the calculation.

This global nature of the calculation is also the reason why there is a focus on problems that are difficult to break down into elementary calculations. The travelling salesman problem is a great example: we need to consider all cities and all distances to minimize overall travel length.

To finish off the circuit, we could add a measurement to each qubit:
"""

# ╔═╡ fae3bb78-928c-4eab-b697-842a50b617f3
plot(chain(2, put(1:2=>circuit), Measure(2)))

# ╔═╡ 686c4aa3-f083-4d1d-9a5f-d911d96dc6e8
md"""
Finally, we can plot the statistics:
"""

# ╔═╡ c3919ddd-1a0e-4f41-9363-4545b10fdfcb
let q = zero_state(2)
	essentials.plot_hist(measure(q |> circuit, nshots=100))
end

# ╔═╡ 063fb7af-043b-4d2e-bbcb-b5809841eb7c
md"""
As we have seen before, 01 and 10 never appear.
"""

# ╔═╡ 84591132-dfda-4093-b135-37320bd2887f
md"""
# Compilation

The circuit is the way to describe a quantum algorithm. It may also contain some arbitrary single or two-qubit unitary and controlled versions thereof. A quantum compiler should be able to decompose these into elementary gates.
"""

# ╔═╡ 85385053-491c-437c-b98d-825206bb50c8
md"""
For instance, in Qiskit, you can get access to the general unitary using the $u3$ gate

``u3(\theta, \phi, \lambda) = \begin{pmatrix} \cos(\theta/2) & -e^{i\lambda}\sin(\theta/2) \\ e^{i\phi}\sin(\theta/2) & e^{i\lambda+i\phi}\cos(\theta/2) \end{pmatrix}.``

The compiler decomposes it into an actual gate sequence.
"""

# ╔═╡ 3b371926-1391-4737-a91e-d5aa62ea444e
md"""
This is one task of a quantum compiler. The next one is to translate the gates given in the circuit to the gates implemented in the hardware or the simulator. In the table above, we defined many gates, but a well-chosen set of three is sufficient for universality. For engineering constraints, typically one minimal set of universal gates is implemented in the hardware. It depends on the physical architecture which three.

At this point, the number of gates applied is probably already increasing: the decomposition of unitary will create many gates and the translation of gates is also likely to add more gates. An additional problem is the topology of the qubits: in some implementations not all qubit registers are connected to each other. The most popular implementation is superconducting qubits, which are manufactured on silicon chips just like any digital device you have. Since this is a quintessentially two dimensional layout, most qubits on the chip will not be connected. Here is an example topology of eight qubits on a superconducting quantum computer:
"""

# ╔═╡ 7996fb0c-89b7-11eb-012a-81d6a97709dc
LocalResource("../figures/eight_qubits.svg", :alt => "8-qubit topology", :width => 300)

# ╔═╡ 7ad8f61e-89b7-11eb-2af6-89c67e20e090
md"""
If we want to perform a two-qubit operations between two qubits that are not neighbouring, we have to perform SWAP operations to switch the qubit states between registers. A SWAP consists of three CNOT gates in a sequence.

The total number of gates at the end of the compilation reflects the true requirement of the hardware. *Circuit depth* is the number of time steps required to execute the circuit, assuming that gates acting on distinct qubits can operate in parallel. On current and near-term quantum computers, we want circuits to be shallow, otherwise decoherence or other forms of noise destroy our calculations.

We have to emphasize that the compilation depends on the backend. On the simulator, physical constraints do not apply. If we compile the circuit above, its depth will not increase:
"""

# ╔═╡ dc2b0813-946d-424b-b3e6-d90c5f12fa27
chain(2, put(1:2=>circuit), Measure(2))

# ╔═╡ 9e4553fb-613d-4399-b52e-bcfba4d8d5c2
md"""
# References

[1] M. Fingerhuth, T. Babej, P. Wittek. (2018). [Open source software in quantum computing](https://doi.org/10.1371/journal.pone.0208561). *PLOS ONE* 13(12):e0208561.
"""

# ╔═╡ Cell order:
# ╟─457000db-5d94-47aa-8370-5d795c96e307
# ╟─dcb19828-89b7-11eb-11f9-b354cb1e45c5
# ╟─de0c2aa8-89b7-11eb-00ff-6d3a927ac4ab
# ╠═ab1f2d3a-4a37-48eb-8f05-e12263660b7a
# ╟─0d4a96ff-ef37-41a3-bcdc-c67211a9ac53
# ╠═fae3bb78-928c-4eab-b697-842a50b617f3
# ╟─686c4aa3-f083-4d1d-9a5f-d911d96dc6e8
# ╠═bc3234fa-89b6-11eb-0b7e-f13e208505b3
# ╠═c3919ddd-1a0e-4f41-9363-4545b10fdfcb
# ╟─063fb7af-043b-4d2e-bbcb-b5809841eb7c
# ╟─84591132-dfda-4093-b135-37320bd2887f
# ╟─85385053-491c-437c-b98d-825206bb50c8
# ╟─3b371926-1391-4737-a91e-d5aa62ea444e
# ╟─7996fb0c-89b7-11eb-012a-81d6a97709dc
# ╟─7ad8f61e-89b7-11eb-2af6-89c67e20e090
# ╠═dc2b0813-946d-424b-b3e6-d90c5f12fa27
# ╟─9e4553fb-613d-4399-b52e-bcfba4d8d5c2

ext: jl

max_stars_repo_name: VarLad/qml-mooc

max_issues_repo_name: VarLad/qml-mooc

max_forks_repo_name: VarLad/qml-mooc


text:
@testset "Construction" begin
    @test isa(ReferenceSequence(dna""), ReferenceSequence)
    @test isa(ReferenceSequence(dna"ACGTN"^30), ReferenceSequence)
    @test isa(ReferenceSequence("ACGTN"^30), ReferenceSequence)
    @test_throws Exception ReferenceSequence(dna"ACGRT")
    @test_throws Exception ReferenceSequence("ACGRT")
end

@testset "Conversion" begin
    seq = LongDNASeq(random_dna(100))
    refseq = ReferenceSequence(seq)
    @test refseq == seq
    @test LongDNASeq(refseq) == seq
end


text:
Qx, x = FlintQQ["x"]
K1, a1 = NumberField(x^6 - x^5 + x^4 - x^3 + x^2 - x + 1, "a") # totally complex
K2, a2 = NumberField(x^6 - x^5 - 7*x^4 + 2*x^3 + 7*x^2 - 2*x - 1, "a") # totally real
K3, a3 = NumberField(x^6 - x^5 - x^4 + 4*x^3 + 3*x^2 - 1, "a") # signature (2, 2)

@testset "Totally real/complex" begin
  @test @inferred istotally_complex(K1)
  @test @inferred !istotally_complex(K2)
  @test @inferred !istotally_complex(K3)

  @test @inferred !istotally_real(K1)
  @test @inferred istotally_real(K2)
  @test @inferred !istotally_real(K3)
end

@testset "conjugates_arb" begin
  CC = AcbField(64)

  co = @inferred conjugates(10001 * a1 + 1, 256)

  @test overlaps(co[1], CC("[9011.589647892093681 +/- 9.30e-16]", "[4339.271274914698763 +/- 9.19e-16]"))
  @test Hecke.radiuslttwopower(co[1], -256)
  @test overlaps(co[2], CC("[-6234.521508389194039 +/- 8.09e-16]", "[7819.096656162766117 +/- 5.63e-16]"))
  @test Hecke.radiuslttwopower(co[2], -256)
  @test overlaps(co[3], CC("[2226.431860497100357 +/- 4.55e-16]", "[9750.25404973041789 +/- 4.36e-15]"))
  @test Hecke.radiuslttwopower(co[3], -256)
  @test overlaps(co[4], CC("[9011.589647892093681 +/- 9.30e-16]", "[-4339.271274914698763 +/- 9.19e-16]"))
  @test Hecke.radiuslttwopower(co[4], -256)
  @test overlaps(co[5], CC("[-6234.521508389194039 +/- 8.09e-16]", "[-7819.096656162766117 +/- 5.63e-16]"))
  @test Hecke.radiuslttwopower(co[5], -256)
  @test overlaps(co[6], CC("[2226.431860497100357 +/- 4.55e-16]", "[-9750.25404973041789 +/- 4.36e-15]"))
  @test Hecke.radiuslttwopower(co[6], -256)

  RR = ArbField(64)

  co = @inferred conjugates_real(a3, 32)
  @test length(co) == 2

  @test Hecke.radiuslttwopower(co[1], -32)
  @test Hecke.overlaps(co[1], RR("[-1.221417721 +/- 5.71e-10]"))
  @test Hecke.radiuslttwopower(co[2], -32)
  @test Hecke.overlaps(co[2], RR("[0.4665 +/- 4.52e-5]"))

  co = @inferred conjugates_complex(a3, 32)
  @test length(co) == 2
  @test Hecke.radiuslttwopower(co[1], -32)
  @test Hecke.overlaps(co[1], CC("[-0.542287022 +/- 3.58e-10]", "[0.460349889 +/- 4.60e-10]" ))
  @test Hecke.radiuslttwopower(co[2], -32)
  @test Hecke.overlaps(co[2], CC("[1.419725855 +/- 5.85e-10]", "[1.205211655 +/- 7.03e-10]"))

  colog = @inferred conjugates_log(a1, 16)
  @test length(colog) == 3
  @test contains(colog[1], 0)
  @test contains(colog[2], 0)
  @test contains(colog[3], 0)

  mink = @inferred minkowski_map(a3, 32)
  @test length(mink) == 6
  co = conjugates(a3, 32)
  sqrt2 = sqrt(RR(2))
  @test overlaps(mink[1], real(co[1]))
  @test overlaps(mink[2], real(co[2]))
  @test overlaps(mink[3], sqrt2 * real(co[3]))
  @test overlaps(mink[4], sqrt2 * imag(co[3]))
  @test overlaps(mink[5], sqrt2 * real(co[4]))
  @test overlaps(mink[6], sqrt2 * imag(co[4]))
end

@testset "T2" begin
  t = @inferred t2(a1, 64)
  @test contains(t, 6)
  @test Hecke.radiuslttwopower(t, -64)

  t = @inferred t2(a3^2, 128)
  @test overlaps(t, parent(t)("[26.8413130281669208224194865303660830716653210471312375675880809411649241305119684287180747356640152133866595977425833731119286780918091464992340699076157 +/- 6.82e-152]"))
end

@testset "signs" begin
  s = @inferred Hecke._signs(a3)
  @test s == Int[-1, 1]

  s = @inferred Hecke._signs(a3^2)
  @test s == Int[1, 1]

  s = @inferred Hecke._signs(a1)
  @test s == Int[]

  s = @inferred Hecke._signs(a2)
  @test s == Int[-1, -1, -1, 1, 1, 1]

  @test_throws ErrorException Hecke._signs(0*a1)
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Auxiliary functions related to coordinate transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export FK5, ITRF, PEF, TOD, MOD, GCRF, J2000, TEME, CIRS, TIRS

# Auxiliary functions to define the conversion models.
@inline FK5()   = Val{:FK5}

# Auxiliary functions to define the reference frames.
@inline ITRF()  = Val{:ITRF}
@inline PEF()   = Val{:PEF}
@inline TOD()   = Val{:TOD}
@inline MOD()   = Val{:MOD}
@inline GCRF()  = Val{:GCRF}
@inline J2000() = Val{:J2000}
@inline TEME()  = Val{:TEME}
@inline CIRS()  = Val{:CIRS}
@inline TIRS()  = Val{:TIRS}


text:
# In this tutorial, we present the use of Stopping to specify an "unrelated"
# stopping criterion. A typical example is the so-called case of "overfitting",
# where the optimization process is actually designed to approximate another problem.

using ADNLPModels, LinearAlgebra, NLPModels, Plots, Printf, Random, Stopping

# In this tutorial, we will use the classical steepest descent method with an Armijo line-search.

import Stopping.armijo
function armijo(xk, dk, fk, slope, f)
  t = 1.0
  fk_new = f(xk + dk)
  while f(xk + t * dk) > fk + 1.0e-4 * t * slope
    t /= 1.5
    fk_new = f(xk + t * dk)
  end
  return t, fk_new
end

function steepest_descent(stp :: NLPStopping)

  xk = stp.current_state.x
  fk, gk = objgrad(stp.pb, xk)

  OK = update_and_start!(stp, fx = fk, gx = gk)

  @printf "%2s %9s %7s %7s %7s\n" "k" "fk" "||∇f(x)||" "t" "λ"
  @printf "%2d %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score)
  while !OK
    dk = - gk
    slope = dot(dk, gk)
    t, fk = armijo(xk, dk, fk, slope, x->obj(stp.pb, x))
    xk += t * dk
    fk, gk = objgrad(stp.pb, xk)

    OK = update_and_stop!(stp, x = xk, fx = fk, gx = gk)

    @printf "%2d %9.2e %7.1e %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score) t slope
  end
  return stp
end

# We also generate randomly some data.

Random.seed!(1234)
m, n = 50, 10
A  = rand(m, n)
b  = A * ones(n)
D  = vcat(A, A) + vcat(zeros(m,n), rand(m,n))
Db = vcat(b, b)
rperm = shuffle(1:2m)
Dr = D[rperm,:]
Dbr = Db[rperm,:]

# We initialize two different problems evaluating respectively a training set, and a test set.

f(x, A, b, λ) = norm(A * x - b)^2 + λ * norm(x)^2
train_pb = ADNLPModel(x -> f(x, Dr[1:m,:], Dbr[1:m], 1e-2), zeros(n))
test_pb = ADNLPModel(x -> f(x, Dr[m+1:2m,:], Dbr[m+1:2m], 0.0), zeros(n))

# In this tutorial, the motivation is to use the solver on the `test_pb`, and track the efficiency of the computed solution on the `train_pb`.
# We specialize the stopping structure to store the objective function of both problems in the `stopping_user_struct`, and evaluate them in the `stp.meta.user_check_func!`.

train_check(stp, b) = begin
  stp.stopping_user_struct[:test_obj][stp.meta.nb_of_stop+1] = obj(stp.stopping_user_struct[:test], stp.current_state.x)
  stp.stopping_user_struct[:train_obj][stp.meta.nb_of_stop+1] = obj(stp.pb, stp.current_state.x)
end
train_obj = NaN * ones(101)
test_obj  = NaN * ones(101)
train_stp = NLPStopping(train_pb, user_struct = Dict(:test => test_pb, :train_obj => train_obj, :test_obj => test_obj), user_check_func! = train_check, max_iter = 100)

# We now run the steepest descent algorithm with the `train_stp` stopping.

steepest_descent(train_stp)

# On the following plot in logarithmic scale, we can see that after a number of iterations, the progress made by the solver are no longer improving the other problem.

#hcat(log.(train_stp.stopping_user_struct[:train_obj]), log.(train_stp.stopping_user_struct[:test_obj]))
plot(log.(train_stp.stopping_user_struct[:train_obj]), label=["train obj"])
plot!(log.(train_stp.stopping_user_struct[:test_obj]), label=["test obj"], title="overfitting")

# ![](../assets/overfitting/overfitting.png)

# To overcome this issue, one possibility is to stop the solver when the second is no longer being minimized. The only modification is to set the entry `meta.stopbyuser` to `true`.

train_check_opt(stp, b) = begin
  k = stp.meta.nb_of_stop
  stp.stopping_user_struct[:test_obj][k+1] = obj(stp.stopping_user_struct[:test], stp.current_state.x)
  stp.stopping_user_struct[:train_obj][k+1] = obj(stp.pb, stp.current_state.x)
  diff = k!=0 && stp.stopping_user_struct[:test_obj][k+1] - stp.stopping_user_struct[:train_obj][k+1] > 10
  inc  = k!=0 && stp.stopping_user_struct[:test_obj][k+1] > stp.stopping_user_struct[:test_obj][k]
  if diff && inc
    stp.meta.stopbyuser = true
  end
end
train_stp.meta.user_check_func! = train_check_opt
reset!(train_stp.pb)
reinit!(train_stp, rstate=true, x = zeros(10))
steepest_descent(train_stp)

nb_iter = train_stp.meta.nb_of_stop
hcat(log.(train_stp.stopping_user_struct[:train_obj][1:nb_iter+1]), log.(train_stp.stopping_user_struct[:test_obj][1:nb_iter+1]))
plot(log.(train_stp.stopping_user_struct[:train_obj][1:nb_iter+1]), label=["train obj"])
plot!(log.(train_stp.stopping_user_struct[:test_obj][1:nb_iter+1]), label=["test obj"], title="overfitting")

# ![](../assets/overfitting/overfitting_cut.png)


text:
function pack(X, U)
    z = zeros(num_var)
    for t = 1:T
        z[x_idx[t]] = X[t]
        t == T && continue
        z[u_idx[t]] = U[t]
    end
    return z
end

function unpack(z)
    x = [z[x_idx[t]] for t = 1:T]
    u = [z[u_idx[t]] for t = 1:T-1]
    return x, u
end

function primal_bounds(x1, xT)
    zl = -Inf * ones(num_var)
    zu = Inf * ones(num_var)

	zl[x_idx[1]] = x1
	zu[x_idx[1]] = x1

	# zl[x_idx[T]] = xT
	# zu[x_idx[T]] = xT

    return zl, zu
end

function constraint_bounds()
    cl = zeros(num_con)
    cu = zeros(num_con)
    return cl, cu
end

function eval_objective(prob::MOIGeneralProblem, z)
	J = 0.0
	for t = 1:T
		if t < T
			x = z[x_idx[t]]
			u = z[u_idx[t]]
			J += objective(x, u, t)
		else
			x = z[x_idx[t]]
			J += objective(x, nothing, t)
		end
	end
	# parameter regularization
	θ = z[p_idx]
	J += 1.0 * θ' * θ
	return J
end

function eval_objective_gradient!(∇J, z, prob::MOIGeneralProblem)
    ∇J .= 0.0
	for t = 1:T
		if t < T
			x = z[x_idx[t]]
			u = z[u_idx[t]]
			obj_x(a) = objective(a, u, t)
			obj_u(a) = objective(x, a, t)
			∇J[x_idx[t]] = ForwardDiff.gradient(obj_x, x)
			∇J[u_idx[t]] = ForwardDiff.gradient(obj_u, u)
		else
			x = z[x_idx[t]]
			obj_xT(a) = objective(a, nothing, t)
			∇J[x_idx[t]] = ForwardDiff.gradient(obj_xT, x)
		end
	end
	θ = z[p_idx]
	∇J[p_idx] = 2.0 * θ

    return nothing
end

function eval_constraint!(c, z, prob::MOIGeneralProblem)
	# dynamics
    for t = 1:T-1
		x = z[x_idx[t]]
		u = z[u_idx[t]]
		y = z[x_idx[t+1]]
		c[n_idx[t]] = dynamics(y, x, u, t)
	end

	# policy
	θ = z[p_idx]
	for t = 1:T-1
		x = z[x_idx[t]]
		u = z[u_idx[t]]
		c[num_dyn .+ (m_idx[t])] = u - policy(x, θ)
	end
    return nothing
end

function eval_constraint_jacobian!(∇c, z, prob::MOIGeneralProblem)
	shift = 0

    for t = 1:T-1
        x = view(z, x_idx[t])
        u = view(z, u_idx[t])
        y = view(z, x_idx[t + 1])

        dyn_x(a) = dynamics(y, a, u, t)
        dyn_u(a) = dynamics(y, x, a, t)
        dyn_y(a) = dynamics(a, x, u, t)

        r_idx = n_idx[t]

        s = n * n
        ∇c[shift .+ (1:s)] = vec(ForwardDiff.jacobian(dyn_x, x))
        shift += s

        s = n * m
        ∇c[shift .+ (1:s)] = vec(ForwardDiff.jacobian(dyn_u, u))
        shift += s

        s = n * n
        ∇c[shift .+ (1:s)] = vec(ForwardDiff.jacobian(dyn_y, y))
        shift += s
    end

	# policy
	θ = z[p_idx]
	for t = 1:T-1
		x = z[x_idx[t]]
		u = z[u_idx[t]]

		px(a) = policy(a, θ)
		pθ(a) = policy(x, a)

		r_idx = num_dyn .+ (m_idx[t])

		s = m * m
		∇c[shift .+ (1:s)] = vec(Diagonal(ones(m)))
		shift += s

		s = m * n
		∇c[shift .+ (1:s)] = -vec(ForwardDiff.jacobian(px, x))
		shift += s

		s = m * p
		∇c[shift .+ (1:s)] = -vec(ForwardDiff.jacobian(pθ, θ))
		shift += s
	end

    return nothing
end

function sparsity_jacobian(prob::MOIGeneralProblem)

	row = []
	col = []

	for t = 1:T-1
		r_idx = n_idx[t]
		row_col!(row, col, r_idx, x_idx[t])
		row_col!(row, col, r_idx, u_idx[t])
		row_col!(row, col, r_idx, x_idx[t + 1])
	end

	for t = 1:T-1
		r_idx = num_dyn .+ (m_idx[t])
		row_col!(row, col, r_idx, u_idx[t])
		row_col!(row, col, r_idx, x_idx[t])
		row_col!(row, col, r_idx, p_idx)
	end

	return collect(zip(row, col))
end

## Trajectory Optimization

# time
T = 11
h = 1.0

# dynamics
n = 2
m = 1

A = SMatrix{n,n}([1.0 1.0; 0.0 1.0])
B = SMatrix{n,m}([0.0; 1.0])

function dynamics(y, x, u, t)
	y - (A * x + B * u)
end

# objective
function objective(x, u, t)
	J = 0.0
	if t == T
		Q = Diagonal(100.0 * ones(n))
		J += x' * Q  * x
	else
		Q = Diagonal(ones(n))
		R = Diagonal(0.1 * ones(m))
		J += x' * Q  * x
		J += transpose(u) * R * u
	end
	return J
end

# policy
function policy(x, θ)
	M = reshape(θ, m, n)
	M * x
end

p = m * n

# trajectory indices
x_idx = [(t - 1) * (n + m) .+ (1:n) for t = 1:T]
u_idx = [(t - 1) * (n + m) + n .+ (1:m) for t = 1:T-1]
n_idx = [(t - 1) * n .+ (1:n) for t = 1:T-1]
m_idx = [(t - 1) * m .+ (1:m) for t = 1:T-1]
p_idx = n * T + m * (T - 1) .+ (1:p)

# MOI setup
num_x = n * T
num_u = m * (T - 1)
num_var = num_x + num_u + p

num_dyn = n * (T - 1)
num_policy = m * (T - 1)
num_con = num_dyn + num_policy

# initialization
x1 = [1.0; 0.0]
xT = [0.0; 0.0]
X = linear_interpolation(x1, xT, T)
U = [0.01 * randn(m) for t = 1:T-1]

z = [pack(X, U); 0.01 * randn(p)]

prob = MOIGeneralProblem(num_var, num_con, primal_bounds(x1, xT), constraint_bounds())
prob_moi = moi_problem(prob)

z_sol, info = solve(prob_moi, z,
		tol = 1.0e-3,
		c_tol = 1.0e-3,
		max_iter = 1000,
		nlp = :ipopt)

x_sol, u_sol = unpack(z_sol)
θ_sol = z_sol[num_x + num_u .+ (1:p)]
M = reshape(θ_sol, m, n)

using Plots

plot(hcat(x_sol...)')
plot(hcat(u_sol..., u_sol[end])', linetype = :steppost)

## simulate policy
x1_sim = copy(x1) + 0.01 * randn(n)

x_hist = [x1_sim]
u_hist = []
for t = 1:T-1
	push!(u_hist, M * x_hist[end])
	# push!(u_hist, u_sol[t])

	push!(x_hist, A * x_hist[end] + B * u_hist[end])
end

plot(hcat(x_hist...)')
plot(hcat(u_hist..., u_hist[end])', linetype = :steppost)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
const ADF = AbstractDataField

# Two-dimensional fields
@inline Base.getindex( r::ADF{Nothing, Y, Z},    i, j, k) where {Y, Z} = @inbounds r.data[1, j, k]
@inline Base.setindex!(r::ADF{Nothing, Y, Z}, d, i, j, k) where {Y, Z} = @inbounds r.data[1, j, k] = d
@inline Base.getindex( r::ADF{X, Nothing, Z},    i, j, k) where {X, Z} = @inbounds r.data[i, 1, k]
@inline Base.setindex!(r::ADF{X, Nothing, Z}, d, i, j, k) where {X, Z} = @inbounds r.data[i, 1, k] = d
@inline Base.getindex( r::ADF{X, Y, Nothing},    i, j, k) where {X, Y} = @inbounds r.data[i, j, 1]
@inline Base.setindex!(r::ADF{X, Y, Nothing}, d, i, j, k) where {X, Y} = @inbounds r.data[i, j, 1] = d

@inline Base.getindex( r::ADF{Nothing, Y, Z},    j, k) where {Y, Z} = @inbounds r.data[1, j, k]
@inline Base.setindex!(r::ADF{Nothing, Y, Z}, d, j, k) where {Y, Z} = @inbounds r.data[1, j, k] = d
@inline Base.getindex( r::ADF{X, Nothing, Z},    i, k) where {X, Z} = @inbounds r.data[i, 1, k]
@inline Base.setindex!(r::ADF{X, Nothing, Z}, d, i, k) where {X, Z} = @inbounds r.data[i, 1, k] = d
@inline Base.getindex( r::ADF{X, Y, Nothing},    i, j) where {X, Y} = @inbounds r.data[i, j, 1]
@inline Base.setindex!(r::ADF{X, Y, Nothing}, d, i, j) where {X, Y} = @inbounds r.data[i, j, 1] = d

# One-dimensional fields
@inline Base.getindex( r::ADF{X, Nothing, Nothing},    i, j, k) where X = @inbounds r.data[i, 1, 1]
@inline Base.setindex!(r::ADF{X, Nothing, Nothing}, d, i, j, k) where X = @inbounds r.data[i, 1, 1] = d
@inline Base.getindex( r::ADF{Nothing, Y, Nothing},    i, j, k) where Y = @inbounds r.data[1, j, 1]
@inline Base.setindex!(r::ADF{Nothing, Y, Nothing}, d, i, j, k) where Y = @inbounds r.data[1, j, 1] = d
@inline Base.getindex( r::ADF{Nothing, Nothing, Z},    i, j, k) where Z = @inbounds r.data[1, 1, k]
@inline Base.setindex!(r::ADF{Nothing, Nothing, Z}, d, i, j, k) where Z = @inbounds r.data[1, 1, k] = d

@inline Base.getindex( r::ADF{X, Nothing, Nothing},    i) where X = @inbounds r.data[i, 1, 1]
@inline Base.setindex!(r::ADF{X, Nothing, Nothing}, d, i) where X = @inbounds r.data[i, 1, 1] = d
@inline Base.getindex( r::ADF{Nothing, Y, Nothing},    j) where Y = @inbounds r.data[1, j, 1]
@inline Base.setindex!(r::ADF{Nothing, Y, Nothing}, d, j) where Y = @inbounds r.data[1, j, 1] = d
@inline Base.getindex( r::ADF{Nothing, Nothing, Z},    k) where Z = @inbounds r.data[1, 1, k]
@inline Base.setindex!(r::ADF{Nothing, Nothing, Z}, d, k) where Z = @inbounds r.data[1, 1, k] = d

# Zero-dimensional fields
@inline Base.getindex( r::ADF{Nothing, Nothing, Nothing},    i, j, k) = @inbounds r.data[1, 1, 1]
@inline Base.setindex!(r::ADF{Nothing, Nothing, Nothing}, d, i, j, k) = @inbounds r.data[1, 1, 1] = d

@inline Base.getindex( r::ADF{Nothing, Nothing, Nothing}   ) = @inbounds r.data[1, 1, 1]
@inline Base.setindex!(r::ADF{Nothing, Nothing, Nothing}, d) = @inbounds r.data[1, 1, 1] = d


text:
@testset "Solvers" begin
  @testset "CookieCutter" begin
    problem = SimulationProblem(RegularGrid{Float64}(100,100), (:facies => Int, :property => Float64), 3)
    solver = CookieCutter(DummySimSolver(:facies => NamedTuple()),
                          [0 => DummySimSolver(), 1 => DummySimSolver()])

    Random.seed!(1234)
    solution = solve(problem, solver)

    if visualtests
      gr(size=(800,600))
      @plottest plot(solution) joinpath(datadir,"CookieCutter.png") !istravis
    end
  end

  @testset "SeqSim" begin
    sdomain = RegularGrid{Float64}(100, 100)
    problem = SimulationProblem(sdomain, :var => Float64, 3)
    solver = SeqSim(:var => (estimator=DummyEstimator(),
                             neighborhood=BallNeighborhood(10.),
                             minneighbors=1, maxneighbors=10,
                             marginal=Normal(), path=LinearPath(sdomain)))

    Random.seed!(1234)
    solution = solve(problem, solver)

    if visualtests
      gr(size=(900,300))
      @plottest plot(solution) joinpath(datadir,"SeqSim.png") !istravis
    end
  end
end


text:
#Define SunType and object
struct SunType <: abstractStar end
Sun = Star(SunType)

#Data, as previously used in AstroDynBase
naif_id(::Type{SunType}) = 10
mu(::Type{SunType}) = 1.32712440041e11
mean_radius(::Type{SunType}) = 696000.0
polar_radius(::Type{SunType}) = 696000.0
equatorial_radius(::Type{SunType}) = 696000.0
parent(::Type{SunType}) = SSB

#Extra data
alpha(::Type{SunType}) = deg2rad.([286.13, 0., 0.])
delta(::Type{SunType}) = deg2rad.([63.87, 0., 0.])
omega(::Type{SunType}) = deg2rad.([84.176, 14.18440, 0.])

alpha_nut_prec(::Type{SunType}) = 0.
delta_nut_prec(::Type{SunType}) = 0.
omega_nut_prec(::Type{SunType}) = 0.
theta0_nut_prec(::Type{SunType}) = 0.
theta1_nut_prec(::Type{SunType}) = 0.


text:
# # using Revise
using Test
using Logging
using LarSurf
using SparseArrays
# using Pio3d
# using BenchmarkTools

interactive = false
if interactive
    using Plasm
end
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation
# Logging.configure(level==Logging.Debug)

# include("../src/LarSurf.jl")
# include("../src/block.jl")

@testset "test smoothing" begin
    # xystep = 1
    # zstep = 1
    # threshold = 4000;
    # pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
    # datap = Pio3d.read3d(pth)
    # data3d = datap["data3d"]
    # segmentation = data3d .> threshold

    segmentation = LarSurf.generate_minecraft_kidney(30)

    block_size = [5,5,5]
    basicmodel, Flin, larmodel = LarSurf.get_surface_grid_per_block(segmentation, block_size; return_all=true)
    someV, topology = basicmodel
    FV = topology[1]
    bigV = someV
    # @info "smoothing_FV test" bigV FV FV[1]
    # Quads -> smoothing -> triangulation
    Vqs = LarSurf.Smoothing.smoothing_FV(bigV, FV, 0.6)
    @test typeof(Vqs) <: Array

    # Plasm.View((Vqs * 100, [FV]))

    FVtri = LarSurf.triangulate_quads(FV)
    if interactive
        Plasm.View((Vqs * 100, [FVtri]))
    end
    # Smoothi accoarding to trinagle faces
    # Quads -> triangulation -> Smoothing
    Vts = LarSurf.Smoothing.smoothing_FV(bigV, FVtri, 0.6)
    @test typeof(Vts) <: Array
    @test size(Vts,1) == 3

    for i = 1:5
        Vts = LarSurf.Smoothing.smoothing_FV(Vts, FVtri, 0.6)
    end
    if interactive
        Plasm.View((Vqs * 100, [FVtri]))
    end
    # taubin
    Vts = LarSurf.Smoothing.smoothing_FV(bigV, FVtri, 0.6)
    for i = 1:5
        Vts = LarSurf.Smoothing.smoothing_FV(Vts, FVtri, 0.6)
        Vts = LarSurf.Smoothing.smoothing_FV(Vts, FVtri, -0.2)
    end
    if interactive
        Plasm.View((Vqs * 100, [FVtri]))
    end
    @test typeof(Vts) <: Array
    @test size(Vts,1) == 3
end

@testset "test iterative smoothing" begin
    segmentation = LarSurf.generate_minecraft_kidney(30)

    block_size = [5,5,5]
    basicmodel = LarSurf.get_surface_grid_per_block(segmentation, block_size; return_all=false)
    someV, topology = basicmodel
    FV = topology[1]
    bigV = someV
    # @info "smoothing_FV test" bigV FV FV[1]
    # Quads -> smoothing -> triangulation
    Vqs = LarSurf.Smoothing.smoothing_FV(bigV, FV, 0.6, 5)
    @test typeof(Vqs) <: Array
    @test size(Vqs,1) == 3
end
@testset "test iterative smoothing small" begin
    # segmentation = LarSurf.generate_minecraft_kidney(6)
    # segmentation = zeros(np)
    #
    # Vertex IDs first is on level 2 in Z-axis and second is on level 3
    #         6,5   15,1
    #  4         *------*
    #            |      |
    #       11,12|  3,14|     9,10
    #  3         *------*------*
    #            |      |      |
    #        1,16|  13,8|   4,7|
    #  2         *------*------*
    #
    #            2      3      4

    seg = zeros(Int8, 4,4,4)
    seg[2,2,2] = 1
    seg[2,3,2] = 1
    seg[3,2,2] = 1
    segmentation = seg .> 0

    block_size = [2,2,2]
    basicmodel = LarSurf.get_surface_grid_per_block(segmentation, block_size; return_all=false)
    someV, topology = basicmodel
    FV = topology[1]
    bigV = someV
    # @info "smoothing_FV test" bigV FV FV[1]
    # Quads -> smoothing -> triangulation
    testEV = LarSurf.Smoothing.get_EV_quads(FV)
    Vqs = LarSurf.Smoothing.smoothing_FV(bigV, FV, 0.6, 1)
    @test typeof(Vqs) <: Array
    @test size(Vqs,1) == 3
end

@testset "test iterative taubin smoothing" begin
    segmentation = LarSurf.generate_minecraft_kidney(30)

    block_size = [5,5,5]
    basicmodel = LarSurf.get_surface_grid_per_block(segmentation, block_size; return_all=false)
    someV, topology = basicmodel
    FV = topology[1]
    bigV = someV
    # @info "smoothing_FV test" bigV FV FV[1]
    # Quads -> smoothing -> triangulation
    Vqs = LarSurf.Smoothing.smoothing_FV_taubin(bigV, FV, 0.6, -0.4, 5)
    @test typeof(Vqs) <: Array
    @test size(Vqs,1) == 3
end

@testset "test get_EV" begin
    larmodel = LarSurf.Lar.cuboidGrid([1,1,2], true)
    V, topology = larmodel
    VV, EVorig, FV, CV = topology
    deleteat!(FV,10)
    # FVsurf = cat([FV[1:9], FV[11:end]])
    EVset = LarSurf.Smoothing.get_EV_quads(FV)
    EVdouble = LarSurf.Smoothing.get_EV_quads2(FV)
    @test length(EVset[1]) == 2
    @test length(EVdouble[1]) == 2
    # @test EVset::Array{Array,1}

    # Plasm.View((V, [EVset, FV]))
    # Plasm.View((V, [EV, FV]))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1FlockerVolumeSource <: SwaggerModel
    datasetName::Any # spec type: Union{ Nothing, String } # spec name: datasetName
    datasetUUID::Any # spec type: Union{ Nothing, String } # spec name: datasetUUID

    function IoK8sApiCoreV1FlockerVolumeSource(;datasetName=nothing, datasetUUID=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1FlockerVolumeSource, Symbol("datasetName"), datasetName)
        setfield!(o, Symbol("datasetName"), datasetName)
        validate_property(IoK8sApiCoreV1FlockerVolumeSource, Symbol("datasetUUID"), datasetUUID)
        setfield!(o, Symbol("datasetUUID"), datasetUUID)
        o
    end
end # type IoK8sApiCoreV1FlockerVolumeSource

const _property_map_IoK8sApiCoreV1FlockerVolumeSource = Dict{Symbol,Symbol}(Symbol("datasetName")=>Symbol("datasetName"), Symbol("datasetUUID")=>Symbol("datasetUUID"))
const _property_types_IoK8sApiCoreV1FlockerVolumeSource = Dict{Symbol,String}(Symbol("datasetName")=>"String", Symbol("datasetUUID")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1FlockerVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1FlockerVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1FlockerVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1FlockerVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1FlockerVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1FlockerVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1FlockerVolumeSource)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1FlockerVolumeSource }, name::Symbol, val)
end


text:
struct PlotList # TODO: Is PlotList needed?
    plots::Dict
end
function PlotList()
    PlotList(Dict())
end

Base.show(io::IO, mm::MIME"text/html", p::PlotList) =
    show(io, mm, "$(Plots.backend()) with $(length(p.plots)) plots, named $(keys(p.plots))")
Base.show(io::IO, mm::MIME"text/plain", p::PlotList) =
    show(io, mm, "$(Plots.backend()) with $(length(p.plots)) plots, named $(keys(p.plots))")

# the fundamental struct for plotting
struct PGData
    data::Dict{Symbol, DataFrames.DataFrame}
    time::Union{StepRange{Dates.DateTime}, Vector{Dates.DateTime}}
end

#### Generation Names ####
function get_generation_variable_names(
    results::IS.Results;
    names::Union{Nothing, Vector{Symbol}} = nothing,
)
    existing_names = PSI.get_existing_variables(results)
    names = isnothing(names) ? existing_names : [n for n in names if n ∈ existing_names]
    filter_names = Vector{Symbol}()
    for name in names
        name_string = string(name)
        if any(startswith.(name_string, SUPPORTEDVARPREFIX))
            if any(endswith.(name_string, all_subtypes(PSY.Generator)))
                push!(filter_names, name)
            end
        elseif any(name .== keys(SLACKVARS))
            push!(filter_names, name)
        end
    end
    return filter_names
end

function get_storage_variable_names(
    results::IS.Results;
    names::Union{Nothing, Vector{Symbol}} = nothing,
)
    existing_names = PSI.get_existing_variables(results)
    names = isnothing(names) ? existing_names : [n for n in names if n ∈ existing_names]
    filter_names = Vector{Symbol}()
    for name in names
        name_string = string(name)
        if any(startswith.(name_string, SUPPORTEDVARPREFIX))
            if any(endswith.(name_string, all_subtypes(PSY.Storage)))
                push!(filter_names, name)
            end
        end
    end
    return filter_names
end

function get_generation_parameter_names(
    results::IS.Results;
    names::Union{Nothing, Vector{Symbol}} = nothing,
)
    existing_names = PSI.get_existing_parameters(results)
    names = isnothing(names) ? existing_names : [n for n in names if n ∈ existing_names]
    filter_names = Vector{Symbol}()
    for name in names
        name_string = string(name)
        if any(startswith.(name_string, SUPPORTEDPARAMPREFIX))
            if any(occursin.(all_subtypes(PSY.Generator), name_string))
                push!(filter_names, name)
            end
        end
    end
    return filter_names
end

#### Load Names ####
function get_load_variable_names(
    results::IS.Results;
    names::Union{Nothing, Vector{Symbol}} = nothing,
)
    existing_names = PSI.get_existing_variables(results)
    names = isnothing(names) ? existing_names : [n for n in names if n ∈ existing_names]
    filter_names = Vector{Symbol}()
    for name in names
        name_string = string(name)
        if any(startswith.(name_string, SUPPORTEDVARPREFIX))
            if any(occursin.(SUPPORTEDLOADPARAMS, name_string))
                push!(filter_names, name)
            end
        end
    end
    return filter_names
end

function get_load_parameter_names(
    results::IS.Results;
    names::Union{Nothing, Vector{Symbol}} = nothing,
)
    existing_names = PSI.get_existing_parameters(results)
    names = isnothing(names) ? existing_names : [n for n in names if n ∈ existing_names]
    filter_names = Vector{Symbol}()
    for name in names
        name_string = string(name)
        if any(startswith.(name_string, SUPPORTEDPARAMPREFIX))
            if any(occursin.(SUPPORTEDLOADPARAMS, name_string))
                push!(filter_names, name)
            end
        end
    end
    return filter_names
end

#### Service Names ####
function get_service_variable_names(
    results::IS.Results;
    names::Union{Nothing, Vector{Symbol}} = nothing,
)
    existing_names = PSI.get_existing_variables(results)
    names = isnothing(names) ? existing_names : [n for n in names if n ∈ existing_names]
    filter_names = Vector{Symbol}()
    for name in names
        name_string = string(name)
        if any(endswith.(name_string, SUPPORTEDSERVICEPARAMS))
            push!(filter_names, name)
        end
    end
    return filter_names
end

function get_service_parameter_names(
    results::IS.Results;
    names::Union{Nothing, Vector{Symbol}} = nothing,
)
    existing_names = PSI.get_existing_parameters(results)
    names = isnothing(names) ? existing_names : [n for n in names if n ∈ existing_names]
    filter_names = Vector{Symbol}()
    for name in names
        name_string = string(name)
        if any(endswith.(name_string, SUPPORTEDSERVICEPARAMS))
            push!(filter_names, name)
        end
    end
    return filter_names
end

#### get data ###

function _get_matching_param(var_name)
    param_name = Symbol(
        join([SUPPORTEDPARAMPREFIX, split(string(var_name), PSI_NAME_DELIMITER)[end]]),
    )
    return param_name
end

function _get_matching_var(param_name)
    var_name = Symbol(
        join(
            [
                s for s in
                vcat(split.(split(string(param_name), "max_active_power"), "_")...) if
                !isempty(s)
            ],
            PSI_NAME_DELIMITER,
        ),
    )
    return var_name
end

no_datetime(df::DataFrames.DataFrame) = df[:, propertynames(df) .!== :DateTime]

function add_fixed_parameters!(
    variables::Dict{Symbol, DataFrames.DataFrame},
    parameters::Dict{Symbol, DataFrames.DataFrame},
)
    # fixed output should be added to plots when there exists a parameter of the form
    # :P__max_active_power__* but there is no corresponding :P__* variable
    for (param_name, param) in parameters
        var_name = _get_matching_var(param_name)
        if !haskey(variables, var_name)
            mult = any(occursin.(NEGATIVE_PARAMETERS, string(param_name))) ? -1.0 : 1.0
            variables[var_name] = param
            variables[var_name][:, propertynames(param) .!== :DateTime] .*= mult
        end
    end
end

# finds the parameters corresponding to variables for curtailment calculations
function _curtailment_parameters(parameters::Vector{Symbol}, variables::Vector{Symbol})
    curtailment_parameters =
        Vector{NamedTuple{(:parameter, :variable), Tuple{Symbol, Symbol}}}()
    for var in variables
        if any(
            length.(
                collect.(eachmatch.(Regex(CURTAILMENTPARAMSUFFIX), string.(parameters))),
            ) .> 1,
        )
            prefix = CURTAILMENTPARAMPREFIX
        else
            prefix = SUPPORTEDPARAMPREFIX
        end
        var_param = Symbol(
            join([prefix, split(string(var), "_")[end], CURTAILMENTPARAMSUFFIX], "_"),
        )
        if var_param in parameters
            push!(curtailment_parameters, (parameter = var_param, variable = var))
        end
    end
    return curtailment_parameters
end

function _filter_curtailment!(
    variable_values::Dict,
    parameter_values::Dict,
    curtailment_parameters::Vector{
        NamedTuple{(:parameter, :variable), Tuple{Symbol, Symbol}},
    },
)
    for curtailment in curtailment_parameters
        if !haskey(variable_values, curtailment.variable)
            variable_values[curtailment.variable] = parameter_values[curtailment.parameter]
        else
            curt =
                parameter_values[curtailment.parameter] .-
                variable_values[curtailment.variable]
            if haskey(variable_values, :Curtailment)
                variable_values[:Curtailment] =
                    hcat(variable_values[:Curtailment], no_datetime(curt))
            else
                variable_values[:Curtailment] = curt
            end
        end
    end
end

function get_generation_data(results::R; kwargs...) where {R <: PSI.PSIResults}
    initial_time = get(kwargs, :initial_time, nothing)
    len = get(kwargs, :horizon, get(kwargs, :len, nothing))
    names = get(kwargs, :names, nothing)
    curtailment = get(kwargs, :curtailment, true)
    storage = get(kwargs, :storage, true)

    if curtailment && !isnothing(names)
        @warn "Cannot guarantee curtailment calculations with specified names"
    end

    var_names = get_generation_variable_names(results; names = names)

    if storage
        var_names = vcat(var_names, get_storage_variable_names(results; names))
    end

    param_names = get_generation_parameter_names(results; names = names)

    variables = PSI.read_realized_variables(
        results;
        names = var_names,
        initial_time = initial_time,
        len = len,
    )
    parameters = PSI.read_realized_parameters(
        results;
        names = param_names,
        initial_time = initial_time,
        len = len,
    )

    add_fixed_parameters!(variables, parameters)

    if curtailment
        curtailment_parameters = _curtailment_parameters(param_names, var_names)
        _filter_curtailment!(variables, parameters, curtailment_parameters)
    end

    timestamps =
        PSI.get_realized_timestamps(results; initial_time = initial_time, len = len)
    return PGData(variables, timestamps)
end

function get_load_data(results::R; kwargs...) where {R <: PSI.PSIResults}
    initial_time = get(kwargs, :initial_time, nothing)
    len = get(kwargs, :horizon, get(kwargs, :len, nothing))
    names = get(kwargs, :names, nothing)

    var_names = get_load_variable_names(results; names = names)
    param_names = get_load_parameter_names(results; names = names)

    variables = PSI.read_realized_variables(
        results;
        names = var_names,
        initial_time = initial_time,
        len = len,
    )
    parameters = PSI.read_realized_parameters(
        results;
        names = param_names,
        initial_time = initial_time,
        len = len,
    )

    add_fixed_parameters!(variables, parameters)

    timestamps =
        PSI.get_realized_timestamps(results; initial_time = initial_time, len = len)
    return PGData(variables, timestamps)
end

################################### INPUT DEMAND #################################

function _get_loads(system::PSY.System, bus::PSY.Bus)
    return [
        load for load in PSY.get_components(PSY.PowerLoad, system, PSY.get_available) if
        PSY.get_bus(load) == bus
    ]
end
function _get_loads(system::PSY.System, agg::T) where {T <: PSY.AggregationTopology}
    return PSY.get_components_in_aggregation_topology(PSY.PowerLoad, system, agg)
end
function _get_loads(system::PSY.System, load::PSY.PowerLoad)
    return [load]
end
function _get_loads(system::PSY.System, sys::PSY.System)
    return PSY.get_components(PSY.PowerLoad, system, PSY.get_available)
end

get_base_power(system::PSY.System) = PSY.get_base_power(system)
get_base_power(results::PSI.SimulationProblemResults) = IS.get_base_power(results)
get_base_power(results::PSI.ProblemResults) = results.base_power

function get_load_data(
    system::PSY.System;
    aggregation::Union{
        Type{PSY.PowerLoad},
        Type{PSY.Bus},
        Type{PSY.System},
        Type{<:PSY.AggregationTopology},
    } = PSY.PowerLoad,
    kwargs...,
)
    aggregation_components =
        aggregation == PSY.System ? [system] : PSY.get_components(aggregation, system)
    if isempty(aggregation_components)
        throw(ArgumentError("System does not have type $aggregation."))
    end
    horizon = get(kwargs, :horizon, PSY.get_forecast_horizon(system))
    initial_time = get(kwargs, :initial_time, PSY.get_forecast_initial_timestamp(system))
    parameters = Dict{Symbol, DataFrames.DataFrame}()
    PSY.set_units_base_system!(system, "SYSTEM_BASE")
    for agg in aggregation_components
        loads = _get_loads(system, agg)
        length(loads) == 0 && continue
        colname = aggregation == PSY.System ? "System" : PSY.get_name(agg)
        load_values = DataFrames.DataFrame()
        for load in loads
            f = PSY.get_time_series_values( # TODO: this isn't applying the scaling factors
                PSY.Deterministic,
                load,
                "max_active_power",
                start_time = initial_time,
                len = horizon,
            )
            load_values[:, PSY.get_name(load)] = f
        end
        parameters[Symbol(colname)] = load_values
    end
    time_range =
        range(initial_time, step = PSY.get_time_series_resolution(system), length = horizon)

    return PGData(parameters, time_range)
end

function get_service_data(results::R; kwargs...) where {R <: PSI.PSIResults}
    initial_time = get(kwargs, :initial_time, nothing)
    len = get(kwargs, :horizon, get(kwargs, :len, nothing))
    names = get(kwargs, :names, nothing)

    var_names = get_service_variable_names(results; names = names)

    variables = PSI.read_realized_variables(
        results;
        names = var_names,
        initial_time = initial_time,
        len = len,
    )
    timestamps =
        PSI.get_realized_timestamps(results; initial_time = initial_time, len = len)

    return PGData(variables, timestamps)
end

#### result combination and aggregation ####

"""
aggregates and combines data into single DataFrame

# Example

```julia
PG.combine_categories(gen_uc.data)
```

"""
function combine_categories(
    data::Union{Dict{Symbol, DataFrames.DataFrame}, Dict{String, DataFrames.DataFrame}};
    names::Union{Vector{String}, Vector{Symbol}, Nothing} = nothing,
    aggregate::Union{Function, Nothing} = nothing,
)
    aggregate = isnothing(aggregate) ? x -> sum(x, dims = 2) : aggregate
    names = isnothing(names) ? keys(data) : names
    values = []
    keep_names = []
    for k in names
        if !isempty(data[k])
            push!(values, aggregate(Matrix(no_datetime(data[k]))))
            push!(keep_names, k)
        end
    end
    data = hcat(values...)
    keep_names = string.(keep_names)
    isempty(data) && return DataFrames.DataFrame()
    return DataFrames.DataFrame(data, keep_names)
end

"""
Re-categorizes data according to an aggregation dictionary
* makes no guarantee of complete data collection *

# Example

```julia
aggregation = PG.make_fuel_dictionary(results_uc.system)
PG.categorize_data(gen_uc.data, aggregation)
```

"""
function categorize_data(
    data::Dict{Symbol, DataFrames.DataFrame},
    aggregation::Dict;
    curtailment = true,
    slacks = true,
)
    category_dataframes = Dict{String, DataFrames.DataFrame}()
    var_types = Dict(zip(last.(split.(string.(keys(data)), "_")), keys(data)))
    for (category, list) in aggregation
        category_df = DataFrames.DataFrame()
        for tuple in list
            if haskey(var_types, tuple[1])
                category_data = data[var_types[tuple[1]]]
                colname =
                    typeof(names(category_data)[1]) == String ? "$(tuple[2])" :
                    Symbol(tuple[2])
                DataFrames.insertcols!(
                    category_df,
                    (colname => category_data[:, colname]),
                    makeunique = true,
                )
            end
        end
        category_dataframes[string(category)] = category_df
    end
    if curtailment && haskey(data, :Curtailment)
        category_dataframes["Curtailment"] = data[:Curtailment]
    end
    for (slack, slack_name) in SLACKVARS
        if slacks && haskey(data, slack)
            category_dataframes[slack_name] = data[slack]
        end
    end

    return category_dataframes
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Kriging(var₁=>param₁, var₂=>param₂, ...)

A polyalgorithm Kriging estimation solver.

Each pair `var=>param` specifies the `KrigingParam` `param`
for the Kriging variable `var`. In order to avoid boilerplate
code, the constructor expects pairs of `Symbol` and `NamedTuple`
instead.

## Parameters

* `variogram` - Variogram model (default to `GaussianVariogram()`)
* `mean`      - Simple Kriging mean
* `degree`    - Universal Kriging degree
* `drifts`    - External Drift Kriging drift functions

Latter options override former options. For example, by specifying
`drifts`, the user is telling the algorithm to ignore `degree` and
`mean`. If no option is specified, Ordinary Kriging is used by
default with the `variogram` only.

* `minneighbors` - Minimum number of neighbors (default to `1`)
* `maxneighbors` - Maximum number of neighbors (default to `10`)
* `neighborhood` - Search neighborhood (default to `nothing`)
* `distance`     - Distance used to find nearest neighbors (default to `Euclidean()`)

The `maxneighbors` option can be used to perform approximate Kriging
with a subset of data points per estimation location. If `maxneighbors`
is set to `nothing`, global Kriging is performed instead. Two neighborhood
search methods are available depending on the value of `neighborhood`:

* If a `neighborhood` is provided, local Kriging is performed by sliding
  the `neighborhood` in the domain.

* If `neighborhood` is not provided, the Kriging system is built using
  `maxneighbors` nearest neighbors according to a `distance`.

## Examples

Solve the variable `:var₁` with Simple Kriging by specifying
the `mean`, and the variable `:var₂` with Universal Kriging
by specifying the `degree` and the `variogram` model.

```julia
julia> Kriging(
  :var₁ => (mean=1.,),
  :var₂ => (degree=1, variogram=SphericalVariogram(range=20.))
)
```

Solve all variables of the problem with the default parameters
(i.e. Ordinary Kriging with unit Gaussian variogram):

```julia
julia> Kriging()
```
"""
@estimsolver Kriging begin
  @param variogram = GaussianVariogram()
  @param mean = nothing
  @param degree = nothing
  @param drifts = nothing
  @param minneighbors = 1
  @param maxneighbors = 10
  @param neighborhood = nothing
  @param distance = Euclidean()
end

function preprocess(problem::EstimationProblem, solver::Kriging)
  # retrieve problem info
  pdomain = domain(problem)
  pdata   = data(problem)
  ndata   = nelements(pdata)

  # result of preprocessing
  preproc = Dict{Symbol,NamedTuple}()

  for covars in covariables(problem, solver)
    for var in covars.names
      # get user parameters
      varparams = covars.params[(var,)]

      # determine which Kriging variant to use
      if varparams.drifts ≠ nothing
        estimator = ExternalDriftKriging(varparams.variogram, varparams.drifts)
      elseif varparams.degree ≠ nothing
        estimator = UniversalKriging(varparams.variogram, varparams.degree, embeddim(pdomain))
      elseif varparams.mean ≠ nothing
        estimator = SimpleKriging(varparams.variogram, varparams.mean)
      else
        estimator = OrdinaryKriging(varparams.variogram)
      end

      # determine minimum/maximum number of neighbors
      minneighbors = varparams.minneighbors
      maxneighbors = varparams.maxneighbors

      # determine neighborhood search method
      if varparams.maxneighbors ≠ nothing
        # upper bound in maxneighbors
        maxneighbors > ndata && (maxneighbors = ndata)
        if varparams.neighborhood ≠ nothing
          # local search with a neighborhood
          neigh = varparams.neighborhood

          if neigh isa MetricBall
            bsearcher = KBallSearch(pdata, maxneighbors, neigh)
          else
            searcher  = BallSearch(pdata, neigh)
            bsearcher = BoundedSearch(searcher, maxneighbors)
          end
        else
          # nearest neighbor search with a distance
          distance  = varparams.distance
          bsearcher = KNearestSearch(pdata, maxneighbors, metric=distance)
        end
      else
        # use all data points as neighbors
        bsearcher = nothing
      end

      # save preprocessed input
      preproc[var] = (estimator=estimator,
                      minneighbors=minneighbors,
                      maxneighbors=maxneighbors,
                      bsearcher=bsearcher)
    end
  end

  preproc
end

function solve(problem::EstimationProblem, solver::Kriging)
  # preprocess user input
  preproc = preprocess(problem, solver)

  # results for each variable
  μs = []; σs = []
  for var in name.(variables(problem))
    if preproc[var].maxneighbors ≠ nothing
      # perform Kriging with reduced number of neighbors
      varμ, varσ = solve_approx(problem, var, preproc)
    else
      # perform Kriging with all data points as neighbors
      varμ, varσ = solve_exact(problem, var, preproc)
    end

    push!(μs, var => varμ)
    push!(σs, Symbol(var,"_variance") => varσ)
  end

  georef((; μs..., σs...), domain(problem))
end

function solve_approx(problem::EstimationProblem, var::Symbol, preproc)
    # retrieve problem info
    pdata = data(problem)
    pdomain = domain(problem)

    mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))

    # unpack preprocessed parameters
    estimator, minneighbors, maxneighbors, bsearcher = preproc[var]

    # determine value type
    V = mactypeof[var]

    # pre-allocate memory for result
    varμ = Vector{V}(undef, nelements(pdomain))
    varσ = Vector{V}(undef, nelements(pdomain))

    # pre-allocate memory for neighbors
    neighbors = Vector{Int}(undef, maxneighbors)

    # estimation loop
    for location in traverse(pdomain, LinearPath())
      pₒ = centroid(pdomain, location)

      # find neighbors with previously estimated values
      nneigh = search!(neighbors, pₒ, bsearcher)

      # skip location in there are too few neighbors
      if nneigh < minneighbors
        varμ[location] = NaN
        varσ[location] = NaN
      else
        # final set of neighbors
        nview = view(neighbors, 1:nneigh)

        # view neighborhood with data
        𝒟 = view(pdata, nview)

        # fit estimator to data
        krig = fit(estimator, 𝒟)

        # retrieve element at location
        uₒ = pdomain[location]

        # save mean and variance
        μ, σ² = predict(krig, var, uₒ)

        varμ[location] = μ
        varσ[location] = σ²
      end
    end

    varμ, varσ
end

function solve_exact(problem::EstimationProblem, var::Symbol, preproc)
    # retrieve problem info
    pdata = data(problem)
    pdomain = domain(problem)

    mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))

    # unpack preprocessed parameters
    estimator, minneighbors, maxneighbors, bsearcher = preproc[var]

    # determine value type
    V = mactypeof[var]

    # pre-allocate memory for result
    varμ = Vector{V}(undef, nelements(pdomain))
    varσ = Vector{V}(undef, nelements(pdomain))

    # retrieve non-missing data
    locs = findall(!ismissing, pdata[var])
    𝒟 = view(pdata, locs)

    # fit estimator once
    krig = fit(estimator, 𝒟)

    # predict at all locations
    for location in traverse(pdomain, LinearPath())
      uₒ = pdomain[location]

      μ, σ² = predict(krig, var, uₒ)

      varμ[location] = μ
      varσ[location] = σ²
    end

    varμ, varσ
end


text:
@testset "Test How to State NLP" begin
  ###############################################################################
  #
  # The data used through the algorithmic process in the Stopping framework
  # are stored in a State.
  # We illustrate here the NLPAtX which is a specialization of the State for
  # non-linear programming.
  #
  ###############################################################################
  #using Test, NLPModels, Stopping

  include("../test-stopping/rosenbrock.jl")
  #Formulate the problem with NLPModels
  x0 = ones(6)
  y0 = ones(1)
  c(x) = [x[1] - x[2]]
  lcon = [0.0]
  ucon = [0.0]

  #We can create a NLPAtX for constrained optimization.
  #Here we provide y0 = [1.0]
  #Note that the default value is [0.0]
  #Warning: see https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/autodiff_model.jl
  #for the proper way of defining an ADNLPModel
  meta = NLPModelMeta(
    6,
    x0 = x0,
    lvar = zeros(6),
    uvar = Inf * ones(6),
    ncon = 1,
    y0 = [0.0],
    lcon = lcon,
    ucon = ucon,
  )
  nlp = ADNLPModel(meta, Counters(), ADNLPModels.ForwardDiffAD(6, 1, rosenbrock, x0), rosenbrock, c)
  #We can create a NLPAtX for bounds-constrained optimization:
  meta = NLPModelMeta(6, x0 = x0, lvar = zeros(6), uvar = Inf * ones(6))
  nlp2 =
    ADNLPModel(meta, Counters(), ADNLPModels.ForwardDiffAD(6, rosenbrock, x0), rosenbrock, x -> [])
  #We can create a NLPAtX for unconstrained optimization:
  nlp3 = ADNLPModel(x -> rosenbrock(x), x0)

  ###############################################################################
  #I. Initialize a NLPAtX:
  #
  #There are two main constructor for the States.
  #The unconstrained:
  state_unc = NLPAtX(x0)
  #The constrained:
  state_con = NLPAtX(x0, y0)

  #By default, all the values in the State are set to nothing except x and lambda
  #In the unconstrained case lambda is a vector of length 0
  @test !(state_unc.lambda == nothing)
  #From the default initialization, all the other entries are void:
  @test state_unc.mu == [] && state_con.mu == []
  @test isnan(state_unc.fx) && isnan(state_con.fx)

  #Note that the constructor proceeds to a size checking on gx, Hx, mu, cx, Jx.
  #It returns an error if this test fails.
  try
    NLPAtX(x0, Jx = ones(1, 1))
    @test false
  catch
    #printstyled("NLPAtX(x0, Jx = ones(1,1)) is invalid as length(lambda)=0\n")
    @test true
  end

  ###############################################################################
  #II. Update the entries
  #
  #At the creation of a NLPAtX, keyword arguments populate the state:
  state_bnd = NLPAtX(x0, mu = zeros(6))
  @test state_bnd.mu == zeros(6) #initialize multipliers with bounds constraints

  #The NLPAtX has two functions: update! and reinit!
  #The update! has the same behavior as in the GenericState:
  update!(state_bnd, fx = 1.0, blah = 1) #update! ignores unnecessary keywords
  @test state_bnd.mu == zeros(6) && state_bnd.fx == 1.0 && state_bnd.x == x0

  #reinit! by default reuse x and lambda and reset all the entries at their
  #default values (void or empty Counters):
  reinit!(state_bnd, mu = ones(6))
  @test state_bnd.mu == ones(6) && isnan(state_bnd.fx)
  @test state_bnd.x == x0 && state_bnd.lambda == zeros(0)
  #However, we can specify both entries
  reinit!(state_bnd, 2 * ones(6), zeros(0))
  @test state_bnd.x == 2 * ones(6) && state_bnd.lambda == zeros(0)
  @test state_bnd.mu == []

  ###############################################################################
  #III. Domain Error
  #Similar to the GenericState we can use _domain_check to verify there are no NaN
  @test Stopping._domain_check(state_bnd) == false
  update!(state_bnd, mu = [NaN])
  @test Stopping._domain_check(state_bnd) == true

  ###############################################################################
  #IV. Use the NLPAtX
  #
  #For algorithmic use, it might be conveninent to fill in all the entries of then
  #State. In this case, we can use the Stopping:
  stop = NLPStopping(nlp, state_unc, optimality_check = (x, y) -> unconstrained_check(x, y))
  #Note that the fill_in! can receive known informations via keywords.
  #If we don't want to store the hessian matrix, we turn the keyword
  #matrix_info as false.
  fill_in!(stop, x0, matrix_info = false)

  #printstyled("Hx has not been updated: ",stop.current_state.Hx == nothing,"\n")
  @test stop.current_state.Hx == zeros(0, 0)

  # We can now use the updated step in the algorithmic procedure
  @test start!(stop)
end


text:
@testset "atari" begin
    atari_env_names = ReinforcementLearningEnvironments.list_atari_rom_names()
    atari_env_names = filter(x -> x ∉ ["pacman", "surround"], atari_env_names)

    @testset "basic tests" begin
        for x in atari_env_names
            env = AtariEnv(x)
            # RLBase.test_interfaces!(env)
            RLBase.test_runnable!(env)
        end
    end

    @testset "seed" begin
        env = AtariEnv(; name = "pong", seed = 456)
        old_states = []
        actions = [rand(action_space(env)) for i in 1:10, j in 1:100]

        for i in 1:10
            for j in 1:100
                env(actions[i, j])
                push!(old_states, copy(state(env)))
            end
            reset!(env)
        end

        env = AtariEnv(; name = "pong", seed = 456)
        new_states = []
        for i in 1:10
            for j in 1:100
                env(actions[i, j])
                push!(new_states, copy(state(env)))
            end
            reset!(env)
        end

        @test old_states == new_states
    end

    @testset "frame_skip" begin
        env = AtariEnv(; name = "pong", frame_skip = 4, seed = 456)
        states = []
        actions = [rand(action_space(env)) for _ in 1:100]

        for i in 1:100
            env(actions[i])
            push!(states, copy(state(env)))
        end

        env = AtariEnv(; name = "pong", frame_skip = 1, seed = 456)
        for i in 1:100
            env(actions[i])
            env(actions[i])
            env(actions[i])
            s1 = copy(state(env))
            env(actions[i])
            s2 = copy(state(env))
            @test states[i] == max.(s1, s2)
        end
    end

    @testset "repeat_action_probability" begin
        env = AtariEnv(; name = "pong", repeat_action_probability = 1.0, seed = 456)
        states = []
        actions = [rand(action_space(env)) for _ in 1:100]
        for i in 1:100
            env(actions[i])
            push!(states, copy(state(env)))
        end

        env = AtariEnv(; name = "pong", repeat_action_probability = 1.0, seed = 456)
        for i in 1:100
            env(actions[1])
            @test states[i] == state(env)
        end
    end

    @testset "max_num_frames_per_episode" begin
        for i in 1:10
            env = AtariEnv(; name = "pong", max_num_frames_per_episode = i, seed = 456)
            for _ in 1:i
                env(1)
            end
            @test true == is_terminated(env)
        end
    end
end


text:
"""
    observe_from_history!(model::DriverModel, roadway::Roadway, trajdata::Vector{<:EntityScene}, egoid, start::Int, stop::Int)

Given a prerecorded trajectory `trajdata`, run the observe function of a driver model for the scenes between `start` and `stop` for the vehicle of id `egoid`.
The ego vehicle does not take any actions, it just observe the scenes,
"""
function observe_from_history!(
        model::DriverModel,
        roadway::Roadway,
        trajdata::Vector{<:EntityScene},
        egoid,
        start::Int = 1,
        stop::Int = length(trajdata))
    reset_hidden_state!(model)

    for i=start:stop
        observe!(model, trajdata[i], roadway, egoid)
    end

    return model
end

"""
    maximum_entities(trajdata::Vector{<:EntityScene})
return the maximum number of entities present in a given scene of the trajectory.
"""
function maximum_entities(trajdata::Vector{<:EntityScene})
    return maximum(capacity, trajdata)
end

"""
    simulate_from_history(model::DriverModel, roadway::Roadway, trajdata::Vector{Scene{E}}, egoid, timestep::Float64,
                          start::Int = 1, stop::Int = length(trajdata);
                          rng::AbstractRNG = Random.GLOBAL_RNG) where {E<:Entity}

Replay the given trajectory except for the entity `egoid` which follows the given driver model.
See `simulate_from_history!` if you want to provide a container for the results or log actions.
"""
function simulate_from_history(
        model::DriverModel,
        roadway::Roadway,
        trajdata::Vector{Scene{E}},
        egoid,
        timestep::Float64,
        start::Int = 1,
        stop::Int = length(trajdata);
        rng::AbstractRNG = Random.GLOBAL_RNG
        ) where {E<:Entity}
    scenes = [Scene(E, maximum_entities(trajdata)) for i=1:(stop - start + 1)]
    n = simulate_from_history!(model, roadway, trajdata, egoid, timestep,
                               start, stop, scenes,
                               rng=rng)
    return scenes[1:(n+1)]
end

"""
    simulate_from_history!(model::DriverModel, roadway::Roadway, trajdata::Vector{Scene{E}}, egoid, timestep::Float64,
                          start::Int, stop::Int, scenes::Vector{Scene{E}};
                          actions::Union{Nothing, Vector{Scene{A}}} = nothing, rng::AbstractRNG = Random.GLOBAL_RNG) where {E<:Entity}

Replay the given trajectory except for the entity `egoid` which follows the given driver model.
The resulting trajectory is stored in `scenes`, the actions of the ego vehicle are stored in `actions`.

"""
function simulate_from_history!(
    model::DriverModel,
    roadway::Roadway,
    trajdata::Vector{Scene{E}},
    egoid,
    timestep::Float64,
    start::Int,
    stop::Int,
    scenes::Vector{Scene{E}};
    actions::Union{Nothing, Vector{Scene{A}}} = nothing,
    rng::AbstractRNG = Random.GLOBAL_RNG
    ) where {E<:Entity, A<:EntityAction}

    # run model (unsure why it is needed, it was in the old code )
    observe_from_history!(model, roadway, trajdata, egoid, start, stop)

    copyto!(scenes[1], trajdata[start])
    for tick=1:(stop - start)

        empty!(scenes[tick + 1])
        if (actions !== nothing) empty!(actions[tick]) end

        ego = get_by_id(scenes[tick], egoid)
        observe!(model, scenes[tick], roadway, egoid)
        a = rand(rng, model)

        ego_state_p = propagate(ego, a, roadway, timestep)

        copyto!(scenes[tick+1], trajdata[start+tick])
        egoind = findfirst(egoid, scenes[tick+1])
        scenes[tick+1][egoind] = Entity(ego_state_p, ego.def, egoid)

        if (actions !== nothing) push!(actions[tick], EntityAction(a, egoid)) end

    end
    return (stop - start)
end