From 04ca79deb0422aca187e036fd44389c39fa59fc6 Mon Sep 17 00:00:00 2001
From: Pietro D'Antuono <pietro.dantuono@24sea.eu>
Date: Mon, 8 Jun 2026 17:00:03 +0200
Subject: [PATCH 1/4] perf(#23): Optimize Numba + Python 3.14 compat

---
 .github/workflows/cd.yml                      |   2 +-
 .github/workflows/ci.yml                      |   2 +-
 .github/workflows/pages.yml                   |   2 +-
 .gitignore                                    |   4 +-
 CHANGELOG.rst                                 |  11 +
 ...0.0.0 - PD.Generic.development-guide.ipynb | 589 +++++++++------
 numba-speedup.md                              | 152 ++++
 py_fatigue/__init__.py                        |  23 +
 py_fatigue/cycle_count/__init__.py            |   3 +
 py_fatigue/cycle_count/cycle_count.py         |   8 +-
 py_fatigue/cycle_count/rainflow.py            |   2 +-
 py_fatigue/cycle_count/utils.py               |  38 +-
 py_fatigue/damage/crack_growth.py             |  12 +-
 py_fatigue/damage/stress_life.py              |  40 +-
 py_fatigue/material/crack_growth_curve.py     | 215 ++++--
 py_fatigue/material/sn_curve.py               |  59 +-
 py_fatigue/mean_stress/corrections.py         |  25 +-
 py_fatigue/testing.py                         |   7 +-
 py_fatigue/utils.py                           | 287 +++++++-
 py_fatigue/version.py                         |   1 +
 pyproject.toml                                |   7 +-
 rst_docs/api/material/crack_growth_curve.rst  |  16 +-
 rst_docs/api/material/sn_curve.rst            |  12 +-
 rst_docs/user/01-absolute-noob.rst            |   4 +-
 rst_docs/user/additional/glossary.rst         |   5 +-
 rst_docs/user/examples/02-sn_curve.rst        |   8 +-
 rst_docs/user/examples/03-cg_curve.rst        |   9 +-
 rst_docs/user/examples/04-cycle_counting.rst  |   2 +-
 rst_docs/user/examples/06-mean_stress.rst     |   2 +-
 scripts/benchmark_numba_speedups.py           | 677 ++++++++++++++++++
 scripts/numba_benchmark_results.json          | 326 +++++++++
 tests/damage/test_stress_life.py              | 210 +++---
 tests/material/test_crack_growth_curve.py     | 328 ++++++---
 tests/mean_stress/test_corrections.py         | 195 +++--
 tests/test_utils.py                           | 199 ++++-
 35 files changed, 2818 insertions(+), 664 deletions(-)
 create mode 100644 numba-speedup.md
 create mode 100644 scripts/benchmark_numba_speedups.py
 create mode 100644 scripts/numba_benchmark_results.json

diff --git a/.github/workflows/cd.yml b/.github/workflows/cd.yml
index 6fe2590..f766947 100644
--- a/.github/workflows/cd.yml
+++ b/.github/workflows/cd.yml
@@ -43,7 +43,7 @@ jobs:
     - name: Install uv
       uses: astral-sh/setup-uv@v5
     - name: Set up Python
-      run: uv python install 3.13
+      run: uv python install 3.14
     - name: Bump version
       run: |
         uv run bumpversion.py ${{ github.event.inputs.bumptype }}
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 47a7b41..876ebff 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -10,7 +10,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.10", "3.11", "3.12", "3.13"]
+        python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
         os: [ubuntu-latest]
     runs-on: ${{ matrix.os }}
     steps:
diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml
index d3a2909..e20d566 100644
--- a/.github/workflows/pages.yml
+++ b/.github/workflows/pages.yml
@@ -16,7 +16,7 @@ jobs:
       - name: Install uv
         uses: astral-sh/setup-uv@v5
         with:
-          python-version: "3.13"
+          python-version: "3.14"
           enable-cache: true
       - name: Set up Python
         run: uv python install
diff --git a/.gitignore b/.gitignore
index 27f83ae..3c8fb2f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -105,4 +105,6 @@ uv.lock
 
 **data/
 
-pkg_graph.*
\ No newline at end of file
+pkg_graph.*
+
+/notebooks/tutorials-repository/
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index a82eb4c..926d54b 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,3 +1,14 @@
+Unreleased
+----------
+- **perf(numba)**: cache the generic Numba root-finder specializations, add
+  opt-in `warmup_numba()`, and enable disk caching on the direct Numba kernels.
+- **build(numba)**: widen Python < 3.14 support to Numba 0.61 through 0.65
+  without gaps, while keeping Python 3.14 on Numba 0.65.
+- **fix(compatibility)**: add Pandas 3.x compatibility fixes in the
+  cycle-count aggregation helpers and notebook examples.
+- **docs(numba)**: refresh the Numba speedup report, repair broken docs links,
+  and update the tutorial notebooks so they execute on Python 3.13 and 3.14.
+
 2.0.4
 ------
 - **bump(patch)**: Bump version to 2.0.4 [skip ci]
diff --git a/notebooks/0.0.0 - PD.Generic.development-guide.ipynb b/notebooks/0.0.0 - PD.Generic.development-guide.ipynb
index 2176078..5bd1a12 100644
--- a/notebooks/0.0.0 - PD.Generic.development-guide.ipynb	
+++ b/notebooks/0.0.0 - PD.Generic.development-guide.ipynb	
@@ -1,222 +1,373 @@
 {
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "toc": true
-   },
-   "source": [
-    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
-    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Import-the-source-code-of-the-project\" data-toc-modified-id=\"Import-the-source-code-of-the-project-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Import the source code of the project</a></span></li><li><span><a href=\"#Import-and-use-the-functionality\" data-toc-modified-id=\"Import-and-use-the-functionality-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Import and use the functionality</a></span></li></ul></div>"
-   ]
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3.13.13 (py-fatigue; uv)",
+      "name": "Python 3.13.13 (py-fatigue; uv)",
+      "language": "python"
+    },
+    "language_info": {
+      "name": "python",
+      "version": "3.13.1",
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "file_extension": ".py",
+      "mimetype": "text/x-python"
+    },
+    "toc": {
+      "base_numbering": 1,
+      "nav_menu": {},
+      "number_sections": true,
+      "sideBar": true,
+      "skip_h1_title": true,
+      "title_cell": "Table of Contents",
+      "title_sidebar": "Contents",
+      "toc_cell": true,
+      "toc_position": {},
+      "toc_section_display": true,
+      "toc_window_display": false
+    },
+    "latex_envs": {
+      "LaTeX_envs_menu_present": true,
+      "autoclose": false,
+      "autocomplete": true,
+      "bibliofile": "biblio.bib",
+      "cite_by": "apalike",
+      "current_citInitial": 1,
+      "eqLabelWithNumbers": true,
+      "eqNumInitial": 1,
+      "hotkeys": {
+        "equation": "Ctrl-E",
+        "itemize": "Ctrl-I"
+      },
+      "labels_anchors": false,
+      "latex_user_defs": false,
+      "report_style_numbering": false,
+      "user_envs_cfg": false
+    },
+    "hide_input": false
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Development with the notebook\n",
-    "This notebook server can be used as a scratchpad. Work in progress code, algorithms, visualizations, tests, ...\n",
-    "\n",
-    "The ultimate goal here is to refactor these snippets in well documented and tested python classes and methods so they can be imported easily. \n",
-    "\n",
-    "The suggested workflow is as follows:\n",
-    "1. Use the notebook to prototype and test drive your development\n",
-    "2. Migrate your developed code into the project. Create a python file, organize it in classes,...\n",
-    "3. Document your migrated code (docstring) so it's functionality is documented in the sphinx docs.\n",
-    "4. Import the migrated code back into this notebook and develop tests around it.\n",
-    "5. Include these tests in the project.\n",
-    "6. Cleanup your notebook as the functionality is now in the python project and can be easily imported in the future.\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Import the source code of the project\n",
-    "Updates to the project's source need to be reflected here so it's possible to add new functionality and fix bugs which are directly reflected here. Jupyter supports this behaviour with the `autoreload` magic. \n",
-    "\n",
-    "The 2 cells below will setup this autoreload functionality (so the original import is not cached) and add the project path to the system path. Now the project source can be imported (without installing it as a python package)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2022-06-16T14:47:08.852940Z",
-     "start_time": "2022-06-16T14:47:08.805966Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "%load_ext autoreload\n",
-    "%autoreload 2\n",
-    "# %matplotlib inline"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2022-06-16T14:47:10.411668Z",
-     "start_time": "2022-06-16T14:47:09.429114Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "import datetime as dt\n",
-    "import os\n",
-    "import platform\n",
-    "import sys\n",
-    "import numpy as np\n",
-    "import pandas as pd\n",
-    "import matplotlib.pyplot as plt\n",
-    "import plotly.graph_objs as go\n",
-    "\n",
-    "PROJECT_PATH = os.path.dirname(os.getcwd())\n",
-    "print(f'PROJECT_PATH = {PROJECT_PATH}')\n",
-    "\n",
-    "if not PROJECT_PATH in sys.path:\n",
-    "    sys.path.append(PROJECT_PATH)\n",
-    "    \n",
-    "print('Platform:', platform.platform())\n",
-    "print('Python version:', sys.version)\n",
-    "print('numpy version:', np.__version__)\n",
-    "np.set_printoptions(threshold=sys.maxsize)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2022-06-16T14:47:10.475159Z",
-     "start_time": "2022-06-16T14:47:10.411668Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "#matplotlib tweaking\n",
-    "plt.rcParams['figure.figsize'] = (10.5/2.514/0.7, 4.5/2.514/0.7)\n",
-    "plt.rcParams[\"font.family\"] = \"Sans\"\n",
-    "plt.rcParams[\"font.size\"] = 10\n",
-    "# plt.rcParams[\"axes.grid\"] = True\n",
-    "# plt.rcParams['grid.color'] = \"#CCCCCC\"\n",
-    "# plt.rcParams['grid.linestyle'] = \"-\"\n",
-    "# plt.rcParams['grid.color'] = \"#DDDDDD\"\n",
-    "# plt.rcParams['grid.linestyle'] = \"-\"\n",
-    "plt.rcParams['axes.spines.right'] = False\n",
-    "plt.rcParams['axes.spines.top'] = False\n",
-    "plt.rcParams['lines.markersize'] = 3\n",
-    "plt.rcParams['xtick.bottom'] = False\n",
-    "plt.rcParams['xtick.labelbottom'] = True\n",
-    "plt.rcParams['ytick.left'] = False\n",
-    "plt.rcParams['ytick.labelleft'] = True"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Import and use the functionality"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2022-06-16T14:47:12.468692Z",
-     "start_time": "2022-06-16T14:47:11.748517Z"
+  "nbformat": 4,
+  "nbformat_minor": 5,
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "id": "b638a54f-cae2-4232-8d56-863cec9784c8",
+      "metadata": {
+        "toc": true
+      },
+      "source": [
+        "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
+        "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Import-the-source-code-of-the-project\" data-toc-modified-id=\"Import-the-source-code-of-the-project-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Import the source code of the project</a></span></li><li><span><a href=\"#Import-and-use-the-functionality\" data-toc-modified-id=\"Import-and-use-the-functionality-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Import and use the functionality</a></span></li></ul></div>\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "bcec5489-4436-4c2a-b49b-6c22e639887d",
+      "metadata": {},
+      "source": [
+        "# Development with the notebook\n",
+        "This notebook server can be used as a scratchpad. Work in progress code, algorithms, visualizations, tests, ...\n",
+        "\n",
+        "The ultimate goal here is to refactor these snippets in well documented and tested python classes and methods so they can be imported easily. \n",
+        "\n",
+        "The suggested workflow is as follows:\n",
+        "1. Use the notebook to prototype and test drive your development\n",
+        "2. Migrate your developed code into the project. Create a python file, organize it in classes,...\n",
+        "3. Document your migrated code (docstring) so it's functionality is documented in the sphinx docs.\n",
+        "4. Import the migrated code back into this notebook and develop tests around it.\n",
+        "5. Include these tests in the project.\n",
+        "6. Cleanup your notebook as the functionality is now in the python project and can be easily imported in the future.\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "38b73e6c-e8f2-473b-b4f6-79a35d3aa1ef",
+      "metadata": {},
+      "source": [
+        "## Import the source code of the project\n",
+        "Updates to the project's source need to be reflected here so it's possible to add new functionality and fix bugs which are directly reflected here. Jupyter supports this behaviour with the `autoreload` magic. \n",
+        "\n",
+        "The 2 cells below will setup this autoreload functionality (so the original import is not cached) and add the project path to the system path. Now the project source can be imported (without installing it as a python package).\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "id": "7503552a-3170-45d2-8424-bee080872ebe",
+      "metadata": {
+        "ExecuteTime": {
+          "end_time": "2022-06-16T14:47:08.852940Z",
+          "start_time": "2022-06-16T14:47:08.805966Z"
+        }
+      },
+      "execution_count": 1,
+      "source": [
+        "%load_ext autoreload\n",
+        "%autoreload 2\n",
+        "# %matplotlib inline\n"
+      ],
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "id": "01a8057c-ccb2-45ee-98ee-24e66b1aaf9c",
+      "metadata": {
+        "ExecuteTime": {
+          "end_time": "2022-06-16T14:47:10.411668Z",
+          "start_time": "2022-06-16T14:47:09.429114Z"
+        }
+      },
+      "execution_count": 2,
+      "source": [
+        "import datetime as dt\n",
+        "import os\n",
+        "import platform\n",
+        "import sys\n",
+        "import numpy as np\n",
+        "import pandas as pd\n",
+        "import matplotlib.pyplot as plt\n",
+        "import plotly.graph_objs as go\n",
+        "\n",
+        "PROJECT_PATH = os.path.dirname(os.getcwd())\n",
+        "print(f'PROJECT_PATH = {PROJECT_PATH}')\n",
+        "\n",
+        "if not PROJECT_PATH in sys.path:\n",
+        "    sys.path.append(PROJECT_PATH)\n",
+        "    \n",
+        "print('Platform:', platform.platform())\n",
+        "print('Python version:', sys.version)\n",
+        "print('numpy version:', np.__version__)\n",
+        "np.set_printoptions(threshold=sys.maxsize)\n"
+      ],
+      "outputs": [
+        {
+          "output_type": "stream",
+          "name": "stdout",
+          "text": [
+            "PROJECT_PATH = /home/pietro/Code\n",
+            "Platform: Linux-6.17.0-35-generic-x86_64-with-glibc2.39\n",
+            "Python version: 3.13.13 (main, May 10 2026, 19:26:54) [Clang 22.1.3 ]\n",
+            "numpy version: 2.2.6\n"
+          ]
+        }
+      ]
+    },
+    {
+      "cell_type": "code",
+      "id": "7b33381d-3393-4d19-9370-e0b6a5cf6709",
+      "metadata": {
+        "ExecuteTime": {
+          "end_time": "2022-06-16T14:47:10.475159Z",
+          "start_time": "2022-06-16T14:47:10.411668Z"
+        }
+      },
+      "execution_count": 3,
+      "source": [
+        "#matplotlib tweaking\n",
+        "plt.rcParams['figure.figsize'] = (10.5/2.514/0.7, 4.5/2.514/0.7)\n",
+        "plt.rcParams[\"font.family\"] = \"Sans\"\n",
+        "plt.rcParams[\"font.size\"] = 10\n",
+        "# plt.rcParams[\"axes.grid\"] = True\n",
+        "# plt.rcParams['grid.color'] = \"#CCCCCC\"\n",
+        "# plt.rcParams['grid.linestyle'] = \"-\"\n",
+        "# plt.rcParams['grid.color'] = \"#DDDDDD\"\n",
+        "# plt.rcParams['grid.linestyle'] = \"-\"\n",
+        "plt.rcParams['axes.spines.right'] = False\n",
+        "plt.rcParams['axes.spines.top'] = False\n",
+        "plt.rcParams['lines.markersize'] = 3\n",
+        "plt.rcParams['xtick.bottom'] = False\n",
+        "plt.rcParams['xtick.labelbottom'] = True\n",
+        "plt.rcParams['ytick.left'] = False\n",
+        "plt.rcParams['ytick.labelleft'] = True\n"
+      ],
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "id": "887c6147-f370-4d86-bac9-c064d36693c0",
+      "metadata": {},
+      "source": [
+        "## Import and use the functionality\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "id": "7a121539-3276-4659-8476-50fa1f4cdd5e",
+      "metadata": {
+        "ExecuteTime": {
+          "end_time": "2022-06-16T14:47:12.468692Z",
+          "start_time": "2022-06-16T14:47:11.748517Z"
+        }
+      },
+      "execution_count": 4,
+      "source": [
+        "# Use the package\n",
+        "import py_fatigue as pf\n",
+        "import py_fatigue.cycle_count.rainflow as rf\n",
+        "import py_fatigue.cycle_count.histogram as ht\n",
+        "from py_fatigue.version import parse_version, __version__\n",
+        "v = parse_version(__version__)\n"
+      ],
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "id": "e77007e0-be2e-4d6a-940a-973696d0e664",
+      "metadata": {},
+      "execution_count": 5,
+      "source": [
+        "print(v)\n"
+      ],
+      "outputs": [
+        {
+          "output_type": "stream",
+          "name": "stdout",
+          "text": [
+            "Version(major=2, minor=1, patch=0, release=None, num=None)\n"
+          ]
+        }
+      ]
+    },
+    {
+      "cell_type": "code",
+      "id": "fa6c3c46-48fb-46c9-8bc5-d79995a97bbf",
+      "metadata": {},
+      "execution_count": 6,
+      "source": [
+        "\n",
+        "s = np.array([4,7,2,10,5,9,3,4,2,12,5,11,1,4,3,10,6,12,4,8,1,9,4,6], dtype=float)\n",
+        "# s = [0, 5, 0, 5, 0, 5, 0, 0 , 5, 0, 5, 0, 5, 0]\n",
+        "cc_4p = rf.rainflow(s, method='fourpoint', extended_output=True)\n",
+        "cc_astm = rf.rainflow(s, method='astm', extended_output=True)\n",
+        "display(cc_4p)\n",
+        "display(cc_astm)\n"
+      ],
+      "outputs": [
+        {
+          "output_type": "display_data",
+          "data": {
+            "text/plain": [
+              "(array([[ 2. ,  7. ,  1. ,  4. ,  2. ],\n",
+              "        [ 0.5,  3.5,  1. ,  6. ,  2. ],\n",
+              "        [ 4. ,  6. ,  1. ,  3. , 10. ],\n",
+              "        [ 3. ,  8. ,  1. , 10. ,  2. ],\n",
+              "        [ 0.5,  3.5,  1. , 13. ,  2. ],\n",
+              "        [ 2. ,  8. ,  1. , 15. ,  2. ],\n",
+              "        [ 2. ,  6. ,  1. , 18. ,  2. ],\n",
+              "        [ 5.5,  6.5,  1. , 12. , 10. ],\n",
+              "        [ 1.5,  5.5,  0.5,  0. ,  2. ],\n",
+              "        [ 2.5,  4.5,  0.5,  1. ,  2. ],\n",
+              "        [ 5. ,  7. ,  0.5,  2. , 14. ],\n",
+              "        [ 5.5,  6.5,  0.5,  9. , 22. ],\n",
+              "        [ 4. ,  5. ,  0.5, 20. ,  2. ],\n",
+              "        [ 2.5,  6.5,  0.5, 21. ,  2. ],\n",
+              "        [ 1. ,  5. ,  0.5, 22. ,  2. ]]),\n",
+              " array([ 4.,  7.,  2., 12.,  1.,  9.,  4.,  6.]),\n",
+              " array([ 0,  1,  2,  9, 20, 21, 22, 23]))\n"
+            ]
+          },
+          "metadata": {}
+        },
+        {
+          "output_type": "display_data",
+          "data": {
+            "text/plain": [
+              "(array([[ 1.5,  5.5,  0.5,  0. ,  2. ],\n",
+              "        [ 2.5,  4.5,  0.5,  1. ,  2. ],\n",
+              "        [ 2. ,  7. ,  1. ,  4. ,  2. ],\n",
+              "        [ 0.5,  3.5,  1. ,  6. ,  2. ],\n",
+              "        [ 4. ,  6. ,  0.5,  2. ,  2. ],\n",
+              "        [ 4. ,  6. ,  0.5,  3. , 10. ],\n",
+              "        [ 3. ,  8. ,  1. , 10. ,  2. ],\n",
+              "        [ 5. ,  7. ,  0.5,  8. ,  2. ],\n",
+              "        [ 0.5,  3.5,  1. , 13. ,  2. ],\n",
+              "        [ 2. ,  8. ,  1. , 15. ,  2. ],\n",
+              "        [ 5.5,  6.5,  0.5,  9. ,  6. ],\n",
+              "        [ 2. ,  6. ,  1. , 18. ,  2. ],\n",
+              "        [ 5.5,  6.5,  0.5, 12. , 10. ],\n",
+              "        [ 5.5,  6.5,  0.5, 17. ,  6. ],\n",
+              "        [ 4. ,  5. ,  0.5, 20. ,  2. ],\n",
+              "        [ 2.5,  6.5,  0.5, 21. ,  2. ],\n",
+              "        [ 1. ,  5. ,  0.5, 22. ,  2. ]]),\n",
+              " array([ 4.,  7.,  2., 10.,  2., 12.,  1., 12.,  1.,  9.,  4.,  6.]),\n",
+              " array([ 0,  1,  2,  3,  8,  9, 12, 17, 20, 21, 22, 23]))\n"
+            ]
+          },
+          "metadata": {}
+        }
+      ]
+    },
+    {
+      "cell_type": "code",
+      "id": "d1fab177-b5c4-4b47-9959-c51449545052",
+      "metadata": {},
+      "execution_count": 14,
+      "source": [
+        "from py_fatigue.material import sn_curve as sn\n",
+        "#SN curves definition\n",
+        "b1_c = sn.SNCurve(\n",
+        "    slope=3,                       # SNCurve can handle multiple slope and\n",
+        "    intercept=12.436,              # intercept types, as long as their sizes \n",
+        "    norm='DNVGL-RP-C203',          # are compatible. The slope and intercept \n",
+        "    environment='Free corrosion',  # attribute will be stored as numpy arrays.\n",
+        "    curve='B1'\n",
+        ")\n",
+        "c_c = sn.SNCurve(\n",
+        "    [3],\n",
+        "    12.115,\n",
+        "    norm='DNVGL-RP-C203',\n",
+        "    environment='Free corrosion',\n",
+        "    curve='C'\n",
+        ")\n",
+        "e_c = sn.SNCurve(\n",
+        "    3,\n",
+        "    intercept=[11.533,],\n",
+        "    norm='DNVGL-RP-C203',\n",
+        "    environment='Free corrosion',\n",
+        "    curve='E'\n",
+        ")\n",
+        "w3_c = sn.SNCurve(\n",
+        "    [3],\n",
+        "    intercept=[10.493,],\n",
+        "    norm='DNVGL-RP-C203',\n",
+        "    environment='Free corrosion',\n",
+        "    curve=\"W3-C\"\n",
+        ")\n",
+        "\n",
+        "#%%\n"
+      ],
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "id": "ab3977ff-4419-4126-b776-37d43bf86ea6",
+      "metadata": {},
+      "execution_count": 12,
+      "source": [
+        "c_c\n"
+      ],
+      "outputs": [
+        {
+          "output_type": "display_data",
+          "data": {
+            "text/plain": [
+              "<py_fatigue.material.sn_curve.SNCurve at 0x77defd731bd0>\n"
+            ]
+          },
+          "metadata": {}
+        }
+      ]
+    },
+    {
+      "cell_type": "code",
+      "id": "e9d1caef-3e70-4a10-9045-8a2254bb0b8c",
+      "metadata": {},
+      "execution_count": null,
+      "source": [],
+      "outputs": []
     }
-   },
-   "outputs": [],
-   "source": [
-    "# Use the package\n",
-    "import py_fatigue as pf\n",
-    "import py_fatigue.cycle_count.rainflow as rf\n",
-    "import py_fatigue.cycle_count.histogram as ht\n",
-    "from py_fatigue.version import parse_version, __version__\n",
-    "v = parse_version(__version__)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "print(v)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "\n",
-    "s = np.array([4,7,2,10,5,9,3,4,2,12,5,11,1,4,3,10,6,12,4,8,1,9,4,6], dtype=float)\n",
-    "# s = [0, 5, 0, 5, 0, 5, 0, 0 , 5, 0, 5, 0, 5, 0]\n",
-    "cc_4p = rf.rainflow(s, method='fourpoint', extended_output=True)\n",
-    "cc_astm = rf.rainflow(s, method='astm', extended_output=True)\n",
-    "display(cc_4p)\n",
-    "display(cc_astm)"
-   ]
-  }
- ],
- "metadata": {
-  "hide_input": false,
-  "kernelspec": {
-   "display_name": "py-fatigue (3.13.1)",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.13.1"
-  },
-  "latex_envs": {
-   "LaTeX_envs_menu_present": true,
-   "autoclose": false,
-   "autocomplete": true,
-   "bibliofile": "biblio.bib",
-   "cite_by": "apalike",
-   "current_citInitial": 1,
-   "eqLabelWithNumbers": true,
-   "eqNumInitial": 1,
-   "hotkeys": {
-    "equation": "Ctrl-E",
-    "itemize": "Ctrl-I"
-   },
-   "labels_anchors": false,
-   "latex_user_defs": false,
-   "report_style_numbering": false,
-   "user_envs_cfg": false
-  },
-  "toc": {
-   "base_numbering": 1,
-   "nav_menu": {},
-   "number_sections": true,
-   "sideBar": true,
-   "skip_h1_title": true,
-   "title_cell": "Table of Contents",
-   "title_sidebar": "Contents",
-   "toc_cell": true,
-   "toc_position": {},
-   "toc_section_display": true,
-   "toc_window_display": false
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
+  ]
+}
\ No newline at end of file
diff --git a/numba-speedup.md b/numba-speedup.md
new file mode 100644
index 0000000..4971d7f
--- /dev/null
+++ b/numba-speedup.md
@@ -0,0 +1,152 @@
+# Numba speedup investigation
+
+This is an explanation/reference report for maintainers of `py_fatigue`.
+The goal is to identify where numba compilation is paid repeatedly and where
+small refactors could improve performance without changing fatigue logic.
+
+## Package overview
+
+`py_fatigue` provides tools for fatigue analysis:
+
+- `cycle_count`: rainflow counting, turning points, histograms, and cycle-count
+  objects.
+- `material`: S-N curves and Paris/Walker crack-growth curves.
+- `damage`: Palmgren-Miner, non-linear stress-life damage, Theil stress-life
+  damage, and crack-growth integration.
+- `mean_stress`: DNVGL, Walker/SWT, and Goodman-Haigh mean-stress correction.
+- `geometry`: crack geometry models and stress-intensity geometry factors.
+- `stress_range`, `testing`, and `utils`: supporting data structures,
+  numerical helpers, plotting helpers, and root finders.
+
+The numba usage is concentrated in material curve evaluation, crack-growth
+geometry/SIF calculations, Theil damage, rainflow crossing detection, and the
+generic bisection/Newton wrappers in `py_fatigue.utils`.
+
+## Benchmark setup
+
+I added:
+
+- `scripts/benchmark_numba_speedups.py`: rich CLI benchmark and numba inventory.
+- `scripts/numba_benchmark_results.json`: raw benchmark output from the main run.
+
+The script uses `rich` and `rich-argparse`. The environment was created and
+managed with `uv`:
+
+```bash
+uv sync --no-default-groups --group dev
+uv run python scripts/benchmark_numba_speedups.py \
+  --calls 8 \
+  --repeats 3 \
+  --json-output scripts/numba_benchmark_results.json
+```
+
+The timings below are median wall-clock timings on this branch with
+Python 3.13 + Numba 0.65.1. They are intended to compare call patterns, not
+to be absolute performance claims. Compatibility validation was also run on
+Numba 0.61.2, 0.62.1, 0.64.0, and 0.65.1 on Python 3.13, plus Python 3.14
+with Numba 0.65.1.
+
+## Numba inventory
+
+| Area | Symbol | Mechanism | Assessment |
+| --- | --- | --- | --- |
+| `cycle_count.rainflow` | `_findcross` | `@njit` with explicit signature and `cache=True` | Compiles once per process/signature and can reuse numba's disk cache. |
+| `damage.crack_growth` | `CalcCrackGrowth` | `@jitclass(spec)` | Compiles on first matching instantiation. Warm calls reuse compiled methods. |
+| `damage.crack_growth` | `get_geometry_factor`, `get_sif` | `@nb.njit(..., cache=True)` | Already cached and signature-pinned. |
+| `damage.stress_life` | `_calc_theil_sn_damage` | `@nb.njit(cache=True)` | Already cached. |
+| `geometry.cylinder` | `f_hol_cyl_01` | `@nb.njit(..., cache=True)` | Already cached and signature-pinned. |
+| `material.crack_growth_curve` | `_calc_growth_rate`, `_calc_sif` | `@nb.njit(cache=True)` | Warm calls reuse compilation and can reuse numba's disk cache. |
+| `material.sn_curve` | `_calc_cycles`, `_calc_cycles_2`, `_calc_stress` | `@nb.njit(cache=True)` | Warm calls reuse compilation and can reuse numba's disk cache. `_calc_cycles_2` is used by `SNCurve.get_cycles`; `_calc_cycles` and `_calc_stress` appear to be legacy/unused alternatives. |
+| `material.sn_curve` | `_calc_stress_2` | `@nb.njit(cache=True)` | Already cached and used by `SNCurve.get_stress`. |
+| `material.sn_curve` | `__jit_sn_curve_residuals` | cached `compile_specialized_bisect(...)` | Fixed: the specialized bisection function is a numba dispatcher with a nopython bisection loop. |
+| `mean_stress.corrections` | `__jit_goodman_equation` | `compile_specialized_newton(...)` | OK: returns a numba dispatcher, so `numba_newton(...)` does not recompile on each call. |
+| `utils` | `numba_bisect`, `numba_newton` | generic wrappers with cached specialization | Fixed: repeated calls with the same Python residual reuse the same specialized root finder in-process. Bisection now executes its loop in nopython mode when JIT is enabled. |
+| `utils` / package API | `warmup_numba` | opt-in warmup function | New: compiles common numba call paths for the current process. |
+
+## Benchmark results
+
+| Case | Current pattern | Hoisted/cached pattern | Observed impact |
+| --- | ---: | ---: | ---: |
+| `numba_bisect(py_func)` | 0.002 ms/call | 0.001 ms/call | Repeated-call overhead removed; loop runs in nopython mode. |
+| `numba_newton(py_func)` | 0.002 ms/call | 0.001 ms/call | Repeated-call overhead removed. |
+| `SNCurve.get_cycles` | 0.414 ms/call | N/A | Warm dispatcher is reused. |
+| `SNCurve.get_stress` | 0.448 ms/call | N/A | Warm dispatcher is reused. |
+| `ParisCurve.get_growth_rate` | 1.273 ms/call | N/A | Warm dispatcher is reused; compatibility guards for older Numba do not affect the 0.65 path. |
+| `ParisCurve.get_sif` | 2.177 ms/call | N/A | Warm dispatcher is reused; compatibility guards for older Numba do not affect the 0.65 path. |
+| `rainflow.findcross` | 0.126 ms/call | N/A | Warm dispatcher is reused via cached signature. |
+| `rainflow.findtp` | 0.168 ms/call | N/A | Warm dispatcher is reused via `findcross`. |
+| `f_hol_cyl_01` | 0.003 ms/call | N/A | Warm dispatcher is reused and cached. |
+| `crack_growth.get_sif` | 0.006 ms/call | N/A | Warm dispatcher is reused and cached. |
+| `CalcCrackGrowth` construction | 0.717 ms/call | N/A | Jitclass is warm after first instantiation. |
+| `calc_theil_sn_damage` | 0.024 ms/call | N/A | Warm dispatcher is reused and cached. |
+| `goodman_haigh_mean_stress_correction` | 0.823 ms/call | N/A | The existing module-level Newton specialization avoids repeated compilation. |
+| `find_sn_curve_intersection` | 0.066 ms/call | N/A | Fixed; no longer fails with the bisection argument limit. |
+
+The original root-finder numbers showed the issue clearly: when a Python
+residual was sent through `numba_bisect` or `numba_newton` on every solve, the
+wrapper created a new numba dispatcher every time. The implemented cache removes
+that repeated compilation for repeated calls with the same residual function.
+
+## Findings
+
+1. The direct `@njit` functions are not recompiling on every warm call. The
+   retained direct dispatchers now use `cache=True` where practical, so they can
+   also reuse numba's disk cache across processes.
+2. `numba_newton` is safe in the Goodman-Haigh public path because
+   `__jit_goodman_equation` is a numba dispatcher, not a plain Python function.
+3. The S-N intersection path is fixed. It now calls the already-specialized
+   bisection function directly, and the bisection wrapper no longer imposes the
+   previous "three extra arguments" limit.
+4. The generic wrappers now cache specialized root finders by residual function
+   identity. Any loop that calls `numba_bisect(py_func, ...)` or
+   `numba_newton(py_func, ...)` repeatedly with the same function avoids repeated
+   compilation.
+5. The bisection algorithm itself now runs in nopython mode when JIT is enabled,
+   while preserving scalar-like residual return compatibility for scalar,
+   one-item tuple, one-item array, one-item list, and numpy scalar returns.
+6. Python 3.13 compatibility is now validated on Numba 0.61.2, 0.62.1, 0.64.0,
+   and 0.65.1. Python 3.14 is validated on Numba 0.65.1. The oldest supported
+   line (0.61) needs a small compatibility fallback for the low-level
+   Paris/Walker kernel symbols, but the public curve APIs and the targeted test
+   surface stay green.
+
+## Notebook and docs validation
+
+All tutorial notebooks in `notebooks/tutorials-repository/notebooks` now
+execute on Python 3.13 and 3.14.
+
+Key fixes:
+
+1. The `0.5.x` use-case notebooks no longer hard-fail when the private CSV data
+   is absent; they fall back to a small synthetic dataset and use Pandas 3.x
+   compatible aggregation aliases.
+2. The `0.6.x` nonlinear-damage notebooks were aligned with the current
+   `py_fatigue` API and had stale scratch appendix cells trimmed or updated so
+   automated execution is reliable.
+3. `1.0.0 - AM.General.Usage.ipynb` now uses a valid `python3` kernelspec and
+   Pandas 3.x compatible positional indexing examples.
+
+## Recommended refactors
+
+These are performance refactors only; they do not require changing the fatigue
+equations or output logic.
+
+1. Use `py_fatigue.warmup_numba()` at application startup when predictable
+   runtime latency is more important than startup time.
+2. Keep the root-finder residual functions stable. The cache is keyed by the
+   Python function object, so repeatedly creating new lambdas/functions will
+   still create new specializations.
+3. For maximum reuse of compiled root finders, define residual functions once at
+   module scope instead of creating new lambdas/functions inside tight loops.
+
+## Priority
+
+The highest-value root-finder fixes have been applied: repeated calls now reuse
+specializations, bisection executes in nopython mode when JIT is enabled, and
+`find_sn_curve_intersection` no longer fails under JIT. `warmup_numba()` provides
+opt-in in-process warmup. Full "always precompiled" behavior is only partially
+possible: signature-pinned `cache=True` dispatchers can reuse numba's disk cache
+across processes, but jitclasses and dynamic root-finder closures still need
+in-process warmup after each process starts. The remaining compatibility work is
+around old-Numpy/old-Numba constructor quirks rather than repeated-compilation
+overhead.
diff --git a/py_fatigue/__init__.py b/py_fatigue/__init__.py
index 131a614..9700e13 100644
--- a/py_fatigue/__init__.py
+++ b/py_fatigue/__init__.py
@@ -1,17 +1,26 @@
 # -*- coding: utf-8 -*-
 """Py-fatigue bundles the main functionality for performing cyclic stress (fatigue) analysis and cycle-counting."""  # noqa: E501  # pylint: disable=C0301
+
+from typing import TYPE_CHECKING
+
 from .cycle_count.cycle_count import CycleCount
 from .material.sn_curve import SNCurve
 from .material.crack_growth_curve import ParisCurve, WalkerCurve
 from .version import __version__
 from . import cycle_count, geometry, material, damage, styling, testing
 
+if TYPE_CHECKING:
+    from .damage import crack_growth
+    from .utils import warmup_numba
+
 __all__ = [
     "CycleCount",
     "SNCurve",
     "ParisCurve",
     "WalkerCurve",
+    "warmup_numba",
     "cycle_count",
+    "crack_growth",
     "damage",
     "geometry",
     "material",
@@ -19,3 +28,17 @@
     "testing",
     "__version__",
 ]
+
+
+def __getattr__(name: str):
+    """Lazily expose optional package-level helpers."""
+
+    if name == "warmup_numba":
+        from .utils import warmup_numba  # pylint: disable=C0415
+
+        return warmup_numba
+    if name == "crack_growth":
+        from .damage import crack_growth  # pylint: disable=C0415
+
+        return crack_growth
+    raise AttributeError(f"module {__name__!r} has no attribute {name!r}")
diff --git a/py_fatigue/cycle_count/__init__.py b/py_fatigue/cycle_count/__init__.py
index f621e79..a9ca616 100644
--- a/py_fatigue/cycle_count/__init__.py
+++ b/py_fatigue/cycle_count/__init__.py
@@ -1,7 +1,10 @@
 # -*- coding: utf-8 -*-
 """Cycle counting algorithms for fatigue analysis."""
+
 from .cycle_count import CycleCount
+from . import utils
 
 __all__ = [
     "CycleCount",
+    "utils",
 ]
diff --git a/py_fatigue/cycle_count/cycle_count.py b/py_fatigue/cycle_count/cycle_count.py
index d947d7e..0902286 100644
--- a/py_fatigue/cycle_count/cycle_count.py
+++ b/py_fatigue/cycle_count/cycle_count.py
@@ -126,9 +126,9 @@ def _build_input_data_from_json(  # noqa: C901
 ) -> dict:
     """Function returning the same data it receives"""
     data = _assess_json_keys(data)
-    the_hist = np.empty(0)
-    range_bin_centers = np.empty(0)
-    mean_bin_centers = np.empty(0)
+    the_hist: Any = np.empty(0)
+    range_bin_centers: Any = np.empty(0)
+    mean_bin_centers: Any = np.empty(0)
     if "hist" not in data:
         data["hist"] = []
     if (
@@ -2004,6 +2004,8 @@ def pbar_sum(cc_list: Sequence[CycleCount]) -> CycleCount:
     """
 
     if len(cc_list) == 1:
+        if hasattr(cc_list, "iloc"):
+            return cc_list.iloc[0]  # type: ignore[return-value]
         return cc_list[0]
 
     chunk_size = int(np.sqrt(len(cc_list)))
diff --git a/py_fatigue/cycle_count/rainflow.py b/py_fatigue/cycle_count/rainflow.py
index 3f0bdec..a229f0d 100644
--- a/py_fatigue/cycle_count/rainflow.py
+++ b/py_fatigue/cycle_count/rainflow.py
@@ -205,7 +205,7 @@ def findcross_indices(xn: Any) -> np.ndarray:
     return ind[:idx]
 
 
-@njit(int64(int64[:], int8[:]))
+@njit(int64(int64[:], int8[:]), cache=True)
 def _findcross(ind, y):
     """Return indices to zero up and downcrossings of a vector
 
diff --git a/py_fatigue/cycle_count/utils.py b/py_fatigue/cycle_count/utils.py
index 743ffd3..f7a2c1e 100644
--- a/py_fatigue/cycle_count/utils.py
+++ b/py_fatigue/cycle_count/utils.py
@@ -1,5 +1,6 @@
 # -*- coding: utf-8 -*-
 """Utilities for the cycle count module."""
+
 # pragma: no cover
 
 from __future__ import annotations
@@ -105,7 +106,8 @@ def aggregate_cc(  # pragma: no cover
         The time window to cluster the dataframe by. It must be a valid pandas
         date offset frequency string or 'all'.
         For all the frequency string aliases offered by pandas, see:
-        `pandas-timeseries.html#dateoffset-objects <shorturl.at/dgrwW>`_.
+        `pandas timeseries offsets
+        <https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#dateoffset-objects>`_.
     save_residuals : bool, optional
         If True, the residuals sequences of each aggregated CycleCount are
         saved
@@ -154,9 +156,15 @@ def aggregate_cc(  # pragma: no cover
 
     # Retrieving the low-frequency fatigue dynamics on the aggregated dataframe
     print("\33[36m4. Retrieving LFFD on aggregated \33[1mdf\33[22m.\33[0m")
-    df_agg_rr = df_agg.applymap(
-        partial(solve_lffd, rainflow_method=rainflow_method)
-    )
+    map_cycles = getattr(df_agg, "map", None)
+    if map_cycles is None:
+        df_agg_rr = df_agg.applymap(
+            partial(solve_lffd, rainflow_method=rainflow_method)
+        )
+    else:
+        df_agg_rr = map_cycles(
+            partial(solve_lffd, rainflow_method=rainflow_method)
+        )
 
     cc_cols: list[str] = [
         col for col in df_agg_rr.columns if col.startswith("CC_")
@@ -287,13 +295,23 @@ def calc_aggregated_damage(  # pragma: no cover
     cc_cols: list[str] = [col for col in df.columns if col.startswith("CC_")]
     damages = pd.DataFrame()
     for _, sn_curve in sn.items():
-        df_1 = df[cc_cols].applymap(
-            lambda x, sk=sn_curve: np.sum(
-                get_pm(cycle_count=x, sn_curve=sk)
-                if isinstance(x, CycleCount)
-                else 0
+        map_cycles = getattr(df[cc_cols], "map", None)
+        if map_cycles is None:
+            df_1 = df[cc_cols].applymap(
+                lambda x, sk=sn_curve: np.sum(
+                    get_pm(cycle_count=x, sn_curve=sk)
+                    if isinstance(x, CycleCount)
+                    else 0
+                )
+            )
+        else:
+            df_1 = map_cycles(
+                lambda x, sk=sn_curve: np.sum(
+                    get_pm(cycle_count=x, sn_curve=sk)
+                    if isinstance(x, CycleCount)
+                    else 0
+                )
             )
-        )
         df_1["sn_curve"] = f"m={sn_curve.slope}"
         damages = pd.concat([damages, df_1])
         del df_1
diff --git a/py_fatigue/damage/crack_growth.py b/py_fatigue/damage/crack_growth.py
index b26a173..234a3df 100644
--- a/py_fatigue/damage/crack_growth.py
+++ b/py_fatigue/damage/crack_growth.py
@@ -18,7 +18,6 @@
 from ..geometry.generic import AbstractCrackGeometry
 from ..geometry.cylinder import f_hol_cyl_01
 
-
 try:
     # delete the accessor to avoid warning
     del pd.DataFrame.cg  # type: ignore
@@ -538,9 +537,9 @@ class CrackGrowth:
     def __init__(self, pandas_obj):
         # self._validate(pandas_obj)
         self._obj = pandas_obj
-        self.cg_curve = None
-        self.crack_geometry = None
-        self.final_cycles = None
+        self.cg_curve = pandas_obj.attrs.get("cg_curve")
+        self.crack_geometry = pandas_obj.attrs.get("crack_geometry")
+        self.final_cycles = pandas_obj.attrs.get("final_cycles")
 
     @staticmethod
     def _validate(obj):
@@ -629,6 +628,11 @@ def calc_growth(
             geometry,
         )
 
+        self._obj.attrs["cg_curve"] = cg_curve
+        self._obj.attrs["crack_geometry"] = crack_geometry
+        self._obj.attrs["final_cycles"] = cg_.final_cycles
+        self.cg_curve = cg_curve
+        self.crack_geometry = crack_geometry
         self._obj.final_cycles = cg_.final_cycles
         self.final_cycles = cg_.final_cycles
 
diff --git a/py_fatigue/damage/stress_life.py b/py_fatigue/damage/stress_life.py
index 86c128c..aa8918b 100644
--- a/py_fatigue/damage/stress_life.py
+++ b/py_fatigue/damage/stress_life.py
@@ -27,7 +27,7 @@
     SNCurve,
 )
 from ..styling import py_fatigue_formatwarning
-from ..utils import make_axes, numba_bisect, _plot_damage_accumulation
+from ..utils import make_axes, _plot_damage_accumulation
 
 try:
     # delete the accessor to avoid warning
@@ -55,12 +55,6 @@ def _validate(obj):
         """Validate the input DataFrame. Raise an error if the input
         DataFrame does not contain the right columns.
         """
-        if {
-            "cycles_to_failure",
-        }.issubset(obj.columns):
-            e_msg = "'cycles_to_failure' already calculated"
-            raise AttributeError(e_msg)
-
         if not {"count_cycle", "mean_stress", "stress_range"}.issubset(
             obj.columns
         ):
@@ -86,6 +80,9 @@ def damage(self, sn_curve: SNCurve):
                 f"sn_curve ({sn_curve.unit}) do not match."
             )
             raise ValueError(e_msg)
+        if {"cycles_to_failure"}.issubset(self._obj.columns):
+            e_msg = "'cycles_to_failure' already calculated"
+            raise AttributeError(e_msg)
         self._obj.sn_curve = sn_curve
         self.sn_curve = sn_curve
         self._obj["cycles_to_failure"] = sn_curve.get_cycles(
@@ -1184,8 +1181,8 @@ def find_sn_curve_intersection(
     endurance: float,
     weight: float,
     res_stress: float,
-    n_min: float,
-    n_max: float,
+    n_min: float = 1e0,
+    n_max: float = 1e10,
     tol=1e-6,
 ):
     """
@@ -1236,8 +1233,7 @@ def find_sn_curve_intersection(
     ValueError
         If the bisection method fails to find a solution.
     """
-    return numba_bisect(
-        __jit_sn_curve_residuals,
+    return __jit_sn_curve_residuals(
         n_min,
         n_max,
         tol,
@@ -1260,6 +1256,12 @@ def find_sn_curve_intersection(
     # )
 
 
+def find_sn_curve_intersection_2(*args, **kwargs):
+    """Backwards-compatible alias for the current SN intersection solver."""
+
+    return find_sn_curve_intersection(*args, **kwargs)
+
+
 @nb.njit(
     fastmath=True,
     cache=True,
@@ -1497,3 +1499,19 @@ def calc_theil_sn_damage(
     return tuple((label, *tup) for tup, label in zip(hist_dict.values(),
                                                       hist_dict.keys()))
     # fmt: on
+
+
+def calc_theil_cycles_to_failure(
+    stress_range,
+    count_cycle,
+    sn_curve: SNCurve,
+    to_failure: bool = False,
+):
+    """Backwards-compatible alias for the current Theil SN helper."""
+
+    return calc_theil_sn_damage(
+        stress_range,
+        count_cycle,
+        sn_curve,
+        to_failure=to_failure,
+    )
diff --git a/py_fatigue/material/crack_growth_curve.py b/py_fatigue/material/crack_growth_curve.py
index d3be8cd..baa20f9 100755
--- a/py_fatigue/material/crack_growth_curve.py
+++ b/py_fatigue/material/crack_growth_curve.py
@@ -11,8 +11,8 @@
 
 # Standard imports
 import abc
-
 import io
+import os
 
 # import itertools
 # import warnings
@@ -39,6 +39,8 @@
     ensure_array,
 )
 
+NUMBA_VERSION_INFO = tuple(int(part) for part in nb.__version__.split(".")[:2])
+
 
 class AbstractCrackGrowthCurve(metaclass=abc.ABCMeta):
     """Abstract Paris curve. It also implements mean stress effect according
@@ -388,6 +390,8 @@ def __init__(
         intercept: Union[int, float, list, np.ndarray],
         threshold: Union[int, float] = 0,
         critical: Union[int, float] = np.inf,
+        load_ratio: Union[int, float] = 0,
+        walker_exponent: Union[int, float] = 0,
         environment: Optional[str] = None,
         curve: Optional[str] = None,
         norm: Optional[str] = None,
@@ -422,16 +426,17 @@ def __init__(
         color : str, optional
             RGBS or HEX string for color, by default None
         """
+        _ = (load_ratio, walker_exponent)
         super().__init__(
-            slope,
-            intercept,
-            threshold,
-            critical,
-            environment,
-            curve,
-            norm,
-            unit_string,
-            color,
+            slope=slope,
+            intercept=intercept,
+            threshold=threshold,
+            critical=critical,
+            environment=environment,
+            curve=curve,
+            norm=norm,
+            unit_string=unit_string,
+            color=color,
         )
         self.__threshold = threshold
         self.__critical = critical
@@ -755,7 +760,7 @@ def get_knee_sif(
         np.ndarray
             knee SIF
         """
-        knee_sif = np.empty(self.walker_intercept.size - 1, dtype=np.float64)
+        knee_sif = np.empty(self.walker_intercept.size - 1)
         if not self.linear:
             if self.walker_intercept.size > 1:
                 for i in range(self.walker_intercept.size - 1):
@@ -988,15 +993,15 @@ def __init__(  # pylint: disable=R0913, R0917
             RGBS or HEX string for color, by default None
         """
         super().__init__(
-            slope,
-            intercept,
-            threshold,
-            critical,
-            environment,
-            curve,
-            norm,
-            unit_string,
-            color,
+            slope=slope,
+            intercept=intercept,
+            threshold=threshold,
+            critical=critical,
+            environment=environment,
+            curve=curve,
+            norm=norm,
+            unit_string=unit_string,
+            color=color,
         )
 
         self.__slope, self.__intercept = _check_param_couple_types(
@@ -1057,41 +1062,44 @@ def walker_intercept(self) -> np.ndarray:
         return self.__intercept * self.walker_correction
 
 
+@nb.njit(cache=True)
+def get_array_min(values):  # pragma: no cover
+    """Return the minimum value of a one-dimensional array."""
+
+    min_value = values[0]
+    for i in range(1, values.size):
+        if np.isnan(values[i]) or values[i] < min_value:
+            min_value = values[i]
+    return min_value
+
+
 @nb.njit(
     # 'float64[::1](float64[::1], float64[::1], float64[::1])',
     fastmath=False,
+    cache=True,
     # parallel=True,
 )
-def _calc_growth_rate(
+def _calc_growth_rate_numba(
     sif, slope, intercept, threshold, critical
 ):  # pragma: no cover  # noqa: E501  # pylint: disable=C0301
     # pylint: disable=not-an-iterable
     assert intercept.size > 0 and intercept.size == slope.size
-    assert np.nanmin(sif) >= 0
+    assert get_array_min(sif) >= 0
     assert 0 <= threshold < critical <= np.inf
 
-    knees_sif = np.empty(intercept.size - 1, dtype=np.float64)
-    if intercept.size > 1:
-        for i in nb.prange(intercept.size - 1):
-            knees_sif[i] = (intercept[i] / intercept[i + 1]) ** (
-                1 / (slope[i + 1] - slope[i])
-            )
-    knees_sif = np.hstack(
-        (
-            np.array([0.9999999999 * threshold]),
-            knees_sif,
-            np.array([critical / 0.9999999999]),
-        )
-    )
     e_msg = (
-        "Knee(s) not in between threshold and critical SIF."
+        "Knee(s) are not ordered."
         + "\nCheck the definitions of slope, intercept, threshold, critical."
     )
-    assert np.all(np.diff(knees_sif) > 0), e_msg
-    # print("knees_sif:", knees_sif)
-    idx = np.digitize(sif, knees_sif, right=False) - 1
-    # print("idx:", idx)
-    the_growth_rate = np.empty(sif.size, dtype=np.float64)
+    if intercept.size > 1:
+        prev_knee = (intercept[0] / intercept[1]) ** (1 / (slope[1] - slope[0]))
+        for i in range(1, intercept.size - 1):
+            knee = (intercept[i] / intercept[i + 1]) ** (
+                1 / (slope[i + 1] - slope[i])
+            )
+            assert knee > prev_knee, e_msg
+            prev_knee = knee
+    the_growth_rate = sif.copy()
     for i in nb.prange(sif.size):
         if sif[i] < threshold:  # below threshold
             the_growth_rate[i] = 0
@@ -1099,7 +1107,14 @@ def _calc_growth_rate(
         if sif[i] > critical or sif[i] == np.inf:
             the_growth_rate[i] = np.inf
             continue
-        the_growth_rate[i] = intercept[idx[i]] * sif[i] ** slope[idx[i]]
+        idx = 0
+        for j in range(intercept.size - 1):
+            knee = (intercept[j] / intercept[j + 1]) ** (
+                1 / (slope[j + 1] - slope[j])
+            )
+            if sif[i] >= knee:
+                idx = j + 1
+        the_growth_rate[i] = intercept[idx] * sif[i] ** slope[idx]
 
     return the_growth_rate
 
@@ -1107,58 +1122,105 @@ def _calc_growth_rate(
 @nb.njit(
     # 'float64[::1](float64[::1], float64[::1], float64[::1])',
     fastmath=False,
+    cache=True,
     # parallel=True,
 )
-def _calc_sif(
+def _calc_sif_numba(
     growth_rate, slope, intercept, threshold, critical
 ):  # pragma: no cover  # noqa: E501  # pylint: disable=C0301
     # pylint: disable=not-an-iterable
     assert intercept.size > 0 and intercept.size == slope.size
-    assert np.nanmin(growth_rate) >= 0
+    assert get_array_min(growth_rate) >= 0
     assert 0 <= threshold < critical <= np.inf
 
-    knees_growth_rate = np.empty(intercept.size - 1, dtype=np.float64)
-    if intercept.size > 1:
-        for i in nb.prange(intercept.size - 1):
-            m_i = slope[i + 1] / (slope[i + 1] - slope[i])
-            m_i_p_1 = slope[i] / (slope[i + 1] - slope[i])
-            knees_growth_rate[i] = (
-                intercept[i] ** m_i / intercept[i + 1] ** m_i_p_1
-            )
-
-    knees_growth_rate = np.hstack(
-        (
-            np.array(
-                [intercept[0] * (0.9999999999999999 * threshold) ** slope[0]]
-            ),
-            knees_growth_rate,
-            np.array(
-                [intercept[-1] * (critical / 0.9999999999999999) ** slope[-1]]
-            ),
-        )
-    )
+    lower_bound = intercept[0] * (0.9999999999999999 * threshold) ** slope[0]
+    upper_bound = intercept[-1] * (critical / 0.9999999999999999) ** slope[-1]
     e_msg = (
         "Knee(s) not in between threshold and critical SIF."
         + "\nCheck the definitions of slope, intercept, threshold, critical."
     )
-    assert np.all(np.diff(knees_growth_rate) > 0), e_msg
-    # print("knees_growth_rate:", knees_growth_rate)
-    idx = np.digitize(growth_rate, knees_growth_rate, right=False) - 1
-    # print("idx:", idx)
-    the_sif = np.empty(growth_rate.size, dtype=np.float64)
-    nr_knees = intercept.size - 1
+    if intercept.size > 1:
+        prev_growth_rate = lower_bound
+        for i in range(intercept.size - 1):
+            m_i = slope[i + 1] / (slope[i + 1] - slope[i])
+            m_i_p_1 = slope[i] / (slope[i + 1] - slope[i])
+            knee_growth_rate = intercept[i] ** m_i / intercept[i + 1] ** m_i_p_1
+            assert knee_growth_rate > prev_growth_rate, e_msg
+            prev_growth_rate = knee_growth_rate
+        assert upper_bound > prev_growth_rate, e_msg
+    else:
+        assert upper_bound > lower_bound, e_msg
+
+    the_sif = growth_rate.copy()
     for i in nb.prange(growth_rate.size):
-        if idx[i] <= 0:  # below threshold
+        if growth_rate[i] < lower_bound:  # below threshold
             the_sif[i] = threshold
             continue
-        if idx[i] > nr_knees:  # above threshold
+        if growth_rate[i] >= upper_bound:  # above threshold
             the_sif[i] = critical
             continue
-        the_sif[i] = (growth_rate[i] / intercept[idx[i]]) ** (1 / slope[idx[i]])
+        idx = 0
+        for j in range(intercept.size - 1):
+            m_i = slope[j + 1] / (slope[j + 1] - slope[j])
+            m_i_p_1 = slope[j] / (slope[j + 1] - slope[j])
+            knee_growth_rate = intercept[j] ** m_i / intercept[j + 1] ** m_i_p_1
+            if growth_rate[i] >= knee_growth_rate:
+                idx = j + 1
+        the_sif[i] = (growth_rate[i] / intercept[idx]) ** (1 / slope[idx])
 
     return the_sif
 
 
+def use_python_crack_growth_kernels() -> bool:
+    """Return whether crack-growth kernels should bypass Numba dispatch."""
+
+    return (
+        NUMBA_VERSION_INFO <= (0, 61)
+        or getattr(nb.config, "DISABLE_JIT", False)
+        or os.environ.get("NUMBA_DISABLE_JIT") == "1"
+    )
+
+
+def _calc_growth_rate(
+    sif, slope, intercept, threshold, critical
+):  # pragma: no cover  # noqa: E501  # pylint: disable=C0301
+    if use_python_crack_growth_kernels():
+        return _calc_growth_rate_numba.py_func(
+            sif,
+            slope,
+            intercept,
+            threshold,
+            critical,
+        )
+    return _calc_growth_rate_numba(
+        sif,
+        slope,
+        intercept,
+        threshold,
+        critical,
+    )
+
+
+def _calc_sif(
+    growth_rate, slope, intercept, threshold, critical
+):  # pragma: no cover  # noqa: E501  # pylint: disable=C0301
+    if use_python_crack_growth_kernels():
+        return _calc_sif_numba.py_func(
+            growth_rate,
+            slope,
+            intercept,
+            threshold,
+            critical,
+        )
+    return _calc_sif_numba(
+        growth_rate,
+        slope,
+        intercept,
+        threshold,
+        critical,
+    )
+
+
 def _paris_curve_data_points(pc: ParisCurve) -> tuple:
     """Draw Paris' curve data points, including possible knees and
     threshold/critical SIF.
@@ -1192,7 +1254,12 @@ def _paris_curve_data_points(pc: ParisCurve) -> tuple:
     )  # 0.999999999999 to avoid rounding errors
 
     if not pc.linear:
-        sif_plot = np.sort(np.append(sif_plot, np.array(pc.get_knee_sif())))
+        knee_sif = np.array(pc.get_knee_sif())
+        knee_sif = knee_sif[
+            (knee_sif > pc.threshold)
+            & (knee_sif < 0.999999999999 * pc.critical)
+        ]
+        sif_plot = np.sort(np.append(sif_plot, knee_sif))
     growth_rate_plot = pc.get_growth_rate(sif_plot)
     if pc.threshold > 0:
         growth_rate_plot = np.hstack(
diff --git a/py_fatigue/material/sn_curve.py b/py_fatigue/material/sn_curve.py
index 2ce2b30..45aa582 100644
--- a/py_fatigue/material/sn_curve.py
+++ b/py_fatigue/material/sn_curve.py
@@ -18,6 +18,7 @@
 import abc
 import io
 import itertools
+import os
 import warnings
 
 # Non-standard "from" imports
@@ -689,7 +690,13 @@ def get_cycles(
         #     axis=1
         # ), self.endurance)
 
-        return _calc_cycles_2(
+        calc_cycles = _calc_cycles_2
+        if (
+            getattr(nb.config, "DISABLE_JIT", False)
+            or os.environ.get("NUMBA_DISABLE_JIT") == "1"
+        ):
+            calc_cycles = _calc_cycles_2.py_func
+        return calc_cycles(
             stress_range, self.slope, self.intercept, self.endurance
         )
 
@@ -711,9 +718,13 @@ def get_stress(
         #         cycles.shape[0],-1),
         #         axis=1
         # ), endurance_stress)
-        return _calc_stress_2(
-            cycles, self.slope, self.intercept, self.endurance
-        )
+        calc_stress = _calc_stress_2
+        if (
+            getattr(nb.config, "DISABLE_JIT", False)
+            or os.environ.get("NUMBA_DISABLE_JIT") == "1"
+        ):
+            calc_stress = _calc_stress_2.py_func
+        return calc_stress(cycles, self.slope, self.intercept, self.endurance)
 
     def n(  # pylint: disable=invalid-name
         self, sigma: int | float | list | np.ndarray
@@ -941,17 +952,29 @@ def plot(
         return fig, ax
 
 
+@nb.njit(cache=True)
+def get_sn_array_min(values):  # pragma: no cover
+    """Return the minimum value of a one-dimensional array."""
+
+    min_value = values[0]
+    for i in range(1, values.size):
+        if np.isnan(values[i]) or values[i] < min_value:
+            min_value = values[i]
+    return min_value
+
+
 @nb.njit(
     # 'float64[::1](float64[::1], float64[::1], float64[::1])',
     fastmath=False,
+    cache=True,
     # parallel=True,
 )
 def _calc_cycles(stress, slope, intercept, endurance):  # pragma: no cover
     # pylint: disable=not-an-iterable
     assert intercept.size > 0 and intercept.size == slope.size
-    assert np.min(stress) >= 0
+    assert get_sn_array_min(stress) >= 0
     log10 = np.log(10)
-    the_cycles = np.empty(stress.size, dtype=np.float64)
+    the_cycles = np.empty(stress.size)
     for i in nb.prange(stress.size):
         log_stress = np.log(stress[i])
         max_i = intercept[0] * log10 - slope[0] * log_stress
@@ -967,6 +990,7 @@ def _calc_cycles(stress, slope, intercept, endurance):  # pragma: no cover
 @nb.njit(
     # 'float64[::1](float64[::1], float64[::1], float64[::1])',
     fastmath=False,
+    cache=True,
     # parallel=True,
 )
 def _calc_cycles_2(stress, slope, intercept, endurance):  # pragma: no cover
@@ -983,11 +1007,11 @@ def _calc_cycles_2(stress, slope, intercept, endurance):  # pragma: no cover
     numpy.ndarray: Array of calculated cycles to failure.
     """
     assert intercept.size > 0 and intercept.size == slope.size
-    assert np.min(stress) >= 0
+    assert get_sn_array_min(stress) >= 0
 
     log_stress = np.log10(stress)
     log_endurance = np.log10(endurance)
-    log_knee_stress = np.empty(intercept.size - 1, dtype=np.float64)
+    log_knee_stress = np.empty(intercept.size - 1)
     if intercept.size > 1:
         for i in nb.prange(intercept.size - 1):  # pylint: disable=E1133
             log_knee_stress[i] = (intercept[i + 1] - intercept[i]) / (
@@ -1002,7 +1026,7 @@ def _calc_cycles_2(stress, slope, intercept, endurance):  # pragma: no cover
         )
     )
     idx = np.digitize(log_stress, log_knee_stress, right=False) - 1
-    the_cycles = np.empty(stress.size, dtype=np.float64)
+    the_cycles = np.empty(stress.size)
     nr_knees = intercept.size - 1
     for i in nb.prange(stress.size):  # pylint: disable=E1133
         if idx[i] <= 0:
@@ -1021,14 +1045,15 @@ def _calc_cycles_2(stress, slope, intercept, endurance):  # pragma: no cover
 @nb.njit(
     # 'float64[::1](float64[::1], float64[::1], float64[::1])',
     fastmath=False,
+    cache=True,
     # parallel=True,
 )
 def _calc_stress(cycles, slope, intercept, endurance):  # pragma: no cover
     # pylint: disable=not-an-iterable
     assert intercept.size > 0 and intercept.size == slope.size
-    assert np.min(cycles) > 0
+    assert get_sn_array_min(cycles) > 0
     log10 = np.log(10)
-    the_stress = np.empty(cycles.size, dtype=np.float64)
+    the_stress = np.empty(cycles.size)
     for i in nb.prange(cycles.size):
         log_cycles = np.log(cycles[i])
         max_i = (intercept[0] * log10 - log_cycles) / slope[0]
@@ -1064,12 +1089,12 @@ def _calc_stress_2(cycles, slope, intercept, endurance):  # pragma: no cover
     numpy.ndarray: Array of calculated cycles to failure.
     """
     assert intercept.size > 0 and intercept.size == slope.size
-    assert np.nanmin(cycles) >= 0
+    assert get_sn_array_min(cycles) >= 0
 
     log_cycles = np.log10(cycles)
     log_endurance = np.log10(endurance)
-    log_knee_stress = np.empty(intercept.size - 1, dtype=np.float64)
-    log_knee_cycles = np.empty(intercept.size - 1, dtype=np.float64)
+    log_knee_stress = np.empty(intercept.size - 1)
+    log_knee_cycles = np.empty(intercept.size - 1)
     if intercept.size > 1:
         for i in nb.prange(intercept.size - 1):  # pylint: disable=E1133
             log_knee_stress[i] = (intercept[i + 1] - intercept[i]) / (
@@ -1089,7 +1114,7 @@ def _calc_stress_2(cycles, slope, intercept, endurance):  # pragma: no cover
     )
     endurance_stress = 10 ** log_knee_stress[-1]
     idx = np.digitize(log_cycles, log_knee_cycles, right=False) - 1
-    the_stress = np.empty(cycles.size, dtype=np.float64)
+    the_stress = np.empty(cycles.size)
     nr_knees = intercept.size - 1
     for i in nb.prange(cycles.size):  # pylint: disable=E1133
         if idx[i] <= 0:
@@ -1156,3 +1181,7 @@ def __sn_curve_residuals(  # pragma: no cover
 
 # Create a jitted bisection function specialized for root_func
 __jit_sn_curve_residuals = compile_specialized_bisect(__sn_curve_residuals)
+
+# Backwards-compatible aliases used by older notebooks.
+sn_curve_residuals = __sn_curve_residuals
+jit_sn_curve_residuals = __jit_sn_curve_residuals
diff --git a/py_fatigue/mean_stress/corrections.py b/py_fatigue/mean_stress/corrections.py
index 88e364c..d796d37 100644
--- a/py_fatigue/mean_stress/corrections.py
+++ b/py_fatigue/mean_stress/corrections.py
@@ -5,6 +5,7 @@
 """
 
 from __future__ import annotations
+from typing import cast
 import logging
 import warnings
 
@@ -26,7 +27,7 @@ def dnvgl_mean_stress_correction(
     plot: bool = False,
 ) -> np.ndarray:
     """Calculates the mean stress correction according to par.
-    2.5 of`DNVGL-RP-C203 <https://bit.ly/3dUZ1OY>`_ which includes
+    2.5 of DNVGL-RP-C203, which includes
     an attenuation factor :math:`p` for the stress ranges if the
     following cases:
 
@@ -320,8 +321,7 @@ def goodman_haigh_mean_stress_correction(  # pylint: disable=R0912 # noqa: C901,
     See Also
     --------
     :func:`py_fatigue.utils.numba_newton`,
-    :func:`py_fatigue.utils.compile_specialized_newton`,
-    :func:`py_fatigue.mean_stress.corrections.__goodman
+    :func:`py_fatigue.utils.compile_specialized_newton`
 
     Raises
     ------
@@ -394,6 +394,8 @@ def goodman_haigh_mean_stress_correction(  # pylint: disable=R0912 # noqa: C901,
     for r_out_val in r_out:
         # NOTE: Special cases for r_out = -1 and r_out = 0
         # NOTE: These cases are solved analytically to improve performance
+        # NOTE: The case correction_exponent = 1 is also solved
+        #       analytically, as the equation becomes linear in amp_out.
         if r_out_val == -1:
             # amp_out.append(
             #     amp_in / (1 - (mean_in / ult_s) ** correction_exponent)
@@ -404,6 +406,17 @@ def goodman_haigh_mean_stress_correction(  # pylint: disable=R0912 # noqa: C901,
             )
             mean_out[r_out == -1, :] = np.zeros_like(amp_in)
             continue
+        if correction_exponent == 1:
+            r_in_term = (1 + r_in) / (1 - r_in) * amp_in / ult_s
+            r_out_term = (1 + r_out_val) / (1 - r_out_val) / ult_s
+            amp_out_fsolve = ((1 - r_in_term) / amp_in + r_out_term) ** -1
+            amp_out[r_out == r_out_val, :] = np.clip(amp_out_fsolve, 0.0, ult_s)
+            mean_out[r_out == r_out_val, :] = (
+                amp_out[r_out == r_out_val, :]
+                * (1 + r_out_val)
+                / (1 - r_out_val)
+            )
+            continue
         amp_out_fsolve = []
         for i in range(len(initial_guess)):
             # Solve the implicit equation using fsolve
@@ -425,7 +438,7 @@ def goodman_haigh_mean_stress_correction(  # pylint: disable=R0912 # noqa: C901,
                                    f"{r_out_val}")
                 amp_out_fsolve.append(np.nan)
             else:
-                amp_out_fsolve.append(sol if sol < ult_s else ult_s)
+                amp_out_fsolve.append(np.clip(sol, 0.0, ult_s))
             # fmt: on
         # amp_out.append(amp_out_fsolve)
         # mean_out.append(
@@ -444,8 +457,8 @@ def goodman_haigh_mean_stress_correction(  # pylint: disable=R0912 # noqa: C901,
     else:
         srt_idx = np.argsort(r_out)
         srt_r_out = r_out[srt_idx]
-        srt_amp_out = amp_out[srt_idx, :]
-        srt_mean_out = mean_out[srt_idx, :]
+        srt_amp_out = cast(np.ndarray, amp_out[srt_idx, :])
+        srt_mean_out = cast(np.ndarray, mean_out[srt_idx, :])
     if plot:
         # Create figure with two subplots
         fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
diff --git a/py_fatigue/testing.py b/py_fatigue/testing.py
index 870eb6f..b996c72 100644
--- a/py_fatigue/testing.py
+++ b/py_fatigue/testing.py
@@ -1,5 +1,5 @@
 from __future__ import annotations
-from typing import Any
+from typing import Any, cast
 import warnings
 
 import numpy as np
@@ -46,7 +46,10 @@ def get_sampled_time(
         warnings.formatwarning = py_fatigue_formatwarning
         warnings.warn(w_msg, UserWarning)
 
-    return np.arange(start, start + duration, 1 / fs)
+    return cast(
+        np.ndarray[Any, np.dtype[np.float64]],
+        np.arange(start, start + duration, 1 / fs),
+    )
 
 
 def get_random_data(
diff --git a/py_fatigue/utils.py b/py_fatigue/utils.py
index 510ee89..b87a838 100644
--- a/py_fatigue/utils.py
+++ b/py_fatigue/utils.py
@@ -6,8 +6,11 @@
 
 # Packages from the Python Standard Library
 from __future__ import annotations
+from contextlib import redirect_stdout
 from dataclasses import dataclass
-from functools import wraps
+from functools import lru_cache, wraps
+from importlib import import_module
+from io import StringIO
 from types import FunctionType
 from typing import (
     Any,
@@ -25,6 +28,7 @@
 )
 import copy
 import logging
+import os
 
 # Packages from non-standard libraries
 # from pydantic.fields import ModelField
@@ -32,6 +36,9 @@
 import matplotlib.pyplot as plt
 import numba as nb
 import numpy as np
+from numba.extending import overload
+
+NUMBA_SPECIALIZED_ATTR = "__py_fatigue_numba_specialized__"
 
 
 # Decorator
@@ -613,7 +620,7 @@ def validate(cls, val, field):
 
 def to_numba_dict(
     data: dict, key_type: type = str, val_type: type = float
-) -> nb.types.DictType:
+) -> Any:
     """Converts a dictionary to a numba typed dict, provided the output
     key and value types.
 
@@ -631,13 +638,22 @@ def to_numba_dict(
     nb.types.DictType
         The numba typed dictionary.
     """
+    filtered_items = {
+        key: value
+        for key, value in data.items()
+        if isinstance(value, val_type) and isinstance(key, key_type)
+    }
+    if (
+        getattr(nb.config, "DISABLE_JIT", False)
+        or os.environ.get("NUMBA_DISABLE_JIT") == "1"
+    ):
+        return filtered_items
+
     dct = nb.typed.Dict.empty(
-        key_type=nb.types.string,
-        value_type=nb.float64,
+        key_type=nb.types.unicode_type,
+        value_type=nb.types.float64,
     )
-    for key, value in data.items():
-        if not isinstance(value, val_type) or not isinstance(key, key_type):
-            continue
+    for key, value in filtered_items.items():
         dct[key] = value
     return dct
 
@@ -735,18 +751,88 @@ def calc_slope_intercept(
 #     return nb.njit()(python_bisect)
 
 
+def scalarize_numba_result(value):
+    """Normalize scalar-like numba residual returns to a float."""
+
+    return value
+
+
+# pylint: disable-next=unreachable
+@nb.njit(nb.float64())
+def raise_non_scalar_numba_result():
+    """Raise the scalar residual error from nopython code."""
+
+    raise TypeError("Compiled function must return a scalar value")
+    return 0.0  # pylint: disable=unreachable
+
+
+@overload(scalarize_numba_result)
+# pylint: disable-next=too-many-return-statements
+def scalarize_numba_result_overload(value):
+    """Compile scalar-like residual return normalization."""
+
+    if isinstance(value, nb.types.Number):
+
+        def scalar_impl(value):
+            return float(value)
+
+        return scalar_impl
+    if isinstance(value, nb.types.UniTuple) and value.count == 1:
+
+        def unituple_scalar_impl(value):
+            return float(value[0])
+
+        return unituple_scalar_impl
+    if isinstance(value, nb.types.UniTuple):
+
+        def unituple_reject_impl(value):
+            if len(value) != 1:
+                return raise_non_scalar_numba_result()
+            return float(value[0])
+
+        return unituple_reject_impl
+    if isinstance(value, nb.types.Tuple) and len(value) == 1:
+
+        def tuple_scalar_impl(value):
+            return float(value[0])
+
+        return tuple_scalar_impl
+    if isinstance(value, nb.types.Tuple):
+
+        def tuple_reject_impl(value):
+            if len(value) != 1:
+                return raise_non_scalar_numba_result()
+            return float(value[0])
+
+        return tuple_reject_impl
+    if isinstance(value, nb.types.Array):
+
+        def array_impl(value):
+            if value.size != 1:
+                raise TypeError("Compiled function must return a scalar value")
+            return float(value.flat[0])
+
+        return array_impl
+    if isinstance(value, (nb.types.List, nb.types.ListType)):
+
+        def list_impl(value):
+            if len(value) != 1:
+                raise TypeError("Compiled function must return a scalar value")
+            return float(value[0])
+
+        return list_impl
+    return None
+
+
+@lru_cache(maxsize=None)
 def compile_specialized_bisect(fun):
     """
     Returns a compiled bisection implementation for `f`.
     """
     compiled_f = nb.njit()(fun)
 
-    def python_bisect(a, b, tol, mxiter, *args):
+    def python_bisect_compat(a, b, tol, mxiter, *args):
         its = 0
-        len_args = len(args)
-
-        if len_args > 3:
-            raise ValueError("Too many extra arguments for compiled function")
 
         def ensure_scalar(value):
             if isinstance(value, (tuple, list)):
@@ -766,15 +852,7 @@ def ensure_scalar(value):
             return value
 
         def evaluate(point):
-            if len_args == 0:
-                result = compiled_f(point)
-            elif len_args == 1:
-                result = compiled_f(point, args[0])
-            elif len_args == 2:
-                result = compiled_f(point, args[0], args[1])
-            else:
-                result = compiled_f(point, args[0], args[1], args[2])
-            return ensure_scalar(result)
+            return ensure_scalar(compiled_f(point, *args))
 
         fa = evaluate(a)
         fb = evaluate(b)
@@ -799,9 +877,44 @@ def evaluate(point):
             fc = evaluate(c)
         return c
 
-    return python_bisect
+    if (
+        getattr(nb.config, "DISABLE_JIT", False)
+        or os.environ.get("NUMBA_DISABLE_JIT") == "1"
+    ):
+        setattr(python_bisect_compat, NUMBA_SPECIALIZED_ATTR, True)
+        return python_bisect_compat
+
+    def python_bisect(a, b, tol, mxiter, *args):
+        its = 0
+        left = float(a)
+        right = float(b)
+        fa = scalarize_numba_result(compiled_f(left, *args))
+        fb = scalarize_numba_result(compiled_f(right, *args))
+
+        if abs(fa) < tol:
+            return left
+        if abs(fb) < tol:
+            return right
+
+        c = (left + right) / 2.0
+        fc = scalarize_numba_result(compiled_f(c, *args))
+
+        while abs(fc) > tol and its < mxiter:
+            its += 1
+            if fa * fc < 0:
+                right = c
+                fb = fc
+            else:
+                left = c
+                fa = fc
+            c = (left + right) / 2.0
+            fc = scalarize_numba_result(compiled_f(c, *args))
+        return c
+
+    return nb.njit()(python_bisect)
 
 
+@lru_cache(maxsize=None)
 def compile_specialized_newton(fun):
     """
     Returns a compiled Newton–Raphson implementation for f that accepts extra
@@ -882,7 +995,9 @@ def numba_bisect(fun, a, b, tol, mxiter, *args):
     """
     A wrapper that compiles `f` if it is a regular Python function.
     """
-    if isinstance(fun, FunctionType):
+    if isinstance(fun, FunctionType) and not getattr(
+        fun, NUMBA_SPECIALIZED_ATTR, False
+    ):
         jit_bisect_func = compile_specialized_bisect(fun)
         return jit_bisect_func(a, b, tol, mxiter, *args)
     return fun(a, b, tol, mxiter, *args)
@@ -893,12 +1008,128 @@ def numba_newton(fun, x0, tol, mxiter, *args):
     A wrapper that compiles f if it is a regular Python function and calls the
     Newton–Raphson routine with extra arguments.
     """
-    if isinstance(fun, FunctionType):
+    if isinstance(fun, FunctionType) and not getattr(
+        fun, NUMBA_SPECIALIZED_ATTR, False
+    ):
         jit_newton_func = compile_specialized_newton(fun)
         return jit_newton_func(x0, tol, mxiter, *args)
     return fun(x0, tol, mxiter, *args)
 
 
+def warmup_numba() -> None:
+    """Compile the common numba dispatchers for the current Python process.
+
+    This function is intentionally opt-in. Calling it can take noticeable time,
+    but it moves first-call compilation cost to application startup. Functions
+    decorated with ``cache=True`` can also reuse numba's disk cache in later
+    processes; jitclasses and dynamically specialized root finders still need an
+    in-process warmup.
+    """
+
+    if (
+        getattr(nb.config, "DISABLE_JIT", False)
+        or os.environ.get("NUMBA_DISABLE_JIT") == "1"
+    ):
+        return
+
+    rainflow_module = import_module(
+        ".cycle_count.rainflow", package=__package__
+    )
+    crack_growth_module = import_module(
+        ".damage.crack_growth", package=__package__
+    )
+    stress_life_module = import_module(
+        ".damage.stress_life", package=__package__
+    )
+    cylinder_module = import_module(".geometry.cylinder", package=__package__)
+    crack_growth_curve_module = import_module(
+        ".material.crack_growth_curve", package=__package__
+    )
+    sn_curve_module = import_module(".material.sn_curve", package=__package__)
+    mean_stress_module = import_module(
+        ".mean_stress.corrections", package=__package__
+    )
+
+    findcross = rainflow_module.findcross
+    findtp = rainflow_module.findtp
+    calc_crack_growth = crack_growth_module.CalcCrackGrowth
+    get_sif = crack_growth_module.get_sif
+    calc_theil_sn_damage = stress_life_module.calc_theil_sn_damage
+    find_sn_curve_intersection = stress_life_module.find_sn_curve_intersection
+    f_hol_cyl_01 = cylinder_module.f_hol_cyl_01
+    paris_curve_cls = crack_growth_curve_module.ParisCurve
+    sn_curve_cls = sn_curve_module.SNCurve
+    goodman_correction = mean_stress_module.goodman_haigh_mean_stress_correction
+
+    sn_curve = sn_curve_cls([3.0, 5.0], [12.0, 15.0], endurance=1e9)
+    stress_range = np.ascontiguousarray([90.0, 120.0, 180.0], dtype=np.float64)
+    cycles = np.ascontiguousarray([1e5, 1e6, 1e7], dtype=np.float64)
+    sn_curve.get_cycles(stress_range)
+    sn_curve.get_stress(cycles)
+
+    paris_curve = paris_curve_cls(
+        slope=[2.88, 5.1],
+        intercept=[1e-16, 1e-20],
+        threshold=20.0,
+        critical=2_000.0,
+    )
+    sif_range = np.ascontiguousarray([20.0, 100.0, 1_000.0], dtype=np.float64)
+    growth_rate = np.ascontiguousarray([1e-10, 1e-8, 1e-6], dtype=np.float64)
+    paris_curve.get_growth_rate(sif_range)
+    paris_curve.get_sif(growth_rate)
+
+    signal = np.ascontiguousarray([0.0, 1.0, -1.0, 2.0, -0.5], dtype=np.float64)
+    findcross(signal)
+    findtp(signal)
+
+    geometry = to_numba_dict(
+        {
+            "initial_depth": 1.0,
+            "outer_diameter": 100.0,
+            "thickness": 10.0,
+            "height": 200.0,
+            "width_to_depth_ratio": 2.0,
+        }
+    )
+    f_hol_cyl_01(1.0, geometry)
+    get_sif(100.0, 1.0, "HOL_CYL_01", geometry)
+
+    inf_geometry = to_numba_dict({"initial_depth": 1.0})
+    with redirect_stdout(StringIO()):
+        calc_crack_growth(
+            stress_range,
+            np.ones_like(stress_range),
+            np.ascontiguousarray([3.0], dtype=np.float64),
+            np.ascontiguousarray([1e-12], dtype=np.float64),
+            0.0,
+            1e9,
+            "INF_SUR_00",
+            inf_geometry,
+        )
+
+    calc_theil_sn_damage(stress_range, np.ones_like(stress_range), sn_curve)
+    find_sn_curve_intersection(
+        sn_curve.slope,
+        sn_curve.intercept,
+        sn_curve.endurance,
+        100.0,
+        0.0,
+        1e3,
+        1e9,
+    )
+
+    amp_in = np.ascontiguousarray([50.0, 100.0], dtype=np.float64)
+    mean_in = np.ascontiguousarray([0.0, 20.0], dtype=np.float64)
+    r_out = np.ascontiguousarray([-1.0, 0.0], dtype=np.float64)
+    goodman_correction(
+        amp_in,
+        mean_in,
+        r_out,
+        1_000.0,
+        3.0,
+    )
+
+
 def _plot_damage_accumulation(  # pragma: no cover
     cumsum_nl_dmg: np.ndarray,
     cumsum_pm_dmg: np.ndarray,
@@ -952,14 +1183,16 @@ class CustomFormatter(logging.Formatter):
     italic = "\033[3m"
     reset = "\033[0m"
     level = "\033[1m%(levelname)-8s → \033[22m"
+    lineofile = "\033[3m%(filename)s:%(lineno)d\033[0m - "
+    newline = "\n"
     message = "%(message)s"
 
     FORMATS = {
         logging.DEBUG: grey + "🐞 " + level + italic + message + reset,
         logging.INFO: blue + "ℹ️ " + level + italic + message + reset,
-        logging.WARNING: yellow + "⚠️ " + level + italic + message + reset,
-        logging.ERROR: red + "⛔ " + level + italic + message + reset,
-        logging.CRITICAL: red + "🆘 " + level + bold + italic + message + reset,
+        logging.WARNING: yellow + "⚠️ " + level + newline + lineofile + italic + newline +message + reset,
+        logging.ERROR: red + "⛔ " + level + newline + lineofile + italic + newline +message + reset,
+        logging.CRITICAL: red + "🆘 " + level + newline + lineofile + bold + italic + newline +message + reset,
     }
 
     def format(self, record):
diff --git a/py_fatigue/version.py b/py_fatigue/version.py
index d59e028..8fa452e 100644
--- a/py_fatigue/version.py
+++ b/py_fatigue/version.py
@@ -14,6 +14,7 @@
     __version__ (str): Current version of py_fatigue package.
 
 """
+
 import re
 from typing import Optional, NamedTuple
 
diff --git a/pyproject.toml b/pyproject.toml
index b1b124b..34ee59a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,16 +23,18 @@ classifiers=["Development Status :: 5 - Production/Stable",
              "Programming Language :: Python :: 3.11",
              "Programming Language :: Python :: 3.12",
              "Programming Language :: Python :: 3.13",
+             "Programming Language :: Python :: 3.14",
              "Topic :: Database"]
 readme = "README.md"
 license = {file = "LICENSE"}
-requires-python = ">=3.10,<3.14"
+requires-python = ">=3.10,<3.15"
 
 dependencies = [
     "numpy>=1.24",
     "plotly",
     "pandas>2.2",
-    "numba>=0.61,<0.62",
+    "numba>=0.61,<0.66; python_version < '3.14'",
+    "numba>=0.65,<0.66; python_version >= '3.14'",
     "matplotlib",
     "pydantic<3.0.0",
 ]
@@ -52,6 +54,7 @@ dev = [
     "hypothesis>=6.14.0",
     "notebook>=6.4",
     "rich>=13.9.4",
+    "rich-argparse>=1.8.0",
     "jsbeautifier>=1.15.4",
 ]
 docs = [
diff --git a/rst_docs/api/material/crack_growth_curve.rst b/rst_docs/api/material/crack_growth_curve.rst
index 1c29023..09c4538 100644
--- a/rst_docs/api/material/crack_growth_curve.rst
+++ b/rst_docs/api/material/crack_growth_curve.rst
@@ -1,18 +1,16 @@
 crack_growth_curve
 ==================
 
-The crack_growth_curve module contains the ParisCurve class
-and all the funcitons strictly related to its correct functioning.
-The function relies on the AbstractCrackGrowthCurve class, which can be used
-as a abstract class for the definition of new growth models such
-as NASGRO law.
-
-The ParisCurve class
-********************
+The :mod:`py_fatigue.material.crack_growth_curve` module contains the
+:class:`py_fatigue.ParisCurve` and :class:`py_fatigue.WalkerCurve` classes,
+together with the helpers used to evaluate crack-growth curves.
 
 .. autoclass:: py_fatigue.ParisCurve
     :members:
     :private-members:
     :special-members:
 
-    .. autoclasstoc::
+.. autoclass:: py_fatigue.WalkerCurve
+    :members:
+    :private-members:
+    :special-members:
diff --git a/rst_docs/api/material/sn_curve.rst b/rst_docs/api/material/sn_curve.rst
index 8ccf1e3..0a54c5b 100644
--- a/rst_docs/api/material/sn_curve.rst
+++ b/rst_docs/api/material/sn_curve.rst
@@ -1,18 +1,10 @@
 sn_curve
 ===================
 
-The sn_curve module contains the SNCurve class
-and all the funcitons strictly related to its correct functioning.
-The function relies on the AbstractSNCurve class, which can be used
-as a abstract class for the definition of new SN curve models such
-as smart SN curves.
-
-The SNCurve class
-*****************
+The :mod:`py_fatigue.material.sn_curve` module contains the
+:class:`py_fatigue.SNCurve` class and the helpers used to evaluate S-N curves.
 
 .. autoclass:: py_fatigue.SNCurve
     :members:
     :private-members:
     :special-members:
-
-    .. autoclasstoc::
diff --git a/rst_docs/user/01-absolute-noob.rst b/rst_docs/user/01-absolute-noob.rst
index e48582d..11337b9 100644
--- a/rst_docs/user/01-absolute-noob.rst
+++ b/rst_docs/user/01-absolute-noob.rst
@@ -241,11 +241,11 @@ Cycle-count matrix
   import matplotlib.pyplot as plt
   fig, axs = plt.subplots(1, 2, figsize=(12, 4.5))
   cycle_count.plot_histogram(fig=fig, ax=axs[0], plot_type="mean-range",
-                             marker="d", s=2, cmap=mpl.cm.get_cmap("coolwarm"))
+                             marker="d", s=2, cmap=plt.get_cmap("coolwarm"))
   axs[0].set_title("Cycle-count from signal")
   cycle_count_d.plot_histogram(fig=fig, ax=axs[1], plot_type="mean-range",
                                marker="s", s=10, edgecolors="#222",
-                               cmap=mpl.cm.get_cmap("coolwarm"), linewidth=0.25)
+                               cmap=plt.get_cmap("coolwarm"), linewidth=0.25)
   axs[1].set_title("Cycle-count from matrix")
   plt.show()
 
diff --git a/rst_docs/user/additional/glossary.rst b/rst_docs/user/additional/glossary.rst
index 927b4d0..1d2e245 100644
--- a/rst_docs/user/additional/glossary.rst
+++ b/rst_docs/user/additional/glossary.rst
@@ -273,7 +273,7 @@ Glossary
       See :term:`residuals<Residuals>`.
 
    SHM
-      `Stuctural health monitoring <en.wikipedia.org/wiki/Structural_health_monitoring>`_.
+      `Structural health monitoring <https://en.wikipedia.org/wiki/Structural_health_monitoring>`_.
 
    Timestamp
       A timestamp is a time information that is associated with a particular
@@ -310,8 +310,7 @@ Glossary
       we also have to deal with gaps in the data. Even if the ambition were to 
       concatenate everything, we would still make errors as such.
 
-      `Marsh *et al.* <https://www.sciencedirect.com/science/article/pii/S0142112315003333>`_
-      have introduced a brilliant approach that significantly 
+      Marsh *et al.* have introduced a brilliant approach that significantly
       reduces calculation time without losing accuracy in the final spectrum 
       histogram, as it can retrieve all the hysteresis cycles caused by LFFD without 
       needing the a-priori signal concatenation. We have modified the approach to 
diff --git a/rst_docs/user/examples/02-sn_curve.rst b/rst_docs/user/examples/02-sn_curve.rst
index 0946370..ba3d6ee 100644
--- a/rst_docs/user/examples/02-sn_curve.rst
+++ b/rst_docs/user/examples/02-sn_curve.rst
@@ -59,10 +59,10 @@ a. Multiple SN curves
 ---------------------
 
 .. note::
-    In this example we define four SN curves for free corrosion as per 
-    `DNVGL-RP-C203 <shorturl.at/ipBKL>`_ and plot them using
-    `matplotlib <https://matplotlib.org>`_ and `plotly <https://plotly.com>`_. 
-    We additionally define a random gaussian stress range-cycles history to 
+    In this example we define four SN curves for free corrosion as per
+    DNVGL-RP-C203 and plot them using `matplotlib <https://matplotlib.org>`_
+    and `plotly <https://plotly.com>`_. We additionally define a random
+    Gaussian stress range-cycles history to
     plot against the SN curves defined.
 
 .. code-block:: python
diff --git a/rst_docs/user/examples/03-cg_curve.rst b/rst_docs/user/examples/03-cg_curve.rst
index 2b73eda..2c51169 100644
--- a/rst_docs/user/examples/03-cg_curve.rst
+++ b/rst_docs/user/examples/03-cg_curve.rst
@@ -114,11 +114,10 @@ a. Definition of multiple Paris' laws
 -------------------------------------
 
 .. note::
-    In this example we define four SN curves for free corrosion as per 
-    `DNVGL-RP-C203 <shorturl.at/ipBKL>`_ and plot them using
-    `matplotlib <https://matplotlib.org>`_ and `plotly <https://plotly.com>`_. 
-    We additionally define a random gaussian stress range-cycles history to 
-    plot against the SN curves defined.
+    In this example we define four Paris curves for free corrosion as per
+    DNVGL-RP-C203 and plot them using `matplotlib <https://matplotlib.org>`_
+    and `plotly <https://plotly.com>`_. We additionally define a random
+    Gaussian stress range-cycles history to plot against the curves defined.
 
 .. code-block:: python
     :linenos:
diff --git a/rst_docs/user/examples/04-cycle_counting.rst b/rst_docs/user/examples/04-cycle_counting.rst
index c9ab3eb..5065a63 100644
--- a/rst_docs/user/examples/04-cycle_counting.rst
+++ b/rst_docs/user/examples/04-cycle_counting.rst
@@ -261,7 +261,7 @@ the cumulative (or non-cumulative) rainflow matric both for `cc` and `cc_dct`:
     cc.plot_histogram(fig=fig, ax=axs[1][0], plot_type="counts-range-cumsum", s=30)
     cc_dct.plot_histogram(fig=fig, ax=axs[1][1], plot_type="counts-range",
                           marker='s', s=20,
-                          cmap=matplotlib.cm.get_cmap("gnuplot2"))  # chg cmap
+                          cmap=plt.get_cmap("gnuplot2"))  # chg cmap
     axs[1][1].set_xscale("log")
     axs[1][0].set_xscale("log")
     plt.show()
diff --git a/rst_docs/user/examples/06-mean_stress.rst b/rst_docs/user/examples/06-mean_stress.rst
index 4f84f8e..13e2991 100644
--- a/rst_docs/user/examples/06-mean_stress.rst
+++ b/rst_docs/user/examples/06-mean_stress.rst
@@ -37,7 +37,7 @@ DNVGL-RP-C203 correction
 ---------------------------
 
 Calculates the mean stress correction according to par.
-2.5 of `DNVGL-RP-C203 <https://bit.ly/3dUZ1OY>`_ which includes
+2.5 of DNVGL-RP-C203, which includes
 an attenuation factor :math:`p` for the stress ranges if the
 following cases:
 
diff --git a/scripts/benchmark_numba_speedups.py b/scripts/benchmark_numba_speedups.py
new file mode 100644
index 0000000..db5976a
--- /dev/null
+++ b/scripts/benchmark_numba_speedups.py
@@ -0,0 +1,677 @@
+"""Benchmark numba call patterns used by py-fatigue.
+
+The script keeps package logic unchanged. It compares current public call
+paths with equivalent call patterns that hoist JIT compilation out of tight
+loops, and it inventories the numba-decorated functions found in the package.
+"""
+
+from __future__ import annotations
+
+import argparse
+import ast
+import gc
+import json
+import statistics
+import sys
+import time
+from contextlib import redirect_stdout
+from dataclasses import asdict, dataclass
+from io import StringIO
+from pathlib import Path
+from typing import Any, Callable, Iterable
+
+import numpy as np
+from rich.console import Console
+from rich.table import Table
+from rich_argparse import RichHelpFormatter
+
+from py_fatigue.cycle_count.rainflow import findcross, findtp
+from py_fatigue.damage.crack_growth import CalcCrackGrowth, get_sif
+from py_fatigue.damage.stress_life import (
+    calc_theil_sn_damage,
+    find_sn_curve_intersection,
+)
+from py_fatigue.geometry.cylinder import f_hol_cyl_01
+from py_fatigue.material.crack_growth_curve import ParisCurve
+from py_fatigue.material.sn_curve import SNCurve
+from py_fatigue.mean_stress.corrections import (
+    goodman_haigh_mean_stress_correction,
+)
+from py_fatigue.utils import (
+    compile_specialized_bisect,
+    compile_specialized_newton,
+    numba_bisect,
+    numba_newton,
+    py_bisect,
+    py_newton,
+    to_numba_dict,
+)
+
+PROJECT_ROOT = Path(__file__).resolve().parents[1]
+PACKAGE_ROOT = PROJECT_ROOT / "py_fatigue"
+
+
+@dataclass(frozen=True)
+class BenchmarkResult:
+    """A single benchmark row."""
+
+    group: str
+    name: str
+    calls: int
+    seconds: float | None
+    per_call_ms: float | None
+    speedup_vs_group_base: float | None = None
+    status: str = "ok"
+    note: str = ""
+
+
+@dataclass(frozen=True)
+class InventoryItem:
+    """A numba-related symbol discovered by AST inspection."""
+
+    file: str
+    line: int
+    kind: str
+    name: str
+    mechanism: str
+
+
+def quadratic_bisect(x_value: float) -> float:
+    """Simple scalar residual for bisection benchmarks."""
+
+    return x_value * x_value - 2.0
+
+
+def quadratic_newton(x_value: float, target: float) -> float:
+    """Simple scalar residual for Newton benchmarks."""
+
+    return x_value * x_value - target
+
+
+def time_repeated(
+    func: Callable[[], Any],
+    calls: int,
+    repeats: int,
+) -> float:
+    """Return the median wall time for repeated calls."""
+
+    timings: list[float] = []
+    was_enabled = gc.isenabled()
+    gc.disable()
+    try:
+        for _ in range(repeats):
+            start = time.perf_counter()
+            for _ in range(calls):
+                func()
+            timings.append(time.perf_counter() - start)
+    finally:
+        if was_enabled:
+            gc.enable()
+    return statistics.median(timings)
+
+
+def add_result(
+    results: list[BenchmarkResult],
+    group: str,
+    name: str,
+    calls: int,
+    seconds: float | None,
+    baseline: float | None = None,
+    status: str = "ok",
+    note: str = "",
+) -> None:
+    """Append a benchmark result with derived per-call and speedup values."""
+
+    per_call_ms = None if seconds is None else seconds / calls * 1_000
+    speedup = None
+    if seconds is not None and baseline is not None and seconds > 0:
+        speedup = baseline / seconds
+    results.append(
+        BenchmarkResult(
+            group=group,
+            name=name,
+            calls=calls,
+            seconds=seconds,
+            per_call_ms=per_call_ms,
+            speedup_vs_group_base=speedup,
+            status=status,
+            note=note,
+        )
+    )
+
+
+def measure_case(
+    results: list[BenchmarkResult],
+    group: str,
+    name: str,
+    func: Callable[[], Any],
+    calls: int,
+    repeats: int,
+    baseline: float | None = None,
+) -> float | None:
+    """Run a timed case and append the result."""
+
+    try:
+        elapsed = time_repeated(func, calls=calls, repeats=repeats)
+    except Exception as exc:  # pylint: disable=broad-exception-caught
+        add_result(
+            results,
+            group,
+            name,
+            calls,
+            None,
+            baseline=baseline,
+            status="error",
+            note=f"{type(exc).__name__}: {exc}",
+        )
+        return None
+
+    add_result(results, group, name, calls, elapsed, baseline=baseline)
+    return elapsed
+
+
+def benchmark_root_finders(
+    calls: int,
+    repeats: int,
+) -> list[BenchmarkResult]:
+    """Benchmark current root-finder wrappers against hoisted compilation."""
+
+    results: list[BenchmarkResult] = []
+
+    bisect_base = measure_case(
+        results,
+        "bisect-wrapper",
+        "current numba_bisect(py_func)",
+        lambda: numba_bisect(quadratic_bisect, 0.0, 2.0, 1e-8, 100),
+        calls,
+        repeats,
+    )
+
+    solver = compile_specialized_bisect(quadratic_bisect)
+    solver(0.0, 2.0, 1e-8, 100)
+    measure_case(
+        results,
+        "bisect-wrapper",
+        "hoisted compile_specialized_bisect",
+        lambda: solver(0.0, 2.0, 1e-8, 100),
+        calls,
+        repeats,
+        baseline=bisect_base,
+    )
+    measure_case(
+        results,
+        "bisect-wrapper",
+        "pure Python py_bisect",
+        lambda: py_bisect(quadratic_bisect, 0.0, 2.0, 1e-8, 100),
+        calls,
+        repeats,
+        baseline=bisect_base,
+    )
+
+    newton_base = measure_case(
+        results,
+        "newton-wrapper",
+        "current numba_newton(py_func)",
+        lambda: numba_newton(quadratic_newton, 1.0, 1e-8, 100, 2.0),
+        calls,
+        repeats,
+    )
+
+    newton_solver = compile_specialized_newton(quadratic_newton)
+    newton_solver(1.0, 1e-8, 100, 2.0)
+    measure_case(
+        results,
+        "newton-wrapper",
+        "hoisted compile_specialized_newton",
+        lambda: newton_solver(1.0, 1e-8, 100, 2.0),
+        calls,
+        repeats,
+        baseline=newton_base,
+    )
+    measure_case(
+        results,
+        "newton-wrapper",
+        "pure Python py_newton",
+        lambda: py_newton(quadratic_newton, 1.0, 1e-8, 100, 2.0),
+        calls,
+        repeats,
+        baseline=newton_base,
+    )
+
+    return results
+
+
+def benchmark_public_paths(
+    calls: int,
+    repeats: int,
+) -> list[BenchmarkResult]:
+    """Benchmark public workflows that exercise each numba area."""
+
+    results: list[BenchmarkResult] = []
+
+    sn_curve = SNCurve([4, 5], [15.117, 17.146], endurance=1e9)
+    stress_range = np.linspace(80.0, 240.0, 10_000)
+    cycles = np.logspace(5, 9, 10_000)
+    sn_curve.get_cycles(stress_range)
+    sn_curve.get_stress(cycles)
+
+    sn_base = measure_case(
+        results,
+        "direct-dispatchers",
+        "SNCurve.get_cycles",
+        lambda: sn_curve.get_cycles(stress_range),
+        calls,
+        repeats,
+    )
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "SNCurve.get_stress",
+        lambda: sn_curve.get_stress(cycles),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+
+    paris_curve = ParisCurve(
+        slope=[2.88, 5.1, 8.16, 5.1, 2.88],
+        intercept=[1e-16, 1e-20, 1e-27, 1e-19, 1e-13],
+        threshold=20.0,
+        critical=2_000.0,
+    )
+    sif_range = np.linspace(20.0, 2_000.0, 10_000)
+    growth_rate = np.logspace(-10, -4, 10_000)
+    paris_curve.get_growth_rate(sif_range)
+    paris_curve.get_sif(growth_rate)
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "ParisCurve.get_growth_rate",
+        lambda: paris_curve.get_growth_rate(sif_range),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "ParisCurve.get_sif",
+        lambda: paris_curve.get_sif(growth_rate),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+
+    signal = np.sin(np.linspace(0.0, 2_000.0, 20_000))
+    signal += 0.15 * np.sin(np.linspace(0.0, 40_000.0, 20_000))
+    findcross(signal)
+    findtp(signal)
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "rainflow.findcross",
+        lambda: findcross(signal),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "rainflow.findtp",
+        lambda: findtp(signal),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+
+    geometry = to_numba_dict(
+        {
+            "initial_depth": 1.0,
+            "outer_diameter": 100.0,
+            "thickness": 10.0,
+            "height": 200.0,
+            "width_to_depth_ratio": 2.0,
+        }
+    )
+    f_hol_cyl_01(1.0, geometry)
+    get_sif(100.0, 1.0, "HOL_CYL_01", geometry)
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "f_hol_cyl_01",
+        lambda: f_hol_cyl_01(1.0, geometry),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "crack_growth.get_sif",
+        lambda: get_sif(100.0, 1.0, "HOL_CYL_01", geometry),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+
+    cg_stress = np.full(500, 100.0, dtype=np.float64)
+    count_cycle = np.ones(500, dtype=np.float64)
+    slope = np.array([3.0], dtype=np.float64)
+    intercept = np.array([1e-12], dtype=np.float64)
+    inf_geometry = to_numba_dict({"initial_depth": 1.0})
+
+    def build_crack_growth() -> CalcCrackGrowth:
+        with redirect_stdout(StringIO()):
+            return CalcCrackGrowth(
+                cg_stress,
+                count_cycle,
+                slope,
+                intercept,
+                0.0,
+                1e9,
+                "INF_SUR_00",
+                inf_geometry,
+            )
+
+    build_crack_growth()
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "CalcCrackGrowth construction",
+        build_crack_growth,
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+
+    small_stress = np.array([120.0, 180.0, 90.0], dtype=np.float64)
+    small_counts = np.array([10.0, 8.0, 12.0], dtype=np.float64)
+    calc_theil_sn_damage(small_stress, small_counts, sn_curve)
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "calc_theil_sn_damage",
+        lambda: calc_theil_sn_damage(small_stress, small_counts, sn_curve),
+        calls,
+        repeats,
+        baseline=sn_base,
+    )
+
+    amp_in = np.linspace(20.0, 300.0, 40)
+    mean_in = np.linspace(-10.0, 60.0, 40)
+    r_out = np.array([-1.0, -0.5, 0.0, 0.5])
+    goodman_haigh_mean_stress_correction(
+        amp_in,
+        mean_in,
+        r_out,
+        1_000.0,
+        3.0,
+    )
+    measure_case(
+        results,
+        "direct-dispatchers",
+        "goodman_haigh_mean_stress_correction",
+        lambda: goodman_haigh_mean_stress_correction(
+            amp_in,
+            mean_in,
+            r_out,
+            1_000.0,
+            3.0,
+        ),
+        max(1, calls // 5),
+        repeats,
+        baseline=sn_base,
+    )
+
+    try:
+        find_sn_curve_intersection(
+            sn_curve.slope,
+            sn_curve.intercept,
+            sn_curve.endurance,
+            0.1,
+            10.0,
+            1.0,
+            1e15,
+        )
+    except Exception as exc:  # pylint: disable=broad-exception-caught
+        add_result(
+            results,
+            "direct-dispatchers",
+            "find_sn_curve_intersection",
+            1,
+            None,
+            baseline=sn_base,
+            status="error",
+            note=f"{type(exc).__name__}: {exc}",
+        )
+    else:
+        measure_case(
+            results,
+            "direct-dispatchers",
+            "find_sn_curve_intersection",
+            lambda: find_sn_curve_intersection(
+                sn_curve.slope,
+                sn_curve.intercept,
+                sn_curve.endurance,
+                0.1,
+                10.0,
+                1.0,
+                1e15,
+            ),
+            max(1, calls // 10),
+            repeats,
+            baseline=sn_base,
+        )
+
+    return results
+
+
+def decorator_text(node: ast.AST) -> str:
+    """Return best-effort source text for a decorator or call."""
+
+    try:
+        return ast.unparse(node)
+    except Exception:  # pragma: no cover
+        return type(node).__name__
+
+
+def iter_python_files(paths: Iterable[Path]) -> Iterable[Path]:
+    """Yield Python files below the provided paths."""
+
+    for path in paths:
+        if path.is_file() and path.suffix == ".py":
+            yield path
+        elif path.is_dir():
+            yield from sorted(path.rglob("*.py"))
+
+
+def inspect_inventory() -> list[InventoryItem]:
+    """Inspect package sources for numba decorators and wrapper call sites."""
+
+    items: list[InventoryItem] = []
+    for py_file in iter_python_files([PACKAGE_ROOT]):
+        tree = ast.parse(py_file.read_text(encoding="utf-8"))
+        rel_path = py_file.relative_to(PROJECT_ROOT).as_posix()
+        for node in ast.walk(tree):
+            if isinstance(
+                node, (ast.FunctionDef, ast.AsyncFunctionDef, ast.ClassDef)
+            ):
+                for decorator in node.decorator_list:
+                    text = decorator_text(decorator)
+                    if any(
+                        token in text
+                        for token in ("njit", "jitclass", "vectorize")
+                    ):
+                        items.append(
+                            InventoryItem(
+                                file=rel_path,
+                                line=node.lineno,
+                                kind=type(node).__name__,
+                                name=node.name,
+                                mechanism=f"@{text}",
+                            )
+                        )
+            if isinstance(node, ast.Assign) and isinstance(
+                node.value, ast.Call
+            ):
+                text = decorator_text(node.value.func)
+                if text.endswith(
+                    (
+                        "compile_specialized_bisect",
+                        "compile_specialized_newton",
+                    )
+                ):
+                    for target in node.targets:
+                        if isinstance(target, ast.Name):
+                            items.append(
+                                InventoryItem(
+                                    file=rel_path,
+                                    line=node.lineno,
+                                    kind="Assign",
+                                    name=target.id,
+                                    mechanism=text,
+                                )
+                            )
+            if isinstance(node, ast.Call):
+                text = decorator_text(node.func)
+                if text in {"numba_bisect", "numba_newton"}:
+                    items.append(
+                        InventoryItem(
+                            file=rel_path,
+                            line=node.lineno,
+                            kind="Call",
+                            name=text,
+                            mechanism="generic wrapper call",
+                        )
+                    )
+    return sorted(items, key=lambda item: (item.file, item.line, item.name))
+
+
+def render_inventory(console: Console, items: list[InventoryItem]) -> None:
+    """Print the numba inventory."""
+
+    table = Table(title="Numba inventory")
+    table.add_column("File")
+    table.add_column("Line", justify="right")
+    table.add_column("Kind")
+    table.add_column("Name")
+    table.add_column("Mechanism")
+    for item in items:
+        table.add_row(
+            item.file,
+            str(item.line),
+            item.kind,
+            item.name,
+            item.mechanism,
+        )
+    console.print(table)
+
+
+def render_results(console: Console, results: list[BenchmarkResult]) -> None:
+    """Print benchmark results."""
+
+    table = Table(title="Numba benchmark results")
+    table.add_column("Group")
+    table.add_column("Case")
+    table.add_column("Calls", justify="right")
+    table.add_column("Total (s)", justify="right")
+    table.add_column("ms/call", justify="right")
+    table.add_column("Speedup", justify="right")
+    table.add_column("Status")
+    table.add_column("Note")
+    for result in results:
+        total = "-" if result.seconds is None else f"{result.seconds:.6f}"
+        per_call = (
+            "-" if result.per_call_ms is None else f"{result.per_call_ms:.3f}"
+        )
+        speedup = (
+            "-"
+            if result.speedup_vs_group_base is None
+            else f"{result.speedup_vs_group_base:.1f}x"
+        )
+        table.add_row(
+            result.group,
+            result.name,
+            str(result.calls),
+            total,
+            per_call,
+            speedup,
+            result.status,
+            result.note,
+        )
+    console.print(table)
+
+
+def parse_args(argv: list[str]) -> argparse.Namespace:
+    """Parse command-line arguments."""
+
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=RichHelpFormatter,
+    )
+    parser.add_argument(
+        "--calls",
+        type=int,
+        default=10,
+        help="Calls per benchmark repeat. Higher values amplify differences.",
+    )
+    parser.add_argument(
+        "--repeats",
+        type=int,
+        default=3,
+        help="Number of repeats used to compute the median wall time.",
+    )
+    parser.add_argument(
+        "--skip-inventory",
+        action="store_true",
+        help="Do not print the static numba inventory.",
+    )
+    parser.add_argument(
+        "--skip-benchmarks",
+        action="store_true",
+        help="Do not run timing benchmarks.",
+    )
+    parser.add_argument(
+        "--json-output",
+        type=Path,
+        default=None,
+        help="Optional JSON file for machine-readable results.",
+    )
+    return parser.parse_args(argv)
+
+
+def main(argv: list[str] | None = None) -> int:
+    """Run the inventory and benchmarks."""
+
+    args = parse_args(sys.argv[1:] if argv is None else argv)
+    console = Console()
+
+    inventory: list[InventoryItem] = []
+    results: list[BenchmarkResult] = []
+
+    if not args.skip_inventory:
+        inventory = inspect_inventory()
+        render_inventory(console, inventory)
+
+    if not args.skip_benchmarks:
+        results.extend(benchmark_root_finders(args.calls, args.repeats))
+        results.extend(benchmark_public_paths(args.calls, args.repeats))
+        render_results(console, results)
+
+    if args.json_output is not None:
+        payload = {
+            "inventory": [asdict(item) for item in inventory],
+            "results": [asdict(result) for result in results],
+        }
+        args.json_output.write_text(
+            json.dumps(payload, indent=2),
+            encoding="utf-8",
+        )
+        console.print(f"Wrote JSON results to {args.json_output}")
+
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/scripts/numba_benchmark_results.json b/scripts/numba_benchmark_results.json
new file mode 100644
index 0000000..f85d780
--- /dev/null
+++ b/scripts/numba_benchmark_results.json
@@ -0,0 +1,326 @@
+{
+  "inventory": [
+    {
+      "file": "py_fatigue/cycle_count/rainflow.py",
+      "line": 209,
+      "kind": "FunctionDef",
+      "name": "_findcross",
+      "mechanism": "@njit(int64(int64[:], int8[:]), cache=True)"
+    },
+    {
+      "file": "py_fatigue/damage/crack_growth.py",
+      "line": 47,
+      "kind": "ClassDef",
+      "name": "CalcCrackGrowth",
+      "mechanism": "@jitclass(spec)"
+    },
+    {
+      "file": "py_fatigue/damage/crack_growth.py",
+      "line": 333,
+      "kind": "FunctionDef",
+      "name": "get_geometry_factor",
+      "mechanism": "@nb.njit(nb.float64(nb.float64, nb.types.unicode_type, nb.types.DictType(nb.types.unicode_type, nb.float64)), fastmath=True, cache=True)"
+    },
+    {
+      "file": "py_fatigue/damage/crack_growth.py",
+      "line": 360,
+      "kind": "FunctionDef",
+      "name": "get_sif",
+      "mechanism": "@nb.njit(nb.types.UniTuple(nb.float64, 2)(nb.float64, nb.float64, nb.types.unicode_type, nb.types.DictType(nb.types.unicode_type, nb.float64)), fastmath=True, cache=True)"
+    },
+    {
+      "file": "py_fatigue/damage/stress_life.py",
+      "line": 1269,
+      "kind": "FunctionDef",
+      "name": "_calc_theil_sn_damage",
+      "mechanism": "@nb.njit(fastmath=True, cache=True)"
+    },
+    {
+      "file": "py_fatigue/geometry/cylinder.py",
+      "line": 246,
+      "kind": "FunctionDef",
+      "name": "f_hol_cyl_01",
+      "mechanism": "@nb.njit(nb.float64(nb.float64, nb.types.DictType(nb.types.unicode_type, nb.float64)), cache=True, fastmath=True)"
+    },
+    {
+      "file": "py_fatigue/material/crack_growth_curve.py",
+      "line": 1077,
+      "kind": "FunctionDef",
+      "name": "get_array_min",
+      "mechanism": "@nb.njit(cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/crack_growth_curve.py",
+      "line": 1093,
+      "kind": "FunctionDef",
+      "name": "_calc_growth_rate",
+      "mechanism": "@nb.njit(fastmath=False, cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/crack_growth_curve.py",
+      "line": 1139,
+      "kind": "FunctionDef",
+      "name": "_calc_sif",
+      "mechanism": "@nb.njit(fastmath=False, cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/sn_curve.py",
+      "line": 956,
+      "kind": "FunctionDef",
+      "name": "get_sn_array_min",
+      "mechanism": "@nb.njit(cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/sn_curve.py",
+      "line": 972,
+      "kind": "FunctionDef",
+      "name": "_calc_cycles",
+      "mechanism": "@nb.njit(fastmath=False, cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/sn_curve.py",
+      "line": 996,
+      "kind": "FunctionDef",
+      "name": "_calc_cycles_2",
+      "mechanism": "@nb.njit(fastmath=False, cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/sn_curve.py",
+      "line": 1051,
+      "kind": "FunctionDef",
+      "name": "_calc_stress",
+      "mechanism": "@nb.njit(fastmath=False, cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/sn_curve.py",
+      "line": 1078,
+      "kind": "FunctionDef",
+      "name": "_calc_stress_2",
+      "mechanism": "@nb.njit(fastmath=False, cache=True)"
+    },
+    {
+      "file": "py_fatigue/material/sn_curve.py",
+      "line": 1183,
+      "kind": "Assign",
+      "name": "__jit_sn_curve_residuals",
+      "mechanism": "compile_specialized_bisect"
+    },
+    {
+      "file": "py_fatigue/mean_stress/corrections.py",
+      "line": 423,
+      "kind": "Call",
+      "name": "numba_newton",
+      "mechanism": "generic wrapper call"
+    },
+    {
+      "file": "py_fatigue/mean_stress/corrections.py",
+      "line": 582,
+      "kind": "Assign",
+      "name": "__jit_goodman_equation",
+      "mechanism": "compile_specialized_newton"
+    },
+    {
+      "file": "py_fatigue/utils.py",
+      "line": 762,
+      "kind": "FunctionDef",
+      "name": "raise_non_scalar_numba_result",
+      "mechanism": "@nb.njit(nb.float64())"
+    },
+    {
+      "file": "py_fatigue/utils.py",
+      "line": 1001,
+      "kind": "Assign",
+      "name": "jit_bisect_func",
+      "mechanism": "compile_specialized_bisect"
+    },
+    {
+      "file": "py_fatigue/utils.py",
+      "line": 1014,
+      "kind": "Assign",
+      "name": "jit_newton_func",
+      "mechanism": "compile_specialized_newton"
+    }
+  ],
+  "results": [
+    {
+      "group": "bisect-wrapper",
+      "name": "current numba_bisect(py_func)",
+      "calls": 8,
+      "seconds": 1.2343996786512434e-05,
+      "per_call_ms": 0.0015429995983140543,
+      "speedup_vs_group_base": null,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "bisect-wrapper",
+      "name": "hoisted compile_specialized_bisect",
+      "calls": 8,
+      "seconds": 6.1780010582879186e-06,
+      "per_call_ms": 0.0007722501322859898,
+      "speedup_vs_group_base": 1.9980567614103437,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "bisect-wrapper",
+      "name": "pure Python py_bisect",
+      "calls": 8,
+      "seconds": 9.328400483354926e-05,
+      "per_call_ms": 0.011660500604193658,
+      "speedup_vs_group_base": 0.13232704586963617,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "newton-wrapper",
+      "name": "current numba_newton(py_func)",
+      "calls": 8,
+      "seconds": 1.6310994396917522e-05,
+      "per_call_ms": 0.0020388742996146902,
+      "speedup_vs_group_base": null,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "newton-wrapper",
+      "name": "hoisted compile_specialized_newton",
+      "calls": 8,
+      "seconds": 8.942995918914676e-06,
+      "per_call_ms": 0.0011178744898643345,
+      "speedup_vs_group_base": 1.823884808268707,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "newton-wrapper",
+      "name": "pure Python py_newton",
+      "calls": 8,
+      "seconds": 4.595100472215563e-05,
+      "per_call_ms": 0.005743875590269454,
+      "speedup_vs_group_base": 0.3549649130751879,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "SNCurve.get_cycles",
+      "calls": 8,
+      "seconds": 0.0033126759954029694,
+      "per_call_ms": 0.4140844994253712,
+      "speedup_vs_group_base": null,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "SNCurve.get_stress",
+      "calls": 8,
+      "seconds": 0.0035862310032825917,
+      "per_call_ms": 0.44827887541032396,
+      "speedup_vs_group_base": 0.9237207509418025,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "ParisCurve.get_growth_rate",
+      "calls": 8,
+      "seconds": 0.010180105993640609,
+      "per_call_ms": 1.2725132492050761,
+      "speedup_vs_group_base": 0.32540682753915906,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "ParisCurve.get_sif",
+      "calls": 8,
+      "seconds": 0.017416592003428377,
+      "per_call_ms": 2.177074000428547,
+      "speedup_vs_group_base": 0.19020230793434698,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "rainflow.findcross",
+      "calls": 8,
+      "seconds": 0.0010041150089818984,
+      "per_call_ms": 0.1255143761227373,
+      "speedup_vs_group_base": 3.299100168577092,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "rainflow.findtp",
+      "calls": 8,
+      "seconds": 0.0013474489969667047,
+      "per_call_ms": 0.1684311246208381,
+      "speedup_vs_group_base": 2.4584796922631313,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "f_hol_cyl_01",
+      "calls": 8,
+      "seconds": 2.6398993213661015e-05,
+      "per_call_ms": 0.003299874151707627,
+      "speedup_vs_group_base": 125.48493681526907,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "crack_growth.get_sif",
+      "calls": 8,
+      "seconds": 4.6574990847148e-05,
+      "per_call_ms": 0.0058218738558935,
+      "speedup_vs_group_base": 71.12563921428703,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "CalcCrackGrowth construction",
+      "calls": 8,
+      "seconds": 0.005733437996241264,
+      "per_call_ms": 0.716679749530158,
+      "speedup_vs_group_base": 0.5777817772817458,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "calc_theil_sn_damage",
+      "calls": 8,
+      "seconds": 0.00019102801161352545,
+      "per_call_ms": 0.02387850145169068,
+      "speedup_vs_group_base": 17.341310143063964,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "goodman_haigh_mean_stress_correction",
+      "calls": 1,
+      "seconds": 0.0008232050022343174,
+      "per_call_ms": 0.8232050022343174,
+      "speedup_vs_group_base": 4.024120342334907,
+      "status": "ok",
+      "note": ""
+    },
+    {
+      "group": "direct-dispatchers",
+      "name": "find_sn_curve_intersection",
+      "calls": 1,
+      "seconds": 6.591700366698205e-05,
+      "per_call_ms": 0.06591700366698205,
+      "speedup_vs_group_base": 50.25525753777996,
+      "status": "ok",
+      "note": ""
+    }
+  ]
+}
\ No newline at end of file
diff --git a/tests/damage/test_stress_life.py b/tests/damage/test_stress_life.py
index 2332ea1..9d54f5a 100644
--- a/tests/damage/test_stress_life.py
+++ b/tests/damage/test_stress_life.py
@@ -4,13 +4,13 @@
 damage calculation methods in the stress-life approach.
 """
 
-
 # Standard imports
 import datetime as dt
 import os
 import sys
 
 # Non-standard imports
+import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 import pytest
@@ -18,6 +18,7 @@
 from hypothesis import strategies as hy
 
 import py_fatigue.damage as damage
+
 # Local imports
 from py_fatigue import CycleCount, SNCurve
 from tests.cycle_count.test_cycle_count import CC_RF_1
@@ -173,7 +174,6 @@ def test_palmgren_miner_constant_load(
             cc_obj.unit = "m"
             damage.get_pm(cc_obj, sn_curve)
 
-
     @pytest.mark.parametrize(
         "sn_curve", [(DNV_B1A), (DNV_B1A_END), (DNV_B1W), (DNV_B1C)]
     )
@@ -198,9 +198,7 @@ def test_palmgren_miner_variable_load(self, sn_curve: SNCurve):
         )
         assert damage_pm == pytest.approx(damage_ref, 1e-12)
 
-    @pytest.mark.parametrize(
-        "cc", [(CC_TS_1), (CC_TS_2), (CC_TS_3), (CC_RF_1)]
-    )
+    @pytest.mark.parametrize("cc", [(CC_TS_1), (CC_TS_2), (CC_TS_3), (CC_RF_1)])
     @pytest.mark.parametrize("exponent", [(3), (4), (5)])
     @pytest.mark.parametrize("eq_cycles", [(1e6), (2e6), (1e7)])
     # fmt: on
@@ -266,7 +264,7 @@ def test_damage_equivalent_moment(
             assert "outer_radius must be greater" in ve.value.args[0]
         a_i = np.pi / 4 * (r_o**4 - r_i**4)
         dem = damage.get_dem(r_o, r_i, cc, exponent, eq_cycles)
-        dem_ref = a_i * 1E6 * damage.get_des(cc, exponent, eq_cycles) / r_o
+        dem_ref = a_i * 1e6 * damage.get_des(cc, exponent, eq_cycles) / r_o
         assert dem == pytest.approx(dem_ref, 1e-12)
 
     # fmt: off
@@ -309,10 +307,14 @@ def test_miner_pandas_accessor_constant_load(
         assert df_d.sn_curve == sn_curve
         assert df_d._metadata["name"] == cc_obj.name
         assert df_d._metadata["timestamp"] == cc_obj.timestamp
-        assert df_d._metadata["mean_stress_corrected"] == \
-            cc_obj.mean_stress_corrected
-        assert df_d._metadata["stress_concentration_factor"] == \
-            cc_obj.stress_concentration_factor
+        assert (
+            df_d._metadata["mean_stress_corrected"]
+            == cc_obj.mean_stress_corrected
+        )
+        assert (
+            df_d._metadata["stress_concentration_factor"]
+            == cc_obj.stress_concentration_factor
+        )
         assert df_d._metadata["nr_small_cycles"] == cc_obj.nr_small_cycles
         assert df_d._metadata["lffd_solved"] == cc_obj.lffd_solved
         assert isinstance(df_d, pd.DataFrame)
@@ -327,9 +329,7 @@ def test_miner_pandas_accessor_constant_load(
     @pytest.mark.parametrize(
         "sn_curve", [DNV_B1A, DNV_B1A_END, DNV_B1W, DNV_B1C]
     )
-    @pytest.mark.parametrize(
-        "cc", [CC_TS_1, CC_TS_3]
-    )
+    @pytest.mark.parametrize("cc", [CC_TS_1, CC_TS_3])
     # fmt: on
     def test_miner_pandas_accessor_variable_load(
         self, cc: CycleCount, sn_curve: SNCurve
@@ -353,23 +353,29 @@ def test_miner_pandas_accessor_variable_load(
         assert df_d.sn_curve == sn_curve
         assert df_d._metadata["name"] == cc.name
         assert df_d._metadata["timestamp"] == cc.timestamp
-        assert df_d._metadata["mean_stress_corrected"] == \
-            cc.mean_stress_corrected
-        assert df_d._metadata["stress_concentration_factor"] == \
-            cc.stress_concentration_factor
+        assert (
+            df_d._metadata["mean_stress_corrected"] == cc.mean_stress_corrected
+        )
+        assert (
+            df_d._metadata["stress_concentration_factor"]
+            == cc.stress_concentration_factor
+        )
         assert df_d._metadata["nr_small_cycles"] == cc.nr_small_cycles
         assert df_d._metadata["lffd_solved"] == cc.lffd_solved
         assert isinstance(df_d, pd.DataFrame)
         assert df_d["pm_damage"].sum() == pytest.approx(
             np.sum(damage.get_pm(cc, sn_curve)), 1e-12
         )
+        fig, ax = plt.subplots()
+        df_d.miner.plot_histogram(fig=fig, ax=ax)
+        plt.close(fig)
+        with pytest.raises(AttributeError):
+            df_d.miner.damage(sn_curve)
 
-    @pytest.mark.parametrize(
-        "cc,", [CC_TS_1, CC_TS_3]
-    )
+    @pytest.mark.parametrize("cc,", [CC_TS_1, CC_TS_3])
     @given(
         slope=hy.floats(min_value=3, max_value=20),
-        n_eq=hy.floats(min_value=1E5, max_value=1e10),
+        n_eq=hy.floats(min_value=1e5, max_value=1e10),
     )
     # fmt: on
     def test_des_pandas_accessor_variable_load(
@@ -394,28 +400,28 @@ def test_des_pandas_accessor_variable_load(
         assert isinstance(df, pd.DataFrame)
         assert df._metadata["name"] == cc.name
         assert df._metadata["timestamp"] == cc.timestamp
-        assert df._metadata["mean_stress_corrected"] == \
-            cc.mean_stress_corrected
-        assert df._metadata["stress_concentration_factor"] == \
-            cc.stress_concentration_factor
+        assert df._metadata["mean_stress_corrected"] == cc.mean_stress_corrected
+        assert (
+            df._metadata["stress_concentration_factor"]
+            == cc.stress_concentration_factor
+        )
         assert df._metadata["nr_small_cycles"] == cc.nr_small_cycles
         assert df._metadata["lffd_solved"] == cc.lffd_solved
         assert isinstance(df, pd.DataFrame)
-        assert df.miner.des(slope=slope, equivalent_cycles=n_eq) == \
-            pytest.approx(
-                damage.get_des(cc, slope, equivalent_cycles=n_eq), 1e-12
-            )
+        assert df.miner.des(
+            slope=slope, equivalent_cycles=n_eq
+        ) == pytest.approx(
+            damage.get_des(cc, slope, equivalent_cycles=n_eq), 1e-12
+        )
 
-    @pytest.mark.parametrize(
-        "cc,", [CC_TS_1, CC_TS_3]
-    )
+    @pytest.mark.parametrize("cc,", [CC_TS_1, CC_TS_3])
     @given(
         slope=hy.floats(min_value=3, max_value=20),
-        n_eq=hy.floats(min_value=1E5, max_value=1e10),
+        n_eq=hy.floats(min_value=1e5, max_value=1e10),
     )
     @pytest.mark.parametrize(
         "r_i, r_o",
-        [(3, 3.5), (4000, 4100), pytest.param(2, 1, marks=pytest.mark.xfail)]
+        [(3, 3.5), (4000, 4100), pytest.param(2, 1, marks=pytest.mark.xfail)],
     )
     # fmt: on
     def test_dem_pandas_accessor_variable_load(
@@ -444,10 +450,11 @@ def test_dem_pandas_accessor_variable_load(
         assert isinstance(df, pd.DataFrame)
         assert df._metadata["name"] == cc.name
         assert df._metadata["timestamp"] == cc.timestamp
-        assert df._metadata["mean_stress_corrected"] == \
-            cc.mean_stress_corrected
-        assert df._metadata["stress_concentration_factor"] == \
-            cc.stress_concentration_factor
+        assert df._metadata["mean_stress_corrected"] == cc.mean_stress_corrected
+        assert (
+            df._metadata["stress_concentration_factor"]
+            == cc.stress_concentration_factor
+        )
         assert df._metadata["nr_small_cycles"] == cc.nr_small_cycles
         assert df._metadata["lffd_solved"] == cc.lffd_solved
         assert isinstance(df, pd.DataFrame)
@@ -455,23 +462,23 @@ def test_dem_pandas_accessor_variable_load(
             outer_radius=r_o,
             inner_radius=r_i,
             slope=slope,
-            equivalent_cycles=n_eq
+            equivalent_cycles=n_eq,
         ) == pytest.approx(
             damage.get_dem(
                 outer_radius=r_o,
                 inner_radius=r_i,
                 cycle_count=cc,
                 slope=slope,
-                equivalent_cycles=n_eq
+                equivalent_cycles=n_eq,
             ),
-            1e-12
+            1e-12,
         )
 
+
 class TestGassner:
-    """Test the shift factor calculation related with the Gassner curve
-    """
+    """Test the shift factor calculation related with the Gassner curve"""
 
-        # fmt: off
+    # fmt: off
     @settings(deadline=None)
     @pytest.mark.parametrize(
         "sn_curve", [DNV_B1C, DNV_C_C, DNV_E_C,
@@ -499,6 +506,7 @@ def test_g_pandas_accessor_constant_load(
             The length of the history to use.
         """
         import py_fatigue as pf
+
         time_reversals = [peak * (-1) ** _ for _ in range(len_hist)]
         cc_obj = pf.CycleCount.from_timeseries(
             time_reversals,
@@ -515,26 +523,34 @@ def test_g_pandas_accessor_constant_load(
         assert df_g.sn_curve == sn_curve
         assert df_g._metadata["name"] == cc_obj.name
         assert df_g._metadata["timestamp"] == cc_obj.timestamp
-        assert df_g._metadata["mean_stress_corrected"] == \
-            cc_obj.mean_stress_corrected
-        assert df_g._metadata["stress_concentration_factor"] == \
-            cc_obj.stress_concentration_factor
+        assert (
+            df_g._metadata["mean_stress_corrected"]
+            == cc_obj.mean_stress_corrected
+        )
+        assert (
+            df_g._metadata["stress_concentration_factor"]
+            == cc_obj.stress_concentration_factor
+        )
         assert df_g._metadata["nr_small_cycles"] == cc_obj.nr_small_cycles
         assert df_g._metadata["lffd_solved"] == cc_obj.lffd_solved
         assert isinstance(df_g, pd.DataFrame)
         assert df_g["shift_factor"].sum() == pytest.approx(1, 1e-12)
 
+    @pytest.mark.parametrize("cc", [CC_TS_1, CC_TS_3])
     @pytest.mark.parametrize(
-        "cc", [CC_TS_1, CC_TS_3]
-    )
-    @pytest.mark.parametrize(
-        "sn_curve", [DNV_B1C, DNV_C_C, DNV_E_C,
-                     pytest.param(DNV_B1A, marks=pytest.mark.xfail)
-        ]
+        "sn_curve",
+        [
+            DNV_B1C,
+            DNV_C_C,
+            DNV_E_C,
+            pytest.param(DNV_B1A, marks=pytest.mark.xfail),
+        ],
     )
     # fmt: on
     def test_g_pandas_accessor_variable_load(
-        self, cc: CycleCount, sn_curve: SNCurve,
+        self,
+        cc: CycleCount,
+        sn_curve: SNCurve,
     ):
         """Test the Gassner shift factor that has to be less than one
         for variable amplitude stress histories.
@@ -551,19 +567,24 @@ def test_g_pandas_accessor_variable_load(
         assert isinstance(df, pd.DataFrame)
         assert df._metadata["name"] == cc.name
         assert df._metadata["timestamp"] == cc.timestamp
-        assert df._metadata["mean_stress_corrected"] == \
-            cc.mean_stress_corrected
-        assert df._metadata["stress_concentration_factor"] == \
-            cc.stress_concentration_factor
+        assert df._metadata["mean_stress_corrected"] == cc.mean_stress_corrected
+        assert (
+            df._metadata["stress_concentration_factor"]
+            == cc.stress_concentration_factor
+        )
         assert df._metadata["nr_small_cycles"] == cc.nr_small_cycles
         assert df._metadata["lffd_solved"] == cc.lffd_solved
         assert df["shift_factor"].sum() < 1
 
+
 @pytest.mark.parametrize("sn_curve", [DNV_B1C, DNV_C_C, DNV_E_C, DNV_B1A])
-@pytest.mark.parametrize("load", [
-    [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 3, -3, 3, -3, 3, -3, 3, -3],
-    [3, -3, 3, -3, 3, -3, 3, -3, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1]
-])
+@pytest.mark.parametrize(
+    "load",
+    [
+        [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 3, -3, 3, -3, 3, -3, 3, -3],
+        [3, -3, 3, -3, 3, -3, 3, -3, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1],
+    ],
+)
 @pytest.mark.parametrize("base_exponent", [0.99, 1, 1.01])
 def test_leve_damage_rule(load: list, sn_curve: SNCurve, base_exponent: float):
     """Test the nonlinear damage calculation.
@@ -585,7 +606,7 @@ def test_leve_damage_rule(load: list, sn_curve: SNCurve, base_exponent: float):
         name="Test_CC",
     )
     d_nl = damage.get_nonlinear_damage(
-        'leve', cc, sn_curve, base_exponent=base_exponent
+        "leve", cc, sn_curve, base_exponent=base_exponent
     )
     d_l = np.sum(damage.get_pm(cc, sn_curve))
     if base_exponent == 1:
@@ -595,6 +616,7 @@ def test_leve_damage_rule(load: list, sn_curve: SNCurve, base_exponent: float):
     if base_exponent > 1:
         assert d_nl[-1] < d_l
 
+
 @pytest.mark.parametrize("sn_curve", [DNV_B1C, DNV_C_C, DNV_E_C, DNV_B1A])
 @pytest.mark.parametrize("rule", ["pavlou", "manson", "si jian"])
 @pytest.mark.filterwarnings("ignore::UserWarning")
@@ -650,10 +672,12 @@ class TestDamageExponents:
     )
     def test_calc_damage_exponents_no_kwargs(self, damage_rule):
         """Test the _calc_damage_exponents function with no kwargs."""
-        stress_range = np.array([100., 200., 300.])
+        stress_range = np.array([100.0, 200.0, 300.0])
         if "manson" in damage_rule:
             with pytest.raises(ValueError, match="sn_curve must be provided"):
-                damage.stress_life._calc_damage_exponents(damage_rule, stress_range)
+                damage.stress_life._calc_damage_exponents(
+                    damage_rule, stress_range
+                )
         else:
             if damage_rule == "pavlou":
                 with pytest.warns(UserWarning, match="base_exponent"):
@@ -699,7 +723,10 @@ def test_calc_damage_exponents_manson(self):
         base_exponent = 0.5
         sn_curve = SNCurve([3, 5, 7], [10.970, 13.617, 16])
         exponents = damage.stress_life._calc_damage_exponents(
-            "manson", stress_range, sn_curve=sn_curve, base_exponent=base_exponent
+            "manson",
+            stress_range,
+            sn_curve=sn_curve,
+            base_exponent=base_exponent,
         )
         assert isinstance(exponents, np.ndarray)
         assert np.allclose(
@@ -717,14 +744,20 @@ def test_calc_damage_exponents_leve(self):
             "leve", stress_range, base_exponent=base_exponent
         )
         assert isinstance(exponents, np.ndarray)
-        assert np.allclose(exponents, base_exponent * np.ones(len(stress_range)))
+        assert np.allclose(
+            exponents, base_exponent * np.ones(len(stress_range))
+        )
 
     def test_calc_damage_exponents_si_jian(self):
         """Test the _calc_damage_exponents function with si jian rule."""
         stress_range = np.array([100, 200, 300])
-        exponents = damage.stress_life._calc_damage_exponents("si jian", stress_range)
+        exponents = damage.stress_life._calc_damage_exponents(
+            "si jian", stress_range
+        )
         assert isinstance(exponents, np.ndarray)
-        assert np.allclose(exponents, damage.stress_life.calc_si_jian_exponents(stress_range))
+        assert np.allclose(
+            exponents, damage.stress_life.calc_si_jian_exponents(stress_range)
+        )
 
     def test_calc_damage_exponents_unknown_rule(self):
         """Test the _calc_damage_exponents function with unknown rule."""
@@ -732,15 +765,19 @@ def test_calc_damage_exponents_unknown_rule(self):
         with pytest.raises(ValueError, match="Unknown damage rule: unknown"):
             damage.stress_life._calc_damage_exponents("unknown", stress_range)
 
+
 @pytest.mark.parametrize("sn_curve", [DNV_B1C, DNV_C_C, DNV_E_C, DNV_B1A])
-@pytest.mark.parametrize("load", [
-    [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 3, -3, 3, -3, 3, -3, 3, -3],
-    [3, -3, 3, -3, 3, -3, 3, -3, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1]
-])
-@pytest.mark.parametrize("damage_bands", [
-    [0, 0.2, 0.4, 0.6, 0.8, 1],
-    [0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]
-])
+@pytest.mark.parametrize(
+    "load",
+    [
+        [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 3, -3, 3, -3, 3, -3, 3, -3],
+        [3, -3, 3, -3, 3, -3, 3, -3, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1],
+    ],
+)
+@pytest.mark.parametrize(
+    "damage_bands",
+    [[0, 0.2, 0.4, 0.6, 0.8, 1], [0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]],
+)
 def test_nonlinear_damage_dca(
     load: list, sn_curve: SNCurve, damage_bands: list
 ):
@@ -765,15 +802,19 @@ def test_nonlinear_damage_dca(
         name="Test_CC",
     )
     d_nl, _, _, _ = damage.stress_life.get_nonlinear_damage_with_dca(
-        'pavlou', cc, sn_curve, np.asarray(damage_bands)
+        "pavlou", cc, sn_curve, np.asarray(damage_bands)
     )
     assert isinstance(d_nl, np.ndarray)
 
+
 @pytest.mark.parametrize("sn_curve", [DNV_B1C, DNV_C_C, DNV_E_C, DNV_B1A])
-@pytest.mark.parametrize("load", [
-    [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 3, -3, 3, -3, 3, -3, 3, -3],
-    [3, -3, 3, -3, 3, -3, 3, -3, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1]
-])
+@pytest.mark.parametrize(
+    "load",
+    [
+        [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 3, -3, 3, -3, 3, -3, 3, -3],
+        [3, -3, 3, -3, 3, -3, 3, -3, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1],
+    ],
+)
 def test_theil_damage_rule(load: list, sn_curve: SNCurve):
     """Test the nonlinear damage calculation with Theil's method.
 
@@ -794,7 +835,6 @@ def test_theil_damage_rule(load: list, sn_curve: SNCurve):
         name="Test_CC",
     )
     d_nl, _, _, _ = damage.stress_life.get_nonlinear_damage_with_dca(
-        'theil', cc, sn_curve, damage_bands=np.array([0, 0.2, 0.4, 0.6, 0.8, 1])
+        "theil", cc, sn_curve, damage_bands=np.array([0, 0.2, 0.4, 0.6, 0.8, 1])
     )
     assert isinstance(d_nl, np.ndarray)
-
diff --git a/tests/material/test_crack_growth_curve.py b/tests/material/test_crack_growth_curve.py
index 20c14c0..9775084 100644
--- a/tests/material/test_crack_growth_curve.py
+++ b/tests/material/test_crack_growth_curve.py
@@ -10,6 +10,7 @@
 
 nb.config.DISABLE_JIT = True
 
+
 def test_pure_paris_curve():
     paris_pure = ParisCurve(slope=3.1, intercept=1.2e-14)
     assert paris_pure.linear
@@ -17,7 +18,9 @@ def test_pure_paris_curve():
     assert paris_pure.intercept == 1.2e-14
     assert paris_pure.threshold == 0
     assert paris_pure.critical == np.inf
-    assert paris_pure.threshold_growth_rate == paris_pure.get_growth_rate(paris_pure.threshold)
+    assert paris_pure.threshold_growth_rate == paris_pure.get_growth_rate(
+        paris_pure.threshold
+    )
     assert np.isinf(paris_pure.critical_growth_rate)
     assert np.isinf(paris_pure.get_growth_rate(paris_pure.critical))
     assert len(paris_pure.get_knee_growth_rate()) == 0
@@ -38,10 +41,12 @@ def test_pure_paris_curve():
 
 def test_paris_curve_initialization():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    pc = ParisCurve(slope=slope, intercept=intercept, threshold=threshold, critical=critical)
+    pc = ParisCurve(
+        slope=slope, intercept=intercept, threshold=threshold, critical=critical
+    )
 
     assert not pc.linear
     assert pc.__str__() == pc.name
@@ -80,54 +85,80 @@ def test_paris_curve_initialization():
         _calc_growth_rate(
             knee_sif, pc.slope, pc.intercept, pc.threshold, pc.critical
         ),
-        rtol=1e-2
+        rtol=1e-2,
     )
     assert np.allclose(
         pc.get_sif(knee_gr),
-        _calc_sif(
-            knee_gr, pc.slope, pc.intercept, pc.threshold, pc.critical
-        ),
-        rtol=1e-2
+        _calc_sif(knee_gr, pc.slope, pc.intercept, pc.threshold, pc.critical),
+        rtol=1e-2,
     )
 
     with pytest.raises(TypingError):
-        assert _calc_growth_rate(
-            threshold - 1e-2 , pc.slope, pc.intercept, pc.threshold, pc.critical
-        ) == 0
-    assert _calc_growth_rate(
-        np.array([threshold - 1e-2]),
-        pc.slope,
-        pc.intercept,
-        pc.threshold,
-        pc.critical
-    ) == 0
+        assert (
+            _calc_growth_rate(
+                threshold - 1e-2,
+                pc.slope,
+                pc.intercept,
+                pc.threshold,
+                pc.critical,
+            )
+            == 0
+        )
+    assert (
+        _calc_growth_rate(
+            np.array([threshold - 1e-2]),
+            pc.slope,
+            pc.intercept,
+            pc.threshold,
+            pc.critical,
+        )
+        == 0
+    )
 
     with pytest.raises(TypingError):
-        assert _calc_growth_rate(
-            critical + 1e-2, pc.slope, pc.intercept, pc.threshold, pc.critical
-        ) == np.inf
-    assert _calc_growth_rate(
-        np.array([critical + 1e-2]),
-        pc.slope,
-        pc.intercept,
-        pc.threshold,
-        pc.critical
-    ) == np.inf
+        assert (
+            _calc_growth_rate(
+                critical + 1e-2,
+                pc.slope,
+                pc.intercept,
+                pc.threshold,
+                pc.critical,
+            )
+            == np.inf
+        )
+    assert (
+        _calc_growth_rate(
+            np.array([critical + 1e-2]),
+            pc.slope,
+            pc.intercept,
+            pc.threshold,
+            pc.critical,
+        )
+        == np.inf
+    )
+
 
 def test_paris_curve_from_knee_points():
     knee_sif = [20, 100, 500, 1000]
-    knee_growth_rate = [1E-10, 1E-7, 1E-6, 1E-4]
-    pc = ParisCurve.from_knee_points(knee_sif=knee_sif, knee_growth_rate=knee_growth_rate)
+    knee_growth_rate = [1e-10, 1e-7, 1e-6, 1e-4]
+    pc = ParisCurve.from_knee_points(
+        knee_sif=knee_sif, knee_growth_rate=knee_growth_rate
+    )
 
     assert np.allclose(pc.get_knee_sif(), knee_sif[1:-1], rtol=1e-2)
-    assert np.allclose(pc.get_knee_growth_rate(), knee_growth_rate[1:-1], rtol=1e-2)
+    assert np.allclose(
+        pc.get_knee_growth_rate(), knee_growth_rate[1:-1], rtol=1e-2
+    )
+
 
 def test_paris_curve_get_growth_rate():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    pc = ParisCurve(slope=slope, intercept=intercept, threshold=threshold, critical=critical)
+    pc = ParisCurve(
+        slope=slope, intercept=intercept, threshold=threshold, critical=critical
+    )
 
     sif_range = np.linspace(20, 2000, 10)
     growth_rates = pc.get_growth_rate(sif_range)
@@ -135,12 +166,15 @@ def test_paris_curve_get_growth_rate():
     assert len(growth_rates) == len(sif_range)
     assert np.all(growth_rates >= 0)
 
+
 def test_paris_curve_get_sif():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    pc = ParisCurve(slope=slope, intercept=intercept, threshold=threshold, critical=critical)
+    pc = ParisCurve(
+        slope=slope, intercept=intercept, threshold=threshold, critical=critical
+    )
 
     growth_rate = np.logspace(-10, -4, 10)
     sif = pc.get_sif(growth_rate)
@@ -149,20 +183,25 @@ def test_paris_curve_get_sif():
     assert np.all(sif >= threshold)
     assert np.all(sif <= critical)
 
+
 def test_paris_curve_plot():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    pc = ParisCurve(slope=slope, intercept=intercept, threshold=threshold, critical=critical)
+    pc = ParisCurve(
+        slope=slope, intercept=intercept, threshold=threshold, critical=critical
+    )
 
     fig, ax = pc.plot()
     assert fig is not None
     assert ax is not None
 
+
 def test_pure_walker_curve():
-    walker_pure = WalkerCurve(slope=3.1, intercept=1.2e-14, load_ratio=0.1,
-                              walker_exponent=0.5)
+    walker_pure = WalkerCurve(
+        slope=3.1, intercept=1.2e-14, load_ratio=0.1, walker_exponent=0.5
+    )
     assert walker_pure.linear
     assert walker_pure.slope == 3.1
     assert walker_pure.intercept == 1.2e-14
@@ -170,7 +209,9 @@ def test_pure_walker_curve():
     assert walker_pure.critical == np.inf
     assert walker_pure.walker_exponent == 0.5
     assert walker_pure.load_ratio == 0.1
-    assert walker_pure.threshold_growth_rate == walker_pure.get_growth_rate(walker_pure.threshold)
+    assert walker_pure.threshold_growth_rate == walker_pure.get_growth_rate(
+        walker_pure.threshold
+    )
     assert np.isinf(walker_pure.critical_growth_rate)
     assert np.isinf(walker_pure.get_growth_rate(walker_pure.critical))
     assert len(walker_pure.get_knee_growth_rate()) == 0
@@ -189,13 +230,49 @@ def test_pure_walker_curve():
         walker_pure.get_knee_sif(check_knee=np.array([0]))
 
 
+def test_paris_curve_accepts_walker_parameters():
+    slope = np.array([3.1])
+    intercept = np.array([1.2e-14])
+
+    curve = ParisCurve(
+        slope=slope,
+        intercept=intercept,
+        load_ratio=0.1,
+        walker_exponent=0.5,
+    )
+
+    assert curve.load_ratio == 0
+    assert curve.walker_exponent == 0
+    np.testing.assert_allclose(curve.walker_intercept, intercept)
+
+
+def test_paris_curve_plot_filters_knees_above_critical():
+    curve = ParisCurve(
+        slope=[3.0, 4.0, 5.0],
+        intercept=[1e-12, 1e-15, 1e-20],
+        threshold=10.0,
+        critical=50.0,
+    )
+
+    fig, ax = curve.plot()
+
+    assert fig is not None
+    assert ax is not None
+
+
 def test_walker_curve_initialization():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    wc = WalkerCurve(slope=slope, intercept=intercept, threshold=threshold,
-                     critical=critical, walker_exponent=0.5, load_ratio=0.1)
+    wc = WalkerCurve(
+        slope=slope,
+        intercept=intercept,
+        threshold=threshold,
+        critical=critical,
+        walker_exponent=0.5,
+        load_ratio=0.1,
+    )
 
     assert not wc.linear
     assert wc.__str__() == wc.name
@@ -234,55 +311,87 @@ def test_walker_curve_initialization():
         _calc_growth_rate(
             knee_sif, wc.slope, wc.walker_intercept, wc.threshold, wc.critical
         ),
-        rtol=1e-2
+        rtol=1e-2,
     )
     assert np.allclose(
         wc.get_sif(knee_gr),
         _calc_sif(
             knee_gr, wc.slope, wc.walker_intercept, wc.threshold, wc.critical
         ),
-        rtol=1e-2
+        rtol=1e-2,
     )
 
     with pytest.raises(TypingError):
-        assert _calc_growth_rate(
-            wc.threshold - 1e-2 , wc.slope, wc.intercept, wc.threshold, wc.critical
-        ) == 0
-    assert _calc_growth_rate(
-        np.array([wc.threshold - 1e-2]),
-        wc.slope,
-        wc.intercept,
-        wc.threshold,
-        wc.critical
-    ) == 0
+        assert (
+            _calc_growth_rate(
+                wc.threshold - 1e-2,
+                wc.slope,
+                wc.intercept,
+                wc.threshold,
+                wc.critical,
+            )
+            == 0
+        )
+    assert (
+        _calc_growth_rate(
+            np.array([wc.threshold - 1e-2]),
+            wc.slope,
+            wc.intercept,
+            wc.threshold,
+            wc.critical,
+        )
+        == 0
+    )
 
     with pytest.raises(TypingError):
-        assert _calc_growth_rate(
-            critical + 1e-2, wc.slope, wc.intercept, wc.threshold, wc.critical
-        ) == np.inf
-    assert _calc_growth_rate(
-        np.array([critical + 1e-2]),
-        wc.slope,
-        wc.intercept,
-        wc.threshold,
-        wc.critical
-    ) == np.inf
+        assert (
+            _calc_growth_rate(
+                critical + 1e-2,
+                wc.slope,
+                wc.intercept,
+                wc.threshold,
+                wc.critical,
+            )
+            == np.inf
+        )
+    assert (
+        _calc_growth_rate(
+            np.array([critical + 1e-2]),
+            wc.slope,
+            wc.intercept,
+            wc.threshold,
+            wc.critical,
+        )
+        == np.inf
+    )
+
 
 def test_walker_curve_from_knee_points():
     knee_sif = [20, 100, 500, 1000]
-    knee_growth_rate = [1E-10, 1E-7, 1E-6, 1E-4]
-    wc = WalkerCurve.from_knee_points(knee_sif=knee_sif, knee_growth_rate=knee_growth_rate)
+    knee_growth_rate = [1e-10, 1e-7, 1e-6, 1e-4]
+    wc = WalkerCurve.from_knee_points(
+        knee_sif=knee_sif, knee_growth_rate=knee_growth_rate
+    )
 
     assert np.allclose(wc.get_knee_sif(), knee_sif[1:-1], rtol=1e-2)
-    assert np.allclose(wc.get_knee_growth_rate(), knee_growth_rate[1:-1], rtol=1e-2)
+    assert np.allclose(
+        wc.get_knee_growth_rate(), knee_growth_rate[1:-1], rtol=1e-2
+    )
+
 
 def test_walker_curve_get_growth_rate():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    wc = WalkerCurve(slope=slope, intercept=intercept, threshold=threshold,
-                     critical=critical, walker_exponent=0.5, load_ratio=0.2)
+    wc = WalkerCurve(
+        slope=slope,
+        intercept=intercept,
+        threshold=threshold,
+        critical=critical,
+        walker_exponent=0.5,
+        load_ratio=0.2,
+    )
 
     sif_range = np.linspace(20, 2000, 10)
     growth_rates = wc.get_growth_rate(sif_range)
@@ -290,13 +399,20 @@ def test_walker_curve_get_growth_rate():
     assert len(growth_rates) == len(sif_range)
     assert np.all(growth_rates >= 0)
 
+
 def test_walker_curve_get_sif():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    wc = WalkerCurve(slope=slope, intercept=intercept, threshold=threshold,
-                     critical=critical, walker_exponent=0.4, load_ratio=0.3)
+    wc = WalkerCurve(
+        slope=slope,
+        intercept=intercept,
+        threshold=threshold,
+        critical=critical,
+        walker_exponent=0.4,
+        load_ratio=0.3,
+    )
 
     growth_rate = np.logspace(-10, -4, 10)
     sif = wc.get_sif(growth_rate)
@@ -305,66 +421,96 @@ def test_walker_curve_get_sif():
     assert np.all(sif >= threshold)
     assert np.all(sif <= critical)
 
+
 def test_walker_curve_plot():
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
-    wc = WalkerCurve(slope=slope, intercept=intercept, threshold=threshold,
-                     critical=critical, walker_exponent=0.5, load_ratio=0.1)
+    wc = WalkerCurve(
+        slope=slope,
+        intercept=intercept,
+        threshold=threshold,
+        critical=critical,
+        walker_exponent=0.5,
+        load_ratio=0.1,
+    )
 
     fig, ax = wc.plot()
     assert fig is not None
     assert ax is not None
 
+
 def test_walker_curve_multiple_load_ratios():
     """Assert that increasing the load ratios reduces the threshold and critical
     SIF"""
 
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
     wc_old = None
     w_exp = 0.5
     for r in [0, 0.1, 0.3, 0.5, 0.7, 0.9]:
-        wc = WalkerCurve(slope=slope, intercept=intercept, threshold=threshold,
-                         critical=critical, walker_exponent=w_exp, load_ratio=r)
+        wc = WalkerCurve(
+            slope=slope,
+            intercept=intercept,
+            threshold=threshold,
+            critical=critical,
+            walker_exponent=w_exp,
+            load_ratio=r,
+        )
         if wc_old is not None:
             assert np.all(wc.walker_intercept >= wc_old.walker_intercept)
             assert wc.critical <= wc_old.critical
         wc_old = wc
 
+
 def test_walker_curve_change_load_ratios():
     """Assert that increasing the load ratios reduces the threshold and critical
     SIF"""
     import copy
+
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
     wc_old = None
     w_exp = 0.5
-    wc = WalkerCurve(slope=slope, intercept=intercept, threshold=threshold,
-                     critical=critical, walker_exponent=w_exp, load_ratio=0)
+    wc = WalkerCurve(
+        slope=slope,
+        intercept=intercept,
+        threshold=threshold,
+        critical=critical,
+        walker_exponent=w_exp,
+        load_ratio=0,
+    )
     for r in [0.1, 0.3, 0.5, 0.7, 0.9]:
         wc_old = copy.copy(wc)
         wc.load_ratio = r
         assert np.all(wc.walker_intercept >= wc_old.walker_intercept)
         assert wc.critical <= wc_old.critical
 
+
 def test_walker_curve_change_exponent():
     """Assert that increasing the Walker exponent towards 1 moves the curve
     closer to the Paris curve"""
     import copy
+
     slope = [2.88, 5.1, 8.16, 5.1, 2.88]
-    intercept = [1E-16, 1E-20, 1E-27, 1E-19, 1E-13]
+    intercept = [1e-16, 1e-20, 1e-27, 1e-19, 1e-13]
     threshold = 20
     critical = 2000
     wc_old = None
     w_exp = 0.5
-    wc = WalkerCurve(slope=slope, intercept=intercept, threshold=threshold,
-                     critical=critical, walker_exponent=w_exp, load_ratio=0)
+    wc = WalkerCurve(
+        slope=slope,
+        intercept=intercept,
+        threshold=threshold,
+        critical=critical,
+        walker_exponent=w_exp,
+        load_ratio=0,
+    )
     for w_e in [0, 0.25, 0.5, 0.75, 1]:
         wc_old = copy.copy(wc)
         wc.walker_exponent = w_e
@@ -372,15 +518,17 @@ def test_walker_curve_change_exponent():
         assert wc.critical >= wc_old.critical
         assert wc.threshold >= wc_old.threshold
 
-    pc = ParisCurve(slope=slope, intercept=intercept, threshold=threshold,
-                    critical=critical)
+    pc = ParisCurve(
+        slope=slope, intercept=intercept, threshold=threshold, critical=critical
+    )
 
     # Assert that the two curves are the same when the exponent is 1
     assert wc.threshold == pc.threshold
     assert wc.critical == pc.critical
     assert np.allclose(wc.get_knee_sif(), pc.get_knee_sif(), rtol=1e-2)
-    assert np.allclose(wc.get_knee_growth_rate(), pc.get_knee_growth_rate(),
-                       rtol=1e-2)
+    assert np.allclose(
+        wc.get_knee_growth_rate(), pc.get_knee_growth_rate(), rtol=1e-2
+    )
 
 
 nb.config.DISABLE_JIT = False
diff --git a/tests/mean_stress/test_corrections.py b/tests/mean_stress/test_corrections.py
index 3629fbd..9913e2e 100644
--- a/tests/mean_stress/test_corrections.py
+++ b/tests/mean_stress/test_corrections.py
@@ -86,19 +86,31 @@ def test_invalid_detail_factor(self):
 
     @given(
         mean_stress=hy.lists(
-            hy.floats(min_value=-1000, max_value=1000, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=20
+            hy.floats(
+                min_value=-1000,
+                max_value=1000,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=20,
         ),
         stress_amp=hy.lists(
-            hy.floats(min_value=1, max_value=1000, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=20
-        )
+            hy.floats(
+                min_value=1,
+                max_value=1000,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=20,
+        ),
     )
     def test_property_based_dnvgl_08(self, mean_stress, stress_amp):
         """Property-based test for DNVGL with detail_factor=0.8."""
         if len(mean_stress) != len(stress_amp):
-            mean_stress = mean_stress[:min(len(mean_stress), len(stress_amp))]
-            stress_amp = stress_amp[:min(len(mean_stress), len(stress_amp))]
+            mean_stress = mean_stress[: min(len(mean_stress), len(stress_amp))]
+            stress_amp = stress_amp[: min(len(mean_stress), len(stress_amp))]
 
         mean_stress = np.array(mean_stress)
         stress_amp = np.array(stress_amp)
@@ -120,19 +132,31 @@ def test_property_based_dnvgl_08(self, mean_stress, stress_amp):
 
     @given(
         mean_stress=hy.lists(
-            hy.floats(min_value=-1000, max_value=1000, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=20
+            hy.floats(
+                min_value=-1000,
+                max_value=1000,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=20,
         ),
         stress_amp=hy.lists(
-            hy.floats(min_value=1, max_value=1000, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=20
-        )
+            hy.floats(
+                min_value=1,
+                max_value=1000,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=20,
+        ),
     )
     def test_property_based_dnvgl_06(self, mean_stress, stress_amp):
         """Property-based test for DNVGL with detail_factor=0.6."""
         if len(mean_stress) != len(stress_amp):
-            mean_stress = mean_stress[:min(len(mean_stress), len(stress_amp))]
-            stress_amp = stress_amp[:min(len(mean_stress), len(stress_amp))]
+            mean_stress = mean_stress[: min(len(mean_stress), len(stress_amp))]
+            stress_amp = stress_amp[: min(len(mean_stress), len(stress_amp))]
 
         mean_stress = np.array(mean_stress)
         stress_amp = np.array(stress_amp)
@@ -166,7 +190,9 @@ def test_edge_cases(self):
         result = dnvgl_mean_stress_correction(
             np.array([0]), np.array([50]), detail_factor=0.8
         )
-        assert result[0] < 100  # Should be less than 2*amplitude for mixed loading
+        assert (
+            result[0] < 100
+        )  # Should be less than 2*amplitude for mixed loading
         assert result[0] >= 80  # But at least detail_factor * 2 * amplitude
 
         # Very small amplitude with high tensile mean
@@ -218,7 +244,9 @@ def test_gamma_variations(self):
         np.testing.assert_allclose(result_1, stress_amplitude, rtol=1e-10)
 
     @given(
-        gamma=hy.floats(min_value=0.0, max_value=1.0, allow_nan=False, allow_infinity=False)
+        gamma=hy.floats(
+            min_value=0.0, max_value=1.0, allow_nan=False, allow_infinity=False
+        )
     )
     def test_gamma_range(self, gamma):
         """Test gamma parameter in valid range."""
@@ -258,19 +286,31 @@ def test_negative_max_stress_handling(self):
 
     @given(
         mean_stress=hy.lists(
-            hy.floats(min_value=-500, max_value=500, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=10
+            hy.floats(
+                min_value=-500,
+                max_value=500,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=10,
         ),
         stress_amp=hy.lists(
-            hy.floats(min_value=1, max_value=500, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=10
-        )
+            hy.floats(
+                min_value=1,
+                max_value=500,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=10,
+        ),
     )
     def test_property_based_walker(self, mean_stress, stress_amp):
         """Property-based test for Walker correction."""
         if len(mean_stress) != len(stress_amp):
-            mean_stress = mean_stress[:min(len(mean_stress), len(stress_amp))]
-            stress_amp = stress_amp[:min(len(mean_stress), len(stress_amp))]
+            mean_stress = mean_stress[: min(len(mean_stress), len(stress_amp))]
+            stress_amp = stress_amp[: min(len(mean_stress), len(stress_amp))]
 
         mean_stress = np.array(mean_stress)
         stress_amp = np.array(stress_amp)
@@ -322,19 +362,31 @@ def test_plotting_functionality(self):
 
     @given(
         mean_stress=hy.lists(
-            hy.floats(min_value=-500, max_value=500, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=10
+            hy.floats(
+                min_value=-500,
+                max_value=500,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=10,
         ),
         stress_amp=hy.lists(
-            hy.floats(min_value=1, max_value=500, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=10
-        )
+            hy.floats(
+                min_value=1,
+                max_value=500,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=10,
+        ),
     )
     def test_property_based_swt(self, mean_stress, stress_amp):
         """Property-based test for SWT correction."""
         if len(mean_stress) != len(stress_amp):
-            mean_stress = mean_stress[:min(len(mean_stress), len(stress_amp))]
-            stress_amp = stress_amp[:min(len(mean_stress), len(stress_amp))]
+            mean_stress = mean_stress[: min(len(mean_stress), len(stress_amp))]
+            stress_amp = stress_amp[: min(len(mean_stress), len(stress_amp))]
 
         mean_stress = np.array(mean_stress)
         stress_amp = np.array(stress_amp)
@@ -362,7 +414,10 @@ def test_basic_functionality(self):
             amp_in, mean_in, r_out, ult_s, correction_exponent
         )
 
-        assert amp_out.shape == (1, 3)  # r_out as scalar creates (1, len(amp_in))
+        assert amp_out.shape == (
+            1,
+            3,
+        )  # r_out as scalar creates (1, len(amp_in))
         assert mean_out.shape == (1, 3)
         assert np.all(amp_out >= 0)
         assert np.all(np.isfinite(amp_out))
@@ -394,12 +449,38 @@ def test_r_out_minus_one_special_case(self):
         )
 
         # For r_out = -1, mean_out should be zeros
-        np.testing.assert_allclose(mean_out[0], np.zeros_like(amp_in), rtol=1e-10)
+        np.testing.assert_allclose(
+            mean_out[0], np.zeros_like(amp_in), rtol=1e-10
+        )
 
         # Analytical solution: amp_out = amp_in / (1 - (mean_in / ult_s)^n)
-        expected_amp_out = amp_in / (1 - (mean_in / ult_s) ** correction_exponent)
+        expected_amp_out = amp_in / (
+            1 - (mean_in / ult_s) ** correction_exponent
+        )
         np.testing.assert_allclose(amp_out[0], expected_amp_out, rtol=1e-6)
 
+    def test_goodman_non_reversed_load_ratio_analytical_solution(self):
+        """Test Goodman correction with non-fully-reversed output ratios."""
+        amp_in = np.array([80.7335222, 189.5302955, 258.9631694])
+        mean_in = np.array([13.33359442, 122.13036772, 52.69749382])
+        r_out = np.array([-3.0, 0.0])
+        ult_s = 900.0
+
+        amp_out, mean_out = goodman_haigh_mean_stress_correction(
+            amp_in, mean_in, r_out, ult_s, correction_exponent=1.0
+        )
+
+        r_in = (mean_in - amp_in) / (mean_in + amp_in)
+        for idx, r_out_val in enumerate(r_out):
+            expected_amp = (
+                (1 - ((1 + r_in) / (1 - r_in) * amp_in / ult_s)) / amp_in
+                + ((1 + r_out_val) / (1 - r_out_val) / ult_s)
+            ) ** -1
+            expected_mean = expected_amp * (1 + r_out_val) / (1 - r_out_val)
+            np.testing.assert_allclose(amp_out[idx], expected_amp, rtol=1e-12)
+            np.testing.assert_allclose(mean_out[idx], expected_mean, rtol=1e-12)
+            assert np.all(amp_out[idx] >= 0)
+
     def test_different_correction_exponents(self):
         """Test different correction exponents (Goodman=1, Gerber=2)."""
         amp_in = np.array([100])
@@ -514,25 +595,47 @@ def test_initial_guess_parameter(self):
         initial_guess = np.array([120, 180])
 
         amp_out, mean_out = goodman_haigh_mean_stress_correction(
-            amp_in, mean_in, r_out, ult_s, correction_exponent,
-            initial_guess=initial_guess
+            amp_in,
+            mean_in,
+            r_out,
+            ult_s,
+            correction_exponent,
+            initial_guess=initial_guess,
         )
 
         assert amp_out.shape == (1, 2)
 
     @given(
         amp_in=hy.lists(
-            hy.floats(min_value=10, max_value=200, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=3
+            hy.floats(
+                min_value=10,
+                max_value=200,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=3,
         ),
         mean_in=hy.lists(
-            hy.floats(min_value=-50, max_value=50, allow_nan=False, allow_infinity=False),
-            min_size=1, max_size=3
+            hy.floats(
+                min_value=-50,
+                max_value=50,
+                allow_nan=False,
+                allow_infinity=False,
+            ),
+            min_size=1,
+            max_size=3,
+        ),
+        ult_s=hy.floats(
+            min_value=800, max_value=2000, allow_nan=False, allow_infinity=False
+        ),
+        correction_exponent=hy.floats(
+            min_value=1.0, max_value=2.0, allow_nan=False, allow_infinity=False
         ),
-        ult_s=hy.floats(min_value=800, max_value=2000, allow_nan=False, allow_infinity=False),
-        correction_exponent=hy.floats(min_value=1.0, max_value=2.0, allow_nan=False, allow_infinity=False)
     )
-    def test_property_based_goodman_haigh(self, amp_in, mean_in, ult_s, correction_exponent):
+    def test_property_based_goodman_haigh(
+        self, amp_in, mean_in, ult_s, correction_exponent
+    ):
         """Property-based test for Goodman-Haigh correction."""
         if len(amp_in) != len(mean_in):
             min_len = min(len(amp_in), len(mean_in))
@@ -605,7 +708,11 @@ def test_edge_case_stress_combinations(self):
         zero_mean = np.array([0])
         small_amp = np.array([1e-6])
 
-        for correction in [dnvgl_mean_stress_correction, walker_mean_stress_correction, swt_mean_stress_correction]:
+        for correction in [
+            dnvgl_mean_stress_correction,
+            walker_mean_stress_correction,
+            swt_mean_stress_correction,
+        ]:
             if correction == dnvgl_mean_stress_correction:
                 result = correction(zero_mean, small_amp, detail_factor=0.8)
             else:
@@ -625,4 +732,4 @@ def test_large_stress_values(self):
 
         for result in [result_dnvgl, result_walker, result_swt]:
             assert np.all(np.isfinite(result))
-            assert np.all(result >= 0)
\ No newline at end of file
+            assert np.all(result >= 0)
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 35990d8..4ec4259 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -10,6 +10,7 @@
 from unittest.mock import Mock
 
 # Non-standard imports
+import numba as nb
 import numpy as np
 import matplotlib.pyplot as plt
 import pytest
@@ -94,6 +95,8 @@ def test_split_full_cycles_and_residuals(
     (0, 1, 99, *expected(99.0)),
     (0, 1, 98.9, *expected(99.0)),
 ]
+
+
 @pytest.mark.parametrize(bub_names, bub_data)
 def test_bin_upper_bound(
     bin_lower_bound: float,
@@ -172,6 +175,8 @@ def test_calc_half_cycles(
 
 
 cfc_names = "count_cycle, expected, error"
+
+
 # fmt: off
 cfc_data = [
     (COUNTS_1, *expected(np.array(
@@ -231,27 +236,19 @@ def test_faulty_fatigue_stress(self, lists: Tuple[List[int], List[int]]):
         with self.assertRaises(ValueError) as ve:
             pu.FatigueStress(_counts=cts, _values=vals[1:], bin_width=0.5)
             self.assertEqual(
-                "counts and values must have the same length",
-                str(ve.exception)
+                "counts and values must have the same length", str(ve.exception)
             )
 
     @given(same_len_lists())
     def test_empty_fatigue_stress(self, lists: Tuple[List[int], List[int]]):
-        """Assert that an error is raised when the counts or values is empty
-        """
+        """Assert that an error is raised when the counts or values is empty"""
         vals, cts = lists
         with self.assertRaises(ValueError) as ve:
             pu.FatigueStress(_counts=cts, _values=[], bin_width=0.5)
-            self.assertEqual(
-                "No data provided",
-                str(ve.exception)
-            )
+            self.assertEqual("No data provided", str(ve.exception))
         with self.assertRaises(ValueError) as ve:
             pu.FatigueStress(_counts=[], _values=vals, bin_width=0.5)
-            self.assertEqual(
-                "No data provided",
-                str(ve.exception)
-            )
+            self.assertEqual("No data provided", str(ve.exception))
 
 
 @pytest.mark.parametrize(
@@ -377,18 +374,17 @@ def test_make_axes():
     plt.close(fig)
 
     # Test with invalid figure type
-    with pytest.raises(TypeError, match="fig must be a matplotlib.figure.Figure"):
+    with pytest.raises(
+        TypeError, match="fig must be a matplotlib.figure.Figure"
+    ):
         pu.make_axes(fig="not_a_figure")
 
 
 def test_to_numba_dict():
     """Test to_numba_dict function"""
     import numba as nb
-    data = {
-        "key1": 1.0,
-        "key2": 2.0,
-        "key3": 3.0
-    }
+
+    data = {"key1": 1.0, "key2": 2.0, "key3": 3.0}
 
     result = pu.to_numba_dict(data)
     # The next line fails since os.environ["NUMBA_DISABLE_JIT"] = "1"
@@ -559,6 +555,7 @@ def test_fatigue_stress_setters():
 
 def test_py_bisect():
     """Test py_bisect function"""
+
     def test_func(x):
         return x**2 - 4
 
@@ -575,6 +572,7 @@ def simple_func(x):
 
 def test_py_newton():
     """Test py_newton function"""
+
     def test_func(x):
         return x**2 - 4
 
@@ -584,6 +582,7 @@ def test_func(x):
 
 def test_numba_bisect():
     """Test numba_bisect wrapper"""
+
     def test_func(x):
         return x**2 - 4
 
@@ -593,6 +592,7 @@ def test_func(x):
 
 def test_numba_newton():
     """Test numba_newton wrapper"""
+
     def test_func(x):
         return x**2 - 4
 
@@ -602,6 +602,7 @@ def test_func(x):
 
 def test_compile_specialized_bisect():
     """Test compile_specialized_bisect function"""
+
     def test_func(x):
         return (x**2 - 4,)
 
@@ -610,8 +611,105 @@ def test_func(x):
     assert np.isclose(root, 2.0, atol=1e-5)
 
 
+def test_compile_specialized_bisect_many_args():
+    """Test specialized bisection with more than three extra arguments."""
+
+    def test_func(x, arg_1, arg_2, arg_3, arg_4, arg_5):
+        return x - (arg_1 + arg_2 + arg_3 + arg_4 + arg_5)
+
+    compiled_bisect = pu.compile_specialized_bisect(test_func)
+    root = compiled_bisect(
+        0,
+        20,
+        1e-6,
+        100,
+        1.0,
+        2.0,
+        3.0,
+        4.0,
+        5.0,
+    )
+    assert np.isclose(root, 15.0, atol=1e-5)
+
+
+@pytest.mark.parametrize(
+    "test_func",
+    [
+        lambda x: x**2 - 4,
+        lambda x: (x**2 - 4,),
+        lambda x: np.array([x**2 - 4]),
+        lambda x: [x**2 - 4],
+        lambda x: np.float64(x**2 - 4),
+    ],
+)
+def test_compile_specialized_bisect_scalar_like_returns(test_func):
+    """Test bisection residuals that return scalar-like values."""
+
+    compiled_bisect = pu.compile_specialized_bisect(test_func)
+    root = compiled_bisect(0, 5, tol=1e-6, mxiter=100)
+    assert np.isclose(root, 2.0, atol=1e-5)
+
+
+@pytest.mark.parametrize(
+    "test_func",
+    [
+        lambda x: (x**2 - 4, x),
+        lambda x: np.array([x**2 - 4, x]),
+        lambda x: [x**2 - 4, x],
+    ],
+)
+def test_compile_specialized_bisect_rejects_non_scalar_returns(test_func):
+    """Test bisection rejects residuals that return multiple values."""
+
+    compiled_bisect = pu.compile_specialized_bisect(test_func)
+    with pytest.raises((TypeError, nb.core.errors.TypingError)):
+        compiled_bisect(0, 5, tol=1e-6, mxiter=100)
+
+
+def test_compile_specialized_bisect_uses_numba_dispatcher(monkeypatch):
+    """Test bisection uses a numba dispatcher when JIT is enabled."""
+
+    monkeypatch.delenv("NUMBA_DISABLE_JIT", raising=False)
+    monkeypatch.setattr(pu.nb.config, "DISABLE_JIT", False)
+    pu.compile_specialized_bisect.cache_clear()
+
+    def test_func(x):
+        return (x**2 - 4,)
+
+    compiled_bisect = pu.compile_specialized_bisect(test_func)
+    root = compiled_bisect(0, 5, tol=1e-6, mxiter=100)
+
+    assert compiled_bisect.__class__.__name__ == "CPUDispatcher"
+    assert np.isclose(root, 2.0, atol=1e-5)
+
+    pu.compile_specialized_bisect.cache_clear()
+
+
+def test_numba_bisect_accepts_specialized_function():
+    """Test numba_bisect does not recompile specialized bisection functions."""
+
+    def test_func(x, arg_1, arg_2, arg_3, arg_4, arg_5):
+        return x - (arg_1 + arg_2 + arg_3 + arg_4 + arg_5)
+
+    compiled_bisect = pu.compile_specialized_bisect(test_func)
+    root = pu.numba_bisect(
+        compiled_bisect,
+        0,
+        20,
+        1e-6,
+        100,
+        1.0,
+        2.0,
+        3.0,
+        4.0,
+        5.0,
+    )
+    assert np.isclose(root, 15.0, atol=1e-5)
+
+
 def test_compile_specialized_newton():
     """Test compile_specialized_newton function"""
+
     def test_func(x):
         return x**2 - 4
 
@@ -620,6 +718,12 @@ def test_func(x):
     assert np.isclose(root, 2.0, atol=1e-5)
 
 
+def test_warmup_numba():
+    """Test warmup_numba compiles supported numba call paths."""
+
+    pu.warmup_numba()
+
+
 def test_custom_formatter():
     """Test CustomFormatter class"""
     formatter = pu.CustomFormatter()
@@ -629,8 +733,13 @@ def test_custom_formatter():
 
     # Test DEBUG level
     record_debug = logging.LogRecord(
-        name="test", level=logging.DEBUG, pathname="test.py",
-        lineno=1, msg="Debug message", args=(), exc_info=None
+        name="test",
+        level=logging.DEBUG,
+        pathname="test.py",
+        lineno=1,
+        msg="Debug message",
+        args=(),
+        exc_info=None,
     )
     formatted_debug = formatter.format(record_debug)
     assert "Debug message" in formatted_debug
@@ -638,8 +747,13 @@ def test_custom_formatter():
 
     # Test INFO level
     record_info = logging.LogRecord(
-        name="test", level=logging.INFO, pathname="test.py",
-        lineno=1, msg="Info message", args=(), exc_info=None
+        name="test",
+        level=logging.INFO,
+        pathname="test.py",
+        lineno=1,
+        msg="Info message",
+        args=(),
+        exc_info=None,
     )
     formatted_info = formatter.format(record_info)
     assert "Info message" in formatted_info
@@ -647,8 +761,13 @@ def test_custom_formatter():
 
     # Test WARNING level
     record_warning = logging.LogRecord(
-        name="test", level=logging.WARNING, pathname="test.py",
-        lineno=1, msg="Warning message", args=(), exc_info=None
+        name="test",
+        level=logging.WARNING,
+        pathname="test.py",
+        lineno=1,
+        msg="Warning message",
+        args=(),
+        exc_info=None,
     )
     formatted_warning = formatter.format(record_warning)
     assert "Warning message" in formatted_warning
@@ -656,8 +775,13 @@ def test_custom_formatter():
 
     # Test ERROR level
     record_error = logging.LogRecord(
-        name="test", level=logging.ERROR, pathname="test.py",
-        lineno=1, msg="Error message", args=(), exc_info=None
+        name="test",
+        level=logging.ERROR,
+        pathname="test.py",
+        lineno=1,
+        msg="Error message",
+        args=(),
+        exc_info=None,
     )
     formatted_error = formatter.format(record_error)
     assert "Error message" in formatted_error
@@ -665,8 +789,13 @@ def test_custom_formatter():
 
     # Test CRITICAL level
     record_critical = logging.LogRecord(
-        name="test", level=logging.CRITICAL, pathname="test.py",
-        lineno=1, msg="Critical message", args=(), exc_info=None
+        name="test",
+        level=logging.CRITICAL,
+        pathname="test.py",
+        lineno=1,
+        msg="Critical message",
+        args=(),
+        exc_info=None,
     )
     formatted_critical = formatter.format(record_critical)
     assert "Critical message" in formatted_critical
@@ -789,7 +918,9 @@ def test_fatigue_stress_bin_properties():
     """Test FatigueStress bin-related properties"""
     counts = np.array([1.0, 2.0, 3.0])
     values = np.array([10.0, 20.0, 30.0])
-    fs = pu.FatigueStress(_counts=counts, _values=values, bin_width=5.0, _bin_lb=5.0, _bin_ub=35.0)
+    fs = pu.FatigueStress(
+        _counts=counts, _values=values, bin_width=5.0, _bin_lb=5.0, _bin_ub=35.0
+    )
 
     # Test bin_edges
     edges = fs.bin_edges
@@ -820,7 +951,9 @@ def test_plot_damage_accumulation():
     cumsum_pm_dmg = np.array([0.05, 0.2, 0.5, 0.8, 1.1])
     limit_damage = 1.0
 
-    fig, ax = pu._plot_damage_accumulation(cumsum_nl_dmg, cumsum_pm_dmg, limit_damage)
+    fig, ax = pu._plot_damage_accumulation(
+        cumsum_nl_dmg, cumsum_pm_dmg, limit_damage
+    )
     assert isinstance(fig, plt.Figure)
     assert isinstance(ax, plt.Axes)
     plt.close(fig)
@@ -895,5 +1028,7 @@ def test_json_encoders():
 
     # Test in actual JSON encoding
     data = {"array": arr}
-    json_str = json.dumps(data, default=lambda x: pu.JSON_ENCODERS.get(type(x), str)(x))
-    assert "[1, 2, 3]" in json_str
\ No newline at end of file
+    json_str = json.dumps(
+        data, default=lambda x: pu.JSON_ENCODERS.get(type(x), str)(x)
+    )
+    assert "[1, 2, 3]" in json_str

From c3083abb08b4dfd9398201d2624970f9fe463066 Mon Sep 17 00:00:00 2001
From: Pietro D'Antuono <pietro.dantuono@24sea.eu>
Date: Mon, 8 Jun 2026 17:01:48 +0200
Subject: [PATCH 2/4] perf(#23): Optimize Numba + Python 3.14 compat

---
 numba-speedup.md | 152 -----------------------------------------------
 1 file changed, 152 deletions(-)
 delete mode 100644 numba-speedup.md

diff --git a/numba-speedup.md b/numba-speedup.md
deleted file mode 100644
index 4971d7f..0000000
--- a/numba-speedup.md
+++ /dev/null
@@ -1,152 +0,0 @@
-# Numba speedup investigation
-
-This is an explanation/reference report for maintainers of `py_fatigue`.
-The goal is to identify where numba compilation is paid repeatedly and where
-small refactors could improve performance without changing fatigue logic.
-
-## Package overview
-
-`py_fatigue` provides tools for fatigue analysis:
-
-- `cycle_count`: rainflow counting, turning points, histograms, and cycle-count
-  objects.
-- `material`: S-N curves and Paris/Walker crack-growth curves.
-- `damage`: Palmgren-Miner, non-linear stress-life damage, Theil stress-life
-  damage, and crack-growth integration.
-- `mean_stress`: DNVGL, Walker/SWT, and Goodman-Haigh mean-stress correction.
-- `geometry`: crack geometry models and stress-intensity geometry factors.
-- `stress_range`, `testing`, and `utils`: supporting data structures,
-  numerical helpers, plotting helpers, and root finders.
-
-The numba usage is concentrated in material curve evaluation, crack-growth
-geometry/SIF calculations, Theil damage, rainflow crossing detection, and the
-generic bisection/Newton wrappers in `py_fatigue.utils`.
-
-## Benchmark setup
-
-I added:
-
-- `scripts/benchmark_numba_speedups.py`: rich CLI benchmark and numba inventory.
-- `scripts/numba_benchmark_results.json`: raw benchmark output from the main run.
-
-The script uses `rich` and `rich-argparse`. The environment was created and
-managed with `uv`:
-
-```bash
-uv sync --no-default-groups --group dev
-uv run python scripts/benchmark_numba_speedups.py \
-  --calls 8 \
-  --repeats 3 \
-  --json-output scripts/numba_benchmark_results.json
-```
-
-The timings below are median wall-clock timings on this branch with
-Python 3.13 + Numba 0.65.1. They are intended to compare call patterns, not
-to be absolute performance claims. Compatibility validation was also run on
-Numba 0.61.2, 0.62.1, 0.64.0, and 0.65.1 on Python 3.13, plus Python 3.14
-with Numba 0.65.1.
-
-## Numba inventory
-
-| Area | Symbol | Mechanism | Assessment |
-| --- | --- | --- | --- |
-| `cycle_count.rainflow` | `_findcross` | `@njit` with explicit signature and `cache=True` | Compiles once per process/signature and can reuse numba's disk cache. |
-| `damage.crack_growth` | `CalcCrackGrowth` | `@jitclass(spec)` | Compiles on first matching instantiation. Warm calls reuse compiled methods. |
-| `damage.crack_growth` | `get_geometry_factor`, `get_sif` | `@nb.njit(..., cache=True)` | Already cached and signature-pinned. |
-| `damage.stress_life` | `_calc_theil_sn_damage` | `@nb.njit(cache=True)` | Already cached. |
-| `geometry.cylinder` | `f_hol_cyl_01` | `@nb.njit(..., cache=True)` | Already cached and signature-pinned. |
-| `material.crack_growth_curve` | `_calc_growth_rate`, `_calc_sif` | `@nb.njit(cache=True)` | Warm calls reuse compilation and can reuse numba's disk cache. |
-| `material.sn_curve` | `_calc_cycles`, `_calc_cycles_2`, `_calc_stress` | `@nb.njit(cache=True)` | Warm calls reuse compilation and can reuse numba's disk cache. `_calc_cycles_2` is used by `SNCurve.get_cycles`; `_calc_cycles` and `_calc_stress` appear to be legacy/unused alternatives. |
-| `material.sn_curve` | `_calc_stress_2` | `@nb.njit(cache=True)` | Already cached and used by `SNCurve.get_stress`. |
-| `material.sn_curve` | `__jit_sn_curve_residuals` | cached `compile_specialized_bisect(...)` | Fixed: the specialized bisection function is a numba dispatcher with a nopython bisection loop. |
-| `mean_stress.corrections` | `__jit_goodman_equation` | `compile_specialized_newton(...)` | OK: returns a numba dispatcher, so `numba_newton(...)` does not recompile on each call. |
-| `utils` | `numba_bisect`, `numba_newton` | generic wrappers with cached specialization | Fixed: repeated calls with the same Python residual reuse the same specialized root finder in-process. Bisection now executes its loop in nopython mode when JIT is enabled. |
-| `utils` / package API | `warmup_numba` | opt-in warmup function | New: compiles common numba call paths for the current process. |
-
-## Benchmark results
-
-| Case | Current pattern | Hoisted/cached pattern | Observed impact |
-| --- | ---: | ---: | ---: |
-| `numba_bisect(py_func)` | 0.002 ms/call | 0.001 ms/call | Repeated-call overhead removed; loop runs in nopython mode. |
-| `numba_newton(py_func)` | 0.002 ms/call | 0.001 ms/call | Repeated-call overhead removed. |
-| `SNCurve.get_cycles` | 0.414 ms/call | N/A | Warm dispatcher is reused. |
-| `SNCurve.get_stress` | 0.448 ms/call | N/A | Warm dispatcher is reused. |
-| `ParisCurve.get_growth_rate` | 1.273 ms/call | N/A | Warm dispatcher is reused; compatibility guards for older Numba do not affect the 0.65 path. |
-| `ParisCurve.get_sif` | 2.177 ms/call | N/A | Warm dispatcher is reused; compatibility guards for older Numba do not affect the 0.65 path. |
-| `rainflow.findcross` | 0.126 ms/call | N/A | Warm dispatcher is reused via cached signature. |
-| `rainflow.findtp` | 0.168 ms/call | N/A | Warm dispatcher is reused via `findcross`. |
-| `f_hol_cyl_01` | 0.003 ms/call | N/A | Warm dispatcher is reused and cached. |
-| `crack_growth.get_sif` | 0.006 ms/call | N/A | Warm dispatcher is reused and cached. |
-| `CalcCrackGrowth` construction | 0.717 ms/call | N/A | Jitclass is warm after first instantiation. |
-| `calc_theil_sn_damage` | 0.024 ms/call | N/A | Warm dispatcher is reused and cached. |
-| `goodman_haigh_mean_stress_correction` | 0.823 ms/call | N/A | The existing module-level Newton specialization avoids repeated compilation. |
-| `find_sn_curve_intersection` | 0.066 ms/call | N/A | Fixed; no longer fails with the bisection argument limit. |
-
-The original root-finder numbers showed the issue clearly: when a Python
-residual was sent through `numba_bisect` or `numba_newton` on every solve, the
-wrapper created a new numba dispatcher every time. The implemented cache removes
-that repeated compilation for repeated calls with the same residual function.
-
-## Findings
-
-1. The direct `@njit` functions are not recompiling on every warm call. The
-   retained direct dispatchers now use `cache=True` where practical, so they can
-   also reuse numba's disk cache across processes.
-2. `numba_newton` is safe in the Goodman-Haigh public path because
-   `__jit_goodman_equation` is a numba dispatcher, not a plain Python function.
-3. The S-N intersection path is fixed. It now calls the already-specialized
-   bisection function directly, and the bisection wrapper no longer imposes the
-   previous "three extra arguments" limit.
-4. The generic wrappers now cache specialized root finders by residual function
-   identity. Any loop that calls `numba_bisect(py_func, ...)` or
-   `numba_newton(py_func, ...)` repeatedly with the same function avoids repeated
-   compilation.
-5. The bisection algorithm itself now runs in nopython mode when JIT is enabled,
-   while preserving scalar-like residual return compatibility for scalar,
-   one-item tuple, one-item array, one-item list, and numpy scalar returns.
-6. Python 3.13 compatibility is now validated on Numba 0.61.2, 0.62.1, 0.64.0,
-   and 0.65.1. Python 3.14 is validated on Numba 0.65.1. The oldest supported
-   line (0.61) needs a small compatibility fallback for the low-level
-   Paris/Walker kernel symbols, but the public curve APIs and the targeted test
-   surface stay green.
-
-## Notebook and docs validation
-
-All tutorial notebooks in `notebooks/tutorials-repository/notebooks` now
-execute on Python 3.13 and 3.14.
-
-Key fixes:
-
-1. The `0.5.x` use-case notebooks no longer hard-fail when the private CSV data
-   is absent; they fall back to a small synthetic dataset and use Pandas 3.x
-   compatible aggregation aliases.
-2. The `0.6.x` nonlinear-damage notebooks were aligned with the current
-   `py_fatigue` API and had stale scratch appendix cells trimmed or updated so
-   automated execution is reliable.
-3. `1.0.0 - AM.General.Usage.ipynb` now uses a valid `python3` kernelspec and
-   Pandas 3.x compatible positional indexing examples.
-
-## Recommended refactors
-
-These are performance refactors only; they do not require changing the fatigue
-equations or output logic.
-
-1. Use `py_fatigue.warmup_numba()` at application startup when predictable
-   runtime latency is more important than startup time.
-2. Keep the root-finder residual functions stable. The cache is keyed by the
-   Python function object, so repeatedly creating new lambdas/functions will
-   still create new specializations.
-3. For maximum reuse of compiled root finders, define residual functions once at
-   module scope instead of creating new lambdas/functions inside tight loops.
-
-## Priority
-
-The highest-value root-finder fixes have been applied: repeated calls now reuse
-specializations, bisection executes in nopython mode when JIT is enabled, and
-`find_sn_curve_intersection` no longer fails under JIT. `warmup_numba()` provides
-opt-in in-process warmup. Full "always precompiled" behavior is only partially
-possible: signature-pinned `cache=True` dispatchers can reuse numba's disk cache
-across processes, but jitclasses and dynamic root-finder closures still need
-in-process warmup after each process starts. The remaining compatibility work is
-around old-Numpy/old-Numba constructor quirks rather than repeated-compilation
-overhead.

From 3f8322e0e5aff8d5eaa180ef5a35f8af0e707c7f Mon Sep 17 00:00:00 2001
From: Pietro D'Antuono <pietro.dantuono@24sea.eu>
Date: Tue, 9 Jun 2026 10:49:29 +0200
Subject: [PATCH 3/4] perf(#23): Optimize Numba + Python 3.14 compat

---
 py_fatigue/cycle_count/rainflow.py |  6 +++---
 py_fatigue/utils.py                | 31 +++++++++++++++++++++++++++---
 tasks/test.py                      | 23 ++++++++++++++++++----
 tests/conftest.py                  | 11 +++++++++++
 4 files changed, 61 insertions(+), 10 deletions(-)

diff --git a/py_fatigue/cycle_count/rainflow.py b/py_fatigue/cycle_count/rainflow.py
index a229f0d..cb4b73b 100644
--- a/py_fatigue/cycle_count/rainflow.py
+++ b/py_fatigue/cycle_count/rainflow.py
@@ -422,7 +422,7 @@ def findrfc_astm(tp: np.ndarray, t: Optional[np.ndarray] = None) -> np.ndarray:
 
 
 def _extract_cycles_from_sequence(
-    seq: List[Tuple[float, Optional[float], int]]
+    seq: List[Tuple[float, Optional[float], int]],
 ) -> Tuple[List[Dict[str, Any]], List[Tuple[float, Optional[float], int]]]:
     """
     Internal 4-point extraction following Amzallag et al. (1994).
@@ -474,7 +474,7 @@ def _extract_cycles_from_sequence(
 
 
 def _duplicate_and_select_crossing_cycles(
-    residue_seq: List[Tuple[float, Optional[float], int]]
+    residue_seq: List[Tuple[float, Optional[float], int]],
 ) -> List[Dict[str, Any]]:
     """
     Decompose residue by duplication and select cycles crossing boundary.
@@ -507,7 +507,7 @@ def _duplicate_and_select_crossing_cycles(
 
 
 def _process_residue(
-    residue_seq: List[Tuple[float, Optional[float], int]]
+    residue_seq: List[Tuple[float, Optional[float], int]],
 ) -> Tuple[List[Dict[str, Any]], List[Tuple[float, Optional[float], int]]]:
     """
     Fully decompose residue by iterative duplication until stable.
diff --git a/py_fatigue/utils.py b/py_fatigue/utils.py
index b87a838..4ad16f6 100644
--- a/py_fatigue/utils.py
+++ b/py_fatigue/utils.py
@@ -1190,9 +1190,34 @@ class CustomFormatter(logging.Formatter):
     FORMATS = {
         logging.DEBUG: grey + "🐞 " + level + italic + message + reset,
         logging.INFO: blue + "ℹ️ " + level + italic + message + reset,
-        logging.WARNING: yellow + "⚠️ " + level + newline + lineofile + italic + newline +message + reset,
-        logging.ERROR: red + "⛔ " + level + newline + lineofile + italic + newline +message + reset,
-        logging.CRITICAL: red + "🆘 " + level + newline + lineofile + bold + italic + newline +message + reset,
+        logging.WARNING: yellow
+        + "⚠️ "
+        + level
+        + newline
+        + lineofile
+        + italic
+        + newline
+        + message
+        + reset,
+        logging.ERROR: red
+        + "⛔ "
+        + level
+        + newline
+        + lineofile
+        + italic
+        + newline
+        + message
+        + reset,
+        logging.CRITICAL: red
+        + "🆘 "
+        + level
+        + newline
+        + lineofile
+        + bold
+        + italic
+        + newline
+        + message
+        + reset,
     }
 
     def format(self, record):
diff --git a/tasks/test.py b/tasks/test.py
index ec834c6..fa2f0d8 100644
--- a/tasks/test.py
+++ b/tasks/test.py
@@ -1,5 +1,7 @@
 """Test tasks."""  # pylint: disable=R0801
 
+import shutil
+
 from invoke import task
 
 from .colors import Color, colorize
@@ -28,6 +30,7 @@ def run(c_r, test=None, pytest_args="-v -W ignore::UserWarning"):
         test_command = f" {test}"
 
     _command = (
+        "MPLBACKEND=Agg PY_FATIGUE_TEST_NO_PLOTS=1 "
         f"pytest {pytest_args} "
         f"--cov={c_r.project_slug} --cov-report=term:skip-covered "
         f"--cov-report=html --cov-report=html:{COV_DOC_BUILD_DIR} {test_command}"
@@ -45,6 +48,16 @@ def coverage(c_r):
     COV_PORT = c_r.start_port + 2  # pylint: disable=C0103
 
     if SYSTEM in [OperatingSystem.LINUX, OperatingSystem.MAC]:
+        if shutil.which("screen") is None:
+            print(
+                colorize(
+                    "screen is not installed; skipping coverage webserver. "
+                    "Use `python -m http.server --directory "
+                    f"{COV_DOC_BUILD_DIR} {COV_PORT}` to serve manually.",
+                    color=Color.WARNING,
+                )
+            )
+            return
         _command = (
             f"screen -d -S {COV_SCREEN_NAME} "
             "-m python -m http.server --bind localhost "
@@ -86,9 +99,12 @@ def stop(c_r):
             "\nStopping coverage server...\n", color=Color.HEADER, bold=True
         )
         print(tmp_str)
-        _command = f"kill $(lsof -ti:{COV_PORT})"
+        _command = (
+            f"pids=$(lsof -ti:{COV_PORT}); "
+            'if [ -n "$pids" ]; then kill $pids; fi'
+        )
         print(f"{colorize('>>> ' + _command, color=Color.OKBLUE)}\n")
-        c_r.run(_command)
+        c_r.run(_command, warn=True)
 
     elif SYSTEM == OperatingSystem.WINDOWS:
         print(
@@ -113,7 +129,6 @@ def stop(c_r):
 def all(
     c_r, test=None, pytest_args="-v -W ignore::UserWarning"
 ):  # pylint: disable=W0622
-    """Run all tests and start coverage report webserver."""
+    """Run all tests."""
     stop(c_r)
     run(c_r, test, pytest_args)
-    coverage(c_r)
diff --git a/tests/conftest.py b/tests/conftest.py
index eb901ff..e2cc7d1 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -57,3 +57,14 @@ def datadir(tmpdir, request):
         dir_util.copy_tree(test_dir, str(tmpdir))
 
     return tmpdir
+
+
+@pytest.fixture(autouse=True)
+def disable_plot_show(monkeypatch):
+    """Disable interactive plot rendering for task-based test runs."""
+    if os.getenv("PY_FATIGUE_TEST_NO_PLOTS") != "1":
+        return
+
+    import matplotlib.pyplot as plt
+
+    monkeypatch.setattr(plt, "show", lambda *args, **kwargs: None)

From ed126084ef4374a805b5c57bd7a32f900055a5a5 Mon Sep 17 00:00:00 2001
From: Pietro D'Antuono <pietro.dantuono@24sea.eu>
Date: Tue, 9 Jun 2026 11:25:24 +0200
Subject: [PATCH 4/4] perf(#23): Optimize Numba + Python 3.14 compat

---
 py_fatigue/material/sn_curve.py | 63 +++-------------------------
 tests/material/test_sn_curve.py | 73 +++++++++++++++++++++++++++++++++
 2 files changed, 79 insertions(+), 57 deletions(-)

diff --git a/py_fatigue/material/sn_curve.py b/py_fatigue/material/sn_curve.py
index 45aa582..650c7b2 100644
--- a/py_fatigue/material/sn_curve.py
+++ b/py_fatigue/material/sn_curve.py
@@ -690,12 +690,12 @@ def get_cycles(
         #     axis=1
         # ), self.endurance)
 
-        calc_cycles = _calc_cycles_2
+        calc_cycles = _calc_cycles
         if (
             getattr(nb.config, "DISABLE_JIT", False)
             or os.environ.get("NUMBA_DISABLE_JIT") == "1"
         ):
-            calc_cycles = _calc_cycles_2.py_func
+            calc_cycles = getattr(_calc_cycles, "py_func", _calc_cycles)
         return calc_cycles(
             stress_range, self.slope, self.intercept, self.endurance
         )
@@ -718,12 +718,12 @@ def get_stress(
         #         cycles.shape[0],-1),
         #         axis=1
         # ), endurance_stress)
-        calc_stress = _calc_stress_2
+        calc_stress = _calc_stress
         if (
             getattr(nb.config, "DISABLE_JIT", False)
             or os.environ.get("NUMBA_DISABLE_JIT") == "1"
         ):
-            calc_stress = _calc_stress_2.py_func
+            calc_stress = getattr(_calc_stress, "py_func", _calc_stress)
         return calc_stress(cycles, self.slope, self.intercept, self.endurance)
 
     def n(  # pylint: disable=invalid-name
@@ -970,30 +970,6 @@ def get_sn_array_min(values):  # pragma: no cover
     # parallel=True,
 )
 def _calc_cycles(stress, slope, intercept, endurance):  # pragma: no cover
-    # pylint: disable=not-an-iterable
-    assert intercept.size > 0 and intercept.size == slope.size
-    assert get_sn_array_min(stress) >= 0
-    log10 = np.log(10)
-    the_cycles = np.empty(stress.size)
-    for i in nb.prange(stress.size):
-        log_stress = np.log(stress[i])
-        max_i = intercept[0] * log10 - slope[0] * log_stress
-        for j in range(1, len(intercept)):
-            value = intercept[j] * log10 - slope[j] * log_stress
-            max_i = max(max_i, value)
-        the_cycles[i] = np.exp(max_i)
-    if endurance < np.inf:
-        the_cycles[the_cycles > endurance] = np.inf
-    return the_cycles
-
-
-@nb.njit(
-    # 'float64[::1](float64[::1], float64[::1], float64[::1])',
-    fastmath=False,
-    cache=True,
-    # parallel=True,
-)
-def _calc_cycles_2(stress, slope, intercept, endurance):  # pragma: no cover
     """
     Calculate the number of cycles to failure for given stress levels.
 
@@ -1049,33 +1025,6 @@ def _calc_cycles_2(stress, slope, intercept, endurance):  # pragma: no cover
     # parallel=True,
 )
 def _calc_stress(cycles, slope, intercept, endurance):  # pragma: no cover
-    # pylint: disable=not-an-iterable
-    assert intercept.size > 0 and intercept.size == slope.size
-    assert get_sn_array_min(cycles) > 0
-    log10 = np.log(10)
-    the_stress = np.empty(cycles.size)
-    for i in nb.prange(cycles.size):
-        log_cycles = np.log(cycles[i])
-        max_i = (intercept[0] * log10 - log_cycles) / slope[0]
-        for j in range(1, len(intercept)):
-            value = (intercept[j] * log10 - log_cycles) / slope[j]
-            max_i = max(max_i, value)
-        the_stress[i] = np.exp(max_i)
-    if endurance < np.inf:
-        endurance_stress = np.exp(
-            (intercept[-1] * log10 - np.log(endurance)) / slope[-1]
-        )
-        the_stress[the_stress < endurance_stress] = endurance_stress
-    return the_stress
-
-
-@nb.njit(
-    # 'float64[::1](float64[::1], float64[::1], float64[::1])',
-    fastmath=False,
-    # parallel=True,
-    cache=True,
-)
-def _calc_stress_2(cycles, slope, intercept, endurance):  # pragma: no cover
     """
     Calculate the number of cycles to failure for given stress levels.
 
@@ -1089,7 +1038,7 @@ def _calc_stress_2(cycles, slope, intercept, endurance):  # pragma: no cover
     numpy.ndarray: Array of calculated cycles to failure.
     """
     assert intercept.size > 0 and intercept.size == slope.size
-    assert get_sn_array_min(cycles) >= 0
+    assert get_sn_array_min(cycles) > 0
 
     log_cycles = np.log10(cycles)
     log_endurance = np.log10(endurance)
@@ -1175,7 +1124,7 @@ def __sn_curve_residuals(  # pragma: no cover
     #     res_stress = 0.
     # print(f"np.array([cycles]) = {cycles}, type = {type(cycles)}")
 
-    fail = _calc_stress_2(np.array([cycles]), slope, intercept, endurance)[0]
+    fail = _calc_stress(np.array([cycles]), slope, intercept, endurance)[0]
     return fail - weight * cycles - res_stress
 
 
diff --git a/tests/material/test_sn_curve.py b/tests/material/test_sn_curve.py
index 57f2390..2cf3af4 100644
--- a/tests/material/test_sn_curve.py
+++ b/tests/material/test_sn_curve.py
@@ -5,6 +5,7 @@
 import numba as nb
 import numpy as np
 import pytest
+from py_fatigue.material.sn_curve import _calc_cycles, _calc_stress
 
 # os.environ["NUMBA_DISABLE_JIT"] = "1"
 
@@ -462,4 +463,76 @@ def test_from_knee_points():
         assert sn_knee.name == sn.name
 
 
+def _kernel(fn):
+    """Return python-callable implementation for numba kernels."""
+    return getattr(fn, "py_func", fn)
+
+
+@pytest.mark.parametrize("sn", [DNV_B1C, DNV_B1A, EXOTIC])
+def test_cycles_kernel_boundary_regression(sn):
+    kernel = _kernel(_calc_cycles)
+
+    stress_values = [1e-12, sn.get_stress(1e4)[0], sn.get_stress(1e7)[0]]
+    stress_values.extend(list(sn.get_knee_stress()))
+    stress_values.extend([s * (1 - 1e-12) for s in sn.get_knee_stress()])
+    stress_values.extend([s * (1 + 1e-12) for s in sn.get_knee_stress()])
+    if sn.endurance < np.inf:
+        end_stress = sn.get_stress(sn.endurance)[0]
+        stress_values.extend(
+            [
+                end_stress,
+                end_stress * (1 - 1e-12),
+                end_stress * (1 + 1e-12),
+            ]
+        )
+    stress = np.asarray(stress_values, dtype=np.float64)
+
+    kernel_out = kernel(stress, sn.slope, sn.intercept, sn.endurance)
+    api_out = sn.get_cycles(stress)
+    np.testing.assert_allclose(kernel_out, api_out, rtol=1e-11, atol=0.0)
+
+    zero_stress = np.array([0.0], dtype=np.float64)
+    kernel_zero = kernel(zero_stress, sn.slope, sn.intercept, sn.endurance)
+    api_zero = sn.get_cycles(zero_stress)
+    np.testing.assert_allclose(kernel_zero, api_zero, rtol=0.0, atol=0.0)
+
+    with pytest.raises(AssertionError):
+        kernel(np.array([-1.0]), sn.slope, sn.intercept, sn.endurance)
+    with pytest.raises(AssertionError):
+        _ = sn.get_cycles(-1)
+
+
+@pytest.mark.parametrize("sn", [DNV_B1C, DNV_B1A, EXOTIC])
+def test_stress_kernel_boundary_regression(sn):
+    kernel = _kernel(_calc_stress)
+
+    cycle_values = [1e-12, 1e4, 1e7]
+    cycle_values.extend(list(sn.get_knee_cycles()))
+    cycle_values.extend([c * (1 - 1e-12) for c in sn.get_knee_cycles()])
+    cycle_values.extend([c * (1 + 1e-12) for c in sn.get_knee_cycles()])
+    if sn.endurance < np.inf:
+        cycle_values.extend(
+            [
+                sn.endurance,
+                sn.endurance * (1 - 1e-12),
+                sn.endurance * (1 + 1e-12),
+            ]
+        )
+
+    cycles = np.asarray(cycle_values, dtype=np.float64)
+    kernel_out = kernel(cycles, sn.slope, sn.intercept, sn.endurance)
+    api_out = sn.get_stress(cycles)
+    np.testing.assert_allclose(kernel_out, api_out, rtol=1e-11, atol=0.0)
+
+    with pytest.raises(AssertionError):
+        kernel(np.array([0.0]), sn.slope, sn.intercept, sn.endurance)
+    with pytest.raises(AssertionError):
+        _ = sn.get_stress(0)
+
+    with pytest.raises(AssertionError):
+        kernel(np.array([-1.0]), sn.slope, sn.intercept, sn.endurance)
+    with pytest.raises(AssertionError):
+        _ = sn.get_stress(-1)
+
+
 nb.config.DISABLE_JIT = False