Opened 4 years ago
Closed 3 years ago
#22322 closed enhancement (fixed)
allow sympy algorithm in solve
Reported by:  kcrisman  Owned by:  

Priority:  major  Milestone:  sage8.2 
Component:  symbolics  Keywords:  
Cc:  rws, mforets, tmonteil, charpent  Merged in:  
Authors:  Ralf Stephan  Reviewers:  Emmanuel Charpentier 
Report Upstream:  N/A  Work issues:  
Branch:  8e1e240 (Commits, GitHub, GitLab)  Commit:  8e1e240bb9b8d735251847e67deec1fdb32eccba 
Dependencies:  #23990, #24104, #24062  Stopgaps: 
Description
What it says. See e.g. this post for why.
Change History (52)
comment:1 Changed 4 years ago by
 Cc mforets added
comment:2 Changed 4 years ago by
 Cc tmonteil added
comment:3 followup: ↓ 4 Changed 4 years ago by
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.
comment:4 in reply to: ↑ 3 Changed 4 years ago by
Replying to mforets:
solve
(fromexpression.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 passingex.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).
comment:5 Changed 4 years ago by
Correction, doing ex.lhs()ex.rhs()
is not needed if the SympyConverter
is fixed.
comment:6 Changed 3 years ago by
 Dependencies set to #23990
comment:7 Changed 3 years ago by
 Cc charpent added
comment:8 Changed 3 years ago by
 Status changed from new to needs_info
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) <ipythoninput1425de6a4c8f0a> 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.
comment:9 followup: ↓ 10 Changed 3 years ago by
As I explained in #23990 that is the SymPy convention.
comment:10 in reply to: ↑ 9 Changed 3 years ago by
comment:11 Changed 3 years ago by
 Dependencies changed from #23990 to #23990, #24078
 Milestone changed from sage7.6 to sage8.2
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.
comment:12 Changed 3 years ago by
 Dependencies changed from #23990, #24078 to #23990
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).
comment:13 Changed 3 years ago by
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]]
comment:14 followup: ↓ 16 Changed 3 years ago by
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) <ipythoninput8c8fc37b12ad9> 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()
.
comment:15 Changed 3 years ago by
 Dependencies changed from #23990 to #23990, #24080
comment:16 in reply to: ↑ 14 Changed 3 years ago by
 Dependencies changed from #23990, #24080 to #23990
Replying to rws:
So people would need to know that
operands()
gets the solutions. Or we support indexing as alias tooperands()
.
Supporting expression indexing seems tricky. Also this is an interface change. Better move it to a future ticket.
comment:17 Changed 3 years ago by
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
.
comment:18 followups: ↓ 19 ↓ 20 ↓ 26 Changed 3 years ago by
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.)
comment:19 in reply to: ↑ 18 ; followup: ↓ 24 Changed 3 years ago by
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.
comment:20 in reply to: ↑ 18 ; followup: ↓ 22 Changed 3 years ago by
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?
comment:21 Changed 3 years ago by
 Dependencies changed from #23990 to #23990, #24104
comment:22 in reply to: ↑ 20 ; followup: ↓ 23 Changed 3 years ago by
Replying to rws:
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()
andsolve(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.
comment:23 in reply to: ↑ 22 ; followup: ↓ 35 Changed 3 years ago by
Replying to charpent:
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.
comment:24 in reply to: ↑ 19 Changed 3 years ago by
Replying to charpent:
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.
comment:25 Changed 3 years ago by
Caveat of using RealSet
:
sage: pi.n(30) in RealSet.point(pi) False sage: pi.n(50) in RealSet.point(pi) True
comment:26 in reply to: ↑ 18 Changed 3 years ago by
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.
comment:27 Changed 3 years ago by
 Dependencies changed from #23990, #24104 to #23990, #24104, #24062
comment:28 Changed 3 years ago by
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
.
comment:29 Changed 3 years ago by
Ah no, I'm wrong. We do the right translation.
comment:30 Changed 3 years ago by
 Branch set to u/rws/22322
comment:31 Changed 3 years ago by
 Commit set to 3073840e81e526524eea0ce89cb28bc83723c7d5
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 singleexpression 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

comment:32 Changed 3 years ago by
And domain='real'
is recognized too.
comment:33 Changed 3 years ago by
Ah, multivariate doesn't work at the moment.
comment:34 Changed 3 years ago by
 Commit changed from 3073840e81e526524eea0ce89cb28bc83723c7d5 to eac5f8bf8a12da0373bed819bcd060673a2f7aa4
Branch pushed to git repo; I updated commit sha1. New commits:
eac5f8b  22322: fixes

comment:35 in reply to: ↑ 23 Changed 3 years ago by
Replying to rws:
[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...
[ Resnip ... ]
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 nonprofessionalmathematician 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.
comment:36 followup: ↓ 37 Changed 3 years ago by
Stefan,
I'm a bit lost about your various tickets. I suggest that you make a metaticket (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 nolongeruseful parts should be also possible.
comment:37 in reply to: ↑ 36 ; followup: ↓ 38 Changed 3 years ago by
Replying to charpent:
Stefan,
I'm a bit lost about your various tickets. I suggest that you make a metaticket (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 nolongeruseful 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.
comment:38 in reply to: ↑ 37 Changed 3 years ago by
 Status changed from needs_info to positive_review
Replying to rws: [ 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
.
comment:39 Changed 3 years ago by
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...).
comment:40 Changed 3 years ago by
 Status changed from positive_review to needs_work
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.
comment:41 Changed 3 years ago by
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.
comment:42 Changed 3 years ago by
 Branch changed from u/rws/22322 to u/rws/223221
comment:43 Changed 3 years ago by
 Commit changed from eac5f8bf8a12da0373bed819bcd060673a2f7aa4 to 963cc0e6cb2ea5441422f73beb961b882dfec420
 Status changed from needs_work to needs_review
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 223221

86e9c20  24104: move doctests

91cb91f  24104: solve no longer calls itself

0f6e826  Move the single expression solver into standalone function for granularity plus some mild cleanup.

4048e25  Merge branch 'public/symbolcs/merge_solves24104' of git://trac.sagemath.org/sage into 223221

963cc0e  22322: allow sympy algorithm in solve

comment:44 Changed 3 years ago by
 Status changed from needs_review to positive_review
comment:45 Changed 3 years ago by
Reviewer name.
comment:46 Changed 3 years ago by
 Status changed from positive_review to needs_work
comment:47 Changed 3 years ago by
 Reviewers set to Emmanuel Charpentier
@charpent I hope you don't mind me adding your name.
comment:48 Changed 3 years ago by
 Status changed from needs_work to positive_review
comment:49 Changed 3 years ago by
 Status changed from positive_review to needs_work
On OSX:
********************************************************************** File "src/sage/symbolic/relation.py", line 894, in sage.symbolic.relation.solve Failed example: solve([x^2  y^2/exp(x), y1], x, y, algorithm='sympy') Expected: [{y: 1, x: 2*lambert_w(1/2)}] Got: [{x: 2*lambert_w(1/2), y: 1}] ********************************************************************** File "src/sage/symbolic/relation.py", line 902, in sage.symbolic.relation.solve Failed example: solve([x + y + z + t, z  t], x, y, z, t, algorithm='sympy') Expected: [{z: t, x: y}] Got: [{x: y, z: t}] ********************************************************************** File "src/sage/symbolic/relation.py", line 904, in sage.symbolic.relation.solve Failed example: solve([x^2+y+z, y+x^2+z, x+y+z^2], x, y,z, algorithm='sympy') Expected: [{y: (z + 1)*z, x: z}, {y: z^2 + z  1, x: z + 1}] Got: [{x: z, y: (z + 1)*z}, {x: z + 1, y: z^2 + z  1}] ********************************************************************** 1 item had failures: 3 of 112 in sage.symbolic.relation.solve [377 tests, 3 failures, 185.81 s]  sage t long src/sage/symbolic/relation.py # 3 doctests failed 
comment:50 Changed 3 years ago by
 Commit changed from 963cc0e6cb2ea5441422f73beb961b882dfec420 to 8e1e240bb9b8d735251847e67deec1fdb32eccba
Branch pushed to git repo; I updated commit sha1. New commits:
8e1e240  22322: dict orderindependent doctests

comment:51 Changed 3 years ago by
 Status changed from needs_work to positive_review
comment:52 Changed 3 years ago by
 Branch changed from u/rws/223221 to 8e1e240bb9b8d735251847e67deec1fdb32eccba
 Resolution set to fixed
 Status changed from positive_review to closed
AFAIK, also useful for equations where the unkown appears as an exponent, cf. solve(4.94 * 1.062^x == 15, x)