diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml
index 07282413b..a57745b01 100644
--- a/.github/workflows/build_docs.yml
+++ b/.github/workflows/build_docs.yml
@@ -29,6 +29,8 @@ jobs:
           pip install nbsphinx
           pip install sphinx_rtd_theme
           pip install jupyter
+          pip install numpydoc
+          pip install myst-parser
       - name: Install package
         run: |
           pip install -e .
diff --git a/.gitignore b/.gitignore
index 3f4f6d54e..6556c710e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,4 +32,6 @@ venv.bak/
 
 docs/*
 docsrc/.ipynb_checkpoints/*
-docsrc/*.ipynb
\ No newline at end of file
+docsrc/*.ipynb
+docsrc/static/*
+docsrc/README.md
diff --git a/README.md b/README.md
index ac94cc649..652dc9a18 100644
--- a/README.md
+++ b/README.md
@@ -1,16 +1,18 @@
 <div align="center">
-  <img src="https://github.com/fmfn/BayesianOptimization/blob/master/examples/func.png"><br><br>
+  <img src="https://raw.githubusercontent.com/bayesian-optimization/BayesianOptimization/master/static/func.png"><br><br>
 </div>
 
 # Bayesian Optimization
 
-![tests](https://github.com/fmfn/BayesianOptimization/actions/workflows/run_tests.yml/badge.svg)
-[![Codecov](https://codecov.io/github/fmfn/BayesianOptimization/badge.svg?branch=master&service=github)](https://codecov.io/github/fmfn/BayesianOptimization?branch=master)
+![tests](https://github.com/bayesian-optimization/BayesianOptimization/actions/workflows/run_tests.yml/badge.svg)
+[![Codecov](https://codecov.io/github/bayesian-optimization/BayesianOptimization/badge.svg?branch=master&service=github)](https://codecov.io/github/bayesian-optimization/BayesianOptimization?branch=master)
 [![Pypi](https://img.shields.io/pypi/v/bayesian-optimization.svg)](https://pypi.python.org/pypi/bayesian-optimization)
 
 Pure Python implementation of bayesian global optimization with gaussian
 processes.
 
+## Installation
+
 * PyPI (pip):
 
 ```console
@@ -30,25 +32,18 @@ suited for optimization of high cost functions, situations where the balance
 between exploration and exploitation is important.
 
 ## Quick Start
-See below for a quick tour over the basics of the Bayesian Optimization package. More detailed information, other advanced features, and tips on usage/implementation can be found in the [examples](https://github.com/fmfn/BayesianOptimization/tree/master/examples) folder. I suggest that you:
-- Follow the
-[basic tour notebook](https://github.com/fmfn/BayesianOptimization/blob/master/examples/basic-tour.ipynb)
-to learn how to use the package's most important features.
-- Take a look at the
-[advanced tour notebook](https://github.com/fmfn/BayesianOptimization/blob/master/examples/advanced-tour.ipynb)
-to learn how to make the package more flexible, how to deal with categorical parameters, how to use observers, and more.
-- Check out this
-[notebook](https://github.com/fmfn/BayesianOptimization/blob/master/examples/visualization.ipynb)
-with a step by step visualization of how this method works.
-- To understand how to use bayesian optimization when additional constraints are present, see the
-[constrained optimization notebook](https://github.com/fmfn/BayesianOptimization/blob/master/examples/constraints.ipynb).
-- Explore this [notebook](https://github.com/fmfn/BayesianOptimization/blob/master/examples/exploitation_vs_exploration.ipynb)
+See below for a quick tour over the basics of the Bayesian Optimization package. More detailed information, other advanced features, and tips on usage/implementation can be found in the [examples](http://bayesian-optimization.github.io/BayesianOptimization/examples.html) folder. I suggest that you:
+- Follow the [basic tour notebook](http://bayesian-optimization.github.io/BayesianOptimization/basic-tour.html) to learn how to use the package's most important features.
+- Take a look at the [advanced tour notebook](http://bayesian-optimization.github.io/BayesianOptimization/advanced-tour.html) to learn how to make the package more flexible, how to deal with categorical parameters, how to use observers, and more.
+- Check out this [notebook](http://bayesian-optimization.github.io/BayesianOptimization/visualization.html) with a step by step visualization of how this method works.
+- To understand how to use bayesian optimization when additional constraints are present, see the [constrained optimization notebook](http://bayesian-optimization.github.io/BayesianOptimization/constraints.html).
+- Explore this [notebook](http://bayesian-optimization.github.io/BayesianOptimization/exploitation_vs_exploration.html)
 exemplifying the balance between exploration and exploitation and how to
 control it.
-- Go over this [script](https://github.com/fmfn/BayesianOptimization/blob/master/examples/sklearn_example.py)
+- Go over this [script](https://github.com/bayesian-optimization/BayesianOptimization/blob/master/examples/sklearn_example.py)
 for examples of how to tune parameters of Machine Learning models using cross validation and bayesian optimization.
-- Explore the [domain reduction notebook](https://github.com/fmfn/BayesianOptimization/blob/master/examples/domain_reduction.ipynb) to learn more about how search can be sped up by dynamically changing parameters' bounds.
-- Finally, take a look at this [script](https://github.com/fmfn/BayesianOptimization/blob/master/examples/async_optimization.py)
+- Explore the [domain reduction notebook](http://bayesian-optimization.github.io/BayesianOptimization/domain_reduction.html) to learn more about how search can be sped up by dynamically changing parameters' bounds.
+- Finally, take a look at this [script](https://github.com/bayesian-optimization/BayesianOptimization/blob/master/examples/async_optimization.py)
 for ideas on how to implement bayesian optimization in a distributed fashion using this package.
 
 
@@ -56,11 +51,11 @@ for ideas on how to implement bayesian optimization in a distributed fashion usi
 
 Bayesian optimization works by constructing a posterior distribution of functions (gaussian process) that best describes the function you want to optimize. As the number of observations grows, the posterior distribution improves, and the algorithm becomes more certain of which regions in parameter space are worth exploring and which are not, as seen in the picture below.
 
-![BayesianOptimization in action](./examples/bo_example.png)
+![BayesianOptimization in action](./static/bo_example.png)
 
 As you iterate over and over, the algorithm balances its needs of exploration and exploitation taking into account what it knows about the target function. At each step a Gaussian Process is fitted to the known samples (points previously explored), and the posterior distribution, combined with a exploration strategy (such as UCB (Upper Confidence Bound), or EI (Expected Improvement)), are used to determine the next point that should be explored (see the gif below).
 
-![BayesianOptimization in action](./examples/bayesian_optimization.gif)
+![BayesianOptimization in action](./static/bayesian_optimization.gif)
 
 This process is designed to minimize the number of steps required to find a combination of parameters that are close to the optimal combination. To do so, this method uses a proxy optimization problem (finding the maximum of the acquisition function) that, albeit still a hard problem, is cheaper (in the computational sense) and common tools can be employed. Therefore Bayesian Optimization is most adequate for situations where sampling the function to be optimized is a very expensive endeavor. See the references for a proper discussion of this method.
 
@@ -68,10 +63,9 @@ This project is under active development, if you find a bug, or anything that
 needs correction, please let me know.
 
 
-Basic tour of the Bayesian Optimization package
-===============================================
+## Basic tour of the Bayesian Optimization package
 
-## 1. Specifying the function to be optimized
+### 1. Specifying the function to be optimized
 
 This is a function optimization package, therefore the first and most important ingredient is, of course, the function to be optimized.
 
@@ -89,7 +83,7 @@ def black_box_function(x, y):
     return -x ** 2 - (y - 1) ** 2 + 1
 ```
 
-## 2. Getting Started
+### 2. Getting Started
 
 All we need to get started is to instantiate a `BayesianOptimization` object specifying a function to be optimized `f`, and its parameters with their corresponding bounds, `pbounds`. This is a constrained optimization technique, so you must specify the minimum and maximum values that can be probed for each parameter in order for it to work
 
@@ -160,7 +154,7 @@ for i, res in enumerate(optimizer.res):
 ```
 
 
-### 2.1 Changing bounds
+#### 2.1 Changing bounds
 
 During the optimization process you may realize the bounds chosen for some parameters are not adequate. For these situations you can invoke the method `set_bounds` to alter them. You can pass any combination of **existing** parameters and their associated new bounds.
 
@@ -183,17 +177,17 @@ optimizer.maximize(
     |  10       | -1.762    |  1.442    |  0.1735   |
     =================================================
 
-### 2.2 Sequential Domain Reduction
+#### 2.2 Sequential Domain Reduction
 
 Sometimes the initial boundaries specified for a problem are too wide, and adding points to improve the response surface in regions of the solution domain is extraneous. Other times the cost function is very expensive to compute, and minimizing the number of calls is extremely beneficial.
 
 When it's worthwhile to converge on an optimal point quickly rather than try to find the optimal point, contracting the domain around the current optimal value as the search progresses can speed up the search progress considerably. Using the `SequentialDomainReductionTransformer` the bounds of the problem can be panned and zoomed dynamically in an attempt to improve convergence.
 
-![sequential domain reduction](./examples/sdr.png)
+![sequential domain reduction](./static/sdr.png)
 
-An example of using the `SequentialDomainReductionTransformer` is shown in the [domain reduction notebook](https://github.com/fmfn/BayesianOptimization/blob/master/examples/domain_reduction.ipynb). More information about this method can be found in the paper ["On the robustness of a simple domain reduction scheme for simulation‐based optimization"](http://www.truegrid.com/srsm_revised.pdf).
+An example of using the `SequentialDomainReductionTransformer` is shown in the [domain reduction notebook](http://bayesian-optimization.github.io/BayesianOptimization/domain_reduction.html). More information about this method can be found in the paper ["On the robustness of a simple domain reduction scheme for simulation‐based optimization"](http://www.truegrid.com/srsm_revised.pdf).
 
-## 3. Guiding the optimization
+### 3. Guiding the optimization
 
 It is often the case that we have an idea of regions of the parameter space where the maximum of our function might lie. For these situations the `BayesianOptimization` object allows the user to specify points to be probed. By default these will be explored lazily (`lazy=True`), meaning these points will be evaluated only the next time you call `maximize`. This probing process happens before the gaussian process takes over.
 
@@ -221,11 +215,11 @@ optimizer.maximize(init_points=0, n_iter=0)
     =================================================
 
 
-## 4. Saving, loading and restarting
+### 4. Saving, loading and restarting
 
 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.
 
-### 4.1 Saving progress
+#### 4.1 Saving progress
 
 
 ```python
@@ -255,7 +249,7 @@ optimizer.maximize(
 
 By default the previous data in the json file is removed. If you want to keep working with the same logger, the `reset` parameter in `JSONLogger` should be set to False.
 
-### 4.2 Loading progress
+#### 4.2 Loading progress
 
 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.
 
@@ -277,54 +271,71 @@ load_logs(new_optimizer, logs=["./logs.log"]);
 
 ## Next Steps
 
-This introduction covered the most basic functionality of the package. Checkout the [basic-tour](https://github.com/fmfn/BayesianOptimization/blob/master/examples/basic-tour.ipynb) and [advanced-tour](https://github.com/fmfn/BayesianOptimization/blob/master/examples/advanced-tour.ipynb) notebooks in the example folder, where you will find detailed explanations and other more advanced functionality. Also, browse the examples folder for implementation tips and ideas.
-
-Installation
-============
-
-### Installation
-
-The latest release can be obtained by two ways:
-
-* With PyPI (pip):
-
-      pip install bayesian-optimization
+This introduction covered the most basic functionality of the package. Checkout the [basic-tour](http://bayesian-optimization.github.io/BayesianOptimization/basic-tour.html) and [advanced-tour](http://bayesian-optimization.github.io/BayesianOptimization/advanced-tour.html), where you will find detailed explanations and other more advanced functionality. Also, browse the [examples](http://bayesian-optimization.github.io/BayesianOptimization/examples.html) for implementation tips and ideas.
 
-* With conda (from conda-forge channel):
-
-      conda install -c conda-forge bayesian-optimization
-
-The bleeding edge version can be installed with:
-
-    pip install git+https://github.com/fmfn/BayesianOptimization.git
+## Installing from GitHub
+To install directly from master, simply run
+```
+pip install https://github.com/bayesian-optimization/BayesianOptimization/archive/master.zip
+```
+Please note that these builds are not necessarily stable.
 
-If you prefer, you can clone it and run the setup.py file. Use the following
-commands to get a copy from Github and install all dependencies:
+If you prefer, you can also clone the repository and run `setup.py`:
+```
+git clone https://github.com/bayesian-optimization/BayesianOptimization.git
+cd BayesianOptimization
+python setup.py install
+```
+    
 
-    git clone https://github.com/fmfn/BayesianOptimization.git
-    cd BayesianOptimization
-    python setup.py install
+## Minutiae
 
-Citation
-============
+### Citation
 
-If you used this package in your research and is interested in citing it here's how you do it:
+If you used this package in your research, please cite it:
 
 ```
 @Misc{,
     author = {Fernando Nogueira},
     title = {{Bayesian Optimization}: Open source constrained global optimization tool for {Python}},
     year = {2014--},
-    url = " https://github.com/fmfn/BayesianOptimization"
+    url = " https://github.com/bayesian-optimization/BayesianOptimization"
+}
+```
+If you used any of the advanced functionalities, please additionally cite the corresponding publication:
+
+For the `SequentialDomainTransformer`:
+```
+@article{
+    author = {Stander, Nielen and Craig, Kenneth},
+    year = {2002},
+    month = {06},
+    pages = {},
+    title = {On the robustness of a simple domain reduction scheme for simulation-based optimization},
+    volume = {19},
+    journal = {International Journal for Computer-Aided Engineering and Software (Eng. Comput.)},
+    doi = {10.1108/02644400210430190}
+}
+```
+
+For constrained optimization:
+```
+@inproceedings{gardner2014bayesian,
+    title={Bayesian optimization with inequality constraints.},
+    author={Gardner, Jacob R and Kusner, Matt J and Xu, Zhixiang Eddie and Weinberger, Kilian Q and Cunningham, John P},
+    booktitle={ICML},
+    volume={2014},
+    pages={937--945},
+    year={2014}
 }
 ```
 
-# Dependencies
+### Dependencies
 * Numpy
 * Scipy
 * Scikit-learn
  
-# References:
+### References:
 * http://papers.nips.cc/paper/4522-practical-bayesian-optimization-of-machine-learning-algorithms.pdf
 * http://arxiv.org/pdf/1012.2599v1.pdf
 * http://www.gaussianprocess.org/gpml/
diff --git a/bayes_opt/constraint.py b/bayes_opt/constraint.py
index 5ea333525..10142e43b 100644
--- a/bayes_opt/constraint.py
+++ b/bayes_opt/constraint.py
@@ -81,7 +81,7 @@ def eval(self, **kwargs):
         Raises
         ------
         TypeError
-            If the kwargs keys don't match the function argument names.
+            If the kwargs' keys don't match the function argument names.
         """
         try:
             return self.fun(**kwargs)
diff --git a/bayes_opt/domain_reduction.py b/bayes_opt/domain_reduction.py
index f4b86ffc8..f1d31a418 100644
--- a/bayes_opt/domain_reduction.py
+++ b/bayes_opt/domain_reduction.py
@@ -1,26 +1,37 @@
+"""Implement domain transformation.
+
+In particular, this provides a base transformer class and a sequential domain
+reduction transformer as based on Stander and Craig's "On the robustness of a
+simple domain reduction scheme for simulation-based optimization"
+"""
 from typing import Optional, Union, List
 
 import numpy as np
 from .target_space import TargetSpace
+from warnings import warn
 
 
 class DomainTransformer():
-    '''The base transformer class'''
+    """Base class."""
 
-    def __init__(self, **kwargs):
+    def __init__(self, **kwargs) -> None:
+        """To override with specific implementation."""
         pass
 
-    def initialize(self, target_space: TargetSpace):
+    def initialize(self, target_space: TargetSpace) -> None:
+        """To override with specific implementation."""
         raise NotImplementedError
 
-    def transform(self, target_space: TargetSpace):
+    def transform(self, target_space: TargetSpace) -> dict:
+        """To override with specific implementation."""
         raise NotImplementedError
 
 
 class SequentialDomainReductionTransformer(DomainTransformer):
-    """
+    """Reduce the searchable space.
+
     A sequential domain reduction transformer based on the work by Stander, N. and Craig, K:
-    "On the robustness of a simple domain reduction scheme for simulation‐based optimization"
+    "On the robustness of a simple domain reduction scheme for simulation-based optimization"
     """
 
     def __init__(
@@ -36,7 +47,14 @@ def __init__(
         self.minimum_window_value = minimum_window
 
     def initialize(self, target_space: TargetSpace) -> None:
-        """Initialize all of the parameters"""
+        """Initialize all of the parameters.
+
+        Parameters
+        ----------
+        target_space : TargetSpace
+            TargetSpace this DomainTransformer operates on.
+        """
+        # Set the original bounds
         self.original_bounds = np.copy(target_space.bounds)
         self.bounds = [self.original_bounds]
 
@@ -47,6 +65,7 @@ def initialize(self, target_space: TargetSpace) -> None:
         else:
             self.minimum_window = [self.minimum_window_value] * len(target_space.bounds)
 
+        # Set initial values
         self.previous_optimal = np.mean(target_space.bounds, axis=1)
         self.current_optimal = np.mean(target_space.bounds, axis=1)
         self.r = target_space.bounds[:, 1] - target_space.bounds[:, 0]
@@ -72,14 +91,12 @@ def initialize(self, target_space: TargetSpace) -> None:
         self._window_bounds_compatibility(self.original_bounds)
 
     def _update(self, target_space: TargetSpace) -> None:
-
+        """Update contraction rate, window size, and window center."""
         # setting the previous
         self.previous_optimal = self.current_optimal
         self.previous_d = self.current_d
-
-        self.current_optimal = target_space.params[
-            np.argmax(target_space.target)
-        ]
+        
+        self.current_optimal = target_space.params_to_array(target_space.max()['params'])
 
         self.current_d = 2.0 * (self.current_optimal -
                                 self.previous_optimal) / self.r
@@ -97,32 +114,82 @@ def _update(self, target_space: TargetSpace) -> None:
         self.r = self.contraction_rate * self.r
 
     def _trim(self, new_bounds: np.array, global_bounds: np.array) -> np.array:
-        for i, variable in enumerate(new_bounds):
-            if variable[0] < global_bounds[i, 0]:
-                variable[0] = global_bounds[i, 0]
-            if variable[1] > global_bounds[i, 1]:
-                variable[1] = global_bounds[i, 1]
-        for i, entry in enumerate(new_bounds):
-            if entry[0] > entry[1]:
-                new_bounds[i, 0] = entry[1]
-                new_bounds[i, 1] = entry[0]
-            window_width = abs(entry[0] - entry[1])
-            if window_width < self.minimum_window[i]:
-                dw = (self.minimum_window[i] - window_width) / 2.0
-                left_expansion_space = abs(global_bounds[i, 0] - entry[0]) # should be non-positive
-                right_expansion_space = abs(global_bounds[i, 1] - entry[1]) # should be non-negative
-                # conservative
-                dw_l = min(dw, left_expansion_space)
-                dw_r = min(dw, right_expansion_space)
-                # this crawls towards the edge
-                ddw_r = dw_r + max(dw - dw_l, 0)
-                ddw_l = dw_l + max(dw - dw_r, 0)
-                new_bounds[i, 0] -= ddw_l
-                new_bounds[i, 1] += ddw_r
+        """
+        Adjust the new_bounds and verify that they adhere to global_bounds and minimum_window.
+
+        Parameters
+        ----------
+        new_bounds : np.array
+            The proposed new_bounds that (may) need adjustment.
+
+        global_bounds : np.array
+            The maximum allowable bounds for each parameter.
+
+        Returns
+        -------
+        new_bounds : np.array
+            The adjusted bounds after enforcing constraints.
+        """
+        #sort bounds
+        new_bounds = np.sort(new_bounds)
+
+        # Validate each parameter's bounds against the global_bounds
+        for i, pbounds in enumerate(new_bounds):
+            # If the one of the bounds is outside the global bounds, reset the bound to the global bound
+            # This is expected to happen when the window is near the global bounds, no warning is issued
+            if (pbounds[0] < global_bounds[i, 0]):
+                pbounds[0] = global_bounds[i, 0]
+
+            if (pbounds[1] > global_bounds[i, 1]):
+                pbounds[1] = global_bounds[i, 1]
+
+            # If a lower bound is greater than the associated global upper bound, reset it to the global lower bound
+            if (pbounds[0] > global_bounds[i, 1]):
+                pbounds[0] = global_bounds[i, 0]
+                warn("\nDomain Reduction Warning:\n"+
+                    "A parameter's lower bound is greater than the global upper bound."+
+                    "The offensive boundary has been reset."+
+                    "Be cautious of subsequent reductions.", stacklevel=2)
+
+            # If an upper bound is less than the associated global lower bound, reset it to the global upper bound
+            if (pbounds[1] < global_bounds[i, 0]):
+                pbounds[1] = global_bounds[i, 1]
+                warn("\nDomain Reduction Warning:\n"+
+                    "A parameter's lower bound is greater than the global upper bound."+
+                    "The offensive boundary has been reset."+
+                    "Be cautious of subsequent reductions.", stacklevel=2)
+
+        # Adjust new_bounds to ensure they respect the minimum window width for each parameter
+        for i, pbounds in enumerate(new_bounds):
+            current_window_width = abs(pbounds[0] - pbounds[1])
+
+            # If the window width is less than the minimum allowable width, adjust it
+            # Note that when minimum_window < width of the global bounds one side always has more space than required
+            if current_window_width < self.minimum_window[i]:
+                width_deficit = (self.minimum_window[i] - current_window_width) / 2.0
+                available_left_space = abs(global_bounds[i, 0] - pbounds[0])
+                available_right_space = abs(global_bounds[i, 1] - pbounds[1])
+                
+                # determine how much to expand on the left and right
+                expand_left = min(width_deficit, available_left_space)
+                expand_right = min(width_deficit, available_right_space)
+                
+                # calculate the deficit on each side
+                expand_left_deficit = width_deficit - expand_left
+                expand_right_deficit = width_deficit - expand_right
+
+                # shift the deficit to the side with more space
+                adjust_left = expand_left + max(expand_right_deficit, 0)
+                adjust_right = expand_right + max(expand_left_deficit, 0)
+                
+                # adjust the bounds
+                pbounds[0] -= adjust_left
+                pbounds[1] += adjust_right
+
         return new_bounds
 
     def _window_bounds_compatibility(self, global_bounds: np.array) -> bool:
-        """Checks if global bounds are compatible with the minimum window sizes."""
+        """Check if global bounds are compatible with the minimum window sizes."""
         for i, entry in enumerate(global_bounds):
             global_window_width = abs(entry[1] - entry[0])
             if global_window_width < self.minimum_window[i]:
@@ -133,7 +200,7 @@ def _create_bounds(self, parameters: dict, bounds: np.array) -> dict:
         return {param: bounds[i, :] for i, param in enumerate(parameters)}
 
     def transform(self, target_space: TargetSpace) -> dict:
-
+        """Transform the bounds of the target space."""
         self._update(target_space)
 
         new_bounds = np.array(
@@ -143,6 +210,6 @@ def transform(self, target_space: TargetSpace) -> dict:
             ]
         ).T
 
-        self._trim(new_bounds, self.original_bounds)
+        new_bounds = self._trim(new_bounds, self.original_bounds)
         self.bounds.append(new_bounds)
         return self._create_bounds(target_space.keys, new_bounds)
diff --git a/bayes_opt/target_space.py b/bayes_opt/target_space.py
index 195474cb2..bf727c5c2 100644
--- a/bayes_opt/target_space.py
+++ b/bayes_opt/target_space.py
@@ -1,6 +1,7 @@
 """Manages the optimization domain and holds points."""
 import numpy as np
 from colorama import Fore
+from warnings import warn
 from .util import ensure_rng, NotUniqueError
 
 
@@ -38,9 +39,11 @@ class TargetSpace():
     >>>     return p1 + p2
     >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
     >>> space = TargetSpace(target_func, pbounds, random_state=0)
-    >>> x = space.random_points(1)[0]
-    >>> y = space.register_point(x)
-    >>> assert self.max_point()['max_val'] == y
+    >>> x = np.array([4 , 5])
+    >>> y = target_func(x)
+    >>> space.register(x, y)
+    >>> assert self.max()['target'] == 9
+    >>> assert self.max()['params'] == {'p1': 1.0, 'p2': 2.0}
     """
 
     def __init__(self, target_func, pbounds, constraint=None, random_state=None,
@@ -178,6 +181,30 @@ def constraint_values(self):
 
         return self._constraint_values
 
+    @property
+    def mask(self):
+        """Return a boolean array of valid points.
+         
+        Points are valid if they satisfy both the constraint and boundary conditions.
+
+        Returns
+        -------
+        np.ndarray
+        """
+        mask = np.ones_like(self.target, dtype=bool)
+
+        # mask points that don't satisfy the constraint
+        if self._constraint is not None:
+            mask &= self._constraint.allowed(self._constraint_values)
+
+        # mask points that are outside the bounds
+        if self._bounds is not None:
+            within_bounds = np.all((self._bounds[:, 0] <= self._params) & 
+                    (self._params <= self._bounds[:, 1]), axis=1)
+            mask &= within_bounds
+
+        return mask
+
     def params_to_array(self, params):
         """Convert a dict representation of parameters into an array version.
 
@@ -260,13 +287,14 @@ def register(self, params, target, constraint_value=None):
 
         Examples
         --------
+        >>> target_func = lambda p1, p2: p1 + p2
         >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
-        >>> space = TargetSpace(lambda p1, p2: p1 + p2, pbounds)
+        >>> space = TargetSpace(target_func, pbounds)
         >>> len(space)
         0
         >>> x = np.array([0, 0])
         >>> y = 1
-        >>> space.add_observation(x, y)
+        >>> space.register(x, y)
         >>> len(space)
         1
         """
@@ -281,6 +309,11 @@ def register(self, params, target, constraint_value=None):
                 raise NotUniqueError(f'Data point {x} is not unique. You can set'
                                      ' "allow_duplicate_points=True" to avoid this error')
 
+        # if x is not within the bounds of the parameter space, warn the user
+        if self._bounds is not None:
+            if not np.all((self._bounds[:, 0] <= x) & (x <= self._bounds[:, 1])):
+                warn(f'\nData point {x} is outside the bounds of the parameter space. ', stacklevel=2)
+
         self._params = np.concatenate([self._params, x.reshape(1, -1)])
         self._target = np.concatenate([self._target, [target]])
 
@@ -313,6 +346,15 @@ def probe(self, params):
         -------
         y : float
             target function value.
+
+        Example
+        -------
+        >>> target_func = lambda p1, p2: p1 + p2
+        >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
+        >>> space = TargetSpace(target_func, pbounds)
+        >>> space.probe([1, 5])
+        >>> assert self.max()['target'] == 6
+        >>> assert self.max()['params'] == {'p1': 1.0, 'p2': 5.0}
         """
         x = self._as_array(params)
         params = dict(zip(self._keys, x))
@@ -327,19 +369,20 @@ def probe(self, params):
         return target, constraint_value
 
     def random_sample(self):
-        """Create random points within the bounds of the space.
+        """
+        Sample a random point from within the bounds of the space.
 
         Returns
         -------
-        data: np.ndarray
-            [num x dim] array points with dimensions corresponding to `self._keys`
+        data: ndarray
+            [1 x dim] array with dimensions corresponding to `self._keys`
 
         Examples
         --------
         >>> target_func = lambda p1, p2: p1 + p2
         >>> pbounds = {'p1': (0, 1), 'p2': (1, 100)}
         >>> space = TargetSpace(target_func, pbounds, random_state=0)
-        >>> space.random_points(1)
+        >>> space.random_sample()
         array([[ 55.33253689,   0.54488318]])
         """
         data = np.empty((1, self.dim))
@@ -348,61 +391,48 @@ def random_sample(self):
         return data.ravel()
 
     def _target_max(self):
-        """Get maximum target value found.
+        """Get the maximum target value within the current parameter bounds.
         
         If there is a constraint present, the maximum value that fulfills the
-        constraint is returned.
+        constraint within the parameter bounds is returned.
 
-        
         Returns
         -------
-        dict | None
-            The maximum allowed point's target function value.
-            Returns None if there is no (allowed) maximum.
+        max: float
+            The maximum target value.
+        
         """
         if len(self.target) == 0:
             return None
 
-        if self._constraint is None:
-            return self.target.max()
-
-        allowed = self._constraint.allowed(self._constraint_values)
-        if allowed.any():
-            return self.target[allowed].max()
+        if len(self.target[self.mask]) == 0:
+            return None
 
-        return None
+        return self.target[self.mask].max()
 
     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 within the parameter bounds 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.
+        res: dict
+            A dictionary with the keys 'target' and 'params'. The value of
+            'target' is the maximum target value, and the value of 'params' is
+            a dictionary with the parameter names as keys and the parameter
+            values as values.
+        
         """
         target_max = self._target_max()
-
         if target_max is None:
             return None
 
-        if self._constraint is not None:
-            allowed = self._constraint.allowed(self._constraint_values)
-
-            target = self.target[allowed]
-            params = self.params[allowed]
-            constraint_values = self.constraint_values[allowed]
-        else:
-            target = self.target
-            params = self.params
-
-        target_max_idx = np.where(target == target_max)[0][0]
-
-
+        target = self.target[self.mask]
+        params = self.params[self.mask]
+        target_max_idx = np.argmax(target)
+        
         res = {
                 'target': target_max,
                 'params': dict(
@@ -411,6 +441,7 @@ def max(self):
         }
 
         if self._constraint is not None:
+            constraint_values = self.constraint_values[self.mask]
             res['constraint'] = constraint_values[target_max_idx]
 
         return res
@@ -420,9 +451,17 @@ def res(self):
 
         Returns
         -------
-        dict
-            A dictionary containing the target function values, parameters,
-            and, if applicable, constraint function values and allowedness.
+        res: list
+            A list of dictionaries with the keys 'target', 'params', and
+            'constraint'. The value of 'target' is the target value, the value
+            of 'params' is a dictionary with the parameter names as keys and the
+            parameter values as values, and the value of 'constraint' is the
+            constraint fulfillment.
+
+        Notes
+        -----
+        Does not report if points are within the bounds of the parameter space.
+
         """
         if self._constraint is None:
             params = [dict(zip(self.keys, p)) for p in self.params]
diff --git a/bayes_opt/util.py b/bayes_opt/util.py
index 417e307f7..a90f344d7 100644
--- a/bayes_opt/util.py
+++ b/bayes_opt/util.py
@@ -8,7 +8,7 @@
 
 def acq_max(ac, gp, y_max, bounds, random_state, constraint=None, n_warmup=10000, n_iter=10, y_max_params=None):
     """Find the maximum of the acquisition function.
-    
+
     It uses a combination of random sampling (cheap) and the 'L-BFGS-B'
     optimization method. First by sampling `n_warmup` (1e5) points at random,
     and then running L-BFGS-B from `n_iter` (10) random starting points.
@@ -106,7 +106,6 @@ def adjusted_ac(x):
     x_seeds = random_state.uniform(bounds[:, 0], bounds[:, 1],
                                    size=(1+n_iter+int(not y_max_params is None),
                                    bounds.shape[0]))
-
     x_seeds[0] = x_max
     if not y_max_params is None:
         # Add the provided best sample to the seeds so that the optimization
diff --git a/docsrc/Makefile b/docsrc/Makefile
index e0e46fb58..c3886de30 100644
--- a/docsrc/Makefile
+++ b/docsrc/Makefile
@@ -15,6 +15,8 @@ BUILDDIR      = ../docs
 .PHONY: help Makefile
 
 github:
+	@cp ../README.md .
+	@cp -r ../static .
 	@make html
 	@cp -a ../docs/html/. ../docs
 
diff --git a/docsrc/conf.py b/docsrc/conf.py
index b758e36f4..159db2748 100644
--- a/docsrc/conf.py
+++ b/docsrc/conf.py
@@ -42,7 +42,9 @@
     'sphinx.ext.githubpages',
     'nbsphinx',
     'IPython.sphinxext.ipython_console_highlighting',
-    'sphinx.ext.mathjax']
+    'sphinx.ext.mathjax',
+    'myst_parser',
+    'numpydoc']
 
 source_suffix = {
     '.rst': 'restructuredtext',
@@ -63,6 +65,7 @@
 # a list of builtin themes.
 #
 html_theme = 'sphinx_rtd_theme'
+html_favicon = 'func.ico'
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
diff --git a/docsrc/func.ico b/docsrc/func.ico
new file mode 100644
index 000000000..f29470261
Binary files /dev/null and b/docsrc/func.ico differ
diff --git a/docsrc/index.rst b/docsrc/index.rst
index 03e5da765..cfa017774 100644
--- a/docsrc/index.rst
+++ b/docsrc/index.rst
@@ -5,6 +5,7 @@ Bayesian Optimization
    :maxdepth: 2
    :caption: Contents:
 
+   /quickstart
    /examples
    /code_docs
 
diff --git a/docsrc/quickstart.rst b/docsrc/quickstart.rst
new file mode 100644
index 000000000..9f12a4c65
--- /dev/null
+++ b/docsrc/quickstart.rst
@@ -0,0 +1,2 @@
+.. include:: README.md
+   :parser: myst_parser.sphinx_
\ No newline at end of file
diff --git a/examples/domain_reduction.ipynb b/examples/domain_reduction.ipynb
index 4dd1041ad..81ec8f164 100644
--- a/examples/domain_reduction.ipynb
+++ b/examples/domain_reduction.ipynb
@@ -18,7 +18,7 @@
     "\n",
     "**Zoom**: contract the region of interest.\n",
     "\n",
-    "![](sdr.png)\n",
+    "![](../static/sdr.png)\n",
     "\n",
     "\n",
     "## Parameters\n",
diff --git a/examples/bayesian_optimization.gif b/static/bayesian_optimization.gif
similarity index 100%
rename from examples/bayesian_optimization.gif
rename to static/bayesian_optimization.gif
diff --git a/examples/bo_example.png b/static/bo_example.png
similarity index 100%
rename from examples/bo_example.png
rename to static/bo_example.png
diff --git a/examples/func.png b/static/func.png
similarity index 100%
rename from examples/func.png
rename to static/func.png
diff --git a/examples/sdr.png b/static/sdr.png
similarity index 100%
rename from examples/sdr.png
rename to static/sdr.png
diff --git a/tests/test_logs_bounds.log b/tests/test_logs_bounds.log
new file mode 100644
index 000000000..0a7d9ba46
--- /dev/null
+++ b/tests/test_logs_bounds.log
@@ -0,0 +1,5 @@
+{"datetime": {"delta": 0.0, "datetime": "2018-11-25 08:29:25", "elapsed": 0.0}, "params": {"y": 0, "x": 0}, "target": 0}
+{"datetime": {"delta": 0.001301, "datetime": "2018-11-25 08:29:25", "elapsed": 0.001301}, "params": {"y": 1, "x": 1}, "target": 2}
+{"datetime": {"delta": 1.075242, "datetime": "2018-11-25 08:29:26", "elapsed": 1.076543}, "params": {"y": 2, "x": 2}, "target": 4}
+{"datetime": {"delta": 0.239797, "datetime": "2018-11-25 08:29:26", "elapsed": 1.31634}, "params": {"y": 3, "x": 3}, "target": 6}
+{"datetime": {"delta": 0.001301, "datetime": "2018-11-25 08:29:25", "elapsed": 0.001301}, "params": {"y": 1.5, "x": 1.5}, "target": 3}
diff --git a/tests/test_seq_domain_red.py b/tests/test_seq_domain_red.py
index 4b60bf049..f03dc71dc 100644
--- a/tests/test_seq_domain_red.py
+++ b/tests/test_seq_domain_red.py
@@ -68,6 +68,7 @@ def reset(self):
     assert not (standard_optimizer._space.bounds ==
                 mutated_optimizer._space.bounds).any()
 
+
 def test_minimum_window_is_kept():
     bounds_transformer = SequentialDomainReductionTransformer(minimum_window=1.0)
     pbounds = {'x': (-0.5, 0.5), 'y': (-1.0, 0.0)}
@@ -106,6 +107,7 @@ def test_minimum_window_array_is_kept():
     window_widths = np.diff(bounds_transformer.bounds)
     assert np.all(np.isclose(np.squeeze(np.min(window_widths, axis=0)), window_ranges))
 
+
 def test_trimming_bounds():
     """Test if the bounds are trimmed correctly within the bounds"""
     def dummy_function(x1, x2, x3, x4, x5):
@@ -143,4 +145,51 @@ def test_exceeded_bounds():
                 verbose=0,
                 random_state=1,
                 bounds_transformer=bounds_transformer
-            )
\ No newline at end of file
+            )
+
+
+def test_trim_when_both_new_bounds_exceed_global_bounds():
+    """Test if the global bounds are respected when both new bounds for a given parameter
+    are beyond the global bounds."""
+
+    # initialize a bounds transformer
+    bounds_transformer = SequentialDomainReductionTransformer(minimum_window=10)
+    pbounds  = {'x': (-10, 10),'y': (-10, 10)}
+    target_sp = TargetSpace(target_func=black_box_function, pbounds=pbounds)
+    bounds_transformer.initialize(target_sp)
+    global_bounds = np.asarray(list(pbounds.values()))
+
+    def verify_bounds_in_range(new_bounds, global_bounds):
+        """Check if the new bounds are within the global bounds."""
+        test = True
+        for i, pbounds in enumerate(new_bounds):
+            if (pbounds[0] < global_bounds[i, 0] or pbounds[0] > global_bounds[i, 1]):
+                test = False
+            if (pbounds[1] > global_bounds[i, 1] or pbounds[1] < global_bounds[i, 0]):
+                test = False
+        return test
+
+    # test if the sorting of the bounds is correct
+    new_bounds = np.array( [[5, -5], [-10, 10]] )
+    trimmed_bounds = bounds_transformer._trim(new_bounds, global_bounds)
+    assert (trimmed_bounds == np.array( [[-5, 5], [-10, 10]] )).all()
+
+    # test if both (upper/lower) bounds for a parameter exceed the global bounds
+    new_bounds = np.array( [[-50, -20], [20, 50]] )
+    with pytest.warns(UserWarning):
+        trimmed_bounds = bounds_transformer._trim(new_bounds, global_bounds)
+    assert verify_bounds_in_range(trimmed_bounds, global_bounds)
+
+    # test if both (upper/lower) bounds for a parameter exceed the global bounds 
+    # while they are out of order
+    new_bounds = np.array( [[-20, -50], [-10, 10]] )
+    with pytest.warns(UserWarning):
+        trimmed_bounds = bounds_transformer._trim(new_bounds, global_bounds)
+    assert verify_bounds_in_range(trimmed_bounds, global_bounds)
+
+if __name__ == '__main__':
+    r"""
+    CommandLine:
+        python tests/test_seq_domain_red.py
+    """
+    pytest.main([__file__])
\ No newline at end of file
diff --git a/tests/test_target_space.py b/tests/test_target_space.py
index c554dea0e..494b42495 100644
--- a/tests/test_target_space.py
+++ b/tests/test_target_space.py
@@ -74,6 +74,7 @@ def test_as_array():
 
 
 def test_register():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     space = TargetSpace(target_func, PBOUNDS)
 
     assert len(space) == 0
@@ -96,6 +97,7 @@ def test_register():
 
 
 def test_register_with_constraint():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     constraint = ConstraintModel(lambda x: x, -2, 2)
     space = TargetSpace(target_func, PBOUNDS, constraint=constraint)
 
@@ -118,7 +120,16 @@ def test_register_with_constraint():
         space.register(params={"p1": 2, "p2": 2}, target=3)
 
 
+def test_register_point_beyond_bounds():
+    PBOUNDS = {'p1': (0, 1), 'p2': (1, 10)}
+    space = TargetSpace(target_func, PBOUNDS)
+
+    with pytest.warns(UserWarning):
+        space.register(params={"p1": 0.5, "p2": 20}, target=2.5)
+
+
 def test_probe():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     space = TargetSpace(target_func, PBOUNDS, allow_duplicate_points=True)
 
     assert len(space) == 0
@@ -166,12 +177,14 @@ def test_random_sample():
 def test_y_max():
     space = TargetSpace(target_func, PBOUNDS)
     assert space._target_max() == None
-    space.probe(params={"p1": 1, "p2": 2})
-    space.probe(params={"p1": 5, "p2": 1})
+    space.probe(params={"p1": 1, "p2": 7})
+    space.probe(params={"p1": 0.5, "p2": 1})
     space.probe(params={"p1": 0, "p2": 1})
-    assert space._target_max() == 6
+    assert space._target_max() == 8
+
 
 def test_y_max_with_constraint():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     constraint = ConstraintModel(lambda p1, p2: p1-p2, -2, 2)
     space = TargetSpace(target_func, PBOUNDS, constraint)
     assert space._target_max() == None
@@ -181,10 +194,21 @@ def test_y_max_with_constraint():
     assert space._target_max() == 3
 
 
+def test_y_max_within_pbounds():
+    PBOUNDS = {'p1': (0, 2), 'p2': (1, 100)}
+    space = TargetSpace(target_func, PBOUNDS)
+    assert space._target_max() == None
+    space.probe(params={"p1": 1, "p2": 2})
+    space.probe(params={"p1": 0, "p2": 1})
+    with pytest.warns(UserWarning):
+        space.probe(params={"p1": 5, "p2": 1})
+    assert space._target_max() == 3
+
 
 def test_max():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     space = TargetSpace(target_func, PBOUNDS)
-
+    
     assert space.max() == None
     space.probe(params={"p1": 1, "p2": 2})
     space.probe(params={"p1": 5, "p2": 4})
@@ -192,7 +216,9 @@ def test_max():
     space.probe(params={"p1": 1, "p2": 6})
     assert space.max() == {"params": {"p1": 5, "p2": 4}, "target": 9}
 
+
 def test_max_with_constraint():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     constraint = ConstraintModel(lambda p1, p2: p1-p2, -2, 2)
     space = TargetSpace(target_func, PBOUNDS, constraint=constraint)
 
@@ -203,7 +229,9 @@ def test_max_with_constraint():
     space.probe(params={"p1": 1, "p2": 6}) # Unfeasible
     assert space.max() == {"params": {"p1": 2, "p2": 3}, "target": 5, "constraint": -1}
 
+
 def test_max_with_constraint_identical_target_value():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     constraint = ConstraintModel(lambda p1, p2: p1-p2, -2, 2)
     space = TargetSpace(target_func, PBOUNDS, constraint=constraint)
 
@@ -215,7 +243,9 @@ def test_max_with_constraint_identical_target_value():
     space.probe(params={"p1": 1, "p2": 6}) # Unfeasible
     assert space.max() == {"params": {"p1": 2, "p2": 3}, "target": 5, "constraint": -1}
 
+
 def test_res():
+    PBOUNDS = {'p1': (0, 10), 'p2': (1, 100)}
     space = TargetSpace(target_func, PBOUNDS)
 
     assert space.res() == []
diff --git a/tests/test_util.py b/tests/test_util.py
index b3a73668a..89b79c01e 100644
--- a/tests/test_util.py
+++ b/tests/test_util.py
@@ -132,7 +132,7 @@ def f(x, y):
 
     optimizer = BayesianOptimization(
         f=f,
-        pbounds={"x": (-2, 2), "y": (-2, 2)}
+        pbounds={"x": (-200, 200), "y": (-200, 200)}
     )
     assert len(optimizer.space) == 0
 
@@ -150,6 +150,21 @@ def f(x, y):
         load_logs(other_optimizer, [str(test_dir / "test_logs.log")])
 
 
+def test_logs_bounds():
+    def f(x, y):
+        return x + y
+
+    optimizer = BayesianOptimization(
+        f=f,
+        pbounds={"x": (-2, 2), "y": (-2, 2)}
+    )
+
+    with pytest.warns(UserWarning):
+        load_logs(optimizer, [str(test_dir / "test_logs_bounds.log")])
+    
+    assert len(optimizer.space) == 5
+
+
 def test_logs_constraint():
 
     def f(x, y):
@@ -162,7 +177,7 @@ def c(x, y):
 
     optimizer = BayesianOptimization(
         f=f,
-        pbounds={"x": (-2, 2), "y": (-2, 2)},
+        pbounds={"x": (-200, 200), "y": (-200, 200)},
         constraint=constraint
     )
 
@@ -179,5 +194,4 @@ def c(x, y):
     CommandLine:
         python tests/test_target_space.py
     """
-    import pytest
-    pytest.main([__file__])
+    pytest.main([__file__])
\ No newline at end of file