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: sage-8.9
Component: symbolics Keywords: integration
Cc: gh-mwageringel, 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:

Change History (52)

comment:1 Changed 5 months ago by chapoton

  • Branch set to u/chapoton/27958
  • Cc rws slelievre added
  • Commit set to 3ddd815e4fbbef38f3bb2b9b1c75674d20c44cd0
  • Status changed from new to needs_review

New commits:

3ddd815enhance our integration capacities by using maxima, giac and sympy

comment:2 Changed 5 months ago by git

  • Commit changed from 3ddd815e4fbbef38f3bb2b9b1c75674d20c44cd0 to e297e21f5bae8d5b61e7701820bb691857bea647

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

e297e21trac 27958 polishing the details

comment:3 Changed 5 months ago by chapoton

  • Cc slabbe kcrisman gh-mwageringel pbruin jakobkroeker added

comment:4 Changed 5 months ago by chapoton

  • Cc tmonteil added
  • Keywords integration added

comment:5 Changed 5 months ago by chapoton

sometimes sympy takes ages, so that's why it comes last:

sage: integrate(sqrt(1-m*sin(x)^2),x,algorithm='sympy')
... hangs ...

comment:6 Changed 5 months ago by chapoton

  • Status changed from needs_review to needs_work

some failing doctests, work in progress

comment:7 Changed 5 months ago by git

  • Commit changed from e297e21f5bae8d5b61e7701820bb691857bea647 to 4a9fbbb2656d3af8fe015ad00fde43daf71cf0fc

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

4a9fbbbtrac 27958 fixing doctests

comment:8 Changed 5 months ago by chapoton

  • Cc mforets added

This seems to work well enough. Waiting for the patchbots' green lights.

comment:9 Changed 5 months ago by tmonteil

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 ?

Last edited 5 months ago by tmonteil (previous) (diff)

comment:10 Changed 5 months ago by tmonteil

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

  • Commit changed from 4a9fbbb2656d3af8fe015ad00fde43daf71cf0fc to 8eaf4c1c9df9d4097dae5483a911fe71f95f5ab8

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

8eaf4c1trac 27958 some details in doc

comment:12 Changed 5 months ago by chapoton

thx, I have modified this part of the doc

comment:13 Changed 5 months ago by tmonteil

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.

Last edited 5 months ago by tmonteil (previous) (diff)

comment:14 Changed 5 months ago by tmonteil

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 tmonteil

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

  • Commit changed from 8eaf4c1c9df9d4097dae5483a911fe71f95f5ab8 to 406e27278cc325b9fe0b84b448da6cbad537bba8

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

406e272trac 27958 doc polishing

comment:17 Changed 5 months ago by chapoton

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

full green bot, please review

comment:19 Changed 5 months ago by kcrisman

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 mforets

+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 chapoton

as I said in comment:5, sometimes sympy seems to hang forever. Here is one case :

sage: %time integrate(sqrt(1-3*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(1-3*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(1-3*sin(x)^2),x,algorithm='sympy')
NOTHING HAPPENS...

So this behaviour is now inherited by our integration.

One can hope that the users will CTRL-C. Or use a time-out ??

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(1-3*sin(x)^2),x)
CPU times: user 1min 47s, sys: 346 ms, total: 1min 47s
Wall time: 1min 50s
elliptic_e(x, 3)
Last edited 5 months ago by chapoton (previous) (diff)

comment:22 Changed 5 months ago by kcrisman

I'd be careful, because not every user perhaps knows about Ctrl-C, 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 git

  • Commit changed from 406e27278cc325b9fe0b84b448da6cbad537bba8 to 43c3bfa2c30ca2f406ddaa149f412116ec6fca62

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

43c3bfaenhance our integration capacities by using maxima, giac and sympy

comment:24 Changed 5 months ago by chapoton

squashed branch, with more doctests showing the new capacities

comment:25 follow-up: Changed 5 months ago by chapoton

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..
Last edited 5 months ago by chapoton (previous) (diff)

comment:26 in reply to: ↑ 25 Changed 5 months ago by slelievre

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)))
Last edited 5 months ago by slelievre (previous) (diff)

comment:27 Changed 5 months ago by chapoton

So, any comment or review, somebody ?

comment:28 Changed 5 months ago by chapoton

  • Cc vdelecroix vklein added

Anybody out there ? Does somebody care ?

comment:29 Changed 5 months ago by vdelecroix

What happend here

   sage: N(integrate(exp(-x^2)*log(x), x, 17, 42))  # rel tol 2e-12
-  2.5657285006962035e-127
+  -8.88178419700125e-16
Last edited 5 months ago by vdelecroix (previous) (diff)

comment:30 Changed 5 months ago by vdelecroix

More precisely: I don't understand the doctest. The answer is indeed close to 2.5657285006962035e-127 but then it makes no sense to have an error tolerance of e-12.

comment:31 Changed 5 months ago by vdelecroix

Actually, I believe more

sage: ComplexBallField(512).integral(lambda x,_: (-x*x).exp()*x.log(), 17, 42)
[2.565728500561051482917356396130478590014770955402032663e-127 +/- 8.94e-182]

(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.5577642856694726097777e-130 +/- 4.52e-153]

Is the formula provided by sympy even correct?

comment:32 Changed 5 months ago by vdelecroix

Or the symbolic ring does something completely stupid to evaluate hypergeometric functions with real ball input?

Last edited 5 months ago by vdelecroix (previous) (diff)

comment:33 Changed 5 months ago by vdelecroix

  • Status changed from needs_review to needs_info

comment:34 Changed 5 months ago by vdelecroix

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.42e-16]
sage: CBF.integral(lambda x,_: (-x*x).exp()*x.log(), 1, 2)
[0.032612042358610 +/- 5.89e-16]

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.218662695682179212714334629004026381551686783767107654124923915615978e-12 +/- 6.31e-82]
sage: R(f.subs(x=b)) - R(f.subs(x=a))
[2.218664132884571616593899629338237856352627658842503518158550169e-12 +/- 2.32e-76]

The two values start to differ from the 7th digit.

comment:35 Changed 5 months ago by vdelecroix

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 chapoton

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

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 chapoton

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 follow-up: Changed 5 months ago by 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?

comment:40 in reply to: ↑ 39 ; follow-up: Changed 5 months ago by 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?

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 vdelecroix

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?

[...] in ticket #27958, comment 58 [...]

There is no comment number 58 on ticket #27958!

comment:42 follow-up: Changed 5 months ago by zimmerma

sorry that was ticket #23572

comment:43 in reply to: ↑ 42 Changed 5 months ago by vdelecroix

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 git

  • Commit changed from 43c3bfa2c30ca2f406ddaa149f412116ec6fca62 to d4371e3be96a2386f65b57d2e9f36f69eaf71587

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

d4371e3enhance our integration capacities by using maxima, giac and sympy

comment:45 Changed 5 months ago by chapoton

Here is a proposal that keeps the spirit of the modified books' doctests.

comment:46 Changed 5 months ago by chapoton

  • Status changed from needs_info to needs_review

comment:47 Changed 4 months ago by embray

  • Milestone changed from sage-8.8 to sage-8.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 tmonteil

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

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

Last edited 4 months ago by chapoton (previous) (diff)

comment:50 Changed 4 months ago by tmonteil

  • Description modified (diff)
  • Status changed from needs_review to positive_review

comment:51 Changed 4 months ago by tmonteil

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

  • Branch changed from u/chapoton/27958 to d4371e3be96a2386f65b57d2e9f36f69eaf71587
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.