Add support for Integral with lambdify (using scipy or mpmath)

[ehren]

References to other Issues or PRs

Fixes #19641, Fixes #6408, Fixes #4470, Fixes #4471, Fixes #5932

Brief description of what is fixed or changed

Add support for Integral with lambdify (using scipy or mpmath printers)

Other comments

Release Notes

• utilities
  • Add support for Integral with lambdify (using scipy or mpmath)

[sympy-bot]

✅

Hi, I am the SymPy bot (v160). 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:

• utilities
  • Add support for Integral with lambdify (using scipy or mpmath) (#20134 by @ehren)

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

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

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 #19641, Fixes #6408, Fixes #4470, Fixes #4471, Fixes #5932

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

Add support for Integral with lambdify (using scipy or mpmath printers)

#### Other comments


#### 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 -->

* utilities
   * Add support for Integral with lambdify (using scipy or mpmath)

<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[codecov]

Codecov Report

Merging #20134 into master will increase coverage by 7.894%.
The diff coverage is 88.138%.

@@              Coverage Diff              @@
##            master    #20134       +/-   ##
=============================================
+ Coverage   64.466%   72.360%   +7.894%     
=============================================
  Files          671       671               
  Lines       173707    173803       +96     
  Branches     41020     41048       +28     
=============================================
+ Hits        111982    125765    +13783     
+ Misses       55307     42086    -13221     
+ Partials      6418      5952      -466     

[moorepants]

This looks good to me. Nice addition! +1 merge after all CI passes.

[oscarbenjamin]

Does this work with vectorised inputs?

[ehren]

Does this work with vectorised inputs?

If I understand correctly this seems to be the same question asked about the simpler case of constant functions in issue #5642 where it is desired that when a sympy constant function is lambdified but passed a vector parameter the constant result is returned in an array with the same shape as the parameter passed to the lambdified func. So right now this has the same problem. Adapting an example from that issue:

>>> from sympy import *
>>> import numpy
>>> x, y = symbols("x, y")
>>> u = numpy.linspace(-2, 2, 10)
>>> s = lambdify(y, Integral(x**2, (x, 0, y)))
>>> s(1)
0.33333333333333337
>>> s(u)  # fail
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<lambdifygenerated-4>", line 2, in _lambdifygenerated
  File "/Users/ehren/Documents/symdev/envpycharm1/lib/python3.7/site-packages/scipy/integrate/quadpack.py", line 348, in quad
    flip, a, b = b < a, min(a, b), max(a, b)
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
>>> v = numpy.vectorize(s)  # workaround from #5642
>>> v(u)
array([-2.66666667, -1.25468679, -0.45724737, -0.09876543, -0.00365798,
        0.00365798,  0.09876543,  0.45724737,  1.25468679,  2.66666667])

There's also scipy.integrate.quad_vec for vector valued functions but I believe that's a separate issue than the vectorized input problem (and I suppose there's no unevaluated form of vector_integrate in sympy anyway).

[moorepants]

For the vectorized situation, using the numpy.vectorize() function internal to lambdify seems like a workable approach given that scipy.quad() is not a vectorized function and scipy doesn't offer them. The other option is to raise an error for the user if they try to use vector inputs. numpy.vectorized() is a computationally slow option, but at least it works.

[moorepants]

Another option is using trapz, simps, or romb for vectorized input: https://stackoverflow.com/questions/41223186/numpy-vectorization-with-integration

That also begs the question of why quad is chosen versus another algorithm. If a user specifies a replacement function, it should also work. The lambdify api lets you swap out functions.

[moorepants]

This is an aside, but this is a pretty sweet package: https://github.com/nschloe/quadpy (just discovered it myself).

[ehren]

One issue with the user specified replacement functions is they do have to accept the same signature as the sympy version but I guess it's documented that way (so the user would basically have to repeat the parsing code added here for supporting another numerical routine). Just a check that custom modules work with Integral (this would already work in master):

>>> autodoit = [{'Integral': lambda *args: Integral(*args).doit()}]
>>> f = lambdify(y, Integral(x**2, (x, 0, y)), modules=autodoit)
>>> f(1)
1/3

This PR also sweeps under the rug the additional options that can be passed to mpmath and scipy's quad but we likely have that issue with lambdifying other sympy functions.

Incidentally, found a few more duplicate issues (added to references section).

[moorepants]

This PR also sweeps under the rug the additional options that can be passed to mpmath and scipy's quad but we likely have that issue with lambdifying other sympy functions.

I think that's what the modules option is for.

[ehren]

I think that's what the modules option is for.

It might be nice if you could just pass additional keyword arguments to a module function but I think this could be a maintenance burden so it's best to leave it as is.

[oscarbenjamin]

This looks good. Thanks for fixing this long-standing issue.