Opened 5 months ago
Closed 4 months ago
#27958 closed enhancement (fixed)
enhance the integration by using also giac and sympy
Reported by:  chapoton  Owned by:  

Priority:  major  Milestone:  sage8.9 
Component:  symbolics  Keywords:  integration 
Cc:  ghmwageringel, jakobkroeker, kcrisman, pbruin, rws, slabbe, slelievre, tmonteil, vdelecroix, vklein, zimmerma, mforets, mmezzarobba  Merged in:  
Authors:  Frédéric Chapoton  Reviewers:  Thierry Monteil 
Report Upstream:  N/A  Work issues:  
Branch:  d4371e3 (Commits)  Commit:  d4371e3be96a2386f65b57d2e9f36f69eaf71587 
Dependencies:  Stopgaps: 
Description (last modified by )
Change History (52)
comment:1 Changed 5 months ago by
 Branch set to u/chapoton/27958
 Cc rws slelievre added
 Commit set to 3ddd815e4fbbef38f3bb2b9b1c75674d20c44cd0
 Status changed from new to needs_review
comment:2 Changed 5 months ago by
 Commit changed from 3ddd815e4fbbef38f3bb2b9b1c75674d20c44cd0 to e297e21f5bae8d5b61e7701820bb691857bea647
Branch pushed to git repo; I updated commit sha1. New commits:
e297e21  trac 27958 polishing the details

comment:3 Changed 5 months ago by
 Cc slabbe kcrisman ghmwageringel pbruin jakobkroeker added
comment:4 Changed 5 months ago by
 Cc tmonteil added
 Keywords integration added
comment:5 Changed 5 months ago by
sometimes sympy takes ages, so that's why it comes last:
sage: integrate(sqrt(1m*sin(x)^2),x,algorithm='sympy') ... hangs ...
comment:6 Changed 5 months ago by
 Status changed from needs_review to needs_work
some failing doctests, work in progress
comment:7 Changed 5 months ago by
 Commit changed from e297e21f5bae8d5b61e7701820bb691857bea647 to 4a9fbbb2656d3af8fe015ad00fde43daf71cf0fc
Branch pushed to git repo; I updated commit sha1. New commits:
4a9fbbb  trac 27958 fixing doctests

comment:8 Changed 5 months ago by
 Cc mforets added
This seems to work well enough. Waiting for the patchbots' green lights.
comment:9 Changed 5 months ago by
Having more integration algorithms tried by default is a very good move. However, doctests should get more analyzed:
This part of the patch looks weird:
can do, but Sage currently can't do:: +can do, but Sage currently cannot do::  sage: integrate(log(x)*exp(x^2), x) # todo  Mathematica can do this  integrate(e^(x^2)*log(x), x) + sage: integrate(log(x)*exp(x^2), x) + 1/2*sqrt(pi)*erf(x)*log(x)  x*hypergeometric((1/2, 1/2), (3/2, 3/2), x^2)
Even if Sage is not able to perform the nice simplifications, this looks pretty correct (though i am not analyst):
sage: f = 1/2*sqrt(pi)*erf(x)*log(x)  x*hypergeometric((1/2, 1/2), (3/2, 3/2), x^2) sage: g = f.derivative(x) sage: h = g  log(x)*exp(x^2) sage: h.full_simplify() 1/18*(4*x^3*hypergeometric((3/2, 3/2), (5/2, 5/2), x^2)  18*x*hypergeometric((1/2, 1/2), (3/2, 3/2), x^2) + 9*sqrt(pi)*erf(x))/x sage: h.plot(x,0,1)
That function looks like the zero function ?
comment:10 Changed 5 months ago by
 Reviewers set to Thierry Monteil
Actually, mathematica gives a very similar result : https://www.wolframalpha.com/input/?i=integrate+log(x)*exp(x%5E2)
comment:11 Changed 5 months ago by
 Commit changed from 4a9fbbb2656d3af8fe015ad00fde43daf71cf0fc to 8eaf4c1c9df9d4097dae5483a911fe71f95f5ab8
Branch pushed to git repo; I updated commit sha1. New commits:
8eaf4c1  trac 27958 some details in doc

comment:12 Changed 5 months ago by
thx, I have modified this part of the doc
comment:13 Changed 5 months ago by
Also, in
+Todo, Mathematica can do this and gets `\pi^2/15`:: sage: integrate(log(1+sqrt(1+4*x)/2)/x, x, 0, 1) Traceback (most recent call last): ... ValueError: Integral is divergent.
Well, the integral is divergent, as the function behaves like log(3/2)/x
around zero. And i bet mathematica never returned pi^2/15
, maybe the original example was modified along the time without noticing.
comment:14 Changed 5 months ago by
Some other remarks:
 when you write that Sage can "now" compute some integral, it might be good to refer to
:trac:`27958`
so that interested people could see how and where the issue was fixed.  when you add some conversions dict that are not indirectly tested with the examples, it might be good to doctest them explicitely (i am thinking in particular to the
sign
conversion for fricas, but maybe there are other that you added without the need to let the new example work, but only you can track those without effort).  why not also adding
fricas
to the tried integrators ? Though it is optional, it could be nice to have it implicitely used when the user has it installed. This may be part of another ticket though.
comment:15 Changed 5 months ago by
 Cc zimmerma added
Let me CC @zimmerma since some doctests in the Sage book are modified, so that perhaps the example of exp(−x^2 )*ln(x)
might not be good anymore to illustrate numerical integration as a fallback.
comment:16 Changed 5 months ago by
 Commit changed from 8eaf4c1c9df9d4097dae5483a911fe71f95f5ab8 to 406e27278cc325b9fe0b84b448da6cbad537bba8
Branch pushed to git repo; I updated commit sha1. New commits:
406e272  trac 27958 doc polishing

comment:17 Changed 5 months ago by
 Status changed from needs_work to needs_review
Thx, I have made the suggested changes.
Let us keep the possibility of also adding the fricas integrator for another ticket.
comment:18 Changed 5 months ago by
full green bot, please review
comment:19 Changed 5 months ago by
Any timing information? I'd be really concerned about stuff slowing down by trying these additional integrators if there was significant overhead somewhere we didn't know about.
comment:20 Changed 5 months ago by
+1 for this PR.
I'd be really concerned about stuff slowing down by trying these additional integrators
Since it still uses maxima
as its first try, the time will get worse only in those cases were Sage was returning an error or an unevaluated integral (both results are useless, i would say, so one shouldn't care).
comment:21 Changed 5 months ago by
as I said in comment:5, sometimes sympy seems to hang forever. Here is one case :
sage: %time integrate(sqrt(13*sin(x)^2),x) CPU times: user 252 ms, sys: 3.93 ms, total: 256 ms Wall time: 274 ms integrate(sqrt(3*sin(x)^2 + 1), x) sage: %time integrate(sqrt(13*sin(x)^2),x,algorithm='giac') CPU times: user 2.85 ms, sys: 3.85 ms, total: 6.7 ms Wall time: 85.3 ms integrate(sqrt(3*sin(x)^2 + 1), x) sage: %time integrate(sqrt(13*sin(x)^2),x,algorithm='sympy') NOTHING HAPPENS...
So this behaviour is now inherited by our integration.
One can hope that the users will CTRLC. Or use a timeout ??
It seems to me that this is a reasonable prize to pay for having more success, in particular with integrals contaiing abs(), now that we have disabled maxima's abs_integrate
.
EDIT:
in this case, sympy succeeds after not so long!!
sage: %time integrate(sqrt(13*sin(x)^2),x) CPU times: user 1min 47s, sys: 346 ms, total: 1min 47s Wall time: 1min 50s elliptic_e(x, 3)
comment:22 Changed 5 months ago by
I'd be careful, because not every user perhaps knows about CtrlC, especially in a notebook/worksheet situation. But of course if we get a significantly higher proportion of integrals it is good. I just wanted to raise this.
comment:23 Changed 5 months ago by
 Commit changed from 406e27278cc325b9fe0b84b448da6cbad537bba8 to 43c3bfa2c30ca2f406ddaa149f412116ec6fca62
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
43c3bfa  enhance our integration capacities by using maxima, giac and sympy

comment:24 Changed 5 months ago by
squashed branch, with more doctests showing the new capacities
comment:25 followup: ↓ 26 Changed 5 months ago by
and the bot is still fully green.
Note that sometimes it is giac, rather than sympy, that hangs, for example in
integrate(sqrt(x + sqrt(x)), x)
EDIT
In giac itself, integrate(sqrt(x + sqrt(x)),x);
returns an error quickly.
EDIT Does not hang in Python2, only in Python3 !!!
EDIT: Apparently, in sage+py3, the giac interface cannot recover from a crash of Giac..
giac('int(sqrt(x+sqrt(x)),x)') Giac crashed  automatically restarting. HANGS..
comment:26 in reply to: ↑ 25 Changed 5 months ago by
Replying to chapoton:
[sometimes giac hangs, for example in]
integrate(sqrt(x + sqrt(x)), x)
Bernard Parisse made a workaround commit:
and Giac can now evaluate this integral to:
2*(2*((1/6*sqrt(x)+1/24)*sqrt(x)1/16)*sqrt(x+sqrt(x))1/16*ln(abs(2*(sqrt(x+sqrt(x))sqrt(x))1)))
comment:27 Changed 5 months ago by
So, any comment or review, somebody ?
comment:28 Changed 5 months ago by
 Cc vdelecroix vklein added
Anybody out there ? Does somebody care ?
comment:29 Changed 5 months ago by
What happend here
sage: N(integrate(exp(x^2)*log(x), x, 17, 42)) # rel tol 2e12  2.5657285006962035e127 + 8.88178419700125e16
comment:30 Changed 5 months ago by
More precisely: I don't understand the doctest. The answer is indeed close to 2.5657285006962035e127
but then it makes no sense to have an error tolerance of e12
.
comment:31 Changed 5 months ago by
Actually, I believe more
sage: ComplexBallField(512).integral(lambda x,_: (x*x).exp()*x.log(), 17, 42) [2.565728500561051482917356396130478590014770955402032663e127 +/ 8.94e182]
(which is close to what provide gsl
and the initial doctest) rather than
sage: f = integral(exp(x^2)*log(x), (x,17,42), algorithm='sympy') sage: RealBallField(512)(f) [1.5577642856694726097777e130 +/ 4.52e153]
Is the formula provided by sympy even correct?
comment:32 Changed 5 months ago by
Or the symbolic ring does something completely stupid to evaluate hypergeometric functions with real ball input?
comment:33 Changed 5 months ago by
 Status changed from needs_review to needs_info
comment:34 Changed 5 months ago by
According to Mathematica it is correct... also confirmed by this smaller range computation
sage: f = integral(exp(x^2)*log(x), x, algorithm='sympy') sage: RBF(f.subs(x=2))  RBF(f.subs(x=1)) [0.032612042358609 +/ 8.42e16] sage: CBF.integral(lambda x,_: (x*x).exp()*x.log(), 1, 2) [0.032612042358610 +/ 5.89e16]
The conversion to ball fields appear to be wrong also for much smaller values
sage: prec = 256 sage: R = RealBallField(prec) sage: C = ComplexBallField(prec) sage: a = 5; b = 6 sage: C.integral(lambda x,_: (x*x).exp()*x.log(), a, b) [2.218662695682179212714334629004026381551686783767107654124923915615978e12 +/ 6.31e82] sage: R(f.subs(x=b))  R(f.subs(x=a)) [2.218664132884571616593899629338237856352627658842503518158550169e12 +/ 2.32e76]
The two values start to differ from the 7th digit.
comment:35 Changed 5 months ago by
Conclusion: change the doctest to something more relevant such as
sage: integrate(exp(x^2)*log(x), x) 1/2*sqrt(pi)*erf(x)*log(x)  x*hypergeometric((1/2, 1/2), (3/2, 3/2), x^2)
Though this takes 10 secs!!!
Also taking the numerical approximation of a symbolic integral is the last thing you want to do to evaluate numerically an integral. We should not tell users to do that.
comment:36 Changed 5 months ago by
 Cc mmezzarobba added
Anyway, these doctests do no longer make any sense in the context of the book, because sage is now able to integrate symbolically this expression. In the book, this is given as an example where symbolic integration fails..
I would rather remove them from our testsuite..
comment:37 Changed 5 months ago by
I would rather remove them from our test suite..
But didn't they still want a failing doctest somehow, or is this okay by their point of view? (There are so many authors, maybe you are one of them and can speak to that.) It would be good to somewhere have a failing integral, because that tests whether that branch of the code (i.e. the one that says what to do if everything still doesn't work) is not inadvertently no longer working. That said, I haven't been too aggressive on this ticket so I'll let all of you decide on this.
comment:38 Changed 5 months ago by
As a still failing integral, one could use then
integral(cos(ln(cos(x))),x,0,pi/2)
which cannot be done by maxima, giac, sympy or fricas.
comment:39 followup: ↓ 40 Changed 5 months ago by
As one of the authors, that's fine with me, but I'm not the main author of the chapter, and I think the author with the clearest idea of what to do and not to do with these doctests is Paul (already Cc:ed). Paul, if you are reading this, is it okay with you to remove the test?
comment:40 in reply to: ↑ 39 ; followup: ↓ 41 Changed 5 months ago by
Replying to mmezzarobba:
As one of the authors, that's fine with me, but I'm not the main author of the chapter, and I think the author with the clearest idea of what to do and not to do with these doctests is Paul (already Cc:ed). Paul, if you are reading this, is it okay with you to remove the test?
since it was decided in ticket #27958, comment 58 (while the book was already in print by SIAM), not to accept doctests that do not pass in Python3 (while the book was publicly available since a long time, and nobody asked before for the doctests to be compatible with Python 3), the doctests in Sage do no longer correspond to those in the book, and I have decided not to spend any more time on this.
In summary: do whatever you want, those doctests are no more relevant for the french book.
comment:41 in reply to: ↑ 40 Changed 5 months ago by
Replying to zimmerma:
Replying to mmezzarobba:
As one of the authors, that's fine with me, but I'm not the main author of the chapter, and I think the author with the clearest idea of what to do and not to do with these doctests is Paul (already Cc:ed). Paul, if you are reading this, is it okay with you to remove the test?
There is no comment number 58 on ticket #27958!
comment:42 followup: ↓ 43 Changed 5 months ago by
sorry that was ticket #23572
comment:43 in reply to: ↑ 42 Changed 5 months ago by
Replying to zimmerma:
sorry that was ticket #23572
I did not know about what happend with this ticket. It is indeed unfortunate... the switch to Python 3 is a long and painful route.
But the matter under discussion here is of different nature. The single test under discussion is
N(integrate(exp(x^2)*log(x), x, 17, 42))
There is a positive change of behavior since Sage will now symbolically integrate exp(x^2)*log(x)
. As the doctest was used to illustrate the fact that numerical integration was possible on unevaluated integrals it would not make sense anymore. It is rather the book that has to be adapted here (contrarily to #23572 where the doctest framework should have been adapted).
Isn't there a document that explains what has changed between the available version of the book and the current Sage version? I think this would be very useful.
comment:44 Changed 5 months ago by
 Commit changed from 43c3bfa2c30ca2f406ddaa149f412116ec6fca62 to d4371e3be96a2386f65b57d2e9f36f69eaf71587
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
d4371e3  enhance our integration capacities by using maxima, giac and sympy

comment:45 Changed 5 months ago by
Here is a proposal that keeps the spirit of the modified books' doctests.
comment:46 Changed 5 months ago by
 Status changed from needs_info to needs_review
comment:47 Changed 4 months ago by
 Milestone changed from sage8.8 to sage8.9
Moving tickets from the Sage 8.8 milestone that have been actively worked on in the last six months to the next release milestone (optimistically).
comment:48 Changed 4 months ago by
 Status changed from needs_review to needs_work
Regarding builtins, instead of testing that a string is produced:
sage: log_integral(x)._fricas_init_() 'li(x)'
You could (optionally) test that fricas actually handles the function:
sage: log_integral(x)._fricas_() # optional  fricas li(x)
as you did for giac.
comment:49 Changed 4 months ago by
 Status changed from needs_work to needs_review
Well, I prefer the first test, because Fricas is optional.
edit
I could change that if you insist, but this is really not the core of the matter.
comment:50 Changed 4 months ago by
 Description modified (diff)
 Status changed from needs_review to positive_review
comment:51 Changed 4 months ago by
 Description modified (diff)
The thing is that this test does not test anything, you could change with any string not representing a fricas
function, it will work as well.
Anyway, i see that there are other similar fricas
functions that are not tested, so let us postpone that for a dedicated ticket.
comment:52 Changed 4 months ago by
 Branch changed from u/chapoton/27958 to d4371e3be96a2386f65b57d2e9f36f69eaf71587
 Resolution set to fixed
 Status changed from positive_review to closed
New commits:
enhance our integration capacities by using maxima, giac and sympy