diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml
index 43e7d988e..4227baf14 100644
--- a/.github/workflows/build_docs.yml
+++ b/.github/workflows/build_docs.yml
@@ -28,6 +28,8 @@ jobs:
         run: sudo apt-get install -y pandoc
       - name: Install Poetry
         uses: snok/install-poetry@v1
+        with:
+          version: 'latest'
       - name: Install package and test dependencies
         run: |
           poetry install --with dev,nbtools
diff --git a/.github/workflows/format_and_lint.yml b/.github/workflows/format_and_lint.yml
index 7a870481a..a7fffe021 100644
--- a/.github/workflows/format_and_lint.yml
+++ b/.github/workflows/format_and_lint.yml
@@ -19,6 +19,8 @@ jobs:
         python-version: "3.9"
     - name: Install Poetry
       uses: snok/install-poetry@v1
+      with:
+        version: 'latest'
     - name: Install dependencies
       run: |
         poetry install --with dev
diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
index 8501d8cd6..6ae4602ca 100644
--- a/.github/workflows/run_tests.yml
+++ b/.github/workflows/run_tests.yml
@@ -19,6 +19,9 @@ jobs:
       matrix:
         python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
         numpy-version: [">=1.25,<2", ">=2"]
+        exclude:
+          - python-version: "3.13"
+            numpy-version: ">=1.25,<2" # numpy<2 is not supported on Python 3.13
 
     steps:
     - uses: actions/checkout@v3
@@ -26,14 +29,28 @@ jobs:
       uses: actions/setup-python@v3
       with:
         python-version: ${{ matrix.python-version }}
+    - name: Install the latest version of uv
+      uses: astral-sh/setup-uv@v5
+      with:
+        version: "latest"
     - name: Install Poetry
       uses: snok/install-poetry@v1
+      with:
+        version: 'latest'
     - name: Install test dependencies
       run: |
-        poetry add "numpy${{ matrix.numpy-version }}"
-        poetry install --with dev,nbtools
+        poetry self add poetry-plugin-export
+        poetry export -f requirements.txt --with dev,nbtools --without-hashes --output requirements-dev.txt
+        echo "numpy${{ matrix.numpy-version }}" >> constraints.txt
+        uv pip compile requirements-dev.txt --output-file requirements.txt \
+          --python-version ${{ matrix.python-version }} \
+          --override constraints.txt
+        poetry run pip install -r requirements.txt
+        poetry install --only-root
     - name: Test with pytest
       run: |
         poetry run pytest --cov-report xml --cov=bayes_opt/
     - name: Upload coverage to Codecov
-      uses: codecov/codecov-action@v3
+      uses: codecov/codecov-action@v4
+      with:
+          token: ${{ secrets.CODECOV_TOKEN }}
diff --git a/.gitignore b/.gitignore
index 786d95c87..4b159c12f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -36,4 +36,8 @@ docsrc/*.ipynb
 docsrc/static/*
 docsrc/README.md
 
-poetry.lock
\ No newline at end of file
+poetry.lock
+
+# Add log files and optimizer state files to gitignore
+examples/logs.log
+examples/optimizer_state.json
diff --git a/bayes_opt/acquisition.py b/bayes_opt/acquisition.py
index 167bd5dc0..a81a6df13 100644
--- a/bayes_opt/acquisition.py
+++ b/bayes_opt/acquisition.py
@@ -69,6 +69,33 @@ def __init__(self, random_state: int | RandomState | None = None) -> None:
             self.random_state = RandomState()
         self.i = 0
 
+    def _serialize_random_state(self) -> dict | None:
+        """Convert random state to JSON serializable format."""
+        if self.random_state is not None:
+            state = self.random_state.get_state()
+            return {
+                "bit_generator": state[0],
+                "state": state[1].tolist(),  # Convert numpy array to list
+                "pos": state[2],
+                "has_gauss": state[3],
+                "cached_gaussian": state[4],
+            }
+        return None
+
+    def _deserialize_random_state(self, state_dict: dict | None) -> None:
+        """Restore random state from JSON serializable format."""
+        if state_dict is not None:
+            if self.random_state is None:
+                self.random_state = RandomState()
+            state = (
+                state_dict["bit_generator"],
+                np.array(state_dict["state"], dtype=np.uint32),
+                state_dict["pos"],
+                state_dict["has_gauss"],
+                state_dict["cached_gaussian"],
+            )
+            self.random_state.set_state(state)
+
     @abc.abstractmethod
     def base_acq(self, *args: Any, **kwargs: Any) -> NDArray[Float]:
         """Provide access to the base acquisition function."""
@@ -82,6 +109,34 @@ def _fit_gp(self, gp: GaussianProcessRegressor, target_space: TargetSpace) -> No
             if target_space.constraint is not None:
                 target_space.constraint.fit(target_space.params, target_space._constraint_values)
 
+    def get_acquisition_params(self) -> dict[str, Any]:
+        """
+        Get the parameters of the acquisition function.
+
+        Returns
+        -------
+        dict
+            The parameters of the acquisition function.
+        """
+        error_msg = (
+            "Custom AcquisitionFunction subclasses must implement their own get_acquisition_params method."
+        )
+        raise NotImplementedError(error_msg)
+
+    def set_acquisition_params(self, **params) -> None:
+        """
+        Set the parameters of the acquisition function.
+
+        Parameters
+        ----------
+        **params : dict
+            The parameters of the acquisition function.
+        """
+        error_msg = (
+            "Custom AcquisitionFunction subclasses must implement their own set_acquisition_params method."
+        )
+        raise NotImplementedError(error_msg)
+
     def suggest(
         self,
         gp: GaussianProcessRegressor,
@@ -454,6 +509,34 @@ def decay_exploration(self) -> None:
         ):
             self.kappa = self.kappa * self.exploration_decay
 
+    def get_acquisition_params(self) -> dict:
+        """Get the current acquisition function parameters.
+
+        Returns
+        -------
+        dict
+            Dictionary containing the current acquisition function parameters.
+        """
+        return {
+            "kappa": self.kappa,
+            "exploration_decay": self.exploration_decay,
+            "exploration_decay_delay": self.exploration_decay_delay,
+            "random_state": self._serialize_random_state(),
+        }
+
+    def set_acquisition_params(self, params: dict) -> None:
+        """Set the acquisition function parameters.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary containing the acquisition function parameters.
+        """
+        self.kappa = params["kappa"]
+        self.exploration_decay = params["exploration_decay"]
+        self.exploration_decay_delay = params["exploration_decay_delay"]
+        self._deserialize_random_state(params["random_state"])
+
 
 class ProbabilityOfImprovement(AcquisitionFunction):
     r"""Probability of Improvement acqusition function.
@@ -587,6 +670,34 @@ def decay_exploration(self) -> None:
         ):
             self.xi = self.xi * self.exploration_decay
 
+    def get_acquisition_params(self) -> dict:
+        """Get the current acquisition function parameters.
+
+        Returns
+        -------
+        dict
+            Dictionary containing the current acquisition function parameters.
+        """
+        return {
+            "xi": self.xi,
+            "exploration_decay": self.exploration_decay,
+            "exploration_decay_delay": self.exploration_decay_delay,
+            "random_state": self._serialize_random_state(),
+        }
+
+    def set_acquisition_params(self, params: dict) -> None:
+        """Set the acquisition function parameters.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary containing the acquisition function parameters.
+        """
+        self.xi = params["xi"]
+        self.exploration_decay = params["exploration_decay"]
+        self.exploration_decay_delay = params["exploration_decay_delay"]
+        self._deserialize_random_state(params["random_state"])
+
 
 class ExpectedImprovement(AcquisitionFunction):
     r"""Expected Improvement acqusition function.
@@ -728,6 +839,34 @@ def decay_exploration(self) -> None:
         ):
             self.xi = self.xi * self.exploration_decay
 
+    def get_acquisition_params(self) -> dict:
+        """Get the current acquisition function parameters.
+
+        Returns
+        -------
+        dict
+            Dictionary containing the current acquisition function parameters.
+        """
+        return {
+            "xi": self.xi,
+            "exploration_decay": self.exploration_decay,
+            "exploration_decay_delay": self.exploration_decay_delay,
+            "random_state": self._serialize_random_state(),
+        }
+
+    def set_acquisition_params(self, params: dict) -> None:
+        """Set the acquisition function parameters.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary containing the acquisition function parameters.
+        """
+        self.xi = params["xi"]
+        self.exploration_decay = params["exploration_decay"]
+        self.exploration_decay_delay = params["exploration_decay_delay"]
+        self._deserialize_random_state(params["random_state"])
+
 
 class ConstantLiar(AcquisitionFunction):
     """Constant Liar acquisition function.
@@ -918,6 +1057,38 @@ def suggest(
 
         return x_max
 
+    def get_acquisition_params(self) -> dict:
+        """Get the current acquisition function parameters.
+
+        Returns
+        -------
+        dict
+            Dictionary containing the current acquisition function parameters.
+        """
+        return {
+            "dummies": [dummy.tolist() for dummy in self.dummies],
+            "base_acquisition_params": self.base_acquisition.get_acquisition_params(),
+            "strategy": self.strategy,
+            "atol": self.atol,
+            "rtol": self.rtol,
+            "random_state": self._serialize_random_state(),
+        }
+
+    def set_acquisition_params(self, params: dict) -> None:
+        """Set the acquisition function parameters.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary containing the acquisition function parameters.
+        """
+        self.dummies = [np.array(dummy) for dummy in params["dummies"]]
+        self.base_acquisition.set_acquisition_params(params["base_acquisition_params"])
+        self.strategy = params["strategy"]
+        self.atol = params["atol"]
+        self.rtol = params["rtol"]
+        self._deserialize_random_state(params["random_state"])
+
 
 class GPHedge(AcquisitionFunction):
     """GPHedge acquisition function.
@@ -1035,3 +1206,38 @@ def suggest(
         self.previous_candidates = np.array(x_max)
         idx = self._sample_idx_from_softmax_gains()
         return x_max[idx]
+
+    def get_acquisition_params(self) -> dict:
+        """Get the current acquisition function parameters.
+
+        Returns
+        -------
+        dict
+            Dictionary containing the current acquisition function parameters.
+        """
+        return {
+            "base_acquisitions_params": [acq.get_acquisition_params() for acq in self.base_acquisitions],
+            "gains": self.gains.tolist(),
+            "previous_candidates": self.previous_candidates.tolist()
+            if self.previous_candidates is not None
+            else None,
+            "gphedge_random_state": self._serialize_random_state(),
+        }
+
+    def set_acquisition_params(self, params: dict) -> None:
+        """Set the acquisition function parameters.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary containing the acquisition function parameters.
+        """
+        for acq, acq_params in zip(self.base_acquisitions, params["base_acquisitions_params"]):
+            acq.set_acquisition_params(acq_params)
+
+        self.gains = np.array(params["gains"])
+        self.previous_candidates = (
+            np.array(params["previous_candidates"]) if params["previous_candidates"] is not None else None
+        )
+
+        self._deserialize_random_state(params["gphedge_random_state"])
diff --git a/bayes_opt/bayesian_optimization.py b/bayes_opt/bayesian_optimization.py
index d7f2e4035..90c22d647 100644
--- a/bayes_opt/bayesian_optimization.py
+++ b/bayes_opt/bayesian_optimization.py
@@ -6,11 +6,15 @@
 
 from __future__ import annotations
 
+import json
 from collections import deque
+from os import PathLike
+from pathlib import Path
 from typing import TYPE_CHECKING, Any
 from warnings import warn
 
 import numpy as np
+from scipy.optimize import NonlinearConstraint
 from sklearn.gaussian_process import GaussianProcessRegressor
 from sklearn.gaussian_process.kernels import Matern
 
@@ -28,7 +32,6 @@
 
     from numpy.random import RandomState
     from numpy.typing import NDArray
-    from scipy.optimize import NonlinearConstraint
 
     from bayes_opt.acquisition import AcquisitionFunction
     from bayes_opt.constraint import ConstraintModel
@@ -356,3 +359,113 @@ def set_gp_params(self, **params: Any) -> None:
         if "kernel" in params:
             params["kernel"] = wrap_kernel(kernel=params["kernel"], transform=self._space.kernel_transform)
         self._gp.set_params(**params)
+
+    def save_state(self, path: str | PathLike[str]) -> None:
+        """Save complete state for reconstruction of the optimizer.
+
+        Parameters
+        ----------
+        path : str or PathLike
+            Path to save the optimization state
+
+        Raises
+        ------
+        ValueError
+            If attempting to save state before collecting any samples.
+        """
+        if len(self._space) == 0:
+            msg = (
+                "Cannot save optimizer state before collecting any samples. "
+                "Please probe or register at least one point before saving."
+            )
+            raise ValueError(msg)
+
+        random_state = None
+        if self._random_state is not None:
+            state_tuple = self._random_state.get_state()
+            random_state = {
+                "bit_generator": state_tuple[0],
+                "state": state_tuple[1].tolist(),
+                "pos": state_tuple[2],
+                "has_gauss": state_tuple[3],
+                "cached_gaussian": state_tuple[4],
+            }
+
+        # Get constraint values if they exist
+        constraint_values = self._space._constraint_values.tolist() if self.is_constrained else None
+        acquisition_params = self._acquisition_function.get_acquisition_params()
+        state = {
+            "pbounds": {key: self._space._bounds[i].tolist() for i, key in enumerate(self._space.keys)},
+            # Add current transformed bounds if using bounds transformer
+            "transformed_bounds": (self._space.bounds.tolist() if self._bounds_transformer else None),
+            "keys": self._space.keys,
+            "params": np.array(self._space.params).tolist(),
+            "target": self._space.target.tolist(),
+            "constraint_values": constraint_values,
+            "gp_params": {
+                "kernel": self._gp.kernel.get_params(),
+                "alpha": self._gp.alpha,
+                "normalize_y": self._gp.normalize_y,
+                "n_restarts_optimizer": self._gp.n_restarts_optimizer,
+            },
+            "allow_duplicate_points": self._allow_duplicate_points,
+            "verbose": self._verbose,
+            "random_state": random_state,
+            "acquisition_params": acquisition_params,
+        }
+
+        with Path(path).open("w") as f:
+            json.dump(state, f, indent=2)
+
+    def load_state(self, path: str | PathLike[str]) -> None:
+        """Load optimizer state from a JSON file.
+
+        Parameters
+        ----------
+        path : str or PathLike
+            Path to the JSON file containing the optimizer state.
+        """
+        with Path(path).open("r") as file:
+            state = json.load(file)
+
+        params_array = np.asarray(state["params"], dtype=np.float64)
+        target_array = np.asarray(state["target"], dtype=np.float64)
+        constraint_array = (
+            np.array(state["constraint_values"]) if state["constraint_values"] is not None else None
+        )
+
+        for i in range(len(params_array)):
+            params = self._space.array_to_params(params_array[i])
+            target = target_array[i]
+            constraint = constraint_array[i] if constraint_array is not None else None
+            self.register(params=params, target=target, constraint_value=constraint)
+
+        self._acquisition_function.set_acquisition_params(state["acquisition_params"])
+
+        if state.get("transformed_bounds") and self._bounds_transformer:
+            new_bounds = {
+                key: bounds for key, bounds in zip(self._space.keys, np.array(state["transformed_bounds"]))
+            }
+            self._space.set_bounds(new_bounds)
+            self._bounds_transformer.initialize(self._space)
+
+        self._gp.set_params(**state["gp_params"])
+        if isinstance(self._gp.kernel, dict):
+            kernel_params = self._gp.kernel
+            self._gp.kernel = Matern(
+                length_scale=kernel_params["length_scale"],
+                length_scale_bounds=tuple(kernel_params["length_scale_bounds"]),
+                nu=kernel_params["nu"],
+            )
+
+        self._gp.fit(self._space.params, self._space.target)
+
+        if state["random_state"] is not None:
+            random_state_tuple = (
+                state["random_state"]["bit_generator"],
+                np.array(state["random_state"]["state"], dtype=np.uint32),
+                state["random_state"]["pos"],
+                state["random_state"]["has_gauss"],
+                state["random_state"]["cached_gaussian"],
+            )
+            self._random_state.set_state(random_state_tuple)
diff --git a/examples/acquisition_functions.ipynb b/examples/acquisition_functions.ipynb
index 53e693754..b3a84684b 100644
--- a/examples/acquisition_functions.ipynb
+++ b/examples/acquisition_functions.ipynb
@@ -54,14 +54,14 @@
     "        'Gaussian Process and Utility Function After {} Steps'.format(steps),\n",
     "        fontsize=30\n",
     "    )\n",
-    "    \n",
-    "    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) \n",
+    "\n",
+    "    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])\n",
     "    axis = plt.subplot(gs[0])\n",
     "    acq = plt.subplot(gs[1])\n",
-    "    \n",
+    "\n",
     "    x_obs = np.array([[res[\"params\"][\"x\"]] for res in optimizer.res])\n",
     "    y_obs = np.array([res[\"target\"] for res in optimizer.res])\n",
-    "    \n",
+    "\n",
     "    acquisition_function_._fit_gp(optimizer._gp, optimizer._space)\n",
     "    mu, sigma = posterior(optimizer, x)\n",
     "\n",
@@ -69,10 +69,10 @@
     "    axis.plot(x_obs.flatten(), y_obs, 'D', markersize=8, label=u'Observations', color='r')\n",
     "    axis.plot(x, mu, '--', color='k', label='Prediction')\n",
     "\n",
-    "    axis.fill(np.concatenate([x, x[::-1]]), \n",
+    "    axis.fill(np.concatenate([x, x[::-1]]),\n",
     "              np.concatenate([mu - 1.9600 * sigma, (mu + 1.9600 * sigma)[::-1]]),\n",
     "        alpha=.6, fc='c', ec='None', label='95% confidence interval')\n",
-    "    \n",
+    "\n",
     "    axis.set_xlim((-2, 10))\n",
     "    axis.set_ylim((None, None))\n",
     "    axis.set_ylabel('f(x)', fontdict={'size':20})\n",
@@ -82,13 +82,13 @@
     "    x = x.flatten()\n",
     "\n",
     "    acq.plot(x, utility, label='Utility Function', color='purple')\n",
-    "    acq.plot(x[np.argmax(utility)], np.max(utility), '*', markersize=15, \n",
+    "    acq.plot(x[np.argmax(utility)], np.max(utility), '*', markersize=15,\n",
     "             label=u'Next Best Guess', markerfacecolor='gold', markeredgecolor='k', markeredgewidth=1)\n",
     "    acq.set_xlim((-2, 10))\n",
     "    #acq.set_ylim((0, np.max(utility) + 0.5))\n",
     "    acq.set_ylabel('Utility', fontdict={'size':20})\n",
     "    acq.set_xlabel('x', fontdict={'size':20})\n",
-    "    \n",
+    "\n",
     "    axis.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)\n",
     "    acq.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)\n",
     "    return fig, fig.axes"
@@ -110,7 +110,7 @@
     "class GreedyAcquisition(acquisition.AcquisitionFunction):\n",
     "    def __init__(self, random_state=None):\n",
     "        super().__init__(random_state)\n",
-    "    \n",
+    "\n",
     "    def base_acq(self, mean, std):\n",
     "        return mean # disregard std"
    ]
@@ -357,11 +357,17 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
-   "source": []
+   "source": [
+    "### Saving and loading with custom acquisition functions\n",
+    "\n",
+    "If saving and loading is desired functionality, the acquisition function needs to have a `get_acquisition_params` and `set_acquisition_params` method.\n",
+    "\n",
+    "Default acquisition functions have a `get_acquisition_params` and `set_acquisition_params` method that is used to save and load the acquisition function state when the optimizer is being saved and loaded. When writing and using a custom acquisition function, those methods need to also be defined. In order to have perfect replicability, all parameters that the acquisition function depends on need to be saved and loaded. The get_acquisition_params method should return a dictionary containing all acquisition function parameters. The set_acquisition_params method should take that same dictionary containing the acquisition function parameters and update the respective parameters.\n",
+    "\n",
+    "_serialize_random_state and _deserialize_random_state are provided as part of the base class and should be included in the dictionary returned by `get_acquisition_params`. See any of the default acquisition functions for reference."
+   ]
   },
   {
    "cell_type": "code",
diff --git a/examples/basic-tour.ipynb b/examples/basic-tour.ipynb
index 4ecd83296..959782aeb 100644
--- a/examples/basic-tour.ipynb
+++ b/examples/basic-tour.ipynb
@@ -106,11 +106,11 @@
      "text": [
       "|   iter    |  target   |     x     |     y     |\n",
       "-------------------------------------------------\n",
-      "| \u001b[0m1        \u001b[0m | \u001b[0m-7.135   \u001b[0m | \u001b[0m2.834    \u001b[0m | \u001b[0m1.322    \u001b[0m |\n",
-      "| \u001b[0m2        \u001b[0m | \u001b[0m-7.78    \u001b[0m | \u001b[0m2.0      \u001b[0m | \u001b[0m-1.186   \u001b[0m |\n",
-      "| \u001b[95m3        \u001b[0m | \u001b[95m-6.967   \u001b[0m | \u001b[95m2.582    \u001b[0m | \u001b[95m-0.1396  \u001b[0m |\n",
-      "| \u001b[0m4        \u001b[0m | \u001b[0m-12.24   \u001b[0m | \u001b[0m2.015    \u001b[0m | \u001b[0m-2.029   \u001b[0m |\n",
-      "| \u001b[0m5        \u001b[0m | \u001b[0m-18.0    \u001b[0m | \u001b[0m3.302    \u001b[0m | \u001b[0m-1.846   \u001b[0m |\n",
+      "| \u001b[39m1        \u001b[39m | \u001b[39m-7.135   \u001b[39m | \u001b[39m2.8340440\u001b[39m | \u001b[39m1.3219469\u001b[39m |\n",
+      "| \u001b[39m2        \u001b[39m | \u001b[39m-7.78    \u001b[39m | \u001b[39m2.0002287\u001b[39m | \u001b[39m-1.186004\u001b[39m |\n",
+      "| \u001b[39m3        \u001b[39m | \u001b[39m-7.157   \u001b[39m | \u001b[39m2.8375977\u001b[39m | \u001b[39m1.3238498\u001b[39m |\n",
+      "| \u001b[35m4        \u001b[39m | \u001b[35m-6.633   \u001b[39m | \u001b[35m2.7487090\u001b[39m | \u001b[35m1.2790562\u001b[39m |\n",
+      "| \u001b[35m5        \u001b[39m | \u001b[35m-5.751   \u001b[39m | \u001b[35m2.5885326\u001b[39m | \u001b[35m1.2246876\u001b[39m |\n",
       "=================================================\n"
      ]
     }
@@ -138,7 +138,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "{'target': -6.966506881014352, 'params': {'x': 2.5822074982517598, 'y': -0.13961016009280103}}\n"
+      "{'target': np.float64(-5.750985875689304), 'params': {'x': np.float64(2.5885326650623566), 'y': np.float64(1.2246876000015976)}}\n"
      ]
     }
    ],
@@ -163,15 +163,15 @@
      "output_type": "stream",
      "text": [
       "Iteration 0: \n",
-      "\t{'target': -7.135455292718879, 'params': {'x': 2.8340440094051482, 'y': 1.3219469606529486}}\n",
+      "\t{'target': np.float64(-7.135455292718879), 'params': {'x': np.float64(2.8340440094051482), 'y': np.float64(1.3219469606529486)}}\n",
       "Iteration 1: \n",
-      "\t{'target': -7.779531005607566, 'params': {'x': 2.0002287496346898, 'y': -1.1860045642089614}}\n",
+      "\t{'target': np.float64(-7.779531005607566), 'params': {'x': np.float64(2.0002287496346898), 'y': np.float64(-1.1860045642089614)}}\n",
       "Iteration 2: \n",
-      "\t{'target': -6.966506881014352, 'params': {'x': 2.5822074982517598, 'y': -0.13961016009280103}}\n",
+      "\t{'target': np.float64(-7.156839989425082), 'params': {'x': np.float64(2.8375977943744273), 'y': np.float64(1.3238498831039895)}}\n",
       "Iteration 3: \n",
-      "\t{'target': -12.235835023240657, 'params': {'x': 2.0154397119682423, 'y': -2.0288343947238228}}\n",
+      "\t{'target': np.float64(-6.633273772355583), 'params': {'x': np.float64(2.7487090390562576), 'y': np.float64(1.2790562505410115)}}\n",
       "Iteration 4: \n",
-      "\t{'target': -17.99963711795217, 'params': {'x': 3.301617813728339, 'y': -1.8458666395359906}}\n"
+      "\t{'target': np.float64(-5.750985875689304), 'params': {'x': np.float64(2.5885326650623566), 'y': np.float64(1.2246876000015976)}}\n"
      ]
     }
    ],
@@ -209,11 +209,11 @@
      "text": [
       "|   iter    |  target   |     x     |     y     |\n",
       "-------------------------------------------------\n",
-      "| \u001b[95m6        \u001b[0m | \u001b[95m-1.378   \u001b[0m | \u001b[95m1.229    \u001b[0m | \u001b[95m0.06916  \u001b[0m |\n",
-      "| \u001b[95m7        \u001b[0m | \u001b[95m-0.9471  \u001b[0m | \u001b[95m-0.9434  \u001b[0m | \u001b[95m-0.02816 \u001b[0m |\n",
-      "| \u001b[0m8        \u001b[0m | \u001b[0m-3.651   \u001b[0m | \u001b[0m-1.168   \u001b[0m | \u001b[0m2.813    \u001b[0m |\n",
-      "| \u001b[0m9        \u001b[0m | \u001b[0m-17.13   \u001b[0m | \u001b[0m-2.0     \u001b[0m | \u001b[0m-2.758   \u001b[0m |\n",
-      "| \u001b[95m10       \u001b[0m | \u001b[95m0.8427   \u001b[0m | \u001b[95m-0.193   \u001b[0m | \u001b[95m0.6535   \u001b[0m |\n",
+      "| \u001b[35m6        \u001b[39m | \u001b[35m-4.438   \u001b[39m | \u001b[35m2.3269441\u001b[39m | \u001b[35m1.1533794\u001b[39m |\n",
+      "| \u001b[35m7        \u001b[39m | \u001b[35m-2.42    \u001b[39m | \u001b[35m1.8477442\u001b[39m | \u001b[35m0.9230233\u001b[39m |\n",
+      "| \u001b[35m8        \u001b[39m | \u001b[35m-0.2088  \u001b[39m | \u001b[35m1.0781674\u001b[39m | \u001b[35m1.2152869\u001b[39m |\n",
+      "| \u001b[35m9        \u001b[39m | \u001b[35m0.7797   \u001b[39m | \u001b[35m-0.298812\u001b[39m | \u001b[35m1.3619705\u001b[39m |\n",
+      "| \u001b[39m10       \u001b[39m | \u001b[39m-3.391   \u001b[39m | \u001b[39m-0.655060\u001b[39m | \u001b[39m2.9904883\u001b[39m |\n",
       "=================================================\n"
      ]
     }
@@ -295,8 +295,8 @@
      "text": [
       "|   iter    |  target   |     x     |     y     |\n",
       "-------------------------------------------------\n",
-      "| \u001b[0m11       \u001b[0m | \u001b[0m0.66     \u001b[0m | \u001b[0m0.5      \u001b[0m | \u001b[0m0.7      \u001b[0m |\n",
-      "| \u001b[0m12       \u001b[0m | \u001b[0m0.1      \u001b[0m | \u001b[0m-0.3     \u001b[0m | \u001b[0m0.1      \u001b[0m |\n",
+      "| \u001b[39m11       \u001b[39m | \u001b[39m0.66     \u001b[39m | \u001b[39m0.5      \u001b[39m | \u001b[39m0.7      \u001b[39m |\n",
+      "| \u001b[39m12       \u001b[39m | \u001b[39m0.1      \u001b[39m | \u001b[39m-0.3     \u001b[39m | \u001b[39m0.1      \u001b[39m |\n",
       "=================================================\n"
      ]
     }
@@ -309,176 +309,68 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## 4. Saving, loading and restarting\n",
+    "## 4. Saving and loading the optimizer\n",
     "\n",
-    "By default you can follow the progress of your optimization by setting `verbose>0` when instantiating the `BayesianOptimization` object. If you need more control over logging/alerting you will need to use an observer. For more information about observers checkout the advanced tour notebook. Here we will only see how to use the native `JSONLogger` object to save to and load progress from files.\n",
+    "The optimizer state can be saved to a file and loaded from a file. This is useful for continuing an optimization from a previous state, or for analyzing the optimization history without running the optimizer again.\n",
     "\n",
-    "### 4.1 Saving progress"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 14,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from bayes_opt.logger import JSONLogger\n",
-    "from bayes_opt.event import Events"
+    "Note: if you are using your own custom acquisition function, you will need to save and load the acquisition function state as well. This is done by calling the `get_acquisition_params` and `set_acquisition_params` methods of the acquisition function. See the acquisition function documentation for more information."
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The observer paradigm works by:\n",
-    "1. Instantiating an observer object.\n",
-    "2. Tying the observer object to a particular event fired by an optimizer.\n",
-    "\n",
-    "The `BayesianOptimization` object fires a number of internal events during optimization, in particular, every time it probes the function and obtains a new parameter-target combination it will fire an `Events.OPTIMIZATION_STEP` event, which our logger will listen to.\n",
+    "### 4.1 Saving the optimizer state\n",
     "\n",
-    "**Caveat:** The logger will not look back at previously probed points."
+    "The optimizer state can be saved to a file using the `save_state` method.\n",
+    "optimizer.save_state(\"./optimizer_state.json\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
-    "logger = JSONLogger(path=\"./logs.log\")\n",
-    "optimizer.subscribe(Events.OPTIMIZATION_STEP, logger)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 16,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "|   iter    |  target   |     x     |     y     |\n",
-      "-------------------------------------------------\n",
-      "| \u001b[0m13       \u001b[0m | \u001b[0m-12.48   \u001b[0m | \u001b[0m-1.266   \u001b[0m | \u001b[0m-2.446   \u001b[0m |\n",
-      "| \u001b[0m14       \u001b[0m | \u001b[0m-3.854   \u001b[0m | \u001b[0m-1.069   \u001b[0m | \u001b[0m-0.9266  \u001b[0m |\n",
-      "| \u001b[95m15       \u001b[0m | \u001b[95m0.967    \u001b[0m | \u001b[95m0.04749  \u001b[0m | \u001b[95m1.175    \u001b[0m |\n",
-      "| \u001b[95m16       \u001b[0m | \u001b[95m0.9912   \u001b[0m | \u001b[95m0.07374  \u001b[0m | \u001b[95m0.9421   \u001b[0m |\n",
-      "| \u001b[0m17       \u001b[0m | \u001b[0m-3.565   \u001b[0m | \u001b[0m0.7821   \u001b[0m | \u001b[0m2.988    \u001b[0m |\n",
-      "=================================================\n"
-     ]
-    }
-   ],
-   "source": [
-    "optimizer.maximize(\n",
-    "    init_points=2,\n",
-    "    n_iter=3,\n",
-    ")"
+    "optimizer.save_state(\"optimizer_state.json\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### 4.2 Loading progress\n",
+    "## 4.2 Loading the optimizer state\n",
     "\n",
-    "Naturally, if you stored progress you will be able to load that onto a new instance of `BayesianOptimization`. The easiest way to do it is by invoking the `load_logs` function, from the `util` submodule."
+    "To load with a previously saved state, pass the path of your saved state file to the `load_state_path` parameter. Note that if you've changed the bounds of your parameters, you'll need to pass the updated bounds to the new optimizer.\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [],
-   "source": [
-    "from bayes_opt.util import load_logs"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 18,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "0\n"
-     ]
-    }
-   ],
    "source": [
     "new_optimizer = BayesianOptimization(\n",
     "    f=black_box_function,\n",
-    "    pbounds={\"x\": (-2, 2), \"y\": (-2, 2)},\n",
-    "    verbose=2,\n",
-    "    random_state=7,\n",
+    "    pbounds={\"x\": (-2, 3), \"y\": (-3, 3)},\n",
+    "    random_state=1,\n",
+    "    verbose=0\n",
     ")\n",
-    "print(len(new_optimizer.space))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 19,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "load_logs(new_optimizer, logs=[\"./logs.log\"]);"
+    "\n",
+    "new_optimizer.load_state(\"./optimizer_state.json\")\n",
+    "\n",
+    "# Continue optimization\n",
+    "new_optimizer.maximize(\n",
+    "    init_points=0,\n",
+    "    n_iter=5\n",
+    ")"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": 20,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "New optimizer is now aware of 5 points.\n"
-     ]
-    }
-   ],
    "source": [
-    "print(\"New optimizer is now aware of {} points.\".format(len(new_optimizer.space)))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 21,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "|   iter    |  target   |     x     |     y     |\n",
-      "-------------------------------------------------\n",
-      "| \u001b[0m1        \u001b[0m | \u001b[0m-0.6164  \u001b[0m | \u001b[0m-1.271   \u001b[0m | \u001b[0m1.045    \u001b[0m |\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "| \u001b[95m2        \u001b[0m | \u001b[95m0.9988   \u001b[0m | \u001b[95m0.005532 \u001b[0m | \u001b[95m1.034    \u001b[0m |\n",
-      "| \u001b[0m3        \u001b[0m | \u001b[0m0.9757   \u001b[0m | \u001b[0m0.1501   \u001b[0m | \u001b[0m1.043    \u001b[0m |\n",
-      "| \u001b[0m4        \u001b[0m | \u001b[0m0.9819   \u001b[0m | \u001b[0m-0.07993 \u001b[0m | \u001b[0m0.8917   \u001b[0m |\n",
-      "| \u001b[0m5        \u001b[0m | \u001b[0m-3.35    \u001b[0m | \u001b[0m1.976    \u001b[0m | \u001b[0m0.3338   \u001b[0m |\n",
-      "| \u001b[0m6        \u001b[0m | \u001b[0m0.9778   \u001b[0m | \u001b[0m-0.1166  \u001b[0m | \u001b[0m1.093    \u001b[0m |\n",
-      "| \u001b[0m7        \u001b[0m | \u001b[0m0.6605   \u001b[0m | \u001b[0m-0.02272 \u001b[0m | \u001b[0m0.4178   \u001b[0m |\n",
-      "| \u001b[0m8        \u001b[0m | \u001b[0m-0.5044  \u001b[0m | \u001b[0m-0.7146  \u001b[0m | \u001b[0m1.997    \u001b[0m |\n",
-      "| \u001b[95m9        \u001b[0m | \u001b[95m0.9997   \u001b[0m | \u001b[95m0.0084   \u001b[0m | \u001b[95m1.016    \u001b[0m |\n",
-      "| \u001b[0m10       \u001b[0m | \u001b[0m-11.83   \u001b[0m | \u001b[0m1.995    \u001b[0m | \u001b[0m-1.974   \u001b[0m |\n",
-      "=================================================\n"
-     ]
-    }
-   ],
-   "source": [
-    "new_optimizer.maximize(\n",
-    "    init_points=0,\n",
-    "    n_iter=10,\n",
-    ")"
+    "This provides a simpler alternative to the logging system shown in section 4, especially when you want to continue optimization from a previous state."
    ]
   },
   {
@@ -489,6 +381,11 @@
     "\n",
     "This tour should be enough to cover most usage scenarios of this package. If, however, you feel like you need to know more, please checkout the `advanced-tour` notebook. There you will be able to find other, more advanced features of this package that could be what you're looking for. Also, browse the examples folder for implementation tips and ideas."
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
   }
  ],
  "metadata": {
@@ -507,7 +404,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.0"
+   "version": "3.13.1"
   }
  },
  "nbformat": 4,
diff --git a/pyproject.toml b/pyproject.toml
index 884c3ac48..a1ff1d7be 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,24 +1,36 @@
-[tool.poetry]
+[project]
 name = "bayesian-optimization"
 version = "3.0.0b1"
 description = "Bayesian Optimization package"
-authors = ["Fernando Nogueira"]
-license = "MIT"
+authors = [{ name = "Fernando Nogueira", email = "fmfnogueira@gmail.com" }]
+license = { file = "LICENSE" }
 readme = "README.md"
-packages = [{include = "bayes_opt"}]
-
-[tool.poetry.dependencies]
-python = "^3.9"
-scikit-learn = "^1.0.0"
-numpy = ">=1.25"
-scipy =  [
-    {version = "^1.0.0", python = "<3.13"},
-    {version = "^1.14.1", python = ">=3.13"}
+requires-python = ">=3.9,<4.0"
+classifiers = [
+    "License :: OSI Approved :: MIT License",
+    "Programming Language :: Python",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.9",
+    "Programming Language :: Python :: 3.10",
+    "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
+]
+dependencies = [
+    "scikit-learn>=1.0.0,<2.0.0",
+    "numpy>=1.25; python_version<'3.13'",
+    "numpy>=2.1.3; python_version>='3.13'",
+    "scipy>=1.0.0,<2.0.0; python_version<'3.13'",
+    "scipy>=1.14.1,<2.0.0; python_version>='3.13'",
+    "colorama>=0.4.6,<1.0.0",
 ]
-colorama = "^0.4.6"
+
+[tool.poetry]
+requires-poetry = ">=2.0"
+packages = [{ include = "bayes_opt" }]
 
 
-[tool.poetry.group.dev]  # for testing/developing
+[tool.poetry.group.dev] # for testing/developing
 optional = true
 [tool.poetry.group.dev.dependencies]
 pytest = "^8.0.0"
@@ -28,7 +40,7 @@ ruff = "0.6.6"
 pre-commit = "^3.7.1"
 
 
-[tool.poetry.group.nbtools]  # for running/converting notebooks
+[tool.poetry.group.nbtools] # for running/converting notebooks
 optional = true
 [tool.poetry.group.nbtools.dependencies]
 nbformat = "^5.9.2"
@@ -38,17 +50,17 @@ matplotlib = "^3.0"
 nbsphinx = "^0.9.4"
 sphinx-immaterial = "^0.12.0"
 sphinx = [
-    {version = "^7.0.0", python = "<3.10"},
-    {version = "^8.0.0", python = ">=3.10"}
+    { version = "^7.0.0", python = "<3.10" },
+    { version = "^8.0.0", python = ">=3.10" },
 ]
 sphinx-autodoc-typehints = [
-    {version = "^2.3.0", python = "<3.10"},
-    {version = "^2.4.0", python = ">=3.10"}
+    { version = "^2.3.0", python = "<3.10" },
+    { version = "^2.4.0", python = ">=3.10" },
 ]
 
 
 [build-system]
-requires = ["poetry-core"]
+requires = ["poetry-core>=2.0"]
 build-backend = "poetry.core.masonry.api"
 
 [tool.coverage.report]
diff --git a/tests/test_acquisition.py b/tests/test_acquisition.py
index 1191976df..e3d3e51ca 100644
--- a/tests/test_acquisition.py
+++ b/tests/test_acquisition.py
@@ -4,14 +4,38 @@
 
 import numpy as np
 import pytest
+from scipy.optimize import NonlinearConstraint
 from scipy.spatial.distance import pdist
 from sklearn.gaussian_process import GaussianProcessRegressor
 
-from bayes_opt import acquisition, exception
+from bayes_opt import BayesianOptimization, acquisition, exception
+from bayes_opt.acquisition import (
+    ConstantLiar,
+    ExpectedImprovement,
+    GPHedge,
+    ProbabilityOfImprovement,
+    UpperConfidenceBound,
+)
 from bayes_opt.constraint import ConstraintModel
 from bayes_opt.target_space import TargetSpace
 
 
+# Test fixtures
+@pytest.fixture
+def target_func_x_and_y():
+    return lambda x, y: -((x - 1) ** 2) - (y - 2) ** 2
+
+
+@pytest.fixture
+def pbounds():
+    return {"x": (-5, 5), "y": (-5, 5)}
+
+
+@pytest.fixture
+def constraint(constraint_func):
+    return NonlinearConstraint(fun=constraint_func, lb=-1.0, ub=4.0)
+
+
 @pytest.fixture
 def target_func():
     return lambda x: sum(x)
@@ -32,6 +56,11 @@ def target_space(target_func):
     return TargetSpace(target_func=target_func, pbounds={"x": (1, 4), "y": (0, 3.0)})
 
 
+@pytest.fixture
+def constraint_func():
+    return lambda x, y: x + y
+
+
 @pytest.fixture
 def constrained_target_space(target_func):
     constraint_model = ConstraintModel(fun=lambda params: params["x"] + params["y"], lb=0.0, ub=1.0)
@@ -53,6 +82,12 @@ def mock_acq(x: np.ndarray):
     def base_acq(self, mean, std):
         pass
 
+    def get_acquisition_params(self) -> dict:
+        return {}
+
+    def set_acquisition_params(self, params: dict) -> None:
+        pass
+
 
 def test_base_acquisition():
     acq = acquisition.UpperConfidenceBound()
@@ -351,3 +386,255 @@ def test_gphedge_integration(gp, target_space, random_state):
 def test_upper_confidence_bound_invalid_kappa_error(kappa: float):
     with pytest.raises(ValueError, match="kappa must be greater than or equal to 0."):
         acquisition.UpperConfidenceBound(kappa=kappa)
+
+
+def verify_optimizers_match(optimizer1, optimizer2):
+    """Helper function to verify two optimizers match."""
+    assert len(optimizer1.space) == len(optimizer2.space)
+    assert optimizer1.max["target"] == optimizer2.max["target"]
+    assert optimizer1.max["params"] == optimizer2.max["params"]
+
+    np.testing.assert_array_equal(optimizer1.space.params, optimizer2.space.params)
+    np.testing.assert_array_equal(optimizer1.space.target, optimizer2.space.target)
+
+    if optimizer1.is_constrained:
+        np.testing.assert_array_equal(
+            optimizer1.space._constraint_values, optimizer2.space._constraint_values
+        )
+        assert optimizer1.space._constraint.lb == optimizer2.space._constraint.lb
+        assert optimizer1.space._constraint.ub == optimizer2.space._constraint.ub
+
+    rng = np.random.default_rng()
+    assert rng.bit_generator.state["state"]["state"] == rng.bit_generator.state["state"]["state"]
+
+    assert optimizer1._gp.kernel.get_params() == optimizer2._gp.kernel.get_params()
+
+    suggestion1 = optimizer1.suggest()
+    suggestion2 = optimizer2.suggest()
+    assert suggestion1 == suggestion2, f"\nSuggestion 1: {suggestion1}\nSuggestion 2: {suggestion2}"
+
+
+def test_integration_upper_confidence_bound(target_func_x_and_y, pbounds, tmp_path):
+    """Test save/load integration with UpperConfidenceBound acquisition."""
+    acquisition_function = UpperConfidenceBound(kappa=2.576)
+
+    # Create and run first optimizer
+    optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=acquisition_function,
+        random_state=1,
+        verbose=0,
+    )
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    # Save state
+    state_path = tmp_path / "ucb_state.json"
+    optimizer.save_state(state_path)
+
+    # Create new optimizer and load state
+    new_optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=UpperConfidenceBound(kappa=2.576),
+        random_state=1,
+        verbose=0,
+    )
+    new_optimizer.load_state(state_path)
+
+    verify_optimizers_match(optimizer, new_optimizer)
+
+
+def test_integration_probability_improvement(target_func_x_and_y, pbounds, tmp_path):
+    """Test save/load integration with ProbabilityOfImprovement acquisition."""
+    acquisition_function = ProbabilityOfImprovement(xi=0.01)
+
+    optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=acquisition_function,
+        random_state=1,
+        verbose=0,
+    )
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    state_path = tmp_path / "pi_state.json"
+    optimizer.save_state(state_path)
+
+    new_optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=ProbabilityOfImprovement(xi=0.01),
+        random_state=1,
+        verbose=0,
+    )
+    new_optimizer.load_state(state_path)
+
+    verify_optimizers_match(optimizer, new_optimizer)
+
+
+def test_integration_expected_improvement(target_func_x_and_y, pbounds, tmp_path):
+    """Test save/load integration with ExpectedImprovement acquisition."""
+    acquisition_function = ExpectedImprovement(xi=0.01)
+
+    optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=acquisition_function,
+        random_state=1,
+        verbose=0,
+    )
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    state_path = tmp_path / "ei_state.json"
+    optimizer.save_state(state_path)
+
+    new_optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=ExpectedImprovement(xi=0.01),
+        random_state=1,
+        verbose=0,
+    )
+    new_optimizer.load_state(state_path)
+
+    verify_optimizers_match(optimizer, new_optimizer)
+
+
+def test_integration_constant_liar(target_func_x_and_y, pbounds, tmp_path):
+    """Test save/load integration with ConstantLiar acquisition."""
+    base_acq = UpperConfidenceBound(kappa=2.576)
+    acquisition_function = ConstantLiar(base_acquisition=base_acq)
+
+    optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=acquisition_function,
+        random_state=1,
+        verbose=0,
+    )
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    state_path = tmp_path / "cl_state.json"
+    optimizer.save_state(state_path)
+
+    new_optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=ConstantLiar(base_acquisition=UpperConfidenceBound(kappa=2.576)),
+        random_state=1,
+        verbose=0,
+    )
+    new_optimizer.load_state(state_path)
+
+    verify_optimizers_match(optimizer, new_optimizer)
+
+
+def test_integration_gp_hedge(target_func_x_and_y, pbounds, tmp_path):
+    """Test save/load integration with GPHedge acquisition."""
+    base_acquisitions = [
+        UpperConfidenceBound(kappa=2.576),
+        ProbabilityOfImprovement(xi=0.01),
+        ExpectedImprovement(xi=0.01),
+    ]
+    acquisition_function = GPHedge(base_acquisitions=base_acquisitions)
+
+    optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=acquisition_function,
+        random_state=1,
+        verbose=0,
+    )
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    state_path = tmp_path / "gphedge_state.json"
+    optimizer.save_state(state_path)
+
+    new_base_acquisitions = [
+        UpperConfidenceBound(kappa=2.576),
+        ProbabilityOfImprovement(xi=0.01),
+        ExpectedImprovement(xi=0.01),
+    ]
+    new_optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=GPHedge(base_acquisitions=new_base_acquisitions),
+        random_state=1,
+        verbose=0,
+    )
+    new_optimizer.load_state(state_path)
+
+    verify_optimizers_match(optimizer, new_optimizer)
+
+
+def test_integration_constrained(target_func_x_and_y, pbounds, constraint, tmp_path):
+    """Test save/load integration with constraints."""
+    acquisition_function = ExpectedImprovement(xi=0.01)
+
+    optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=acquisition_function,
+        constraint=constraint,
+        random_state=1,
+        verbose=0,
+    )
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    state_path = tmp_path / "constrained_state.json"
+    optimizer.save_state(state_path)
+
+    new_optimizer = BayesianOptimization(
+        f=target_func_x_and_y,
+        pbounds=pbounds,
+        acquisition_function=ExpectedImprovement(xi=0.01),
+        constraint=constraint,
+        random_state=1,
+        verbose=0,
+    )
+    new_optimizer.load_state(state_path)
+
+    verify_optimizers_match(optimizer, new_optimizer)
+
+
+def test_custom_acquisition_without_get_params():
+    """Test that a custom acquisition function without get_acquisition_params raises NotImplementedError."""
+
+    class CustomAcqWithoutGetParams(acquisition.AcquisitionFunction):
+        def __init__(self, random_state=None):
+            super().__init__(random_state=random_state)
+
+        def base_acq(self, mean, std):
+            return mean + std
+
+        def set_acquisition_params(self, params):
+            pass
+
+    acq = CustomAcqWithoutGetParams()
+    with pytest.raises(
+        NotImplementedError,
+        match="Custom AcquisitionFunction subclasses must implement their own get_acquisition_params method",
+    ):
+        acq.get_acquisition_params()
+
+
+def test_custom_acquisition_without_set_params():
+    """Test that a custom acquisition function without set_acquisition_params raises NotImplementedError."""
+
+    class CustomAcqWithoutSetParams(acquisition.AcquisitionFunction):
+        def __init__(self, random_state=None):
+            super().__init__(random_state=random_state)
+
+        def base_acq(self, mean, std):
+            return mean + std
+
+        def get_acquisition_params(self):
+            return {}
+
+    acq = CustomAcqWithoutSetParams()
+    with pytest.raises(
+        NotImplementedError,
+        match="Custom AcquisitionFunction subclasses must implement their own set_acquisition_params method",
+    ):
+        acq.set_acquisition_params(params={})
diff --git a/tests/test_bayesian_optimization.py b/tests/test_bayesian_optimization.py
index 48e1af115..236de2de7 100644
--- a/tests/test_bayesian_optimization.py
+++ b/tests/test_bayesian_optimization.py
@@ -5,13 +5,17 @@
 
 import numpy as np
 import pytest
+from scipy.optimize import NonlinearConstraint
 
 from bayes_opt import BayesianOptimization, acquisition
 from bayes_opt.acquisition import AcquisitionFunction
+from bayes_opt.domain_reduction import SequentialDomainReductionTransformer
 from bayes_opt.event import DEFAULT_EVENTS, Events
 from bayes_opt.exception import NotUniqueError
 from bayes_opt.logger import ScreenLogger
+from bayes_opt.parameter import BayesParameter
 from bayes_opt.target_space import TargetSpace
+from bayes_opt.util import ensure_rng
 
 
 def target_func(**kwargs):
@@ -333,3 +337,312 @@ def test_duplicate_points():
     optimizer.register(params=next_point_to_probe, target=target)
     # and again (should throw warning)
     optimizer.register(params=next_point_to_probe, target=target)
+
+
+def test_save_load_state(tmp_path):
+    """Test saving and loading optimizer state."""
+    # Initialize and run original optimizer
+    optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    # Save state
+    state_path = tmp_path / "optimizer_state.json"
+    optimizer.save_state(state_path)
+
+    # Create new optimizer and load state
+    new_optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+    new_optimizer.load_state(state_path)
+
+    # Test that key properties match
+    assert len(optimizer.space) == len(new_optimizer.space)
+    assert optimizer.max["target"] == new_optimizer.max["target"]
+    assert optimizer.max["params"] == new_optimizer.max["params"]
+    np.testing.assert_array_equal(optimizer.space.params, new_optimizer.space.params)
+    np.testing.assert_array_equal(optimizer.space.target, new_optimizer.space.target)
+
+
+def test_save_load_w_categorical_params(tmp_path):
+    """Test saving and loading optimizer state with categorical parameters."""
+
+    def str_target_func(param1: str, param2: str) -> float:
+        # Simple function that maps strings to numbers
+        value_map = {"low": 1.0, "medium": 2.0, "high": 3.0}
+        return value_map[param1] + value_map[param2]
+
+    str_pbounds = {"param1": ["low", "medium", "high"], "param2": ["low", "medium", "high"]}
+
+    optimizer = BayesianOptimization(f=str_target_func, pbounds=str_pbounds, random_state=1, verbose=0)
+
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    state_path = tmp_path / "optimizer_state.json"
+    optimizer.save_state(state_path)
+
+    new_optimizer = BayesianOptimization(f=str_target_func, pbounds=str_pbounds, random_state=1, verbose=0)
+    new_optimizer.load_state(state_path)
+
+    assert len(optimizer.space) == len(new_optimizer.space)
+    assert optimizer.max["target"] == new_optimizer.max["target"]
+    assert optimizer.max["params"] == new_optimizer.max["params"]
+    for i in range(len(optimizer.space)):
+        assert isinstance(optimizer.res[i]["params"]["param1"], str)
+        assert isinstance(optimizer.res[i]["params"]["param2"], str)
+        assert isinstance(new_optimizer.res[i]["params"]["param1"], str)
+        assert isinstance(new_optimizer.res[i]["params"]["param2"], str)
+        assert optimizer.res[i]["params"] == new_optimizer.res[i]["params"]
+
+
+def test_suggest_point_returns_same_point(tmp_path):
+    """Check that suggest returns same point after save/load."""
+    optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    state_path = tmp_path / "optimizer_state.json"
+    optimizer.save_state(state_path)
+
+    new_optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+    new_optimizer.load_state(state_path)
+
+    # Both optimizers should suggest the same point
+    suggestion1 = optimizer.suggest()
+    suggestion2 = new_optimizer.suggest()
+    assert suggestion1 == suggestion2
+
+
+def test_save_load_random_state(tmp_path):
+    """Test that random state is properly preserved."""
+    # Initialize optimizer
+    optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+
+    # Register a point before saving
+    optimizer.probe(params={"p1": 1, "p2": 2}, lazy=False)
+
+    # Save state
+    state_path = tmp_path / "optimizer_state.json"
+    optimizer.save_state(state_path)
+
+    # Create new optimizer with same configuration
+    new_optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+    new_optimizer.load_state(state_path)
+
+    # Both optimizers should suggest the same point
+    suggestion1 = optimizer.suggest()
+    suggestion2 = new_optimizer.suggest()
+    assert suggestion1 == suggestion2
+
+
+def test_save_load_unused_optimizer(tmp_path):
+    """Test saving and loading optimizer state with unused optimizer."""
+    optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+
+    # Test that saving without samples raises an error
+    with pytest.raises(ValueError, match="Cannot save optimizer state before collecting any samples"):
+        optimizer.save_state(tmp_path / "optimizer_state.json")
+
+    # Add a sample point
+    optimizer.probe(params={"p1": 1, "p2": 2}, lazy=False)
+
+    # Now saving should work
+    optimizer.save_state(tmp_path / "optimizer_state.json")
+
+    new_optimizer = BayesianOptimization(f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0)
+    new_optimizer.load_state(tmp_path / "optimizer_state.json")
+
+    assert len(optimizer.space) == len(new_optimizer.space)
+    assert optimizer.max["target"] == new_optimizer.max["target"]
+    assert optimizer.max["params"] == new_optimizer.max["params"]
+    np.testing.assert_array_equal(optimizer.space.params, new_optimizer.space.params)
+    np.testing.assert_array_equal(optimizer.space.target, new_optimizer.space.target)
+
+    """Test saving and loading optimizer state with constraints."""
+
+    def constraint_func(x: float, y: float) -> float:
+        return x + y  # Simple constraint: sum of parameters should be within bounds
+
+    constraint = NonlinearConstraint(fun=constraint_func, lb=0.0, ub=3.0)
+
+    # Initialize optimizer with constraint
+    optimizer = BayesianOptimization(
+        f=target_func, pbounds={"x": (-1, 3), "y": (0, 5)}, constraint=constraint, random_state=1, verbose=0
+    )
+
+    # Register some points, some that satisfy constraint and some that don't
+    optimizer.register(
+        params={"x": 1.0, "y": 1.0},  # Satisfies constraint: sum = 2.0
+        target=2.0,
+        constraint_value=2.0,
+    )
+    optimizer.register(
+        params={"x": 2.0, "y": 2.0},  # Violates constraint: sum = 4.0
+        target=4.0,
+        constraint_value=4.0,
+    )
+    optimizer.register(
+        params={"x": 0.5, "y": 0.5},  # Satisfies constraint: sum = 1.0
+        target=1.0,
+        constraint_value=1.0,
+    )
+
+    state_path = tmp_path / "optimizer_state.json"
+    optimizer.save_state(state_path)
+
+    new_optimizer = BayesianOptimization(
+        f=target_func, pbounds={"x": (-1, 3), "y": (0, 5)}, constraint=constraint, random_state=1, verbose=0
+    )
+    new_optimizer.load_state(state_path)
+
+    # Test that key properties match
+    assert len(optimizer.space) == len(new_optimizer.space)
+    assert optimizer.max["target"] == new_optimizer.max["target"]
+    assert optimizer.max["params"] == new_optimizer.max["params"]
+    np.testing.assert_array_equal(optimizer.space.params, new_optimizer.space.params)
+    np.testing.assert_array_equal(optimizer.space.target, new_optimizer.space.target)
+
+    # Test that constraint values were properly saved and loaded
+    np.testing.assert_array_equal(optimizer.space._constraint_values, new_optimizer.space._constraint_values)
+
+    # Test that both optimizers suggest the same point (should respect constraints)
+    suggestion1 = optimizer.suggest()
+    suggestion2 = new_optimizer.suggest()
+    assert suggestion1 == suggestion2
+
+    # Verify that suggested point satisfies constraint
+    constraint_value = constraint_func(**suggestion1)
+    assert 0.0 <= constraint_value <= 3.0, "Suggested point violates constraint"
+
+
+def test_save_load_w_domain_reduction(tmp_path):
+    """Test saving and loading optimizer state with domain reduction transformer."""
+    # Initialize optimizer with bounds transformer
+    bounds_transformer = SequentialDomainReductionTransformer()
+    optimizer = BayesianOptimization(
+        f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0, bounds_transformer=bounds_transformer
+    )
+
+    # Run some iterations to trigger domain reduction
+    optimizer.maximize(init_points=2, n_iter=3)
+
+    # Save state
+    state_path = tmp_path / "optimizer_state.json"
+    optimizer.save_state(state_path)
+
+    # Create new optimizer with same configuration
+    new_bounds_transformer = SequentialDomainReductionTransformer()
+    new_optimizer = BayesianOptimization(
+        f=target_func, pbounds=PBOUNDS, random_state=1, verbose=0, bounds_transformer=new_bounds_transformer
+    )
+    new_optimizer.load_state(state_path)
+
+    # Both optimizers should probe the same point
+    point = {"p1": 1.5, "p2": 0.5}
+    probe1 = optimizer.probe(point)
+    probe2 = new_optimizer.probe(point)
+    assert probe1 == probe2
+
+    # Both optimizers should suggest the same point
+    suggestion1 = optimizer.suggest()
+    suggestion2 = new_optimizer.suggest()
+    assert suggestion1 == suggestion2
+
+    # Verify that the transformed bounds match
+    assert optimizer._space.bounds.tolist() == new_optimizer._space.bounds.tolist()
+
+
+def test_save_load_w_custom_parameter(tmp_path):
+    """Test saving and loading optimizer state with custom parameter types."""
+
+    class FixedPerimeterTriangleParameter(BayesParameter):
+        def __init__(self, name: str, bounds, perimeter) -> None:
+            super().__init__(name, bounds)
+            self.perimeter = perimeter
+
+        @property
+        def is_continuous(self):
+            return True
+
+        def random_sample(self, n_samples: int, random_state):
+            random_state = ensure_rng(random_state)
+            samples = []
+            while len(samples) < n_samples:
+                samples_ = random_state.dirichlet(np.ones(3), n_samples)
+                samples_ = samples_ * self.perimeter  # scale samples by perimeter
+
+                samples_ = samples_[
+                    np.all((self.bounds[:, 0] <= samples_) & (samples_ <= self.bounds[:, 1]), axis=-1)
+                ]
+                samples.extend(np.atleast_2d(samples_))
+            return np.array(samples[:n_samples])
+
+        def to_float(self, value):
+            return value
+
+        def to_param(self, value):
+            return value * self.perimeter / sum(value)
+
+        def kernel_transform(self, value):
+            return value * self.perimeter / np.sum(value, axis=-1, keepdims=True)
+
+        def to_string(self, value, str_len: int) -> str:
+            len_each = (str_len - 2) // 3
+            str_ = "|".join([f"{float(np.round(value[i], 4))}"[:len_each] for i in range(3)])
+            return str_.ljust(str_len)
+
+        @property
+        def dim(self):
+            return 3  # as we have three float values, each representing the length of one side.
+
+    def area_of_triangle(sides):
+        a, b, c = sides
+        s = np.sum(sides, axis=-1)  # perimeter
+        return np.sqrt(s * (s - a) * (s - b) * (s - c))
+
+    # Create parameter and bounds
+    param = FixedPerimeterTriangleParameter(
+        name="sides", bounds=np.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]]), perimeter=1.0
+    )
+    pbounds = {"sides": param}
+
+    # Print initial pbounds
+    print("\nOriginal pbounds:")
+    print(pbounds)
+
+    # Initialize first optimizer
+    optimizer = BayesianOptimization(f=area_of_triangle, pbounds=pbounds, random_state=1, verbose=0)
+
+    # Run iterations and immediately save state
+    optimizer.maximize(init_points=2, n_iter=5)
+
+    # Force GP update before saving
+    optimizer._gp.fit(optimizer.space.params, optimizer.space.target)
+
+    # Save state
+    state_path = tmp_path / "optimizer_state.json"
+    optimizer.save_state(state_path)
+
+    # Create new optimizer and load state
+    new_optimizer = BayesianOptimization(f=area_of_triangle, pbounds=pbounds, random_state=1, verbose=0)
+    new_optimizer.load_state(state_path)
+
+    # Test that key properties match
+    assert len(optimizer.space) == len(new_optimizer.space)
+    assert optimizer.max["target"] == new_optimizer.max["target"]
+    np.testing.assert_array_almost_equal(
+        optimizer.max["params"]["sides"], new_optimizer.max["params"]["sides"], decimal=10
+    )
+
+    # Test that all historical data matches
+    for i in range(len(optimizer.space)):
+        np.testing.assert_array_almost_equal(
+            optimizer.space.params[i], new_optimizer.space.params[i], decimal=10
+        )
+        assert optimizer.space.target[i] == new_optimizer.space.target[i]
+        np.testing.assert_array_almost_equal(
+            optimizer.res[i]["params"]["sides"], new_optimizer.res[i]["params"]["sides"], decimal=10
+        )
+        assert optimizer.res[i]["target"] == new_optimizer.res[i]["target"]
+
+    # Test that multiple subsequent suggestions match
+    for _ in range(5):
+        suggestion1 = optimizer.suggest()
+        suggestion2 = new_optimizer.suggest()
+        np.testing.assert_array_almost_equal(suggestion1["sides"], suggestion2["sides"], decimal=10)