Repository: ComPWA/tensorwaves (language: python)

Issue:
Re-implement performance optimization
#357 removes `Model.performance_optimize()`, because this feature was not well tested and has to be reconsidered with the upcoming interface changes. The relevant old code can be found here:

- [`SympyModel.performance_optimize()`](https://github.com/ComPWA/tensorwaves/blob/18fce7df9153bb241bf0d87b0b2df24cf0d87de0/src/tensorwaves/model/sympy.py#L338-L339)
- [`_ConstantSubExpressionSympyModel`](https://github.com/ComPWA/tensorwaves/blob/18fce7df9153bb241bf0d87b0b2df24cf0d87de0/src/tensorwaves/model/sympy.py#L229-L357)
- The locations [where these were used](https://github.com/ComPWA/tensorwaves/blob/908a33f2f46f59335d8c571fda838a26cd4f1641/src/tensorwaves/estimator.py#L80-L87) in `Estimator`

Closes #330 


Discussion / hints:

Fix patch:
diff --git a/codecov.yml b/codecov.yml
index 184eda1a..5bef0c52 100644
--- a/codecov.yml
+++ b/codecov.yml
@@ -1,5 +1,5 @@
 codecov:
-  require_ci_to_pass: yes
+  require_ci_to_pass: true
 
 coverage:
   precision: 2
@@ -10,7 +10,7 @@ coverage:
     project:
       default:
         # basic
-        target: 70% # can't go below this percentage
+        target: 80% # can't go below this percentage
         threshold: 3% # allow drops by this percentage
         base: auto
         # advanced
@@ -35,7 +35,7 @@ coverage:
         only_pulls: false
         flags: null
         paths: null
-    changes: no
+    changes: false
 
 parsers:
   gcov:
diff --git a/docs/index.md b/docs/index.md
index 814eb725..105fe848 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -28,13 +28,15 @@ The decision to use Python only has been further inspired by
 
 :::
 
-TensorWaves is a fitter package that optimizes mathematical models to data
-samples. The models can be any mathematical expression (best expressed with
-{mod}`sympy`) that is then converted to any computational backend. In addition,
-the {mod}`tensorwaves.data` module allows one to generate toy data samples. The
-package is developed in parallel with {doc}`AmpForm <ampform:index>`, which
-implements physics models, but its 'lambdifying' mechanisms can in principle be
-used independently.
+TensorWaves is a fitter package that can fit models to data samples. The models
+can be any mathematical expression (best expressed with {mod}`sympy`) that is
+then expressed to any computational back-end. In addition, the
+{mod}`tensorwaves.data` module allows one to generate data distributions for
+those mathematical expressions.
+
+The package is developed in parallel with {doc}`AmpForm <ampform:index>`, which
+implements physics models, but its mechanisms for creating functions can in
+principle be used independently.
 
 ```{link-button} usage
 ---
diff --git a/docs/usage.ipynb b/docs/usage.ipynb
index 75e277a3..b737bac6 100644
--- a/docs/usage.ipynb
+++ b/docs/usage.ipynb
@@ -55,7 +55,7 @@
    "source": [
     "TensorWaves is a package for fitting general mathematical expressions to data distributions. The fundamentals behind the package are illustrated in {doc}`/usage/basics`.\n",
     "\n",
-    "While general in design, the package is intended for doing {doc}`Partial Wave Analysis <pwa:index>`. First, the {mod}`ampform` package determines which transitions are allowed from some initial state to a final state. It then formulates those transitions mathematically as an {ref}`amplitude model <usage:Construct a model>`. TensorWaves can then {meth}`~.Model.lambdify` this expression to some computational backend. Finally, TensorWaves {ref}`'fits' <usage:Optimize the model>` this model to some data sample. Optionally, a data sample can be {ref}`generated <usage:Generate data sample>` from the model.\n",
+    "While general in design, the package is intended for doing {doc}`Partial Wave Analysis <pwa:index>`. First, the {mod}`qrules` package determines which transitions are allowed from some initial state to a final state. The {mod}`ampform` package then formulates those transitions mathematically as an {ref}`amplitude model <usage:Construct a model>`. TensorWaves can then create an function ({func}`.create_function`) for this expression in some computational back-end. Finally, TensorWaves can {ref}`'fit' <usage:Optimize the model>` this function to a data sample. Optionally, a data sample can be {ref}`generated <usage:Generate data sample>` from the model.\n",
     "\n",
     "This page shows a brief overview of the complete workflow. More info about each step can be found under {ref}`usage:Step-by-step workflow`."
    ]
@@ -94,7 +94,7 @@
     "from tensorwaves.data import generate_data, generate_phsp\n",
     "from tensorwaves.data.transform import HelicityTransformer\n",
     "from tensorwaves.estimator import UnbinnedNLL\n",
-    "from tensorwaves.model import LambdifiedFunction, SympyModel\n",
+    "from tensorwaves.function.sympy import create_function\n",
     "from tensorwaves.optimizer.callbacks import CSVSummary\n",
     "from tensorwaves.optimizer.minuit import Minuit2\n",
     "\n",
@@ -205,11 +205,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sympy_model = SympyModel(\n",
+    "intensity = create_function(\n",
     "    expression=model.expression.doit(),\n",
     "    parameters=model.parameter_defaults,\n",
+    "    backend=\"jax\",\n",
     ")\n",
-    "intensity = LambdifiedFunction(sympy_model, backend=\"jax\")\n",
     "data_converter = HelicityTransformer(model.adapter)\n",
     "reaction_info = model.adapter.reaction_info\n",
     "initial_state_mass = reaction_info.initial_state[-1].mass\n",
@@ -224,6 +224,19 @@
     ")"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "phsp_set = data_converter.transform(phsp_sample)\n",
+    "data_set = data_converter.transform(data_sample)\n",
+    "data_frame = pd.DataFrame(data_set)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -232,7 +245,8 @@
      "source_hidden": true
     },
     "tags": [
-     "hide-cell"
+     "full-width",
+     "hide-input"
     ]
    },
    "outputs": [],
@@ -260,36 +274,9 @@
     "    for i, p in enumerate(intermediate_states):\n",
     "        plt.gca().axvline(\n",
     "            x=p.mass, linestyle=\"dotted\", label=p.name, color=colors[i]\n",
-    "        )"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "tags": []
-   },
-   "outputs": [],
-   "source": [
-    "phsp_set = data_converter.transform(phsp_sample)\n",
-    "data_set = data_converter.transform(data_sample)\n",
-    "data_frame = pd.DataFrame(data_set)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "jupyter": {
-     "source_hidden": true
-    },
-    "tags": [
-     "full-width",
-     "hide-input"
-    ]
-   },
-   "outputs": [],
-   "source": [
+    "        )\n",
+    "\n",
+    "\n",
     "fig, ax = plt.subplots(figsize=(9, 4))\n",
     "data_frame[\"m_12\"].hist(bins=100, alpha=0.5, density=True, ax=ax)\n",
     "indicate_masses()\n",
@@ -341,15 +328,15 @@
     "    data = data_set[variable_name]\n",
     "    phsp = phsp_set[variable_name]\n",
     "    intensities = intensity_model(phsp_set)\n",
-    "    fig, ax = plt.subplots(figsize=(9, 4))\n",
-    "    plt.hist(\n",
+    "    _, ax = plt.subplots(figsize=(9, 4))\n",
+    "    ax.hist(\n",
     "        data,\n",
     "        bins=bins,\n",
     "        alpha=0.5,\n",
     "        label=\"data\",\n",
     "        density=True,\n",
     "    )\n",
-    "    plt.hist(\n",
+    "    ax.hist(\n",
     "        phsp,\n",
     "        weights=intensities,\n",
     "        bins=bins,\n",
@@ -359,7 +346,7 @@
     "        density=True,\n",
     "    )\n",
     "    indicate_masses()\n",
-    "    plt.legend()"
+    "    ax.legend()"
    ]
   },
   {
@@ -381,6 +368,7 @@
     "    \"Gamma_f(0)(980)\": 0.1,\n",
     "    \"Gamma_f(0)(1710)\": 0.2,\n",
     "}\n",
+    "original_parameters = intensity.parameters\n",
     "intensity.update_parameters(initial_parameters)"
    ]
   },
@@ -480,7 +468,7 @@
     "    label = line.get_label()\n",
     "    color = line.get_color()\n",
     "    ax2.axhline(\n",
-    "        y=complex(sympy_model.parameters[label]).real,\n",
+    "        y=complex(original_parameters[label]).real,\n",
     "        color=color,\n",
     "        alpha=0.5,\n",
     "        linestyle=\"dotted\",\n",
diff --git a/docs/usage/analytic-continuation.ipynb b/docs/usage/analytic-continuation.ipynb
index 477ec35b..5322e1bc 100644
--- a/docs/usage/analytic-continuation.ipynb
+++ b/docs/usage/analytic-continuation.ipynb
@@ -72,7 +72,7 @@
     "\n",
     "from tensorwaves.data import generate_data, generate_phsp\n",
     "from tensorwaves.data.transform import HelicityTransformer\n",
-    "from tensorwaves.model import LambdifiedFunction, SympyModel\n",
+    "from tensorwaves.function.sympy import create_function\n",
     "\n",
     "logger = logging.getLogger()\n",
     "logger.setLevel(logging.ERROR)\n",
@@ -193,11 +193,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sympy_model = SympyModel(\n",
+    "intensity = create_function(\n",
     "    expression=model.expression.doit(),\n",
     "    parameters=model.parameter_defaults,\n",
+    "    backend=\"jax\",\n",
     ")\n",
-    "intensity = LambdifiedFunction(sympy_model, backend=\"jax\")\n",
     "data_converter = HelicityTransformer(model.adapter)\n",
     "reaction_info = model.adapter.reaction_info\n",
     "initial_state_mass = reaction_info.initial_state[-1].mass\n",
diff --git a/docs/usage/basics.ipynb b/docs/usage/basics.ipynb
index 3e2cfebc..0edb4287 100644
--- a/docs/usage/basics.ipynb
+++ b/docs/usage/basics.ipynb
@@ -51,14 +51,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The {doc}`/usage` page illustrates how to use {mod}`tensorwaves` in combination with {mod}`ampform`. At core, however, {mod}`tensorwaves` is a package that can 'fit' arbitrary models to a data set using different backends, as well as generate toy data samples.\n",
+    "The {doc}`/usage` page illustrates how to use {mod}`tensorwaves` in combination with {mod}`ampform`. At core, however, {mod}`tensorwaves` is a package that can 'fit' arbitrary mathematical expressions to a data set using different computational back-ends. It can also use those expressions to describe a distribution over which to generate data samples.\n",
     "\n",
-    "These pages illustrate what's going on behind the scenes with some simple 1-dimensional and 2-dimensional models. Since we don't have a data sample to fit a model to, we follow the the same procedure as in the {doc}`/usage` page:\n",
+    "This page illustrate what's going on behind the scenes with some simple 1-dimensional and 2-dimensional expressions. The workflow is just like in {doc}`/usage`, but without using {mod}`ampform`:\n",
     "\n",
-    "1. Formulate a mathematical description with {mod}`sympy`.\n",
-    "2. Generate a 'toy' data sample with hit-and-miss.\n",
-    "3. Convert the model to some backend with {mod}`tensorwaves`.\n",
-    "4. Tweak the parameters and fit the modified to the generated data sample."
+    "1. Formulate a mathematical expression with {mod}`sympy`.\n",
+    "2. Generate a distribution data sample for that expression.\n",
+    "3. Express the expression as a function in some computational back-end.\n",
+    "4. Tweak the {attr}`~.Function.parameters` and fit the {class}`.Function` to the generated distribution."
    ]
   },
   {
@@ -74,7 +74,7 @@
    },
    "outputs": [],
    "source": [
-    "!pip install black"
+    "!pip -qqq install black"
    ]
   },
   {
@@ -99,8 +99,8 @@
     "from sympy.plotting import plot3d\n",
     "\n",
     "from tensorwaves.estimator import UnbinnedNLL\n",
+    "from tensorwaves.function.sympy import create_function\n",
     "from tensorwaves.interface import Function\n",
-    "from tensorwaves.model import LambdifiedFunction, SympyModel\n",
     "from tensorwaves.optimizer.callbacks import CSVSummary\n",
     "from tensorwaves.optimizer.minuit import Minuit2"
    ]
@@ -118,7 +118,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "First of all, we'll formulate a 1-dimensional model with {mod}`sympy`. In this example, we take a sum of [Gaussians](https://en.wikipedia.org/wiki/Gaussian_function) plus some [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution) distribution."
+    "First of all, we'll formulate some expression with {mod}`sympy`. In this example, we take a sum of [Gaussians](https://en.wikipedia.org/wiki/Gaussian_function) plus some [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution) distribution."
    ]
   },
   {
@@ -231,9 +231,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "So far, all we did was using {mod}`sympy` to mathematically formulate a model. We now need to {func}`~sympy.utilities.lambdify.lambdify` that expression to some computational backend, so that we can efficiently generate data and/or optimize the parameters to 'fit' the model to some data sample.\n",
-    "\n",
-    "{mod}`tensorwaves` does this in two intermediate stages. First, the expression and parameter defaults are expressed in terms of a {class}`.SympyModel` (interface: {class}`.Model`) that serves as a template for different computational backends."
+    "So far, all we did was using {mod}`sympy` to symbolically formulate a mathematical expression. We now need to {func}`~sympy.utilities.lambdify.lambdify` that expression to some computational backend, so that we can efficiently generate data and/or optimize the parameters in the function to 'fit' the model to some data sample. TensorWaves can do this with the {func}`.create_function` function:"
    ]
   },
   {
@@ -242,9 +240,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model_1d = SympyModel(\n",
+    "function_1d = create_function(\n",
     "    expression=expression_1d,\n",
     "    parameters=parameter_defaults,\n",
+    "    backend=\"numpy\",\n",
     ")"
    ]
   },
@@ -252,23 +251,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Next, we use this {class}`.SympyModel` as a template to construct a {class}`.LambdifiedFunction` (interface: {class}`.Function`). We choose {mod}`numpy` as computational backend, as we're not doing any optimizing yet (for optimizing it's better to use {obj}`jax <jax.jit>`)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "function_1d = LambdifiedFunction(model_1d, backend=\"numpy\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Internally, the {class}`.LambdifiedFunction` carries some source code that numpy understands. And it indeed looks similar to the expression that we formulated in {ref}`usage/basics:Formulate model`:"
+    "The resulting {class}`.LambdifiedFunction` internally carries some source code that {mod}`numpy` understands. And it indeed looks similar to the expression that we formulated in {ref}`usage/basics:Formulate model`:"
    ]
   },
   {
@@ -288,12 +271,11 @@
     "\n",
     "from black import FileMode, format_str\n",
     "\n",
-    "print(\n",
-    "    format_str(\n",
-    "        inspect.getsource(function_1d._LambdifiedFunction__lambdified_model),\n",
-    "        mode=FileMode(),\n",
-    "    )\n",
-    ")"
+    "src = format_str(\n",
+    "    inspect.getsource(function_1d._LambdifiedFunction__function),\n",
+    "    mode=FileMode(),\n",
+    ")\n",
+    "print(src)"
    ]
   },
   {
@@ -373,7 +355,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "So, we now have a function $f$ of $x$ expressed in some computational backend. This function is to describe a distribution of $x$ values. In the real world, these $x$-values are some observables from a process you measure. But sometimes, it's useful to generate a 'toy' data sample for your model function as well, to try it out."
+    "So, we now have a function $f$ of $x$ expressed in some computational backend. This function is to describe a distribution over $x$. In the real world, $x$ is an observable from a process you measure. But sometimes, it's useful to generate a 'toy' data sample for your model function as well, to try it out."
    ]
   },
   {
@@ -737,7 +719,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "estimator = UnbinnedNLL(model_1d, data, domain_sample, backend=\"jax\")"
+    "estimator = UnbinnedNLL(function_1d, data, domain_sample, backend=\"jax\")"
    ]
   },
   {
@@ -891,11 +873,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model_2d = SympyModel(\n",
+    "function_2d = create_function(\n",
     "    expression=expression_2d,\n",
     "    parameters=parameter_defaults,\n",
-    ")\n",
-    "function_2d = LambdifiedFunction(model_2d, backend=\"numpy\")"
+    "    backend=\"numpy\",\n",
+    ")"
    ]
   },
   {
@@ -1008,7 +990,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "estimator = UnbinnedNLL(model_2d, data, domain_sample, backend=\"jax\")\n",
+    "estimator = UnbinnedNLL(function_2d, data, domain_sample, backend=\"jax\")\n",
     "minuit2 = Minuit2(callback=CSVSummary(\"traceback.csv\"))\n",
     "fit_result = minuit2.optimize(estimator, initial_parameters)\n",
     "fit_result"
diff --git a/docs/usage/faster-lambdify.ipynb b/docs/usage/faster-lambdify.ipynb
index 9f8724f0..9ae9e598 100644
--- a/docs/usage/faster-lambdify.ipynb
+++ b/docs/usage/faster-lambdify.ipynb
@@ -69,11 +69,14 @@
     "from ampform.dynamics.builder import create_relativistic_breit_wigner_with_ff\n",
     "from IPython.display import HTML, SVG\n",
     "\n",
-    "from tensorwaves.model import LambdifiedFunction, SympyModel\n",
-    "from tensorwaves.model.sympy import optimized_lambdify, split_expression\n",
+    "from tensorwaves.function.sympy import (\n",
+    "    create_function,\n",
+    "    fast_lambdify,\n",
+    "    split_expression,\n",
+    ")\n",
     "\n",
-    "logger = logging.getLogger()\n",
-    "logger.setLevel(logging.ERROR)"
+    "LOGGER = logging.getLogger()\n",
+    "LOGGER.setLevel(logging.ERROR)"
    ]
   },
   {
@@ -97,7 +100,7 @@
    "outputs": [],
    "source": [
     "x, y, z = sp.symbols(\"x:z\")\n",
-    "expr = x ** z + 2 * y\n",
+    "expr = x ** z + 2 * y + sp.log(y * z)\n",
     "expr"
    ]
   },
@@ -138,15 +141,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "top_expr, sub_expressions = split_expression(expr, max_complexity=3)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "top_expr, sub_expressions = split_expression(expr, max_complexity=3)\n",
     "top_expr"
    ]
   },
@@ -238,14 +233,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Optimized lambdify"
+    "## Fast lambdify"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Generally, the lambdify time scales exponentially with the size of an expression tree. With larger expression trees, it's therefore much faster to lambdify these sub-expressions separately and to recombine them. TensorWaves offers a function that does this for you: {func}`.optimized_lambdify`. We'll use an {class}`~ampform.helicity.HelicityModel` to illustrate this:"
+    "Generally, the lambdify time scales exponentially with the size of an expression tree. With larger expression trees, it's therefore much faster to lambdify these sub-expressions separately and to recombine them. TensorWaves offers a function that does this for you: {func}`.fast_lambdify`. We'll use an {class}`~ampform.helicity.HelicityModel` to illustrate this:"
    ]
   },
   {
@@ -282,7 +277,7 @@
    "outputs": [],
    "source": [
     "%%time\n",
-    "lambdified_optimized = optimized_lambdify(\n",
+    "split_function = fast_lambdify(\n",
     "    expression,\n",
     "    sorted_symbols,\n",
     "    max_complexity=100,\n",
@@ -303,7 +298,7 @@
    },
    "outputs": [],
    "source": [
-    "lambdified_optimized"
+    "split_function"
    ]
   },
   {
@@ -315,41 +310,37 @@
    "outputs": [],
    "source": [
     "%%time\n",
-    "sp.lambdify(sorted_symbols, expression)"
+    "normal_function = sp.lambdify(sorted_symbols, expression)"
    ]
   },
   {
-   "cell_type": "markdown",
-   "metadata": {},
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "remove-cell"
+    ]
+   },
+   "outputs": [],
    "source": [
-    "## Specifying complexity"
+    "normal_function"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the usually workflow (see {doc}`/usage`), TensorWaves uses SymPy's own {func}`~sympy.utilities.lambdify.lambdify` by default. You can change this behavior with the `max_complexity` argument of {class}`.SympyModel`:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sympy_model = SympyModel(\n",
-    "    expression=model.expression.doit(),\n",
-    "    parameters=model.parameter_defaults,\n",
-    "    max_complexity=100,\n",
-    ")"
+    "## Specifying complexity"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If `max_complexity` is specified (i.e., is not {obj}`None`), {class}`.LambdifiedFunction` uses TensorWaves's {func}`.optimized_lambdify`."
+    "In the {doc}`usual workflow </usage>`, we used the {func}`.create_function` function. By default, this internally calls SymPy's own {func}`~sympy.utilities.lambdify.lambdify` function. But if you specify its `max_complexity` argument, {func}`.create_function` uses TensorWaves's {func}`.fast_lambdify`."
    ]
   },
   {
@@ -359,7 +350,12 @@
    "outputs": [],
    "source": [
     "%%time\n",
-    "intensity = LambdifiedFunction(sympy_model, backend=\"jax\")"
+    "function = create_function(\n",
+    "    expression=model.expression.doit(),\n",
+    "    parameters=model.parameter_defaults,\n",
+    "    max_complexity=100,\n",
+    "    backend=\"numpy\",\n",
+    ")"
    ]
   }
  ],
diff --git a/docs/usage/step1.ipynb b/docs/usage/step1.ipynb
index dd4cc1dc..5ed8b90e 100644
--- a/docs/usage/step1.ipynb
+++ b/docs/usage/step1.ipynb
@@ -51,9 +51,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "TensorWaves requires you to first formulate a mathematical {class}`.Model` that you want to fit to your data set. When this model is an amplitude model, it is most convenient to construct it with {mod}`ampform`.\n",
+    "Whether {doc}`generating data <step2>` or {doc}`fitting a model <step3>`, TensorWaves takes mathematical expressions as input. When that expression is an amplitude model, it is most convenient to formulate it with {mod}`ampform`.\n",
     "\n",
-    "This notebook briefly illustrates how to create such an amplitude model with {mod}`ampform` and how to write it to a recipe file that can be understood by {mod}`tensorwaves`. For more control, have a look at {doc}`the usage guides of AmpForm <ampform:usage>`."
+    "This notebook briefly illustrates how to create such an amplitude model with {mod}`ampform` and how to write it to a recipe file that can be understood by {mod}`tensorwaves`. For more advanced examples, have a look at {doc}`the usage guides of AmpForm <ampform:usage>`."
    ]
   },
   {
@@ -139,7 +139,7 @@
     "1. The coefficients for the different amplitudes are **complex** valued.\n",
     "2. By default there is no dynamics in the model, so it still has to be specified.\n",
     "\n",
-    "We choose to use {func}`~ampform.dynamics.relativistic_breit_wigner_with_ff` as the lineshape for all resonances and use a Blatt-Weisskopf form factor (no dynamics; see {class}`~ampform.dynamics.BlattWeisskopfSquared`) for the production decay. The {meth}`~ampform.helicity.HelicityAmplitudeBuilder.set_dynamics` is a convenience interface for replacing the dynamics for intermediate states."
+    "We choose to use {func}`~ampform.dynamics.relativistic_breit_wigner_with_ff` as the lineshape for all resonances and use a Blatt-Weisskopf form factor ({func}`~ampform.dynamics.builder.create_non_dynamic_with_ff`) for the production decay. The {meth}`~ampform.helicity.HelicityAmplitudeBuilder.set_dynamics` is a convenience interface for replacing the dynamics for intermediate states."
    ]
   },
   {
@@ -183,7 +183,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Finally, we can write the {class}`~ampform.helicity.HelicityModel` to disk via {mod}`pickle`, as well as store the {class}`~qrules.transition.ReactionInfo` as JSON:"
+    "Finally, we can write the {class}`~ampform.helicity.HelicityModel` to disk via {mod}`pickle`. The {class}`~qrules.transition.ReactionInfo` object can be pickled as well, but here, we write it to JSON:"
    ]
   },
   {
diff --git a/docs/usage/step2.ipynb b/docs/usage/step2.ipynb
index 0a610985..e79cbac4 100644
--- a/docs/usage/step2.ipynb
+++ b/docs/usage/step2.ipynb
@@ -44,7 +44,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Step 2: Generate data samples"
+    "# Step 2: Generate data"
    ]
   },
   {
@@ -64,8 +64,10 @@
    "source": [
     "import pickle\n",
     "\n",
+    "from ampform.helicity import HelicityModel\n",
+    "\n",
     "with open(\"helicity_model.pickle\", \"rb\") as model_file:\n",
-    "    model = pickle.load(model_file)"
+    "    model: HelicityModel = pickle.load(model_file)"
    ]
   },
   {
@@ -108,7 +110,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "As opposed to the main {ref}`usage:Generate data sample` of the main usage example page, we will generate a **deterministic** data sample. This can be done by feeding a {class}`.UniformRealNumberGenerator` with a specific {attr}`.UniformRealNumberGenerator.seed` and feeding that generator to the {func}`.generate_phsp` function:"
+    "As opposed to the main {ref}`usage:Generate data sample` of the main usage example page, we will generate a **deterministic** data sample. This can be done by feeding a {class}`.UniformRealNumberGenerator` with a specific {attr}`~.UniformRealNumberGenerator.seed` and feeding that generator to the {func}`.generate_phsp` function:"
    ]
   },
   {
@@ -169,9 +171,7 @@
     "\n",
     "An intensity-based sample is generated with the function {func}`.generate_data`. Its usage is similar to {func}`.generate_phsp`, but now you have to provide a {obj}`.Function` as well as a {obj}`.DataTransformer` that is used to transform the four-momentum phase space sample to a data sample that can be understood by the {obj}`.Function`.\n",
     "\n",
-    "Now, recall that in {doc}`step1`, we used the helicity formalism to mathematically express the reaction in terms of an amplitude model. TensorWaves needs to convert this {obj}`~ampform.helicity.HelicityModel` to an {class}`.Model` object that it can then {meth}`~.Model.lambdify` to a {obj}`.Function` object.\n",
-    "\n",
-    "The {obj}`~ampform.helicity.HelicityModel`  was expressed in terms of {mod}`sympy`, so we express the model as a {class}`.SympyModel` and lambdify it to a {class}`.LambdifiedFunction`:"
+    "Now, recall that in {doc}`step1`, we used the helicity formalism to mathematically express the reaction in terms of an amplitude model. TensorWaves needs to convert this {obj}`~ampform.helicity.HelicityModel` to a {obj}`.Function` object that can perform fast computations. This can be done with {func}`.create_function`:"
    ]
   },
   {
@@ -180,7 +180,7 @@
    "source": [
     ":::{margin}\n",
     "\n",
-    "Here, we make use of {func}`.optimized_lambdify` by specifying `max_complexity`. See {ref}`usage/faster-lambdify:Specifying complexity`.\n",
+    "Here, we make use of {func}`.fast_lambdify` by specifying `max_complexity`. See {ref}`usage/faster-lambdify:Specifying complexity`.\n",
     "\n",
     ":::"
    ]
@@ -195,14 +195,14 @@
    },
    "outputs": [],
    "source": [
-    "from tensorwaves.model import LambdifiedFunction, SympyModel\n",
+    "from tensorwaves.function.sympy import create_function\n",
     "\n",
-    "sympy_model = SympyModel(\n",
+    "intensity = create_function(\n",
     "    expression=model.expression.doit(),\n",
     "    parameters=model.parameter_defaults,\n",
     "    max_complexity=200,\n",
-    ")\n",
-    "%time intensity = LambdifiedFunction(sympy_model, backend=\"numpy\")"
+    "    backend=\"numpy\",\n",
+    ")"
    ]
   },
   {
@@ -220,7 +220,7 @@
    "source": [
     "from tensorwaves.data.transform import HelicityTransformer\n",
     "\n",
-    "data_converter = HelicityTransformer(model.adapter)"
+    "helicity_transformer = HelicityTransformer(model.adapter)"
    ]
   },
   {
@@ -246,7 +246,7 @@
     "    size=10_000,\n",
     "    initial_state_mass=initial_state_mass,\n",
     "    final_state_masses=final_state_masses,\n",
-    "    data_transformer=data_converter,\n",
+    "    data_transformer=helicity_transformer,\n",
     "    intensity=intensity,\n",
     "    random_generator=rng,\n",
     ")\n",
@@ -263,7 +263,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "As before, we use a {class}`.UniformRealNumberGenerator` with a specific {attr}`.UniformRealNumberGenerator.seed` to ensure we get a **deterministic** data sample."
+    "As before, we use a {class}`.UniformRealNumberGenerator` with a specific {attr}`~.UniformRealNumberGenerator.seed` to ensure we get a **deterministic** data sample."
    ]
   },
   {
@@ -286,8 +286,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "phsp_set = data_converter.transform(phsp_sample)\n",
-    "data_set = data_converter.transform(data_sample)\n",
+    "phsp_set = helicity_transformer.transform(phsp_sample)\n",
+    "data_set = helicity_transformer.transform(data_sample)\n",
     "list(data_set)"
    ]
   },
@@ -308,10 +308,7 @@
    },
    "outputs": [],
    "source": [
-    "for (\n",
-    "    state_id,\n",
-    "    particle,\n",
-    ") in model.adapter.reaction_info.final_state.items():\n",
+    "for state_id, particle in reaction_info.final_state.items():\n",
     "    print(f\"ID {state_id}:\", particle.name)"
    ]
   },
@@ -333,7 +330,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The {obj}`~ampform.data.DataSet` can easily be converted it to a {class}`pandas.DataFrame`:"
+    "The {obj}`~ampform.data.DataSet` can easily be converted to a {class}`pandas.DataFrame`:"
    ]
   },
   {
@@ -364,14 +361,16 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "tags": []
+    "tags": [
+     "hide-cell"
+    ]
    },
    "outputs": [],
    "source": [
     "from matplotlib import cm\n",
     "\n",
     "reaction_info = model.adapter.reaction_info\n",
-    "intermediate_states = sorted(\n",
+    "resonances = sorted(\n",
     "    (\n",
     "        p\n",
     "        for p in model.particles\n",
@@ -381,7 +380,7 @@
     "    key=lambda p: p.mass,\n",
     ")\n",
     "\n",
-    "evenly_spaced_interval = np.linspace(0, 1, len(intermediate_states))\n",
+    "evenly_spaced_interval = np.linspace(0, 1, len(resonances))\n",
     "colors = [cm.rainbow(x) for x in evenly_spaced_interval]"
    ]
   },
@@ -399,9 +398,10 @@
     "\n",
     "data_frame[\"m_12\"].hist(bins=100, alpha=0.5, density=True, figsize=(9, 4))\n",
     "plt.xlabel(\"$m$ [GeV]\")\n",
-    "for i, p in enumerate(intermediate_states):\n",
-    "    plt.axvline(x=p.mass, linestyle=\"dotted\", label=p.name, color=colors[i])\n",
-    "plt.legend();"
+    "for p, color in zip(resonances, colors):\n",
+    "    plt.axvline(x=p.mass, linestyle=\"dotted\", label=p.name, color=color)\n",
+    "plt.legend()\n",
+    "plt.show()"
    ]
   },
   {
@@ -426,7 +426,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "TensorWaves currently has no export functionality for data samples, so we just {func}`pickle.dump` these data samples as follows:"
+    "To export the generated data samples, simply {func}`pickle.dump` them as follows:"
    ]
   },
   {
diff --git a/docs/usage/step3.ipynb b/docs/usage/step3.ipynb
index d34bda30..55255573 100644
--- a/docs/usage/step3.ipynb
+++ b/docs/usage/step3.ipynb
@@ -51,9 +51,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "As explained in the {doc}`previous step <step2>`, a {class}`.Model` object and its lambdified {class}`.Function` form behave just like a mathematical function that takes a set of data points as an argument and returns a list of intensities (real numbers). At this stage, we want to optimize the parameters of the intensity model, so that it matches the distribution of our data sample. This is what we call 'performing a fit'.\n",
+    "As explained in the {doc}`previous step <step2>`, a {class}`.Function` can compute a list of intensities (real numbers) for an input `.DataSample`. At this stage, we want to optimize the parameters of this {class}`.Function`, so that it matches the distribution of our data sample. This is what we call 'fitting'.\n",
     "\n",
-    "First, load the relevant data from the previous steps."
+    "First, we load the relevant data from the previous steps."
    ]
   },
   {
@@ -64,7 +64,8 @@
      "source_hidden": true
     },
     "tags": [
-     "hide-input"
+     "hide-input",
+     "remove-output"
     ]
    },
    "outputs": [],
@@ -72,19 +73,23 @@
     "import pickle\n",
     "\n",
     "import qrules\n",
+    "from ampform.helicity import HelicityModel\n",
     "\n",
-    "from tensorwaves.model import SympyModel\n",
+    "from tensorwaves.function.sympy import create_function\n",
     "\n",
     "reaction = qrules.io.load(\"transitions.json\")\n",
     "with open(\"helicity_model.pickle\", \"rb\") as stream:\n",
-    "    model = pickle.load(stream)\n",
+    "    model: HelicityModel = pickle.load(stream)\n",
     "with open(\"data_set.pickle\", \"rb\") as stream:\n",
     "    data_set = pickle.load(stream)\n",
     "with open(\"phsp_set.pickle\", \"rb\") as stream:\n",
     "    phsp_set = pickle.load(stream)\n",
-    "sympy_model = SympyModel(\n",
+    "\n",
+    "function = create_function(\n",
     "    expression=model.expression.doit(),\n",
     "    parameters=model.parameter_defaults,\n",
+    "    backend=\"jax\",\n",
+    "    max_complexity=100,\n",
     ")"
    ]
   },
@@ -99,12 +104,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "To perform a fit, you need to define an {class}`.Estimator`. An estimator is a measure for the discrepancy between the intensity model and the data distribution to which you fit. In PWA, we usually use an *unbinned negative log likelihood estimator*.\n",
+    "To perform a fit, you need to define an {class}`.Estimator`. This is a measure for the discrepancy between the {class}`.Function` and the data distribution to which you fit it. In PWA, we usually use an **unbinned negative log likelihood estimator** ({class}`.UnbinnedNLL`).\n",
     "\n",
-    "Generally, the intensity model is not normalized, but a log likelihood estimator requires a normalized function. This is where the {ref}`phase space dataset <usage/step2:2.1 Generate phase space sample>` comes into play again: the intensity is evaluated separately with the phase space dataset so that the output of the {class}`.Function` can be normalized with it. The phase space sample is therefore a required argument!\n",
+    "Generally, the {class}`.Function` is not normalized with regards to the data sample, while a log likelihood estimator requires a normalized function. This is where the {ref}`phase space dataset <usage/step2:2.1 Generate phase space sample>` comes into play again: the {class}`.Function` is evaluated over the phase space dataset so that its output can be used as a normalization factor.\n",
     "\n",
     "```{margin}\n",
-    "If you want to correct for the efficiency of the detector, you should insert a *detector-reconstructed* phase space sample here.\n",
+    "If you want to correct for the efficiency of the detector, you should use a *detector-reconstructed* phase space sample.\n",
     "```"
    ]
   },
@@ -117,9 +122,9 @@
     "from tensorwaves.estimator import UnbinnedNLL\n",
     "\n",
     "estimator = UnbinnedNLL(\n",
-    "    sympy_model,\n",
-    "    data_set,\n",
-    "    phsp_set,\n",
+    "    function,\n",
+    "    data=data_set,\n",
+    "    phsp=phsp_set,\n",
     "    backend=\"jax\",\n",
     ")"
    ]
@@ -128,7 +133,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Note that the {class}`.UnbinnedNLL` can be expressed with different backends (it creates a {class}`.LambdifiedFunction` internally from the template {class}`.Model`). Here, we use {func}`jax <jax.jit>`, which is turns out to be the fastest backend for this model."
+    "Note that the {class}`.UnbinnedNLL` can be expressed with different backends, because it uses statistical operations like `log` and `mean`. Here, we use {func}`jax <jax.jit>`, which turns out to be the fastest backend for this model."
    ]
   },
   {
@@ -151,9 +156,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Starting the fit itself is quite simple: just create an optimizer instance of your choice, here [Minuit2](https://root.cern.ch/doc/master/Minuit2Page.html), and call its {meth}`~.Minuit2.optimize` method to start the fitting process. The {meth}`~.Minuit2.optimize` method requires a mapping of parameter names to their initial values. **Only the parameters listed in the mapping are optimized.**\n",
+    "Starting the fit itself is quite simple: just create an {mod}`.optimizer` instance of your choice, here [Minuit2](https://root.cern.ch/doc/master/Minuit2Page.html), and call its {meth}`~.Optimizer.optimize` method to start the fitting process. The {meth}`~.Optimizer.optimize` method requires a mapping of parameter names to their initial values. **Only the parameters listed in the mapping are optimized.**\n",
     "\n",
-    "Let's first select a few of the parameters that we saw in {ref}`Step 3.1 <usage/step3:3.1 Define estimator>` and feed them to the optimizer to run the fit. Notice that we modify the parameters slightly to make the fit more interesting (we are running fitting to a data sample that was generated with this very same amplitude model after all)."
+    "Let's first select a few of parameters and give them a small offset with regard to their original value. We give this offset to make the fit more interesting―after all, we are fitting to a data sample that was generated with this very same `.Function`."
    ]
   },
   {
@@ -166,6 +171,7 @@
    },
    "outputs": [],
    "source": [
+    "original_parameters = function.parameters\n",
     "initial_parameters = {\n",
     "    R\"C_{J/\\psi(1S) \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500) \\to \\pi^{0}_{0} \\pi^{0}_{0}}\": 1.0\n",
     "    + 0.0j,\n",
@@ -179,7 +185,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Recall that a {class}`.Function` object computes the intensity for a certain dataset. This can be seen now nicely when we use these intensities as weights on the phase space sample and plot it together with the original dataset. Here, we look at the invariant mass distribution projection of the final states `1` and `2`, which, {ref}`as we saw before <usage/step2:2.3 Visualize kinematic variables>`, are the final state particles $\\pi^0,\\pi^0$.\n",
+    "Recall that a {class}`.Function` object computes the intensity for a certain dataset. This can be seen now nicely when we use these intensities as weights on the phase space sample and plot it together with the original dataset. Here, we look at the invariant mass distribution projection of the final states `1` and `2`, which, {ref}`as we saw before <usage/step2:2.3 Visualize kinematic variables>`, is the final state particle pair $\\pi^0\\pi^0$.\n",
     "\n",
     "Don't forget to use {meth}`~.Function.update_parameters` first!"
    ]
@@ -202,7 +208,7 @@
     "from matplotlib import cm\n",
     "\n",
     "reaction_info = model.adapter.reaction_info\n",
-    "intermediate_states = sorted(\n",
+    "resonances = sorted(\n",
     "    (\n",
     "        p\n",
     "        for p in model.particles\n",
@@ -212,15 +218,18 @@
     "    key=lambda p: p.mass,\n",
     ")\n",
     "\n",
-    "evenly_spaced_interval = np.linspace(0, 1, len(intermediate_states))\n",
+    "evenly_spaced_interval = np.linspace(0, 1, len(resonances))\n",
     "colors = [cm.rainbow(x) for x in evenly_spaced_interval]\n",
     "\n",
     "\n",
     "def indicate_masses():\n",
     "    plt.xlabel(\"$m$ [GeV]\")\n",
-    "    for i, p in enumerate(intermediate_states):\n",
+    "    for color, resonance in zip(colors, resonances):\n",
     "        plt.gca().axvline(\n",
-    "            x=p.mass, linestyle=\"dotted\", label=p.name, color=colors[i]\n",
+    "            x=resonance.mass,\n",
+    "            linestyle=\"dotted\",\n",
+    "            label=resonance.name,\n",
+    "            color=color,\n",
     "        )\n",
     "\n",
     "\n",
@@ -228,22 +237,22 @@
     "    variable_name,\n",
     "    data_set,\n",
     "    phsp_set,\n",
-    "    intensity_model,\n",
+    "    function,\n",
     "    bins=100,\n",
     "):\n",
     "    data = data_set[variable_name]\n",
     "    phsp = phsp_set[variable_name]\n",
-    "    intensities = intensity_model(phsp_set)\n",
-    "    fig, ax = plt.subplots(figsize=(9, 4))\n",
-    "    axis = plt.gca()\n",
-    "    axis.hist(\n",
+    "    intensities = function(phsp_set)\n",
+    "    _, ax = plt.subplots(figsize=(9, 4))\n",
+    "    ax = plt.gca()\n",
+    "    ax.hist(\n",
     "        data,\n",
     "        bins=bins,\n",
     "        alpha=0.5,\n",
     "        label=\"data\",\n",
     "        density=True,\n",
     "    )\n",
-    "    axis.hist(\n",
+    "    ax.hist(\n",
     "        phsp,\n",
     "        weights=intensities,\n",
     "        bins=bins,\n",
@@ -253,18 +262,7 @@
     "        density=True,\n",
     "    )\n",
     "    indicate_masses()\n",
-    "    plt.legend()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from tensorwaves.model import LambdifiedFunction\n",
-    "\n",
-    "intensity = LambdifiedFunction(sympy_model, backend=\"numpy\")"
+    "    ax.legend()"
    ]
   },
   {
@@ -277,8 +275,8 @@
    },
    "outputs": [],
    "source": [
-    "intensity.update_parameters(initial_parameters)\n",
-    "compare_model(\"m_12\", data_set, phsp_set, intensity)"
+    "function.update_parameters(initial_parameters)\n",
+    "compare_model(\"m_12\", data_set, phsp_set, function)"
    ]
   },
   {
@@ -310,27 +308,27 @@
    },
    "outputs": [],
    "source": [
-    "import os\n",
-    "\n",
     "from IPython.display import clear_output\n",
     "\n",
     "from tensorwaves.optimizer.callbacks import Callback\n",
     "\n",
-    "variable_to_plot = \"m_12\"\n",
-    "is_doc = \"READTHEDOCS\" in os.environ or \"EXECUTE_NB\" in os.environ\n",
-    "\n",
     "\n",
     "class PyplotCallback(Callback):\n",
-    "    def __init__(self, step_size=10):\n",
+    "    def __init__(self, variable=\"m_12\", step_size=10):\n",
+    "        self.__variable = variable\n",
     "        self.__step_size = step_size\n",
     "        self.__fig = None\n",
     "        self.__ax = None\n",
     "        self.__latest_parameters = {}\n",
     "\n",
     "    def on_optimize_start(self, logs):\n",
+    "        if STATIC_WEB_PAGE:\n",
+    "            return\n",
     "        self.__fig, self.__ax = plt.subplots(1, figsize=(8, 5))\n",
     "\n",
     "    def on_optimize_end(self, logs):\n",
+    "        if STATIC_WEB_PAGE:\n",
+    "            return\n",
     "        self.update_plot(nbins=200)\n",
     "        self.__ax = None\n",
     "        self.__fig = None\n",
@@ -339,9 +337,9 @@
     "        pass\n",
     "\n",
     "    def on_function_call_end(self, function_call, logs):\n",
-    "        self.__latest_parameters = logs[\"parameters\"]\n",
-    "        if is_doc:\n",
+    "        if STATIC_WEB_PAGE:\n",
     "            return\n",
+    "        self.__latest_parameters = logs[\"parameters\"]\n",
     "        if function_call % self.__step_size != 0:\n",
     "            return\n",
     "        if function_call > 75:\n",
@@ -353,10 +351,10 @@
     "    def update_plot(self, nbins: int):\n",
     "        if self.__fig is None or self.__ax is None:\n",
     "            self.__fig, self.__ax = plt.subplots(1, figsize=(8, 5))\n",
-    "        data = data_set[variable_to_plot]\n",
-    "        phsp = phsp_set[variable_to_plot]\n",
-    "        intensity.update_parameters(self.__latest_parameters)\n",
-    "        intensities = intensity(phsp_set)\n",
+    "        data = data_set[self.__variable]\n",
+    "        phsp = phsp_set[self.__variable]\n",
+    "        function.update_parameters(self.__latest_parameters)\n",
+    "        intensities = function(phsp_set)\n",
     "        self.__ax.cla()\n",
     "        self.__ax.hist(data, bins=nbins, alpha=0.5, label=\"data\", density=True)\n",
     "        self.__ax.hist(\n",
@@ -417,10 +415,10 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Finally, we create an {class}`.Optimizer` to {meth}`~.Minuit2.optimize` the model (which is embedded in the {class}`.Estimator`). Here, we choose the {class}`.Minuit2` optimizer, which is the most common optimizer in high-energy physics. Since we constructed the {class}`.UnbinnedNLL` with {obj}`jax <jax.jit>`, can an optimize the model with analytic gradient over the {class}`.LambdifiedFunction`:\n",
+    "Finally, we create an {class}`.Optimizer` to {meth}`~.Minuit2.optimize` the parameters in the {class}`.Function` (which is embedded in the {class}`.Estimator`). Here, we choose the {class}`.Minuit2` optimizer, which is the most common optimizer in high-energy physics. Since we constructed the {class}`.UnbinnedNLL` with {obj}`jax <jax.jit>`, can an optimize the model with analytic gradient over the {class}`.LambdifiedFunction`:\n",
     "\n",
     "```{margin}\n",
-    "The computation time depends on the complexity of the model, number of data events, the size of the phase space sample, and the number of free parameters. This model is rather small and has but a few free parameters, so the optimization shouldn't take more than a **few seconds**.\n",
+    "The computation time depends on the complexity of the model, the number of data events, the size of the phase space sample, and the number of free parameters. This model is rather small and has but a few free parameters, so the optimization shouldn't take more than a minute.\n",
     "\n",
     "The custom callback example slows down the optimization, but you can remove it to get the raw fit performance.\n",
     "```"
@@ -466,7 +464,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "As can be seen, the values of the optimized parameters in the result are again comparable to the original values that are contained in the model ({attr}`.SympyModel.parameters`):"
+    "As can be seen, the values of the optimized parameters in the {class}`.FitResult` are again comparable to the original parameter values."
    ]
   },
   {
@@ -480,7 +478,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Each of the {mod}`.optimizer`s offer more specific information about the fit result. This information can be accessed with {attr}`.FitResult.specifics`. A common example would be to get the `~iminuit.Minuit.covariance` matrix:"
+    "Each of the {mod}`.optimizer`s offer more specific information about the fit result. This information can be accessed with {attr}`.FitResult.specifics`. A common example would be to get the {attr}`~iminuit.Minuit.covariance` matrix:"
    ]
   },
   {
@@ -490,7 +488,7 @@
    "outputs": [],
    "source": [
     "covariance_matrix = fit_result.specifics.covariance\n",
-    "covariance_matrix"
+    "covariance_matrix.correlation()"
    ]
   },
   {
@@ -557,7 +555,7 @@
     "    print(p)\n",
     "    print(f\"  initial:   {initial_parameters[p]:.3}\")\n",
     "    print(f\"  optimized: {optimized_parameters[p]:.3}\")\n",
-    "    print(f\"  original:  {sympy_model.parameters[p]:.3}\")"
+    "    print(f\"  original:  {original_parameters[p]:.3}\")"
    ]
   },
   {
@@ -584,7 +582,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In {ref}`usage/step3:3.2 Optimize intensity model`, we initialized {obj}`.Minuit2` with some callbacks that are {class}`.Loadable`. Such callback classes offer the possibility to {meth}`.Loadable.load_latest_parameters`, so you can pick up the optimize process in case it crashes or if you pause it. Loading the latest parameters goes as follows:"
+    "In {ref}`usage/step3:3.2 Optimize intensity model`, we initialized {obj}`.Minuit2` with some callbacks that are {class}`.Loadable`. Such callback classes offer the possibility to {meth}`~.Loadable.load_latest_parameters`, so you can pick up the optimize process in case it crashes or if you pause it. Loading the latest parameters goes as follows:"
    ]
   },
   {
@@ -658,7 +656,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Using the same method as above, we renew the parameters of the {class}`.Function` and plot it again over the data distribution."
+    "Using the same method as above, we renew the parameters of the {class}`.Function` and plot it again over the phase space sample."
    ]
   },
   {
@@ -671,8 +669,8 @@
    },
    "outputs": [],
    "source": [
-    "intensity.update_parameters(latest_parameters)\n",
-    "compare_model(\"m_12\", data_set, phsp_set, intensity)"
+    "function.update_parameters(latest_parameters)\n",
+    "compare_model(\"m_12\", data_set, phsp_set, function)"
    ]
   },
   {
@@ -722,7 +720,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Just like in {ref}`usage/step2:2.2 Generate intensity-based sample`, these _intensity components_ can each be expressed in a computational backend. This can be done with the method {meth}`~ampform.helicity.HelicityModel.sum_components` from the original model and creating a new {class}`.LambdifiedFunction` from that expression. After that we, update the parameters with the optimized parameter values we found in {ref}`usage/step3:Perform fit`. The two components in the example below should be the same:"
+    "Just like in {ref}`usage/step2:2.2 Generate intensity-based sample`, these _intensity components_ can each be expressed in a computational backend. This can be done with the method {meth}`~ampform.helicity.HelicityModel.sum_components` from the original model and creating a new {class}`.LambdifiedFunction` from that expression with {func}`.create_function`. After that we, update the parameters with the optimized parameter values we found in {ref}`usage/step3:Perform fit`. The two components in the example below should be the same:"
    ]
   },
   {
@@ -749,11 +747,9 @@
     "        R\" \\to \\pi^{0}_{0} \\pi^{0}_{0}}\",\n",
     "    ]\n",
     ")\n",
-    "from_amplitudes = LambdifiedFunction(\n",
-    "    model=SympyModel(\n",
-    "        expression=added_components.doit(),\n",
-    "        parameters=model.parameter_defaults,\n",
-    "    ),\n",
+    "from_amplitudes = create_function(\n",
+    "    expression=added_components.doit(),\n",
+    "    parameters=model.parameter_defaults,\n",
     "    backend=\"numpy\",\n",
     ")\n",
     "from_amplitudes.update_parameters(latest_parameters)"
@@ -768,11 +764,9 @@
     "added_components = model.sum_components(\n",
     "    components=[R\"I_{J/\\psi(1S)_{+1} \\to \\gamma_{+1} \\pi^{0}_{0} \\pi^{0}_{0}}\"]\n",
     ")\n",
-    "from_intensity = LambdifiedFunction(\n",
-    "    model=SympyModel(\n",
-    "        expression=added_components.doit(),\n",
-    "        parameters=model.parameter_defaults,\n",
-    "    ),\n",
+    "from_intensity = create_function(\n",
+    "    expression=added_components.doit(),\n",
+    "    parameters=model.parameter_defaults,\n",
     "    backend=\"numpy\",\n",
     ")\n",
     "from_intensity.update_parameters(latest_parameters)"
@@ -814,7 +808,7 @@
     "bins = 150\n",
     "ax.hist(\n",
     "    phsp_set[\"m_12\"],\n",
-    "    weights=intensity(phsp_set),\n",
+    "    weights=function(phsp_set),\n",
     "    bins=bins,\n",
     "    alpha=0.2,\n",
     "    label=\"full intensity\",\n",
@@ -829,7 +823,8 @@
     "ax.set_xlim(0.25, 2.5)\n",
     "ax.set_xlabel(R\"$m_{\\pi^0\\pi^0}$ [GeV]\")\n",
     "ax.set_yticks([])\n",
-    "plt.legend();"
+    "plt.legend()\n",
+    "plt.show()"
    ]
   },
   {
@@ -859,11 +854,9 @@
     "logging.basicConfig()\n",
     "logging.getLogger().setLevel(logging.ERROR)\n",
     "intensity_components = [\n",
-    "    LambdifiedFunction(\n",
-    "        SympyModel(\n",
-    "            expression=model.sum_components([c]).doit(),\n",
-    "            parameters=model.parameter_defaults,\n",
-    "        ),\n",
+    "    create_function(\n",
+    "        expression=model.sum_components([c]).doit(),\n",
+    "        parameters=model.parameter_defaults,\n",
     "        backend=\"numpy\",\n",
     "    )\n",
     "    for c in model.components\n",
@@ -895,7 +888,7 @@
     "ax.set_yticks([])\n",
     "ax.legend(\n",
     "    labels=[f\"${c[3:-1]}$\" for c in model.components if c.startswith(\"I\")]\n",
-    ");"
+    ")"
    ]
   },
   {
@@ -989,7 +982,7 @@
     "    label = line.get_label()\n",
     "    color = line.get_color()\n",
     "    ax2.axhline(\n",
-    "        y=complex(sympy_model.parameters[label]).real,\n",
+    "        y=complex(original_parameters[label]).real,\n",
     "        color=color,\n",
     "        alpha=0.5,\n",
     "        linestyle=\"dotted\",\n",
@@ -1072,7 +1065,7 @@
     "    label = line.get_label()\n",
     "    color = line.get_color()\n",
     "    ax2.axhline(\n",
-    "        y=complex(sympy_model.parameters[label]).real,\n",
+    "        y=complex(original_parameters[label]).real,\n",
     "        color=color,\n",
     "        alpha=0.5,\n",
     "        linestyle=\"dotted\",\n",
diff --git a/src/tensorwaves/__init__.py b/src/tensorwaves/__init__.py
index 4abb1c9e..b5b537db 100644
--- a/src/tensorwaves/__init__.py
+++ b/src/tensorwaves/__init__.py
@@ -2,17 +2,18 @@
 
 The `tensorwaves` package contains four main ingredients:
 
-  Model lambdifying (`tensorwaves.model`)
-    ― Convert arbitrary mathematical expressions to an computational backends
-    like `numpy`, `jax <jax.jit>`, and `tensorflow <tf.Tensor>`. This module
-    works best with `sympy` expressions as input.
+  Function creation (`tensorwaves.function`)
+    ― Express arbitrary mathematical expressions as functions in different
+    kinds of computational backends, like `numpy`, `jax <jax.jit>`, and
+    `tensorflow <tf.Tensor>`. This module works best with `sympy` expressions
+    as input.
 
   Data generation (`tensorwaves.data`)
-    ― Generate phase space samples and toy Monte Carlo samples based on an
-    intensity `.Model`.
+    ― Generate phase space samples as well as hit-and-miss Monte Carlo samples
+    for the input mathematical expression.
 
-  Model optimization (`tensorwaves.optimizer` and `~.estimator`)
-    ― Optimize the `.Model` with respect to some data sample and an
+  Optimization (`tensorwaves.optimizer` and `~.estimator`)
+    ― Optimize the `.Function` with respect to some data sample and an
     `.Estimator` (loss function).
 
 The `.interface` module defines how the main classes interact.
@@ -21,8 +22,8 @@
 __all__ = [
     "data",
     "estimator",
-    "model",
+    "function",
     "optimizer",
 ]
 
-from . import data, estimator, model, optimizer
+from . import data, estimator, function, optimizer
diff --git a/src/tensorwaves/estimator.py b/src/tensorwaves/estimator.py
index c84a8865..4bff6952 100644
--- a/src/tensorwaves/estimator.py
+++ b/src/tensorwaves/estimator.py
@@ -2,7 +2,7 @@
 
 All estimators have to implement the `.Estimator` interface.
 """
-from typing import Callable, Dict, Mapping, Optional, Union
+from typing import Callable, Dict, Mapping, Union
 
 import numpy as np
 
@@ -11,10 +11,8 @@
     DataSample,
     Estimator,
     Function,
-    Model,
     ParameterValue,
 )
-from tensorwaves.model import LambdifiedFunction
 
 
 def gradient_creator(
@@ -44,54 +42,30 @@ class UnbinnedNLL(Estimator):  # pylint: disable=too-many-instance-attributes
     """Unbinned negative log likelihood estimator.
 
     Args:
-        model: A model that should be compared to the dataset.
-        dataset: The dataset used for the comparison. The model has to be
+        function: A `.Function` that describes a distribution over a certain
+            domain.
+        data: The dataset used for the comparison. The model has to be
             evaluateable with this dataset.
-        phsp_set: A phase space dataset, which is used for the normalization.
-            The model has to be evaluateable with this dataset. When correcting
-            for the detector efficiency use a phase space sample, that passed
-            the detector reconstruction.
+        phsp: The domain (phase space) over which to execute the function is
+            used for the normalization. When correcting for the detector
+            efficiency, use a phase space sample that passed the detector
+            reconstruction.
+        backend: The computational back-end with which the negative log
+            likelihood should be computed.
 
     """
 
     def __init__(  # pylint: disable=too-many-arguments
         self,
-        model: Union[Function, Model],
-        dataset: DataSample,
-        phsp_dataset: DataSample,
+        function: Function,
+        data: DataSample,
+        phsp: DataSample,
         phsp_volume: float = 1.0,
         backend: Union[str, tuple, dict] = "numpy",
-        use_caching: bool = False,
-        fixed_parameters: Optional[Dict[str, ParameterValue]] = None,
     ) -> None:
-        self.__use_caching = use_caching
-        self.__dataset = {k: np.array(v) for k, v in dataset.items()}
-        self.__phsp_dataset = {k: np.array(v) for k, v in phsp_dataset.items()}
-        if isinstance(model, Function):
-            self.__data_function = model
-            self.__phsp_function = model
-        elif isinstance(model, Model):
-            if self.__use_caching:
-                fixed_data_inputs = dict(self.__dataset)
-                fixed_phsp_inputs = dict(self.__phsp_dataset)
-                if fixed_parameters:
-                    fixed_data_inputs.update(fixed_parameters)  # type: ignore[arg-type]
-                    fixed_phsp_inputs.update(fixed_parameters)  # type: ignore[arg-type]
-                self.__data_function = LambdifiedFunction(
-                    model.performance_optimize(fix_inputs=fixed_data_inputs),
-                    backend,
-                )
-                self.__phsp_function = LambdifiedFunction(
-                    model.performance_optimize(fix_inputs=fixed_phsp_inputs),
-                    backend,
-                )
-            else:
-                self.__data_function = LambdifiedFunction(model, backend)
-                self.__phsp_function = self.__data_function
-        else:
-            raise TypeError(
-                f"{model.__class__} not of type {Function} or {Model}"
-            )
+        self.__data = {k: np.array(v) for k, v in data.items()}
+        self.__phsp = {k: np.array(v) for k, v in phsp.items()}
+        self.__function = function
         self.__gradient = gradient_creator(self.__call__, backend)
 
         self.__mean_function = find_function(backend, "mean")
@@ -101,14 +75,9 @@ def __init__(  # pylint: disable=too-many-arguments
         self.__phsp_volume = phsp_volume
 
     def __call__(self, parameters: Mapping[str, ParameterValue]) -> float:
-        self.__data_function.update_parameters(parameters)
-        if self.__use_caching:
-            self.__phsp_function.update_parameters(parameters)
-            bare_intensities = self.__data_function({})
-            phsp_intensities = self.__phsp_function({})
-        else:
-            bare_intensities = self.__data_function(self.__dataset)
-            phsp_intensities = self.__phsp_function(self.__phsp_dataset)
+        self.__function.update_parameters(parameters)
+        bare_intensities = self.__function(self.__data)
+        phsp_intensities = self.__function(self.__phsp)
         normalization_factor = 1.0 / (
             self.__phsp_volume * self.__mean_function(phsp_intensities)
         )
diff --git a/src/tensorwaves/model/__init__.py b/src/tensorwaves/function/__init__.py
similarity index 53%
rename from src/tensorwaves/model/__init__.py
rename to src/tensorwaves/function/__init__.py
index b909f60c..ea0a299a 100644
--- a/src/tensorwaves/model/__init__.py
+++ b/src/tensorwaves/function/__init__.py
@@ -1,44 +1,38 @@
-"""Evaluateable physics models for amplitude analysis.
+"""Express mathematical expressions in terms of computational functions."""
 
-The `.model` module takes care of lambdifying mathematical expressions to
-computational backends.
-"""
-
-# pyright: reportUnusedImport=false
-from typing import Dict, Mapping, Union
+from typing import Callable, Dict, Mapping, Sequence
 
 import numpy as np
 
-from tensorwaves.interface import DataSample, Function, Model, ParameterValue
-
-from .sympy import SympyModel  # noqa: F401
+from tensorwaves.interface import DataSample, Function, ParameterValue
 
 
 class LambdifiedFunction(Function):
-    """Implements `.Function` based on a `.Model` using {meth}`~.Model.lambdify`."""
+    """Implements `.Function` for a specific computational back-end."""
 
     def __init__(
         self,
-        model: Model,
-        backend: Union[str, tuple, dict] = "numpy",
+        function: Callable[..., np.ndarray],
+        argument_order: Sequence[str],
+        parameters: Mapping[str, ParameterValue],
     ) -> None:
-        self.__lambdified_model = model.lambdify(backend=backend)
-        self.__parameters = model.parameters
-        self.__ordered_args = model.argument_order
+        self.__function = function
+        self.__argument_order = argument_order
+        self.__parameters = dict(parameters)
 
     def __call__(self, dataset: DataSample) -> np.ndarray:
-        return self.__lambdified_model(
+        return self.__function(
             *[
                 dataset[var_name]
                 if var_name in dataset
                 else self.__parameters[var_name]
-                for var_name in self.__ordered_args
+                for var_name in self.__argument_order
             ],
         )
 
     @property
     def parameters(self) -> Dict[str, ParameterValue]:
-        return self.__parameters
+        return dict(self.__parameters)
 
     def update_parameters(
         self, new_parameters: Mapping[str, ParameterValue]
diff --git a/src/tensorwaves/function/sympy.py b/src/tensorwaves/function/sympy.py
new file mode 100644
index 00000000..89df12c7
--- /dev/null
+++ b/src/tensorwaves/function/sympy.py
@@ -0,0 +1,237 @@
+"""Lambdify `sympy` expression trees to a `.Function`."""
+
+import logging
+from typing import (
+    Any,
+    Callable,
+    Dict,
+    List,
+    Mapping,
+    Optional,
+    Sequence,
+    Tuple,
+    Union,
+)
+
+import sympy as sp
+from sympy.printing.numpy import (
+    NumPyPrinter,
+    _numpy_known_constants,
+    _numpy_known_functions,
+)
+from tqdm.auto import tqdm
+
+from tensorwaves._backend import get_backend_modules
+from tensorwaves.interface import ParameterValue
+
+from . import LambdifiedFunction
+
+_jax_known_functions = {
+    k: v.replace("numpy.", "jnp.") for k, v in _numpy_known_functions.items()
+}
+_jax_known_constants = {
+    k: v.replace("numpy.", "jnp.") for k, v in _numpy_known_constants.items()
+}
+
+
+class _JaxPrinter(NumPyPrinter):  # pylint: disable=abstract-method
+    # pylint: disable=invalid-name
+    module_imports = {"jax": {"numpy as jnp"}}
+    _module = "jnp"
+    _kc = _jax_known_constants
+    _kf = _jax_known_functions
+
+    def _print_ComplexSqrt(self, expr: sp.Expr) -> str:  # noqa: N802
+        x = self._print(expr.args[0])
+        return (
+            "jnp.select("
+            f"[jnp.less({x}, 0), True], "
+            f"[1j * jnp.sqrt(-{x}), jnp.sqrt({x})], "
+            "default=jnp.nan)"
+        )
+
+
+def split_expression(
+    expression: sp.Expr,
+    max_complexity: int,
+    min_complexity: int = 0,
+) -> Tuple[sp.Expr, Dict[sp.Symbol, sp.Expr]]:
+    """Split an expression into a 'top expression' and several sub-expressions.
+
+    Replace nodes in the expression tree of a `sympy.Expr
+    <sympy.core.expr.Expr>` that lie within a certain complexity range (see
+    :meth:`~sympy.core.basic.Basic.count_ops`) with symbols and keep a mapping
+    of each to these symbols to the sub-expressions that they replaced.
+
+    .. seealso:: :doc:`/usage/faster-lambdify`
+    """
+    i = 0
+    symbol_mapping: Dict[sp.Symbol, sp.Expr] = {}
+    n_operations = sp.count_ops(expression)
+    if n_operations < max_complexity:
+        return expression, symbol_mapping
+    progress_bar = tqdm(
+        total=n_operations,
+        desc="Splitting expression",
+        unit="node",
+        disable=not _use_progress_bar(),
+    )
+
+    def recursive_split(sub_expression: sp.Expr) -> sp.Expr:
+        nonlocal i
+        for arg in sub_expression.args:
+            complexity = sp.count_ops(arg)
+            if min_complexity < complexity < max_complexity:
+                progress_bar.update(n=complexity)
+                symbol = sp.Symbol(f"f{i}")
+                i += 1
+                symbol_mapping[symbol] = arg
+                sub_expression = sub_expression.xreplace({arg: symbol})
+            else:
+                new_arg = recursive_split(arg)
+                sub_expression = sub_expression.xreplace({arg: new_arg})
+        return sub_expression
+
+    top_expression = recursive_split(expression)
+    remainder = progress_bar.total - progress_bar.n
+    progress_bar.update(n=remainder)  # pylint crashes if total is set directly
+    progress_bar.close()
+    return top_expression, symbol_mapping
+
+
+def create_function(
+    expression: sp.Expr,
+    parameters: Mapping[sp.Symbol, ParameterValue],
+    backend: str,
+    max_complexity: Optional[int] = None,
+    **kwargs: Any,
+) -> LambdifiedFunction:
+    sorted_symbols = sorted(expression.free_symbols, key=lambda s: s.name)
+    if max_complexity is None:
+        lambdified_function = lambdify(
+            expression=expression,
+            symbols=sorted_symbols,
+            backend=backend,
+            **kwargs,
+        )
+    else:
+        lambdified_function = fast_lambdify(
+            expression=expression,
+            symbols=sorted_symbols,
+            backend=backend,
+            max_complexity=max_complexity,
+            **kwargs,
+        )
+    return LambdifiedFunction(
+        function=lambdified_function,
+        argument_order=list(map(str, sorted_symbols)),
+        parameters={
+            symbol.name: value for symbol, value in parameters.items()
+        },
+    )
+
+
+def fast_lambdify(
+    expression: sp.Expr,
+    symbols: Sequence[sp.Symbol],
+    backend: Union[str, tuple, dict],
+    *,
+    min_complexity: int = 0,
+    max_complexity: int,
+    **kwargs: Any,
+) -> Callable:
+    """Speed up :func:`.lambdify` with :func:`.split_expression`.
+
+    .. seealso:: :doc:`/usage/faster-lambdify`
+    """
+    top_expression, sub_expressions = split_expression(
+        expression,
+        min_complexity=min_complexity,
+        max_complexity=max_complexity,
+    )
+    if not sub_expressions:
+        return lambdify(top_expression, symbols, backend, **kwargs)
+
+    sorted_top_symbols = sorted(sub_expressions, key=lambda s: s.name)
+    top_function = lambdify(
+        top_expression, sorted_top_symbols, backend, **kwargs
+    )
+    sub_functions: List[Callable] = []
+    for symbol in tqdm(
+        iterable=sorted_top_symbols,
+        desc="Lambdifying sub-expressions",
+        unit="expr",
+        disable=not _use_progress_bar(),
+    ):
+        sub_expression = sub_expressions[symbol]
+        sub_function = lambdify(sub_expression, symbols, backend, **kwargs)
+        sub_functions.append(sub_function)
+
+    def recombined_function(*args: Any) -> Any:
+        new_args = [sub_function(*args) for sub_function in sub_functions]
+        return top_function(*new_args)
+
+    return recombined_function
+
+
+def lambdify(
+    expression: sp.Expr,
+    symbols: Sequence[sp.Symbol],
+    backend: Union[str, tuple, dict],
+    **kwargs: Any,
+) -> Callable:
+    """A wrapper around :func:`~sympy.utilities.lambdify.lambdify`."""
+    # pylint: disable=import-outside-toplevel,too-many-return-statements
+    def jax_lambdify() -> Callable:
+        import jax
+
+        return jax.jit(
+            sp.lambdify(
+                symbols,
+                expression,
+                modules=modules,
+                printer=_JaxPrinter,
+                **kwargs,
+            )
+        )
+
+    def numba_lambdify() -> Callable:
+        # pylint: disable=import-error
+        import numba
+
+        return numba.jit(
+            sp.lambdify(symbols, expression, modules="numpy", **kwargs),
+            forceobj=True,
+            parallel=True,
+        )
+
+    def tensorflow_lambdify() -> Callable:
+        # pylint: disable=import-error
+        import tensorflow.experimental.numpy as tnp  # pyright: reportMissingImports=false
+
+        return sp.lambdify(symbols, expression, modules=tnp, **kwargs)
+
+    modules = get_backend_modules(backend)
+    if isinstance(backend, str):
+        if backend == "jax":
+            return jax_lambdify()
+        if backend == "numba":
+            return numba_lambdify()
+        if backend in {"tensorflow", "tf"}:
+            return tensorflow_lambdify()
+
+    if isinstance(backend, tuple):
+        if any("jax" in x.__name__ for x in backend):
+            return jax_lambdify()
+        if any("numba" in x.__name__ for x in backend):
+            return numba_lambdify()
+        if any(
+            "tensorflow" in x.__name__ or "tf" in x.__name__ for x in backend
+        ):
+            return tensorflow_lambdify()
+
+    return sp.lambdify(symbols, expression, modules=modules, **kwargs)
+
+
+def _use_progress_bar() -> bool:
+    return logging.getLogger().level <= logging.WARNING
diff --git a/src/tensorwaves/interface.py b/src/tensorwaves/interface.py
index 654e94dc..0250f4b3 100644
--- a/src/tensorwaves/interface.py
+++ b/src/tensorwaves/interface.py
@@ -1,16 +1,7 @@
 """Defines top-level interface of tensorwaves."""
 
 from abc import ABC, abstractmethod
-from typing import (
-    Any,
-    Callable,
-    Dict,
-    FrozenSet,
-    Mapping,
-    Optional,
-    Tuple,
-    Union,
-)
+from typing import Any, Dict, Mapping, Optional, Tuple, Union
 
 import attr
 import numpy as np
@@ -77,41 +68,6 @@ def transform(self, dataset: DataSample) -> DataSample:
         """
 
 
-class Model(ABC):
-    """Interface of a model which can be lambdified into a callable."""
-
-    @abstractmethod
-    def lambdify(self, backend: Union[str, tuple, dict]) -> Callable:
-        """Lambdify the model into a Callable.
-
-        Args:
-          backend: Choice of backend for fast evaluations.
-
-        The arguments of the Callable are union of the variables and parameters.
-        The return value of the Callable is Any. In theory the return type
-        should be a value type depending on the model. Currently, there no
-        typing support is implemented for this.
-        """
-
-    @abstractmethod
-    def performance_optimize(self, fix_inputs: DataSample) -> "Model":
-        """Create a performance optimized model, based on fixed inputs."""
-
-    @property
-    @abstractmethod
-    def parameters(self) -> Dict[str, ParameterValue]:
-        """Get mapping of parameters to suggested initial values."""
-
-    @property
-    @abstractmethod
-    def variables(self) -> FrozenSet[str]:
-        """Expected input variable names."""
-
-    @property
-    def argument_order(self) -> Tuple[str, ...]:
-        """Order of arguments of lambdified function signature."""
-
-
 class Estimator(ABC):
     """Estimator for discrepancy model and data."""
 
diff --git a/src/tensorwaves/model/sympy.py b/src/tensorwaves/model/sympy.py
deleted file mode 100644
index 3cd1e57f..00000000
--- a/src/tensorwaves/model/sympy.py
+++ /dev/null
@@ -1,448 +0,0 @@
-"""Lambdify `sympy` expression trees from a `.Model` to a `.Function`."""
-
-import logging
-from copy import deepcopy
-from typing import (
-    Any,
-    Callable,
-    Dict,
-    FrozenSet,
-    List,
-    Optional,
-    Sequence,
-    Tuple,
-    Union,
-)
-
-import sympy as sp
-from sympy.printing.numpy import (
-    NumPyPrinter,
-    _numpy_known_constants,
-    _numpy_known_functions,
-)
-from tqdm.auto import tqdm
-
-from tensorwaves._backend import get_backend_modules
-from tensorwaves.interface import DataSample, Model, ParameterValue
-
-_jax_known_functions = {
-    k: v.replace("numpy.", "jnp.") for k, v in _numpy_known_functions.items()
-}
-_jax_known_constants = {
-    k: v.replace("numpy.", "jnp.") for k, v in _numpy_known_constants.items()
-}
-
-
-class _JaxPrinter(NumPyPrinter):  # pylint: disable=abstract-method
-    # pylint: disable=invalid-name
-    module_imports = {"jax": {"numpy as jnp"}}
-    _module = "jnp"
-    _kc = _jax_known_constants
-    _kf = _jax_known_functions
-
-    def _print_ComplexSqrt(self, expr: sp.Expr) -> str:  # noqa: N802
-        x = self._print(expr.args[0])
-        return (
-            "jnp.select("
-            f"[jnp.less({x}, 0), True], "
-            f"[1j * jnp.sqrt(-{x}), jnp.sqrt({x})], "
-            "default=jnp.nan)"
-        )
-
-
-def split_expression(
-    expression: sp.Expr,
-    max_complexity: int,
-    min_complexity: int = 0,
-) -> Tuple[sp.Expr, Dict[sp.Symbol, sp.Expr]]:
-    """Split an expression into a 'top expression' and several sub-expressions.
-
-    Replace nodes in the expression tree of a `sympy.Expr
-    <sympy.core.expr.Expr>` that lie within a certain complexity range (see
-    :meth:`~sympy.core.basic.Basic.count_ops`) with symbols and keep a mapping
-    of each to these symbols to the sub-expressions that they replaced.
-
-    .. seealso:: :doc:`/usage/faster-lambdify`
-    """
-    i = 0
-    symbol_mapping: Dict[sp.Symbol, sp.Expr] = {}
-    n_operations = sp.count_ops(expression)
-    if n_operations < max_complexity:
-        return expression, symbol_mapping
-    progress_bar = tqdm(
-        total=n_operations,
-        desc="Splitting expression",
-        unit="node",
-        disable=not _use_progress_bar(),
-    )
-
-    def recursive_split(sub_expression: sp.Expr) -> sp.Expr:
-        nonlocal i
-        for arg in sub_expression.args:
-            complexity = sp.count_ops(arg)
-            if min_complexity < complexity < max_complexity:
-                progress_bar.update(n=complexity)
-                symbol = sp.Symbol(f"f{i}")
-                i += 1
-                symbol_mapping[symbol] = arg
-                sub_expression = sub_expression.xreplace({arg: symbol})
-            else:
-                new_arg = recursive_split(arg)
-                sub_expression = sub_expression.xreplace({arg: new_arg})
-        return sub_expression
-
-    top_expression = recursive_split(expression)
-    remainder = progress_bar.total - progress_bar.n
-    progress_bar.update(n=remainder)  # pylint crashes if total is set directly
-    progress_bar.close()
-    return top_expression, symbol_mapping
-
-
-def optimized_lambdify(
-    expression: sp.Expr,
-    symbols: Sequence[sp.Symbol],
-    backend: Union[str, tuple, dict],
-    *,
-    min_complexity: int = 0,
-    max_complexity: int,
-    **kwargs: Any,
-) -> Callable:
-    """Speed up `~sympy.utilities.lambdify.lambdify` with `.split_expression`.
-
-    .. seealso:: :doc:`/usage/faster-lambdify`
-    """
-    top_expression, sub_expressions = split_expression(
-        expression,
-        min_complexity=min_complexity,
-        max_complexity=max_complexity,
-    )
-    if not sub_expressions:
-        return _backend_lambdify(top_expression, symbols, backend, **kwargs)
-
-    sorted_top_symbols = sorted(sub_expressions, key=lambda s: s.name)
-    top_function = _backend_lambdify(
-        top_expression, sorted_top_symbols, backend, **kwargs
-    )
-    sub_functions: List[Callable] = []
-    for symbol in tqdm(
-        iterable=sorted_top_symbols,
-        desc="Lambdifying sub-expressions",
-        unit="expr",
-        disable=not _use_progress_bar(),
-    ):
-        sub_expression = sub_expressions[symbol]
-        sub_function = _backend_lambdify(
-            sub_expression, symbols, backend, **kwargs
-        )
-        sub_functions.append(sub_function)
-
-    def recombined_function(*args: Any) -> Any:
-        new_args = [sub_function(*args) for sub_function in sub_functions]
-        return top_function(*new_args)
-
-    return recombined_function
-
-
-def _sympy_lambdify(
-    expression: sp.Expr,
-    symbols: Sequence[sp.Symbol],
-    backend: Union[str, tuple, dict],
-    *,
-    max_complexity: Optional[int] = None,
-    **kwargs: Any,
-) -> Callable:
-    if max_complexity is None:
-        return _backend_lambdify(
-            expression=expression,
-            symbols=symbols,
-            backend=backend,
-            **kwargs,
-        )
-    return optimized_lambdify(
-        expression=expression,
-        symbols=symbols,
-        backend=backend,
-        max_complexity=max_complexity,
-        **kwargs,
-    )
-
-
-def _backend_lambdify(
-    expression: sp.Expr,
-    symbols: Sequence[sp.Symbol],
-    backend: Union[str, tuple, dict],
-    **kwargs: Any,
-) -> Callable:
-    """A wrapper around :func:`~sympy.utilities.lambdify.lambdify`."""
-    # pylint: disable=import-outside-toplevel,too-many-return-statements
-    def jax_lambdify() -> Callable:
-        import jax
-
-        return jax.jit(
-            sp.lambdify(
-                symbols,
-                expression,
-                modules=modules,
-                printer=_JaxPrinter,
-                **kwargs,
-            )
-        )
-
-    def numba_lambdify() -> Callable:
-        # pylint: disable=import-error
-        import numba
-
-        return numba.jit(
-            sp.lambdify(symbols, expression, modules="numpy", **kwargs),
-            forceobj=True,
-            parallel=True,
-        )
-
-    def tensorflow_lambdify() -> Callable:
-        # pylint: disable=import-error
-        import tensorflow.experimental.numpy as tnp  # pyright: reportMissingImports=false
-
-        return sp.lambdify(symbols, expression, modules=tnp, **kwargs)
-
-    modules = get_backend_modules(backend)
-    if isinstance(backend, str):
-        if backend == "jax":
-            return jax_lambdify()
-        if backend == "numba":
-            return numba_lambdify()
-        if backend in {"tensorflow", "tf"}:
-            return tensorflow_lambdify()
-
-    if isinstance(backend, tuple):
-        if any("jax" in x.__name__ for x in backend):
-            return jax_lambdify()
-        if any("numba" in x.__name__ for x in backend):
-            return numba_lambdify()
-        if any(
-            "tensorflow" in x.__name__ or "tf" in x.__name__ for x in backend
-        ):
-            return tensorflow_lambdify()
-
-    return sp.lambdify(symbols, expression, modules=modules, **kwargs)
-
-
-class _ConstantSubExpressionSympyModel(Model):
-    """Implements a performance optimized sympy based model.
-
-    Based on which symbols of the sympy expression are declared.
-    """
-
-    # pylint: disable=too-many-instance-attributes
-    def __init__(
-        self,
-        expression: sp.Expr,
-        parameters: Dict[sp.Symbol, ParameterValue],
-        fix_inputs: DataSample,
-    ) -> None:
-        self.__fix_inputs = fix_inputs
-        self.__constant_symbols = set(self.__fix_inputs)
-        self.__constant_sub_expressions: Dict[sp.Symbol, sp.Expr] = {}
-        self.__find_constant_subexpressions(expression)
-        self.__expression = self.__replace_constant_sub_expressions(expression)
-        self.__not_fixed_parameters = {
-            k: v
-            for k, v in parameters.items()
-            if k.name not in self.__constant_symbols
-        }
-        self.__not_fixed_variables: FrozenSet[sp.Symbol] = frozenset(
-            s
-            for s in self.__expression.free_symbols
-            if s.name not in self.parameters
-            and s.name not in self.__constant_symbols
-            and s not in self.__constant_sub_expressions
-        )
-        self.__argument_order = tuple(self.__not_fixed_variables) + tuple(
-            self.__not_fixed_parameters
-        )
-
-    def __find_constant_subexpressions(self, expr: sp.Expr) -> bool:
-        if not expr.args:
-            if (
-                isinstance(expr, sp.Symbol)
-                and expr.name not in self.__constant_symbols
-            ):
-                return False
-            return True
-
-        is_constant = True
-        temp_constant_sub_expression = []
-        for arg in expr.args:
-            if self.__find_constant_subexpressions(arg):
-                if arg.args:
-                    temp_constant_sub_expression.append(arg)
-            else:
-                is_constant = False
-
-        if not is_constant:
-            for sub_expr in temp_constant_sub_expression:
-                placeholder = sp.Symbol(f"cached[{str(sub_expr)}]")
-                self.__constant_sub_expressions[placeholder] = sub_expr
-        return is_constant
-
-    def __replace_constant_sub_expressions(
-        self, expression: sp.Expr
-    ) -> sp.Expr:
-        new_expression = deepcopy(expression)
-        return new_expression.xreplace(
-            {v: k for k, v in self.__constant_sub_expressions.items()}
-        )
-
-    def lambdify(self, backend: Union[str, tuple, dict]) -> Callable:
-        input_symbols = tuple(self.__expression.free_symbols)
-        lambdified_model = _backend_lambdify(
-            expression=self.__expression,
-            symbols=input_symbols,
-            backend=backend,
-        )
-        constant_input_storage = {}
-        for placeholder, sub_expr in self.__constant_sub_expressions.items():
-            temp_lambdify = _backend_lambdify(
-                expression=sub_expr,
-                symbols=tuple(sub_expr.free_symbols),
-                backend=backend,
-            )
-            free_symbol_names = {x.name for x in sub_expr.free_symbols}
-            constant_input_storage[placeholder.name] = temp_lambdify(
-                *(self.__fix_inputs[k] for k in free_symbol_names)
-            )
-
-        input_args: list = []
-        non_fixed_arg_positions = list(range(0, len(self.argument_order)))
-
-        for input_arg in input_symbols:
-            if input_arg in self.__argument_order:
-                non_fixed_arg_positions[
-                    self.__argument_order.index(input_arg)
-                ] = len(input_args)
-                input_args.append(0.0)
-            elif input_arg.name in self.__fix_inputs:
-                input_args.append(self.__fix_inputs[input_arg.name])
-            else:
-                input_args.append(constant_input_storage[input_arg.name])
-
-        def update_args(*args: Tuple[Any, ...]) -> None:
-            for i, x in enumerate(args):
-                input_args[non_fixed_arg_positions[i]] = x
-
-        def wrapper(*args: Tuple[Any, ...]) -> Any:
-            update_args(*args)
-            return lambdified_model(*input_args)
-
-        return wrapper
-
-    def performance_optimize(self, fix_inputs: DataSample) -> "Model":
-        return NotImplemented
-
-    @property
-    def parameters(self) -> Dict[str, ParameterValue]:
-        return {
-            symbol.name: value
-            for symbol, value in self.__not_fixed_parameters.items()
-        }
-
-    @property
-    def variables(self) -> FrozenSet[str]:
-        """Expected input variable names."""
-        return frozenset(
-            {symbol.name for symbol in self.__not_fixed_variables}
-        )
-
-    @property
-    def argument_order(self) -> Tuple[str, ...]:
-        return tuple(x.name for x in self.__argument_order)
-
-
-class SympyModel(Model):
-    r"""Full definition of an arbitrary model based on `sympy`.
-
-    Note that input for particle physics amplitude models are based on
-    four-momenta. However, for reasons of convenience, some models may define
-    and use a distinct set of kinematic variables (e.g. in the helicity
-    formalism: angles :math:`\theta` and :math:`\phi`). In this case, a
-    `.DataTransformer` instance (adapter) is needed to transform four momentum
-    information into the custom set of kinematic variables.
-
-    Args:
-        expression: A sympy expression that contains the complete information
-            of the model based on some inputs. The inputs are defined via the
-            `~sympy.core.basic.Basic.free_symbols` attribute of the
-            `sympy.Expr <sympy.core.expr.Expr>`.
-
-        parameters: Defines which inputs of the model are parameters. The keys
-            represent the parameter set, while the values represent their
-            default values. Consequently, the variables of the model are
-            defined as the intersection of the total input set with the
-            parameter set.
-    """
-
-    def __init__(
-        self,
-        expression: sp.Expr,
-        parameters: Dict[sp.Symbol, ParameterValue],
-        max_complexity: Optional[int] = None,
-    ) -> None:
-        if not all(map(lambda p: isinstance(p, sp.Symbol), parameters)):
-            raise TypeError(f"Not all parameters are of type {sp.Symbol}")
-
-        if not set(parameters) <= set(expression.free_symbols):
-            unused_parameters = set(parameters) - set(expression.free_symbols)
-            logging.warning(
-                f"Parameters {unused_parameters} are defined but do not appear"
-                " in the model!"
-            )
-
-        self.__expression = expression
-        # after .doit() certain symbols like the meson radius can disappear
-        # hence the parameters have to be shrunk to this space
-        self.__parameters = {
-            k: v
-            for k, v in parameters.items()
-            if k in self.__expression.free_symbols
-        }
-        self.__variables: FrozenSet[sp.Symbol] = frozenset(
-            s
-            for s in self.__expression.free_symbols
-            if s.name not in self.parameters
-        )
-        self.__argument_order = tuple(self.__variables) + tuple(
-            self.__parameters
-        )
-        self.max_complexity = max_complexity
-
-    def lambdify(self, backend: Union[str, tuple, dict]) -> Callable:
-        """Lambdify the model using `~sympy.utilities.lambdify.lambdify`."""
-        return _sympy_lambdify(
-            expression=self.__expression,
-            symbols=self.__argument_order,
-            backend=backend,
-            max_complexity=self.max_complexity,
-        )
-
-    def performance_optimize(self, fix_inputs: DataSample) -> "Model":
-        return _ConstantSubExpressionSympyModel(
-            self.__expression, self.__parameters, fix_inputs
-        )
-
-    @property
-    def parameters(self) -> Dict[str, ParameterValue]:
-        return {
-            symbol.name: value for symbol, value in self.__parameters.items()
-        }
-
-    @property
-    def variables(self) -> FrozenSet[str]:
-        """Expected input variable names."""
-        return frozenset({symbol.name for symbol in self.__variables})
-
-    @property
-    def argument_order(self) -> Tuple[str, ...]:
-        return tuple(x.name for x in self.__argument_order)
-
-
-def _use_progress_bar() -> bool:
-    return logging.getLogger().level <= logging.WARNING
diff --git a/src/tensorwaves/optimizer/__init__.py b/src/tensorwaves/optimizer/__init__.py
index a5a9d242..f9aba424 100644
--- a/src/tensorwaves/optimizer/__init__.py
+++ b/src/tensorwaves/optimizer/__init__.py
@@ -1,7 +1,7 @@
 """Collection of fit optimizers.
 
-Optimizers 'fit' a `.Model` to some data sample with regard to `.Estimator` as
-loss function.
+Optimizers 'fit' a `.Function` to some data sample with regard to `.Estimator`
+as loss function.
 """
 
 __all__ = [
diff --git a/tox.ini b/tox.ini
index 70c25b80..3c5ecc42 100644
--- a/tox.ini
+++ b/tox.ini
@@ -66,6 +66,8 @@ commands =
         --re-ignore docs/.*\.yml \
         --re-ignore docs/_build/.* \
         --re-ignore docs/api/.* \
+        --re-ignore docs/usage/logs/.* \
+        --re-ignore docs/usage/sub_expr_f? \
         --open-browser \
         docs/ docs/_build/html
 
@@ -124,7 +126,7 @@ allowlist_externals =
     pytest
 commands =
     pytest {posargs} \
-        --cov-fail-under=90 \
+        --cov-fail-under=80 \
         --cov-report=html \
         --cov-report=xml \
         --cov=qrules


Test patch:
diff --git a/tests/integration/ampform/conftest.py b/tests/integration/ampform/conftest.py
index adf25de0..b9fced42 100644
--- a/tests/integration/ampform/conftest.py
+++ b/tests/integration/ampform/conftest.py
@@ -1,4 +1,6 @@
 # pylint: disable=redefined-outer-name
+from typing import Tuple
+
 import ampform
 import pytest
 import qrules
@@ -9,8 +11,9 @@
 from tensorwaves.data import generate_data, generate_phsp
 from tensorwaves.data.phasespace import TFUniformRealNumberGenerator
 from tensorwaves.data.transform import HelicityTransformer
+from tensorwaves.function import LambdifiedFunction
+from tensorwaves.function.sympy import create_function
 from tensorwaves.interface import DataSample, DataTransformer
-from tensorwaves.model import LambdifiedFunction, SympyModel
 
 
 @pytest.fixture(scope="session", params=["canonical", "helicity"])
@@ -43,20 +46,6 @@ def helicity_model(reaction: qrules.ReactionInfo) -> HelicityModel:
     return model_builder.formulate()
 
 
-@pytest.fixture(scope="session", params=["lambdify", "optimized_lambdify"])
-def sympy_model(
-    helicity_model: HelicityModel, request: SubRequest
-) -> SympyModel:
-    max_complexity = None
-    if request.param == "optimized_lambdify":
-        max_complexity = 200
-    return SympyModel(
-        expression=helicity_model.expression.doit(),
-        parameters=helicity_model.parameter_defaults,
-        max_complexity=max_complexity,
-    )
-
-
 @pytest.fixture(scope="session")
 def kinematics(helicity_model: HelicityModel) -> DataTransformer:
     return HelicityTransformer(helicity_model.adapter)
@@ -85,17 +74,43 @@ def phsp_set(
     return kinematics.transform(phsp_sample)
 
 
-@pytest.fixture(scope="session")
-def intensity(sympy_model: SympyModel) -> LambdifiedFunction:
-    return LambdifiedFunction(sympy_model, backend="numpy")
+@pytest.fixture(
+    scope="session",
+    params=[
+        ("jax", False),
+        ("jax", True),
+    ],
+    ids=[
+        "jax-normal",
+        "jax-fast",
+    ],
+)
+def function_fixture(
+    helicity_model: HelicityModel, request: SubRequest
+) -> Tuple[LambdifiedFunction, str]:
+    backend, fast_lambdify = request.param
+    if fast_lambdify:
+        max_complexity = None
+        lambdify_type = f"{backend}-fast"
+    else:
+        max_complexity = 200
+        lambdify_type = f"{backend}-normal"
+    function = create_function(
+        expression=helicity_model.expression.doit(),
+        parameters=helicity_model.parameter_defaults,
+        max_complexity=max_complexity,
+        backend=backend,
+    )
+    return function, lambdify_type
 
 
 @pytest.fixture(scope="session")
 def data_sample(
     reaction: qrules.ReactionInfo,
     kinematics: DataTransformer,
-    intensity: LambdifiedFunction,
+    function_fixture: Tuple[LambdifiedFunction, str],
 ) -> DataSample:
+    function, _ = function_fixture
     n_events = int(1e4)
     initial_state = reaction.initial_state
     final_state = reaction.final_state
@@ -105,7 +120,7 @@ def data_sample(
         initial_state_mass=initial_state[-1].mass,
         final_state_masses={i: p.mass for i, p in final_state.items()},
         data_transformer=kinematics,
-        intensity=intensity,
+        intensity=function,
         random_generator=rng,
     )
     assert all(map(lambda v: len(v) == n_events, sample.values()))
diff --git a/tests/integration/ampform/test_model.py b/tests/integration/ampform/test_model.py
deleted file mode 100644
index c70ac516..00000000
--- a/tests/integration/ampform/test_model.py
+++ /dev/null
@@ -1,40 +0,0 @@
-# pylint: disable=line-too-long, no-self-use, redefined-outer-name
-import qrules
-
-from tensorwaves.model import SympyModel
-
-
-class TestSympyModel:
-    def test_parameters(
-        self, reaction: qrules.ReactionInfo, sympy_model: SympyModel
-    ):
-        if reaction.formalism == "canonical-helicity":
-            assert set(sympy_model.parameters) == {
-                R"C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(500) \gamma;"
-                R" f_{0}(500)"
-                R" \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}}",
-                R"C_{J/\psi(1S) \xrightarrow[S=1]{L=0} f_{0}(980) \gamma;"
-                R" f_{0}(980)"
-                R" \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}}",
-                R"C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(500) \gamma;"
-                R" f_{0}(500)"
-                R" \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}}",
-                R"C_{J/\psi(1S) \xrightarrow[S=1]{L=2} f_{0}(980) \gamma;"
-                R" f_{0}(980)"
-                R" \xrightarrow[S=0]{L=0} \pi^{0} \pi^{0}}",
-                "Gamma_f(0)(500)",
-                "Gamma_f(0)(980)",
-                "m_f(0)(500)",
-                "m_f(0)(980)",
-            }
-        else:
-            assert set(sympy_model.parameters) == {
-                R"C_{J/\psi(1S) \to f_{0}(980)_{0} \gamma_{+1}; f_{0}(980) \to"
-                R" \pi^{0}_{0} \pi^{0}_{0}}",
-                R"C_{J/\psi(1S) \to f_{0}(500)_{0} \gamma_{+1}; f_{0}(500) \to"
-                R" \pi^{0}_{0} \pi^{0}_{0}}",
-                "m_f(0)(980)",
-                "Gamma_f(0)(980)",
-                "m_f(0)(500)",
-                "Gamma_f(0)(500)",
-            }
diff --git a/tests/integration/ampform/test_optimize.py b/tests/integration/ampform/test_optimize.py
index 4738353d..54dc06b5 100644
--- a/tests/integration/ampform/test_optimize.py
+++ b/tests/integration/ampform/test_optimize.py
@@ -1,12 +1,13 @@
 # pylint: disable=redefined-outer-name
 from pathlib import Path
+from typing import Tuple
 
 import pytest
 import qrules
 
 from tensorwaves.estimator import UnbinnedNLL
+from tensorwaves.function import LambdifiedFunction
 from tensorwaves.interface import DataSample
-from tensorwaves.model.sympy import SympyModel
 from tensorwaves.optimizer.callbacks import (
     CallbackList,
     CSVSummary,
@@ -17,10 +18,13 @@
 
 @pytest.fixture(scope="session")
 def estimator(
-    sympy_model: SympyModel, data_set: DataSample, phsp_set: DataSample
+    function_fixture: Tuple[LambdifiedFunction, str],
+    data_set: DataSample,
+    phsp_set: DataSample,
 ) -> UnbinnedNLL:
+    function, _ = function_fixture
     return UnbinnedNLL(
-        sympy_model,
+        function,
         dict(data_set),
         dict(phsp_set),
         backend="jax",
@@ -29,15 +33,12 @@ def estimator(
 
 def test_fit_and_callbacks(  # pylint: disable=too-many-locals
     estimator: UnbinnedNLL,
+    function_fixture: Tuple[LambdifiedFunction, str],
     output_dir: Path,
     reaction: qrules.ReactionInfo,
-    sympy_model: SympyModel,
 ):
+    _, lambdify_type = function_fixture
     formalism_alias = reaction.formalism[:4]
-    if sympy_model.max_complexity is None:
-        lambdify_type = "normal"
-    else:
-        lambdify_type = "optimized"
     filename = output_dir / f"fit_result_{formalism_alias}_{lambdify_type}"
     csv_file = f"{filename}.csv"
     yml_file = f"{filename}.yml"
diff --git a/tests/unit/model/__init__.py b/tests/unit/function/__init__.py
similarity index 100%
rename from tests/unit/model/__init__.py
rename to tests/unit/function/__init__.py
diff --git a/tests/unit/function/test_function.py b/tests/unit/function/test_function.py
new file mode 100644
index 00000000..fb2f1ed4
--- /dev/null
+++ b/tests/unit/function/test_function.py
@@ -0,0 +1,50 @@
+# pylint: disable=no-self-use, redefined-outer-name
+import numpy as np
+import pytest
+import sympy as sp
+
+from tensorwaves.function import LambdifiedFunction
+from tensorwaves.function.sympy import create_function
+from tensorwaves.interface import DataSample, Function
+
+
+class TestLambdifiedFunction:
+    @pytest.fixture(scope="module")
+    def function(self) -> LambdifiedFunction:
+        c_1, c_2, c_3, c_4 = sp.symbols("c_(1:5)")
+        x = sp.Symbol("x", complex_twice=True)
+        parameters = {
+            c_1: 1 + 1j,
+            c_2: -1 + 1j,
+            c_3: 1 - 1j,
+            c_4: -1 - 1j,
+        }
+        expression = (
+            c_1 * sp.sqrt(x) / x
+            + c_2
+            * sp.exp(-sp.Rational(1, 2) * ((x - 2) / sp.Rational(1, 2)) ** 2)
+            + c_3 * (x ** 2 - 3 * x)
+            + c_4
+        )
+        expression = sp.simplify(sp.conjugate(expression) * expression)
+        return create_function(expression, parameters, backend="numpy")
+
+    @pytest.mark.parametrize(
+        ("test_data", "expected_results"),
+        [
+            (
+                {"x": np.array([0.5, 1.0, 1.5, 2.0, 2.5])},
+                [3.52394, 9.11931, 16.3869, 18.1716, 7.16359],
+            ),
+        ],
+    )
+    def test_call(
+        self,
+        function: Function,
+        test_data: DataSample,
+        expected_results: np.ndarray,
+    ):
+        results = function(test_data)
+        np.testing.assert_array_almost_equal(
+            results, expected_results, decimal=4
+        )
diff --git a/tests/unit/function/test_sympy.py b/tests/unit/function/test_sympy.py
new file mode 100644
index 00000000..07109517
--- /dev/null
+++ b/tests/unit/function/test_sympy.py
@@ -0,0 +1,60 @@
+# cspell:ignore lambdifygenerated
+import sys
+
+import numpy as np
+import pytest
+import sympy as sp
+
+from tensorwaves.function.sympy import fast_lambdify, split_expression
+
+
+def create_expression(x, y, z):
+    return x ** z + 2 * y
+
+
+@pytest.mark.parametrize("backend", ["jax", "math", "numpy", "tf"])
+@pytest.mark.parametrize("max_complexity", [2, 3, 4])
+def test_fast_lambdify(backend: str, max_complexity: int):
+    x, y, z = sp.symbols("x y z")
+    expression = create_expression(x, y, z)
+    function = fast_lambdify(
+        expression,
+        symbols=[x, y, z],
+        max_complexity=max_complexity,
+        backend=backend,
+    )
+
+    func_repr = str(function)
+    if max_complexity <= 3:
+        assert func_repr.startswith("<function fast_lambdify.<locals>")
+    else:
+        repr_start = "<function _lambdifygenerated"
+        if backend == "jax":
+            if sys.version_info >= (3, 7):
+                repr_start = "<CompiledFunction of " + repr_start
+            else:
+                repr_start = "<CompiledFunction object at 0x"
+        assert func_repr.startswith(repr_start)
+
+    data = (
+        np.array([1, 2]),
+        np.array([1, np.e]),
+        np.array([1, 2]),
+    )
+    output = function(*data)
+    expected = create_expression(*data)
+    assert pytest.approx(output) == expected
+
+
+def test_split_expression():
+    x, y, z = sp.symbols("x y z")
+    expression = create_expression(x, y, z)
+    top_expr, sub_expressions = split_expression(expression, max_complexity=3)
+    assert top_expr.free_symbols == set(sub_expressions)
+    assert expression == top_expr.xreplace(sub_expressions)
+
+    sub_symbols = sorted(top_expr.free_symbols, key=lambda s: s.name)
+    assert len(sub_symbols) == 2
+    f1, f2 = tuple(sub_symbols)
+    assert sub_expressions[f1] == x ** z
+    assert sub_expressions[f2] == 2 * y
diff --git a/tests/unit/model/test_sympy.py b/tests/unit/model/test_sympy.py
deleted file mode 100644
index aff2280f..00000000
--- a/tests/unit/model/test_sympy.py
+++ /dev/null
@@ -1,106 +0,0 @@
-# cspell:ignore lambdifygenerated
-# pylint: disable=invalid-name, no-self-use, redefined-outer-name
-import sys
-
-import numpy as np
-import pytest
-import sympy as sp
-
-from tensorwaves.interface import DataSample, Function
-from tensorwaves.model import LambdifiedFunction, SympyModel
-from tensorwaves.model.sympy import optimized_lambdify, split_expression
-
-
-class TestLambdifiedFunction:
-    @pytest.fixture(scope="module")
-    def function(self) -> LambdifiedFunction:
-        c_1, c_2, c_3, c_4 = sp.symbols("c_(1:5)")
-        x = sp.Symbol("x", complex_twice=True)
-        parameters = {
-            c_1: 1 + 1j,
-            c_2: -1 + 1j,
-            c_3: 1 - 1j,
-            c_4: -1 - 1j,
-        }
-        expression = (
-            c_1 * sp.sqrt(x) / x
-            + c_2
-            * sp.exp(-sp.Rational(1, 2) * ((x - 2) / sp.Rational(1, 2)) ** 2)
-            + c_3 * (x ** 2 - 3 * x)
-            + c_4
-        )
-        expression = sp.simplify(sp.conjugate(expression) * expression)
-        model = SympyModel(expression=expression, parameters=parameters)
-        return LambdifiedFunction(model, "numpy")
-
-    @pytest.mark.parametrize(
-        ("test_data", "expected_results"),
-        [
-            (
-                {"x": np.array([0.5, 1.0, 1.5, 2.0, 2.5])},
-                [3.52394, 9.11931, 16.3869, 18.1716, 7.16359],
-            ),
-        ],
-    )
-    def test_call(
-        self,
-        function: Function,
-        test_data: DataSample,
-        expected_results: np.ndarray,
-    ):
-        results = function(test_data)
-        np.testing.assert_array_almost_equal(
-            results, expected_results, decimal=4
-        )
-
-
-def create_expression(x, y, z):
-    return x ** z + 2 * y
-
-
-@pytest.mark.parametrize("backend", ["jax", "math", "numpy", "tf"])
-@pytest.mark.parametrize("max_complexity", [2, 3, 4])
-def test_optimized_lambdify(backend: str, max_complexity: int):
-    x, y, z = sp.symbols("x y z")
-    expression = create_expression(x, y, z)
-    function = optimized_lambdify(
-        expression,
-        symbols=[x, y, z],
-        max_complexity=max_complexity,
-        backend=backend,
-    )
-
-    func_repr = str(function)
-    if max_complexity <= 3:
-        assert func_repr.startswith("<function optimized_lambdify.<locals>")
-    else:
-        repr_start = "<function _lambdifygenerated"
-        if backend == "jax":
-            if sys.version_info >= (3, 7):
-                repr_start = "<CompiledFunction of " + repr_start
-            else:
-                repr_start = "<CompiledFunction object at 0x"
-        assert func_repr.startswith(repr_start)
-
-    data = (
-        np.array([1, 2]),
-        np.array([1, np.e]),
-        np.array([1, 2]),
-    )
-    output = function(*data)
-    expected = create_expression(*data)
-    assert pytest.approx(output) == expected
-
-
-def test_split_expression():
-    x, y, z = sp.symbols("x y z")
-    expression = create_expression(x, y, z)
-    top_expr, sub_expressions = split_expression(expression, max_complexity=3)
-    assert top_expr.free_symbols == set(sub_expressions)
-    assert expression == top_expr.xreplace(sub_expressions)
-
-    sub_symbols = sorted(top_expr.free_symbols, key=lambda s: s.name)
-    assert len(sub_symbols) == 2
-    f1, f2 = tuple(sub_symbols)
-    assert sub_expressions[f1] == x ** z
-    assert sub_expressions[f2] == 2 * y
diff --git a/tests/unit/test_estimator.py b/tests/unit/test_estimator.py
index 00621a84..6944f66f 100644
--- a/tests/unit/test_estimator.py
+++ b/tests/unit/test_estimator.py
@@ -7,19 +7,21 @@
 import sympy as sp
 
 from tensorwaves.estimator import UnbinnedNLL
+from tensorwaves.function import LambdifiedFunction
+from tensorwaves.function.sympy import create_function
 from tensorwaves.interface import DataSample, ParameterValue
-from tensorwaves.model import SympyModel
 from tensorwaves.optimizer.minuit import Minuit2
 
 
-def gaussian(mu_: float, sigma_: float) -> SympyModel:
+def gaussian(mu_: float, sigma_: float) -> LambdifiedFunction:
     x, mu, sigma = sp.symbols("x, mu, sigma")
-    return SympyModel(
-        expression=(sp.exp(-(((x - mu) / sigma) ** 2) / 2)),
+    return create_function(
+        expression=sp.exp(-(((x - mu) / sigma) ** 2) / 2),
         parameters={
             mu: mu_,
             sigma: sigma_,
         },
+        backend="numpy",
     )
 
 
@@ -30,7 +32,7 @@ def gaussian_sum(
     a_2: float,
     mu_2: float,
     sigma_2: float,
-) -> SympyModel:
+) -> LambdifiedFunction:
     x, a1, mu1, sigma1, a2, mu2, sigma2 = sp.symbols(
         "x, a1, mu1, sigma1, a2, mu2, sigma2"
     )
@@ -45,7 +47,7 @@ def gaussian_sum(
         * sp.exp(-(((x - mu2) / sigma2) ** 2) / 2)
     )
 
-    return SympyModel(
+    return create_function(
         expression=gaussian1 + gaussian2,
         parameters={
             a1: a_1,
@@ -55,6 +57,7 @@ def gaussian_sum(
             mu2: mu_2,
             sigma2: sigma_2,
         },
+        backend="numpy",
     )
 
 
@@ -66,23 +69,23 @@ def phsp_dataset() -> DataSample:
     }
 
 
-__np_rng = np.random.default_rng(12345)
+NUMPY_RNG = np.random.default_rng(12345)
 
 
 @pytest.mark.parametrize(
-    ("model", "dataset", "true_params"),
+    ("function", "dataset", "true_params"),
     [
         (
             gaussian(1.0, 0.1),
             {
-                "x": __np_rng.normal(0.5, 0.1, 1000),
+                "x": NUMPY_RNG.normal(0.5, 0.1, 1000),
             },
             {"mu": 0.5},
         ),
         (
             gaussian(1.0, 0.1),
             {
-                "x": __np_rng.normal(0.5, 0.3, 1000),
+                "x": NUMPY_RNG.normal(0.5, 0.3, 1000),
             },
             {"mu": 0.5, "sigma": 0.3},
         ),
@@ -90,12 +93,12 @@ def phsp_dataset() -> DataSample:
             gaussian_sum(1.0, 1.0, 0.1, 2.0, 2.0, 0.3),
             {
                 "x": np.append(
-                    __np_rng.normal(
+                    NUMPY_RNG.normal(
                         1.0,
                         0.1,
                         2000,
                     ),
-                    __np_rng.normal(
+                    NUMPY_RNG.normal(
                         2.0,
                         0.3,
                         1000,
@@ -110,12 +113,12 @@ def phsp_dataset() -> DataSample:
             gaussian_sum(1.0, 1.0, 0.1, 1.0, 2.0, 0.3),
             {
                 "x": np.append(
-                    __np_rng.normal(
+                    NUMPY_RNG.normal(
                         0.9,
                         0.3,
                         1000,
                     ),
-                    __np_rng.normal(
+                    NUMPY_RNG.normal(
                         2.5,
                         0.1,
                         1000,
@@ -128,12 +131,12 @@ def phsp_dataset() -> DataSample:
             gaussian_sum(1.0, 1.0, 0.1, 2.0, 2.5, 0.3),
             {
                 "x": np.append(
-                    __np_rng.normal(
+                    NUMPY_RNG.normal(
                         0.9,
                         0.3,
                         2000,
                     ),
-                    __np_rng.normal(
+                    NUMPY_RNG.normal(
                         2.5,
                         0.1,
                         1000,
@@ -145,13 +148,13 @@ def phsp_dataset() -> DataSample:
     ],
 )
 def test_sympy_unbinned_nll(
-    model: SympyModel,
+    function: LambdifiedFunction,
     dataset: DataSample,
     true_params: Dict[str, ParameterValue],
     phsp_dataset: DataSample,
 ):
     estimator = UnbinnedNLL(
-        model,
+        function,
         dataset,
         phsp_dataset,
         phsp_volume=6.0,