From 79d8ac4ca7daf2a5ca133e6f5eabf87ba7c649ab Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Sun, 21 Jan 2024 12:17:02 +0100 Subject: [PATCH 1/5] MAINT: unfreeze TR-003 --- docs/conf.py | 1 - docs/report/003.ipynb | 185 ++++-------------------------------------- 2 files changed, 18 insertions(+), 168 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 5252c9d3..fe2c6890 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -29,7 +29,6 @@ def get_nb_exclusion_patterns() -> list[str]: "report/000*", "report/001*", "report/002*", - "report/003*", "report/005*", "report/008*", "report/009*", diff --git a/docs/report/003.ipynb b/docs/report/003.ipynb index a6f6df96..43ed553a 100644 --- a/docs/report/003.ipynb +++ b/docs/report/003.ipynb @@ -179,25 +179,10 @@ }, "tags": [ "hide-input", - "remove-input", - "keep_output" + "remove-input" ] }, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle q_a = \\sqrt[\\mathrm{c}]{q_a^{2}} = \\begin{cases} i \\sqrt{- q_a^{2}} & \\text{for}\\: q_a^{2} < 0 \\\\\\sqrt{q_a^{2}} & \\text{otherwise} \\end{cases}$" - ], - "text/plain": [ - "" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "q_squared_symbol = sp.Symbol(\"q_a^{2}\", real=True)\n", "q_a_expr = ComplexSqrt(q_squared_symbol)\n", @@ -245,25 +230,10 @@ "execution_count": null, "metadata": { "tags": [ - "full-width", - "keep_output" + "full-width" ] }, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle \\frac{1}{16 \\pi^{2}} \\left(\\frac{2 \\sqrt[\\mathrm{c}]{q^2\\left(s\\right)}}{\\sqrt{s}} \\log{\\left(\\frac{m_{1}^{2} + m_{2}^{2} + 2 \\sqrt{s} \\sqrt[\\mathrm{c}]{q^2\\left(s\\right)} - s}{2 m_{1} m_{2}} \\right)} - \\left(m_{1}^{2} - m_{2}^{2}\\right) \\left(- \\frac{1}{\\left(m_{1} + m_{2}\\right)^{2}} + \\frac{1}{s}\\right) \\log{\\left(\\frac{m_{1}}{m_{2}} \\right)}\\right)$" - ], - "text/plain": [ - "(1/(16*pi**2))*((2*ComplexSqrt(BreakupMomentumSquared(s, m1, m2))/sqrt(s))*log((m1**2 + m2**2 + 2*sqrt(s)*ComplexSqrt(BreakupMomentumSquared(s, m1, m2)) - s)/(2*m1*m2)) - (m1**2 - m2**2)*(-1/(m1 + m2)**2 + 1/s)*log(m1/m2))" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "s, m1, m2 = sp.symbols(\"s m1 m2\", real=True)\n", "chew_mandelstam_s_wave_expr = chew_mandelstam_s_wave(s, m1, m2)\n", @@ -296,25 +266,10 @@ "source_hidden": true }, "tags": [ - "remove-input", - "keep_output" + "remove-input" ] }, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle s \\to s^{\\prime} + 10^{- \\epsilon} i$" - ], - "text/plain": [ - "" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "Math(Rf\"{sp.latex(s)} \\to {sp.latex(s_plus)}\")" ] @@ -563,26 +518,8 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [ - "keep_output" - ] - }, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle \\frac{B_{L}^2\\left(\\frac{q^2\\left(s^{\\prime}\\right)}{q_{0}^{2}}\\right) \\rho\\left(s^{\\prime}\\right)}{\\left(s^{\\prime} - \\left(m_{1} + m_{2}\\right)^{2}\\right) \\left(- i \\epsilon - s + s^{\\prime}\\right)}$" - ], - "text/plain": [ - "BlattWeisskopfSquared(L, BreakupMomentumSquared(s^{\\prime}, m1, m2)/q0**2)*PhaseSpaceFactor(s^{\\prime}, m1, m2)/((s^{\\prime} - (m1 + m2)**2)*(-I*epsilon - s + s^{\\prime}))" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "metadata": {}, + "outputs": [], "source": [ "q0 = sp.Symbol(\"q0\", real=True)\n", "L = sp.Symbol(\"L\", integer=True, positive=True)\n", @@ -648,20 +585,8 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [ - "keep_output" - ] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Computing for L ∈ [1, 2, 3]\n" - ] - } - ], + "metadata": {}, + "outputs": [], "source": [ "s_domain = np.linspace(s_min, s_max, num=50)\n", "max_L = 3\n", @@ -687,21 +612,8 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [ - "keep_output" - ] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU times: user 2.61 s, sys: 7.39 ms, total: 2.62 s\n", - "Wall time: 2.62 s\n" - ] - } - ], + "metadata": {}, + "outputs": [], "source": [ "%%time\n", "s_thr_val = float(s_thr.subs({m1: m1_val, m2: m2_val}))\n", @@ -863,26 +775,8 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [ - "keep_output" - ] - }, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle \\frac{s - \\left(m_{1} + m_{2}\\right)^{2}}{\\pi} \\int\\limits_{\\left(m_{1} + m_{2}\\right)^{2}}^{\\infty} \\frac{B_{L}^2\\left(q^2\\left(x\\right)\\right) \\rho\\left(x\\right)}{\\left(x - \\left(m_{1} + m_{2}\\right)^{2}\\right) \\left(- i \\epsilon - s + x\\right)}\\, dx$" - ], - "text/plain": [ - "((s - (m1 + m2)**2)/pi)*Integral(BlattWeisskopfSquared(L, BreakupMomentumSquared(x, m1, m2))*PhaseSpaceFactor(x, m1, m2)/((x - (m1 + m2)**2)*(-I*epsilon - s + x)), (x, (m1 + m2)**2, oo))" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "metadata": {}, + "outputs": [], "source": [ "def dispersion_integral(\n", " s,\n", @@ -963,38 +857,10 @@ "source_hidden": true }, "tags": [ - "hide-input", - "keep_output" + "hide-input" ] }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "def _lambdifygenerated(s, m1, m2, epsilon):\n", - " x0 = pi ** (-1.0)\n", - " x1 = (m1 + m2) ** 2\n", - " x2 = -x1\n", - " return (\n", - " x0\n", - " * (s + x2)\n", - " * quadpy_quad(\n", - " lambda x: (1 / 16)\n", - " * x0\n", - " * sqrt((x + x2) * (x - (m1 - m2) ** 2) / x)\n", - " / (sqrt(x) * (x + x2) * (-1j * epsilon - s + x)),\n", - " x1,\n", - " PINF,\n", - " epsabs=0.0001,\n", - " epsrel=0.0001,\n", - " limit=50,\n", - " )[0]\n", - " )\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "src = inspect.getsource(integral_func_s_wave.pyfunc)\n", "src = black.format_str(src, mode=black.FileMode())\n", @@ -1012,23 +878,8 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "tags": [ - "keep_output" - ] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU times: user 5.13 s, sys: 0 ns, total: 5.13 s\n", - "Wall time: 5.13 s\n", - "CPU times: user 2.41 s, sys: 0 ns, total: 2.41 s\n", - "Wall time: 2.41 s\n" - ] - } - ], + "metadata": {}, + "outputs": [], "source": [ "s_values = np.linspace(-0.15, 1.4, num=200)\n", "%time s_wave_values = integral_func_s_wave(s_values, m1_val, m2_val, epsilon=1e-5)\n", From 32046e141f7cf64ae949b59d738e3df2e98f85c5 Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Sun, 21 Jan 2024 13:36:32 +0100 Subject: [PATCH 2/5] ENH: compute integral with `scipy.integrate.quad_vec()` --- .cspell.json | 3 + docs/report/003.ipynb | 758 ++++++++++++++++++++++++++++-------------- pyproject.toml | 1 + 3 files changed, 515 insertions(+), 247 deletions(-) diff --git a/.cspell.json b/.cspell.json index db9c1e14..4e5ecba1 100644 --- a/.cspell.json +++ b/.cspell.json @@ -169,6 +169,8 @@ "einsum", "elif", "endswith", + "epsabs", + "epsrel", "eqnarray", "evaluatable", "expertsystem", @@ -269,6 +271,7 @@ "pvalues", "py's", "pycode", + "pyfunc", "pygments", "pylance", "pypi", diff --git a/docs/report/003.ipynb b/docs/report/003.ipynb index 43ed553a..5a26a66f 100644 --- a/docs/report/003.ipynb +++ b/docs/report/003.ipynb @@ -42,12 +42,12 @@ }, "source": [ "::::{margin}\n", - ":::{card} Chew-Mandelstam S-wave and dispersion integrals\n", + ":::{card} Chew-Mandelstam dispersion integrals\n", "TR-003\n", "^^^\n", - "Section {ref}`report/003:S-wave` has been implemented in [ampform#265](https://github.com/ComPWA/ampform/issues/265).\n", + "This report formulates a symbolic dispersion integral to approach the left-hand cut in the form factor for arbitrary angular momentum. The integral is evaluated with SciPy's [`quad_vec`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad_vec.html).\n", "+++\n", - "_To be implemented_\n", + "🚧 [ampform#265](https://github.com/ComPWA/ampform/issues/265)\n", ":::\n", "::::" ] @@ -59,54 +59,47 @@ "# Chew-Mandelstam" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ - "%pip install -q ampform==0.13.3 black==22.3.0 matplotlib==3.5.1 mpl-interactions==0.20.2 ndim==0.1.7 numpy==1.22.3 qrules==0.9.7 quadpy==0.16.15 scipy==1.8.0 sympy==1.10.1 x21==0.2.16" + "%pip install -q ampform==0.14.8 black==23.12.1 mpl-interactions==0.24.1 scipy==1.12.0" ] }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "This report is an attempt formulate the Chew-Mandelstam function described in {pdg-review}`2021; Resonances; pp.13` (Section 50.3.5) with [SymPy](https://docs.sympy.org), so that it can be implemented in [AmpForm](https://ampform.rtfd.io)." - ] - }, - { - "cell_type": "code", - "execution_count": null, "metadata": { - "tags": [ - "remove-input" - ] + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] }, - "outputs": [], "source": [ - "%matplotlib widget" + "This report is an attempt formulate the Chew-Mandelstam function described in [PDG2023, §50.3.3 (Resonances), pp.14–15](https://pdg.lbl.gov/2023/reviews/rpp2023-rev-resonances.pdf#page=15) with [SymPy](https://docs.sympy.org), so that it can be implemented in [AmpForm](https://ampform.rtfd.io)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ "hide-cell" ] @@ -114,26 +107,31 @@ "outputs": [], "source": [ "%config InlineBackend.figure_formats = ['svg']\n", + "%matplotlib widget\n", "import inspect\n", "import warnings\n", "from functools import partial\n", + "from typing import Any\n", "\n", "import black\n", "import matplotlib.pyplot as plt\n", "import mpl_interactions.ipyplot as iplt\n", "import numpy as np\n", "import qrules\n", - "import quadpy\n", - "import symplot\n", "import sympy as sp\n", "from ampform.dynamics import (\n", " BlattWeisskopfSquared,\n", " BreakupMomentumSquared,\n", - " ComplexSqrt,\n", " PhaseSpaceFactor,\n", " PhaseSpaceFactorComplex,\n", ")\n", - "from IPython.display import Math\n", + "from ampform.io import aslatex\n", + "from ampform.sympy import unevaluated\n", + "from ampform.sympy.math import ComplexSqrt\n", + "from IPython.display import SVG, Math, display\n", + "from ipywidgets import FloatSlider\n", + "from scipy.integrate import quad_vec\n", + "from tqdm.auto import tqdm\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "PDG = qrules.load_pdg()" @@ -142,6 +140,10 @@ { "cell_type": "markdown", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "source": [ @@ -150,9 +152,15 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "As can be seen in Eq. (50.40) on {pdg-review}`2021; Resonances; p.13`, the Chew-Mandelstam function $\\Sigma_a$ for a particle $a$ decaying to particles $1, 2$ has a simple form for angular momentum $L=0$ ($S$-wave):\n", + "As can be seen in Eq. (50.44) on [PDG2023, §Resonances, p.15](https://pdg.lbl.gov/2023/reviews/rpp2023-rev-resonances.pdf#page=15), the Chew-Mandelstam function $\\Sigma_a$ for a particle $a$ decaying to particles $1, 2$ has a simple form for angular momentum $L=0$ ($S$-wave):\n", "\n", ":::{math}\n", ":class: full-width\n", @@ -174,136 +182,181 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "hide-input", - "remove-input" + "hide-input" ] }, "outputs": [], "source": [ - "q_squared_symbol = sp.Symbol(\"q_a^{2}\", real=True)\n", - "q_a_expr = ComplexSqrt(q_squared_symbol)\n", - "Math(f\"q_a = {sp.latex(q_a_expr)} = {sp.latex(q_a_expr.evaluate())}\")" + "z = sp.Symbol(\"z\")\n", + "sqrt_expr = ComplexSqrt(z)\n", + "Math(aslatex({sqrt_expr: sqrt_expr.get_definition()}))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { - "tags": [] + "editable": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input", + "scroll-input", + "full-width" + ] }, "outputs": [], "source": [ - "def breakup_momentum(s, m1, m2):\n", - " q_squared = BreakupMomentumSquared(s, m1, m2)\n", - " return ComplexSqrt(q_squared)\n", + "@unevaluated\n", + "class BreakupMomentum(sp.Expr):\n", + " s: Any\n", + " m1: Any\n", + " m2: Any\n", + " _latex_repr_ = R\"q_a\\left({s}\\right)\"\n", "\n", + " def evaluate(self) -> sp.Expr:\n", + " s, m1, m2 = self.args\n", + " q_squared = BreakupMomentumSquared(s, m1, m2)\n", + " return ComplexSqrt(q_squared)\n", "\n", - "def chew_mandelstam_s_wave(s, m1, m2):\n", - " # evaluate=False in order to keep same style as PDG\n", - " q = breakup_momentum(s, m1, m2)\n", - " left_term = sp.Mul(\n", - " 2 * q / sp.sqrt(s),\n", - " sp.log((m1**2 + m2**2 - s + 2 * sp.sqrt(s) * q) / (2 * m1 * m2)),\n", - " evaluate=False,\n", - " )\n", - " right_term = (m1**2 - m2**2) * (1 / s - 1 / (m1 + m2) ** 2) * sp.log(m1 / m2)\n", - " return sp.Mul(\n", - " 1 / (16 * sp.pi**2),\n", - " left_term - right_term,\n", - " evaluate=False,\n", - " )" + "\n", + "@unevaluated\n", + "class ChewMandelstamSWave(sp.Expr):\n", + " s: Any\n", + " m1: Any\n", + " m2: Any\n", + " _latex_repr_ = R\"\\Sigma\\left({s}\\right)\"\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " # evaluate=False in order to keep same style as PDG\n", + " s, m1, m2 = self.args\n", + " q = BreakupMomentum(s, m1, m2)\n", + " left_term = sp.Mul(\n", + " 2 * q / sp.sqrt(s),\n", + " sp.log((m1**2 + m2**2 - s + 2 * sp.sqrt(s) * q) / (2 * m1 * m2)),\n", + " evaluate=False,\n", + " )\n", + " right_term = (m1**2 - m2**2) * (1 / s - 1 / (m1 + m2) ** 2) * sp.log(m1 / m2)\n", + " return sp.Mul(\n", + " 1 / (16 * sp.pi**2),\n", + " left_term - right_term,\n", + " evaluate=False,\n", + " )\n", + "\n", + "\n", + "s, m1, m2 = sp.symbols(\"s m1 m2 \")\n", + "cm_expr = ChewMandelstamSWave(s, m1, m2)\n", + "q_expr = BreakupMomentum(s, m1, m2)\n", + "Math(aslatex({e: e.doit(deep=False) for e in [cm_expr, q_expr]}))" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "To check whether this implementation is correct, let's plug some {class}`~sympy.core.symbol.Symbol`s into this function and compare it to Eq. (50.40) on {pdg-review}`2021; Resonances; p.13`:" + "It should be noted that this equation is not well-defined along the real axis, that is, for $\\mathrm{Im}(s) = 0$. For this reason, we split $s$ into a real part $s'$ with a small imaginary offset (the PDG indicates this with $s+0i$). We parametrized this imaginary offset with $\\epsilon$, and for the interactive plot, we do so with a power of $10$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [ - "full-width" + "hide-input" ] }, "outputs": [], "source": [ - "s, m1, m2 = sp.symbols(\"s m1 m2\", real=True)\n", - "chew_mandelstam_s_wave_expr = chew_mandelstam_s_wave(s, m1, m2)\n", - "chew_mandelstam_s_wave_expr" + "s_prime = sp.Symbol(R\"s^{\\prime}\", real=True)\n", + "epsilon = sp.Symbol(\"epsilon\", positive=True)\n", + "s_plus = s_prime + sp.I * sp.Pow(10, -epsilon)\n", + "Math(Rf\"{sp.latex(s)} \\to {sp.latex(s_plus)}\")" ] }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "It should be noted that this equation is not well-defined along the real axis, that is, for $\\mathrm{Im}(s) = 0$. For this reason, we split $s$ into a real part $s'$ with a small imaginary offset (the PDG indicates this with $s+0i$). We parametrized this imaginary offset with $\\epsilon$, and for the interactive plot, we do so with a power of $10$:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, "source": [ - "epsilon = sp.Symbol(\"epsilon\", positive=True)\n", - "s_prime = sp.Symbol(R\"s^{\\prime}\", real=True)\n", - "s_plus = s_prime + sp.I * sp.Pow(10, -epsilon)" + "We are now ready to use express the symbolic expressions above as a numerical function. For comparison, we will plot the Chew-Mandelstam function for $S$-waves next to AmpForm's {class}`~ampform.dynamics.phasespace.PhaseSpaceFactorComplex`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "remove-input" + "hide-input" ] }, "outputs": [], "source": [ - "Math(Rf\"{sp.latex(s)} \\to {sp.latex(s_plus)}\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We are now ready to use [`mpl_interactions`](https://mpl-interactions.rtfd.io) and AmpForm's {mod}`symplot` to visualize this function:" + "rho_expr = PhaseSpaceFactorComplex(s, m1, m2)\n", + "Math(aslatex({rho_expr: rho_expr.doit(deep=False)}))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "outputs": [], "source": [ - "chew_mandelstam_s_wave_prime = chew_mandelstam_s_wave_expr.subs(s, s_plus)\n", - "np_chew_mandelstam_s_wave, sliders = symplot.prepare_sliders(\n", - " expression=chew_mandelstam_s_wave_prime.doit(),\n", - " plot_symbol=s_prime,\n", - ")\n", - "np_phase_space_factor = sp.lambdify(\n", - " args=(s_prime, m1, m2, epsilon),\n", - " expr=PhaseSpaceFactorComplex(s_plus, m1, m2).doit(),\n", - " modules=\"numpy\",\n", - ")" + "symbols = (s_prime, m1, m2, epsilon)\n", + "cm_func = sp.lambdify(symbols, cm_expr.doit().subs(s, s_plus))\n", + "rho_func = sp.lambdify(symbols, rho_expr.doit().subs(s, s_plus))" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "As starting values for the interactive plot, we assume $\\pi\\eta$ scattering (just like in the PDG section) and use their masses as values for $m_1$ and $m_1$, respectively." ] @@ -311,44 +364,65 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Define sliders and plot range" + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-cell" + ] + }, "outputs": [], "source": [ "s_min, s_max = -0.15, 1.4\n", + "plot_domain = np.linspace(s_min, s_max, num=500)\n", + "\n", "m1_val = PDG[\"pi0\"].mass\n", "m2_val = PDG[\"eta\"].mass\n", - "\n", - "plot_domain = np.linspace(s_min, s_max, 500)\n", - "sliders.set_ranges(\n", - " m1=(0, 2, 200),\n", - " m2=(0, 2, 200),\n", - " epsilon=(1, 12),\n", - ")\n", - "sliders.set_values(\n", - " m1=m1_val,\n", - " m2=m2_val,\n", - " epsilon=4,\n", + "sliders = dict(\n", + " m1=FloatSlider(\n", + " description=str(m1),\n", + " min=0,\n", + " max=1.2,\n", + " value=m1_val,\n", + " ),\n", + " m2=FloatSlider(\n", + " description=str(m2),\n", + " min=0,\n", + " max=1.2,\n", + " value=m2_val,\n", + " ),\n", + " epsilon=FloatSlider(\n", + " description=str(epsilon),\n", + " min=0,\n", + " max=8,\n", + " value=4,\n", + " ),\n", ")" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For comparison, we plot the Chew-Mandelstam function for $S$-waves next to AmpForm's {class}`~ampform.dynamics.phasespace.PhaseSpaceFactorComplex`. Have a look at the resulting plots and compare to Figure 50.4 on {pdg-review}`2021; Resonances; p.12`." - ] - }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "full-width", "remove-output", - "hide-cell" + "scroll-input", + "hide-input" ] }, "outputs": [], @@ -362,6 +436,7 @@ "imag_style = {\"label\": \"Imag part\", \"c\": \"red\"}\n", "threshold_style = {\"label\": R\"$s_\\mathrm{thr}$\", \"c\": \"grey\", \"linewidth\": 0.5}\n", "\n", + "xlim = (plot_domain.min(), plot_domain.max())\n", "ylim = (-1, +1)\n", "y_factor = 16 * np.pi\n", "controls = iplt.axvline(\n", @@ -378,22 +453,20 @@ ")\n", "iplt.plot(\n", " plot_domain,\n", - " lambda *args, **kwargs: (\n", - " y_factor * 1j * np_phase_space_factor(*args, **kwargs)\n", - " ).real,\n", + " lambda *args, **kwargs: (y_factor * 1j * rho_func(*args, **kwargs)).real,\n", " controls=controls,\n", " ylim=ylim,\n", + " xlim=xlim,\n", " alpha=0.7,\n", " ax=ax1,\n", " **real_style,\n", ")\n", "iplt.plot(\n", " plot_domain,\n", - " lambda *args, **kwargs: (\n", - " y_factor * 1j * np_phase_space_factor(*args, **kwargs)\n", - " ).imag,\n", + " lambda *args, **kwargs: (y_factor * 1j * rho_func(*args, **kwargs)).imag,\n", " controls=controls,\n", " ylim=ylim,\n", + " xlim=xlim,\n", " alpha=0.7,\n", " ax=ax1,\n", " **imag_style,\n", @@ -401,20 +474,20 @@ "\n", "iplt.plot(\n", " plot_domain,\n", - " lambda *args, **kwargs: y_factor\n", - " * np_chew_mandelstam_s_wave(*args, **kwargs).real,\n", + " lambda *args, **kwargs: y_factor * cm_func(*args, **kwargs).real,\n", " controls=controls,\n", " ylim=ylim,\n", + " xlim=xlim,\n", " alpha=0.7,\n", " ax=ax2,\n", " **real_style,\n", ")\n", "iplt.plot(\n", " plot_domain,\n", - " lambda *args, **kwargs: y_factor\n", - " * np_chew_mandelstam_s_wave(*args, **kwargs).imag,\n", + " lambda *args, **kwargs: y_factor * cm_func(*args, **kwargs).imag,\n", " controls=controls,\n", " ylim=ylim,\n", + " xlim=xlim,\n", " alpha=0.7,\n", " ax=ax2,\n", " **imag_style,\n", @@ -424,6 +497,7 @@ " ax.legend(loc=\"lower right\")\n", " ax.set_xticks(np.arange(0, 1.21, 0.3))\n", " ax.set_yticks(np.arange(-1, 1.1, 0.5))\n", + " ax.set_xlim()\n", " ax.set_xlabel(\"$s$ (GeV$^2$)\")\n", "\n", "ax1.set_ylabel(R\"$16\\pi \\; i\\rho(s)$\")\n", @@ -437,30 +511,49 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "remove-cell" + "remove-input", + "full-width" ] }, "outputs": [], "source": [ - "plt.savefig(\"003-chew-mandelstam-s-wave.svg\")" + "if STATIC_WEB_PAGE:\n", + " filename = \"003/chew-mandelstam-s-wave.svg\"\n", + " os.makedirs(os.path.dirname(filename), exist_ok=True)\n", + " fig.savefig(filename)\n", + " display(*sliders.values(), SVG(filename))" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "```{figure} https://user-images.githubusercontent.com/29308176/164984924-764a9558-6afd-46a9-8f24-8cc92ce1bc49.svg\n", - ":class: full-width\n", - "```" + ":::{tip}\n", + "Compare the plots above with Figure 50.6 on [PDG2023, §Resonances, p.16](https://pdg.lbl.gov/2023/reviews/rpp2023-rev-resonances.pdf#page=15).\n", + ":::" ] }, { "cell_type": "markdown", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "source": [ @@ -469,9 +562,15 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "For higher angular momenta, the PDG notes that one has to compute the dispersion integral given by Eq. (50.41) on {pdg-review}`2021; Resonances; p.13`:\n", + "For higher angular momenta, the PDG notes that one has to compute the dispersion integral given by Eq. (50.44) on [PDG2023, §Resonances, p.15](https://pdg.lbl.gov/2023/reviews/rpp2023-rev-resonances.pdf#page=15):\n", "\n", "$$\n", "\\Sigma_a(s+0i) =\n", @@ -482,35 +581,76 @@ " (s' - s_{\\mathrm{thr}_a})(s'-s-i0)\n", " }\n", " \\mathop{}\\!\\mathrm{d}s'\n", - "$$ (dispersion-integral)" + "$$ (dispersion-integral)\n", + "\n", + "Equation {eq}`chew-mandelstam` is the analytic solution for $L=0$." ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "Equation {eq}`chew-mandelstam` is the analytic solution for $L=0$.\n", - "\n", - "From Equations (50.26-27) on {pdg-review}`2021; Resonances; p.9`, it can be deduced that the function $n_a^2$ is the same as AmpForm's {class}`~ampform.dynamics.BlattWeisskopfSquared` (note that this function is normalized, whereas the PDG's $F_j$ function has $1$ in the nominator). Furthermore, the PDG seems to suggest that $z = q_a/q_0$, but this is an unconventional choice and is probably a mistake. For this reason, we simply use {class}`~ampform.dynamics.BlattWeisskopfSquared` for the definition of $n_a^2$:" + "From Equations (50.33–34) on [PDG2023, §Resonances, p.12](https://pdg.lbl.gov/2023/reviews/rpp2023-rev-resonances.pdf#page=12), it can be deduced that the function $n_a^2$ is the same as AmpForm's {class}`~ampform.dynamics.BlattWeisskopfSquared` (note that this function is normalized, whereas the PDG's $F_j$ function has $1$ in the nominator). For this reason, we simply use {class}`~ampform.dynamics.BlattWeisskopfSquared` for the definition of $n_a^2$:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input", + "full-width" + ] + }, "outputs": [], "source": [ - "def na2(s, m1, m2, L, q0):\n", - " q_squared = BreakupMomentumSquared(s, m1, m2)\n", - " return BlattWeisskopfSquared(\n", - " z=q_squared / (q0**2),\n", - " angular_momentum=L,\n", - " )" + "@unevaluated\n", + "class FormFactor(sp.Expr):\n", + " s: Any\n", + " m1: Any\n", + " m2: Any\n", + " L: Any\n", + " q0: Any = 1\n", + " _latex_repr_ = R\"n_a^2\\left({s}\\right)\"\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " s, m1, m2, L, q0 = self.args\n", + " q_squared = BreakupMomentumSquared(s, m1, m2)\n", + " return BlattWeisskopfSquared(\n", + " z=q_squared / (q0**2),\n", + " angular_momentum=L,\n", + " )\n", + "\n", + "\n", + "L = sp.Symbol(\"L\", integer=True, positive=True)\n", + "q0 = sp.Symbol(\"q0\", real=True)\n", + "na_expr = FormFactor(s, m1, m2, L, q0)\n", + "bl_expr = BlattWeisskopfSquared(z, L)\n", + "Math(aslatex({e: e.doit(deep=False) for e in [na_expr, bl_expr]}))" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "For $\\rho_a$, we use AmpForm's {class}`~ampform.dynamics.phasespace.PhaseSpaceFactor`:" ] @@ -518,124 +658,183 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, "outputs": [], "source": [ - "q0 = sp.Symbol(\"q0\", real=True)\n", - "L = sp.Symbol(\"L\", integer=True, positive=True)\n", - "s_thr = (m1 + m2) ** 2\n", - "integrand = (PhaseSpaceFactor(s_prime, m1, m2) * na2(s_prime, m1, m2, L, q0)) / (\n", - " (s_prime - s_thr) * (s_prime - s - epsilon * sp.I)\n", - ")\n", - "integrand" + "Math(aslatex({rho_expr: rho_expr.doit(deep=False)}))" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "Next, we {func}`~sympy.utilities.lambdify.lambdify` this integrand to a {mod}`numpy` expression so that we can integrate it efficiently:" + "The symbolic integrand is then formulated as:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, "outputs": [], "source": [ - "np_integrand = sp.lambdify(\n", - " args=(s_prime, s, L, epsilon, m1, m2, q0),\n", - " expr=integrand.doit(),\n", - " modules=\"numpy\",\n", - ")" + "s_thr = (m1 + m2) ** 2\n", + "integrand = (\n", + " PhaseSpaceFactor(s_prime, m1, m2) * FormFactor(s_prime, m1, m2, L, q0)\n", + ") / ((s_prime - s_thr) * (s_prime - s - epsilon * sp.I))\n", + "integrand" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "As discussed in [TR-016](016.ipynb), {func}`scipy.integrate.quad` cannot integrate over complex-valued functions, while [`quadpy`](https://github.com/sigma-py/quadpy) runs into trouble with vectorized input to the integrand. The following function, from {ref}`report/016:Vectorized input` offers a quick solution:" + "Next, we {func}`~sympy.utilities.lambdify.lambdify` this integrand to a {mod}`numpy` expression so that we can integrate it efficiently:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ - "@np.vectorize\n", - "def vectorized_quad(func, a, b, **func_kwargs):\n", - " values, _ = quadpy.quad(partial(func, **func_kwargs), a, b)\n", - " return values" + "integrand_func = sp.lambdify(\n", + " args=(s_prime, s, L, epsilon, m1, m2, q0),\n", + " expr=integrand.doit(),\n", + " modules=\"numpy\",\n", + ")" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ ":::{note}\n", - "\n", - "Integrals can be expressed with SymPy, with some caveats. See {ref}`report/016:SymPy integral`.\n", - "\n", - ":::\n", - "\n", - "Now, for comparison, we compute this integral for a few values of $L>0$:" + "Integrals can be expressed symbolically with SymPy, with some caveats. See {ref}`report/016:SymPy integral`.\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "As discussed in [TR-016](016.ipynb), {func}`scipy.integrate.quad` cannot integrate over complex-valued functions, but {func}`scipy.integrate.quad_vec` can. For comparison, we now compute this integral for a few values of $L>0$:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ "s_domain = np.linspace(s_min, s_max, num=50)\n", "max_L = 3\n", - "l_values = list(range(1, max_L + 1))\n", - "print(\"Computing for L ∈\", l_values)" + "l_values = list(range(1, max_L + 1))" ] }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "It is handy to store the resulting values of each dispersion integral in a {obj}`dict` with $L$ as keys:" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "```{autolink-skip}\n", - "```" + "It is handy to store the numerical results of each dispersion integral in a {obj}`dict` with $L$ as keys:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, "outputs": [], "source": [ - "%%time\n", "s_thr_val = float(s_thr.subs({m1: m1_val, m2: m2_val}))\n", "integral_values = {\n", - " l_val: vectorized_quad(\n", - " np_integrand,\n", + " l_val: quad_vec(\n", + " lambda x: integrand_func(\n", + " x,\n", + " s=s_domain,\n", + " L=l_val,\n", + " epsilon=1e-3,\n", + " m1=m1_val,\n", + " m2=m2_val,\n", + " q0=1.0,\n", + " ),\n", " a=s_thr_val,\n", " b=np.inf,\n", - " s=s_domain,\n", - " L=l_val,\n", - " epsilon=1e-3,\n", - " m1=m1_val,\n", - " m2=m2_val,\n", - " q0=1.0,\n", - " )\n", - " for l_val in l_values\n", + " )[0]\n", + " for l_val in tqdm(l_values, desc=\"Evaluating integrals\")\n", "}" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "Finally, as can be seen from Eq. {eq}`dispersion-integral`, the resulting values from the integral have to be shifted with a factor $\\frac{s-s_{\\mathrm{thr}_a}}{\\pi}$ to get $\\Sigma_a$. We also scale the values with $16\\pi$ so that it can be compared with the plot generated in {ref}`report/003:S-wave`." ] @@ -644,6 +843,10 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "outputs": [], @@ -659,11 +862,16 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "remove-cell" + "remove-output", + "hide-input" ] }, "outputs": [], @@ -691,34 +899,35 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "remove-cell" + "remove-input" ] }, "outputs": [], "source": [ - "plt.savefig(\"003-chew-mandelstam-l-non-zero.svg\")" + "if STATIC_WEB_PAGE:\n", + " filename = \"003/chew-mandelstam-l-non-zero.svg\"\n", + " os.makedirs(os.path.dirname(filename), exist_ok=True)\n", + " fig.savefig(filename)\n", + " display(SVG(filename))" ] }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "![Chew-Mandelstam for higher angular momenta](https://user-images.githubusercontent.com/29308176/164985017-7600941e-0481-4282-8d9b-0680f720e6ef.svg)\n", - "\n", - ":::{note}\n", - "\n", - "In {ref}`report/003:SymPy expressions` we'll see that the dispersion integral indeed reproduces the same shape as the analytic expression from {ref}`report/003:S-wave`.\n", - "\n", - ":::" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "## SymPy expressions" ] @@ -734,8 +943,12 @@ "cell_type": "code", "execution_count": null, "metadata": { - "jupyter": { - "source_hidden": true + "editable": true, + "mystnb": { + "code_prompt_show": "Symbolic integral definition" + }, + "slideshow": { + "slide_type": "" }, "tags": [ "hide-cell" @@ -760,13 +973,14 @@ " if len(limits) != 1:\n", " msg = f\"Cannot handle {len(limits)}-dimensional integrals\"\n", " raise ValueError(msg)\n", - " integrate = \"quadpy_quad\"\n", - " printer.module_imports[\"quadpy\"].update({f\"quad as {integrate}\"})\n", - " limit_str = \"{}, {}\".format(*tuple(map(printer._print, limits[0])))\n", + " integrate = \"quad_vec\"\n", + " printer.module_imports[\"scipy.integrate\"].add(integrate)\n", + " a, b = map(printer._print, limits[0])\n", " args = \", \".join(map(printer._print, integration_vars))\n", " expr = printer._print(self.args[0])\n", " return (\n", - " f\"{integrate}(lambda {args}: {expr}, {limit_str},\"\n", + " f\"{integrate}(lambda {args}: {expr},\"\n", + " f\" a={a}, b={b},\"\n", " f\" epsabs={self.abs_tolerance}, epsrel={self.abs_tolerance},\"\n", " f\" limit={self.limit})[0]\"\n", " )" @@ -775,7 +989,21 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Definition of the symbolic dispersion integral" + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, "outputs": [], "source": [ "def dispersion_integral(\n", @@ -812,11 +1040,17 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ ":::{warning}\n", "\n", - "We have to keep track of the integration variable ($s'$ in Equation {eq}`dispersion-integral`), so that we don't run into trouble if we use {func}`~sympy.utilities.lambdify.lambdify` with common sub-expressions. The problem is that the integration variable _should not_ be extracted as a common sub-expression, otherwise the lambdified [`quadpy.quad()`](https://github.com/sigma-py/quadpy) expression cannot handle vectorized input.\n", + "We have to keep track of the integration variable ($s'$ in Equation {eq}`dispersion-integral`), so that we don't run into trouble if we use {func}`~sympy.utilities.lambdify.lambdify` with common sub-expressions. The problem is that the integration variable _should not_ be extracted as a common sub-expression, otherwise the lambdified {func}`scipy.integrate.quad_vec` expression cannot handle vectorized input.\n", "\n", ":::\n", "\n", @@ -827,49 +1061,63 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "outputs": [], "source": [ "UnevaluatableIntegral.abs_tolerance = 1e-4\n", "UnevaluatableIntegral.rel_tolerance = 1e-4\n", - "integral_func_s_wave = sp.lambdify(\n", + "integral_s_wave_func = sp.lambdify(\n", " [s, m1, m2, epsilon],\n", " integral_expr.subs(L, 0).doit(),\n", " # integration symbol should not be extracted as common sub-expression!\n", " cse=partial(sp.cse, ignore=[x], list=False),\n", ")\n", - "integral_func_s_wave = np.vectorize(integral_func_s_wave)\n", + "integral_s_wave_func = np.vectorize(integral_s_wave_func)\n", "\n", - "integral_func_p_wave = sp.lambdify(\n", + "integral_p_wave_func = sp.lambdify(\n", " [s, m1, m2, epsilon],\n", " integral_expr.subs(L, 1).doit(),\n", " cse=partial(sp.cse, ignore=[x], list=False),\n", ")\n", - "integral_func_p_wave = np.vectorize(integral_func_p_wave)" + "integral_p_wave_func = np.vectorize(integral_p_wave_func)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ - "src = inspect.getsource(integral_func_s_wave.pyfunc)\n", + "src = inspect.getsource(integral_s_wave_func.pyfunc)\n", "src = black.format_str(src, mode=black.FileMode())\n", "print(src)" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "```{autolink-skip}\n", "```" @@ -878,17 +1126,28 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, "outputs": [], "source": [ "s_values = np.linspace(-0.15, 1.4, num=200)\n", - "%time s_wave_values = integral_func_s_wave(s_values, m1_val, m2_val, epsilon=1e-5)\n", - "%time p_wave_values = integral_func_p_wave(s_values, m1_val, m2_val, epsilon=1e-5)" + "%time s_wave_values = integral_s_wave_func(s_values, m1_val, m2_val, epsilon=1e-5)\n", + "%time p_wave_values = integral_p_wave_func(s_values, m1_val, m2_val, epsilon=1e-5)" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "Note that the dispersion integral for $L=0$ indeed reproduces the same shape as in {ref}`report/003:S-wave`!" ] @@ -897,9 +1156,13 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ "hide-input", "remove-output" @@ -940,23 +1203,24 @@ "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, "jupyter": { "source_hidden": true }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "remove-cell" + "remove-input" ] }, "outputs": [], "source": [ - "plt.savefig(\"003-symbolic-chew-mandelstam.svg\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "![Symbolic Chew-Mandelstam plots](https://user-images.githubusercontent.com/29308176/164984984-dfe73d4c-e604-4d06-b4e1-50be117a57e3.svg)" + "if STATIC_WEB_PAGE:\n", + " filename = \"003/symbolic-chew-mandelstam.svg\"\n", + " os.makedirs(os.path.dirname(filename), exist_ok=True)\n", + " fig.savefig(filename)\n", + " display(SVG(filename))" ] } ], @@ -979,7 +1243,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.18" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index d50daa6e..4be7a194 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -238,6 +238,7 @@ target-version = "py37" task-tags = ["cspell"] [tool.ruff.per-file-ignores] +"**/003.ipynb" = ["B023"] "**/024.ipynb" = ["E731", "E741", "PGH001", "S307"] "**/025.ipynb" = ["E731"] "*.ipynb" = [ From 4fbf679d015469fc4955ac45c4f7a0b3c30039db Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Sun, 21 Jan 2024 11:19:45 +0100 Subject: [PATCH 3/5] MAINT: unfreeze TR-016 --- docs/conf.py | 10 +++- docs/report/016.ipynb | 135 ++++-------------------------------------- 2 files changed, 19 insertions(+), 126 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index fe2c6890..79f95bad 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -32,7 +32,15 @@ def get_nb_exclusion_patterns() -> list[str]: "report/005*", "report/008*", "report/009*", - "report/01*", + "report/010*", + "report/011*", + "report/012*", + "report/013*", + "report/014*", + "report/015*", + "report/017*", + "report/018*", + "report/019*", "report/020*", "report/021*", "report/022*", diff --git a/docs/report/016.ipynb b/docs/report/016.ipynb index 59793e7d..b69eb314 100644 --- a/docs/report/016.ipynb +++ b/docs/report/016.ipynb @@ -75,15 +75,7 @@ "remove-cell" ] }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], + "outputs": [], "source": [ "%pip install -q ndim==0.1.7 quadpy==0.16.15 scipy==1.8.0 sympy==1.10.1 x21==0.2.16" ] @@ -110,21 +102,7 @@ "raises-exception" ] }, - "outputs": [ - { - "ename": "TypeError", - "evalue": "can't convert complex to float", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", - "Input \u001b[0;32mIn [3]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mintegrand\u001b[39m(x):\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m x \u001b[38;5;241m*\u001b[39m (x \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39mj)\n\u001b[0;32m----> 8\u001b[0m \u001b[43mquad\u001b[49m\u001b[43m(\u001b[49m\u001b[43mintegrand\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m0.0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m2.0\u001b[39;49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/miniconda3/envs/compwa-org/lib/python3.8/site-packages/scipy/integrate/_quadpack_py.py:351\u001b[0m, in \u001b[0;36mquad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)\u001b[0m\n\u001b[1;32m 348\u001b[0m flip, a, b \u001b[38;5;241m=\u001b[39m b \u001b[38;5;241m<\u001b[39m a, \u001b[38;5;28mmin\u001b[39m(a, b), \u001b[38;5;28mmax\u001b[39m(a, b)\n\u001b[1;32m 350\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m weight \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 351\u001b[0m retval \u001b[38;5;241m=\u001b[39m \u001b[43m_quad\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfull_output\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mepsabs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mepsrel\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlimit\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 352\u001b[0m \u001b[43m \u001b[49m\u001b[43mpoints\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 353\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 354\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m points \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n", - "File \u001b[0;32m~/miniconda3/envs/compwa-org/lib/python3.8/site-packages/scipy/integrate/_quadpack_py.py:463\u001b[0m, in \u001b[0;36m_quad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points)\u001b[0m\n\u001b[1;32m 461\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m points \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m 462\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m infbounds \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m--> 463\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_quadpack\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_qagse\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43mfull_output\u001b[49m\u001b[43m,\u001b[49m\u001b[43mepsabs\u001b[49m\u001b[43m,\u001b[49m\u001b[43mepsrel\u001b[49m\u001b[43m,\u001b[49m\u001b[43mlimit\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 464\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 465\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m _quadpack\u001b[38;5;241m.\u001b[39m_qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)\n", - "\u001b[0;31mTypeError\u001b[0m: can't convert complex to float" - ] - } - ], + "outputs": [], "source": [ "from scipy.integrate import quad\n", "\n", @@ -178,18 +156,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "((2.666666666666667+2j), (8.765121169122355e-28+2.220446049250313e-14j))" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "complex_integrate(integrand, 0.0, 2.0)" ] @@ -223,18 +190,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "((2.6666666666666665+2.0000000000000004j), 2.0082667671941473e-19)" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "import quadpy\n", "\n", @@ -280,24 +236,7 @@ "raises-exception" ] }, - "outputs": [ - { - "ename": "ValueError", - "evalue": "operands could not be broadcast together with shapes (21,) (10,) ", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", - "Input \u001b[0;32mIn [7]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m s_prime \u001b[38;5;241m*\u001b[39m (s_prime \u001b[38;5;241m+\u001b[39m s \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39mj)\n\u001b[1;32m 8\u001b[0m s_array \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mlinspace(\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m1\u001b[39m, num\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m10\u001b[39m)\n\u001b[0;32m----> 9\u001b[0m \u001b[43mquadpy\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mquad\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 10\u001b[0m \u001b[43m \u001b[49m\u001b[43mpartial\u001b[49m\u001b[43m(\u001b[49m\u001b[43mparametrized_func\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43ms\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43ms_array\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 11\u001b[0m \u001b[43m \u001b[49m\u001b[43ma\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0.0\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 12\u001b[0m \u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m2.0\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 13\u001b[0m \u001b[43m)\u001b[49m\n", - "File \u001b[0;32m:55\u001b[0m, in \u001b[0;36mquad\u001b[0;34m(f, a, b, args, epsabs, epsrel, limit)\u001b[0m\n", - "File \u001b[0;32m:45\u001b[0m, in \u001b[0;36mintegrate_adaptive\u001b[0;34m(f, intervals, eps_abs, eps_rel, criteria_connection, kronrod_degree, minimum_interval_length, max_num_subintervals, dot, domain_shape, range_shape)\u001b[0m\n", - "File \u001b[0;32m:129\u001b[0m, in \u001b[0;36m_gauss_kronrod_integrate\u001b[0;34m(k, f, intervals, dot, domain_shape, range_shape)\u001b[0m\n", - "File \u001b[0;32m:49\u001b[0m, in \u001b[0;36mg\u001b[0;34m(x)\u001b[0m\n", - "Input \u001b[0;32mIn [7]\u001b[0m, in \u001b[0;36mparametrized_func\u001b[0;34m(s_prime, s)\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mparametrized_func\u001b[39m(s_prime, s):\n\u001b[0;32m----> 5\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m s_prime \u001b[38;5;241m*\u001b[39m (\u001b[43ms_prime\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43ms\u001b[49m \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39mj)\n", - "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (21,) (10,) " - ] - } - ], + "outputs": [], "source": [ "from functools import partial\n", "\n", @@ -325,23 +264,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(array([0.66666667+2.j, 1.11111111+2.j, 1.55555556+2.j, 2. +2.j,\n", - " 2.44444444+2.j, 2.88888889+2.j, 3.33333333+2.j, 3.77777778+2.j,\n", - " 4.22222222+2.j, 4.66666667+2.j]),\n", - " array([1.94765926e-19, 2.69476631e-19, 2.24127752e-19, 2.79100064e-19,\n", - " 2.67216263e-19, 1.43065895e-19, 3.08910645e-19, 3.62329394e-19,\n", - " 4.86795288e-19, 2.21702097e-19]))" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "from functools import partial\n", "\n", @@ -426,21 +349,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle \\int\\limits_{a}^{b} s_{p} \\left(s + s_{p} + i\\right)\\, ds_{p}$" - ], - "text/plain": [ - "Integral(s_p*(s + s_p + I), (s_p, a, b))" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "s, s_prime, a, b = sp.symbols(\"s s_p a b\")\n", "integral_expr: sp.Expr = UnevaluatableIntegral(\n", @@ -461,17 +370,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "def _lambdifygenerated(s, a, b):\n", - " return quadpy_quad(lambda s_p: s_p*(s + s_p + 1j), a, b, epsabs=1e-05, epsrel=1e-05, limit=50)[0]\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "import inspect\n", "\n", @@ -493,20 +392,7 @@ "metadata": { "tags": [] }, - "outputs": [ - { - "data": { - "text/plain": [ - "array([0.66666667+2.j, 1.11111111+2.j, 1.55555556+2.j, 2. +2.j,\n", - " 2.44444444+2.j, 2.88888889+2.j, 3.33333333+2.j, 3.77777778+2.j,\n", - " 4.22222222+2.j, 4.66666667+2.j])" - ] - }, - "execution_count": null, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "vec_integral_func = np.vectorize(integral_func)\n", "vec_integral_func(s_array, a=0.0, b=2.0)" @@ -528,7 +414,6 @@ "colab": { "toc_visible": true }, - "keep_output": true, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", @@ -544,7 +429,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.18" + "version": "3.10.13" } }, "nbformat": 4, From 408580d16132649785416cf4019da10fb574df90 Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Sun, 21 Jan 2024 18:39:23 +0100 Subject: [PATCH 4/5] DOC: implement `UnevaluatedIntegral` with `quad_vec()` --- .cspell.json | 2 + docs/report/016.ipynb | 553 ++++++++++++++++++++++++++++++++---------- 2 files changed, 421 insertions(+), 134 deletions(-) diff --git a/.cspell.json b/.cspell.json index 4e5ecba1..62376a85 100644 --- a/.cspell.json +++ b/.cspell.json @@ -62,6 +62,8 @@ "conda", "Dalitz", "defaultdict", + "dummifies", + "dummify", "eval", "flatté", "functools", diff --git a/docs/report/016.ipynb b/docs/report/016.ipynb index b69eb314..c0fe31e4 100644 --- a/docs/report/016.ipynb +++ b/docs/report/016.ipynb @@ -35,16 +35,20 @@ { "cell_type": "markdown", "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [ "sympy" ] }, "source": [ "::::{margin}\n", - ":::{card} Complex integral\n", + ":::{card} Symbolic integral\n", "TR-016\n", "^^^\n", - "As noted in [TR-003](./003.ipynb), {func}`scipy.integrate.quad` cannot handle complex integrals. In addition, one can get into trouble with vectorized input ({obj}`numpy.array`s) on a lambdified {class}`sympy.Integral `. This report discusses both problems and proposes some solutions.\n", + "This report investigates how to formulate a symbolic integral that correctly evaluates to\n", "+++\n", "To be implemented\n", ":::\n", @@ -53,60 +57,121 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "# Complex integral" + "# Lambdifying a symbolic integral" ] }, { - "cell_type": "markdown", - "metadata": {}, + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "remove-cell" + ] + }, + "outputs": [], "source": [ - "\n", - "\n", - "" + "%pip install -q ampform==0.14.8 black==23.12.1 scipy==1.12.0 sympy==1.12" ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Import Python libraries" + }, + "slideshow": { + "slide_type": "" + }, "tags": [ - "remove-cell" + "hide-input" ] }, "outputs": [], "source": [ - "%pip install -q ndim==0.1.7 quadpy==0.16.15 scipy==1.8.0 sympy==1.10.1 x21==0.2.16" + "import inspect\n", + "\n", + "import black\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import sympy as sp\n", + "from ampform.io import aslatex\n", + "from ampform.sympy import unevaluated\n", + "from IPython.display import Markdown, Math\n", + "from scipy.integrate import quad, quad_vec\n", + "from sympy.printing.pycode import _unpack_integral_limits" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "## Complex integration" + "## Numerical integration" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "SciPy cannot integrate complex functions:" + "(quad)=\n", + "### SciPy's `quad()` function" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "SciPy's {func}`scipy.integrate.quad` cannot integrate complex-valued functions:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [ "raises-exception" ] }, "outputs": [], "source": [ - "from scipy.integrate import quad\n", - "\n", - "\n", "def integrand(x):\n", " return x * (x + 1j)\n", "\n", @@ -116,14 +181,13 @@ }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Split real and imaginary integral" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ "A [proposed solution](https://stackoverflow.com/a/5966088) 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:" ] @@ -131,39 +195,42 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ - "import numpy as np\n", - "\n", - "\n", "def complex_integrate(func, a, b, **quad_kwargs):\n", " def real_func(x):\n", - " return np.real(func(x))\n", + " return func(x).real\n", "\n", " def imag_func(x):\n", - " return np.imag(func(x))\n", + " return func(x).imag\n", "\n", " real_integral, real_integral_err = quad(real_func, a, b, **quad_kwargs)\n", " imag_integral, imag_integral_err = quad(imag_func, a, b, **quad_kwargs)\n", " return (\n", " real_integral + 1j * imag_integral,\n", " real_integral_err**2 + 1j * imag_integral_err,\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + " )\n", + "\n", + "\n", "complex_integrate(integrand, 0.0, 2.0)" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ ":::{warning}\n", "\n", @@ -174,70 +241,121 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "### Integrate with `quadpy`" + "(quad_vec)=\n", + "### SciPy's `quad_vec()` function" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "[Alternatively](https://stackoverflow.com/a/42866568), one could use [`quadpy`](https://github.com/sigma-py/quadpy), which essentially does the same as in {ref}`report/016:Split real and imaginary integral`, but can also (to a large degree) handle vectorized input and properly handles uncertainties." + "The easiest solution, however, seems to be {func}`scipy.integrate.quad_vec`:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ - "import quadpy\n", - "\n", - "quadpy.quad(integrand, a=0.0, b=2.0)" + "quad_vec(integrand, 0.0, 2.0)" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - ":::{note}\n", - "\n", - "One may need to play around with the tolerance if the function is not smooth, see [sigma-py/quadpy#255](https://github.com/sigma-py/quadpy/issues/255).\n", - "\n", - ":::\n", - "\n", - ":::{tip}\n", + "This has the added benefit that it can handle functions that return arrays:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "def gaussian(x, mu, sigma):\n", + " return np.exp(-((x - mu) ** 2) / (2 * sigma**2)) / (sigma * np.sqrt(2 * np.pi))\n", "\n", - "[`quadpy`](https://github.com/sigma-py/quadpy) raises exceptions with {obj}`ModuleNotFoundError`s that are a bit unreadable. They are caused by the fact that [`orthopy`](https://pypi.org/project/orthopy) and [`ndim`](https://pypi.org/project/ndim) need to be installed separately.\n", "\n", - ":::" + "mu_values = np.linspace(-2, +3, num=10)\n", + "result, _ = quad_vec(lambda x: gaussian(x, mu_values, sigma=0.5), 0, 2.0)\n", + "result" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "## Vectorized input" + "### Integrate with `quadpy`" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "The dispersion integral from Eq. {eq}`dispersion-integral` in [TR-003](003.ipynb) features a variable $s$ that is an argument to the function $\\Sigma_a$. This becomes a problem when $s$ gets vectorized (in this case: gets an event-wise {obj}`numpy.array` of invariant masses). Here's a simplified version of the problem:" + ":::{warning}\n", + "`quadpy` now requires a license. The examples below are only shown for documentation purposes.\n", + ":::" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "tags": [ - "raises-exception" - ] + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] }, - "outputs": [], "source": [ + "[Alternatively](https://stackoverflow.com/a/42866568), one could use [`quadpy`](https://github.com/sigma-py/quadpy), which essentially does the same as in [`quad()`](#quad), but can also (to a large degree) handle vectorized input and properly handles uncertainties. For example:\n", + "\n", + "```python\n", "from functools import partial\n", "\n", "\n", @@ -250,72 +368,119 @@ " partial(parametrized_func, s=s_array),\n", " a=0.0,\n", " b=2.0,\n", - ")" + ")\n", + "```" ] }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "The way out seems to be to **vectorize the [`quadpy.quad()`](https://github.com/sigma-py/quadpy) call itself** and forward the function arguments through {func}`functools.partial`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "from functools import partial\n", - "\n", + ":::{note}\n", "\n", - "@np.vectorize\n", - "def vectorized_quad(func, a, b, **func_kwargs):\n", - " return quadpy.quad(partial(func, **func_kwargs), a, b)\n", + "One may need to play around with the tolerance if the function is not smooth, see [sigma-py/quadpy#255](https://github.com/sigma-py/quadpy/issues/255).\n", "\n", + ":::\n", "\n", - "vectorized_quad(parametrized_func, a=0.0, b=2.0, s=s_array)" + ":::{tip}\n", + "\n", + "[`quadpy`](https://github.com/sigma-py/quadpy) raises exceptions with {obj}`ModuleNotFoundError`s that are a bit unreadable. They are caused by the fact that [`orthopy`](https://pypi.org/project/orthopy) and [`ndim`](https://pypi.org/project/ndim) need to be installed separately.\n", + ":::" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "Note, however, that this becomes difficult to implement as {func}`~sympy.utilities.lambdify.lambdify` output for a {class}`sympy.Integral `. An attempt at this is made in [TR-003](003.ipynb)." + "## SymPy integral" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "## SymPy integral" + "The dispersion integral from Eq. {eq}`dispersion-integral` in **[TR-003](003.ipynb)** 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()`](#quad_vec) can handle this well though." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "raises-exception" + ] + }, + "outputs": [], + "source": [ + "def parametrized_func(s_prime, s):\n", + " return s_prime * (s_prime + s + 1j)\n", + "\n", + "\n", + "s_array = np.linspace(-1, +1, num=10)\n", + "quad_vec(\n", + " lambda x: parametrized_func(x, s=s_array),\n", + " a=0.0,\n", + " b=2.0,\n", + ")" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "There is no good way to write integrals as [SymPy](https://docs.sympy.org) expressions that correctly {func}`~sympy.utilities.lambdify.lambdify` to a **vectorized** integral that handles complex values. Here is a first step however. Note that this integral expression class derives from {class}`sympy.Integral ` and:\n", + "We now attempt to design [SymPy](https://docs.sympy.org) 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 ` and that:\n", "\n", - "1. overwrites its {meth}`~sympy.core.basic.Basic.doit` method, so that the integral cannot be evaluated by SymPy.\n", - "2. provides a custom NumPy printer method (see [TR-001](001.ipynb)) that lambdifies this expression node to [`quadpy.quad()`](https://github.com/sigma-py/quadpy).\n", - "3. adds class variables that can affect the behavior of [`quadpy.quad()`](https://github.com/sigma-py/quadpy)." + "1. overwrites its {meth}`~sympy.core.basic.Basic.doit` method, so that the integral cannot be evaluated by SymPy,\n", + "2. provides a custom NumPy printer method (see **[TR-001](001.ipynb)**) that lambdifies this expression node to [`quadpy.quad()`](https://github.com/sigma-py/quadpy),\n", + "4. adds class variables that allow configuring [`quad_vec()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad_vec.html),\n", + "5. dummifies the integration variable in case it is not a valid Python variable name." ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ - "import sympy as sp\n", - "from sympy.printing.pycode import _unpack_integral_limits\n", - "\n", - "\n", "class UnevaluatableIntegral(sp.Integral):\n", " abs_tolerance = 1e-5\n", " rel_tolerance = 1e-5\n", " limit = 50\n", + " dummify = True\n", "\n", " def doit(self, **hints):\n", " args = [arg.doit(**hints) for arg in self.args]\n", @@ -323,16 +488,21 @@ "\n", " def _numpycode(self, printer, *args):\n", " integration_vars, limits = _unpack_integral_limits(self)\n", - " if len(limits) != 1:\n", + " if len(limits) != 1 or len(integration_vars) != 1:\n", " msg = f\"Cannot handle {len(limits)}-dimensional integrals\"\n", " raise ValueError(msg)\n", - " integrate = \"quadpy_quad\"\n", - " printer.module_imports[\"quadpy\"].update({f\"quad as {integrate}\"})\n", - " limit_str = \"{}, {}\".format(*tuple(map(printer._print, limits[0])))\n", - " args = \", \".join(map(printer._print, integration_vars))\n", - " expr = printer._print(self.args[0])\n", + " x = integration_vars[0]\n", + " a, b = limits[0]\n", + " expr = self.args[0]\n", + " if self.dummify:\n", + " dummy = sp.Dummy()\n", + " expr = expr.xreplace({x: dummy})\n", + " x = dummy\n", + " integrate_func = \"quad_vec\"\n", + " printer.module_imports[\"scipy.integrate\"].add(integrate_func)\n", " return (\n", - " f\"{integrate}(lambda {args}: {expr}, {limit_str},\"\n", + " f\"{integrate_func}(lambda {printer._print(x)}: {printer._print(expr)},\"\n", + " f\" {printer._print(a)}, {printer._print(b)},\"\n", " f\" epsabs={self.abs_tolerance}, epsrel={self.abs_tolerance},\"\n", " f\" limit={self.limit})[0]\"\n", " )" @@ -340,72 +510,187 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "To test whether this works, we write the expression from {ref}`report/016:Vectorized input` as a {class}`sympy.Expr `. Note that the integral indeed does not evaluate when calling {meth}`~sympy.core.basic.Basic.doit`:" + "To test whether this works, test this integral expression on another {func}`~ampform.sympy.unevaluated` expression:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, "outputs": [], "source": [ - "s, s_prime, a, b = sp.symbols(\"s s_p a b\")\n", - "integral_expr: sp.Expr = UnevaluatableIntegral(\n", - " s_prime * (s_prime + s + sp.I),\n", - " (s_prime, a, b),\n", - ")\n", - "integral_expr.doit()" + "@unevaluated\n", + "class MyFunction(sp.Expr):\n", + " x: sp.Symbol\n", + " omega1: sp.Symbol\n", + " omega2: sp.Symbol\n", + " phi1: sp.Symbol\n", + " phi2: sp.Symbol\n", + " _latex_repr_ = R\"f\\left({x}\\right)\"\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " x, omega1, omega2, phi1, phi2 = self.args\n", + " return sp.sin(omega1 * x + phi2) + sp.sin(omega2 * x + phi2)\n", + "\n", + "\n", + "x, omega1, omega2, phi1, phi2 = sp.symbols(\"x omega1 omega2 phi1 phi2\")\n", + "expr = MyFunction(x, omega1, omega2, phi1, phi2)\n", + "Math(aslatex({expr: expr.doit(deep=False)}))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "w, a, b = sp.symbols(\"w a b\")\n", + "fourier_expr = UnevaluatableIntegral(expr * sp.exp(-sp.I * w * x), (x, a, b))\n", + "fourier_expr" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ - "Indeed the expression correctly lambdifies correctly:" + "Indeed the expression correctly lambdifies correctly, despite the {meth}`~sympy.core.basic.Basic.doit` call:" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "outputs": [], "source": [ - "import inspect\n", - "\n", - "integral_func = sp.lambdify([s, a, b], integral_expr)\n", - "src = inspect.getsource(integral_func)\n", - "print(src)" + "func = sp.lambdify([x, omega1, omega2, phi1, phi2], expr.doit())\n", + "fourier_func = sp.lambdify(\n", + " [w, omega1, omega2, phi1, phi2, a, b], fourier_expr.doit()\n", + ")" ] }, { - "cell_type": "markdown", - "metadata": {}, + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], "source": [ - "Note, however, that the lambdified function has to be vectorized before it can handle {obj}`numpy.array`s:" + "src = inspect.getsource(fourier_func)\n", + "src = f\"\"\"```python\n", + "{black.format_str(src, mode=black.FileMode()).strip()}\n", + "```\"\"\"\n", + "Markdown(src)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, "tags": [] }, "outputs": [], "source": [ - "vec_integral_func = np.vectorize(integral_func)\n", - "vec_integral_func(s_array, a=0.0, b=2.0)" + "domain = np.linspace(-7, +7, num=500)\n", + "parameters = dict(\n", + " omega1=1.2,\n", + " omega2=2.3,\n", + " phi1=-1.2,\n", + " phi2=+0.4,\n", + ")\n", + "func_output = func(domain, **parameters)\n", + "fourier_output = fourier_func(domain, **parameters, a=-10, b=+10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "editable": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "%config InlineBackend.figure_formats = ['svg']\n", + "fig, ax = plt.subplots()\n", + "ax.set_xlabel(\"$x,w$\")\n", + "ax.plot(domain, func_output, label=\"$f(x)$\")\n", + "ax.plot(domain, fourier_output.real, label=R\"$\\mathrm{Re}\\,F(w)$\")\n", + "ax.plot(domain, fourier_output.imag, label=R\"$\\mathrm{Im}\\,F(w)$\")\n", + "ax.legend()\n", + "plt.show()" ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, "source": [ ":::{tip}\n", - "\n", - "For a more complicated and challenging expression, see {ref}`report/003:SymPy expressions` in [TR-003](003.ipynb).\n", - "\n", + "See how this integral expression class is applied to the phase space factor in **[TR-003](003.ipynb)**.\n", ":::" ] } From 1f9f9f4cb0ca9a7268d7d8db3216794787c9efef Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Sun, 21 Jan 2024 18:46:12 +0100 Subject: [PATCH 5/5] ENH: support dummy variables in TR-003 --- docs/report/003.ipynb | 46 +++++++++++++++++++++++++------------------ 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/docs/report/003.ipynb b/docs/report/003.ipynb index 5a26a66f..0f3ad90f 100644 --- a/docs/report/003.ipynb +++ b/docs/report/003.ipynb @@ -97,6 +97,9 @@ "jupyter": { "source_hidden": true }, + "mystnb": { + "code_prompt_show": "Import Python libraries" + }, "slideshow": { "slide_type": "" }, @@ -128,9 +131,10 @@ "from ampform.io import aslatex\n", "from ampform.sympy import unevaluated\n", "from ampform.sympy.math import ComplexSqrt\n", - "from IPython.display import SVG, Math, display\n", + "from IPython.display import SVG, Markdown, Math, display\n", "from ipywidgets import FloatSlider\n", "from scipy.integrate import quad_vec\n", + "from sympy.printing.pycode import _unpack_integral_limits\n", "from tqdm.auto import tqdm\n", "\n", "warnings.filterwarnings(\"ignore\")\n", @@ -956,13 +960,11 @@ }, "outputs": [], "source": [ - "from sympy.printing.pycode import _unpack_integral_limits\n", - "\n", - "\n", "class UnevaluatableIntegral(sp.Integral):\n", " abs_tolerance = 1e-5\n", " rel_tolerance = 1e-5\n", " limit = 50\n", + " dummify = True\n", "\n", " def doit(self, **hints):\n", " args = [arg.doit(**hints) for arg in self.args]\n", @@ -970,17 +972,21 @@ "\n", " def _numpycode(self, printer, *args):\n", " integration_vars, limits = _unpack_integral_limits(self)\n", - " if len(limits) != 1:\n", + " if len(limits) != 1 or len(integration_vars) != 1:\n", " msg = f\"Cannot handle {len(limits)}-dimensional integrals\"\n", " raise ValueError(msg)\n", - " integrate = \"quad_vec\"\n", - " printer.module_imports[\"scipy.integrate\"].add(integrate)\n", - " a, b = map(printer._print, limits[0])\n", - " args = \", \".join(map(printer._print, integration_vars))\n", - " expr = printer._print(self.args[0])\n", + " x = integration_vars[0]\n", + " a, b = limits[0]\n", + " expr = self.args[0]\n", + " if self.dummify:\n", + " dummy = sp.Dummy()\n", + " expr = expr.xreplace({x: dummy})\n", + " x = dummy\n", + " integrate_func = \"quad_vec\"\n", + " printer.module_imports[\"scipy.integrate\"].add(integrate_func)\n", " return (\n", - " f\"{integrate}(lambda {args}: {expr},\"\n", - " f\" a={a}, b={b},\"\n", + " f\"{integrate_func}(lambda {printer._print(x)}: {printer._print(expr)},\"\n", + " f\" {printer._print(a)}, {printer._print(b)},\"\n", " f\" epsabs={self.abs_tolerance}, epsrel={self.abs_tolerance},\"\n", " f\" limit={self.limit})[0]\"\n", " )" @@ -1033,8 +1039,7 @@ " )\n", "\n", "\n", - "x = sp.Symbol(\"x\", real=True)\n", - "integral_expr = dispersion_integral(s, m1, m2, angular_momentum=L, s_prime=x)\n", + "integral_expr = dispersion_integral(s, m1, m2, angular_momentum=L, s_prime=s_prime)\n", "integral_expr" ] }, @@ -1075,14 +1080,14 @@ " [s, m1, m2, epsilon],\n", " integral_expr.subs(L, 0).doit(),\n", " # integration symbol should not be extracted as common sub-expression!\n", - " cse=partial(sp.cse, ignore=[x], list=False),\n", + " cse=partial(sp.cse, ignore=[s_prime], list=False),\n", ")\n", "integral_s_wave_func = np.vectorize(integral_s_wave_func)\n", "\n", "integral_p_wave_func = sp.lambdify(\n", " [s, m1, m2, epsilon],\n", " integral_expr.subs(L, 1).doit(),\n", - " cse=partial(sp.cse, ignore=[x], list=False),\n", + " cse=partial(sp.cse, ignore=[s_prime], list=False),\n", ")\n", "integral_p_wave_func = np.vectorize(integral_p_wave_func)" ] @@ -1105,8 +1110,10 @@ "outputs": [], "source": [ "src = inspect.getsource(integral_s_wave_func.pyfunc)\n", - "src = black.format_str(src, mode=black.FileMode())\n", - "print(src)" + "src = f\"\"\"```python\n", + "{black.format_str(src, mode=black.FileMode()).strip()}\n", + "```\"\"\"\n", + "Markdown(src)" ] }, { @@ -1130,7 +1137,8 @@ "editable": true, "slideshow": { "slide_type": "" - } + }, + "tags": [] }, "outputs": [], "source": [