search

chkwon / PATHSolver.jl

provides a Julia wrapper for the PATH Solver for solving mixed complementarity problems
http://pages.cs.wisc.edu/%7Eferris/path.html
MIT License
50 stars 15 forks source link

solve_mcp(): NNZ or NNZ+1 #77

Closed sdirkse67 closed 1 year ago

sdirkse67 commented 1 year ago

In the code for solve_mcp() is this line:

nnz = Int(dnnz + 1)

This essentially adds an extra/unused Cint or Cdouble to the vectors of length nnz that are passed to J(). It doesn't really matter if J() is always resetting nnz to the actual count, but it's a little confusing or offputting to call solve_mcp(...,nnz=8) but get nnz=9 in the call to J().

There might have been a purpose to this bumping of NNZ with PATH 4.x, but with PATH 5.x I don't see any need for this.

Of course, if you change it, anybody who was running with solve_mcp(...,nnz=actualNNZ-1) is going to be unhappy. And there is surely at least one such user . . .

odow commented 1 year ago

Hmm. I don't remember why I did this... Let me take a look.

We set the actual nnz I the callback:

https://github.com/chkwon/PATHSolver.jl/blob/606c5f271a6f3ac3d00ce3660ec83a93cc87841a/src/C_API.jl#L243

And there is surely at least one such user . . .

I'm not actually sure how many people use PATH directly from Julia. You'd think they would have opened an issue, because this is arguably a bug.

JuliaHub tells me that there are only 67 unique users in the last month. But that includes the holidays so its likely a little higher:

image

I also audited every public instance of PATHSolver.solve_mcp, https://juliahub.com/ui/Search?q=PATHSolver.solve_mcp&type=code, and I don't see anyone futzing with the nnz count so I think we're okay to change.

chkwon commented 1 year ago

@sdirkse67 I have a question for nnz.

Complementarity.jl is one of those that use solve_mcp(). I wasn't sure how to set nnz by the way. Please see here: https://github.com/chkwon/Complementarity.jl/blob/v0.9.0/src/mcp.jl#L220-L228

In the above, I put a "sufficiently large" number for nnz in a strange way. I remember that when I entered a tight number, it didn't work properly. What would be the best way to set nnz? Should I just use the default value as set by @odow (which is length(lb)^2)?

odow commented 1 year ago

I remember that when I entered a tight number, it didn't work properly.

This might have been PATH 4.x? Or it might have been due to the Julia<->C conversion. Since Julia is 1-indexed, we often need to add an extra element to arrays when passing to 0-indexed C programs which use 1-based indexing (GLPK is a killer for this).

I obviously added the +1 for a similar reason, although I don't have notes on why.

I updated a PR here: https://github.com/chkwon/PATHSolver.jl/pull/78

sdirkse67 commented 1 year ago

Looking back at some of the earliest PATH driver examples I see a similar, seemingly arbitrary bump in NNZ, i.e. NNZ = NNZ + 1. So it's quite possible @odow got this from these examples. It's also possible that using a tight/proper NNZ count was exposing some bugs that are hidden when using an over-estimate, again as @odow mentions. But that should not be necessary. Assuming that Complementarity.jl provides a Jacobian with a constant structure, I would just pass the actual count of the nonzeros in this Jacobian to solve_mcp().

sdirkse67 commented 1 year ago

Recently @lassepe sent me an example driver that passes NNZ-1 to solve_mcp() so that the expected NNZ value appears in J(). But it should be no issue for him to pass NNZ instead of NNZ-1. It should be easy enough for him to do this independently of which version/revision of PATHSolver.jl he is using.

odow commented 1 year ago

Looking back at some of the earliest PATH driver examples I see a similar, seemingly arbitrary bump in NNZ, i.e. NNZ = NNZ + 1. So it's quite possible @odow got this from these examples

Yes. This could sound right. I think I wrote the original wrapper while visiting Michael in Madison, and I was using a mix of the PATH 4 examples online and a custom PATH 5 binary that he sent me before you released it officially.

Recently @lassepe sent me an example driver that passes NNZ-1 to solve_mcp() so that the expected NNZ value appears in J()

Ah. I had found only this usage: https://github.com/lassepe/DifferentiableTrajectoryOptimization.jl/blob/77a3b6e639b5f3ee7a39b30012c55e5d5514d0f5/src/mcp_solver.jl#L100

@lassepe want to chime in on this "bug" fix? https://github.com/chkwon/PATHSolver.jl/pull/78

lassepe commented 1 year ago

Yes, in another code snippet that I shared with @sdirkse67 for the purpose of chasing down #70 I noticed this discrepancy and subtracted 1 manually. I forgot about it and thus never filed an issue. I'd be happy to have it fixed here and adjust the code on my end. I don't think many people will be affected by this since passing nnz - 1 seems uncommon and is not strictly necessary (I just did it in hopes of finding #70). Even if someone was bitten by this bugfix, it's easy enough to set a compat bound in their Project.toml.

Also, according to JuliaHub, PATHSolver only has 4 direct (registered) dependents

https://juliahub.com/ui/Packages/PATHSolver/xyGrF/1.3.0?page=2

TensorGames.jl and DifferentiableTrajectoryOptimization.jl are co-maintained by me and will be unaffected by this change. Complementarity.jl is maintained by you and Dynare.jl uses the C-interface but does not correct nnz manually and thus should be unaffected.