diff --git a/.github/workflows/codecov.yml b/.github/workflows/codecov.yml
index a3cfffaa..eb288b2f 100644
--- a/.github/workflows/codecov.yml
+++ b/.github/workflows/codecov.yml
@@ -15,10 +15,10 @@ jobs:
 
       - uses: actions/checkout@v4
 
-      - name: Set up Python 3.7
+      - name: Set up Python 3.10
         uses: actions/setup-python@v5
         with:
-          python-version: 3.7
+          python-version: "3.10"
 
       - name: Install dependencies
         run: |
diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml
index de5a03cc..c26d80ec 100644
--- a/.github/workflows/format.yml
+++ b/.github/workflows/format.yml
@@ -11,10 +11,10 @@ jobs:
     runs-on: ubuntu-latest
     steps:
     - uses: actions/checkout@v4
-    - name: Set up Python 3.7
+    - name: Set up Python 3.10
       uses: actions/setup-python@v5
       with:
-        python-version: 3.7
+        python-version: "3.10"
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 86230f1a..20d77d37 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8, 3.9]
+        python-version: ["3.9", "3.10", "3.11"]
 
     steps:
     - uses: actions/checkout@v4
@@ -40,15 +40,15 @@ jobs:
 
     strategy:
       matrix:
-        os: [ubuntu-latest, macos-latest, windows-latest]
+        os: [ubuntu-latest, macos-latest]
     runs-on: ${{ matrix.os }}
 
     steps:
     - uses: actions/checkout@v4
-    - name: Set up Python 3.7
+    - name: Set up Python 3.10
       uses: actions/setup-python@v5
       with:
-        python-version: 3.7
+        python-version: "3.10"
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
diff --git a/emukit/sensitivity/monte_carlo/monte_carlo_sensitivity.py b/emukit/sensitivity/monte_carlo/monte_carlo_sensitivity.py
index 2c0ae664..31a79be2 100644
--- a/emukit/sensitivity/monte_carlo/monte_carlo_sensitivity.py
+++ b/emukit/sensitivity/monte_carlo/monte_carlo_sensitivity.py
@@ -51,9 +51,7 @@ def saltelli_estimators(
         Saltelli estimators of the total mean and variance
         """
 
-        variable_main_variance = (
-            sum(f_main_sample * f_new_fixing_sample) / (num_monte_carlo_points - 1) - total_mean**2
-        )
+        variable_main_variance = sum(f_main_sample * f_new_fixing_sample) / (num_monte_carlo_points - 1) - total_mean**2
         variable_total_variance = (
             total_variance - sum(f_fixing_sample * f_new_fixing_sample) / (num_monte_carlo_points - 1) + total_mean**2
         )
diff --git a/integration_tests/emukit/quadrature/test_wsabil_loop.py b/integration_tests/emukit/quadrature/test_wsabil_loop.py
index 9f2772e3..f8f4a9f5 100644
--- a/integration_tests/emukit/quadrature/test_wsabil_loop.py
+++ b/integration_tests/emukit/quadrature/test_wsabil_loop.py
@@ -6,7 +6,6 @@
 import numpy as np
 import pytest
 from numpy.testing import assert_array_equal
-from pytest_lazyfixture import lazy_fixture
 
 from emukit.core.loop.user_function import UserFunctionWrapper
 from emukit.model_wrappers.gpy_quadrature_wrappers import BaseGaussianProcessGPy, QuadratureRBFGaussianMeasure, RBFGPy
@@ -58,12 +57,12 @@ def loop_fixed(wsabil_fixed):
     return emukit_loop, Y.shape[0], X, Y
 
 
-wsabi_test_list = [lazy_fixture("loop_adapt"), lazy_fixture("loop_fixed")]
+wsabi_test_list = ["loop_adapt", "loop_fixed"]
 
 
 @pytest.mark.parametrize("loop", wsabi_test_list)
-def test_wsabil_loop(loop):
-    emukit_loop, init_size, _, _ = loop
+def test_wsabil_loop(loop, request):
+    emukit_loop, init_size, _, _ = request.getfixturevalue(loop)
     num_iter = 5
 
     emukit_loop.run_loop(user_function=UserFunctionWrapper(func), stopping_condition=num_iter)
@@ -73,8 +72,8 @@ def test_wsabil_loop(loop):
 
 
 @pytest.mark.parametrize("loop", wsabi_test_list)
-def test_wsabil_loop_initial_state(loop):
-    emukit_loop, _, x_init, y_init = loop
+def test_wsabil_loop_initial_state(loop, request):
+    emukit_loop, _, x_init, y_init = request.getfixturevalue(loop)
 
     assert_array_equal(emukit_loop.loop_state.X, x_init)
     assert_array_equal(emukit_loop.loop_state.Y, y_init)
diff --git a/requirements/requirements.txt b/requirements/requirements.txt
index a4fdda6c..c1d11415 100644
--- a/requirements/requirements.txt
+++ b/requirements/requirements.txt
@@ -1,10 +1,10 @@
 setuptools>=40.8.0
 # numpy deprecated np.bool in 1.24, and it is still used in GPy
 # until GPy is updated we need to avoid 1.24 or higher
-numpy>=1.14.5,<1.24
+numpy>=1.23
 # This is unfortunate - we don't need matplotlib
 # but until GPy and GPyOpt get their dependencies straight
 # we need GPy's plotting extra to ensure smooth installation
-GPy[plotting]>=1.10.0
+GPy[plotting]>=1.13.0
 emcee>=2.2.1
 scipy>=1.1.0
diff --git a/requirements/test_requirements.txt b/requirements/test_requirements.txt
index 24746d33..3d340200 100644
--- a/requirements/test_requirements.txt
+++ b/requirements/test_requirements.txt
@@ -6,7 +6,6 @@ black
 pytest>=3.5.1
 pytest-cov>=2.5.1
 mock>=2.0.0
-pytest-lazy-fixture>=0.4.2
 
 # For Latin design
 PyDOE>=0.3.0
diff --git a/setup.py b/setup.py
index b1ef58d2..04603ae5 100644
--- a/setup.py
+++ b/setup.py
@@ -31,7 +31,7 @@
     include_package_data=True,
     install_requires=requires,
     extras_require={"benchmarking": ["matplotlib"]},
-    python_requires=">=3",
+    python_requires=">=3.9",
     license="Apache License 2.0",
     classifiers=(
         # https://pypi.org/pypi?%3Aaction=list_classifiers
@@ -40,9 +40,6 @@
         "Intended Audience :: Education",
         "Intended Audience :: Science/Research",
         "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.4",
-        "Programming Language :: Python :: 3.5",
-        "Programming Language :: Python :: 3.6",
         "License :: OSI Approved :: Apache Software License",
         "Operating System :: OS Independent",
         "Topic :: Scientific/Engineering :: Artificial Intelligence",
diff --git a/tests/emukit/model_wrappers/test_gpy_wrappers_quadrature.py b/tests/emukit/model_wrappers/test_gpy_wrappers_quadrature.py
index 4d6381fa..558188b4 100644
--- a/tests/emukit/model_wrappers/test_gpy_wrappers_quadrature.py
+++ b/tests/emukit/model_wrappers/test_gpy_wrappers_quadrature.py
@@ -7,7 +7,6 @@
 import GPy
 import numpy as np
 import pytest
-from pytest_lazyfixture import lazy_fixture
 
 from emukit.model_wrappers.gpy_quadrature_wrappers import (
     BaseGaussianProcessGPy,
@@ -204,21 +203,22 @@ def wrapper_matern52_2(dim1, gpy_matern52):
 
 
 gpy_test_list = [
-    lazy_fixture("wrapper_rbf_1"),
-    lazy_fixture("wrapper_rbf_2"),
-    lazy_fixture("wrapper_brownian_1"),
-    lazy_fixture("wrapper_brownian_2"),
-    lazy_fixture("wrapper_matern12_1"),
-    lazy_fixture("wrapper_matern12_2"),
-    lazy_fixture("wrapper_matern32_1"),
-    lazy_fixture("wrapper_matern32_2"),
-    lazy_fixture("wrapper_matern52_1"),
-    lazy_fixture("wrapper_matern52_2"),
+    "wrapper_rbf_1",
+    "wrapper_rbf_2",
+    "wrapper_brownian_1",
+    "wrapper_brownian_2",
+    "wrapper_matern12_1",
+    "wrapper_matern12_2",
+    "wrapper_matern32_1",
+    "wrapper_matern32_2",
+    "wrapper_matern52_1",
+    "wrapper_matern52_2",
 ]
 
 
-@pytest.mark.parametrize("wrapper", gpy_test_list)
-def test_create_emukit_model_from_gpy_model_types(wrapper):
+@pytest.mark.parametrize("wrapper_name", gpy_test_list)
+def test_create_emukit_model_from_gpy_model_types(wrapper_name, request):
+    wrapper = request.getfixturevalue(wrapper_name)
     gpy_model = GPy.models.GPRegression(kernel=wrapper["gpy_kernel"], X=wrapper["data"][0], Y=wrapper["data"][1])
     emukit_gp = create_emukit_model_from_gpy_model(gpy_model=gpy_model, measure=wrapper["measure"])
 
@@ -252,7 +252,7 @@ def test_create_emukit_model_from_gpy_model_raises_warns():
         create_emukit_model_from_gpy_model(gpy_model=gpy_model, integral_bounds=bounds, measure=measure)
 
 
-def test_base_gp_gpy_raises(gpy_prodbrownian):
+def test_base_gp_gpy_raises():
     incompatible_offset = -3
 
     n_dim = 2
diff --git a/tests/emukit/multi_fidelity/test_kernels.py b/tests/emukit/multi_fidelity/test_kernels.py
index c0776519..4fb29ffc 100644
--- a/tests/emukit/multi_fidelity/test_kernels.py
+++ b/tests/emukit/multi_fidelity/test_kernels.py
@@ -10,7 +10,7 @@
 """
 import GPy
 import numpy as np
-from GPy.testing.kernel_tests import check_kernel_gradient_functions
+from GPy.testing.test_kernel import check_kernel_gradient_functions
 
 import emukit.multi_fidelity
 
diff --git a/tests/emukit/quadrature/test_measures.py b/tests/emukit/quadrature/test_measures.py
index 9b16205e..b4b7dff0 100644
--- a/tests/emukit/quadrature/test_measures.py
+++ b/tests/emukit/quadrature/test_measures.py
@@ -9,7 +9,6 @@
 
 import numpy as np
 import pytest
-from pytest_lazyfixture import lazy_fixture
 from utils import check_grad
 
 from emukit.quadrature.measures import BoxDomain, GaussianMeasure, LebesgueMeasure
@@ -77,10 +76,10 @@ def gauss_measure():
 
 
 measure_test_list = [
-    lazy_fixture("lebesgue_measure"),
-    lazy_fixture("lebesgue_measure_normalized"),
-    lazy_fixture("gauss_iso_measure"),
-    lazy_fixture("gauss_measure"),
+    DataLebesgueMeasure(),
+    DataLebesgueNormalizedMeasure(),
+    DataGaussIsoMeasure(),
+    DataGaussMeasure(),
 ]
 
 
diff --git a/tests/emukit/quadrature/test_quadrature_acquisitions.py b/tests/emukit/quadrature/test_quadrature_acquisitions.py
index ffc40130..6c5a4425 100644
--- a/tests/emukit/quadrature/test_quadrature_acquisitions.py
+++ b/tests/emukit/quadrature/test_quadrature_acquisitions.py
@@ -8,7 +8,6 @@
 import GPy
 import numpy as np
 import pytest
-from pytest_lazyfixture import lazy_fixture
 from utils import check_grad
 
 from emukit.model_wrappers.gpy_quadrature_wrappers import BaseGaussianProcessGPy, RBFGPy
@@ -58,50 +57,49 @@ def model_gaussian(gpy_model):
 
 
 model_test_list = [
-    lazy_fixture("model_gaussian"),
-    lazy_fixture("model_lebesgue"),
-    lazy_fixture("model_lebesgue_normalized"),
+    "model_gaussian",
+    "model_lebesgue",
+    "model_lebesgue_normalized",
 ]
 
-
-@pytest.fixture(params=model_test_list)
-def model_test_list_fixture(request):
-    return request.param
-
-
 # === acquisition fixtures start here
 
 
 @pytest.fixture
-def mutual_information(model_test_list_fixture):
-    return MutualInformation(model_test_list_fixture)
+def mutual_information():
+    return lambda model: MutualInformation(model)
 
 
 @pytest.fixture
-def squared_correlation(model_test_list_fixture):
-    return SquaredCorrelation(model_test_list_fixture)
+def squared_correlation():
+    return lambda model: SquaredCorrelation(model)
 
 
 @pytest.fixture
-def integral_variance_reduction(model_test_list_fixture):
-    return IntegralVarianceReduction(model_test_list_fixture)
+def integral_variance_reduction():
+    return lambda model: IntegralVarianceReduction(model)
 
 
 @pytest.fixture
-def uncertainty_sampling(model_test_list_fixture):
-    return UncertaintySampling(model_test_list_fixture)
+def uncertainty_sampling():
+    return lambda model: UncertaintySampling(model)
 
 
-acquisitions_test_list = [
-    lazy_fixture("mutual_information"),
-    lazy_fixture("squared_correlation"),
-    lazy_fixture("integral_variance_reduction"),
-    lazy_fixture("uncertainty_sampling"),
+acquisition_test_list = [
+    "mutual_information",
+    "squared_correlation",
+    "integral_variance_reduction",
+    "uncertainty_sampling",
 ]
 
 
-@pytest.mark.parametrize("aq", acquisitions_test_list)
-def test_quadrature_acquisition_shapes(aq):
+@pytest.mark.parametrize("model_name", model_test_list)
+@pytest.mark.parametrize("aq_name", acquisition_test_list)
+def test_quadrature_acquisition_shapes(model_name, aq_name, request):
+    model = request.getfixturevalue(model_name)
+    aq_factory = request.getfixturevalue(aq_name)
+    aq = aq_factory(model)
+
     x = np.array([[-1, 1], [0, 0], [-2, 0.1]])
 
     # value
@@ -114,8 +112,13 @@ def test_quadrature_acquisition_shapes(aq):
     assert res[1].shape == (3, 2)
 
 
-@pytest.mark.parametrize("aq", acquisitions_test_list)
-def test_quadrature_acquisition_gradient_values(aq):
+@pytest.mark.parametrize("model_name", model_test_list)
+@pytest.mark.parametrize("aq_name", acquisition_test_list)
+def test_quadrature_acquisition_gradient_values(model_name, aq_name, request):
+    model = request.getfixturevalue(model_name)
+    aq_factory = request.getfixturevalue(aq_name)
+    aq = aq_factory(model)
+
     func = lambda x: aq.evaluate(x)[:, 0]
     dfunc = lambda x: aq.evaluate_with_gradients(x)[1].T
     check_grad(func, dfunc, in_shape=(3, 2), bounds=aq.model.X.shape[1] * [(-3, 3)])
diff --git a/tests/emukit/quadrature/test_quadrature_kernels.py b/tests/emukit/quadrature/test_quadrature_kernels.py
index a486524e..f92505ec 100644
--- a/tests/emukit/quadrature/test_quadrature_kernels.py
+++ b/tests/emukit/quadrature/test_quadrature_kernels.py
@@ -7,7 +7,6 @@
 import GPy
 import numpy as np
 import pytest
-from pytest_lazyfixture import lazy_fixture
 from utils import check_grad, sample_uniform
 
 from emukit.model_wrappers.gpy_quadrature_wrappers import (
@@ -297,25 +296,25 @@ def lebesgue_normalized_qprodbrownian():
 
 
 gaussian_embeddings_test_dict = {
-    "qrbf": lazy_fixture("gaussian_qrbf"),
+    "qrbf": "gaussian_qrbf",
 }
 
 lebesgue_embeddings_test_dict = {
-    "qrbf": lazy_fixture("lebesgue_qrbf"),
-    "qmatern12": lazy_fixture("lebesgue_qmatern12"),
-    "qmatern32": lazy_fixture("lebesgue_qmatern32"),
-    "qmatern52": lazy_fixture("lebesgue_qmatern52"),
-    "qbrownian": lazy_fixture("lebesgue_qbrownian"),
-    "qprodbrownian": lazy_fixture("lebesgue_qprodbrownian"),
+    "qrbf": "lebesgue_qrbf",
+    "qmatern12": "lebesgue_qmatern12",
+    "qmatern32": "lebesgue_qmatern32",
+    "qmatern52": "lebesgue_qmatern52",
+    "qbrownian": "lebesgue_qbrownian",
+    "qprodbrownian": "lebesgue_qprodbrownian",
 }
 
 lebesgue_normalized_embeddings_test_dict = {
-    "qrbf": lazy_fixture("lebesgue_normalized_qrbf"),
-    "qmatern12": lazy_fixture("lebesgue_normalized_qmatern12"),
-    "qmatern32": lazy_fixture("lebesgue_normalized_qmatern32"),
-    "qmatern52": lazy_fixture("lebesgue_normalized_qmatern52"),
-    "qbrownian": lazy_fixture("lebesgue_normalized_qbrownian"),
-    "qprodbrownian": lazy_fixture("lebesgue_normalized_qprodbrownian"),
+    "qrbf": "lebesgue_normalized_qrbf",
+    "qmatern12": "lebesgue_normalized_qmatern12",
+    "qmatern32": "lebesgue_normalized_qmatern32",
+    "qmatern52": "lebesgue_normalized_qmatern52",
+    "qbrownian": "lebesgue_normalized_qbrownian",
+    "qprodbrownian": "lebesgue_normalized_qprodbrownian",
 }
 
 embeddings_test_list = (
@@ -326,8 +325,8 @@ def lebesgue_normalized_qprodbrownian():
 
 
 @pytest.mark.parametrize("kernel_embedding", embeddings_test_list)
-def test_qkernel_shapes(kernel_embedding):
-    emukit_qkernel, x1, x2, N, M, D, _ = kernel_embedding
+def test_qkernel_shapes(kernel_embedding, request):
+    emukit_qkernel, x1, x2, N, M, D, _ = request.getfixturevalue(kernel_embedding)
 
     # kernel shapes
     assert emukit_qkernel.K(x1, x2).shape == (N, M)
@@ -360,13 +359,13 @@ def test_qkernel_shapes(kernel_embedding):
         (lebesgue_normalized_embeddings_test_dict["qprodbrownian"], [5.199166451419295, 5.204923979414083]),
     ],
 )
-def test_qkernel_qKq(kernel_embedding, interval):
+def test_qkernel_qKq(kernel_embedding, interval, request):
     # To test the integral value of the kernel embedding, we check if it lies in some confidence interval.
     # These intervals were computed as follows: The kernel emukit_qkernel.qK was integrated by
     # simple random sampling with 1e6 samples. This was done 100 times. The intervals show mean\pm 3 std of the 100
     # integrals obtained by sampling. There might be a small chance that the true integrals lies outside the
     # specified intervals.
-    emukit_qkernel = kernel_embedding[0]
+    emukit_qkernel = request.getfixturevalue(kernel_embedding)[0]
     qKq = emukit_qkernel.qKq()
     assert interval[0] < qKq < interval[1]
 
@@ -519,9 +518,9 @@ def test_qkernel_qKq(kernel_embedding, interval):
         ),
     ],
 )
-def test_qkernel_qK(kernel_embedding, intervals):
+def test_qkernel_qK(kernel_embedding, intervals, request):
     # See test_qkernel_qKq on how the intervals were computed.
-    emukit_qkernel, _, x2, _, _, _, _ = kernel_embedding
+    emukit_qkernel, _, x2, _, _, _, _ = request.getfixturevalue(kernel_embedding)
     qK = emukit_qkernel.qK(x2)[0, :]
     for i in range(4):
         assert intervals[i, 0] < qK[i] < intervals[i, 1]
@@ -531,8 +530,8 @@ def test_qkernel_qK(kernel_embedding, intervals):
 
 
 @pytest.mark.parametrize("kernel_embedding", embeddings_test_list)
-def test_qkernel_gradient_shapes(kernel_embedding):
-    emukit_qkernel, x1, x2, N, M, D, _ = kernel_embedding
+def test_qkernel_gradient_shapes(kernel_embedding, request):
+    emukit_qkernel, x1, x2, N, M, D, _ = request.getfixturevalue(kernel_embedding)
 
     # gradient of kernel
     assert emukit_qkernel.dK_dx1(x1, x2).shape == (D, N, M)
@@ -545,8 +544,8 @@ def test_qkernel_gradient_shapes(kernel_embedding):
 
 
 @pytest.mark.parametrize("kernel_embedding", embeddings_test_list)
-def test_qkernel_gradient_values(kernel_embedding):
-    emukit_qkernel, x1, x2, N, M, D, dat_bounds = kernel_embedding
+def test_qkernel_gradient_values(kernel_embedding, request):
+    emukit_qkernel, x1, x2, N, M, D, dat_bounds = request.getfixturevalue(kernel_embedding)
 
     np.random.seed(42)
     x1 = sample_uniform(in_shape=(N, D), bounds=dat_bounds)
diff --git a/tests/emukit/quadrature/test_quadrature_models.py b/tests/emukit/quadrature/test_quadrature_models.py
index 9c4b1f5d..cfcf8d59 100644
--- a/tests/emukit/quadrature/test_quadrature_models.py
+++ b/tests/emukit/quadrature/test_quadrature_models.py
@@ -9,7 +9,6 @@
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose
-from pytest_lazyfixture import lazy_fixture
 from utils import check_grad
 
 from emukit.model_wrappers.gpy_quadrature_wrappers import BaseGaussianProcessGPy, RBFGPy
@@ -46,32 +45,6 @@ def get_base_gp():
     return BaseGaussianProcessGPy(kern=qrbf, gpy_model=gpy_model), dat
 
 
-def get_vanilla_bq_model():
-    base_gp, dat = get_base_gp()
-    return VanillaBayesianQuadrature(base_gp=base_gp, X=dat.X, Y=dat.Y)
-
-
-def get_bounded_bq_lower():
-    base_gp, dat = get_base_gp()
-    return BoundedBayesianQuadrature(base_gp=base_gp, X=dat.X, Y=dat.Y, lower_bound=dat.bound_lower)
-
-
-def get_bounded_bq_upper():
-    base_gp, dat = get_base_gp()
-    return BoundedBayesianQuadrature(base_gp=base_gp, X=dat.X, Y=dat.Y, upper_bound=dat.bound_upper)
-
-
-def get_wsabil_adapt():
-    base_gp, dat = get_base_gp()
-    return WSABIL(base_gp=base_gp, X=dat.X, Y=dat.Y, adapt_alpha=True)
-
-
-def get_wsabil_fixed():
-    base_gp, dat = get_base_gp()
-    wsabil = WSABIL(base_gp=base_gp, X=dat.X, Y=dat.Y, adapt_alpha=False)
-    return wsabil
-
-
 # === fixtures start here
 @pytest.fixture
 def data():
@@ -80,7 +53,8 @@ def data():
 
 @pytest.fixture
 def gpy_model():
-    return get_gpy_model()
+    model, _ = get_gpy_model()
+    return model
 
 
 @pytest.fixture
@@ -90,41 +64,47 @@ def base_gp():
 
 @pytest.fixture
 def vanilla_bq():
-    return get_vanilla_bq_model()
+    base_gp, dat = get_base_gp()
+    return VanillaBayesianQuadrature(base_gp=base_gp, X=dat.X, Y=dat.Y)
 
 
 @pytest.fixture
 def bounded_bq_lower():
-    return get_bounded_bq_lower()
+    base_gp, dat = get_base_gp()
+    return BoundedBayesianQuadrature(base_gp=base_gp, X=dat.X, Y=dat.Y, lower_bound=dat.bound_lower)
 
 
 @pytest.fixture
 def bounded_bq_upper():
-    return get_bounded_bq_upper()
+    base_gp, dat = get_base_gp()
+    return BoundedBayesianQuadrature(base_gp=base_gp, X=dat.X, Y=dat.Y, upper_bound=dat.bound_upper)
 
 
 @pytest.fixture
 def wsabil_adapt():
-    return get_wsabil_adapt()
+    base_gp, dat = get_base_gp()
+    return WSABIL(base_gp=base_gp, X=dat.X, Y=dat.Y, adapt_alpha=True)
 
 
 @pytest.fixture
 def wsabil_fixed():
-    return get_wsabil_fixed()
+    base_gp, dat = get_base_gp()
+    wsabil = WSABIL(base_gp=base_gp, X=dat.X, Y=dat.Y, adapt_alpha=False)
+    return wsabil
 
 
 vanilla_test_list = [
-    lazy_fixture("vanilla_bq"),
+    "vanilla_bq",
 ]
 
 bounded_test_list = [
-    lazy_fixture("bounded_bq_lower"),
-    lazy_fixture("bounded_bq_upper"),
+    "bounded_bq_lower",
+    "bounded_bq_upper",
 ]
 
 wsabi_test_list = [
-    lazy_fixture("wsabil_adapt"),
-    lazy_fixture("wsabil_fixed"),
+    "wsabil_adapt",
+    "wsabil_fixed",
 ]
 
 all_models_test_list = vanilla_test_list + bounded_test_list + wsabi_test_list
@@ -133,16 +113,18 @@ def wsabil_fixed():
 # === tests shared by all warped models start here
 
 
-@pytest.mark.parametrize("model", all_models_test_list)
-def test_warped_model_data(model, data):
+@pytest.mark.parametrize("model_name", all_models_test_list)
+def test_warped_model_data(model_name, data, request):
+    model = request.getfixturevalue(model_name)
     ABS_TOL = 1e-5
     REL_TOL = 1e-6
     assert_allclose(model.X, data.X, rtol=REL_TOL, atol=ABS_TOL)
     assert_allclose(model.Y, data.Y, rtol=REL_TOL, atol=ABS_TOL)
 
 
-@pytest.mark.parametrize("model", all_models_test_list)
-def test_warped_model_shapes(model):
+@pytest.mark.parametrize("model_name", all_models_test_list)
+def test_warped_model_shapes(model_name, request):
+    model = request.getfixturevalue(model_name)
     x = np.array([[-1, 1], [0, 0], [-2, 0.1], [-3, 4]])
     Y = np.array([[1], [2], [3]])
 
@@ -195,8 +177,9 @@ def test_warped_model_shapes(model):
     assert res[1].shape == (N, M)
 
 
-@pytest.mark.parametrize("model", all_models_test_list)
-def test_warped_model_transforms(model):
+@pytest.mark.parametrize("model_name", all_models_test_list)
+def test_warped_model_transforms(model_name, request):
+    model = request.getfixturevalue(model_name)
     Y = np.array([[1], [2], [3]])
     ABS_TOL = 1e-5
     REL_TOL = 1e-6
@@ -212,8 +195,9 @@ def test_warped_model_transforms(model):
     assert_allclose(model.inverse_transform(Y2), Y1, rtol=REL_TOL, atol=ABS_TOL)
 
 
-@pytest.mark.parametrize("model", all_models_test_list)
-def test_warped_model_gradient_values(model, data):
+@pytest.mark.parametrize("model_name", all_models_test_list)
+def test_warped_model_gradient_values(model_name, data, request):
+    model = request.getfixturevalue(model_name)
     # gradient of mean
     func = lambda z: model.predict(z)[0][:, 0]
     dfunc = lambda z: model.get_prediction_gradients(z)[0].T
@@ -226,16 +210,16 @@ def test_warped_model_gradient_values(model, data):
 
 
 @pytest.mark.parametrize(
-    "model,interval",
+    "model_name,interval",
     [
-        (vanilla_test_list[0], [0.5956279650321574, 0.6000811779371775]),
-        (bounded_test_list[0], [0.8383067891004425, 0.8417905366769567]),
-        (bounded_test_list[1], [2.977651803340788, 2.981939540780773]),
-        (wsabi_test_list[0], [1.0571955335349208, 1.0601420159245922]),
-        (wsabi_test_list[1], [0.47610638476406725, 0.48068140048609603]),
+        ("vanilla_bq", [0.5956279650321574, 0.6000811779371775]),
+        ("bounded_bq_lower", [0.8383067891004425, 0.8417905366769567]),
+        ("bounded_bq_upper", [2.977651803340788, 2.981939540780773]),
+        ("wsabil_adapt", [1.0571955335349208, 1.0601420159245922]),
+        ("wsabil_fixed", [0.47610638476406725, 0.48068140048609603]),
     ],
 )
-def test_warped_model_integrate_mean(model, interval):
+def test_warped_model_integrate_mean(model_name, interval, request):
     # Both outputs of the model.intgerate() method are analytic integrals.
     # To test their values we check if they lie in the confidence interva of an MC estimator.
     # These intervals were computed as follows: the mean model.predict (first argument) was integrated by
@@ -243,21 +227,23 @@ def test_warped_model_integrate_mean(model, interval):
     # times. The intervals show mean\pm 3 std of the 100 integrals obtained by sampling. There might be a very small
     # chance that the true integrals lies outside the specified intervals.
     # See file "ground_truth_integrals_methods.py" for details.
+    model = request.getfixturevalue(model_name)
     res = model.integrate()[0]
     assert interval[0] < res < interval[1]
 
 
 @pytest.mark.parametrize(
-    "model,interval",
+    "model_name,interval",
     [
-        (vanilla_test_list[0], [0.09859906877945852, 0.11181285735935843]),
-        (bounded_test_list[0], None),
-        (bounded_test_list[1], None),
-        (wsabi_test_list[0], None),
-        (wsabi_test_list[1], None),
+        ("vanilla_bq", [0.09859906877945852, 0.11181285735935843]),
+        ("bounded_bq_lower", None),
+        ("bounded_bq_upper", None),
+        ("wsabil_adapt", None),
+        ("wsabil_fixed", None),
     ],
 )
-def test_warped_model_integrate_variance(model, interval):
+def test_warped_model_integrate_variance(model_name, interval, request):
+    model = request.getfixturevalue(model_name)
     # See test_warped_model_integrate_mean on how the intervals were computed
     res = model.integrate()[1]
 
@@ -284,7 +270,6 @@ def test_bounded_bq_correct_bound(data, bounded_bq_lower, bounded_bq_upper):
 
 
 def test_bounded_bq_raises(gpy_model):
-    gpy_model, _ = gpy_model
     measure = LebesgueMeasure.from_bounds(gpy_model.X.shape[1] * [(0, 1)], normalized=False)
     qrbf = QuadratureRBFLebesgueMeasure(RBFGPy(gpy_model.kern), measure=measure)
     base_gp_wrong_kernel = BaseGaussianProcessGPy(kern=qrbf, gpy_model=gpy_model)
diff --git a/tests/emukit/quadrature/test_warpings.py b/tests/emukit/quadrature/test_warpings.py
index 3c1999f5..db107df7 100644
--- a/tests/emukit/quadrature/test_warpings.py
+++ b/tests/emukit/quadrature/test_warpings.py
@@ -7,15 +7,10 @@
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose
-from pytest_lazyfixture import lazy_fixture
 
 from emukit.quadrature.methods.warpings import IdentityWarping, SquareRootWarping
 
 
-def create_fixture_parameters():
-    return [pytest.param(lazy_fixture(warping.name), id=warping.name) for warping in warpings]
-
-
 @pytest.fixture
 def identity_warping():
     return IdentityWarping()
@@ -33,11 +28,10 @@ def inverted_squarerroot_warping():
     return SquareRootWarping(offset=offset, is_inverted=True)
 
 
-warpings_tuple = namedtuple("WarpingTest", ["name"])
 warpings = [
-    warpings_tuple("identity_warping"),
-    warpings_tuple("squarerroot_warping"),
-    warpings_tuple("inverted_squarerroot_warping"),
+    "identity_warping",
+    "squarerroot_warping",
+    "inverted_squarerroot_warping",
 ]
 
 
@@ -45,15 +39,17 @@ def inverted_squarerroot_warping():
 ATOL = 1e-6
 
 
-@pytest.mark.parametrize("warping", create_fixture_parameters())
-def test_warping_shapes(warping):
+@pytest.mark.parametrize("warping_name", warpings)
+def test_warping_shapes(warping_name, request):
+    warping = request.getfixturevalue(warping_name)
     Y = np.ones([5, 1])
     assert warping.transform(Y).shape == Y.shape
     assert warping.inverse_transform(Y).shape == Y.shape
 
 
-@pytest.mark.parametrize("warping", create_fixture_parameters())
-def test_warping_values(warping):
+@pytest.mark.parametrize("warping_name", warpings)
+def test_warping_values(warping_name, request):
+    warping = request.getfixturevalue(warping_name)
     np.random.seed(42)
     Y = np.random.rand(5, 1)
 
diff --git a/tests/emukit/test_acquisitions.py b/tests/emukit/test_acquisitions.py
index 83a763c0..e8792414 100644
--- a/tests/emukit/test_acquisitions.py
+++ b/tests/emukit/test_acquisitions.py
@@ -9,7 +9,6 @@
 
 import numpy as np
 import pytest
-import pytest_lazyfixture
 from scipy.optimize import check_grad
 
 from emukit.bayesian_optimization.acquisitions import (
@@ -160,7 +159,7 @@ def MUMBO_acquisition(gpy_model):
 
 # Helpers for creating parameterized fixtures
 def create_acquisition_fixture_parameters():
-    return [pytest.param(pytest_lazyfixture.lazy_fixture(acq.name), id=acq.name) for acq in acquisition_tests]
+    return [acq.name for acq in acquisition_tests]
 
 
 def create_gradient_acquisition_fixtures():
@@ -168,15 +167,14 @@ def create_gradient_acquisition_fixtures():
     parameters = []
     for acquisition in acquisition_tests:
         if acquisition.has_gradients:
-            acquisition_name = acquisition.name
-            lazy_fixture = pytest_lazyfixture.lazy_fixture(acquisition.name)
-            parameters.append(pytest.param(lazy_fixture, acquisition.rmse_gradient_tolerance, id=acquisition_name))
+            parameters.append(pytest.param(acquisition.name, acquisition.rmse_gradient_tolerance, id=acquisition.name))
     return parameters
 
 
 # Tests
-@pytest.mark.parametrize("acquisition", create_acquisition_fixture_parameters())
-def test_acquisition_evaluate_shape(acquisition, n_dims):
+@pytest.mark.parametrize("acquisition_name", create_acquisition_fixture_parameters())
+def test_acquisition_evaluate_shape(acquisition_name, n_dims, request):
+    acquisition = request.getfixturevalue(acquisition_name)
     x = np.random.rand(1, n_dims)
     acquisition_value = acquisition.evaluate(x)
     assert acquisition_value.shape == (1, 1)
@@ -186,8 +184,9 @@ def test_acquisition_evaluate_shape(acquisition, n_dims):
     assert acquisition_value.shape == (10, 1)
 
 
-@pytest.mark.parametrize(("acquisition", "tol"), create_gradient_acquisition_fixtures())
-def test_acquisition_gradient_computation(acquisition, n_dims, tol):
+@pytest.mark.parametrize(("acquisition_name", "tol"), create_gradient_acquisition_fixtures())
+def test_acquisition_gradient_computation(acquisition_name, n_dims, tol, request):
+    acquisition = request.getfixturevalue(acquisition_name)
     rng = np.random.RandomState(43)
     x_test = rng.rand(10, n_dims)
 
@@ -199,8 +198,9 @@ def test_acquisition_gradient_computation(acquisition, n_dims, tol):
         assert err < tol
 
 
-@pytest.mark.parametrize(("acquisition", "tol"), create_gradient_acquisition_fixtures())
-def test_acquisition_gradient_shapes(acquisition, n_dims, tol):
+@pytest.mark.parametrize(("acquisition_name", "tol"), create_gradient_acquisition_fixtures())
+def test_acquisition_gradient_shapes(acquisition_name, n_dims, tol, request):
+    acquisition = request.getfixturevalue(acquisition_name)
     rng = np.random.RandomState(43)
 
     x_test = rng.rand(1, n_dims)