allow sympy algorithm in solve

[kcrisman]

What it says. See e.g. this post for why.

Depends on #23990 Depends on #24104 Depends on #24062

CC: @rwst @mforets @sagetrac-tmonteil @EmmanuelCharpentier

Component: symbolics

Author: Ralf Stephan

Branch/Commit: 8e1e240

Reviewer: Emmanuel Charpentier

Issue created by migration from https://trac.sagemath.org/ticket/22322

[c22b6800-ec0b-4cbf-94c4-0a74aecc2093]

comment:1

AFAIK, also useful for equations where the unkown appears as an exponent, cf. solve(4.94 * 1.062^x == 15, x)

[c22b6800-ec0b-4cbf-94c4-0a74aecc2093]

comment:3

solve (from expression.pyx) passes a relational expression to Maxima's solve, but this doesn't seem to be available in sympy:

sage: (x==1)._maxima_()
_SAGE_VAR_x=1
sage: (x==1)._sympy_()
Traceback (most recent call last)
...
NotImplementedError: relation

do you have some idea on how to tackle this? otherwise, what about passing ex.lhs()-ex.rhs() to sympy's solve.

[rwst]

comment:4

Replying to @mforets:

solve (from expression.pyx) passes a relational expression to Maxima's solve, but this doesn't seem to be available in sympy: ... do you have some idea on how to tackle this? otherwise, what about passing ex.lhs()-ex.rhs() to sympy's solve.

Yes, that would be it for the equalities. For inequalities SympyConverter.relation() is needed (atm the superclass's relation is called).

[rwst]

comment:5

Correction, doing ex.lhs()-ex.rhs() is not needed if the SympyConverter is fixed.

[rwst]

Dependencies: #23990

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:8

Hmmm...

sage: var("a,b")
(a, b)
sage: sympy.sympify(a==b)._sage_()
a == b
sage: sympy.sympify(a==b).sage()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-14-25de6a4c8f0a> in <module>()
----> 1 sympy.sympify(a==b).sage()

AttributeError: 'Equality' object has no attribute 'sage'

Is that the expected behaviour ? Or do you plan to enhance the .sage() method for various objects ?

It is not obvious to me why (e. g.) mathematica objets can (sometimes...) be converted back to sage via the .sage() (exposed) method, whereas Sympy objects have to use the ._sage_() (unexposed method.

[rwst]

comment:9

As I explained in #23990 that is the SymPy convention.

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:10

Replying to @rwst:

As I explained in #23990 that is the SymPy convention.

Wups, bad wrong ticket. I thought that my comment was cancelled ; it was not.

Apologies...

[rwst]

Changed dependencies from #23990 to #23990, #24078

[rwst]

comment:11

Now that (with #23990) relations are translated from/to SymPy there are still things needed. For example we need to translate any assumptions on variables (and maybe anon. functions?) that were made using assume and var('x', domain=...). This is not only relevant for this ticket, any SymPy operation may access the SymPy knowledge base any time. So I think this should be handled the same way as with Maxima, i.e. at the time the assume/var calls are made. This is #24078.

[rwst]

comment:12

Note that SymPy's solve does not respect the global assumptions database:

In [12]: global_assumptions
Out[12]: AssumptionsContext([Q.is_true(x > 2), Q.positive(x)])

In [13]: solve(x<3)
Out[13]: -∞ < x ∧ x < 3

In [15]: solve((x>2,x<3))
Out[15]: 2 < x ∧ x < 3

So the dependency is irrelevant (at the moment).

[rwst]

Changed dependencies from #23990, #24078 to #23990

[rwst]

comment:13

Our solve doesn't either, so it all is localized in the solve arguments:

sage: assume(x>2)
sage: solve(x<3,x)
[[x < 3]]
sage: solve((x>2,x<3),x)
[[2 < x, x < 3]]

[rwst]

comment:14

Our solve's output is quite opaque IMHO:

sage: solve((x^2<1),x)
[[x > -1, x < 1]]
sage: solve((x^2>1),x)
[[x < -1], [x > 1]]

From a glance, which is AND, which is OR? So I would like to use And(...) and Or(...) as in SymPy but:

sage: Or = function('Or')
sage: ex = Or(x<0, x>1)
sage: ex
Or(x < 0, x > 1)
sage: ex.args()
(x,)
sage: ex.operands()
[x < 0, x > 1]

sage: ex[0]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-8-c8fc37b12ad9> in <module>()
----> 1 ex[Integer(0)]

TypeError: 'sage.symbolic.expression.Expression' object does not support indexing

So people would need to know that operands() gets the solutions. Or we support indexing as alias to operands().

[rwst]

Changed dependencies from #23990 to #23990, #24080

[rwst]

comment:16

Replying to @rwst:

So people would need to know that operands() gets the solutions. Or we support indexing as alias to operands().

Supporting expression indexing seems tricky. Also this is an interface change. Better move it to a future ticket.

[rwst]

Changed dependencies from #23990, #24080 to #23990

[rwst]

comment:17

I cannot believe devs maintained two different versions of solve in symbolic/expression.pyx and symbolic/relation.py. expression.pyx should only be an interface to a function elsewhere, preferably in solve.py.

[c22b6800-ec0b-4cbf-94c4-0a74aecc2093]

comment:18

Note that SymPy's solve does not respect the global assumptions database:

shall we consider interfacing with SymPy's solveset? it is said to supersede solve sooner or later.

i find that the interface (equation, variable, domain) is very clear. isn't it preferable than relying on global assumptions? (or to fallback into the global assumptions if the domain is not provided.)

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:19

Replying to @mforets:

Note that SymPy's solve does not respect the global assumptions database:

shall we consider interfacing with SymPy's solveset? it is said to supersede solve sooner or later.

Fly in the ointment : solve can solve systems of equations, whereas solveset (which is, indeed, cleaner on more than one aspect), solves only equations.

In order to solve system, we would have to map solveset to all the equations in the system (each for one variable (which one ?)), and return the intersection of the results, which is currently problematic.

i find that the interface (equation, variable, domain) is very clear. isn't it preferable than relying on global assumptions? (or to fallback into the global assumptions if the domain is not provided.)

The domain is only a part of the possible assumptions : e. g., maxima often asks if some expression is positive, if such variable is integer, etc... We need to maintain a "knowledge base" about outr variables as least as rich as our current (= Maxima's) facts base.

A few remarks :

1) I note that Sympy's knowledge base is richer than ours, BTW.

2) We might use profitably something analogous to Mathematica's Assuming statement, which does something (roughly )equivalent to :

def assuming(condset, statement):
    try:
        OldAssumptions=assumptions()
        Assumptions=OldAssumptions.copy()
        Assumptions.extend(condset)
        forget(assumptions())
        assume(Assumptions)
        result=<do statement, somehow>
    finally:
        forget(assumptions())
        assume(OldAssumptions)

The wasp in the ointment is of course that statement is evaluated before being passed to Assuming : Python has no macroes...

I think that assuming could be implemented somehow as a decorator, but I'm still to understand that.

[rwst]

comment:20

Replying to @mforets:

i find that the interface (equation, variable, domain) is very clear. isn't it preferable than relying on global assumptions?

I have abandoned the idea of respecting assumptions after I saw that SymPy doesn't either. We don't have to reproduce Maxima results so we're good.

solve can solve systems of equations, whereas solveset (which is, indeed, cleaner on more than one aspect), solves only equations.

If they replace it will be unified. And I think Sage's ex.solve() and solve(ex) should be unified too, if for no other reason than for better maintenance.

In order to solve system, we would have to map solveset to all the equations in the system (each for one variable (which one ?)), and return the intersection of the results, which is currently problematic.

It's always problematic because you need to decide ordering of expressions without variables, and that can take time because of the possible need to increase precision. The question is, should we support SymPy's solve or immediately use solveset exclusively?

We might use profitably something analogous to Mathematica's Assuming

Please review #10035.

We need to maintain a "knowledge base" about outr variables as least as rich as our current (= Maxima's) facts base.

Isn't that the same as a local context?

[rwst]

Changed dependencies from #23990 to #23990, #24104

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:22

Replying to @rwst:

Replying to @mforets:

And me, it seems...

i find that the interface (equation, variable, domain) is very clear. isn't it preferable than relying on global assumptions?

I have abandoned the idea of respecting assumptions after I saw that SymPy doesn't either. We don't have to reproduce Maxima results so we're good.

That's pushing the can down the road : in order to get the solutions, we're bound to solve a system comprising the solutions obtained without the assumptions and the assumptions themselves, which might be unsolvable (e. g. cubic equation with a solution in casus irreducibilis + constraints showing that the roots are real)...

solve can solve systems of equations, whereas solveset (which is, indeed, cleaner on more than one aspect), solves only equations.

If they replace it will be unified. And I think Sage's ex.solve() and solve(ex) should be unified too, if for no other reason than for better maintenance.

Indeed. If...

In order to solve system, we would have to map solveset to all the equations in the system (each for one variable (which one ?)), and return the intersection of the results, which is currently problematic.

It's always problematic because you need to decide ordering of expressions without variables, and that can take time because of the possible need to increase precision. The question is, should we support SymPy's solve or immediately use solveset exclusively?

IMHO, both. But that might be a lot of work.. Depends on your prognosis on the time necessary for the Sympy team to coax solveset in solving systems...

We might use profitably something analogous to Mathematica's Assuming

Please review #10035.

Done (positively). But that is a possible solution only for function calls supporting hold. solve isn't one of them.

We need to maintain a "knowledge base" about outr variables as least as rich as our current (= Maxima's) facts base.

Isn't that the same as a local context?

I do not understand that (yet). Care to amplify ?

And, BTW, thank you for the huge work you are doing on Sympy, which might turn out very important for Sage.

[rwst]

comment:23

Replying to @EmmanuelCharpentier:

Replying to @mforets: Please review #10035.

Done (positively). But that is a possible solution only for function calls supporting hold. solve isn't one of them.

It motivates me to create a context infrastructure where you can say:

sage: mycontext = context(x>0, hold)
sage: with mycontext:
sage:     solve(...)

or

sage: mycontext.start()
sage: ...
sage: solve(...)
sage: mycontext.stop()

We need to maintain a "knowledge base" about outr variables as least as rich as our current (= Maxima's) facts base.

Isn't that the same as a local context?

I do not understand that (yet). Care to amplify ?

I gave the first answer in order to give an idea.

And, BTW, thank you for the huge work you are doing on Sympy, which might turn out very important for Sage.

I remember reading in devel that previous developers hoped so too, also because at that time development on Pynac stopped and practically no one worked on symbolics in Sage. While that's no longer the case, solvers and integration need much help from SymPy, Fricas, and Giac.

[rwst]

comment:24

Replying to @EmmanuelCharpentier:

In order to solve system, we would have to map solveset to all the equations in the system (each for one variable (which one ?)), and return the intersection of the results, which is currently problematic.

Maybe we should not worry about doing intersections, convert to RealSet, intersect, convert back. RealSet uses RealInterval for comparison and I guess the problems have been thought about by its developers. But what to convert the RealSet to in the end? For starters the output of solve at the moment, i.e., lists of lists. The conversion RealSet-->lists might be of broader interest.

[rwst]

comment:25

Caveat of using RealSet:

sage: pi.n(30) in RealSet.point(pi)
False
sage: pi.n(50) in RealSet.point(pi)
True

[rwst]

comment:26

Replying to @mforets:

i find that the interface (equation, variable, domain) is very clear. isn't it preferable than relying on global assumptions? (or to fallback into the global assumptions if the domain is not provided.)

It is. Propagating the solveset interface to our solve would mean at least supporting the option domain='real' with complex the default.

[rwst]

Changed dependencies from #23990, #24104 to #23990, #24104, #24062

[rwst]

comment:28

SymPy seems to make a clear distinction between the operators in x>0 and Gt(x,0). I think it's wrong that (x>0)._sympy_() returns x>0.

[rwst]

comment:29

Ah no, I'm wrong. We do the right translation.

[rwst]

Branch: u/rws/22322

[rwst]

comment:31

This branch works quite well and can be used and reviewed. It lacks doctests and documentation, which I'll add next. But please experiment with it. We use solveset with single-expression univariate tasks, solve for the rest. The output is unified to Sage's list "format".

---

Last 10 new commits:

7899ea9	20204: convert SymPy abstract functions
b261ec3	23923: fix doctest
5010acf	Merge branch 'u/rws/23923' of git://trac.sagemath.org/sage into t/20204/20204
14b1c52	20204: more doctests
9c4119a	16801: Conversion of psi(x,y) to/from SymPy
d5c60e3	24062: pkg version/chksum
1cbc23e	24062: remove obsolote patches; adapt remaining
94f359f	24062: doctest fixes
cf5c81a	Merge branch 'u/rws/24062' of git://trac.sagemath.org/sage into t22322
3073840	22322: sympy algorithm in solve

[rwst]

Commit: 3073840

[rwst]

Author: Ralf Stephan

[rwst]

comment:32

And domain='real' is recognized too.

[rwst]

comment:33

Ah, multivariate doesn't work at the moment.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 3073840 to eac5f8b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

eac5f8b

22322: fixes

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:35

Replying to @rwst:

[Snip ... ]

It motivates me to create a context infrastructure where you can say:

sage: mycontext = context(x>0, hold)
sage: with mycontext:
sage:     solve(...)

or

sage: mycontext.start()
sage: ...
sage: solve(...)
sage: mycontext.stop()

Maybe #24119 might be relevant ? It doesn't medles with evaluation (or not) of the body of the with, it just brackets this body with temporary updates to the facts database. Not the problem you are trying to solve, but might be handy...

[ Re-snip ... ]

I remember reading in devel that previous developers hoped so too, also because at that time development on Pynac stopped and practically no one worked on symbolics in Sage. While that's no longer the case, solvers and integration need much help from SymPy, Fricas, and Giac.

In my (limited) experience, the ability od Sage to call "at will" SymPy, fricas, giac (and even maple and mathematica, but shhh...) is extremely useful for practical work. Add Sage\TeX to the mixture (the notebook is nice for prottyping, but I'm sold on "literate programming" for any serious publication work...), and you get a tool whose usefulness surpasses those of any of its components.

IMHO, Sage needs a lot of development in the "engineering/scientific/undergrad maths" department, which is not, to say the least, the core interest of its core developers. I see the work in this direction as an investment : if we make a tool useful to those non-professional-mathematician users,I would be quite surprised that none of them would dive in Sage development in these departments, thus relieving the algebraists, geometers and number theorist who are the core developers of Sage of the infrastructure and maintenance work.

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:36

Stefan,

I'm a bit lost about your various tickets. I suggest that you make a meta-ticket (say, "Work in progress on SymPy") in which you merge the various pieces you want tested (with their dependencies) : that would ease the retrieval. Reverting no-longer-useful parts should be also possible.

[rwst]

comment:37

Replying to @EmmanuelCharpentier:

Stefan,

I'm a bit lost about your various tickets. I suggest that you make a meta-ticket (say, "Work in progress on SymPy") in which you merge the various pieces you want tested (with their dependencies) : that would ease the retrieval.

Why? Just check out the branch. It contains everything.

Reverting no-longer-useful parts should be also possible.

Can you explain what you are trying to do? We don't need to use the ticket for this, that's also what zulip.sagemath.org is for, just message me please.

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:38

Replying to @rwst: [ Snip ... ]

Why? Just check out the branch. It contains everything.

Aha ! I didn't understand that you had merged the necessary tickets.

Anyway : on top of 8.1.beta9 + #23980+2#4026+#24119, passes ptestlong with no errors whatsoever.

It needs doc, however. And some options do not seem to work on multivariate systems :

sage: solve(x^3+1==0,x,domain="real") ## domain="real" does not work in Maxima's solver...
[x == 1/2*I*sqrt(3)*(-1)^(1/3) - 1/2*(-1)^(1/3), x == -1/2*I*sqrt(3)*(-1)^(1/3) - 1/2*(-1)^(1/3), x == (-1)^(1/3)]
sage: solve(x^3+1==0,x,algorithm="sympy",domain="real") # ...but works with Sympys's solver ...
[x == -1]
sage: solve([x^2+y^2 == 1, y^2 == x^3 + x + 1], [x, y],algorithm="sympy",domain=
....: "real") ### ... but not for systems.

[{y: -1, x: 0},
 {y: 1, x: 0},
 {y: -sqrt(-1/2*I*sqrt(3) + 3/2), x: -1/2*I*sqrt(3) - 1/2},
 {y: sqrt(-1/2*I*sqrt(3) + 3/2), x: -1/2*I*sqrt(3) - 1/2},
 {y: -sqrt(1/2*I*sqrt(3) + 3/2), x: 1/2*I*sqrt(3) - 1/2},
 {y: sqrt(1/2*I*sqrt(3) + 3/2), x: 1/2*I*sqrt(3) - 1/2}]

Notwithstandig that, positive_review.

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:39

Just a thought : I see that the patchbots fail to build the package, for lack of source tarball. It could be useful to add a pointer to it in the ticket description (ISTR that thos worked for #24026...).

[rwst]

comment:40

Thanks. I have to do more Work, at least the merge of #24104 changes and necessary documentation. The tarball will be no longer necessary with the next release because the SymPy upgrade is in there.

[rwst]

comment:41

We'll need to work on #22024 next because solve(x^5 + 3*x^3 + 7, x, algorithm='sympy') will fail. But not for this ticket.

[rwst]

Changed branch from u/rws/22322 to u/rws/22322-1

[rwst]

Changed commit from eac5f8b to 963cc0e

[rwst]

comment:43

This is a new branch because of differences in the dependencies. I also added the missing translation of LambertW. Another problem is that SymPy's solveset is not always better than its solve esp. with nonlinear equations, compare e.g. x+2**x==0. This needs more work later.

Please review, these should be the last big changes.

---

New commits:

269c2de	Merge branch 'develop' into 22322-1
86e9c20	24104: move doctests
91cb91f	24104: solve no longer calls itself
0f6e826	Move the single expression solver into stand-alone function for granularity plus some mild cleanup.
4048e25	Merge branch 'public/symbolcs/merge_solves-24104' of git://trac.sagemath.org/sage into 22322-1
963cc0e	22322: allow sympy algorithm in solve

[7822f248-ba56-45f1-ab3d-4de7482bdf9f]

comment:44

On top of 8.1.beta9 + #23980+2#4026+#24119, passes ptestlong with no errors whatsoever.

[tscrim]

comment:45

Reviewer name.