{autolink-concat}
::::{margin} :::{card} Symbolic integral TR-016 ^^^ This report investigates how to formulate a symbolic integral that correctly evaluates to +++ To be implemented ::: ::::
%pip install -q ampform==0.14.10 black==23.12.1 scipy==1.12.0 sympy==1.12
import inspect
import black
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from ampform.io import aslatex
from ampform.sympy import unevaluated
from IPython.display import Markdown, Math
from scipy.integrate import quad, quad_vec
from sympy.printing.pycode import _unpack_integral_limits
(quad)=
quad()
function¶SciPy's {func}scipy.integrate.quad
cannot integrate complex-valued functions:
def integrand(x):
return x * (x + 1j)
quad(integrand, 0.0, 2.0)
A proposed solution is to wrap the {func}~scipy.integrate.quad
function in a special integrate function that integrates the real and imaginary part of a function separately:
def complex_integrate(func, a, b, **quad_kwargs):
def real_func(x):
return func(x).real
def imag_func(x):
return func(x).imag
real_integral, real_integral_err = quad(real_func, a, b, **quad_kwargs)
imag_integral, imag_integral_err = quad(imag_func, a, b, **quad_kwargs)
return (
real_integral + 1j * imag_integral,
real_integral_err**2 + 1j * imag_integral_err,
)
complex_integrate(integrand, 0.0, 2.0)
:::{warning}
The handling of uncertainties is incorrect.
:::
(quad_vec)=
quad_vec()
function¶The easiest solution, however, seems to be {func}scipy.integrate.quad_vec
:
quad_vec(integrand, 0.0, 2.0)
This has the added benefit that it can handle functions that return arrays:
def gaussian(x, mu, sigma):
return np.exp(-((x - mu) ** 2) / (2 * sigma**2)) / (sigma * np.sqrt(2 * np.pi))
mu_values = np.linspace(-2, +3, num=10)
result, _ = quad_vec(lambda x: gaussian(x, mu_values, sigma=0.5), 0, 2.0)
result
quadpy
¶:::{warning}
quadpy
now requires a license. The examples below are only shown for documentation purposes.
:::
Alternatively, one could use quadpy
, which essentially does the same as in quad()
, but can also (to a large degree) handle vectorized input and properly handles uncertainties. For example:
from functools import partial
def parametrized_func(s_prime, s):
return s_prime * (s_prime + s + 1j)
s_array = np.linspace(-1, 1, num=10)
quadpy.quad(
partial(parametrized_func, s=s_array),
a=0.0,
b=2.0,
)
:::{note}
One may need to play around with the tolerance if the function is not smooth, see sigma-py/quadpy#255.
:::
:::{tip}
quadpy
raises exceptions with {obj}ModuleNotFoundError
s that are a bit unreadable. They are caused by the fact that orthopy
and ndim
need to be installed separately.
:::
The dispersion integral from Eq. {eq}dispersion-integral
in TR-003 features a variable $s$ that is an argument to the function $\Sigma_a$. This becomes a challenge when $s$ gets vectorized (in this case: gets an event-wise {obj}numpy.array
of invariant masses). It seems that quad_vec()
can handle this well though.
def parametrized_func(s_prime, s):
return s_prime * (s_prime + s + 1j)
s_array = np.linspace(-1, +1, num=10)
quad_vec(
lambda x: parametrized_func(x, s=s_array),
a=0.0,
b=2.0,
)
We now attempt to design SymPy expression classes that correctly {func}~sympy.utilities.lambdify.lambdify
using this vectorized numerical integral for handles complex values. Note that this integral expression class derives from {class}sympy.Integral <sympy.integrals.integrals.Integral>
and that:
~sympy.core.basic.Basic.doit
method, so that the integral cannot be evaluated by SymPy,quadpy.quad()
,quad_vec()
,class UnevaluatableIntegral(sp.Integral):
abs_tolerance = 1e-5
rel_tolerance = 1e-5
limit = 50
dummify = True
def doit(self, **hints):
args = [arg.doit(**hints) for arg in self.args]
return self.func(*args)
def _numpycode(self, printer, *args):
integration_vars, limits = _unpack_integral_limits(self)
if len(limits) != 1 or len(integration_vars) != 1:
msg = f"Cannot handle {len(limits)}-dimensional integrals"
raise ValueError(msg)
x = integration_vars[0]
a, b = limits[0]
expr = self.args[0]
if self.dummify:
dummy = sp.Dummy()
expr = expr.xreplace({x: dummy})
x = dummy
integrate_func = "quad_vec"
printer.module_imports["scipy.integrate"].add(integrate_func)
return (
f"{integrate_func}(lambda {printer._print(x)}: {printer._print(expr)},"
f" {printer._print(a)}, {printer._print(b)},"
f" epsabs={self.abs_tolerance}, epsrel={self.abs_tolerance},"
f" limit={self.limit})[0]"
)
To test whether this works, test this integral expression on another {func}~ampform.sympy.unevaluated
expression:
@unevaluated
class MyFunction(sp.Expr):
x: sp.Symbol
omega1: sp.Symbol
omega2: sp.Symbol
phi1: sp.Symbol
phi2: sp.Symbol
_latex_repr_ = R"f\left({x}\right)"
def evaluate(self) -> sp.Expr:
x, omega1, omega2, phi1, phi2 = self.args
return sp.sin(omega1 * x + phi1) + sp.sin(omega2 * x + phi2)
x, omega1, omega2, phi1, phi2 = sp.symbols("x omega1 omega2 phi1 phi2")
expr = MyFunction(x, omega1, omega2, phi1, phi2)
Math(aslatex({expr: expr.doit(deep=False)}))
w, a, b = sp.symbols("w a b")
fourier_expr = UnevaluatableIntegral(expr * sp.exp(-sp.I * w * x), (x, a, b))
fourier_expr
Indeed the expression correctly lambdifies correctly, despite the {meth}~sympy.core.basic.Basic.doit
call:
func = sp.lambdify([x, omega1, omega2, phi1, phi2], expr.doit())
fourier_func = sp.lambdify([w, omega1, omega2, phi1, phi2, a, b], fourier_expr.doit())
src = inspect.getsource(fourier_func)
src = f"""```python
{black.format_str(src, mode=black.FileMode()).strip()}
```"""
Markdown(src)
domain = np.linspace(-7, +7, num=500)
parameters = dict(
omega1=1.2,
omega2=2.3,
phi1=-1.2,
phi2=+0.4,
)
func_output = func(domain, **parameters)
fourier_output = fourier_func(domain, **parameters, a=-10, b=+10)
%config InlineBackend.figure_formats = ['svg']
fig, ax = plt.subplots()
ax.set_xlabel("$x,w$")
ax.plot(domain, func_output, label="$f(x)$")
ax.plot(domain, fourier_output.real, label=R"$\mathrm{Re}\,F(w)$")
ax.plot(domain, fourier_output.imag, label=R"$\mathrm{Im}\,F(w)$")
ax.legend()
plt.show()
:::{tip} See how this integral expression class is applied to the phase space factor in TR-003. :::