I have a loss function L(u1,...,un) that takes the form
L(u) = Integral( (I**2 + J**2) * h, (t,c,d) ) //h=h(t), I=I(t,u)
where
I = Integral( f, (x,a,b) ) // f=f(x,t,u)
J = Integral( g, (x,a,b) ) // g=g(x,t,u)
I want to numerically minimize L with scipy, hence I need to lambdify the expression.
However at this point in time lambdify does not natively support translating integrals. With some tricks one can get it to work with single parametric integrals, see Lambdify A Parametric Integral. However I don't see how the proposed solution could possibly be extended to this generalised case.
One idea that in principle should work is the following:
Take the computational graph defining L. Recursively, starting from the leaves, replace each symbolic operation with the corresponding numerical operation, expressed as a lambda function. However this would lead to an immense nesting of lambda function, which I suspect has a very bad influence on performance.
Ideally we would want to arrive at the same result as one would by hand crafting:
L = lambda u: quad(lambda t: (quad(lambda x: f,a,b)[0]**2
+ quad(lambda x: g,a,b)[0]**2)*h, c, d)[0]
MWE: (using code from old thread)
from sympy import *
from scipy.integrate import quad
import numpy as np
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]
x,a,b,t = sp.symbols('x a b t')
f = (x/t-a)**2
g = (x/t-b)**2
h = exp(-t)
I = Integral( f**2,(x,0,1))
J = Integral( g**2,(x,0,1))
K = Integral( (I**2 + J**2)*h,(t,1,+oo))
F = lambdify( (a,b), K, modules=['sympy',{"Integral": integral_as_quad}])
L = lambda a,b: quad(lambda t: (quad(lambda x: (x/t-a)**2,0,1)[0]**2
+ quad(lambda x: (x/t-b)**2,0,1)[0]**2)*np.exp(-t), 1, np.inf)[0]
display(F(1,1))
display(type(F(1,1)))
display(L(1,1))
