Address non-feasibility in solveset - current code breaks solveset #19189
Conversation
|
✅ Hi, I am the SymPy bot (v160). I'm here to help you write a release notes entry. Please read the guide on how to write release notes. Your release notes are in good order. Here is what the release notes will look like:
This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7. Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. Click here to see the pull request description that was parsed.
Update The release notes on the wiki have been updated. |
Codecov Report
@@ Coverage Diff @@
## master #19189 +/- ##
=============================================
- Coverage 75.666% 75.613% -0.053%
=============================================
Files 651 652 +1
Lines 169336 169572 +236
Branches 39970 40010 +40
=============================================
+ Hits 128131 128220 +89
- Misses 35605 35736 +131
- Partials 5600 5616 +16 |
|
Thanks for this. Can you add an example that demonstrates the bug to the tests? Maybe it would be simpler just to check if |
|
Thank you for the quick response. I have changed the feasible flag into a None assignment and check, as per your suggestion. I've also added the above example as a test, and confirmed that it fails on the master branch and passes with the new code. Thanks very much |
|
Is the given solution correct? If so it would be good to assert that we get that solution in the tests e.g.: You can use The tests have failed because of the use of unicode characters. It is possible to make them pass by adding an encoding header and whitelisting this file. I think it would be better in this case just to use ASCII for the symbol names because it makes it easier for people to edit the file. |
|
Thanks - I can use ascii characters. Actually, the solutions are incorrect. Note for example the equation However, I can see no way this change could possibly have introduced this bug. I've observed similar issues before this change as well. How should this be handled? As you say, it seems incorrect to assert the solution is something it is not. Thanks very much. |
I think the best way is to add a test asserting that the solution is not what is currently returned. That test will obviously fail but you can mark the test as XFAIL. That way if something changes in future and the test begins to pass someone will be alerted to the XPASS and can investigate (the solution might still be incorrect but someone can at least look). You could also open a github issue to describe the incorrect answer and then put a url for the issue in a code comment of the XFAIL test so anyone looking at it can see the discussion. |
|
If you can come up with an example where sympy actually gets the correct result then that's also good :) I see that for the example above There can be differences in the way that I have discussed ideas about this in #16861 |
|
The Groebner basis for the system has each individual unknown as a generator: In [18]: system
Out[18]:
[-λ0⋅ψ00 - λ1⋅ψ01 + μ0 + ψ00 + ψ01, -λ0⋅ψ10 - λ1⋅ψ11 + μ1, -λ2⋅ψ00 - λ3⋅ψ01 + ψ00 + ψ01, -λ2⋅ψ10 - λ3⋅ψ11 + μ3, -λ0⋅ϕ00 - λ2⋅ϕ10 + ϕ00 + ϕ10,
-λ1⋅ϕ00 - λ3⋅ϕ10 + ϕ00 + ϕ10, -λ0⋅ϕ01 - λ2⋅ϕ11, -λ1⋅ϕ01 - λ3⋅ϕ11, -α00 + ψ00⋅ϕ00 + ψ10⋅ϕ01, -α01 + ψ01⋅ϕ00 + ψ11⋅ϕ01, -α10 + ψ00⋅ϕ10 + ψ10⋅ϕ
11, -α11 + ψ01⋅ϕ10 + ψ11⋅ϕ11, -μ0⋅ϕ00, -μ1⋅ϕ01, -μ2⋅ϕ10, -μ3⋅ϕ11, -μ4⋅ψ00, -μ5⋅ψ01, -μ6⋅ψ10, -μ7⋅ψ11, μ2, μ4, μ5, μ6, μ7]
In [19]: solvefor
Out[19]: [ϕ00, ϕ01, ϕ10, ϕ11, ψ00, ψ01, ψ10, ψ11, μ0, μ1, μ3, λ0, λ1, λ2, λ3]
In [20]: groebner(system, solvefor)
Out[20]:
GroebnerBasis([1], ϕ00, ϕ01, ϕ10, ϕ11, ψ00, ψ01, ψ10, ψ11, μ0, μ1, μ3, λ0, λ1, λ2, λ3, domain=ℤ[α00, α01, α10, α11, μ2, μ4, μ5, μ6, μ7], orde
r=lex)I think this implies that the only possible solution would be one in which every unknown is zero in which case the system becomes: In [21]: rep = {x: 0 for x in solvefor}
In [22]: for e in system: pprint(e.subs(rep))
0
0
0
0
0
0
0
0
-α00
-α01
-α10
-α11
0
0
0
0
0
0
0
0
μ2
μ4
μ5
μ6
μ7So maybe the only possible solution is the trivial solution and it only exists if these parameters are all zero. |
|
Thank you, that is very interesting. If I follow, this implies that were I to solve for the system all the variables used in it, I should get a solution where all variables turn out to be zeros. This is not quite what I see. With nonlinsolve I get the following error: Using solve, with dict=True, I obtain the following solutions. This leaves me unsure of which assertion to make. For completeness, the code used to obtain these results is below, with the only change between the first and second being solve -> nonlinsolve: import sympy as sy
symbols = sy.symbols(
'α00, α01, α10, α11, λ0, λ1, λ2, λ3, μ0, μ1, μ2, μ3, μ4, μ5, μ6, μ7, ψ00, ψ01, ψ10, ψ11, ϕ00, ϕ01, ϕ10, ϕ11'
)
α00, α01, α10, α11, λ0, λ1, λ2, λ3, μ0, μ1, μ2, μ3, μ4, μ5, μ6, μ7, ψ00, ψ01, ψ10, ψ11, ϕ00, ϕ01, ϕ10, ϕ11 = symbols
system = [
-λ0 * ψ00 - λ1 * ψ01 + μ0 + ψ00 + ψ01,
-λ0 * ψ10 - λ1 * ψ11 + μ1,
-λ2 * ψ00 - λ3 * ψ01 + ψ00 + ψ01,
-λ2 * ψ10 - λ3 * ψ11 + μ3,
-λ0 * ϕ00 - λ2 * ϕ10 + ϕ00 + ϕ10,
-λ1 * ϕ00 - λ3 * ϕ10 + ϕ00 + ϕ10,
-λ0 * ϕ01 - λ2 * ϕ11,
-λ1 * ϕ01 - λ3 * ϕ11,
-α00 + ψ00 * ϕ00 + ψ10 * ϕ01,
-α01 + ψ01 * ϕ00 + ψ11 * ϕ01,
-α10 + ψ00 * ϕ10 + ψ10 * ϕ11,
-α11 + ψ01 * ϕ10 + ψ11 * ϕ11,
-μ0 * ϕ00,
-μ1 * ϕ01,
-μ2 * ϕ10,
-μ3 * ϕ11,
-μ4 * ψ00,
-μ5 * ψ01,
-μ6 * ψ10,
-μ7 * ψ11,
μ2,
μ4,
μ5,
μ6,
μ7
]
solution = sy.solve(system, symbols, dict=True)
sy.pprint(solution)Thanks very much. |
|
Okay let's work this manually a bit.
That leaves us with: The equations at the bottom imply 8 possibilities (i.e. For this example nonlinsolve gives a different error. With solve we get which implies all phi and psi have to be zero but the lambdas can be anything. Indeed it's not hard to see that that is a valid family of solutions. It looks like solve has missed a solution family though. We could also have all the lambdas be zero and then we get a system with the solution which represents a different infinite family of solutions. Obviously I was wrong in my interpretation of the groebner basis... |
|
Thanks for working that out; I appreciate it! That certainly clarifies things for this system. Two concluding questions: i) To complete this pull request, what would you suggest for the assertion in the test, given all of the above? ii) Could you point me to any materials about whether there are classes of systems for which the solvers in sympy are guaranteed to be sound and complete? Similarly, do you have suggestions for reading on the theory behind solvers of this sort in general? Thanks very much |
|
For now I would just make a test that asserts that we don't get the current incorrect output. Then if anything happens to change so that we start to see a different result the test can be updated. It's important to make sure that the comments on the test explain why you are doing that and clearly link to the discussion on github. Ideally create a new issue so that discussion can continue there after this PR is merged. I can't say whether the solvers are guaranteed to be sound and complete in any particular case. It would be impossible to prove that kind of thing anyway given that there is always the possibility of bugs. When it comes to multivariate systems of polynomials I'm not even sure which classes are guaranteed to be solvable anyway. I imagine your example here is fully solvable but I'm sure it's not hard to construct examples that can't possibly be solvable in any normal kind of closed form. It is intended that |
|
This PR just needs an XFAIL test asserting that nonlinsolve does not return the solution that it currently returns i.e. where |
|
Thanks, and apologies for the delay. |
|
@oscarbenjamin Is this good to go? |
|
Yes. Thanks! |
Brief description of what is fixed or changed
The
_solve_using_known_valuesfunction calls theadd_intersection_complementfunction with potential solutions that may be infeasible, such that the feasible region for a possible parameter might be the emptyset.The current code adds a "result" that is an empty dictionary to the result list. This causes key-errors downstream, when code in
_solve_using_known_valuestries to index into that empty dictionary.This pull request changes the code such that if any symbol has an empty feasible region in a solution, that solution is recognized as non-viable, and, it is "removed" from the results list (instead of being added as an empty dictionary).
Other comments
Here is a example of the bug:
Under the old code, this would throw the following:
Under the new code, the output is:
Thanks very much