Lambdifying a parametric integral leads to integrand symbol not being defined

[notchxu]

using QuadGK

@syms w, t

function integral_with_quad(expr, lims)
    v, a, b = lims
    return quadgk(lambdify(expr, (v,)), a, b)
end

ex = sympy.Integral(t^2, (t,0,w))
fn = lambdify(ex, fns=Dict("Integral" => integral_with_quad))
fn(5)

Attempting to evaluate fn(5) leads to error UndefVarError: t not defined. However, this should not be the case, since the only free symbol in ex is w. The same error occurs if I don't try to replace Integral with the quadgk numerical.

Attempting to fix this by replacing the final two lines with

fn = lambdify(ex, (w,t), fns=Dict("Integral" => integral_with_quad))
fn(5,t)

leads to error MethodError: no method matching iterate(::Symbol) with corresponding stack trace

Stacktrace:
 [1] expr_to_function(body::Expr, vars::Symbol)
   @ SymPy lambdify.jl:242
 [2] lambdify(ex::Sym, vars::Symbol; fns::Dict{Any, Any}, values::Dict{Any, Any}, use_julia_code::Bool, invoke_latest::Bool)
   @ SymPy lambdify.jl:217
 [3] lambdify(ex::Sym, vars::Symbol)
   @ SymPy lambdify.jl:216
 [4] integral_as_quad(f::Sym, lims::Tuple{Symbol, Int64, Symbol})
   @ Main .\REPL[18]:3

[jverzani]

Yeah, the issue is that though t doesn't appear in free_symbols, it does show up when turning the sympy object into an expression:

julia> SymPy.convert_expr(ex)
:(Integral(t ^ 2, (:t, 0, :w)))

julia> free_symbols(ex)
1-element Vector{Sym}:
 w

I'd guess this will only work if you call integrate (which calls doit() on the unevaluated integral expression to actually remove the dummy variable.)

[suburbski]

doit() may solve the issue in this specific case, but you can run into the same problem when the antiderivative has a case distinction which can't be resolved without concrete values. In the following example, x will not appear in free_symbols either, and therefore it yields an error.

x, a, b = symbols("x,a,b", real=true)

res = SymPy.integrate(abs(sin(2π * x)), (x, a, b))
I = SymPy.lambdify(res; use_julia_code=true)
@show res
@show SymPy.free_symbols(res)
@show I(0, 1)

res = Piecewise((Integral(Piecewise((1.0*sin(6.28318530717959*x), sin(6.28318530717959*x) >= 0), (-1.0*sin(6.28318530717959*x), True)), (x, a, b)), a < b), (-Integral(Piecewise((1.0*sin(6.28318530717959*x), sin(6.28318530717959*x) >= 0), (-1.0*sin(6.28318530717959*x), True)), (x, b, a)), True))
SymPy.free_symbols(res) = Sym{PyCall.PyObject}[a, b]
ERROR: LoadError: UndefVarError: `x` not defined

So I've tried to be smart and add x to the final result, which will make it appear in free_symbols again. However, this still doesn't work.

x, a, b = symbols("x,a,b", real=true)

res = x + SymPy.integrate(abs(sin(2π * x)), (x, a, b))
I = SymPy.lambdify(res; use_julia_code=true)
@show res
@show SymPy.free_symbols(res)
@show I(0, 1, 0)

res = x + Piecewise((Integral(Piecewise((1.0*sin(6.28318530717959*x), sin(6.28318530717959*x) >= 0), (-1.0*sin(6.28318530717959*x), True)), (x, a, b)), a < b), (-Integral(Piecewise((1.0*sin(6.28318530717959*x), sin(6.28318530717959*x) >= 0), (-1.0*sin(6.28318530717959*x), True)), (x, b, a)), True))
SymPy.free_symbols(res) = Sym{PyCall.PyObject}[a, b, x]
ERROR: LoadError: UndefVarError: `True` not defined

How do you get around this problem?

[jverzani]

That is a tricky one. It isn't supported in the underlying sympy library either:

julia> sympy.o.julia_code(res.o)
"# Not supported in Julia:\n# Integral\n# Integral\n((a < b) ? (Integral(Piecewise((1.0*sin(6.28318530717959*x), sin(6.28318530717959*x) >= 0), (-1.0*sin(6.28318530717959*x), True)), (x, a, b))) : (-Integral(Piecewise((1.0*sin(6.28318530717959*x), sin(6.28318530717959*x) >= 0), (-1.0*sin(6.28318530717959*x), True)), (x, b, a))))"

The issue is the x shouldn't be part of the conditions, but it is. I'm surprised it isn't in free_symbols given how it prints,

I was hoping this issue would be the one that got me to follow the lead of Symbolics and more flexibly allow the arguments to be specified, but it is some other issue I don't know where to go with,