diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
index 3a809789e..0551a364e 100644
--- a/.github/workflows/run_tests.yml
+++ b/.github/workflows/run_tests.yml
@@ -7,7 +7,6 @@ on:
   push:
     branches: [ "master" ]
   pull_request:
-    branches: [ "master" ]
 
 permissions:
   contents: read
diff --git a/bayes_opt/constraint.py b/bayes_opt/constraint.py
index a31142310..5ea333525 100644
--- a/bayes_opt/constraint.py
+++ b/bayes_opt/constraint.py
@@ -1,3 +1,4 @@
+"""Constraint handling."""
 import numpy as np
 from sklearn.gaussian_process.kernels import Matern
 from sklearn.gaussian_process import GaussianProcessRegressor
@@ -5,42 +6,35 @@
 
 
 class ConstraintModel():
-    """
+    """Model constraints using GPRs.
+
     This class takes the function to optimize as well as the parameters bounds
     in order to find which values for the parameters yield the maximum value
     using bayesian optimization.
 
     Parameters
     ----------
-    fun: function
-        Constraint function. If multiple constraints are handled, this should
-        return a numpy.ndarray of appropriate size.
-
-    lb: numeric or numpy.ndarray
-        Upper limit(s) for the constraints. The return value of `fun` should
-        have exactly this shape.
-
-    ub: numeric or numpy.ndarray
-        Upper limit(s) for the constraints. The return value of `fun` should
-        have exactly this shape.
-
-    random_state: int or numpy.random.RandomState, optional(default=None)
-        If the value is an integer, it is used as the seed for creating a
-        numpy.random.RandomState. Otherwise the random state provided is used.
-        When set to None, an unseeded random state is generated.
-
-    Note
-    ----
-    In case of multiple constraints, this model assumes conditional
-    independence. This means that for each constraint, the probability of
-    fulfillment is the cdf of a univariate Gaussian. The overall probability
-    is a simply the product of the individual probabilities.
+    fun : None or Callable -> float or np.ndarray
+        The constraint function. Should be float-valued or array-valued (if
+        multiple constraints are present). Needs to take the same parameters
+        as the optimization target with the same argument names.
+
+    lb : float or np.ndarray
+        The lower bound on the constraints. Should have the same
+        dimensionality as the return value of the constraint function.
+
+    ub : float or np.ndarray
+        The upper bound on the constraints. Should have the same
+        dimensionality as the return value of the constraint function.
+
+    random_state : np.random.RandomState or int or None, default=None
+        Random state to use.
     """
 
     def __init__(self, fun, lb, ub, random_state=None):
         self.fun = fun
 
-        self._lb = np.atleast_1d(lb)        
+        self._lb = np.atleast_1d(lb)
         self._ub = np.atleast_1d(ub)
 
         if np.any(self._lb >= self._ub):
@@ -58,19 +52,36 @@ def __init__(self, fun, lb, ub, random_state=None):
 
     @property
     def lb(self):
+        """Return lower bounds."""
         return self._lb
 
     @property
     def ub(self):
+        """Return upper bounds."""
         return self._ub
-    
+
     @property
     def model(self):
+        """Return GP regressors of the constraint function."""
         return self._model
 
     def eval(self, **kwargs):
-        """
-        Evaluates the constraint function.
+        """Evaluate the constraint function.
+
+        Parameters
+        ----------
+        **kwargs :
+            Function arguments to evaluate the constraint function on.
+
+
+        Returns
+        -------
+        Value of the constraint function.
+
+        Raises
+        ------
+        TypeError
+            If the kwargs keys don't match the function argument names.
         """
         try:
             return self.fun(**kwargs)
@@ -85,8 +96,19 @@ def eval(self, **kwargs):
             raise
 
     def fit(self, X, Y):
-        """
-        Fits internal GaussianProcessRegressor's to the data.
+        """Fit internal GPRs to the data.
+
+        Parameters
+        ----------
+        X :
+            Parameters of the constraint function.
+        Y :
+            Values of the constraint function.
+
+
+        Returns
+        -------
+        None
         """
         if len(self._model) == 1:
             self._model[0].fit(X, Y)
@@ -95,13 +117,39 @@ def fit(self, X, Y):
                 gp.fit(X, Y[:, i])
 
     def predict(self, X):
-        """
-        Returns the probability that the constraint is fulfilled at `X` based
-        on the internal Gaussian Process Regressors.
+        r"""Calculate the probability that the constraint is fulfilled at `X`.
+
+        Note that this does not try to approximate the values of the
+        constraint function (for this, see `ConstraintModel.approx()`.), but
+        probability that the constraint function is fulfilled. That is, this
+        function calculates
+
+        .. math::
+            p = \text{Pr}\left\{c^{\text{low}} \leq \tilde{c}(x) \leq
+                c^{\text{up}} \right\} = \int_{c^{\text{low}}}^{c^{\text{up}}}
+                \mathcal{N}(c, \mu(x), \sigma^2(x)) \, dc.
+
+        with :math:`\mu(x)`, :math:`\sigma^2(x)` the mean and variance at
+        :math:`x` as given by the GP and :math:`c^{\text{low}}`,
+        :math:`c^{\text{up}}` the lower and upper bounds of the constraint
+        respectively.
+
+        In case of multiple constraints, we assume conditional independence.
+        This means we calculate the probability of constraint fulfilment
+        individually, with the joint probability given as their product.
+
+        Parameters
+        ----------
+        X : np.ndarray of shape (n_samples, n_features)
+            Parameters for which to predict the probability of constraint
+            fulfilment.
+            
+
+        Returns
+        -------
+        np.ndarray of shape (n_samples,)
+            Probability of constraint fulfilment.
 
-        Note that this does not try to approximate the values of the constraint
-        function, but probability that the constraint function is fulfilled.
-        For the former, see `ConstraintModel.approx()`.
         """
         X_shape = X.shape
         X = X.reshape((-1, self._model[0].n_features_in_))
@@ -114,33 +162,53 @@ def predict(self, X):
                             if self._lb[0] != np.inf else np.array([1]))
             result = p_upper - p_lower
             return result.reshape(X_shape[:-1])
-        else:
-            result = np.ones(X.shape[0])
-            for j, gp in enumerate(self._model):
-                y_mean, y_std = gp.predict(X, return_std=True)
-                p_lower = (norm(loc=y_mean, scale=y_std).cdf(self._lb[j])
-                           if self._lb[j] != -np.inf else np.array([0]))
-                p_upper = (norm(loc=y_mean, scale=y_std).cdf(self._ub[j])
-                           if self._lb[j] != np.inf else np.array([1]))
-                result = result * (p_upper - p_lower)
-            return result.reshape(X_shape[:-1])
+
+        result = np.ones(X.shape[0])
+        for j, gp in enumerate(self._model):
+            y_mean, y_std = gp.predict(X, return_std=True)
+            p_lower = (norm(loc=y_mean, scale=y_std).cdf(self._lb[j])
+                        if self._lb[j] != -np.inf else np.array([0]))
+            p_upper = (norm(loc=y_mean, scale=y_std).cdf(self._ub[j])
+                        if self._lb[j] != np.inf else np.array([1]))
+            result = result * (p_upper - p_lower)
+        return result.reshape(X_shape[:-1])
 
     def approx(self, X):
         """
-        Returns the approximation of the constraint function using the internal
-        Gaussian Process Regressors.
+        Approximate the constraint function using the internal GPR model.
+
+        Parameters
+        ----------
+        X : np.ndarray of shape (n_samples, n_features)
+            Parameters for which to estimate the constraint function value.
+
+        Returns
+        -------
+        np.ndarray of shape (n_samples, n_constraints)
+            Constraint function value estimates.
         """
         X_shape = X.shape
         X = X.reshape((-1, self._model[0].n_features_in_))
         if len(self._model) == 1:
             return self._model[0].predict(X).reshape(X_shape[:-1])
-        else:
-            result = np.column_stack([gp.predict(X) for gp in self._model])
-            return result.reshape(X_shape[:-1] + (len(self._lb), ))
+
+        result = np.column_stack([gp.predict(X) for gp in self._model])
+        return result.reshape(X_shape[:-1] + (len(self._lb), ))
 
     def allowed(self, constraint_values):
-        """
-        Checks whether `constraint_values` are below the specified limits.
+        """Check whether `constraint_values` fulfills the specified limits.
+
+        Parameters
+        ----------
+        constraint_values : np.ndarray of shape (n_samples, n_constraints)
+            The values of the constraint function.
+            
+
+        Returns
+        -------
+            np.ndarrray of shape (n_samples,)
+            Specifying wheter the constraints are fulfilled.
+
         """
         if self._lb.size == 1:
             return (np.less_equal(self._lb, constraint_values)
diff --git a/bayes_opt/target_space.py b/bayes_opt/target_space.py
index dea99f976..195474cb2 100644
--- a/bayes_opt/target_space.py
+++ b/bayes_opt/target_space.py
@@ -1,20 +1,39 @@
+"""Manages the optimization domain and holds points."""
 import numpy as np
-from .util import ensure_rng, NotUniqueError
 from colorama import Fore
+from .util import ensure_rng, NotUniqueError
 
 
 def _hashable(x):
-    """ ensure that an point is hashable by a python dict """
+    """Ensure that a point is hashable by a python dict."""
     return tuple(map(float, x))
 
 
-class TargetSpace(object):
-    """
-    Holds the param-space coordinates (X) and target values (Y)
-    Allows for constant-time appends while ensuring no duplicates are added
+class TargetSpace():
+    """Holds the param-space coordinates (X) and target values (Y).
+
+    Allows for constant-time appends.
+
+    Parameters
+    ----------
+    target_func : function
+        Function to be maximized.
+
+    pbounds : dict
+        Dictionary with parameters names as keys and a tuple with minimum
+        and maximum values.
+
+    random_state : int, RandomState, or None
+        optionally specify a seed for a random number generator
 
-    Example
-    -------
+    allow_duplicate_points: bool, optional (default=False)
+        If True, the optimizer will allow duplicate points to be registered.
+        This behavior may be desired in high noise situations where repeatedly probing
+        the same point will give different answers. In other situations, the acquisition
+        may occasionally generate a duplicate point.
+
+    Examples
+    --------
     >>> def target_func(p1, p2):
     >>>     return p1 + p2
     >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
@@ -26,25 +45,6 @@ class TargetSpace(object):
 
     def __init__(self, target_func, pbounds, constraint=None, random_state=None,
                  allow_duplicate_points=False):
-        """
-        Parameters
-        ----------
-        target_func : function
-            Function to be maximized.
-
-        pbounds : dict
-            Dictionary with parameters names as keys and a tuple with minimum
-            and maximum values.
-
-        random_state : int, RandomState, or None
-            optionally specify a seed for a random number generator
-
-        allow_duplicate_points: bool, optional (default=False)
-            If True, the optimizer will allow duplicate points to be registered.
-            This behavior may be desired in high noise situations where repeatedly probing
-            the same point will give different answers. In other situations, the acquisition
-            may occasionally generate a duplicate point.
-        """
         self.random_state = ensure_rng(random_state)
         self._allow_duplicate_points = allow_duplicate_points
         self.n_duplicate_points = 0
@@ -62,7 +62,7 @@ def __init__(self, target_func, pbounds, constraint=None, random_state=None,
 
         # preallocated memory for X and Y points
         self._params = np.empty(shape=(0, self.dim))
-        self._target = np.empty(shape=(0))
+        self._target = np.empty(shape=(0,))
 
         # keep track of unique points we have seen so far
         self._cache = {}
@@ -77,63 +77,149 @@ def __init__(self, target_func, pbounds, constraint=None, random_state=None,
                 self._constraint_values = np.empty(shape=(0, constraint.lb.size), dtype=float)
 
     def __contains__(self, x):
+        """Check if this parameter has already been registered.
+        
+        Returns
+        -------
+        bool
+        """
         return _hashable(x) in self._cache
 
     def __len__(self):
+        """Return number of observations registered.
+        
+        Returns
+        -------
+        int
+        """
         assert len(self._params) == len(self._target)
         return len(self._target)
 
     @property
     def empty(self):
+        """Check if anything has been registered.
+        
+        Returns
+        -------
+        bool
+        """
         return len(self) == 0
 
     @property
     def params(self):
+        """Get the parameter values registered to this TargetSpace.
+        
+        Returns
+        -------
+        np.ndarray
+        """
         return self._params
 
     @property
     def target(self):
+        """Get the target function values registered to this TargetSpace.
+        
+        Returns
+        -------
+        np.ndarray
+        """
         return self._target
 
     @property
     def dim(self):
+        """Get the number of parameter names.
+        
+        Returns
+        -------
+        int
+        """
         return len(self._keys)
 
     @property
     def keys(self):
+        """Get the keys (or parameter names).
+        
+        Returns
+        -------
+        list of str
+        """
         return self._keys
 
     @property
     def bounds(self):
+        """Get the bounds of this TargetSpace.
+        
+        Returns
+        -------
+        np.ndarray
+        """
         return self._bounds
 
     @property
     def constraint(self):
+        """Get the constraint model.
+        
+        Returns
+        -------
+        ConstraintModel
+        """
         return self._constraint
 
     @property
     def constraint_values(self):
-        if self._constraint is not None:
-            return self._constraint_values
+        """Get the constraint values registered to this TargetSpace.
+        
+        Returns
+        -------
+        np.ndarray
+        """
+        if self._constraint is None:
+            raise AttributeError("TargetSpace belongs to an unconstrained optimization")
+
+        return self._constraint_values
 
     def params_to_array(self, params):
+        """Convert a dict representation of parameters into an array version.
+
+        Parameters
+        ----------
+        params : dict
+            a single point, with len(x) == self.dim.
+
+        Returns
+        -------
+        np.ndarray
+            Representation of the parameters as dictionary.
+        """
         try:
             assert set(params) == set(self.keys)
-        except AssertionError:
+        except AssertionError as e:
             raise ValueError(
-                "Parameters' keys ({}) do ".format(sorted(params)) +
-                "not match the expected set of keys ({}).".format(self.keys)
-            )
+                f"Parameters' keys ({sorted(params)}) do " +
+                f"not match the expected set of keys ({self.keys})."
+            ) from e
         return np.asarray([params[key] for key in self.keys])
 
     def array_to_params(self, x):
+        """Convert an array representation of parameters into a dict version.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            a single point, with len(x) == self.dim.
+
+        Returns
+        -------
+        dict
+            Representation of the parameters as dictionary.
+        """
         try:
             assert len(x) == len(self.keys)
-        except AssertionError:
+        except AssertionError as e:
             raise ValueError(
-                "Size of array ({}) is different than the ".format(len(x)) +
-                "expected number of parameters ({}).".format(len(self.keys))
-            )
+                f"Size of array ({len(x)}) is different than the " +
+                f"expected number of parameters ({len(self.keys)})."
+            ) from e
         return dict(zip(self.keys, x))
 
     def _as_array(self, x):
@@ -145,20 +231,20 @@ def _as_array(self, x):
         x = x.ravel()
         try:
             assert x.size == self.dim
-        except AssertionError:
+        except AssertionError as e:
             raise ValueError(
-                "Size of array ({}) is different than the ".format(len(x)) +
-                "expected number of parameters ({}).".format(len(self.keys)))
+                f"Size of array ({len(x)}) is different than the " +
+                f"expected number of parameters ({len(self.keys)})."
+            ) from e
         return x
 
     def register(self, params, target, constraint_value=None):
-        """
-        Append a point and its target value to the known data.
+        """Append a point and its target value to the known data.
 
         Parameters
         ----------
-        x : ndarray
-            a single point, with len(x) == self.dim
+        x : np.ndarray
+            a single point, with len(x) == self.dim.
 
         y : float
             target function value
@@ -172,8 +258,8 @@ def register(self, params, target, constraint_value=None):
         -----
         runs in amortized constant time
 
-        Example
-        -------
+        Examples
+        --------
         >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
         >>> space = TargetSpace(lambda p1, p2: p1 + p2, pbounds)
         >>> len(space)
@@ -188,11 +274,12 @@ def register(self, params, target, constraint_value=None):
         if x in self:
             if self._allow_duplicate_points:
                 self.n_duplicate_points = self.n_duplicate_points + 1
-                print(Fore.RED + f'Data point {x} is not unique. {self.n_duplicate_points} duplicates registered.'
-                              f' Continuing ...')
+
+                print(Fore.RED + f'Data point {x} is not unique. {self.n_duplicate_points}'
+                              ' duplicates registered. Continuing ...')
             else:
-                raise NotUniqueError(f'Data point {x} is not unique. You can set "allow_duplicate_points=True" to '
-                                     f'avoid this error')
+                raise NotUniqueError(f'Data point {x} is not unique. You can set'
+                                     ' "allow_duplicate_points=True" to avoid this error')
 
         self._params = np.concatenate([self._params, x.reshape(1, -1)])
         self._target = np.concatenate([self._target, [target]])
@@ -211,17 +298,15 @@ def register(self, params, target, constraint_value=None):
                                                       [constraint_value]])
 
     def probe(self, params):
-        """
-        Evaluates a single point x, to obtain the value y and then records them
-        as observations.
+        """Evaluate the target function on a point and register the result.
 
         Notes
         -----
-        If x has been previously seen returns a cached value of y.
+        If `params` has been previously seen returns a cached value of `y`.
 
         Parameters
         ----------
-        x : ndarray
+        params : np.ndarray
             a single point, with len(x) == self.dim
 
         Returns
@@ -236,22 +321,21 @@ def probe(self, params):
         if self._constraint is None:
             self.register(x, target)
             return target
-        else:
-            constraint_value = self._constraint.eval(**params)
-            self.register(x, target, constraint_value)
-            return target, constraint_value
+
+        constraint_value = self._constraint.eval(**params)
+        self.register(x, target, constraint_value)
+        return target, constraint_value
 
     def random_sample(self):
-        """
-        Creates random points within the bounds of the space.
+        """Create random points within the bounds of the space.
 
         Returns
-        ----------
-        data: ndarray
+        -------
+        data: np.ndarray
             [num x dim] array points with dimensions corresponding to `self._keys`
 
-        Example
-        -------
+        Examples
+        --------
         >>> target_func = lambda p1, p2: p1 + p2
         >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
         >>> space = TargetSpace(target_func, pbounds, random_state=0)
@@ -267,7 +351,15 @@ def _target_max(self):
         """Get maximum target value found.
         
         If there is a constraint present, the maximum value that fulfills the
-        constraint is returned."""
+        constraint is returned.
+
+        
+        Returns
+        -------
+        dict | None
+            The maximum allowed point's target function value.
+            Returns None if there is no (allowed) maximum.
+        """
         if len(self.target) == 0:
             return None
 
@@ -284,7 +376,15 @@ def max(self):
         """Get maximum target value found and corresponding parameters.
         
         If there is a constraint present, the maximum value that fulfills the
-        constraint is returned."""
+        constraint is returned.
+
+        Returns
+        -------
+        dict | None
+            A dictionary containing the maximum allowed point's target function
+            value, parameters, and, if applicable, constraint function value.
+            Returns None if there is no (allowed) maximum.
+        """
         target_max = self._target_max()
 
         if target_max is None:
@@ -299,10 +399,9 @@ def max(self):
         else:
             target = self.target
             params = self.params
-            constraint_values = self.constraint_values
-        
+
         target_max_idx = np.where(target == target_max)[0][0]
-        
+
 
         res = {
                 'target': target_max,
@@ -318,6 +417,12 @@ def max(self):
 
     def res(self):
         """Get all target values and constraint fulfillment for all parameters.
+
+        Returns
+        -------
+        dict
+            A dictionary containing the target function values, parameters,
+            and, if applicable, constraint function values and allowedness.
         """
         if self._constraint is None:
             params = [dict(zip(self.keys, p)) for p in self.params]
@@ -326,27 +431,26 @@ def res(self):
                 {"target": target, "params": param}
                 for target, param in zip(self.target, params)
             ]
-        else:
-            params = [dict(zip(self.keys, p)) for p in self.params]
 
-            return [
-                {
-                    "target": target,
-                    "constraint": constraint_value,
-                    "params": param,
-                    "allowed": allowed
-                }
-                for target, constraint_value, param, allowed in zip(
-                    self.target,
-                    self._constraint_values,
-                    params,
-                    self._constraint.allowed(self._constraint_values)
-                )
-            ]
+        params = [dict(zip(self.keys, p)) for p in self.params]
+
+        return [
+            {
+                "target": target,
+                "constraint": constraint_value,
+                "params": param,
+                "allowed": allowed
+            }
+            for target, constraint_value, param, allowed in zip(
+                self.target,
+                self._constraint_values,
+                params,
+                self._constraint.allowed(self._constraint_values)
+            )
+        ]
 
     def set_bounds(self, new_bounds):
-        """
-        A method that allows changing the lower and upper searching bounds
+        """Change the lower and upper search bounds.
 
         Parameters
         ----------