From 301085fa1ba5d77051889c5c00a9a880541ab39a Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Fri, 7 Aug 2026 10:28:05 +0200
Subject: [PATCH 1/3] DX: support composable agent configurations

---
 .gitignore   | 4 ++--
 docs/conf.py | 2 ++
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/.gitignore b/.gitignore
index 5900756a..7541f658 100644
--- a/.gitignore
+++ b/.gitignore
@@ -42,8 +42,8 @@ pyvenv*/
 /.agents
 /.claude/
 /.codex
-/AGENTS.md
-/CLAUDE.md
+AGENTS.md
+CLAUDE.md
 
 # Exceptions
 !.cspell.json
diff --git a/docs/conf.py b/docs/conf.py
index af4b2f80..a3d3808e 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -107,6 +107,8 @@ def get_tensorflow_url() -> str:
     "**.ipynb_checkpoints",
     "*build",
     "adr*",
+    "AGENTS.md",
+    "CLAUDE.md",
     "tests",
 ]
 extensions = [

From e6fa1a2b4c52fc93d4b1edf7d6b06e3666fd6ed9 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Wed, 5 Aug 2026 21:54:41 +0200
Subject: [PATCH 2/3] BREAK: replace update_parameters() with pure evaluation
 API

---
 docs/amplitude-analysis.ipynb              | 35 +++++++++++-----------
 docs/usage.ipynb                           | 34 ++++++++++++++-------
 docs/usage/basics.ipynb                    | 12 ++++----
 docs/usage/binned-fit.ipynb                | 16 +++++++---
 docs/usage/chi-squared.ipynb               |  5 ++--
 docs/usage/unbinned-fit.ipynb              |  6 ++--
 src/tensorwaves/estimator.py               |  8 ++---
 src/tensorwaves/function/__init__.py       | 27 +++++++++++++----
 src/tensorwaves/function/sympy/__init__.py |  6 ++--
 src/tensorwaves/interface.py               | 33 +++++++++++++++-----
 tests/function/test_function.py            | 34 ++++++++++++++++-----
 tests/test_estimator.py                    |  3 ++
 12 files changed, 147 insertions(+), 72 deletions(-)

diff --git a/docs/amplitude-analysis.ipynb b/docs/amplitude-analysis.ipynb
index e5889ba9..5447b28f 100644
--- a/docs/amplitude-analysis.ipynb
+++ b/docs/amplitude-analysis.ipynb
@@ -1395,7 +1395,7 @@
     "\n",
     "Let's have a look at our [first guess for the parameter values](#determine-free-parameters). Recall that a {class}`.ParametrizedFunction` object computes the intensity for a certain {obj}`.DataSample`. This can be seen nicely when we use these intensities as weights on the phase space sample and plot it together with the original data sample. Here, we look at the invariant mass distribution projection of the final states `1` and `2`, which, [as we saw before](compwa-step-2.3), is the final state particle pair $\\pi^0\\pi^0$.\n",
     "\n",
-    "Don't forget to use {meth}`~.ParametrizedFunction.update_parameters` first!"
+    "Don't forget to first create a function with these initial parameter values using {meth}`~.ParametrizedFunction.with_parameters`!"
    ]
   },
   {
@@ -1473,8 +1473,8 @@
    "outputs": [],
    "source": [
     "original_parameters = optimized_function.parameters\n",
-    "optimized_function.update_parameters(initial_parameters)\n",
-    "compare_model(\"m_12\", data_real, phsp_real, optimized_function)"
+    "initial_function = optimized_function.with_parameters(initial_parameters)\n",
+    "compare_model(\"m_12\", data_real, phsp_real, initial_function)"
    ]
   },
   {
@@ -1622,7 +1622,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Using the same method as above, we renew the parameters of the {class}`.ParametrizedFunction` and plot it again over the phase space sample."
+    "Using the same method as above, we create a new {class}`.ParametrizedFunction` with the optimized parameter values and plot it again over the phase space sample."
    ]
   },
   {
@@ -1631,7 +1631,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "optimized_function.update_parameters(fit_result.parameter_values)\n",
+    "optimized_function = optimized_function.with_parameters(fit_result.parameter_values)\n",
     "compare_model(\"m_12\", data_real, phsp_real, optimized_function)"
    ]
   },
@@ -1752,22 +1752,21 @@
     "    input_data: DataSample,\n",
     "    resonances: list[str],\n",
     "):\n",
-    "    original_parameters = dict(func.parameters)\n",
     "    negative_lookahead = f\"(?!{'|'.join(map(re.escape, resonances))})\"\n",
     "    # https://regex101.com/r/WrgGyD/1\n",
     "    pattern = rf\"^(\\\\mathcal{{H}}|C_)({negative_lookahead}.)*$\"\n",
-    "    set_parameters_to_zero(func, pattern)\n",
-    "    array = func(input_data)\n",
-    "    func.update_parameters(original_parameters)\n",
-    "    return array\n",
-    "\n",
-    "\n",
-    "def set_parameters_to_zero(func: ParametrizedFunction, name_pattern: str) -> None:\n",
-    "    new_parameters = dict(func.parameters)\n",
-    "    for par_name in func.parameters:\n",
-    "        if re.match(name_pattern, par_name) is not None:\n",
-    "            new_parameters[par_name] = 0\n",
-    "    func.update_parameters(new_parameters)"
+    "    zeroed_parameters = get_zeroed_parameters(func, pattern)\n",
+    "    return func(input_data, zeroed_parameters)\n",
+    "\n",
+    "\n",
+    "def get_zeroed_parameters(\n",
+    "    func: ParametrizedFunction, name_pattern: str\n",
+    ") -> dict[str, complex]:\n",
+    "    return {\n",
+    "        par_name: 0\n",
+    "        for par_name in func.parameters\n",
+    "        if re.match(name_pattern, par_name) is not None\n",
+    "    }"
    ]
   },
   {
diff --git a/docs/usage.ipynb b/docs/usage.ipynb
index 4ebebbb3..60d6a757 100644
--- a/docs/usage.ipynb
+++ b/docs/usage.ipynb
@@ -144,9 +144,12 @@
     "bin_values, bin_edges, _ = ax.hist(data[\"x\"], bins=50, alpha=0.7, label=\"data\")\n",
     "x_values = (bin_edges[1:] + bin_edges[:-1]) / 2\n",
     "y_values = bin_values\n",
-    "function.update_parameters(initial_parameters)\n",
     "lines = ax.plot(\n",
-    "    x_values, function({\"x\": x_values}), c=\"red\", linewidth=2, label=\"model\"\n",
+    "    x_values,\n",
+    "    function({\"x\": x_values}, initial_parameters),\n",
+    "    c=\"red\",\n",
+    "    linewidth=2,\n",
+    "    label=\"model\",\n",
     ")\n",
     "ax.legend(loc=\"upper right\")\n",
     "plt.show()"
@@ -201,6 +204,7 @@
     "class FitAnimation(Callback):\n",
     "    def __init__(self, data, function, x_values, output_file, estimated_iterations=140):\n",
     "        self.__function = function\n",
+    "        self.__parameters = dict(function.parameters)\n",
     "        self.__fig, (self.__ax1, self.__ax2) = plt.subplots(\n",
     "            nrows=2, figsize=(7, 7), tight_layout=True\n",
     "        )\n",
@@ -208,7 +212,7 @@
     "        self.__ax1.hist(data[\"x\"], bins=50, alpha=0.7, label=\"data\")\n",
     "        self.__line = self.__ax1.plot(\n",
     "            x_values,\n",
-    "            function({\"x\": x_values}),\n",
+    "            function({\"x\": x_values}, self.__parameters),\n",
     "            c=\"red\",\n",
     "            linewidth=2,\n",
     "            label=\"model\",\n",
@@ -217,12 +221,12 @@
     "\n",
     "        self.__par_lines = [\n",
     "            self.__ax2.plot(0, value, label=par)[0]\n",
-    "            for par, value in function.parameters.items()\n",
+    "            for par, value in self.__parameters.items()\n",
     "        ]\n",
     "        self.__ax2.set_xlim(0, estimated_iterations)\n",
     "        self.__ax2.set_title(\"Parameter values\")\n",
     "        self.__ax2.legend(\n",
-    "            [f\"${sp.latex(sp.Symbol(par_name))}$\" for par_name in function.parameters],\n",
+    "            [f\"${sp.latex(sp.Symbol(par_name))}$\" for par_name in self.__parameters],\n",
     "            loc=\"upper right\",\n",
     "        )\n",
     "\n",
@@ -230,33 +234,41 @@
     "        self.__writer.setup(self.__fig, outfile=output_file)\n",
     "\n",
     "    def on_optimize_start(self, logs):\n",
+    "        self._update_parameters(logs)\n",
     "        self._update_plot()\n",
     "\n",
     "    def on_optimize_end(self, logs):\n",
+    "        self._update_parameters(logs)\n",
     "        self._update_plot()\n",
     "        self.__writer.finish()\n",
     "\n",
     "    def on_iteration_end(self, iteration, logs):\n",
+    "        self._update_parameters(logs)\n",
     "        self._update_plot()\n",
     "        self.__writer.finish()\n",
     "\n",
     "    def on_function_call_end(self, function_call, logs):\n",
+    "        self._update_parameters(logs)\n",
     "        self._update_plot()\n",
     "\n",
+    "    def _update_parameters(self, logs):\n",
+    "        if logs is not None:\n",
+    "            self.__parameters.update(logs[\"parameters\"])\n",
+    "\n",
     "    def _update_plot(self):\n",
     "        self._update_parametrization_plot()\n",
     "        self._update_traceback()\n",
     "        self.__writer.grab_frame()\n",
     "\n",
     "    def _update_parametrization_plot(self):\n",
-    "        title = self._render_parameters(self.__function.parameters)\n",
+    "        title = self._render_parameters(self.__parameters)\n",
     "        self.__ax1.set_title(title)\n",
-    "        self.__line.set_ydata(self.__function({\"x\": x_values}))\n",
+    "        self.__line.set_ydata(self.__function({\"x\": x_values}, self.__parameters))\n",
     "\n",
     "    def _update_traceback(self):\n",
     "        for line in self.__par_lines:\n",
     "            par_name = line.get_label()\n",
-    "            new_value = function.parameters[par_name]\n",
+    "            new_value = self.__parameters[par_name]\n",
     "            x = line.get_xdata()\n",
     "            x = [*x, x[-1] + 1]\n",
     "            y = [*line.get_ydata(), new_value]\n",
@@ -730,13 +742,13 @@
     ")\n",
     "def plot(dphi, k_r, k_phi, sigma):\n",
     "    global color_mesh, X, Y\n",
-    "    polar_function.update_parameters({\n",
+    "    parameters = {\n",
     "        R\"\\Delta\\phi\": dphi,\n",
     "        \"k_r\": k_r,\n",
     "        \"k_phi\": k_phi,\n",
     "        \"sigma\": sigma,\n",
-    "    })\n",
-    "    Z = polar_function(polar_domain)\n",
+    "    }\n",
+    "    Z = polar_function(polar_domain, parameters)\n",
     "    if color_mesh is not None:\n",
     "        color_mesh.remove()\n",
     "    color_mesh = ax_interactive.pcolormesh(X, Y, Z, cmap=\"coolwarm\")"
diff --git a/docs/usage/basics.ipynb b/docs/usage/basics.ipynb
index ae16c143..51635e3a 100644
--- a/docs/usage/basics.ipynb
+++ b/docs/usage/basics.ipynb
@@ -589,7 +589,7 @@
    "source": [
     "For the rest, the procedure is really just the same as that sketched in {ref}`compwa-step-3`.\n",
     "\n",
-    "We tweak the parameters a bit, then use {meth}`.ParametrizedBackendFunction.update_parameters` to change the function..."
+    "We tweak the parameters a bit, then use {meth}`.ParametrizedBackendFunction.with_parameters` to create a new function with these parameter values..."
    ]
   },
   {
@@ -606,7 +606,7 @@
     "    \"sigma_0\": 0.4,\n",
     "    \"sigma_1\": 0.4,\n",
     "}\n",
-    "function_1d.update_parameters(initial_parameters)"
+    "function_1d = function_1d.with_parameters(initial_parameters)"
    ]
   },
   {
@@ -730,7 +730,7 @@
    "outputs": [],
    "source": [
     "optimized_parameters = fit_result.parameter_values\n",
-    "function_1d.update_parameters(optimized_parameters)"
+    "function_1d = function_1d.with_parameters(optimized_parameters)"
    ]
   },
   {
@@ -1008,7 +1008,7 @@
     "    \"sigma_0\": 0.4,\n",
     "    \"sigma_1\": 0.4,\n",
     "}\n",
-    "function_2d.update_parameters(initial_parameters)"
+    "function_2d = function_2d.with_parameters(initial_parameters)"
    ]
   },
   {
@@ -1148,7 +1148,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If we update the parameters in the {class}`.ParametrizedFunction` with the optimized parameter values found by the {class}`.Optimizer`, we can compare the data distribution with the function."
+    "If we create a new {class}`.ParametrizedFunction` with the optimized parameter values found by the {class}`.Optimizer`, we can compare the data distribution with the function."
    ]
   },
   {
@@ -1158,7 +1158,7 @@
    "outputs": [],
    "source": [
     "optimized_parameters = fit_result.parameter_values\n",
-    "function_2d.update_parameters(optimized_parameters)"
+    "function_2d = function_2d.with_parameters(optimized_parameters)"
    ]
   },
   {
diff --git a/docs/usage/binned-fit.ipynb b/docs/usage/binned-fit.ipynb
index 7d858ff0..6dbf719e 100644
--- a/docs/usage/binned-fit.ipynb
+++ b/docs/usage/binned-fit.ipynb
@@ -126,11 +126,15 @@
    },
    "outputs": [],
    "source": [
-    "function.update_parameters(initial_parameters)\n",
     "fig, ax = plt.subplots(figsize=(8, 5))\n",
     "ax.set_xlabel(\"$x$\")\n",
     "ax.hist(x_distribution, bins=n_bins, label=\"Data distribution\")\n",
-    "ax.plot(x_values, function({\"x\": x_values}), label=\"Initial fit model\", c=\"red\")\n",
+    "ax.plot(\n",
+    "    x_values,\n",
+    "    function({\"x\": x_values}, initial_parameters),\n",
+    "    label=\"Initial fit model\",\n",
+    "    c=\"red\",\n",
+    ")\n",
     "ax.legend()\n",
     "plt.show()"
    ]
@@ -198,11 +202,15 @@
    },
    "outputs": [],
    "source": [
-    "function.update_parameters(fit_result.parameter_values)\n",
     "fig, ax = plt.subplots(figsize=(8, 5))\n",
     "ax.set_xlabel(\"$x$\")\n",
     "ax.hist(x_distribution, bins=n_bins, label=\"Data distribution\")\n",
-    "ax.plot(x_values, function({\"x\": x_values}), label=\"Optimized model\", c=\"red\")\n",
+    "ax.plot(\n",
+    "    x_values,\n",
+    "    function({\"x\": x_values}, fit_result.parameter_values),\n",
+    "    label=\"Optimized model\",\n",
+    "    c=\"red\",\n",
+    ")\n",
     "ax.legend()\n",
     "plt.show()"
    ]
diff --git a/docs/usage/chi-squared.ipynb b/docs/usage/chi-squared.ipynb
index 72fa45f0..65f8af9c 100644
--- a/docs/usage/chi-squared.ipynb
+++ b/docs/usage/chi-squared.ipynb
@@ -98,7 +98,7 @@
    "source": [
     "original_parameters = function.parameters\n",
     "initial_parameters = {\"a\": -25, \"b\": 1.5, \"c\": 2.6}\n",
-    "function.update_parameters(initial_parameters)"
+    "function = function.with_parameters(initial_parameters)"
    ]
   },
   {
@@ -207,7 +207,8 @@
    },
    "outputs": [],
    "source": [
-    "compare_model(function, x_values, observed_y)"
+    "optimized_function = function.with_parameters(fit_result.parameter_values)\n",
+    "compare_model(optimized_function, x_values, observed_y)"
    ]
   }
  ],
diff --git a/docs/usage/unbinned-fit.ipynb b/docs/usage/unbinned-fit.ipynb
index a59b8126..8273d405 100644
--- a/docs/usage/unbinned-fit.ipynb
+++ b/docs/usage/unbinned-fit.ipynb
@@ -182,8 +182,7 @@
     "Y = np.linspace(*ylim, bins_y)\n",
     "X, Y = np.meshgrid(X, Y)\n",
     "\n",
-    "function.update_parameters(initial_parameters)\n",
-    "Z = function({\"x\": X, \"y\": Y})\n",
+    "Z = function({\"x\": X, \"y\": Y}, initial_parameters)\n",
     "\n",
     "fig, (ax1, ax2) = plt.subplots(figsize=(8, 7), nrows=2, sharex=True, tight_layout=True)\n",
     "ax1.set_title(\"Data distribution\")\n",
@@ -258,8 +257,7 @@
    },
    "outputs": [],
    "source": [
-    "function.update_parameters(fit_result.parameter_values)\n",
-    "Z = function({\"x\": X, \"y\": Y})\n",
+    "Z = function({\"x\": X, \"y\": Y}, fit_result.parameter_values)\n",
     "\n",
     "fig, (ax1, ax2) = plt.subplots(figsize=(8, 7), nrows=2, sharex=True, tight_layout=True)\n",
     "ax1.set_title(\"Data distribution\")\n",
diff --git a/src/tensorwaves/estimator.py b/src/tensorwaves/estimator.py
index 3e3ef68e..dbc3f28c 100644
--- a/src/tensorwaves/estimator.py
+++ b/src/tensorwaves/estimator.py
@@ -146,8 +146,7 @@ def __init__(
         self.__sum = find_function("sum", backend)
 
     def __call__(self, parameters: Mapping[str, ParameterValue]) -> float:
-        self.__function.update_parameters(parameters)
-        computed_values = self.__function(self.__domain)
+        computed_values = self.__function(self.__domain, parameters)
         chi_squared = self.__weights * (computed_values - self.__observed_values) ** 2
         return self.__sum(chi_squared)
 
@@ -213,9 +212,8 @@ def __init__(
         self.__phsp_volume = phsp_volume
 
     def __call__(self, parameters: Mapping[str, ParameterValue]) -> float:
-        self.__function.update_parameters(parameters)
-        data_intensities = self.__function(self.__data)
-        phsp_intensities = self.__function(self.__phsp)
+        data_intensities = self.__function(self.__data, parameters)
+        phsp_intensities = self.__function(self.__phsp, parameters)
         if self.__phsp_weights is not None:
             phsp_intensities *= self.__phsp_weights
         normalization_integral = self.__phsp_volume * self.__mean(phsp_intensities)
diff --git a/src/tensorwaves/function/__init__.py b/src/tensorwaves/function/__init__.py
index 4b008013..325e2e4b 100644
--- a/src/tensorwaves/function/__init__.py
+++ b/src/tensorwaves/function/__init__.py
@@ -106,8 +106,12 @@ def __init__(
         self.__function = PositionalArgumentFunction(function, argument_order)
         self.__parameters = dict(parameters)
 
-    def __call__(self, data: DataSample) -> np.ndarray:
-        extended_data = {**data, **self.__parameters}
+    def __call__(
+        self,
+        data: DataSample,
+        parameters: Mapping[str, ParameterValue] | None = None,
+    ) -> np.ndarray:
+        extended_data = {**data, **self.__merge_parameters(parameters)}
         return self.__function(extended_data)  # ty:ignore[invalid-argument-type]
 
     @property
@@ -122,8 +126,21 @@ def argument_order(self) -> tuple[str, ...]:
     def parameters(self) -> dict[str, ParameterValue]:
         return dict(self.__parameters)
 
-    def update_parameters(self, new_parameters: Mapping[str, ParameterValue]) -> None:
-        over_defined = set(new_parameters) - set(self.__parameters)
+    def with_parameters(
+        self, parameters: Mapping[str, ParameterValue]
+    ) -> ParametrizedBackendFunction:
+        return ParametrizedBackendFunction(
+            function=self.function,
+            argument_order=self.argument_order,
+            parameters=self.__merge_parameters(parameters),
+        )
+
+    def __merge_parameters(
+        self, parameters: Mapping[str, ParameterValue] | None
+    ) -> dict[str, ParameterValue]:
+        if parameters is None:
+            return self.__parameters
+        over_defined = set(parameters) - set(self.__parameters)
         if over_defined:
             sep = "\n    "
             parameter_listing = f"{sep}".join(sorted(self.__parameters))
@@ -132,7 +149,7 @@ def update_parameters(self, new_parameters: Mapping[str, ParameterValue]) -> Non
                 f" Expecting one of:{sep}{parameter_listing}"
             )
             raise ValueError(msg)
-        self.__parameters.update(new_parameters)
+        return {**self.__parameters, **parameters}
 
 
 def get_source_code(function: Function) -> str:
diff --git a/src/tensorwaves/function/sympy/__init__.py b/src/tensorwaves/function/sympy/__init__.py
index 28bcad82..d11243e7 100644
--- a/src/tensorwaves/function/sympy/__init__.py
+++ b/src/tensorwaves/function/sympy/__init__.py
@@ -119,8 +119,10 @@ def create_parametrized_function(  # ruff:ignore[too-many-arguments]
         ... )
         >>> array = np.linspace(0, 1, num=5)
         >>> data = {"x": array, "y": array}
-        >>> function.update_parameters({"b": 1})
-        >>> function(data).tolist()
+        >>> function(data, {"b": 1}).tolist()
+        [0.0, 0.0, 0.0, 0.0, 0.0]
+        >>> function_b1 = function.with_parameters({"b": 1})
+        >>> function_b1(data).tolist()
         [0.0, 0.0, 0.0, 0.0, 0.0]
     """
     expression = _substitute_matrix_elements(expression)
diff --git a/src/tensorwaves/interface.py b/src/tensorwaves/interface.py
index 6a97b784..85c5e740 100644
--- a/src/tensorwaves/interface.py
+++ b/src/tensorwaves/interface.py
@@ -43,26 +43,43 @@ def __call__(self, data: InputType) -> OutputType: ...
 
 
 class ParametrizedFunction(Function[InputType, OutputType]):
-    """Interface of a callable function.
+    """Interface of a callable function with parameters.
 
     A `ParametrizedFunction` identifies certain variables in a mathematical expression
     as **parameters**. Remaining variables are considered **domain variables**. Domain
-    variables are the argument of the evaluation (see
-    :func:`~ParametrizedFunction.__call__`), while the parameters are controlled via
-    :attr:`parameters` (getter) and :meth:`update_parameters` (setter). This mechanism
-    is especially important for an `Estimator`.
+    variables are the first argument of the evaluation (see
+    :func:`~ParametrizedFunction.__call__`), while parameter values can be passed as
+    the second argument. Parameter values that are not provided at the call fall back
+    to the default values in :attr:`parameters`. A `ParametrizedFunction` is
+    immutable: a call never affects later calls, which makes it thread-safe and safe
+    to trace for JIT compilers like :code:`jax.jit`. Use :meth:`with_parameters` to
+    create a new function with different default parameter values.
 
     .. automethod:: __call__
     """
 
+    @abstractmethod
+    def __call__(
+        self,
+        data: InputType,
+        parameters: Mapping[str, ParameterValue] | None = None,
+    ) -> OutputType:
+        """Evaluate the function over :code:`data` for these parameter values.
+
+        Given parameter values are merged with the defaults in :attr:`parameters` for
+        this evaluation only.
+        """
+
     @property
     @abstractmethod
     def parameters(self) -> dict[str, ParameterValue]:
-        """`dict` of parameters."""
+        """`dict` of default parameter values."""
 
     @abstractmethod
-    def update_parameters(self, new_parameters: Mapping[str, ParameterValue]) -> None:
-        """Update the collection of parameters."""
+    def with_parameters(
+        self, parameters: Mapping[str, ParameterValue]
+    ) -> ParametrizedFunction[InputType, OutputType]:
+        """Create a new function with updated default parameter values."""
 
 
 class DataTransformer(Function[DataSample, DataSample]):
diff --git a/tests/function/test_function.py b/tests/function/test_function.py
index d25e66f1..25015636 100644
--- a/tests/function/test_function.py
+++ b/tests/function/test_function.py
@@ -58,23 +58,43 @@ def test_call(
     def test_function(self, function: ParametrizedBackendFunction):
         assert callable(function.function)
 
-    def test_update_parameter(self):
-        initial_parameter_values = {"a": 1, "b": 1}
+    def test_call_with_parameters(self):
+        initial_parameter_values = {"a": 1.0, "b": 2.0}
         func = ParametrizedBackendFunction(
             lambda a, b, x: a * x + b,
             argument_order=("a", "b", "x"),
             parameters=initial_parameter_values,
         )
+        data: DataSample = {"x": np.array([0.0, 1.0, 2.0])}
+        np.testing.assert_array_equal(func(data), [2.0, 3.0, 4.0])
+        np.testing.assert_array_equal(func(data, {"a": -1.0}), [2.0, 1.0, 0.0])
         with pytest.raises(
             ValueError,
             match=r"^Parameters {'c'} do not exist in function arguments\.",
         ):
-            func.update_parameters({"a": 2, "c": 1})
+            func(data, {"a": 2.0, "c": 1.0})
         assert func.parameters == initial_parameter_values
-        new_parameter_values = {"a": 2, "b": 2}
-        func.update_parameters(new_parameter_values)
-        assert func.parameters == new_parameter_values
-        assert new_parameter_values != initial_parameter_values
+        np.testing.assert_array_equal(func(data), [2.0, 3.0, 4.0])
+
+    def test_with_parameters(self):
+        initial_parameter_values = {"a": 1.0, "b": 2.0}
+        func = ParametrizedBackendFunction(
+            lambda a, b, x: a * x + b,
+            argument_order=("a", "b", "x"),
+            parameters=initial_parameter_values,
+        )
+        new_func = func.with_parameters({"a": 2.0})
+        assert new_func is not func
+        assert new_func.parameters == {"a": 2.0, "b": 2.0}
+        assert new_func.function is func.function
+        assert func.parameters == initial_parameter_values
+        data: DataSample = {"x": np.array([0.0, 1.0, 2.0])}
+        np.testing.assert_array_equal(new_func(data), [2.0, 4.0, 6.0])
+        with pytest.raises(
+            ValueError,
+            match=r"^Parameters {'c'} do not exist in function arguments\.",
+        ):
+            func.with_parameters({"c": 1.0})
 
 
 class TestPositionalArgumentFunction:
diff --git a/tests/test_estimator.py b/tests/test_estimator.py
index 6f41f77f..ef21a52e 100644
--- a/tests/test_estimator.py
+++ b/tests/test_estimator.py
@@ -32,6 +32,7 @@ def test_call(self, backend):
         assert estimator({}) == 0
         assert estimator({"b": 2}) == 5.0
         assert estimator({"a": 1, "b": 2}) == 14.0
+        assert function.parameters == {"a": 0, "b": 1}, "estimator call is not pure"
         estimator = ChiSquared(
             function,
             x_data,
@@ -213,6 +214,7 @@ def test_sympy_unbinned_nll(
     true_params: dict[str, ParameterValue],
     phsp: DataSample,
 ):
+    original_parameters = function.parameters
     estimator = UnbinnedNLL(
         function,
         data,
@@ -224,6 +226,7 @@ def test_sympy_unbinned_nll(
         estimator,
         initial_parameters=true_params,
     )
+    assert function.parameters == original_parameters, "optimize() is not pure"
 
     par_values = fit_result.parameter_values
     par_errors = fit_result.parameter_errors

From b818ba1848acacb8c0694763bf7de6ab259e6987 Mon Sep 17 00:00:00 2001
From: GitHub <noreply@github.com>
Date: Fri, 7 Aug 2026 13:48:39 +0000
Subject: [PATCH 3/3] MAINT: implement updates from formatters

---
 docs/amplitude-analysis.ipynb | 2 +-
 docs/usage.ipynb              | 2 +-
 docs/usage/basics.ipynb       | 2 +-
 docs/usage/binned-fit.ipynb   | 2 +-
 docs/usage/chi-squared.ipynb  | 2 +-
 docs/usage/unbinned-fit.ipynb | 2 +-
 6 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/docs/amplitude-analysis.ipynb b/docs/amplitude-analysis.ipynb
index 5447b28f..0c59c994 100644
--- a/docs/amplitude-analysis.ipynb
+++ b/docs/amplitude-analysis.ipynb
@@ -1907,7 +1907,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.13"
+   "version": "3.13.14"
   }
  },
  "nbformat": 4,
diff --git a/docs/usage.ipynb b/docs/usage.ipynb
index 60d6a757..6a462c0a 100644
--- a/docs/usage.ipynb
+++ b/docs/usage.ipynb
@@ -805,7 +805,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.12"
+   "version": "3.13.14"
   }
  },
  "nbformat": 4,
diff --git a/docs/usage/basics.ipynb b/docs/usage/basics.ipynb
index 51635e3a..defae220 100644
--- a/docs/usage/basics.ipynb
+++ b/docs/usage/basics.ipynb
@@ -1263,7 +1263,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.13"
+   "version": "3.13.14"
   }
  },
  "nbformat": 4,
diff --git a/docs/usage/binned-fit.ipynb b/docs/usage/binned-fit.ipynb
index 6dbf719e..a003e130 100644
--- a/docs/usage/binned-fit.ipynb
+++ b/docs/usage/binned-fit.ipynb
@@ -278,7 +278,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.12"
+   "version": "3.13.14"
   }
  },
  "nbformat": 4,
diff --git a/docs/usage/chi-squared.ipynb b/docs/usage/chi-squared.ipynb
index 65f8af9c..9b66b1ab 100644
--- a/docs/usage/chi-squared.ipynb
+++ b/docs/usage/chi-squared.ipynb
@@ -231,7 +231,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.12"
+   "version": "3.13.14"
   }
  },
  "nbformat": 4,
diff --git a/docs/usage/unbinned-fit.ipynb b/docs/usage/unbinned-fit.ipynb
index 8273d405..ae0de6c4 100644
--- a/docs/usage/unbinned-fit.ipynb
+++ b/docs/usage/unbinned-fit.ipynb
@@ -290,7 +290,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.12"
+   "version": "3.13.14"
   }
  },
  "nbformat": 4,