From 08ade55f9ab95a29b4638e575632fc55e7a58447 Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Fri, 26 Jan 2024 16:42:52 +0100
Subject: [PATCH 01/12] add unevaluatableintegral definition

https://compwa.github.io/report/016\#sympy-integral
---
 .cspell.json                  |  4 ++++
 src/ampform/sympy/__init__.py | 33 +++++++++++++++++++++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/.cspell.json b/.cspell.json
index 83f931cec..9bc386b8a 100644
--- a/.cspell.json
+++ b/.cspell.json
@@ -80,9 +80,12 @@
     "dtype",
     "dummified",
     "dummifies",
+    "dummify",
     "einsum",
     "elif",
     "epem",
+    "epsabs",
+    "epsrel",
     "eqnarray",
     "eval",
     "evalf",
@@ -273,6 +276,7 @@
     "tensorwaves",
     "toctree",
     "topness",
+    "unevaluatable",
     "unitarity",
     "venv",
     "weisskopf",
diff --git a/src/ampform/sympy/__init__.py b/src/ampform/sympy/__init__.py
index e0c9db2de..c42c15be9 100644
--- a/src/ampform/sympy/__init__.py
+++ b/src/ampform/sympy/__init__.py
@@ -27,6 +27,7 @@
 import sympy as sp
 from sympy.printing.conventions import split_super_sub
 from sympy.printing.precedence import PRECEDENCE
+from sympy.printing.pycode import _unpack_integral_limits
 
 from ._decorator import (
     ExprClass,  # noqa: F401  # pyright: ignore[reportUnusedImport]
@@ -350,3 +351,35 @@ def _warn_about_unsafe_hash():
     """
     message = dedent(message).replace("\n", " ").strip()
     _LOGGER.warning(message)
+
+
+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]"
+        )

From 1c1ad626f3a72230c733437f3b39e11378535176 Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Fri, 26 Jan 2024 16:58:59 +0100
Subject: [PATCH 02/12] add a test file example

---
 tests/sympy/test_integral.py | 2 ++
 1 file changed, 2 insertions(+)
 create mode 100644 tests/sympy/test_integral.py

diff --git a/tests/sympy/test_integral.py b/tests/sympy/test_integral.py
new file mode 100644
index 000000000..733bc9d2a
--- /dev/null
+++ b/tests/sympy/test_integral.py
@@ -0,0 +1,2 @@
+def test_my_function():
+    assert 1 + 1 == 2

From e3c372ab0ad490bb15da11e7ef64e11ef835fa43 Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Wed, 31 Jan 2024 12:11:56 +0100
Subject: [PATCH 03/12] write simple test real value integral

---
 pyproject.toml               |  1 +
 tests/sympy/test_integral.py | 13 +++++++++++--
 2 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index df518f736..6f38c8f68 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -117,6 +117,7 @@ test = [
     "pytest-cov",
     "pytest-profiling",
     "pytest-xdist",
+    "scipy",
 ]
 viz = ["graphviz"]
 
diff --git a/tests/sympy/test_integral.py b/tests/sympy/test_integral.py
index 733bc9d2a..ce425201d 100644
--- a/tests/sympy/test_integral.py
+++ b/tests/sympy/test_integral.py
@@ -1,2 +1,11 @@
-def test_my_function():
-    assert 1 + 1 == 2
+import sympy as sp
+
+from ampform.sympy import UnevaluatableIntegral
+
+
+class TestUnevaluatableIntegral:
+    def test_real_value_function(self):
+        x = sp.symbols("x")
+        integral_expr = UnevaluatableIntegral(x**2, (x, 1, 3))
+        func = sp.lambdify(args=[], expr=integral_expr)
+        assert func() == 26 / 3

From ae20f10ad8281cc94d4324719eae35e2d4509d8d Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Wed, 31 Jan 2024 13:42:36 +0100
Subject: [PATCH 04/12] write simple test array value parameter function

---
 tests/sympy/test_integral.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/tests/sympy/test_integral.py b/tests/sympy/test_integral.py
index ce425201d..282035750 100644
--- a/tests/sympy/test_integral.py
+++ b/tests/sympy/test_integral.py
@@ -1,3 +1,4 @@
+import pytest
 import sympy as sp
 
 from ampform.sympy import UnevaluatableIntegral
@@ -9,3 +10,10 @@ def test_real_value_function(self):
         integral_expr = UnevaluatableIntegral(x**2, (x, 1, 3))
         func = sp.lambdify(args=[], expr=integral_expr)
         assert func() == 26 / 3
+
+    def test_array_value_parameter_function(self):
+        x, p = sp.symbols("x,p")
+        integral_expr = UnevaluatableIntegral(x**p, (x, 1, 3))
+        func = sp.lambdify(args=[p], expr=integral_expr)
+        assert func(p=2) == 26 / 3
+        assert pytest.approx(func(p=1)) == 4

From b189f47ed7505a2c7c5569e80b564b530a7cfbcc Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Wed, 31 Jan 2024 13:53:52 +0100
Subject: [PATCH 05/12] parametrize test

---
 .vscode/settings.json        |  4 ++--
 tests/sympy/test_integral.py | 12 +++++++++---
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/.vscode/settings.json b/.vscode/settings.json
index e74f30354..70c643c6b 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -63,9 +63,9 @@
   "ruff.enable": true,
   "ruff.organizeImports": true,
   "search.exclude": {
-    "typings/**": true,
     "**/tests/**/__init__.py": true,
-    ".constraints/*.txt": true
+    ".constraints/*.txt": true,
+    "typings/**": true
   },
   "yaml.schemas": {
     "https://citation-file-format.github.io/1.2.0/schema.json": "CITATION.cff",
diff --git a/tests/sympy/test_integral.py b/tests/sympy/test_integral.py
index 282035750..49f01d913 100644
--- a/tests/sympy/test_integral.py
+++ b/tests/sympy/test_integral.py
@@ -11,9 +11,15 @@ def test_real_value_function(self):
         func = sp.lambdify(args=[], expr=integral_expr)
         assert func() == 26 / 3
 
-    def test_array_value_parameter_function(self):
+    @pytest.mark.parametrize(
+        "p_value,expected",
+        [
+            (2, 26 / 3),
+            (1, 4),
+        ],
+    )
+    def test_array_value_parameter_function(self, p_value, expected):
         x, p = sp.symbols("x,p")
         integral_expr = UnevaluatableIntegral(x**p, (x, 1, 3))
         func = sp.lambdify(args=[p], expr=integral_expr)
-        assert func(p=2) == 26 / 3
-        assert pytest.approx(func(p=1)) == 4
+        assert pytest.approx(func(p=p_value)) == expected

From 69b926ae6913102dcef37e10d00533b823170b44 Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Wed, 31 Jan 2024 14:38:53 +0100
Subject: [PATCH 06/12] add complex function and array input

---
 tests/sympy/test_integral.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/tests/sympy/test_integral.py b/tests/sympy/test_integral.py
index 49f01d913..8d66423b9 100644
--- a/tests/sympy/test_integral.py
+++ b/tests/sympy/test_integral.py
@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 import sympy as sp
 
@@ -16,9 +17,14 @@ def test_real_value_function(self):
         [
             (2, 26 / 3),
             (1, 4),
+            (1j, (1 / 2 - 1j / 2) * (-1 + 3 ** (1 + 1j))),
+            (
+                np.array([0, 0.5, 1, 2]),
+                np.array([2, 2 * 3 ** (1 / 2) - 2 / 3, 4, 8 + 2 / 3]),
+            ),
         ],
     )
-    def test_array_value_parameter_function(self, p_value, expected):
+    def test_vectorized_parameter_function(self, p_value, expected):
         x, p = sp.symbols("x,p")
         integral_expr = UnevaluatableIntegral(x**p, (x, 1, 3))
         func = sp.lambdify(args=[p], expr=integral_expr)

From 5b6a95d532eadbc840fc789e2b3b4eff23562049 Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Wed, 31 Jan 2024 14:48:29 +0100
Subject: [PATCH 07/12] revert setting change

---
 .vscode/settings.json | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.vscode/settings.json b/.vscode/settings.json
index 70c643c6b..e74f30354 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -63,9 +63,9 @@
   "ruff.enable": true,
   "ruff.organizeImports": true,
   "search.exclude": {
+    "typings/**": true,
     "**/tests/**/__init__.py": true,
-    ".constraints/*.txt": true,
-    "typings/**": true
+    ".constraints/*.txt": true
   },
   "yaml.schemas": {
     "https://citation-file-format.github.io/1.2.0/schema.json": "CITATION.cff",

From df35f6c3027ba14bdf5394cb240922fde0356e15 Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Mon, 5 Feb 2024 17:49:37 +0100
Subject: [PATCH 08/12] add integral example section

---
 .cspell.json           |   1 +
 docs/usage/sympy.ipynb | 128 ++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 128 insertions(+), 1 deletion(-)

diff --git a/.cspell.json b/.cspell.json
index 9bc386b8a..09a5a24ff 100644
--- a/.cspell.json
+++ b/.cspell.json
@@ -206,6 +206,7 @@
     "xlim",
     "xreplace",
     "xticks",
+    "xytext",
     "yaxis",
     "ylabel",
     "ylim",
diff --git a/docs/usage/sympy.ipynb b/docs/usage/sympy.ipynb
index 80369b35d..8c9d07907 100644
--- a/docs/usage/sympy.ipynb
+++ b/docs/usage/sympy.ipynb
@@ -213,6 +213,132 @@
     "Math(aslatex({e: e.doit() for e in exprs}))"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Integral"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Integral over $x^p$ with symbolic range from a to b\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "import sympy as sp\n",
+    "\n",
+    "from ampform.sympy import UnevaluatableIntegral\n",
+    "\n",
+    "x, p, a, b = sp.symbols(\"x,p,a,b\")\n",
+    "integral_expr = UnevaluatableIntegral(x**p, (x, a, b))\n",
+    "integral_expr.doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "integral_func = sp.lambdify(args=[p, a, b], expr=integral_expr)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "value_a = 1.2\n",
+    "value_b = 3.6\n",
+    "value_p = np.array([0.4, 0.6, 0.8])\n",
+    "\n",
+    "areas = integral_func(value_p, value_a, value_b)\n",
+    "areas"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input",
+     "scroll-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "x_area = np.linspace(value_a, value_b, num=100)\n",
+    "x_line = np.linspace(0, 4, num=100)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.set_xlabel(\"$x$\")\n",
+    "ax.set_ylabel(\"$x^p$\")\n",
+    "ax.plot(x_line, x_line ** value_p[0], label=\"$p=0.4$\", color=\"C0\")\n",
+    "ax.fill_between(x_area, x_area ** value_p[0], alpha=1, color=\"C0\")\n",
+    "\n",
+    "ax.plot(x_line, x_line ** value_p[1], label=\"$p=0.6$\", color=\"C1\")\n",
+    "ax.fill_between(x_area, x_area ** value_p[1], alpha=0.6, color=\"C1\")\n",
+    "\n",
+    "ax.plot(x_line, x_line ** value_p[2], label=\"$p=0.8$\", color=\"C2\")\n",
+    "ax.fill_between(x_area, x_area ** value_p[2], alpha=0.7, color=\"C2\")\n",
+    "\n",
+    "ax.text(\n",
+    "    (value_a + value_b) / 2,\n",
+    "    ((value_a ** value_p[0] + value_b ** value_p[0]) / 2) * 0.5,\n",
+    "    \"Area\",\n",
+    "    horizontalalignment=\"center\",\n",
+    "    verticalalignment=\"center\",\n",
+    ")\n",
+    "ax.annotate(\n",
+    "    \"a\", (value_a, 0.08), textcoords=\"offset points\", xytext=(0, -15), ha=\"center\"\n",
+    ")\n",
+    "ax.annotate(\n",
+    "    \"b\", (value_b, 0.08), textcoords=\"offset points\", xytext=(0, -15), ha=\"center\"\n",
+    ")\n",
+    "\n",
+    "ax.legend()\n",
+    "plt.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -275,7 +401,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,

From a726ce696feff527dd920d2bebb67a5c157790a5 Mon Sep 17 00:00:00 2001
From: Vitor Shen <17490173+shenvitor@users.noreply.github.com>
Date: Tue, 6 Feb 2024 18:41:25 +0100
Subject: [PATCH 09/12] add compare sp.Integral and UnevaluatableIntegral

---
 docs/usage/sympy.ipynb | 144 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 144 insertions(+)

diff --git a/docs/usage/sympy.ipynb b/docs/usage/sympy.ipynb
index 8c9d07907..1f07d6d04 100644
--- a/docs/usage/sympy.ipynb
+++ b/docs/usage/sympy.ipynb
@@ -339,6 +339,150 @@
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There are numerous functions whose integrals cannot be expressed in terms of known special functions and for which no closed-form solutions have been identified yet. Here are a few examples:\n",
+    "\n",
+    "1. **Integral Involving Complex Exponential and Polynomials:**\n",
+    "   $$ \\int e^{x + e^{-x}} dx $$\n",
+    "   This integral combines an exponential function with a nested exponential term, leading to a complexity that resists simplification in terms of known functions.\n",
+    "\n",
+    "2. **Integral of a Rational Function Times an Exponential:**\n",
+    "   $$ \\int \\frac{e^x}{x^2 + 1} dx $$\n",
+    "   This integral involves the exponential function and a non-trivial rational function, which together do not yield to standard integration techniques.\n",
+    "\n",
+    "3. **Integral of an Exponential of a Trigonometric Function:**\n",
+    "   $$ \\int e^{\\sin x} dx $$\n",
+    "   The combination of an exponential function and a sinusoidal function is another case where known methods of integration in terms of elementary or special functions fail.\n",
+    "\n",
+    "4. **Product of a Polynomial and a Trigonometric Function:**\n",
+    "   $$ \\int x^4 \\sin(x^2) dx $$\n",
+    "   This integral combines a polynomial with a trigonometric function in a non-standard way, making it resistant to analytical methods.\n",
+    "\n",
+    "5. **Integral of a Logarithmic Function Inside an Exponential:**\n",
+    "   $$ \\int e^{\\ln(x)/x} dx $$\n",
+    "   The combination of a logarithmic function inside an exponential function, as seen here, creates a complex relationship not easily resolved by known functions.\n",
+    "\n",
+    "In the context of hadron physics and high-energy physics, encountering such challenging integrals is not uncommon. They can arise in theoretical models, complex scattering problems, or in the analysis of experimental data. In such cases, physicists often resort to numerical integration, computational approximations, or the development of new special functions or methods tailored to the problem at hand. These integrals, while not solvable in a traditional sense, can still provide significant insights into the behaviors and properties of the systems being studied."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Some example test to compare ``sp.Integral`` and ``UnevaluatableIntegral``"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "UnevaluatableIntegral(sp.exp(-(x**4)), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp.Integral(sp.exp(-(x**4)), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "UnevaluatableIntegral(sp.sin(x) / x, (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp.Integral(sp.sin(x) / x, (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "UnevaluatableIntegral(sp.exp(x + sp.exp(-x)), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp.Integral(sp.exp(x + sp.exp(-x)), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "UnevaluatableIntegral(sp.exp(sp.log(x) / x), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp.Integral(sp.exp(sp.log(x) / x), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "UnevaluatableIntegral(sp.exp(sp.sin(x)), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp.Integral(sp.exp(sp.sin(x)), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "UnevaluatableIntegral(sp.exp(x) / (x**2 + 1), (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp.Integral(sp.exp(x) / (x**2 + 1), (x, a, b)).doit()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

From 1239fe72e4d41ba0fb982ea6886772f720345577 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Fri, 9 Feb 2024 10:40:10 +0100
Subject: [PATCH 10/12] DX: avoid IPython debug warnings

---
 .cspell.json    | 1 +
 environment.yml | 1 +
 2 files changed, 2 insertions(+)

diff --git a/.cspell.json b/.cspell.json
index 09a5a24ff..53c4c1fac 100644
--- a/.cspell.json
+++ b/.cspell.json
@@ -42,6 +42,7 @@
     "Dalitzplot",
     "MAINT",
     "Minkowski",
+    "PYDEVD",
     "adrs",
     "aitchison",
     "arange",
diff --git a/environment.yml b/environment.yml
index fa1316535..8fb56b271 100644
--- a/environment.yml
+++ b/environment.yml
@@ -8,5 +8,6 @@ dependencies:
   - pip:
       - -c .constraints/py3.11.txt -e .[dev]
 variables:
+  PYDEVD_DISABLE_FILE_VALIDATION: 1
   PRETTIER_LEGACY_CLI: "1"
   PYTHONHASHSEED: 0

From df08529011968dbd16ed06562055a898f2dc4b26 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Sun, 11 Feb 2024 12:28:18 +0100
Subject: [PATCH 11/12] DOC: explain integral examples

---
 docs/usage/sympy.ipynb | 247 ++++++++++++++---------------------------
 1 file changed, 86 insertions(+), 161 deletions(-)

diff --git a/docs/usage/sympy.ipynb b/docs/usage/sympy.ipynb
index 1f07d6d04..bc35b5cbf 100644
--- a/docs/usage/sympy.ipynb
+++ b/docs/usage/sympy.ipynb
@@ -217,14 +217,37 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Integral"
+    "## Numerical integrals"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Integral over $x^p$ with symbolic range from a to b\n"
+    "In hadron physics and high-energy physics, it often happens that models contain integrals that do not have an analytical solution.. They can arise in theoretical models, complex scattering problems, or in the analysis of experimental data. In such cases, we need to resort to numerical integrations.\n",
+    "\n",
+    "SymPy provides the [`sympy.Integral`](https://docs.sympy.org/latest/modules/integrals/integrals.html#sympy.integrals.integrals.Integral) class, but this does not give us control over whether or not we want to avoid integrating the class analytically. An example of such an analytically unsolvable integral is shown below. Note that the integral does not evaluate despite the `doit()` call."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sympy as sp\n",
+    "\n",
+    "x, a, b = sp.symbols(\"x a b\")\n",
+    "p = sp.Symbol(\"p\", positive=True)\n",
+    "integral_expr = sp.Integral(sp.exp(x) / (x**p + 1), (x, a, b))\n",
+    "integral_expr.doit()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For amplitude models that contain such integrals that should not be solved analytically, AmpForm provides the {class}`.UnevaluatableIntegral` class. It functions in the same way as [`sympy.Integral`](https://docs.sympy.org/latest/modules/integrals/integrals.html#sympy.integrals.integrals.Integral), but prevents the class from evaluating at all, even if the integral can be solved analytically."
    ]
   },
   {
@@ -239,13 +262,27 @@
    },
    "outputs": [],
    "source": [
-    "import sympy as sp\n",
-    "\n",
     "from ampform.sympy import UnevaluatableIntegral\n",
     "\n",
-    "x, p, a, b = sp.symbols(\"x,p,a,b\")\n",
-    "integral_expr = UnevaluatableIntegral(x**p, (x, a, b))\n",
-    "integral_expr.doit()"
+    "UnevaluatableIntegral(x**p, (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sp.Integral(x**p, (x, a, b)).doit()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This allows {class}`.UnevaluatableIntegral` to serve as a placeholder in expression trees that we call `doit` on when lambdifying to a numerical function. The resulting numerical function takes **complex-valued** and **multidimensional arrays** as function arguments.\n",
+    "\n",
+    "In the following, we see an example where the parameter $p$ inside the integral gets an array as input."
    ]
   },
   {
@@ -260,6 +297,7 @@
    },
    "outputs": [],
    "source": [
+    "integral_expr = UnevaluatableIntegral(sp.exp(x) / (x**p + 1), (x, a, b))\n",
     "integral_func = sp.lambdify(args=[p, a, b], expr=integral_expr)"
    ]
   },
@@ -277,11 +315,11 @@
    "source": [
     "import numpy as np\n",
     "\n",
-    "value_a = 1.2\n",
-    "value_b = 3.6\n",
-    "value_p = np.array([0.4, 0.6, 0.8])\n",
+    "a_val = 1.2\n",
+    "b_val = 3.6\n",
+    "p_array = np.array([0.4, 0.6, 0.8])\n",
     "\n",
-    "areas = integral_func(value_p, value_a, value_b)\n",
+    "areas = integral_func(p_array, a_val, b_val)\n",
     "areas"
    ]
   },
@@ -306,34 +344,27 @@
     "%config InlineBackend.figure_formats = ['svg']\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
-    "x_area = np.linspace(value_a, value_b, num=100)\n",
+    "x_area = np.linspace(a_val, b_val, num=100)\n",
     "x_line = np.linspace(0, 4, num=100)\n",
     "\n",
     "fig, ax = plt.subplots()\n",
     "ax.set_xlabel(\"$x$\")\n",
     "ax.set_ylabel(\"$x^p$\")\n",
-    "ax.plot(x_line, x_line ** value_p[0], label=\"$p=0.4$\", color=\"C0\")\n",
-    "ax.fill_between(x_area, x_area ** value_p[0], alpha=1, color=\"C0\")\n",
     "\n",
-    "ax.plot(x_line, x_line ** value_p[1], label=\"$p=0.6$\", color=\"C1\")\n",
-    "ax.fill_between(x_area, x_area ** value_p[1], alpha=0.6, color=\"C1\")\n",
-    "\n",
-    "ax.plot(x_line, x_line ** value_p[2], label=\"$p=0.8$\", color=\"C2\")\n",
-    "ax.fill_between(x_area, x_area ** value_p[2], alpha=0.7, color=\"C2\")\n",
+    "for i, p_val in enumerate(p_array):\n",
+    "    ax.plot(x_line, x_line**p_val, label=f\"$p={p_val}$\", c=f\"C{i}\")\n",
+    "    ax.fill_between(x_area, x_area**p_val, alpha=(0.7 - i * 0.2), color=\"C0\")\n",
     "\n",
     "ax.text(\n",
-    "    (value_a + value_b) / 2,\n",
-    "    ((value_a ** value_p[0] + value_b ** value_p[0]) / 2) * 0.5,\n",
-    "    \"Area\",\n",
+    "    x=(a_val + b_val) / 2,\n",
+    "    y=((a_val ** p_array[0] + b_val ** p_array[0]) / 2) * 0.5,\n",
+    "    s=\"Area\",\n",
     "    horizontalalignment=\"center\",\n",
     "    verticalalignment=\"center\",\n",
     ")\n",
-    "ax.annotate(\n",
-    "    \"a\", (value_a, 0.08), textcoords=\"offset points\", xytext=(0, -15), ha=\"center\"\n",
-    ")\n",
-    "ax.annotate(\n",
-    "    \"b\", (value_b, 0.08), textcoords=\"offset points\", xytext=(0, -15), ha=\"center\"\n",
-    ")\n",
+    "text_kwargs = dict(ha=\"center\", textcoords=\"offset points\", xytext=(0, -15))\n",
+    "ax.annotate(\"a\", (a_val, 0.08), **text_kwargs)\n",
+    "ax.annotate(\"b\", (b_val, 0.08), **text_kwargs)\n",
     "\n",
     "ax.legend()\n",
     "plt.show()"
@@ -341,151 +372,45 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "There are numerous functions whose integrals cannot be expressed in terms of known special functions and for which no closed-form solutions have been identified yet. Here are a few examples:\n",
-    "\n",
-    "1. **Integral Involving Complex Exponential and Polynomials:**\n",
-    "   $$ \\int e^{x + e^{-x}} dx $$\n",
-    "   This integral combines an exponential function with a nested exponential term, leading to a complexity that resists simplification in terms of known functions.\n",
-    "\n",
-    "2. **Integral of a Rational Function Times an Exponential:**\n",
-    "   $$ \\int \\frac{e^x}{x^2 + 1} dx $$\n",
-    "   This integral involves the exponential function and a non-trivial rational function, which together do not yield to standard integration techniques.\n",
-    "\n",
-    "3. **Integral of an Exponential of a Trigonometric Function:**\n",
-    "   $$ \\int e^{\\sin x} dx $$\n",
-    "   The combination of an exponential function and a sinusoidal function is another case where known methods of integration in terms of elementary or special functions fail.\n",
-    "\n",
-    "4. **Product of a Polynomial and a Trigonometric Function:**\n",
-    "   $$ \\int x^4 \\sin(x^2) dx $$\n",
-    "   This integral combines a polynomial with a trigonometric function in a non-standard way, making it resistant to analytical methods.\n",
-    "\n",
-    "5. **Integral of a Logarithmic Function Inside an Exponential:**\n",
-    "   $$ \\int e^{\\ln(x)/x} dx $$\n",
-    "   The combination of a logarithmic function inside an exponential function, as seen here, creates a complex relationship not easily resolved by known functions.\n",
-    "\n",
-    "In the context of hadron physics and high-energy physics, encountering such challenging integrals is not uncommon. They can arise in theoretical models, complex scattering problems, or in the analysis of experimental data. In such cases, physicists often resort to numerical integration, computational approximations, or the development of new special functions or methods tailored to the problem at hand. These integrals, while not solvable in a traditional sense, can still provide significant insights into the behaviors and properties of the systems being studied."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Some example test to compare ``sp.Integral`` and ``UnevaluatableIntegral``"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "UnevaluatableIntegral(sp.exp(-(x**4)), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sp.Integral(sp.exp(-(x**4)), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "UnevaluatableIntegral(sp.sin(x) / x, (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sp.Integral(sp.sin(x) / x, (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "UnevaluatableIntegral(sp.exp(x + sp.exp(-x)), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sp.Integral(sp.exp(x + sp.exp(-x)), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "UnevaluatableIntegral(sp.exp(sp.log(x) / x), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sp.Integral(sp.exp(sp.log(x) / x), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "UnevaluatableIntegral(sp.exp(sp.sin(x)), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sp.Integral(sp.exp(sp.sin(x)), (x, a, b)).doit()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "source": [
-    "UnevaluatableIntegral(sp.exp(x) / (x**2 + 1), (x, a, b)).doit()"
+    "The arrays can be complex-valued as well. This is particularly useful when calculating dispersion integrals (see **[TR-003](https://compwa.github.io/report/003#general-dispersion-integral)**)."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "outputs": [],
    "source": [
-    "sp.Integral(sp.exp(x) / (x**2 + 1), (x, a, b)).doit()"
+    "integral_func(\n",
+    "    p=np.array([1.5 - 8.6j, -4.6 + 5.5j]),\n",
+    "    a=a_val,\n",
+    "    b=b_val,\n",
+    ")"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": []
+   },
    "source": [
     "## Summations"
    ]

From dbb7ba473163a32a0e77fac4def8c7fc65c2082e Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Sun, 11 Feb 2024 12:41:19 +0100
Subject: [PATCH 12/12] ENH: emit warning if SciPy is not installed

---
 pyproject.toml                | 14 +++++++++++---
 src/ampform/sympy/__init__.py | 13 +++++++++++++
 2 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 6f38c8f68..0a346585f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -50,7 +50,10 @@ name = "ampform"
 requires-python = ">=3.7"
 
 [project.optional-dependencies]
-all = ["ampform[viz]"]
+all = [
+    "ampform[scipy]",
+    "ampform[viz]",
+]
 dev = [
     "ampform[all]",
     "ampform[doc]",
@@ -62,7 +65,7 @@ dev = [
 ]
 doc = [
     "Sphinx >=3",
-    "ampform[viz]",
+    "ampform[all]",
     "black",
     "ipympl",
     "matplotlib",
@@ -102,6 +105,7 @@ mypy = [
     "mypy >=0.730",
     "sphinx-api-relink >=0.0.3",
 ]
+scipy = ["scipy"]
 sty = [
     "ampform[format]",
     "ampform[lint]",
@@ -109,6 +113,7 @@ sty = [
     "pre-commit >=1.4.0",
 ]
 test = [
+    "ampform[scipy]",
     "black",
     "ipywidgets", # symplot
     "nbmake",
@@ -117,7 +122,6 @@ test = [
     "pytest-cov",
     "pytest-profiling",
     "pytest-xdist",
-    "scipy",
 ]
 viz = ["graphviz"]
 
@@ -185,6 +189,10 @@ warn_unused_configs = true
 ignore_missing_imports = true
 module = ["graphviz.*"]
 
+[[tool.mypy.overrides]]
+ignore_missing_imports = true
+module = ["scipy.*"]
+
 [[tool.mypy.overrides]]
 ignore_missing_imports = true
 module = ["ipywidgets.*"]
diff --git a/src/ampform/sympy/__init__.py b/src/ampform/sympy/__init__.py
index c42c15be9..169b6baf7 100644
--- a/src/ampform/sympy/__init__.py
+++ b/src/ampform/sympy/__init__.py
@@ -19,6 +19,7 @@
 import os
 import pickle
 import re
+import warnings
 from abc import abstractmethod
 from os.path import abspath, dirname, expanduser
 from textwrap import dedent
@@ -364,6 +365,7 @@ def doit(self, **hints):
         return self.func(*args)
 
     def _numpycode(self, printer, *args):
+        _warn_if_scipy_not_installed()
         integration_vars, limits = _unpack_integral_limits(self)
         if len(limits) != 1 or len(integration_vars) != 1:
             msg = f"Cannot handle {len(limits)}-dimensional integrals"
@@ -383,3 +385,14 @@ def _numpycode(self, printer, *args):
             f" epsabs={self.abs_tolerance}, epsrel={self.abs_tolerance},"
             f" limit={self.limit})[0]"
         )
+
+
+def _warn_if_scipy_not_installed() -> None:
+    try:
+        import scipy  # noqa: F401  # pyright: ignore[reportUnusedImport, reportMissingImports]
+    except ImportError:
+        warnings.warn(
+            "Scipy is not installed. Install with 'pip install scipy' or with 'pip"
+            " install ampform[scipy]'",
+            stacklevel=1,
+        )