#27958 closed enhancement (fixed)
enhance the integration by using also giac and sympy
Reported by:  Frédéric Chapoton  Owned by:  

Priority:  major  Milestone:  sage8.9 
Component:  symbolics  Keywords:  integration 
Cc:  Markus Wageringel, Jakob Kroeker, KarlDieter Crisman, Peter Bruin, Ralf Stephan, Sébastien Labbé, Samuel Lelièvre, Thierry Monteil, Vincent Delecroix, vklein, Paul Zimmermann, Marcelo Forets, Marc Mezzarobba  Merged in:  
Authors:  Frédéric Chapoton  Reviewers:  Thierry Monteil 
Report Upstream:  N/A  Work issues:  
Branch:  d4371e3 (Commits, GitHub, GitLab)  Commit:  
Dependencies:  Stopgaps: 
Description (last modified by )
Change History (53)
comment:1 Changed 3 years ago by
Branch:  → u/chapoton/27958 

Cc:  Ralf Stephan Samuel Lelièvre added 
Commit:  → 3ddd815e4fbbef38f3bb2b9b1c75674d20c44cd0 
Status:  new → needs_review 
comment:2 Changed 3 years ago by
Commit:  3ddd815e4fbbef38f3bb2b9b1c75674d20c44cd0 → e297e21f5bae8d5b61e7701820bb691857bea647 

Branch pushed to git repo; I updated commit sha1. New commits:
e297e21  trac 27958 polishing the details

comment:3 Changed 3 years ago by
Cc:  Sébastien Labbé KarlDieter Crisman Markus Wageringel Peter Bruin Jakob Kroeker added 

comment:4 Changed 3 years ago by
Cc:  Thierry Monteil added 

Keywords:  integration added 
comment:5 Changed 3 years 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 3 years ago by
Status:  needs_review → needs_work 

some failing doctests, work in progress
comment:7 Changed 3 years ago by
Commit:  e297e21f5bae8d5b61e7701820bb691857bea647 → 4a9fbbb2656d3af8fe015ad00fde43daf71cf0fc 

Branch pushed to git repo; I updated commit sha1. New commits:
4a9fbbb  trac 27958 fixing doctests

comment:8 Changed 3 years ago by
Cc:  Marcelo Forets added 

This seems to work well enough. Waiting for the patchbots' green lights.
comment:9 Changed 3 years 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 3 years ago by
Reviewers:  → Thierry Monteil 

Actually, mathematica gives a very similar result : https://www.wolframalpha.com/input/?i=integrate+log(x)*exp(x%5E2)
comment:11 Changed 3 years ago by
Commit:  4a9fbbb2656d3af8fe015ad00fde43daf71cf0fc → 8eaf4c1c9df9d4097dae5483a911fe71f95f5ab8 

Branch pushed to git repo; I updated commit sha1. New commits:
8eaf4c1  trac 27958 some details in doc

comment:13 Changed 3 years 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 3 years 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 3 years ago by
Cc:  Paul Zimmermann 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 3 years ago by
Commit:  8eaf4c1c9df9d4097dae5483a911fe71f95f5ab8 → 406e27278cc325b9fe0b84b448da6cbad537bba8 

Branch pushed to git repo; I updated commit sha1. New commits:
406e272  trac 27958 doc polishing

comment:17 Changed 3 years ago by
Status:  needs_work → needs_review 

Thx, I have made the suggested changes.
Let us keep the possibility of also adding the fricas integrator for another ticket.
comment:19 Changed 3 years 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 3 years 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 3 years 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 3 years 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 3 years ago by
Commit:  406e27278cc325b9fe0b84b448da6cbad537bba8 → 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:25 followup: 26 Changed 3 years 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 Changed 3 years 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:28 Changed 3 years ago by
Cc:  Vincent Delecroix vklein added 

Anybody out there ? Does somebody care ?
comment:29 Changed 3 years 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 3 years 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 3 years 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 3 years ago by
Or the symbolic ring does something completely stupid to evaluate hypergeometric functions with real ball input?
comment:33 Changed 3 years ago by
Status:  needs_review → needs_info 

comment:34 Changed 3 years 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 3 years 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 3 years ago by
Cc:  Marc Mezzarobba 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 3 years 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 3 years 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 3 years 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 followup: 41 Changed 3 years 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 Changed 3 years 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:43 Changed 3 years 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 3 years ago by
Commit:  43c3bfa2c30ca2f406ddaa149f412116ec6fca62 → 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 3 years ago by
Here is a proposal that keeps the spirit of the modified books' doctests.
comment:46 Changed 3 years ago by
Status:  needs_info → needs_review 

comment:47 Changed 3 years ago by
Milestone:  sage8.8 → 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 3 years ago by
Status:  needs_review → 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 3 years ago by
Status:  needs_work → 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 3 years ago by
Description:  modified (diff) 

Status:  needs_review → positive_review 
comment:51 Changed 3 years 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 3 years ago by
Branch:  u/chapoton/27958 → d4371e3be96a2386f65b57d2e9f36f69eaf71587 

Resolution:  → fixed 
Status:  positive_review → closed 
comment:53 Changed 3 years ago by
Commit:  d4371e3be96a2386f65b57d2e9f36f69eaf71587 

See #28964 for a possible follow up.
New commits:
enhance our integration capacities by using maxima, giac and sympy