Opened 3 years ago

Closed 22 months ago

#22322 closed enhancement (fixed)

allow sympy algorithm in solve

Reported by: kcrisman Owned by:
Priority: major Milestone: sage-8.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) Commit: 8e1e240bb9b8d735251847e67deec1fdb32eccba
Dependencies: #23990, #24104, #24062 Stopgaps:

Description

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

Change History (52)

comment:1 Changed 2 years ago by mforets

  • Cc mforets added

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

comment:2 Changed 2 years ago by tmonteil

  • Cc tmonteil added

comment:3 follow-up: Changed 2 years ago by mforets

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 2 years ago by rws

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).

comment:5 Changed 2 years ago by rws

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

comment:6 Changed 2 years ago by rws

  • Dependencies set to #23990

comment:7 Changed 2 years ago by charpent

  • Cc charpent added

comment:8 Changed 2 years ago by charpent

  • 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)
<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.

comment:9 follow-up: Changed 2 years ago by rws

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

comment:10 in reply to: ↑ 9 Changed 2 years ago by charpent

Replying to rws:

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...

Last edited 2 years ago by charpent (previous) (diff)

comment:11 Changed 2 years ago by rws

  • Dependencies changed from #23990 to #23990, #24078
  • Milestone changed from sage-7.6 to sage-8.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 2 years ago by rws

  • 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 2 years ago by rws

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 follow-up: Changed 2 years ago by rws

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().

comment:15 Changed 2 years ago by rws

  • Dependencies changed from #23990 to #23990, #24080

comment:16 in reply to: ↑ 14 Changed 2 years ago by rws

  • 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 to operands().

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

comment:17 Changed 2 years ago by rws

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 follow-ups: Changed 2 years ago by 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.

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 ; follow-up: Changed 2 years ago by charpent

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 ; follow-up: Changed 2 years ago by rws

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 2 years ago by rws

  • Dependencies changed from #23990 to #23990, #24104

comment:22 in reply to: ↑ 20 ; follow-up: Changed 2 years ago by charpent

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() 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.

comment:23 in reply to: ↑ 22 ; follow-up: Changed 2 years ago by rws

Replying to charpent:

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.

comment:24 in reply to: ↑ 19 Changed 2 years ago by rws

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 2 years ago by rws

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 2 years ago by rws

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 2 years ago by rws

  • Dependencies changed from #23990, #24104 to #23990, #24104, #24062

comment:28 Changed 2 years ago by rws

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 2 years ago by rws

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

comment:30 Changed 2 years ago by rws

  • Branch set to u/rws/22322

comment:31 Changed 2 years ago by rws

  • Authors set to Ralf Stephan
  • 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 single-expression univariate tasks, solve for the rest. The output is unified to Sage's list "format".


Last 10 new commits:

7899ea920204: convert SymPy abstract functions
b261ec323923: fix doctest
5010acfMerge branch 'u/rws/23923' of git://trac.sagemath.org/sage into t/20204/20204
14b1c5220204: more doctests
9c4119a16801: Conversion of psi(x,y) to/from SymPy
d5c60e324062: pkg version/chksum
1cbc23e24062: remove obsolote patches; adapt remaining
94f359f24062: doctest fixes
cf5c81aMerge branch 'u/rws/24062' of git://trac.sagemath.org/sage into t22322
307384022322: sympy algorithm in solve

comment:32 Changed 2 years ago by rws

And domain='real' is recognized too.

comment:33 Changed 2 years ago by rws

Ah, multivariate doesn't work at the moment.

comment:34 Changed 2 years ago by git

  • Commit changed from 3073840e81e526524eea0ce89cb28bc83723c7d5 to eac5f8bf8a12da0373bed819bcd060673a2f7aa4

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

eac5f8b22322: fixes

comment:35 in reply to: ↑ 23 Changed 2 years ago by charpent

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...

[ 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.

comment:36 follow-up: Changed 2 years ago by charpent

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.

comment:37 in reply to: ↑ 36 ; follow-up: Changed 2 years ago by rws

Replying to charpent:

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.

comment:38 in reply to: ↑ 37 Changed 2 years ago by charpent

  • 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 2 years ago by charpent

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 2 years ago by rws

  • 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 2 years ago by rws

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 2 years ago by rws

  • Branch changed from u/rws/22322 to u/rws/22322-1

comment:43 Changed 2 years ago by rws

  • 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:

269c2deMerge branch 'develop' into 22322-1
86e9c2024104: move doctests
91cb91f24104: solve no longer calls itself
0f6e826Move the single expression solver into stand-alone function for granularity plus some mild cleanup.
4048e25Merge branch 'public/symbolcs/merge_solves-24104' of git://trac.sagemath.org/sage into 22322-1
963cc0e22322: allow sympy algorithm in solve

comment:44 Changed 2 years ago by charpent

  • Status changed from needs_review to positive_review

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

comment:45 Changed 2 years ago by tscrim

Reviewer name.

comment:46 Changed 2 years ago by vbraun

  • Status changed from positive_review to needs_work

comment:47 Changed 2 years ago by rws

  • Reviewers set to Emmanuel Charpentier

@charpent I hope you don't mind me adding your name.

comment:48 Changed 2 years ago by rws

  • Status changed from needs_work to positive_review

comment:49 Changed 23 months ago by vbraun

  • 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), y-1], 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 22 months ago by git

  • Commit changed from 963cc0e6cb2ea5441422f73beb961b882dfec420 to 8e1e240bb9b8d735251847e67deec1fdb32eccba

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

8e1e24022322: dict order-independent doctests

comment:51 Changed 22 months ago by rws

  • Status changed from needs_work to positive_review

comment:52 Changed 22 months ago by vbraun

  • Branch changed from u/rws/22322-1 to 8e1e240bb9b8d735251847e67deec1fdb32eccba
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.