Opened 6 years ago

Closed 5 years ago

# allow sympy algorithm in solve

Reported by: Owned by: kcrisman major sage-8.2 symbolics rws, mforets, tmonteil, charpent Ralf Stephan Emmanuel Charpentier N/A 8e1e240 8e1e240bb9b8d735251847e67deec1fdb32eccba #23990, #24104, #24062

### Description

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

### comment:1 Changed 5 years ago by mforets

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

### comment:3 follow-up: ↓ 4 Changed 5 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 5 years ago by rws

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

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

### comment:6 Changed 5 years ago by rws

• Dependencies set to #23990

### comment:8 Changed 5 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: ↓ 10 Changed 5 years ago by rws

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

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

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

Version 0, edited 5 years ago by charpent (next)

### comment:11 Changed 5 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 5 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 5 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: ↓ 16 Changed 5 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)
----> 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 5 years ago by rws

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

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

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

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 5 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: ↓ 19 ↓ 20 ↓ 26 Changed 5 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: ↓ 24 Changed 5 years ago by charpent

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: ↓ 22 Changed 5 years ago by rws

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

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

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

### comment:22 in reply to: ↑ 20 ; follow-up: ↓ 23 Changed 5 years ago by charpent

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

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: ↓ 35 Changed 5 years ago by rws

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

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

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

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

### comment:28 Changed 5 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 5 years ago by rws

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

### comment:30 Changed 5 years ago by rws

• Branch set to u/rws/22322

### comment:31 Changed 5 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:

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

And `domain='real'` is recognized too.

### comment:33 Changed 5 years ago by rws

Ah, multivariate doesn't work at the moment.

### comment:34 Changed 5 years ago by git

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

[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: ↓ 37 Changed 5 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: ↓ 38 Changed 5 years ago by rws

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

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

### comment:43 Changed 5 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:

 ​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`

### comment:44 Changed 5 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.

Reviewer name.

### comment:46 Changed 5 years ago by vbraun

• Status changed from positive_review to needs_work

### comment:47 Changed 5 years ago by rws

• Reviewers set to Emmanuel Charpentier

### comment:48 Changed 5 years ago by rws

• Status changed from needs_work to positive_review

### comment:49 Changed 5 years 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}]
**********************************************************************
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 5 years ago by git

• Commit changed from 963cc0e6cb2ea5441422f73beb961b882dfec420 to 8e1e240bb9b8d735251847e67deec1fdb32eccba

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

 ​8e1e240 `22322: dict order-independent doctests`

### comment:51 Changed 5 years ago by rws

• Status changed from needs_work to positive_review

### comment:52 Changed 5 years 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.