Lambdify A Parametric Integral

Question

I have the following issue: I want to lambdify a sympy expression containing parametric integrals like Integral(tanh(a*x),(x,0,1)). I tried to do a manual implementation like here.

What we want is essentially that the integral gets converted to something like:

lambda theta: quad(lambda x: g(x,theta), a,b)[0]

where

g = sp.lambdify((x,param), f, modules='numpy'))

Consider the following MWE:

import sympy as sp
import numpy as np
from scipy.integrate import quad

def integral_as_quad(function, limits):
    x, a, b = limits
    param = function.free_symbols - {x}
    f = sp.lambdify( (x,*param), function, modules='numpy')
    return lambda y: quad(lambda x: f(x,y), a,b)[0]

a, x = sp.symbols('a,x')
I = sp.Integral(sp.tanh(a*x),(x,0,1))
K = integral_as_quad(sp.tanh(a*x),(x,0,1))
L = sp.lambdify(a, I, modules=['numpy', {'Integral':integral_as_quad}] )

Then calling for example K(1) returns the correct value. However L(1) gives

AttributeError: 'Mul' object has no attribute 'tanh'

Does anyone have an idea how to fix this?

NOTE: Doing it manually is no option, since the expressions I deal with are way more complicated and may contain several different integrals. So I really need to get the lambdify working.

score 3 · Accepted Answer · answered Jul 04 '18 at 16:24

I think returning a lambda from integral_as_quad cannot work, because this lambda will never be called, as the Integral object in SymPy is not callable. Instead, the parameter tuple can be passed to quad via its args argument. Another change I made is in the outer lambdification, replacing

modules=['numpy', {'Integral':integral_as_quad}]

with

modules=[{'Integral': integral_as_quad}, 'sympy']

The idea is that at this stage we don't need NumPy functions yet, we just want to replace the Integral by our callable. The order of modules list matters: the dictionary comes first to prevent SymPy from keeping Integral as an Integral.

Now L(1) returns the correct amount.

import sympy as sp
import numpy as np
from scipy.integrate import quad

def integral_as_quad(function, limits):
    x, a, b = limits
    param = tuple(function.free_symbols - {x})
    f = sp.lambdify((x, *param), function, modules=['numpy'])
    return quad(f, a, b, args=param)[0]

a, x = sp.symbols('a,x')
I = sp.Integral(sp.tanh(a*x), (x,0,1))
L = sp.lambdify(a, I, modules=[{'Integral': integral_as_quad}, 'sympy'])

Thank you so much, you are my saviour! And by adding `a = sp.lambdify( (), a, 'numpy')` `b = sp.lambdify( (), b, 'numpy')` `return quad(f, a(), b(), args=param)[0]` we can even deal with Integral boundaries that are `oo` or depend on `param`!!!! — Hyperplane, Jul 04 '18 at 17:15
Do you think it is possible to generalize this approach to the case when one wants to lamdify a general expression that may contain any number of nested parametric integrals? — Hyperplane, Jul 23 '18 at 12:16
I don't know. Post a new question about nested integrals with a minimal example that's not covered by what's written here. — , Jul 23 '18 at 13:01

score 0 · Answer 2 · answered Jul 04 '18 at 14:37

So one possible workaround I have found, but I am unhappy with because it is too slow for my application, is the following:

def to_lambda(expr, param):
    # Preprocessing
    expr = expr.evalf()
    f = sp.lambdify([param], expr, modules='sympy')
    fun = lambda x: np.array(np.array(f(x).evalf()), dtype='float64')
    return fun

So first, expr gets cast into a lambda function using sympy-functions, e.g. we have

f = lambda a: Integral(tanh(a*x),(x,0,1))

and then we use sympy's internal integrator via evalf() (slow!).

Also, don't ask me why theres the double np.array, if one puts the dtype='float64' into the first one, then it returns TypeError: __array__() takes 1 positional argument but 2 were given

Lambdify A Parametric Integral

2 Answers2

Linked