Regression

assets/jitr_logo.ipynb

@@ -149,7 +149,7 @@
     "            eta\n",
     "            * r\n",
     "            * sc.spherical_jn(l, eta * r)\n",
-    "            * sc.sph_harm_y(l,0, np.pi / 2 + theta, 0)\n",
+    "            * sc.sph_harm_y(l, 0, np.pi / 2 + theta, 0)\n",
     "        )\n",
     "        ** 2\n",
     "    )\n",

docs/regression-tests.md

@@ -0,0 +1,47 @@
+# Regression tests
+
+End-to-end tests compare `jitr` against committed outputs from Frescox and
+TALYS. Reference CSVs are committed; neither external code is required in CI.
+
+## Running
+
+```bash
+uv run pytest tests/regression/
+```
+
+## Frescox cases (F1–F8)
+
+Eight elastic-scattering cases against LLNL
+[Frescox](https://github.com/LLNL/Frescox), all for p/n + `78Ni`:
+
+| Case | Projectile | E\_lab (MeV) | Source deck |
+|------|-----------|-------------|-------------|
+| F1   | p | 6.9   | `B1-example-el.out` (block 1) |
+| F2   | p | 11.0  | `B1-example-el.out` (block 2) |
+| F3   | p | 49.35 | `B1-example-el.out` (block 3) |
+| F4   | p | 100   | `B1-high-el.in` (block 1) |
+| F5   | p | 200   | `B1-high-el.in` (block 2) |
+| F6   | n | 49.35 | `B1_n-high-el.in` (block 1) |
+| F7   | n | 100   | `B1_n-high-el.in` (block 2) |
+| F8   | n | 200   | `B1_n-high-el.in` (block 3) |
+
+## TALYS cases (T1–T4)
+
+Four JLM/JLMB elastic-scattering cases against TALYS 2.2 for `120Sn` at
+10 MeV:
+
+| Case | Projectile | `jlmmode` | Notes |
+|------|-----------|-----------|-------|
+| T1   | n | 0 | standard JLMB normalization |
+| T2   | n | 2 | stronger JLMB normalization |
+| T3   | p | 0 | Coulomb code path |
+| T4   | p | 2 | Coulomb + stronger normalization |
+
+```{toctree}
+:maxdepth: 1
+:caption: Harness details
+
+../tests/regression/README
+../tests/regression/frescox/README
+../tests/regression/talys/README
+```

docs/tests.md

@@ -1,5 +1,13 @@
 # Tests
 
+```{toctree}
+:hidden:
+:maxdepth: 1
+:titlesonly:
+
+regression-tests
+```
+
 ## Set up the development environment
 
 Clone the repository and install the full development environment:
@@ -41,6 +49,12 @@ The example notebooks are tested with `pytest` and `nbval`:
 uv run --group examples pytest --nbval-lax examples/notebooks/
 ```
 
+## Run the regression tests
+
+End-to-end regression tests compare `jitr` against committed outputs from
+Frescox and TALYS.  See [Regression tests](regression-tests.md) for the full
+layout, the landed cases, and regeneration instructions.
+
 ## Browse the published examples
 
 The curated notebook subset that appears on the documentation site is

examples/notebooks/tabulated_density_demo.ipynb

@@ -110,13 +110,21 @@
     "rho_p_bskg3, rho_n_bskg3 = density.densities(A, Z, r, model=\"bskg3\")\n",
     "\n",
     "table = density.density_table(A, Z, model=\"d1m\")\n",
-    "proton_number = 4.0 * np.pi * np.trapezoid(\n",
-    "    table.radial_grid**2 * table.proton_density_grid,\n",
-    "    table.radial_grid,\n",
+    "proton_number = (\n",
+    "    4.0\n",
+    "    * np.pi\n",
+    "    * np.trapezoid(\n",
+    "        table.radial_grid**2 * table.proton_density_grid,\n",
+    "        table.radial_grid,\n",
+    "    )\n",
     ")\n",
-    "neutron_number = 4.0 * np.pi * np.trapezoid(\n",
-    "    table.radial_grid**2 * table.neutron_density_grid,\n",
-    "    table.radial_grid,\n",
+    "neutron_number = (\n",
+    "    4.0\n",
+    "    * np.pi\n",
+    "    * np.trapezoid(\n",
+    "        table.radial_grid**2 * table.neutron_density_grid,\n",
+    "        table.radial_grid,\n",
+    "    )\n",
     ")\n",
     "\n",
     "assert np.isclose(proton_number, Z, rtol=1e-4)\n",
@@ -143,23 +151,27 @@
    "source": [
     "fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharey=True)\n",
     "\n",
-    "axes[0].plot(r, rho_p_d1m, color='r', label=\"p\", linewidth=2)\n",
-    "axes[0].plot(r, rho_n_d1m, color='b', label=\"n\", linewidth=2)\n",
-    "axes[0].plot(r, rho_p_d1m + rho_n_d1m, color='g', label=\"matter\", linestyle=\"--\", linewidth=2)\n",
+    "axes[0].plot(r, rho_p_d1m, color=\"r\", label=\"p\", linewidth=2)\n",
+    "axes[0].plot(r, rho_n_d1m, color=\"b\", label=\"n\", linewidth=2)\n",
+    "axes[0].plot(\n",
+    "    r, rho_p_d1m + rho_n_d1m, color=\"g\", label=\"matter\", linestyle=\"--\", linewidth=2\n",
+    ")\n",
     "axes[0].set_title(\"D1M\")\n",
     "axes[0].set_xlabel(r\"$r$ [fm]\")\n",
     "axes[0].set_ylabel(r\"$\\rho(r)$ [fm$^{-3}$]\")\n",
     "axes[0].legend()\n",
     "\n",
-    "axes[1].plot(r, rho_p_bskg3, color='r', label=\"p\", linewidth=2)\n",
-    "axes[1].plot(r, rho_n_bskg3, color='b', label=\"n\", linewidth=2)\n",
-    "axes[1].plot(r, rho_p_bskg3 + rho_n_bskg3, color='g', label=\"total\", linestyle=\"--\", linewidth=2)\n",
+    "axes[1].plot(r, rho_p_bskg3, color=\"r\", label=\"p\", linewidth=2)\n",
+    "axes[1].plot(r, rho_n_bskg3, color=\"b\", label=\"n\", linewidth=2)\n",
+    "axes[1].plot(\n",
+    "    r, rho_p_bskg3 + rho_n_bskg3, color=\"g\", label=\"total\", linestyle=\"--\", linewidth=2\n",
+    ")\n",
     "axes[1].set_title(\"BSKG3\")\n",
     "axes[1].set_xlabel(r\"$r$ [fm]\")\n",
     "axes[1].legend()\n",
     "\n",
     "fig.suptitle(r\"$^{16}$O\", fontsize=20)\n",
-    "fig.tight_layout()\n"
+    "fig.tight_layout()"
    ]
   },
   {
@@ -183,9 +195,9 @@
    "outputs": [],
    "source": [
     "def norm_and_rms_radii(rho, R):\n",
-    "    m0 = 4*np.pi * np.trapezoid(rho * R**2, R)\n",
-    "    m2 = 4*np.pi * np.trapezoid(rho * R**4, R)\n",
-    "    return m0, np.sqrt(m2 / m0)     "
+    "    m0 = 4 * np.pi * np.trapezoid(rho * R**2, R)\n",
+    "    m2 = 4 * np.pi * np.trapezoid(rho * R**4, R)\n",
+    "    return m0, np.sqrt(m2 / m0)"
    ]
   },
   {
@@ -198,18 +210,19 @@
     "def print_density_stats(rho_p, rho_n, R):\n",
     "    N, rms_n = norm_and_rms_radii(rho_n, R)\n",
     "    Z, rms_p = norm_and_rms_radii(rho_p, R)\n",
-    "    Z, rms_m = norm_and_rms_radii(rho_p+rho_n,R)\n",
+    "    Z, rms_m = norm_and_rms_radii(rho_p + rho_n, R)\n",
     "    print(f\"N: {N:1.2f}\")\n",
     "    print(f\"Z: {Z:1.2f}\")\n",
     "    print(f\"rms_n: {rms_n:1.2f} fm\")\n",
     "    print(f\"rms_p: {rms_p:1.2f} fm\")\n",
     "    print(f\"rms_m: {rms_m:1.2f} fm\")\n",
     "    print(f\"nskin: {rms_n - rms_p:1.3f} fm\")\n",
     "\n",
+    "\n",
     "def plot_densities(ax, rho_p, rho_n, R):\n",
-    "    ax.plot(R, rho_n, color='b', label=\"n\", linewidth=2)\n",
-    "    ax.plot(R, rho_p, color='r', label=\"p\", linewidth=2)\n",
-    "    ax.plot(R, rho_n + rho_p, color='g', label=\"total\", linewidth=2)"
+    "    ax.plot(R, rho_n, color=\"b\", label=\"n\", linewidth=2)\n",
+    "    ax.plot(R, rho_p, color=\"r\", label=\"p\", linewidth=2)\n",
+    "    ax.plot(R, rho_n + rho_p, color=\"g\", label=\"total\", linewidth=2)"
    ]
   },
   {
@@ -261,7 +274,7 @@
     "print(\"208-Pb\")\n",
     "print(\"\\n2pf:\")\n",
     "print(\"=======================\")\n",
-    "print_density_stats( np.array(rho_p_2pf(R)), np.array(rho_n_2pf(R)), R)\n",
+    "print_density_stats(np.array(rho_p_2pf(R)), np.array(rho_n_2pf(R)), R)\n",
     "\n",
     "print(\"\\nd1m:\")\n",
     "print(\"=======================\")\n",
@@ -290,8 +303,8 @@
     }
    ],
    "source": [
-    "fig, axes = plt.subplots( 3, 1, figsize=(6, 7), sharex=True)\n",
-    "fig.suptitle(r\"$^{208}$Pb\", fontsize=20 )\n",
+    "fig, axes = plt.subplots(3, 1, figsize=(6, 7), sharex=True)\n",
+    "fig.suptitle(r\"$^{208}$Pb\", fontsize=20)\n",
     "plot_densities(axes[0], np.array(rho_p_2pf(R)), np.array(rho_n_2pf(R)), R)\n",
     "plot_densities(axes[1], rho_p_d1m, rho_n_d1m, R)\n",
     "plot_densities(axes[2], rho_p_bskg3, rho_n_bskg3, R)\n",
@@ -307,9 +320,9 @@
     "axes[1].set_title(r\"D1M\")\n",
     "axes[2].set_title(r\"BSKG3\")\n",
     "\n",
-    "axes[0].set_xlim([-0.1,11])\n",
+    "axes[0].set_xlim([-0.1, 11])\n",
     "plt.tight_layout()\n",
-    "plt.show()\n"
+    "plt.show()"
    ]
   },
   {

pyproject.toml

@@ -48,6 +48,7 @@ examples = [
     "jupyter",
     "pytest",
     "matplotlib",
+    "black[jupyter]>=26.3.1",
 ]
 
 lint = [

src/jitr/folding/folding.py

@@ -48,6 +48,20 @@ def __init__(self, r_max: float = 20.0, n_quad: int = 200) -> None:
         nodes, weights = leggauss(self.n_quad)
         self.r_q = np.asarray(0.5 * self.r_max * (nodes + 1.0), dtype=float)
         self.w_q = np.asarray(0.5 * self.r_max * weights, dtype=float)
+        self._D = self._build_diff_matrix()
+
+    def _build_diff_matrix(self) -> FloatArray:
+        r = self.r_q
+        w = self.w_q
+        idx = np.arange(len(r))
+        lam = np.sqrt(r * (self.r_max - r) * w)
+        signs = (-1.0) ** (idx[:, None] + idx[None, :])
+        lam_ratio = lam[None, :] / lam[:, None]
+        r_diff = r[:, None] - r[None, :]
+        with np.errstate(divide="ignore", invalid="ignore"):
+            D = np.where(r_diff != 0.0, signs * lam_ratio / r_diff, 0.0)
+        np.fill_diagonal(D, -D.sum(axis=1))
+        return D
 
     def interp_to_quad(self, r_grid: ArrayLike, f_grid: ArrayLike) -> FloatArray:
         """Interpolate tabulated data onto the quadrature grid.
@@ -74,6 +88,29 @@ def integrate(self, f_q: ArrayLike) -> float:
 
         return float(np.sum(self.w_q * np.asarray(f_q, float)))
 
+    def differentiate(self, f_q: ArrayLike) -> FloatArray:
+        """Differentiate a quantity sampled on the quadrature grid.
+
+        Uses the Lagrange–Gauss-Legendre differentiation matrix (Baye 2015).
+        Exact for polynomials of degree ≤ ``n_quad`` − 1; requires no
+        grid-uniformity assumption.
+
+        The off-diagonal elements are derived from the Lagrange-Legendre
+        barycentric weights λ̃_j = √(r_j (r_max − r_j) w_j):
+
+            D_{ij} = (−1)^(i+j) λ̃_j / (λ̃_i (r_i − r_j))   i ≠ j
+
+        The diagonal is set by the row-sum condition ∑_j D_{ij} = 0,
+        which holds because d/dr [∑_j l_j(r)] = 0 for any Lagrange basis.
+
+        Args:
+            f_q: Function values sampled on ``self.r_q``.
+
+        Returns:
+            Derivative df/dr sampled on ``self.r_q``.
+        """
+        return self._D @ np.asarray(f_q, dtype=float)
+
     def Z_from_density(self, rho_q: ArrayLike) -> float:
         """Compute the particle number implied by a spherical density.