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,

[clouds56]

SymPy.SymPyCore.walk_expression(ex, fns=Dict("Integral" => integral_with_quad))

the problem here is why it returns

:((integral_with_quad)(t ^ 3, (:t, 0, :w)))

instead of

:((integral_with_quad)(t ^ 3, (t, 0, w)))

The walk_expression/convert_expr doesn't involves free symbols.

Shouldn't it be Expr(:tuple, :t, 0, :w)?
related codes here:
https://github.com/jverzani/SymPyCore.jl/blame/7bc3029fcba002adff08ff1cecb9fb1f07a514fb/src/lambdify.jl#L253

     elseif fn == "Tuple"
-        return walk_expression.(args(ex), values=values, fns=fns)
+        return Expr(:tuple, walk_expression.(args(ex), values=values, fns=fns)...)

I'm not sure where some other place would use this fn == "Tuple"

[jverzani]

Thanks, I'll try and get to this soon.

[jverzani]

Thanks. With your change (jverzani/SymPyCore.jl#79) this almost works as expected. The next issue is the dummy variable t gets caught in the expression leading to an error unless the expression is invoked in the right context. This pattern seems to work now with your fix:

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))
fn0 = lambdify(ex, fns=Dict("Integral" => integral_with_quad), invoke_latest=false)
fn = eval(fn0)
fn(5)  # (41.66666666666666, 0.0)

[clouds56]

Giving that SymPy accept things like Integral(t^2, (t, 0, t)) where inner t differs from outer t, it might be tricky how walk_expression would handle bounded variables.
We might have to add some support for Integral in walk_expression, that converts Integral(t^2, (t, 0, t)) to things like

:(function var"##236"(t)
      begin
              local var"##235" = Sym("##235")
              (integral_with_quad)((var"##235" ^ 2, (var"##235", 0, t))
          end
  end)

considering it would call Sym every time the function would be invoked, it might be even better to gen code like (according to JuliaLang/julia#20353)

let var"##235" = Sym("##235")
  function var"##236"(t)
       (integral_with_quad)((var"##235" ^ 2, (var"##235", 0, t))
  end
end

While the latter might introduce some breaking change to walk_expression/convert_expr

I'm glad to help here.