Reuse pole detection/reintegration logic from heurisch_wrapper in meijerint #20270

ehren · 2020-10-16T03:35:33Z

References to other Issues or PRs

Brief description of what is fixed or changed

Moves the code in heurisch_wrapper to a shared function, pole_reintegrate that can also be used by meijerint (via a new meijerint_indefinite_wrapper.

Other comments

reused heurisch_wrapper logic originally from #1622

Release Notes

integrals
- meijerg=True option of integrate now returns a Piecewise handling cases where the indefinite integral returned by meijerint_indefinite could have a zero denominator. This affects any uses of the integrals module that rely on meijerint. This behavior can be disabled using the integrate option conds='none'.

…ntegrate (since we're already checking for poles in the denominator); this avoids unnecessary calls to denoms/solve.

…tion that can be shared by heurisch and meijerg

…of a nan by solve when handling a Piecewise (see also sympygh-20269). Because we're calling solve on an equation derived from the integration result anyway, it should be safe to ignore any special cases arrising from poles in the denominator of this integration result.

sympy-bot · 2020-10-16T03:35:43Z

✅

Hi, I am the SymPy bot (v161). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

integrals
- meijerg=True option of integrate now returns a Piecewise handling cases where the indefinite integral returned by meijerint_indefinite could have a zero denominator. This affects any uses of the integrals module that rely on meijerint. This behavior can be disabled using the integrate option conds='none'. (#20270 by @ehren)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.8.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes #5949

#### Brief description of what is fixed or changed

Moves the code in `heurisch_wrapper` to a shared function, `pole_reintegrate` that can also be used by `meijerint` (via a new `meijerint_indefinite_wrapper`.

#### Other comments

reused heurisch_wrapper logic originally from #1622

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->

* integrals
  * `meijerg=True` option of `integrate` now returns a `Piecewise` handling cases where the indefinite integral returned by `meijerint_indefinite` could have a zero denominator. This affects any uses of the `integrals` module that rely on `meijerint`. This behavior can be disabled using the `integrate` option `conds='none'`.

<!-- END RELEASE NOTES -->

ehren · 2020-10-16T03:43:42Z

Regarding the change to ode:

When solving this:

sympy/sympy/solvers/ode/tests/test_ode.py

Line 241 in 87a885b

eq3 = (Eq(diff(x(t),t), y(t)*x(t)), Eq(diff(y(t),t), x(t)**3))

This integral comes up:

Integral(1/(u*sqrt(C1 + 6*u**3)), u).doit()

which after these changes results in

Piecewise((-2*asinh(sqrt(6)*sqrt(C1)/(6*u**(3/2)))/(3*sqrt(C1)), Ne(C1, 0)), (-sqrt(6)/(9*sqrt(u**3)), True))

Then calling solve on an equation derived from this:

solve(-C2 - t - 3*Piecewise((-2*asinh(sqrt(6)*sqrt(C1)/(6*u**(3/2)))/(3*sqrt(C1)), Ne(C1, 0)), (-sqrt(6)/(9*sqrt(u**3)), True)))

results in the solution to the ode involving nan

[Piecewise((6**(2/3)/(6*(sinh(sqrt(C1)*(C2 + t)/2)/sqrt(C1))**(2/3)), Ne(C1, 0)), (nan, True)), Piecewise((2**(1/3)*3**(2/3)*((C2 + t)**(-2))**(1/3)/3, Eq(C1, 0)), (nan, True)), Piecewise((-2**(1/3)*3**(2/3)*((C2 + t)**(-2))**(1/3)/6 - 2**(1/3)*3**(1/6)*I*((C2 + t)**(-2))**(1/3)/2, Eq(C1, 0)), (nan, True)), Piecewise((-2**(1/3)*3**(2/3)*((C2 + t)**(-2))**(1/3)/6 + 2**(1/3)*3**(1/6)*I*((C2 + t)**(-2))**(1/3)/2, Eq(C1, 0)), (nan, True))]

I considered just changing the definition of the constant to C1 = Symbol("C1", zero=False) but this could rule out valid trivial solutions. I think because the result is then used in an equation passed to solve ignoring the Piecewise solution of integrate is perfectly safe.

asmeurer · 2020-10-16T08:12:09Z

I don't think it's a good idea to put assumptions on the constants introduced by the ode solver. That would be very opaque to end-users.

asmeurer · 2020-10-16T08:12:50Z

Should we be using this wrapper in manualintegrate too?

jksuom · 2020-10-16T09:06:58Z

How about implementing a wrapper like this in the main code? Different integrators could be used in special cases.

asmeurer · 2020-10-16T18:10:08Z

sympy/integrals/meijerint.py

+    >>> from sympy import sin
+    >>> from sympy.abc import x, n, k
+    >>> meijerint_indefinite_wrapper(sin(x**(n*k + 1)), x)
+    Piecewise((x*x**(k*n)*x**(k*n/(k*n + 1))*x**(1/(k*n + 1))*gamma(1/2 + 1/(2*(k*n + 1)))*hyper((1/2 + 1/(2*(k*n + 1)),), (3/2, 3/2 + 1/(2*(k*n + 1))), -x**2*x**(2*k*n)/4)/(2*k*n*gamma(3/2 + 1/(2*(k*n + 1))) + 2*gamma(3/2 + 1/(2*(k*n + 1)))), Ne(k, -1/n)), (x*sin(1), True))


I believe you can line wrap this and the doctest will still pass.

asmeurer · 2020-10-16T18:16:17Z

sympy/integrals/heurisch.py

+
+    def reintegrator(f, x):
+        expr = heurisch(f, x, rewrite, hints, mappings, retries, degree_offset,
+                        unnecessary_permutations, _try_heurisch)


I'm not sure why any of the integration algorithms should reuse only that algorithm for the recursive integration, rather than integrate(). Does anything break for heurisch if we change this? Heurisch in particular is the slowest algorithm, which is always tried last, so it's better to avoid using it if we can. It's also unclear to me that all these arguments to heurisch should be the same for the new substituted integrand.

The only place where it might make sense to keep the same algorithm would be Risch, if we can be sure the "simplified" integral is also of the form that Risch can handle (most likely this is the case, though I'd have to look at the specific places where the divisions happen to be sure). Even then, though, Risch has many NotImplementedError holes, so it still might not make sense there.

asmeurer · 2020-10-16T18:18:29Z

sympy/integrals/integrals.py

+
+
+def pole_reintegrate(f, x, antideriv,
+                     reintegrator=lambda f, x: integrate(f, x, conds="none")):


Another concern is that if someone passes a specific algorithm to integrate, like integrate(meijerg=True), should the reintegrator also only use that algorithm? I'm unsure. It does seem to make sense for the flag to be passed through, but it could lead to the confusing situation where integrate(f, x) is computed with, say, meijerg, but integrate(f, x, meijerg=True) doesn't produce the same result, because meijerg isn't what is actually used by the default reintegrator.

asmeurer · 2020-10-16T18:19:55Z

sympy/integrals/integrals.py

+    Takes an integrand f with integration variable x with an initially
+    computed antiderivative and checks for poles in the denominator. For each
+    of these poles, the integral is reevaluated, and the final integration
+    result is given in terms of a Piecewise.


Can you explain the arguments here, or give an example?

asmeurer · 2020-10-16T18:23:55Z

I'm not sure if you already explained this, but what is the advantage of having a function that checks for denominators rather than hooking into the specific place(s) in the integration code where a division happens? It seems to me that the latter would be preferred if it is possible, as it would avoid unnecessarily splitting out cases where a constant also appears in the denominator of the integrand. I see that the helper function here tries to handle those cases, but it still requires extra work to detect them (if I am reading the code correctly, there is a solve() call for every denominator).

asmeurer · 2020-10-16T18:25:47Z

sympy/integrals/integrals.py

+    slns0 = []
+    for d in denoms(f):
+        try:
+            slns0 += solve(d, dict=True, exclude=(x,))


Why do we have to solve here too? Can't we filter out the common denominators from f without calling solve? That would save a lot of unnecessary work if f has a lot of complicated symbolic constant expressions in denominators, which are unrelated to the integration variable. It's also possible for such expressions to not be solvable, in which case they would be filtered if I am reading this code correctly.

This comes from #2112 but we can just check if d.subs(already_computed_solution) == 0 but still looking at optimizations. Have a few changes coming regarding the comments.

ehren · 2020-10-20T01:45:25Z

How about implementing a wrapper like this in the main code? Different integrators could be used in special cases.

Should we be using this wrapper in manualintegrate too?

I experimented with wrapping the entire _eval_integral but found this is unnecessary with manualintegrate. I seem to have a fix for meijerg that follows similar logic to manualintegrate and doesn't require solve or denoms etc. Still have a few lingering test cases to fix but will send an update.

czgdp1807 · 2020-11-20T17:03:08Z

@ehren Some tests seem to fail. Would you like to debug these? Please let us know if you are still working on it. Thanks for your contributions.

ehren added 6 commits October 15, 2020 17:25

reuse initial result in heurisch wrapper

5173739

use conds='none' in heurisch_wrapper when falling back to top-level i…

0787873

…ntegrate (since we're already checking for poles in the denominator); this avoids unnecessary calls to denoms/solve.

move code in heurish_wrapper (originally from sympygh-1622) to a func…

f33dd63

…tion that can be shared by heurisch and meijerg

add meijerint_indefinite_wrapper similar to heurisch_wrapper

dd049f0

update existing tests + new test for issue

979058f

ehren requested a review from oscarbenjamin as a code owner October 16, 2020 03:35

This was referenced Oct 16, 2020

could solve(Piecewise((x, Eq(y, 0)), (1/x, True)), x, y) avoid nan? #20269

Open

[don't merge] handle a potential zero denominator in a substitution in _meijerint_indefinite_1 #20227

Closed

flake8: remove duplicate import

4709634

asmeurer reviewed Oct 16, 2020

View reviewed changes

oscarbenjamin added the integrals label Oct 17, 2020

github-actions bot added the Merge conflict label Oct 26, 2022

ehren mentioned this pull request Apr 15, 2024

IndexError: Index out of range: calling integrate in sympy 1.12 #26504

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Reuse pole detection/reintegration logic from heurisch_wrapper in meijerint #20270

Reuse pole detection/reintegration logic from heurisch_wrapper in meijerint #20270

ehren commented Oct 16, 2020

sympy-bot commented Oct 16, 2020

ehren commented Oct 16, 2020 •

edited

asmeurer commented Oct 16, 2020

asmeurer commented Oct 16, 2020

jksuom commented Oct 16, 2020

asmeurer Oct 16, 2020

asmeurer Oct 16, 2020

asmeurer Oct 16, 2020

asmeurer Oct 16, 2020

asmeurer commented Oct 16, 2020

asmeurer Oct 16, 2020

ehren Oct 17, 2020

ehren commented Oct 20, 2020

czgdp1807 commented Nov 20, 2020



		def pole_reintegrate(f, x, antideriv,
		reintegrator=lambda f, x: integrate(f, x, conds="none")):

Reuse pole detection/reintegration logic from heurisch_wrapper in meijerint #20270

Are you sure you want to change the base?

Reuse pole detection/reintegration logic from heurisch_wrapper in meijerint #20270

Conversation

ehren commented Oct 16, 2020

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

sympy-bot commented Oct 16, 2020

ehren commented Oct 16, 2020 • edited

asmeurer commented Oct 16, 2020

asmeurer commented Oct 16, 2020

jksuom commented Oct 16, 2020

asmeurer Oct 16, 2020

Choose a reason for hiding this comment

asmeurer Oct 16, 2020

Choose a reason for hiding this comment

asmeurer Oct 16, 2020

Choose a reason for hiding this comment

asmeurer Oct 16, 2020

Choose a reason for hiding this comment

asmeurer commented Oct 16, 2020

asmeurer Oct 16, 2020

Choose a reason for hiding this comment

ehren Oct 17, 2020

Choose a reason for hiding this comment

ehren commented Oct 20, 2020

czgdp1807 commented Nov 20, 2020

ehren commented Oct 16, 2020 •

edited