Update mcp.jl

[solliolli]

Came across a weird bug when helping a friend who has very little experience on Julia or complementarity modelling. If for some reason, the last variable defined in the model does not appear in any of the complementarity expressions/mappings, that will cause a bug with the sparse Jacobian where the dimensions of that matrix become smaller than the number of constraints. This will then lead to a BoundsError later in the code. A similar bug would happen if the last complementarity constraint had a constant value for the expression/mapping. My proposed fix is to force the correct dimensions for the sparse Jacobian matrix.

Minimal working example modified from the README example:

using Complementarity, JuMP

m = MCPModel()

@variable(m, x >= 0)
@variable(m, y >= 0)

@mapping(m, F1, x+2)
@mapping(m, F2, 2)

@complementarity(m, F1, x)
@complementarity(m, F2, y)

status = solveMCP(m)

@show result_value(x)
@show result_value(y)

[chkwon]

Thanks!

[andyYaoR]

Sorry I am new to Julia.
The above example does not run and raises boundserror, and I have the same error with my own model. It seems the issue is related to variable not appear in the constraint. Can you please take a look?
Many thanks.

[chkwon]

@andyYaoR Thanks.

The above example produces the following error:

Path 5.0.03 (Fri Jun 26 09:58:07 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris
ERROR: BoundsError: attempt to access 2-element Vector{Int64} at index [3]
Stacktrace:
  [1] getindex
    @ ./array.jl:924 [inlined]
  [2] (::Complementarity.var"#jacobian_callback#9"{Complementarity.var"#j_eval#8"{MathOptInterface.Nonlinear.Evaluator{MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator}}})(n::Int32, nnz::Int32, z::Vector{Float64}, col_start::Vector{Int32}, col_len::Vector{Int32}, row::Vector{Int32}, data::Vector{Float64})
    @ Complementarity ~/.julia/packages/Complementarity/NbsQM/src/mcp.jl:178
  [3] _c_jacobian_evaluation(id_ptr::Ptr{Nothing}, n::Int32, x_ptr::Ptr{Float64}, wantf::Int32, f_ptr::Ptr{Float64}, nnz_ptr::Ptr{Int32}, col_ptr::Ptr{Int32}, len_ptr::Ptr{Int32}, row_ptr::Ptr{Int32}, data_ptr::Ptr{Float64})
    @ PATHSolver ~/.julia/packages/PATHSolver/09Pix/src/C_API.jl:242
  [4] c_api_Path_Solve
    @ ~/.julia/packages/PATHSolver/09Pix/src/C_API.jl:586 [inlined]
  [5] solve_mcp(F::Complementarity.var"#function_callback#7"{MathOptInterface.Nonlinear.Evaluator{MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator}}, J::Complementarity.var"#jacobian_callback#9"{Complementarity.var"#j_eval#8"{MathOptInterface.Nonlinear.Evaluator{MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator}}}, lb::Vector{Float64}, ub::Vector{Float64}, z::Vector{Float64}; nnz::Int64, variable_names::Vector{String}, constraint_names::Vector{String}, silent::Bool, generate_output::Int64, use_start::Bool, use_basics::Bool, jacobian_structure_constant::Bool, jacobian_data_contiguous::Bool, jacobian_linear_elements::Vector{Int64}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PATHSolver ~/.julia/packages/PATHSolver/09Pix/src/C_API.jl:813
  [6] _solve_path!(m::Model; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Complementarity ~/.julia/packages/Complementarity/NbsQM/src/mcp.jl:232
  [7] _solve_path!
    @ ~/.julia/packages/Complementarity/NbsQM/src/mcp.jl:128 [inlined]
  [8] solveMCP(m::Model; solver::Symbol, method::Symbol, linear::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Complementarity ~/.julia/packages/Complementarity/NbsQM/src/mcp.jl:102
  [9] solveMCP(m::Model)
    @ Complementarity ~/.julia/packages/Complementarity/NbsQM/src/mcp.jl:96
 [10] top-level scope
    @ REPL[9]:1

@solliolli any idea?

[solliolli]

Hi! Looking at it quickly, I think commit 4d0be88 might have undone my fix to the sparse Jacobian representation. The problem is, I'm not very familiar with this side of MOI, but if I understand correctly, the line

Complementarity.jl/src/mcp.jl

Line 154 in cd4fa31

nc = length(J)

now sets the dimensions of the Jacobian matrix to be equal to the number of nonzero elements in the Jacobian. Once again, that breaks things in the unlikely scenario that our Jacobian is so extremely sparse that, for example, a 2x2 matrix has one nonzero element.

[solliolli]

The value nc should be the number of constraints in the model, and to fix the issue, I would need to figure out how to extract that number. I'm a bit busy these days, it might take a while for me to familiarize myself with these evaluators, but if anyone feels like they want to try fixing it, please go ahead.

[solliolli]

It seems like replacing the line I mentioned with nc = length(d.ordered_constraints) might fix the problem, but I did not test this further than checking that it removes the error you got.