solvers: modify symbols ordering in _solve_system

[AaronStiff]

References to other Issues or PRs

Fixes #18208

Brief description of what is fixed or changed

solve was solving a system of equations in a different, albeit valid, way to linsolve. Modified the way solvers._solve_system sorts symbols to preserve order and keep solve and linsolve in sync.

Other comments

Will be adding tests for other solvers to test the same system of equations which caused the issue.

Release Notes

• solvers
  • Modified symbols sorting in solvers._solve_system to ensure solve and linsolve solve the same way

[sympy-bot]

✅

Hi, I am the SymPy bot (v161). 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:

• solvers
  • Modified symbols sorting in solvers._solve_system to ensure solve and linsolve solve the same way (#20928 by @AaronStiff)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.8.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #18208 

#### Brief description of what is fixed or changed
`solve` was solving a system of equations in a different, albeit valid, way to `linsolve`. Modified the way `solvers._solve_system` sorts symbols to preserve order and keep `solve` and `linsolve` in sync.

#### Other comments
Will be adding tests for other solvers to test the same system of equations which caused the issue.

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
  * Modified symbols sorting in solvers._solve_system to ensure solve and linsolve solve the same way
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[oscarbenjamin]

The release note should be for solvers rather than polys.

I think it would be better to put these tests together in one place so that the code for the system and solution is not duplicated and to emphasise the fact that we are asserting that all solvers give the same results.

[AaronStiff]

Oh, whoops. Sure thing.

Just put them all in test_solvers then?

[oscarbenjamin]

Just put them all in test_solvers then?

I would put them wherever the tests for linsolve are.

[AaronStiff]

In test_solveset. Okay

[AaronStiff]

Why is test_elliptic_integrals failing? When I test it on my own machine, it passes. I noticed on some of the past tests that vector.frame.py was also failing, but not this time round.

[oscarbenjamin]

The test is a random test so it can potentially fail for some inputs but not necessarily in a reproducible way.

[AaronStiff]

I see. How are failed random tests like this one usually treated? Will the issue have to be debugged somehow?

[oscarbenjamin]

It's probably not possible to reproduce this. At the top of the test output you can see the random seed as well as the PYTHONHASHSEED which controls set ordering. In principle you can set those seeds when running the tests like:

PYTHONHASHSEED=454480227 bin/test sympy/functions/special/tests/test_elliptic_integrals.py -k test_P --seed=2487161

However the random seed is not set at the start of each separate test so to actually reproduce this you would need to run the exact same tests that bin/test did like:

PYTHONHASHSEED=454480227 bin/test --seed=2487161 --split=1/2

That runs half of the test suite and the failing test won't come up for some time. In between other tests may or may not be skipped due to the presence or absence of different libraries or Python versions etc which might alter the random seed by the time you get to the one that failed.

Without knowing what the seed was I just tried running this random test many times:

from sympy.testing.randtest import (test_derivative_numerically as td,
                                      random_complex_number as randcplx)
from sympy.functions.special.elliptic_integrals import elliptic_pi as P
from sympy import srepr

from sympy.abc import m

for i in range(10000):
    print(i)
    rx, ry = randcplx(), randcplx()
    assert td(P(rx, ry, m), m), str((srepr(rx), srepr(ry)))

I think I let it go about 5000 times and didn't see it fail once. That's why we really need to know the seed when something like this happens.

There are two approaches that would help with this sort of thing in future (ideally we should do both of these):

1. A randomised test like this should print out a complete representation of any random values on failure like:

    x = Float(random())
    assert f(x) < x, srepr(x)
   

2. The test runner should reset the random seed before each individual test so that rerunning an individual test with the same --seed value is more likely to reproduce the random failure.

If you feel like fixing either of those then you can push up a commit to do that. Otherwise you can close and reopen the PR which will rerun the tests and I expect that they will pass.

[AaronStiff]

I think I'll just close and reopen. I'm still working out why nonlinsolve is returning a slightly different solution to linsolve.

[oscarbenjamin]

Looks like the same test failed again. Let's see what happens after another commit is pushed...

[AaronStiff]

Sure. I just added a blank line.

[AaronStiff]

@oscarbenjamin, nonlinsolve is returning slightly different results to linsolve, but not in the same way solve was:

linsolve:
FiniteSet((38 - x3, x3 - 10, 23 - x3, x3, 12 - x7, x7 + 6, 16 - x7, x7, 8, 20, 2, 5, 1, 6, 1, 21, 12, 20, -y11 + y9 + 2, y11 - y9 + 21, -y11 - y7 + y9 + 24, y11 + y7 - y9 - 3, 33 - y7, y7, 27 - y9, y9, 27 - y11, y11))

nonlinsolve:
FiniteSet((38 - x3, x3 - 10, 23 - x3, x3, 12 - x7, x7 + 6, 16 - x7, x7, 8, 20, 2, 5, 1, 6, 1, 21, 12, 20, -y5 + y7 - 1, y5 - y7 + 24, 21 - y5, y5, 33 - y7, y7, 27 - y9, y9, -y5 + y7 - y9 + 24, y5 - y7 + y9 + 3))

It seems only equations with y terms are changed.

The problem is that I don't think nonlinsolve uses the same process as any of the other 3 solvers I've added tests for. It doesn't solve or choose free symbols using RREF. Is it still something to try to harmonize, or should it be left out for now? I'm kind of at a dead end at the moment.

[oscarbenjamin]

The problem is that I don't think nonlinsolve uses the same process as any of the other 3 solvers I've added tests for. It doesn't solve or choose free symbols using RREF.

No, I see that now. I think that probably the substitution function when applied to a linear system is more or less equivalent to Gaussian elimination but it doesn't necessarily do the operations in the same order so it can lead to different symbols being given here.

Probably nonlinsolve should pass any subsystems of linear equations to linsolve but that's a separate issue.

Is it still something to try to harmonize, or should it be left out for now? I'm kind of at a dead end at the moment.

For now I would say add a test along with the others showing the different output that nonlinsolve gives but add code comments clearly explaining that ideally nonlinsolve would give the same outputs as the others. Then if the output from nonlinsolve changes later (e.g. because it is made to use linsolve) the test will fail but the comment will make it clear that it's okay to change the test.

[AaronStiff]

Makes sense!

[AaronStiff]

Actually above I would use brackets like:

    (x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15,
    y0, y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11) = vars

Then again:

    x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15 == symbols('x:16')
    y0, y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11 = symbols('y:12')

In that order? I get a TypeError saying it can't unpack a non-iterable.

The equations would be easier to read if you split a new line after each comma like:

   eqs = [
        x0 + x1 + x2 + x3 - 51,
        x0 + x1 + x4 + x5 - 46,
        x2 + x3 + x6 + x7 - 39,
        ...

Sure thing.

[AaronStiff]

There's no need for partsol to be a separate list

Got it.

[AaronStiff]

linear_eq_to_matrix was already called above. I suggest to do it once like

A, b = linear_eq_to_matrix(eqs+partsol, vars)
assert linsolve(matrix, vars) == linsolve_expected
assert A.gauss_jordan ...

Yes, that simplifies it.

[AaronStiff]

This message should be clearer: the solution set is equivalent and is also correct. However we would prefer to use the same symbols as parameters for the solution to the underdetermined system in all cases if possible. So what we want is a solution that is not just equivalent but also given in the same form.

Sure. I kind of just copied what you wrote and modified it a bit. :)

[AaronStiff]

This is the same as linsolve_expected. Can we derive one from the other?

Which is the same? nonlinsolve_expected is different.

[oscarbenjamin]

This is the same as linsolve_expected. Can we derive one from the other?

Which is the same? nonlinsolve_expected is different.

It looks like I accidentally commented on the commit rather than the diff in the PR. I think I was referring to solve_expected.

[oscarbenjamin]

When squashing commits in a PR it's best just to do it at the end (if needed at all). While working through it gets confusing if the commits commented on disappear and it's also helpful to be able to see any test failures etc - those disappear once you force push.

[AaronStiff]

Oh, okay. Sorry about that.

I think it could probably be derived, yes.

[AaronStiff]

That force push was just some spaces I forgot to add. :)

[oscarbenjamin]

Okay, looks good to me. I can't merge this because it's marked as a "draft" though.

[AaronStiff]

Great. I'll update to a regular PR.

[oscarbenjamin]

The tests failed again due to #20933

[oscarbenjamin]

Looks good