(optionally) Sanity-check user functions about ipopt.jl HOT 4 CLOSED

function eval_f(x::Vector{Float64})
    # !!! The x[5] is out-of-bounds !!!
    return x[1] * x[4] * (x[1] + x[2] + x[3]) + x[5]
end

function eval_g(x::Vector{Float64}, g::Vector{Float64})
    g[1] = x[1]   * x[2]   * x[3]   * x[4]
    g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
    grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
    grad_f[2] = x[1] * x[4]
    grad_f[3] = x[1] * x[4] + 1
    grad_f[4] = x[1] * (x[1] + x[2] + x[3])
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
    if mode == :Structure
        # Constraint (row) 1
        rows[1] = 1; cols[1] = 1
        rows[2] = 1; cols[2] = 2
        rows[3] = 1; cols[3] = 3
        rows[4] = 1; cols[4] = 4
        # Constraint (row) 2
        rows[5] = 2; cols[5] = 1
        rows[6] = 2; cols[6] = 2
        rows[7] = 2; cols[7] = 3
        rows[8] = 2; cols[8] = 4
    else
        # Constraint (row) 1
        values[1] = x[2]*x[3]*x[4]  # 1,1
        values[2] = x[1]*x[3]*x[4]  # 1,2
        values[3] = x[1]*x[2]*x[4]  # 1,3
        values[4] = x[1]*x[2]*x[3]  # 1,4
        # Constraint (row) 2
        values[5] = 2*x[1]  # 2,1
        values[6] = 2*x[2]  # 2,2
        values[7] = 2*x[3]  # 2,3
        values[8] = 2*x[4]  # 2,4
    end
end

function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
    if mode == :Structure
        # Symmetric matrix, fill the lower left triangle only
        idx = 1
        for row = 1:4
            for col = 1:row
            rows[idx] = row
            cols[idx] = col
            idx += 1
            end
        end
    else
        # Again, only lower left triangle
        # Objective
        values[1] = obj_factor * (2*x[4])  # 1,1
        values[2] = obj_factor * (  x[4])  # 2,1
        values[3] = 0                      # 2,2
        values[4] = obj_factor * (  x[4])  # 3,1
        values[5] = 0                      # 3,2
        values[6] = 0                      # 3,3
        values[7] = obj_factor * (2*x[1] + x[2] + x[3])  # 4,1
        values[8] = obj_factor * (  x[1])  # 4,2
        values[9] = obj_factor * (  x[1])  # 4,3
        values[10] = 0                     # 4,4

        # First constraint
        values[2] += lambda[1] * (x[3] * x[4])  # 2,1
        values[4] += lambda[1] * (x[2] * x[4])  # 3,1
        values[5] += lambda[1] * (x[1] * x[4])  # 3,2
        values[7] += lambda[1] * (x[2] * x[3])  # 4,1
        values[8] += lambda[1] * (x[1] * x[3])  # 4,2
        values[9] += lambda[1] * (x[1] * x[2])  # 4,3

        # Second constraint
        values[1]  += lambda[2] * 2  # 1,1
        values[3]  += lambda[2] * 2  # 2,2
        values[6]  += lambda[2] * 2  # 3,3
        values[10] += lambda[2] * 2  # 4,4
    end
end

prob = createProblem(
    4, 
    [1.0, 1.0, 1.0, 1.0], 
    [5.0, 5.0, 5.0, 5.0], 
    2, 
    [25.0, 40.0], 
    [Inf, 40.0], 
    8, 
    10,
    eval_f, 
    eval_g,
    eval_grad_f, 
    eval_jac_g, 
    eval_h,
)

prob.x = [1.0, 5.0, 5.0, 1.0]
solvestat = solveProblem(prob)

julia> solvestat = solveProblem(prob)
This is Ipopt version 3.13.2, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        4
Number of nonzeros in inequality constraint Jacobian.:        4
Number of nonzeros in Lagrangian Hessian.............:       10

Total number of variables............................:        4
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        4
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

ERROR: BoundsError: attempt to access 4-element Array{Float64,1} at index [5]
Stacktrace:
 [1] getindex at ./array.jl:809 [inlined]
 [2] eval_f(::Array{Float64,1}) at ./REPL[38]:2
 [3] eval_f_wrapper(::Int32, ::Ptr{Float64}, ::Int32, ::Ptr{Float64}, ::Ptr{Nothing}) at /Users/oscar/.julia/dev/Ipopt/src/Ipopt.jl:143
 [4] solveProblem(::IpoptProblem) at /Users/oscar/.julia/dev/Ipopt/src/Ipopt.jl:361
 [5] top-level scope at REPL[45]:1