Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use the faster libgiac interface for "giac" integration #38686

Open
wants to merge 9 commits into
base: develop
Choose a base branch
from

Conversation

orlitzky
Copy link
Contributor

While working on #38668, I noticed that all of our examples that use giac for integration do so with algorithm='giac'. This uses the pexpect interface that is much slower (and simply sloppier) than the library interface that would be used with algorithm='libgiac'.

To fix this, I could just change those examples to algorithm='libgiac', but since the libgiac interface is better in every way, that seems kind of rude to me. Most users aren't going to know what lib-anything is, and if they just want to use giac, they're going to try algorithm='giac'. So instead this PR makes algorithm='giac' do the right thing, and use the better interface. I've left the algorithm='libgiac' alone for now to avoid breaking any code.

Afterwards I have removed the pexpect giac integrator entirely, because there's no reason to use it, and no other code in sagelib uses it.

Copy link

github-actions bot commented Sep 20, 2024

Documentation preview for this PR (built with commit 8db21a2; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

@orlitzky
Copy link
Contributor Author

AFAICS the CI failure is an unrelated timeout:

2024-09-20T15:19:09.1450225Z ----------------------------------------------------------------------
2024-09-20T15:19:09.1451853Z sage -t --long --random-seed=286735480429121101562228604801325644303 src/sage/data_structures/bounded_integer_sequences.pyx  # Timed out
2024-09-20T15:19:09.1453353Z ----------------------------------------------------------------------

@mkoeppe
Copy link

mkoeppe commented Sep 29, 2024

This looks like a good idea, thanks.

@orlitzky
Copy link
Contributor Author

orlitzky commented Oct 2, 2024

I didn't see the label change, thanks for the review.

@orlitzky orlitzky mentioned this pull request Oct 2, 2024
2 tasks
@vbraun
Copy link
Member

vbraun commented Oct 4, 2024

I'm getting

sage -t --long --warn-long 45.2 --random-seed=123 src/sage/symbolic/integration/integral.py
**********************************************************************
File "src/sage/symbolic/integration/integral.py", line 64, in sage.symbolic.integration.integral.IndefiniteIntegral.__init__
Failed example:
    integrate(Ex, x, algorithm='giac')  # long time
Expected:
    4*(-2*x^(1/3) + 1)^(3/4) + 6*arctan((-2*x^(1/3) + 1)^(1/4)) - 3*log((-2*x^(1/3) + 1)^(1/4) + 1) + 3*log(abs((-2*x^(1/3) + 1)^(1/4) - 1))
Got:
    4*(-2*x^(1/3) + 1)^(3/4) + 6*arctan((-2*x^(1/3) + 1)^(1/4)) - 3*log(abs((-2*x^(1/3) + 1)^(1/4) + 1)) + 3*log(abs((-2*x^(1/3) + 1)^(1/4) - 1))
**********************************************************************
File "src/sage/symbolic/integration/integral.py", line 209, in sage.symbolic.integration.integral.DefiniteIntegral.__init__
Failed example:
    integral(1/max_symbolic(x, 1)**2, x, 0, oo, algorithm='giac')
Exception raised:
    Traceback (most recent call last):
      File "/home/release/Sage/src/sage/doctest/forker.py", line 714, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/release/Sage/src/sage/doctest/forker.py", line 1144, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.symbolic.integration.integral.DefiniteIntegral.__init__[4]>", line 1, in <module>
        integral(Integer(1)/max_symbolic(x, Integer(1))**Integer(2), x, Integer(0), oo, algorithm='giac')
      File "/home/release/Sage/src/sage/misc/functional.py", line 785, in integral
        return x.integral(*args, **kwds)
               ^^^^^^^^^^^^^^^^^^^^^^^^^
      File "sage/symbolic/expression.pyx", line 13222, in sage.symbolic.expression.Expression.integral
        return integral(self, *args, **kwds)
      File "/home/release/Sage/src/sage/symbolic/integration/integral.py", line 1073, in integrate
        return integrator(expression, v, a, b)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      File "/home/release/Sage/src/sage/symbolic/integration/external.py", line 259, in libgiac_integrator
        result = libgiac.integrate(Pygen(expression), v, a, b)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      File "sage/libs/giac/giac.pyx", line 1934, in sage.libs.giac.giac.GiacFunction.__call__
        return Pygen.__call__(self, *args)
      File "sage/libs/giac/giac.pyx", line 1109, in sage.libs.giac.giac.Pygen.__call__
        raise
      File "sage/libs/giac/giac.pyx", line 1099, in sage.libs.giac.giac.Pygen.__call__
        result = self.gptr[0](right.gptr[0], context_ptr)
    RuntimeError: Unable to check sign: 1>Infinity
**********************************************************************
File "src/sage/symbolic/integration/integral.py", line 697, in sage.symbolic.integration.integral.integrate
Failed example:
    integrate(abs(cos(x)), x, 0, 2*pi, algorithm='giac')
Expected:
    4
Got:
    Warning, integration of abs or sign assumes constant sign by intervals (correct if the argument is real):
    Check [abs(cos(sageVARx))]
    Discontinuities at zeroes of cos(sageVARx) were not checked
    4
**********************************************************************

…egration

The library interface to libgiac is much more efficient than the
pexpect one, but the name algorithm="giac" is more attractive
(especially to newcomers) than algorithm="libgiac". We should just
Do The Right Thing and use libgiac when "giac" is requested.
The pexpect integrator is unused now that the "giac" integrator uses
libgiac. We also add a few "# needs sage.libs.giac" tags to the tests
that use the libgiac integrator, because they obviously will fail when
libgiac is not there.
One example in this file integrates with algorithm='giac', which can
fail if sage.libs.giac is not available to perform the integration.
One example in this file integrates with algorithm='giac', which can
fail if sage.libs.giac is not available to perform the integration.
…tegration

A few examples in this file integrate with algorithm='giac', which can
fail if sage.libs.giac is not available to perform the integration.
There's one test in this file for an integral where the "libgiac"
algorithm prints a warning, but the old "giac" algorithm did not. Now
that algorithm="giac" uses libgiac, we expect the warning in that case
too.
One "giac" integral in this file has changed its output slightly now
that the "giac" algorithm uses libgiac as opposed to the pexpect
interface. The difference between the new and old expected outputs
full_simplify() to zero, so this one is nothing to worry about.
When integrating to/from plus/minus infinity, Giac can get stuck:

  sage: integral(1/max_symbolic(x, 1)**2, x, 0, oo, algorithm='giac')
  ...
  RuntimeError: Unable to check sign: 1>Infinity

This is apparently due to the type of infinity that is used for
integration. There are (at least?) two kinds:

  sage: from sage.libs.giac import libgiac
  sage: libgiac(Infinity)
  Infinity
  sage: libgiac(Infinity._giac_())
  +infinity

And, for integration, I guess we want the second kind, because it
allows the integration to proceed with a warning. This commit calls
_giac_() on the endpoints of a definite integral before passing them
to libgiac.integrate(). It fixes the example above, and we update the
corresponding doctest.
With the recent change to the handling of infinity in the libgiac
integrator, a warning about singular points has vanished.
@orlitzky
Copy link
Contributor Author

orlitzky commented Oct 6, 2024

Expected:
4*(-2x^(1/3) + 1)^(3/4) + 6arctan((-2x^(1/3) + 1)^(1/4)) - 3log((-2x^(1/3) + 1)^(1/4) + 1) + 3log(abs((-2x^(1/3) + 1)^(1/4) - 1))
Got:
4
(-2x^(1/3) + 1)^(3/4) + 6arctan((-2x^(1/3) + 1)^(1/4)) - 3log(abs((-2x^(1/3) + 1)^(1/4) + 1)) + 3log(abs((-2*x^(1/3) + 1)^(1/4) - 1))

These are the same anwer. I updated the expected output.

Expected:
4
Got:
Warning, integration of abs or sign assumes constant sign by intervals (correct if the argument is real):
Check [abs(cos(sageVARx))]
Discontinuities at zeroes of cos(sageVARx) were not checked
4

Nothing too scary here. This doctest already looks for the warning output on the previous line, where integrate is used without an algorithm argument. I tweaked the description a bit and added some ....

File "src/sage/symbolic/integration/integral.py", line 209, in sage.symbolic.integration.integral.DefiniteIntegral.init
Failed example:
integral(1/max_symbolic(x, 1)**2, x, 0, oo, algorithm='giac')
...
RuntimeError: Unable to check sign: 1>Infinity

This one on the other hand is a treat. It's reproducible using libgiac_integrator, but not through integrate(), which (if maxima fails) calls libgiac_integrator with exactly the same arguments. My guess is hidden state within sage.libs.giac.

In any case, I was able to work around it. Don't ask me what the difference is, but giac likes +infinity better than Infinity. You get the latter passing Sage's infinity directly to libgiac, and the former by calling _giac_() on Sage's infinity. So I added a few lines to the integrator to do that when the endpoints are plus/minus infinity.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants