diff --git a/.isort.cfg b/.isort.cfg
new file mode 100644
index 00000000..e38a4fd5
--- /dev/null
+++ b/.isort.cfg
@@ -0,0 +1,4 @@
+[settings]
+line_length=79
+multi_line_output=3
+known_third_party = h5py,matplotlib,numpy,pytest,scipy,setuptools
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 00000000..b267469d
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,30 @@
+repos:
+- repo: https://github.com/pre-commit/pre-commit-hooks
+  rev: v2.3.0
+  hooks:
+  - id: trailing-whitespace
+  - id: end-of-file-fixer
+  - id: debug-statements
+
+- repo: https://github.com/asottile/seed-isort-config
+  rev: v1.9.2
+  hooks:
+  - id: seed-isort-config
+    args: [--application-directories=src]
+
+- repo: https://github.com/pre-commit/mirrors-isort
+  rev: v4.3.20
+  hooks:
+  - id: isort
+    args: []
+
+- repo: https://github.com/psf/black
+  rev: 19.3b0
+  hooks:
+  - id: black
+    language_version: python3.6
+
+- repo: https://github.com/dfm/black_nbconvert
+  rev: v0.1.1
+  hooks:
+  - id: black_nbconvert
diff --git a/MANIFEST.in b/MANIFEST.in
index b4b20c74..624fcb9f 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1 +1 @@
-include README.rst LICENSE HISTORY.rst AUTHORS.rst
+include LICENSE *.rst
diff --git a/docs/conf.py b/docs/conf.py
index b025772f..e192a6be 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -3,11 +3,12 @@
 import os
 import sys
 
+from emcee import __version__  # NOQA
+
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
-sys.path.insert(0, os.path.abspath(".."))
-from emcee import __version__  # NOQA
+sys.path.insert(0, os.path.join(os.path.abspath(".."), "src"))
 
 extensions = [
     "sphinx.ext.autodoc",
diff --git a/document/plots/oned.py b/document/plots/oned.py
index cf5da2f7..f9860463 100644
--- a/document/plots/oned.py
+++ b/document/plots/oned.py
@@ -1,15 +1,15 @@
 import os
 import sys
 import time
+from multiprocessing import Pool
 
-import numpy as np
-import matplotlib.pyplot as pl
 import h5py
+import matplotlib.pyplot as pl
+import numpy as np
 
-from multiprocessing import Pool
+import emcee
 
 sys.path.append(os.path.abspath(os.path.join(__file__, "..", "..", "..")))
-import emcee
 
 # import acor
 
@@ -20,8 +20,9 @@ def lnprobfn(p, icov):
 
 def random_cov(ndim, dof=1):
     v = np.random.randn(ndim * (ndim + dof)).reshape((ndim + dof, ndim))
-    return (sum([np.outer(v[i], v[i]) for i in range(ndim + dof)])
-            / (ndim + dof))
+    return sum([np.outer(v[i], v[i]) for i in range(ndim + dof)]) / (
+        ndim + dof
+    )
 
 
 _rngs = {}
@@ -31,33 +32,34 @@ def _worker(args):
     i, outfn, nsteps = args
 
     pid = os.getpid()
-    _random = _rngs.get(pid, np.random.RandomState(int(int(pid)
-        + time.time())))
+    _random = _rngs.get(
+        pid, np.random.RandomState(int(int(pid) + time.time()))
+    )
     _rngs[pid] = _random
 
     ndim = int(np.ceil(2 ** (7 * _random.rand())))
     nwalkers = 2 * ndim + 2
     # nwalkers += nwalkers % 2
-    print ndim, nwalkers
+    print(ndim, nwalkers)
 
     cov = random_cov(ndim)
     icov = np.linalg.inv(cov)
 
-    ens_samp = emcee.EnsembleSampler(nwalkers, ndim, lnprobfn,
-            args=[icov])
+    ens_samp = emcee.EnsembleSampler(nwalkers, ndim, lnprobfn, args=[icov])
     ens_samp.random_state = _random.get_state()
-    pos, lnprob, state = ens_samp.run_mcmc(np.random.randn(nwalkers * ndim)
-            .reshape([nwalkers, ndim]), nsteps)
+    pos, lnprob, state = ens_samp.run_mcmc(
+        np.random.randn(nwalkers * ndim).reshape([nwalkers, ndim]), nsteps
+    )
 
     proposal = np.diag(cov.diagonal())
-    mh_samp = emcee.MHSampler(proposal, ndim, lnprobfn,
-            args=[icov])
+    mh_samp = emcee.MHSampler(proposal, ndim, lnprobfn, args=[icov])
     mh_samp.random_state = state
     mh_samp.run_mcmc(np.random.randn(ndim), nsteps)
 
     f = h5py.File(outfn)
-    f["data"][i, :] = np.array([ndim, np.mean(ens_samp.acor),
-                                np.mean(mh_samp.acor)])
+    f["data"][i, :] = np.array(
+        [ndim, np.mean(ens_samp.acor), np.mean(mh_samp.acor)]
+    )
     f.close()
 
 
@@ -67,7 +69,7 @@ def oned():
     nthreads = 2
 
     outfn = os.path.join(os.path.split(__file__)[0], "gauss_scaling.h5")
-    print outfn
+    print(outfn)
     f = h5py.File(outfn, "w")
     f.create_dataset("data", (niter, 3), "f")
     f.close()
diff --git a/emcee/__init__.py b/emcee/__init__.py
deleted file mode 100644
index d7c45b6a..00000000
--- a/emcee/__init__.py
+++ /dev/null
@@ -1,32 +0,0 @@
-# -*- coding: utf-8 -*-
-
-from __future__ import print_function, absolute_import
-
-__version__ = "3.0.0"
-__bibtex__ = """
-@article{emcee,
-   author = {{Foreman-Mackey}, D. and {Hogg}, D.~W. and {Lang}, D. and {Goodman}, J.},
-    title = {emcee: The MCMC Hammer},
-  journal = {PASP},
-     year = 2013,
-   volume = 125,
-    pages = {306-312},
-   eprint = {1202.3665},
-      doi = {10.1086/670067}
-}
-"""  # NOQA
-
-try:
-    __EMCEE_SETUP__
-except NameError:
-    __EMCEE_SETUP__ = False
-
-if not __EMCEE_SETUP__:
-    from .ensemble import EnsembleSampler
-    from .state import State
-
-    from . import moves
-    from . import autocorr
-    from . import backends
-
-    __all__ = ["EnsembleSampler", "State", "moves", "autocorr", "backends"]
diff --git a/emcee/model.py b/emcee/model.py
deleted file mode 100644
index d9c67952..00000000
--- a/emcee/model.py
+++ /dev/null
@@ -1,12 +0,0 @@
-# -*- coding: utf-8 -*-
-
-from __future__ import division, print_function
-
-__all__ = ["Model"]
-
-from collections import namedtuple
-
-Model = namedtuple(
-    "Model",
-    ("log_prob_fn", "compute_log_prob_fn", "map_fn", "random")
-)
diff --git a/emcee/tests/integration/test_gaussian.py b/emcee/tests/integration/test_gaussian.py
deleted file mode 100644
index 6a11b0d1..00000000
--- a/emcee/tests/integration/test_gaussian.py
+++ /dev/null
@@ -1,62 +0,0 @@
-# -*- coding: utf-8 -*-
-
-from __future__ import division, print_function
-
-import pytest
-import numpy as np
-from itertools import product
-
-from emcee import moves
-
-from .test_proposal import _test_normal, _test_uniform
-
-__all__ = ["test_normal_gaussian", "test_uniform_gaussian",
-           "test_normal_gaussian_nd"]
-
-
-@pytest.mark.parametrize("mode,factor", product(
-    ["vector"], [None, 2.0, 5.0]
-))
-def test_normal_gaussian(mode, factor, **kwargs):
-    _test_normal(moves.GaussianMove(0.5, mode=mode, factor=factor), **kwargs)
-
-
-@pytest.mark.parametrize("mode,factor", product(
-    ["vector", "random", "sequential"], [None, 2.0]
-))
-def test_normal_gaussian_nd(mode, factor, **kwargs):
-    ndim = 3
-    kwargs["nsteps"] = 8000
-
-    # Isotropic.
-    _test_normal(moves.GaussianMove(0.5, factor=factor, mode=mode), ndim=ndim,
-                 **kwargs)
-
-    # Axis-aligned.
-    _test_normal(moves.GaussianMove(0.5 * np.ones(ndim), factor=factor,
-                                    mode=mode), ndim=ndim, **kwargs)
-    with pytest.raises(ValueError):
-        _test_normal(moves.GaussianMove(0.5 * np.ones(ndim-1), factor=factor,
-                                        mode=mode), ndim=ndim,
-                     **kwargs)
-
-    # Full matrix.
-    if mode == "vector":
-        _test_normal(moves.GaussianMove(np.diag(0.5 * np.ones(ndim)),
-                                        factor=factor, mode=mode), ndim=ndim,
-                     **kwargs)
-        with pytest.raises(ValueError):
-            _test_normal(moves.GaussianMove(np.diag(0.5 * np.ones(ndim-1))),
-                         ndim=ndim, **kwargs)
-    else:
-        with pytest.raises(ValueError):
-            _test_normal(moves.GaussianMove(np.diag(0.5 * np.ones(ndim)),
-                                            factor=factor, mode=mode),
-                         ndim=ndim, **kwargs)
-
-
-@pytest.mark.parametrize("mode,factor", product(
-    ["vector"], [None, 2.0, 5.0]
-))
-def test_uniform_gaussian(mode, factor, **kwargs):
-    _test_uniform(moves.GaussianMove(0.5, factor=factor, mode=mode), **kwargs)
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..dd98b06c
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,20 @@
+[build-system]
+requires = ["setuptools>=40.6.0", "wheel", "setuptools_scm"]
+build-backend = "setuptools.build_meta"
+
+[tool.black]
+line-length = 79
+exclude = '''
+/(
+    \.eggs
+  | \.git
+  | \.hg
+  | \.mypy_cache
+  | \.tox
+  | \.venv
+  | _build
+  | buck-out
+  | build
+  | dist
+)/
+'''
diff --git a/setup.py b/setup.py
index c0f992d8..1dd9d11b 100755
--- a/setup.py
+++ b/setup.py
@@ -1,48 +1,72 @@
-#! /usr/bin/env python
-# -*- coding: utf-8 -*-
+#!/usr/bin/env python
 
+# Inspired by:
+# https://hynek.me/articles/sharing-your-labor-of-love-pypi-quick-and-dirty/
+
+import codecs
 import os
-import sys
-from setuptools import setup
-
-if sys.argv[-1] == "publish":
-    os.system("python setup.py sdist; twine upload dist/*")
-    sys.exit()
-
-# Hackishly inject a constant into builtins to enable importing of the
-# package before all the dependencies are built.
-if sys.version_info[0] < 3:
-    import __builtin__ as builtins
-else:
-    import builtins
-builtins.__EMCEE_SETUP__ = True
-import emcee  # NOQA
-
-setup(
-    name="emcee",
-    version=emcee.__version__,
-    author="Daniel Foreman-Mackey",
-    author_email="foreman.mackey@gmail.com",
-    packages=[
-        "emcee",
-        "emcee.moves", "emcee.backends",
-        "emcee.tests", "emcee.tests.unit", "emcee.tests.integration",
-    ],
-    url="http://emcee.readthedocs.io",
-    license="MIT",
-    description=("The Python ensemble sampling toolkit for affine-invariant "
-                 "MCMC"),
-    long_description=open("README.rst").read(),
-    package_data={"": ["LICENSE", "AUTHORS.rst"]},
-    include_package_data=True,
-    install_requires=["numpy"],
-    classifiers=[
-        "Development Status :: 5 - Production/Stable",
-        "Intended Audience :: Developers",
-        "Intended Audience :: Science/Research",
-        "License :: OSI Approved :: MIT License",
-        "Operating System :: OS Independent",
-        "Programming Language :: Python",
-    ],
-    zip_safe=True,
-)
+import re
+
+from setuptools import find_packages, setup
+
+# PROJECT SPECIFIC
+
+NAME = "emcee"
+PACKAGES = find_packages(where="src")
+META_PATH = os.path.join("src", "emcee", "__init__.py")
+CLASSIFIERS = [
+    "Development Status :: 5 - Production/Stable",
+    "Intended Audience :: Developers",
+    "Intended Audience :: Science/Research",
+    "License :: OSI Approved :: MIT License",
+    "Operating System :: OS Independent",
+    "Programming Language :: Python",
+]
+INSTALL_REQUIRES = ["numpy"]
+
+# END PROJECT SPECIFIC
+
+
+HERE = os.path.dirname(os.path.realpath(__file__))
+
+
+def read(*parts):
+    with codecs.open(os.path.join(HERE, *parts), "rb", "utf-8") as f:
+        return f.read()
+
+
+def find_meta(meta, meta_file=read(META_PATH)):
+    meta_match = re.search(
+        r"^__{meta}__ = ['\"]([^'\"]*)['\"]".format(meta=meta), meta_file, re.M
+    )
+    if meta_match:
+        return meta_match.group(1)
+    raise RuntimeError("Unable to find __{meta}__ string.".format(meta=meta))
+
+
+if __name__ == "__main__":
+    setup(
+        name=NAME,
+        use_scm_version={
+            "write_to": os.path.join(
+                "src", NAME, "{0}_version.py".format(NAME)
+            ),
+            "write_to_template": '__version__ = "{version}"\n',
+        },
+        author=find_meta("author"),
+        author_email=find_meta("email"),
+        maintainer=find_meta("author"),
+        maintainer_email=find_meta("email"),
+        url=find_meta("uri"),
+        license=find_meta("license"),
+        description=find_meta("description"),
+        long_description=read("README.rst"),
+        long_description_content_type="text/x-rst",
+        packages=PACKAGES,
+        package_dir={"": "src"},
+        include_package_data=True,
+        install_requires=INSTALL_REQUIRES,
+        classifiers=CLASSIFIERS,
+        zip_safe=False,
+        options={"bdist_wheel": {"universal": "1"}},
+    )
diff --git a/src/emcee/__init__.py b/src/emcee/__init__.py
new file mode 100644
index 00000000..66cc5579
--- /dev/null
+++ b/src/emcee/__init__.py
@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+
+__bibtex__ = """
+@article{emcee,
+   author = {{Foreman-Mackey}, D. and {Hogg}, D.~W. and {Lang}, D. and {Goodman}, J.},
+    title = {emcee: The MCMC Hammer},
+  journal = {PASP},
+     year = 2013,
+   volume = 125,
+    pages = {306-312},
+   eprint = {1202.3665},
+      doi = {10.1086/670067}
+}
+"""
+__uri__ = "https://emcee.readthedocs.io"
+__author__ = "Daniel Foreman-Mackey"
+__email__ = "foreman.mackey@gmail.com"
+__license__ = "MIT"
+__description__ = (
+    "The Python ensemble sampling toolkit for affine-invariant MCMC"
+)
+
+
+from . import autocorr, backends, moves
+from .emcee_version import __version__
+from .ensemble import EnsembleSampler
+from .state import State
+
+__all__ = [
+    "EnsembleSampler",
+    "State",
+    "moves",
+    "autocorr",
+    "backends",
+    "__version__",
+]
diff --git a/emcee/autocorr.py b/src/emcee/autocorr.py
similarity index 93%
rename from emcee/autocorr.py
rename to src/emcee/autocorr.py
index 1b85a5ef..b662840c 100644
--- a/emcee/autocorr.py
+++ b/src/emcee/autocorr.py
@@ -1,7 +1,5 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import logging
 
 import numpy as np
@@ -33,8 +31,8 @@ def function_1d(x):
     n = next_pow_two(len(x))
 
     # Compute the FFT and then (from that) the auto-correlation function
-    f = np.fft.fft(x - np.mean(x), n=2*n)
-    acf = np.fft.ifft(f * np.conjugate(f))[:len(x)].real
+    f = np.fft.fft(x - np.mean(x), n=2 * n)
+    acf = np.fft.ifft(f * np.conjugate(f))[: len(x)].real
     acf /= acf[0]
     return acf
 
@@ -93,7 +91,7 @@ def integrated_time(x, c=5, tol=50, quiet=False):
         for k in range(n_w):
             f += function_1d(x[:, k, d])
         f /= n_w
-        taus = 2.0*np.cumsum(f)-1.0
+        taus = 2.0 * np.cumsum(f) - 1.0
         windows[d] = auto_window(taus, c)
         tau_est[d] = taus[windows[d]]
 
@@ -107,7 +105,7 @@ def integrated_time(x, c=5, tol=50, quiet=False):
             "autocorrelation time for {1} parameter(s). Use this estimate "
             "with caution and run a longer chain!\n"
         ).format(tol, np.sum(flag))
-        msg += "N/{0} = {1:.0f};\ntau: {2}".format(tol, n_t/tol, tau_est)
+        msg += "N/{0} = {1:.0f};\ntau: {2}".format(tol, n_t / tol, tau_est)
         if not quiet:
             raise AutocorrError(tau_est, msg)
         logging.warning(msg)
@@ -122,6 +120,7 @@ class AutocorrError(Exception):
     ``tau`` attribute of this exception.
 
     """
+
     def __init__(self, tau, *args, **kwargs):
         self.tau = tau
         super(AutocorrError, self).__init__(*args, **kwargs)
diff --git a/emcee/backends/__init__.py b/src/emcee/backends/__init__.py
similarity index 68%
rename from emcee/backends/__init__.py
rename to src/emcee/backends/__init__.py
index 0a8fa16a..d1d0d275 100644
--- a/emcee/backends/__init__.py
+++ b/src/emcee/backends/__init__.py
@@ -1,16 +1,10 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-__all__ = [
-    "Backend",
-    "HDFBackend", "TempHDFBackend",
-    "get_test_backends",
-]
-
 from .backend import Backend
 from .hdf import HDFBackend, TempHDFBackend
 
+__all__ = ["Backend", "HDFBackend", "TempHDFBackend", "get_test_backends"]
+
 
 def get_test_backends():
     backends = [Backend]
diff --git a/emcee/backends/backend.py b/src/emcee/backends/backend.py
similarity index 85%
rename from emcee/backends/backend.py
rename to src/emcee/backends/backend.py
index c33ad123..e7808374 100644
--- a/emcee/backends/backend.py
+++ b/src/emcee/backends/backend.py
@@ -1,14 +1,12 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-__all__ = ["Backend"]
-
 import numpy as np
 
 from .. import autocorr
 from ..state import State
 
+__all__ = ["Backend"]
+
 
 class Backend(object):
     """A simple default backend that stores the chain in memory"""
@@ -40,14 +38,16 @@ def has_blobs(self):
 
     def get_value(self, name, flat=False, thin=1, discard=0):
         if self.iteration <= 0:
-            raise AttributeError("you must run the sampler with "
-                                 "'store == True' before accessing the "
-                                 "results")
+            raise AttributeError(
+                "you must run the sampler with "
+                "'store == True' before accessing the "
+                "results"
+            )
 
         if name == "blobs" and not self.has_blobs():
             return None
 
-        v = getattr(self, name)[discard+thin-1:self.iteration:thin]
+        v = getattr(self, name)[discard + thin - 1 : self.iteration : thin]
         if flat:
             s = list(v.shape[1:])
             s[0] = np.prod(v.shape[:2])
@@ -108,16 +108,18 @@ def get_log_prob(self, **kwargs):
     def get_last_sample(self):
         """Access the most recent sample in the chain"""
         if (not self.initialized) or self.iteration <= 0:
-            raise AttributeError("you must run the sampler with "
-                                 "'store == True' before accessing the "
-                                 "results")
+            raise AttributeError(
+                "you must run the sampler with "
+                "'store == True' before accessing the "
+                "results"
+            )
         it = self.iteration
-        blobs = self.get_blobs(discard=it-1)
+        blobs = self.get_blobs(discard=it - 1)
         if blobs is not None:
             blobs = blobs[0]
         return State(
-            self.get_chain(discard=it-1)[0],
-            log_prob=self.get_log_prob(discard=it-1)[0],
+            self.get_chain(discard=it - 1)[0],
+            log_prob=self.get_log_prob(discard=it - 1)[0],
             blobs=blobs,
             random_state=self.random_state,
         )
@@ -184,21 +186,27 @@ def _check(self, state, accepted):
         nwalkers, ndim = self.shape
         has_blobs = self.has_blobs()
         if state.coords.shape != (nwalkers, ndim):
-            raise ValueError("invalid coordinate dimensions; expected {0}"
-                             .format((nwalkers, ndim)))
-        if state.log_prob.shape != (nwalkers, ):
-            raise ValueError("invalid log probability size; expected {0}"
-                             .format(nwalkers))
+            raise ValueError(
+                "invalid coordinate dimensions; expected {0}".format(
+                    (nwalkers, ndim)
+                )
+            )
+        if state.log_prob.shape != (nwalkers,):
+            raise ValueError(
+                "invalid log probability size; expected {0}".format(nwalkers)
+            )
         if state.blobs is not None and not has_blobs:
             raise ValueError("unexpected blobs")
         if state.blobs is None and has_blobs:
             raise ValueError("expected blobs, but none were given")
         if state.blobs is not None and len(state.blobs) != nwalkers:
-            raise ValueError("invalid blobs size; expected {0}"
-                             .format(nwalkers))
-        if accepted.shape != (nwalkers, ):
-            raise ValueError("invalid acceptance size; expected {0}"
-                             .format(nwalkers))
+            raise ValueError(
+                "invalid blobs size; expected {0}".format(nwalkers)
+            )
+        if accepted.shape != (nwalkers,):
+            raise ValueError(
+                "invalid acceptance size; expected {0}".format(nwalkers)
+            )
 
     def save_step(self, state, accepted):
         """Save a step to the backend
diff --git a/emcee/backends/hdf.py b/src/emcee/backends/hdf.py
similarity index 81%
rename from emcee/backends/hdf.py
rename to src/emcee/backends/hdf.py
index ad244d84..27ed8ca0 100644
--- a/emcee/backends/hdf.py
+++ b/src/emcee/backends/hdf.py
@@ -1,22 +1,21 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-__all__ = ["HDFBackend", "TempHDFBackend"]
-
 import os
 from tempfile import NamedTemporaryFile
 
 import numpy as np
 
+from .. import __version__
+from .backend import Backend
+
+__all__ = ["HDFBackend", "TempHDFBackend"]
+
+
 try:
     import h5py
 except ImportError:
     h5py = None
 
-from .backend import Backend
-from .. import __version__
-
 
 class HDFBackend(Backend):
     """A backend that stores the chain in an HDF5 file using h5py
@@ -33,6 +32,7 @@ class HDFBackend(Backend):
             ``RuntimeError`` if the file is opened with write access.
 
     """
+
     def __init__(self, filename, name="mcmc", read_only=False):
         if h5py is None:
             raise ImportError("you must install 'h5py' to use the HDFBackend")
@@ -52,9 +52,11 @@ def initialized(self):
 
     def open(self, mode="r"):
         if self.read_only and mode != "r":
-            raise RuntimeError("The backend has been loaded in read-only "
-                               "mode. Set `read_only = False` to make "
-                               "changes.")
+            raise RuntimeError(
+                "The backend has been loaded in read-only "
+                "mode. Set `read_only = False` to make "
+                "changes."
+            )
         return h5py.File(self.filename, mode)
 
     def reset(self, nwalkers, ndim):
@@ -76,14 +78,18 @@ def reset(self, nwalkers, ndim):
             g.attrs["has_blobs"] = False
             g.attrs["iteration"] = 0
             g.create_dataset("accepted", data=np.zeros(nwalkers))
-            g.create_dataset("chain",
-                             (0, nwalkers, ndim),
-                             maxshape=(None, nwalkers, ndim),
-                             dtype=np.float64)
-            g.create_dataset("log_prob",
-                             (0, nwalkers),
-                             maxshape=(None, nwalkers),
-                             dtype=np.float64)
+            g.create_dataset(
+                "chain",
+                (0, nwalkers, ndim),
+                maxshape=(None, nwalkers, ndim),
+                dtype=np.float64,
+            )
+            g.create_dataset(
+                "log_prob",
+                (0, nwalkers),
+                maxshape=(None, nwalkers),
+                dtype=np.float64,
+            )
 
     def has_blobs(self):
         with self.open() as f:
@@ -91,21 +97,25 @@ def has_blobs(self):
 
     def get_value(self, name, flat=False, thin=1, discard=0):
         if not self.initialized:
-            raise AttributeError("You must run the sampler with "
-                                 "'store == True' before accessing the "
-                                 "results")
+            raise AttributeError(
+                "You must run the sampler with "
+                "'store == True' before accessing the "
+                "results"
+            )
         with self.open() as f:
             g = f[self.name]
             iteration = g.attrs["iteration"]
             if iteration <= 0:
-                raise AttributeError("You must run the sampler with "
-                                     "'store == True' before accessing the "
-                                     "results")
+                raise AttributeError(
+                    "You must run the sampler with "
+                    "'store == True' before accessing the "
+                    "results"
+                )
 
             if name == "blobs" and not g.attrs["has_blobs"]:
                 return None
 
-            v = g[name][discard+thin-1:self.iteration:thin]
+            v = g[name][discard + thin - 1 : self.iteration : thin]
             if flat:
                 s = list(v.shape[1:])
                 s[0] = np.prod(v.shape[:2])
@@ -159,9 +169,12 @@ def grow(self, ngrow, blobs):
                 if not has_blobs:
                     nwalkers = g.attrs["nwalkers"]
                     dt = np.dtype((blobs[0].dtype, blobs[0].shape))
-                    g.create_dataset("blobs", (ntot, nwalkers),
-                                     maxshape=(None, nwalkers),
-                                     dtype=dt)
+                    g.create_dataset(
+                        "blobs",
+                        (ntot, nwalkers),
+                        maxshape=(None, nwalkers),
+                        dtype=dt,
+                    )
                 else:
                     g["blobs"].resize(ntot, axis=0)
                 g.attrs["has_blobs"] = True
@@ -194,7 +207,6 @@ def save_step(self, state, accepted):
 
 
 class TempHDFBackend(object):
-
     def __enter__(self):
         f = NamedTemporaryFile("w", delete=False)
         f.close()
diff --git a/emcee/ensemble.py b/src/emcee/ensemble.py
similarity index 99%
rename from emcee/ensemble.py
rename to src/emcee/ensemble.py
index 612a8be6..9929e59f 100644
--- a/emcee/ensemble.py
+++ b/src/emcee/ensemble.py
@@ -1,6 +1,15 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
+import warnings
+
+import numpy as np
+
+from .backends import Backend
+from .model import Model
+from .moves import StretchMove
+from .pbar import get_progress_bar
+from .state import State
+from .utils import deprecated, deprecation_warning
 
 __all__ = ["EnsembleSampler"]
 
@@ -10,16 +19,6 @@
     # for py2.7, will be an Exception in 3.8
     from collections import Iterable
 
-import numpy as np
-import warnings
-
-from .state import State
-from .model import Model
-from .backends import Backend
-from .moves import StretchMove
-from .pbar import get_progress_bar
-from .utils import deprecated, deprecation_warning
-
 
 class EnsembleSampler(object):
     """An ensemble MCMC sampler
@@ -245,8 +244,8 @@ def sample(
         if np.shape(state.coords) != (self.nwalkers, self.ndim):
             raise ValueError("incompatible input dimensions")
         if (not skip_initial_state_check) and np.linalg.cond(
-             np.atleast_2d(np.cov(state.coords, rowvar=False))
-            ) > 1e8:
+            np.atleast_2d(np.cov(state.coords, rowvar=False))
+        ) > 1e8:
             warnings.warn(
                 "Initial state is not linearly independent and it will not "
                 "allow a full exploration of parameter space",
@@ -475,6 +474,7 @@ def flatchain(self):  # pragma: no cover
     def lnprobability(self):  # pragma: no cover
         log_prob = self.get_log_prob()
         return np.swapaxes(log_prob, 0, 1)
+
     @property
     @deprecated("get_log_prob(flat=True)")
     def flatlnprobability(self):  # pragma: no cover
diff --git a/emcee/interruptible_pool.py b/src/emcee/interruptible_pool.py
similarity index 77%
rename from emcee/interruptible_pool.py
rename to src/emcee/interruptible_pool.py
index 9e15dfe5..ea0b425d 100644
--- a/emcee/interruptible_pool.py
+++ b/src/emcee/interruptible_pool.py
@@ -1,8 +1,6 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-__all__ = ["InterruptiblePool"]
-
 # The standard library now has an interruptible pool
 from multiprocessing.pool import Pool as InterruptiblePool
+
+__all__ = ["InterruptiblePool"]
diff --git a/src/emcee/model.py b/src/emcee/model.py
new file mode 100644
index 00000000..8ba8eff2
--- /dev/null
+++ b/src/emcee/model.py
@@ -0,0 +1,10 @@
+# -*- coding: utf-8 -*-
+
+from collections import namedtuple
+
+__all__ = ["Model"]
+
+
+Model = namedtuple(
+    "Model", ("log_prob_fn", "compute_log_prob_fn", "map_fn", "random")
+)
diff --git a/emcee/moves/__init__.py b/src/emcee/moves/__init__.py
similarity index 64%
rename from emcee/moves/__init__.py
rename to src/emcee/moves/__init__.py
index 3125bc31..94670341 100644
--- a/emcee/moves/__init__.py
+++ b/src/emcee/moves/__init__.py
@@ -1,23 +1,23 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-from .move import Move
-
-from .mh import MHMove
+from .de import DEMove
+from .de_snooker import DESnookerMove
 from .gaussian import GaussianMove
-
+from .kde import KDEMove
+from .mh import MHMove
+from .move import Move
 from .red_blue import RedBlueMove
 from .stretch import StretchMove
 from .walk import WalkMove
-from .kde import KDEMove
-
-from .de import DEMove
-from .de_snooker import DESnookerMove
 
 __all__ = [
     "Move",
-    "MHMove", "GaussianMove",
-    "RedBlueMove", "StretchMove", "WalkMove", "KDEMove",
-    "DEMove", "DESnookerMove",
+    "MHMove",
+    "GaussianMove",
+    "RedBlueMove",
+    "StretchMove",
+    "WalkMove",
+    "KDEMove",
+    "DEMove",
+    "DESnookerMove",
 ]
diff --git a/emcee/moves/de.py b/src/emcee/moves/de.py
similarity index 97%
rename from emcee/moves/de.py
rename to src/emcee/moves/de.py
index 9fdc75b7..f0b6470d 100644
--- a/emcee/moves/de.py
+++ b/src/emcee/moves/de.py
@@ -1,8 +1,7 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
+
 from .red_blue import RedBlueMove
 
 __all__ = ["DEMove"]
diff --git a/emcee/moves/de_snooker.py b/src/emcee/moves/de_snooker.py
similarity index 92%
rename from emcee/moves/de_snooker.py
rename to src/emcee/moves/de_snooker.py
index 5cbe3ede..0b53efa4 100644
--- a/emcee/moves/de_snooker.py
+++ b/src/emcee/moves/de_snooker.py
@@ -1,8 +1,7 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
+
 from .red_blue import RedBlueMove
 
 __all__ = ["DESnookerMove"]
@@ -23,6 +22,7 @@ class DESnookerMove(RedBlueMove):
             reference.
 
     """
+
     def __init__(self, gammas=1.7, **kwargs):
         self.gammas = gammas
         kwargs["nsplits"] = 4
@@ -42,5 +42,5 @@ def get_proposal(self, s, c, random):
             norm = np.linalg.norm(delta)
             u = delta / np.sqrt(norm)
             q[i] = s[i] + u * self.gammas * (np.dot(u, z1) - np.dot(u, z2))
-            metropolis[i] = np.log(np.linalg.norm(q[i]-z)) - np.log(norm)
+            metropolis[i] = np.log(np.linalg.norm(q[i] - z)) - np.log(norm)
         return q, 0.5 * (ndim - 1.0) * metropolis
diff --git a/emcee/moves/gaussian.py b/src/emcee/moves/gaussian.py
similarity index 92%
rename from emcee/moves/gaussian.py
rename to src/emcee/moves/gaussian.py
index 22788cd9..eeafe7ef 100644
--- a/emcee/moves/gaussian.py
+++ b/src/emcee/moves/gaussian.py
@@ -1,7 +1,5 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
 
 from .mh import MHMove
@@ -32,6 +30,7 @@ class GaussianMove(MHMove):
             the other arguments are inconsistent.
 
     """
+
     def __init__(self, cov, mode="vector", factor=None):
         # Parse the proposal type.
         try:
@@ -75,9 +74,12 @@ def __init__(self, scale, factor, mode):
             self._log_factor = np.log(factor)
 
         if mode not in self.allowed_modes:
-            raise ValueError(("'{0}' is not a recognized mode. "
-                              "Please select from: {1}")
-                             .format(mode, self.allowed_modes))
+            raise ValueError(
+                (
+                    "'{0}' is not a recognized mode. "
+                    "Please select from: {1}"
+                ).format(mode, self.allowed_modes)
+            )
         self.mode = mode
 
     def get_factor(self, rng):
@@ -104,7 +106,6 @@ def __call__(self, x0, rng):
 
 
 class _diagonal_proposal(_isotropic_proposal):
-
     def get_updated_vector(self, rng, x0):
         return x0 + self.get_factor(rng) * self.scale * rng.randn(*(x0.shape))
 
@@ -115,4 +116,5 @@ class _proposal(_isotropic_proposal):
 
     def get_updated_vector(self, rng, x0):
         return x0 + self.get_factor(rng) * rng.multivariate_normal(
-            np.zeros(len(self.scale)), self.scale)
+            np.zeros(len(self.scale)), self.scale
+        )
diff --git a/emcee/moves/kde.py b/src/emcee/moves/kde.py
similarity index 86%
rename from emcee/moves/kde.py
rename to src/emcee/moves/kde.py
index a39b50de..0afddfba 100644
--- a/emcee/moves/kde.py
+++ b/src/emcee/moves/kde.py
@@ -1,15 +1,14 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
 
+from .red_blue import RedBlueMove
+
 try:
     from scipy.stats import gaussian_kde
 except ImportError:
     gaussian_kde = None
 
-from .red_blue import RedBlueMove
 
 __all__ = ["KDEMove"]
 
@@ -27,10 +26,12 @@ class KDEMove(RedBlueMove):
             for allowed values.
 
     """
+
     def __init__(self, bw_method=None, **kwargs):
         if gaussian_kde is None:
-            raise ImportError("you need scipy.stats.gaussian_kde to use the "
-                              "KDEMove")
+            raise ImportError(
+                "you need scipy.stats.gaussian_kde to use the " "KDEMove"
+            )
         self.bw_method = bw_method
         super(KDEMove, self).__init__(**kwargs)
 
diff --git a/emcee/moves/mh.py b/src/emcee/moves/mh.py
similarity index 97%
rename from emcee/moves/mh.py
rename to src/emcee/moves/mh.py
index 2d5da96a..4b190498 100644
--- a/emcee/moves/mh.py
+++ b/src/emcee/moves/mh.py
@@ -1,11 +1,9 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
 
-from .move import Move
 from ..state import State
+from .move import Move
 
 __all__ = ["MHMove"]
 
diff --git a/emcee/moves/move.py b/src/emcee/moves/move.py
similarity index 86%
rename from emcee/moves/move.py
rename to src/emcee/moves/move.py
index f77bea12..9a0663bc 100644
--- a/emcee/moves/move.py
+++ b/src/emcee/moves/move.py
@@ -1,22 +1,15 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
 
 __all__ = ["Move"]
 
 
 class Move(object):
-
     def tune(self, state, accepted):
         pass
 
-    def update(self,
-               old_state,
-               new_state,
-               accepted,
-               subset=None):
+    def update(self, old_state, new_state, accepted, subset=None):
         """Update a given subset of the ensemble with an accepted proposal
 
         Args:
@@ -45,7 +38,8 @@ def update(self,
                 raise ValueError(
                     "If you start sampling with a given log_prob, "
                     "you also need to provide the current list of "
-                    "blobs at that position.")
+                    "blobs at that position."
+                )
             old_state.blobs[m1] = new_state.blobs[m2]
 
         return old_state
diff --git a/emcee/moves/red_blue.py b/src/emcee/moves/red_blue.py
similarity index 83%
rename from emcee/moves/red_blue.py
rename to src/emcee/moves/red_blue.py
index 172b8e49..ad82fdad 100644
--- a/emcee/moves/red_blue.py
+++ b/src/emcee/moves/red_blue.py
@@ -1,11 +1,9 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
 
-from .move import Move
 from ..state import State
+from .move import Move
 
 __all__ = ["RedBlueMove"]
 
@@ -35,10 +33,10 @@ class RedBlueMove(Move):
             to @dstndstn for this wonderful terminology.
 
     """
-    def __init__(self,
-                 nsplits=2,
-                 randomize_split=True,
-                 live_dangerously=False):
+
+    def __init__(
+        self, nsplits=2, randomize_split=True, live_dangerously=False
+    ):
         self.nsplits = int(nsplits)
         self.live_dangerously = live_dangerously
         self.randomize_split = randomize_split
@@ -47,8 +45,9 @@ def setup(self, coords):
         pass
 
     def get_proposal(self, sample, complement, random):
-        raise NotImplementedError("The proposal must be implemented by "
-                                  "subclasses")
+        raise NotImplementedError(
+            "The proposal must be implemented by " "subclasses"
+        )
 
     def propose(self, model, state):
         """Use the move to generate a proposal and compute the acceptance
@@ -64,9 +63,11 @@ def propose(self, model, state):
         # Check that the dimensions are compatible.
         nwalkers, ndim = state.coords.shape
         if nwalkers < 2 * ndim and not self.live_dangerously:
-            raise RuntimeError("It is unadvisable to use a red-blue move "
-                               "with fewer walkers than twice the number of "
-                               "dimensions.")
+            raise RuntimeError(
+                "It is unadvisable to use a red-blue move "
+                "with fewer walkers than twice the number of "
+                "dimensions."
+            )
 
         # Run any move-specific setup.
         self.setup(state.coords)
@@ -83,7 +84,7 @@ def propose(self, model, state):
             # Get the two halves of the ensemble.
             sets = [state.coords[inds == j] for j in range(self.nsplits)]
             s = sets[split]
-            c = sets[:split] + sets[split+1:]
+            c = sets[:split] + sets[split + 1 :]
 
             # Get the move-specific proposal.
             q, factors = self.get_proposal(s, c, model.random)
@@ -92,8 +93,9 @@ def propose(self, model, state):
             new_log_probs, new_blobs = model.compute_log_prob_fn(q)
 
             # Loop over the walkers and update them accordingly.
-            for i, (j, f, nlp) in enumerate(zip(
-                    all_inds[S1], factors, new_log_probs)):
+            for i, (j, f, nlp) in enumerate(
+                zip(all_inds[S1], factors, new_log_probs)
+            ):
                 lnpdiff = f + nlp - state.log_prob[j]
                 if lnpdiff > np.log(model.random.rand()):
                     accepted[j] = True
diff --git a/emcee/moves/stretch.py b/src/emcee/moves/stretch.py
similarity index 83%
rename from emcee/moves/stretch.py
rename to src/emcee/moves/stretch.py
index 865f485b..a6aed5c4 100644
--- a/emcee/moves/stretch.py
+++ b/src/emcee/moves/stretch.py
@@ -1,7 +1,5 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
 
 from .red_blue import RedBlueMove
@@ -20,6 +18,7 @@ class StretchMove(RedBlueMove):
         The stretch scale parameter. (default: ``2.0``)
 
     """
+
     def __init__(self, a=2.0, **kwargs):
         self.a = a
         super(StretchMove, self).__init__(**kwargs)
@@ -28,7 +27,7 @@ def get_proposal(self, s, c, random):
         c = np.concatenate(c, axis=0)
         Ns, Nc = len(s), len(c)
         ndim = s.shape[1]
-        zz = ((self.a - 1.) * random.rand(Ns) + 1) ** 2. / self.a
-        factors = (ndim - 1.) * np.log(zz)
+        zz = ((self.a - 1.0) * random.rand(Ns) + 1) ** 2.0 / self.a
+        factors = (ndim - 1.0) * np.log(zz)
         rint = random.randint(Nc, size=(Ns,))
         return c[rint] - (c[rint] - s) * zz[:, None], factors
diff --git a/emcee/moves/walk.py b/src/emcee/moves/walk.py
similarity index 95%
rename from emcee/moves/walk.py
rename to src/emcee/moves/walk.py
index f5323703..a0913f62 100644
--- a/emcee/moves/walk.py
+++ b/src/emcee/moves/walk.py
@@ -1,8 +1,7 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
+
 from .red_blue import RedBlueMove
 
 __all__ = ["WalkMove"]
@@ -20,6 +19,7 @@ class WalkMove(RedBlueMove):
         walkers in the complement.
 
     """
+
     def __init__(self, s=None, **kwargs):
         self.s = s
         super(WalkMove, self).__init__(**kwargs)
diff --git a/emcee/mpi_pool.py b/src/emcee/mpi_pool.py
similarity index 88%
rename from emcee/mpi_pool.py
rename to src/emcee/mpi_pool.py
index a8c7fea5..be728c63 100644
--- a/emcee/mpi_pool.py
+++ b/src/emcee/mpi_pool.py
@@ -1,13 +1,11 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import print_function, absolute_import
 
 try:
     from schwimmbad import MPIPool
 except ImportError:
 
     class MPIPool(object):
-
         def __init__(self, *args, **kwargs):
             raise ImportError(
                 "The MPIPool from emcee has been forked to "
@@ -15,4 +13,5 @@ def __init__(self, *args, **kwargs):
                 "please install that package to continue using the MPIPool"
             )
 
+
 __all__ = ["MPIPool"]
diff --git a/emcee/pbar.py b/src/emcee/pbar.py
similarity index 96%
rename from emcee/pbar.py
rename to src/emcee/pbar.py
index 99ef5a54..07085aa4 100644
--- a/emcee/pbar.py
+++ b/src/emcee/pbar.py
@@ -1,10 +1,9 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
+import logging
 
 __all__ = ["get_progress_bar"]
 
-import logging
 
 try:
     import tqdm
@@ -54,4 +53,3 @@ def get_progress_bar(display, total):
                 return getattr(tqdm, "tqdm_" + display)(total=total)
 
     return _NoOpPBar()
-
diff --git a/emcee/ptsampler.py b/src/emcee/ptsampler.py
similarity index 88%
rename from emcee/ptsampler.py
rename to src/emcee/ptsampler.py
index ac4638f2..05040920 100644
--- a/emcee/ptsampler.py
+++ b/src/emcee/ptsampler.py
@@ -1,13 +1,11 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import print_function, absolute_import
 
 try:
     from ptemcee import Sampler as PTSampler
 except ImportError:
 
     class PTSampler(object):
-
         def __init__(self, *args, **kwargs):
             raise ImportError(
                 "The PTSampler from emcee has been forked to "
@@ -15,4 +13,5 @@ def __init__(self, *args, **kwargs):
                 "please install that package to continue using the PTSampler"
             )
 
+
 __all__ = ["PTSampler"]
diff --git a/emcee/state.py b/src/emcee/state.py
similarity index 87%
rename from emcee/state.py
rename to src/emcee/state.py
index 91ace70b..e6555f3c 100644
--- a/emcee/state.py
+++ b/src/emcee/state.py
@@ -1,11 +1,10 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-__all__ = ["State"]
+from copy import deepcopy
 
 import numpy as np
-from copy import deepcopy
+
+__all__ = ["State"]
 
 
 class State(object):
@@ -28,8 +27,9 @@ class State(object):
 
     __slots__ = "coords", "log_prob", "blobs", "random_state"
 
-    def __init__(self, coords, log_prob=None, blobs=None, random_state=None,
-                 copy=False):
+    def __init__(
+        self, coords, log_prob=None, blobs=None, random_state=None, copy=False
+    ):
         dc = deepcopy if copy else lambda x: x
 
         if hasattr(coords, "coords"):
@@ -52,5 +52,6 @@ def __repr__(self):
     def __iter__(self):
         if self.blobs is None:
             return iter((self.coords, self.log_prob, self.random_state))
-        return iter((self.coords, self.log_prob, self.random_state,
-                     self.blobs))
+        return iter(
+            (self.coords, self.log_prob, self.random_state, self.blobs)
+        )
diff --git a/emcee/tests/__init__.py b/src/emcee/tests/__init__.py
similarity index 100%
rename from emcee/tests/__init__.py
rename to src/emcee/tests/__init__.py
diff --git a/emcee/tests/integration/__init__.py b/src/emcee/tests/integration/__init__.py
similarity index 100%
rename from emcee/tests/integration/__init__.py
rename to src/emcee/tests/integration/__init__.py
diff --git a/emcee/tests/integration/test_de.py b/src/emcee/tests/integration/test_de.py
similarity index 89%
rename from emcee/tests/integration/test_de.py
rename to src/emcee/tests/integration/test_de.py
index 6ce2f0b1..d9c8ec68 100644
--- a/emcee/tests/integration/test_de.py
+++ b/src/emcee/tests/integration/test_de.py
@@ -1,8 +1,7 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 from emcee import moves
+
 from .test_proposal import _test_normal, _test_uniform
 
 __all__ = ["test_normal_de", "test_normal_de_no_gamma", "test_uniform_de"]
diff --git a/emcee/tests/integration/test_de_snooker.py b/src/emcee/tests/integration/test_de_snooker.py
similarity index 88%
rename from emcee/tests/integration/test_de_snooker.py
rename to src/emcee/tests/integration/test_de_snooker.py
index d4afc87f..bdf2d6fe 100644
--- a/emcee/tests/integration/test_de_snooker.py
+++ b/src/emcee/tests/integration/test_de_snooker.py
@@ -1,8 +1,7 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 from emcee import moves
+
 from .test_proposal import _test_normal, _test_uniform
 
 __all__ = ["test_normal_de_snooker", "test_uniform_de_snooker"]
diff --git a/src/emcee/tests/integration/test_gaussian.py b/src/emcee/tests/integration/test_gaussian.py
new file mode 100644
index 00000000..d305d341
--- /dev/null
+++ b/src/emcee/tests/integration/test_gaussian.py
@@ -0,0 +1,79 @@
+# -*- coding: utf-8 -*-
+
+from itertools import product
+
+import numpy as np
+import pytest
+
+from emcee import moves
+
+from .test_proposal import _test_normal, _test_uniform
+
+__all__ = [
+    "test_normal_gaussian",
+    "test_uniform_gaussian",
+    "test_normal_gaussian_nd",
+]
+
+
+@pytest.mark.parametrize("mode,factor", product(["vector"], [None, 2.0, 5.0]))
+def test_normal_gaussian(mode, factor, **kwargs):
+    _test_normal(moves.GaussianMove(0.5, mode=mode, factor=factor), **kwargs)
+
+
+@pytest.mark.parametrize(
+    "mode,factor", product(["vector", "random", "sequential"], [None, 2.0])
+)
+def test_normal_gaussian_nd(mode, factor, **kwargs):
+    ndim = 3
+    kwargs["nsteps"] = 8000
+
+    # Isotropic.
+    _test_normal(
+        moves.GaussianMove(0.5, factor=factor, mode=mode), ndim=ndim, **kwargs
+    )
+
+    # Axis-aligned.
+    _test_normal(
+        moves.GaussianMove(0.5 * np.ones(ndim), factor=factor, mode=mode),
+        ndim=ndim,
+        **kwargs
+    )
+    with pytest.raises(ValueError):
+        _test_normal(
+            moves.GaussianMove(
+                0.5 * np.ones(ndim - 1), factor=factor, mode=mode
+            ),
+            ndim=ndim,
+            **kwargs
+        )
+
+    # Full matrix.
+    if mode == "vector":
+        _test_normal(
+            moves.GaussianMove(
+                np.diag(0.5 * np.ones(ndim)), factor=factor, mode=mode
+            ),
+            ndim=ndim,
+            **kwargs
+        )
+        with pytest.raises(ValueError):
+            _test_normal(
+                moves.GaussianMove(np.diag(0.5 * np.ones(ndim - 1))),
+                ndim=ndim,
+                **kwargs
+            )
+    else:
+        with pytest.raises(ValueError):
+            _test_normal(
+                moves.GaussianMove(
+                    np.diag(0.5 * np.ones(ndim)), factor=factor, mode=mode
+                ),
+                ndim=ndim,
+                **kwargs
+            )
+
+
+@pytest.mark.parametrize("mode,factor", product(["vector"], [None, 2.0, 5.0]))
+def test_uniform_gaussian(mode, factor, **kwargs):
+    _test_uniform(moves.GaussianMove(0.5, factor=factor, mode=mode), **kwargs)
diff --git a/emcee/tests/integration/test_kde.py b/src/emcee/tests/integration/test_kde.py
similarity index 89%
rename from emcee/tests/integration/test_kde.py
rename to src/emcee/tests/integration/test_kde.py
index 647f8aa5..f2142fc2 100644
--- a/emcee/tests/integration/test_kde.py
+++ b/src/emcee/tests/integration/test_kde.py
@@ -1,8 +1,7 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 from emcee import moves
+
 from .test_proposal import _test_normal, _test_uniform
 
 __all__ = ["test_normal_kde", "test_uniform_kde", "test_nsplits_kde"]
diff --git a/emcee/tests/integration/test_proposal.py b/src/emcee/tests/integration/test_proposal.py
similarity index 72%
rename from emcee/tests/integration/test_proposal.py
rename to src/emcee/tests/integration/test_proposal.py
index 2c759a71..87d67444 100644
--- a/emcee/tests/integration/test_proposal.py
+++ b/src/emcee/tests/integration/test_proposal.py
@@ -1,7 +1,5 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
 from scipy import stats
 
@@ -11,11 +9,11 @@
 
 
 def normal_log_prob_blobs(params):
-    return -0.5 * np.sum(params**2), params
+    return -0.5 * np.sum(params ** 2), params
 
 
 def normal_log_prob(params):
-    return -0.5 * np.sum(params**2)
+    return -0.5 * np.sum(params ** 2)
 
 
 def uniform_log_prob(params):
@@ -24,8 +22,16 @@ def uniform_log_prob(params):
     return 0.0
 
 
-def _test_normal(proposal, ndim=1, nwalkers=32, nsteps=2000, seed=1234,
-                 check_acceptance=True, pool=None, blobs=False):
+def _test_normal(
+    proposal,
+    ndim=1,
+    nwalkers=32,
+    nsteps=2000,
+    seed=1234,
+    check_acceptance=True,
+    pool=None,
+    blobs=False,
+):
     # Set up the random number generator.
     np.random.seed(seed)
 
@@ -37,8 +43,9 @@ def _test_normal(proposal, ndim=1, nwalkers=32, nsteps=2000, seed=1234,
     else:
         lp = normal_log_prob
 
-    sampler = emcee.EnsembleSampler(nwalkers, ndim, lp,
-                                    moves=proposal, pool=pool)
+    sampler = emcee.EnsembleSampler(
+        nwalkers, ndim, lp, moves=proposal, pool=pool
+    )
     if hasattr(proposal, "ntune") and proposal.ntune > 0:
         coords = sampler.run_mcmc(coords, proposal.ntune, tune=True)
         sampler.reset()
@@ -47,8 +54,9 @@ def _test_normal(proposal, ndim=1, nwalkers=32, nsteps=2000, seed=1234,
     # Check the acceptance fraction.
     if check_acceptance:
         acc = sampler.acceptance_fraction
-        assert np.all((acc < 0.9) * (acc > 0.1)), \
-            "Invalid acceptance fraction\n{0}".format(acc)
+        assert np.all(
+            (acc < 0.9) * (acc > 0.1)
+        ), "Invalid acceptance fraction\n{0}".format(acc)
 
     # Check the resulting chain using a K-S test and compare to the mean and
     # standard deviation.
@@ -69,14 +77,16 @@ def _test_uniform(proposal, nwalkers=32, nsteps=2000, seed=1234):
     # Initialize the ensemble and proposal.
     coords = np.random.rand(nwalkers, 1)
 
-    sampler = emcee.EnsembleSampler(nwalkers, 1, normal_log_prob,
-                                    moves=proposal)
+    sampler = emcee.EnsembleSampler(
+        nwalkers, 1, normal_log_prob, moves=proposal
+    )
     sampler.run_mcmc(coords, nsteps)
 
     # Check the acceptance fraction.
     acc = sampler.acceptance_fraction
-    assert np.all((acc < 0.9) * (acc > 0.1)), \
-        "Invalid acceptance fraction\n{0}".format(acc)
+    assert np.all(
+        (acc < 0.9) * (acc > 0.1)
+    ), "Invalid acceptance fraction\n{0}".format(acc)
 
     # Check that the resulting chain "fails" the K-S test.
     samps = sampler.get_chain(flat=True)
diff --git a/emcee/tests/integration/test_stretch.py b/src/emcee/tests/integration/test_stretch.py
similarity index 80%
rename from emcee/tests/integration/test_stretch.py
rename to src/emcee/tests/integration/test_stretch.py
index d09289f5..16075a4e 100644
--- a/emcee/tests/integration/test_stretch.py
+++ b/src/emcee/tests/integration/test_stretch.py
@@ -1,13 +1,16 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import pytest
+
 from emcee import moves
+
 from .test_proposal import _test_normal, _test_uniform
 
-__all__ = ["test_normal_stretch", "test_uniform_stretch",
-           "test_nsplits_stretch"]
+__all__ = [
+    "test_normal_stretch",
+    "test_uniform_stretch",
+    "test_nsplits_stretch",
+]
 
 
 @pytest.mark.parametrize("blobs", [True, False])
diff --git a/emcee/tests/integration/test_walk.py b/src/emcee/tests/integration/test_walk.py
similarity index 87%
rename from emcee/tests/integration/test_walk.py
rename to src/emcee/tests/integration/test_walk.py
index d09888ec..8879a62d 100644
--- a/emcee/tests/integration/test_walk.py
+++ b/src/emcee/tests/integration/test_walk.py
@@ -1,8 +1,7 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 from emcee import moves
+
 from .test_proposal import _test_normal, _test_uniform
 
 __all__ = ["test_normal_walk", "test_uniform_walk"]
diff --git a/emcee/tests/unit/__init__.py b/src/emcee/tests/unit/__init__.py
similarity index 100%
rename from emcee/tests/unit/__init__.py
rename to src/emcee/tests/unit/__init__.py
diff --git a/emcee/tests/unit/test_autocorr.py b/src/emcee/tests/unit/test_autocorr.py
similarity index 71%
rename from emcee/tests/unit/test_autocorr.py
rename to src/emcee/tests/unit/test_autocorr.py
index f78dc2ce..9d4365b5 100644
--- a/emcee/tests/unit/test_autocorr.py
+++ b/src/emcee/tests/unit/test_autocorr.py
@@ -1,11 +1,9 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-import pytest
 import numpy as np
+import pytest
 
-from emcee.autocorr import integrated_time, AutocorrError
+from emcee.autocorr import AutocorrError, integrated_time
 
 
 def get_chain(seed=1234, ndim=3, N=100000):
@@ -14,20 +12,20 @@ def get_chain(seed=1234, ndim=3, N=100000):
     x = np.empty((N, ndim))
     x[0] = np.zeros(ndim)
     for i in range(1, N):
-        x[i] = x[i-1] * a + np.random.rand(ndim)
+        x[i] = x[i - 1] * a + np.random.rand(ndim)
     return x
 
 
 def test_1d(seed=1234, ndim=1, N=250000):
     x = get_chain(seed=seed, ndim=ndim, N=N)
     tau = integrated_time(x)
-    assert np.all(np.abs(tau - 19.0) / 19. < 0.2)
+    assert np.all(np.abs(tau - 19.0) / 19.0 < 0.2)
 
 
 def test_nd(seed=1234, ndim=3, N=150000):
     x = get_chain(seed=seed, ndim=ndim, N=N)
     tau = integrated_time(x)
-    assert np.all(np.abs(tau - 19.0) / 19. < 0.2)
+    assert np.all(np.abs(tau - 19.0) / 19.0 < 0.2)
 
 
 def test_too_short(seed=1234, ndim=3, N=100):
@@ -43,7 +41,8 @@ def test_autocorr_multi_works():
 
     # This throws exception unconditionally in buggy impl's
     acls_multi = integrated_time(xs)
-    acls_single = np.array([integrated_time(xs[:, i])
-                            for i in range(xs.shape[1])])
+    acls_single = np.array(
+        [integrated_time(xs[:, i]) for i in range(xs.shape[1])]
+    )
 
     assert np.all(np.abs(acls_multi - acls_single) < 2)
diff --git a/emcee/tests/unit/test_backends.py b/src/emcee/tests/unit/test_backends.py
similarity index 78%
rename from emcee/tests/unit/test_backends.py
rename to src/emcee/tests/unit/test_backends.py
index 392aee6e..8e0ff193 100644
--- a/emcee/tests/unit/test_backends.py
+++ b/src/emcee/tests/unit/test_backends.py
@@ -1,45 +1,48 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import os
 from itertools import product
 
+import h5py
 import numpy as np
-
 import pytest
 
-from emcee import backends, EnsembleSampler
-
-import h5py
+from emcee import EnsembleSampler, backends
 
 __all__ = ["test_backend", "test_reload"]
 
 all_backends = backends.get_test_backends()
 other_backends = all_backends[1:]
-dtypes = [
-    None,
-    [("log_prior", float), ("mean", int)]
-]
+dtypes = [None, [("log_prior", float), ("mean", int)]]
 
 
 def normal_log_prob(params):
-    return -0.5 * np.sum(params**2)
+    return -0.5 * np.sum(params ** 2)
 
 
 def normal_log_prob_blobs(params):
     return normal_log_prob(params), 0.1, int(5)
 
 
-def run_sampler(backend, nwalkers=32, ndim=3, nsteps=25, seed=1234, thin_by=1,
-                dtype=None, blobs=True, lp=None):
+def run_sampler(
+    backend,
+    nwalkers=32,
+    ndim=3,
+    nsteps=25,
+    seed=1234,
+    thin_by=1,
+    dtype=None,
+    blobs=True,
+    lp=None,
+):
     if lp is None:
         lp = normal_log_prob_blobs if blobs else normal_log_prob
     if seed is not None:
         np.random.seed(seed)
     coords = np.random.randn(nwalkers, ndim)
-    sampler = EnsembleSampler(nwalkers, ndim, lp,
-                              backend=backend, blobs_dtype=dtype)
+    sampler = EnsembleSampler(
+        nwalkers, ndim, lp, backend=backend, blobs_dtype=dtype
+    )
     sampler.run_mcmc(coords, nsteps, thin_by=thin_by)
     return sampler
 
@@ -88,8 +91,9 @@ def test_blob_usage_errors(backend):
             run_sampler(be, blobs=True)
 
 
-@pytest.mark.parametrize("backend,dtype,blobs",
-                         product(other_backends, dtypes, [True, False]))
+@pytest.mark.parametrize(
+    "backend,dtype,blobs", product(other_backends, dtypes, [True, False])
+)
 def test_backend(backend, dtype, blobs):
     # Run a sampler with the default backend.
     sampler1 = run_sampler(backends.Backend(), dtype=dtype, blobs=blobs)
@@ -114,9 +118,10 @@ def test_backend(backend, dtype, blobs):
         last2 = sampler2.get_last_sample()
         assert np.allclose(last1.coords, last2.coords)
         assert np.allclose(last1.log_prob, last2.log_prob)
-        assert all(np.allclose(l1, l2)
-                   for l1, l2 in zip(last1.random_state[1:],
-                                     last2.random_state[1:]))
+        assert all(
+            np.allclose(l1, l2)
+            for l1, l2 in zip(last1.random_state[1:], last2.random_state[1:])
+        )
         if blobs:
             _custom_allclose(last1.blobs, last2.blobs)
         else:
@@ -137,15 +142,18 @@ def test_reload(backend, dtype):
         np.random.set_state(state)
 
         # Load the file using a new backend object.
-        backend2 = backends.HDFBackend(backend1.filename, backend1.name,
-                                       read_only=True)
+        backend2 = backends.HDFBackend(
+            backend1.filename, backend1.name, read_only=True
+        )
 
         with pytest.raises(RuntimeError):
             backend2.reset(32, 3)
 
         assert state[0] == backend2.random_state[0]
-        assert all(np.allclose(a, b)
-                   for a, b in zip(state[1:], backend2.random_state[1:]))
+        assert all(
+            np.allclose(a, b)
+            for a, b in zip(state[1:], backend2.random_state[1:])
+        )
 
         # Check all of the components.
         for k in ["chain", "log_prob", "blobs"]:
@@ -157,9 +165,10 @@ def test_reload(backend, dtype):
         last2 = backend2.get_last_sample()
         assert np.allclose(last1.coords, last2.coords)
         assert np.allclose(last1.log_prob, last2.log_prob)
-        assert all(np.allclose(l1, l2)
-                   for l1, l2 in zip(last1.random_state[1:],
-                                     last2.random_state[1:]))
+        assert all(
+            np.allclose(l1, l2)
+            for l1, l2 in zip(last1.random_state[1:], last2.random_state[1:])
+        )
         _custom_allclose(last1.blobs, last2.blobs)
 
         a = backend1.accepted
@@ -188,9 +197,10 @@ def test_restart(backend, dtype):
         last2 = sampler2.get_last_sample()
         assert np.allclose(last1.coords, last2.coords)
         assert np.allclose(last1.log_prob, last2.log_prob)
-        assert all(np.allclose(l1, l2)
-                   for l1, l2 in zip(last1.random_state[1:],
-                                     last2.random_state[1:]))
+        assert all(
+            np.allclose(l1, l2)
+            for l1, l2 in zip(last1.random_state[1:], last2.random_state[1:])
+        )
         _custom_allclose(last1.blobs, last2.blobs)
 
         a = sampler1.acceptance_fraction
@@ -202,12 +212,12 @@ def test_multi_hdf5():
     with backends.TempHDFBackend() as backend1:
         run_sampler(backend1)
 
-        backend2 = backends.HDFBackend(backend1.filename, name='mcmc2')
+        backend2 = backends.HDFBackend(backend1.filename, name="mcmc2")
         run_sampler(backend2)
         chain2 = backend2.get_chain()
 
-        with h5py.File(backend1.filename, 'r') as f:
-            assert set(f.keys()) == {backend1.name, 'mcmc2'}
+        with h5py.File(backend1.filename, "r") as f:
+            assert set(f.keys()) == {backend1.name, "mcmc2"}
 
         backend1.reset(10, 2)
         assert np.allclose(backend2.get_chain(), chain2)
diff --git a/emcee/tests/unit/test_blobs.py b/src/emcee/tests/unit/test_blobs.py
similarity index 63%
rename from emcee/tests/unit/test_blobs.py
rename to src/emcee/tests/unit/test_blobs.py
index 1f3695f3..d069f327 100644
--- a/emcee/tests/unit/test_blobs.py
+++ b/src/emcee/tests/unit/test_blobs.py
@@ -1,35 +1,35 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-import pytest
 import numpy as np
+import pytest
 
-from emcee import backends, EnsembleSampler
+from emcee import EnsembleSampler, backends
 
 __all__ = ["test_blob_shape"]
 
 
 class BlobLogProb(object):
-
     def __init__(self, blob_function):
         self.blob_function = blob_function
 
     def __call__(self, params):
-        return -0.5 * np.sum(params**2), self.blob_function(params)
+        return -0.5 * np.sum(params ** 2), self.blob_function(params)
 
 
 @pytest.mark.parametrize("backend", backends.get_test_backends())
-@pytest.mark.parametrize("blob_spec", [
-    (True, 5, lambda x: np.random.randn(5)),
-    (True, 0, lambda x: np.random.randn()),
-    (False, 2, lambda x: (1.0, np.random.randn(3))),
-    (False, 0, lambda x: "face"),
-    (False, 2, lambda x: (np.random.randn(5), "face")),
-])
+@pytest.mark.parametrize(
+    "blob_spec",
+    [
+        (True, 5, lambda x: np.random.randn(5)),
+        (True, 0, lambda x: np.random.randn()),
+        (False, 2, lambda x: (1.0, np.random.randn(3))),
+        (False, 0, lambda x: "face"),
+        (False, 2, lambda x: (np.random.randn(5), "face")),
+    ],
+)
 def test_blob_shape(backend, blob_spec):
     # HDF backends don't support the object type
-    if backend in (backends.TempHDFBackend, ) and not blob_spec[0]:
+    if backend in (backends.TempHDFBackend,) and not blob_spec[0]:
         return
 
     with backend() as be:
diff --git a/emcee/tests/unit/test_sampler.py b/src/emcee/tests/unit/test_sampler.py
similarity index 65%
rename from emcee/tests/unit/test_sampler.py
rename to src/emcee/tests/unit/test_sampler.py
index ca38a569..6b4cdd30 100644
--- a/emcee/tests/unit/test_sampler.py
+++ b/src/emcee/tests/unit/test_sampler.py
@@ -1,14 +1,12 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import pickle
 from itertools import product
 
-import pytest
 import numpy as np
+import pytest
 
-from emcee import moves, backends, EnsembleSampler
+from emcee import EnsembleSampler, backends, moves
 
 __all__ = ["test_shapes", "test_errors", "test_thin", "test_vectorize"]
 
@@ -16,17 +14,20 @@
 
 
 def normal_log_prob(params):
-    return -0.5 * np.sum(params**2)
-
-
-@pytest.mark.parametrize("backend, moves", product(
-    all_backends,
-    [
-        None,
-        moves.GaussianMove(0.5),
-        [moves.StretchMove(), moves.GaussianMove(0.5)],
-        [(moves.StretchMove(), 0.3), (moves.GaussianMove(0.5), 0.1)],
-    ])
+    return -0.5 * np.sum(params ** 2)
+
+
+@pytest.mark.parametrize(
+    "backend, moves",
+    product(
+        all_backends,
+        [
+            None,
+            moves.GaussianMove(0.5),
+            [moves.StretchMove(), moves.GaussianMove(0.5)],
+            [(moves.StretchMove(), 0.3), (moves.GaussianMove(0.5), 0.1)],
+        ],
+    ),
 )
 def test_shapes(backend, moves, nwalkers=32, ndim=3, nsteps=10, seed=1234):
     # Set up the random number generator.
@@ -35,8 +36,9 @@ def test_shapes(backend, moves, nwalkers=32, ndim=3, nsteps=10, seed=1234):
     with backend() as be:
         # Initialize the ensemble, moves and sampler.
         coords = np.random.randn(nwalkers, ndim)
-        sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob,
-                                  moves=moves, backend=be)
+        sampler = EnsembleSampler(
+            nwalkers, ndim, normal_log_prob, moves=moves, backend=be
+        )
 
         # Run the sampler.
         sampler.run_mcmc(coords, nsteps)
@@ -48,24 +50,38 @@ def test_shapes(backend, moves, nwalkers=32, ndim=3, nsteps=10, seed=1234):
 
         # Check the shapes.
         with pytest.warns(DeprecationWarning):
-            assert sampler.chain.shape == (nwalkers, nsteps, ndim), \
-                "incorrect coordinate dimensions"
+            assert sampler.chain.shape == (
+                nwalkers,
+                nsteps,
+                ndim,
+            ), "incorrect coordinate dimensions"
         with pytest.warns(DeprecationWarning):
-            assert sampler.lnprobability.shape == (nwalkers, nsteps), \
-                "incorrect probability dimensions"
-        assert sampler.get_chain().shape == (nsteps, nwalkers, ndim), \
-            "incorrect coordinate dimensions"
-        assert sampler.get_log_prob().shape == (nsteps, nwalkers), \
-            "incorrect probability dimensions"
-
-        assert sampler.acceptance_fraction.shape == (nwalkers,), \
-            "incorrect acceptance fraction dimensions"
+            assert sampler.lnprobability.shape == (
+                nwalkers,
+                nsteps,
+            ), "incorrect probability dimensions"
+        assert sampler.get_chain().shape == (
+            nsteps,
+            nwalkers,
+            ndim,
+        ), "incorrect coordinate dimensions"
+        assert sampler.get_log_prob().shape == (
+            nsteps,
+            nwalkers,
+        ), "incorrect probability dimensions"
+
+        assert sampler.acceptance_fraction.shape == (
+            nwalkers,
+        ), "incorrect acceptance fraction dimensions"
 
         # Check the shape of the flattened coords.
-        assert sampler.get_chain(flat=True).shape == \
-            (nsteps * nwalkers, ndim), "incorrect coordinate dimensions"
-        assert sampler.get_log_prob(flat=True).shape == \
-            (nsteps*nwalkers,), "incorrect probability dimensions"
+        assert sampler.get_chain(flat=True).shape == (
+            nsteps * nwalkers,
+            ndim,
+        ), "incorrect coordinate dimensions"
+        assert sampler.get_log_prob(flat=True).shape == (
+            nsteps * nwalkers,
+        ), "incorrect probability dimensions"
 
 
 @pytest.mark.parametrize("backend", all_backends)
@@ -76,8 +92,7 @@ def test_errors(backend, nwalkers=32, ndim=3, nsteps=5, seed=1234):
     with backend() as be:
         # Initialize the ensemble, proposal, and sampler.
         coords = np.random.randn(nwalkers, ndim)
-        sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob,
-                                  backend=be)
+        sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob, backend=be)
 
         # Test for not running.
         with pytest.raises(AttributeError):
@@ -94,7 +109,7 @@ def test_errors(backend, nwalkers=32, ndim=3, nsteps=5, seed=1234):
         # ensemble of a different shape.
         sampler.run_mcmc(coords, nsteps, store=False)
 
-        coords2 = np.random.randn(nwalkers, ndim+1)
+        coords2 = np.random.randn(nwalkers, ndim + 1)
         with pytest.raises(ValueError):
             list(sampler.run_mcmc(coords2, nsteps))
 
@@ -105,19 +120,37 @@ def test_errors(backend, nwalkers=32, ndim=3, nsteps=5, seed=1234):
             sampler.run_mcmc(np.ones((nwalkers, ndim)), nsteps)
             assert len(recorded_warnings) == 1
         with pytest.warns(None) as recorded_warnings:
-            sampler.run_mcmc(np.ones((nwalkers, ndim)),
-                    nsteps, skip_initial_state_check=True)
+            sampler.run_mcmc(
+                np.ones((nwalkers, ndim)),
+                nsteps,
+                skip_initial_state_check=True,
+            )
             sampler.run_mcmc(np.random.randn(nwalkers, ndim), nsteps)
             assert len(recorded_warnings) == 0
 
-def run_sampler(backend, nwalkers=32, ndim=3, nsteps=25, seed=1234,
-                thin=None, thin_by=1, progress=False, store=True):
+
+def run_sampler(
+    backend,
+    nwalkers=32,
+    ndim=3,
+    nsteps=25,
+    seed=1234,
+    thin=None,
+    thin_by=1,
+    progress=False,
+    store=True,
+):
     np.random.seed(seed)
     coords = np.random.randn(nwalkers, ndim)
-    sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob,
-                              backend=backend)
-    sampler.run_mcmc(coords, nsteps, thin=thin, thin_by=thin_by,
-                     progress=progress, store=store)
+    sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob, backend=backend)
+    sampler.run_mcmc(
+        coords,
+        nsteps,
+        thin=thin,
+        thin_by=thin_by,
+        progress=progress,
+        store=store,
+    )
     return sampler
 
 
@@ -135,15 +168,16 @@ def test_thin(backend):
         with pytest.warns(DeprecationWarning):
             sampler2 = run_sampler(be, thin=thinby)
         for k in ["get_chain", "get_log_prob"]:
-            a = getattr(sampler1, k)()[thinby-1::thinby]
+            a = getattr(sampler1, k)()[thinby - 1 :: thinby]
             b = getattr(sampler2, k)()
             c = getattr(sampler1, k)(thin=thinby)
             assert np.allclose(a, b), "inconsistent {0}".format(k)
             assert np.allclose(a, c), "inconsistent {0}".format(k)
 
 
-@pytest.mark.parametrize("backend,progress",
-                         product(all_backends, [True, False]))
+@pytest.mark.parametrize(
+    "backend,progress", product(all_backends, [True, False])
+)
 def test_thin_by(backend, progress):
     with backend() as be:
         with pytest.raises(ValueError):
@@ -152,16 +186,17 @@ def test_thin_by(backend, progress):
             run_sampler(be, thin_by=0.1)
         nsteps = 25
         thinby = 3
-        sampler1 = run_sampler(None, nsteps=nsteps*thinby, progress=progress)
-        sampler2 = run_sampler(be, thin_by=thinby, progress=progress,
-                               nsteps=nsteps)
+        sampler1 = run_sampler(None, nsteps=nsteps * thinby, progress=progress)
+        sampler2 = run_sampler(
+            be, thin_by=thinby, progress=progress, nsteps=nsteps
+        )
         for k in ["get_chain", "get_log_prob"]:
-            a = getattr(sampler1, k)()[thinby-1::thinby]
+            a = getattr(sampler1, k)()[thinby - 1 :: thinby]
             b = getattr(sampler2, k)()
             c = getattr(sampler1, k)(thin=thinby)
             assert np.allclose(a, b), "inconsistent {0}".format(k)
             assert np.allclose(a, c), "inconsistent {0}".format(k)
-        assert sampler1.iteration == sampler2.iteration*thinby
+        assert sampler1.iteration == sampler2.iteration * thinby
 
 
 @pytest.mark.parametrize("backend", all_backends)
@@ -181,7 +216,7 @@ def test_restart(backend):
 
 def test_vectorize():
     def lp_vec(p):
-        return -0.5 * np.sum(p**2, axis=1)
+        return -0.5 * np.sum(p ** 2, axis=1)
 
     np.random.seed(42)
     nwalkers, ndim = 32, 3
diff --git a/emcee/tests/unit/test_state.py b/src/emcee/tests/unit/test_state.py
similarity index 92%
rename from emcee/tests/unit/test_state.py
rename to src/emcee/tests/unit/test_state.py
index 22d678aa..44c540e8 100644
--- a/emcee/tests/unit/test_state.py
+++ b/src/emcee/tests/unit/test_state.py
@@ -1,10 +1,9 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import numpy as np
-from emcee.state import State
+
 from emcee import EnsembleSampler
+from emcee.state import State
 
 
 def test_back_compat(seed=1234):
@@ -32,7 +31,7 @@ def test_overwrite(seed=1234):
     np.random.seed(seed)
 
     def ll(x):
-        return -0.5 * np.sum(x**2)
+        return -0.5 * np.sum(x ** 2)
 
     nwalkers = 64
     p0 = np.random.normal(size=(nwalkers, 1))
diff --git a/emcee/tests/unit/test_stretch.py b/src/emcee/tests/unit/test_stretch.py
similarity index 71%
rename from emcee/tests/unit/test_stretch.py
rename to src/emcee/tests/unit/test_stretch.py
index 5dd64c85..a75e0521 100644
--- a/emcee/tests/unit/test_stretch.py
+++ b/src/emcee/tests/unit/test_stretch.py
@@ -1,15 +1,13 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
 import warnings
 
-import pytest
 import numpy as np
+import pytest
 
 from emcee import moves
-from emcee.state import State
 from emcee.model import Model
+from emcee.state import State
 
 __all__ = ["test_live_dangerously"]
 
@@ -19,14 +17,11 @@ def test_live_dangerously(nwalkers=32, nsteps=3000, seed=1234):
 
     # Set up the random number generator.
     np.random.seed(seed)
-    state = State(np.random.randn(nwalkers, 2 * nwalkers),
-                  log_prob=np.random.randn(nwalkers))
-    model = Model(
-        None,
-        lambda x: (np.zeros(len(x)), None),
-        map,
-        np.random
+    state = State(
+        np.random.randn(nwalkers, 2 * nwalkers),
+        log_prob=np.random.randn(nwalkers),
     )
+    model = Model(None, lambda x: (np.zeros(len(x)), None), map, np.random)
     proposal = moves.StretchMove()
 
     # Test to make sure that the error is thrown if there aren't enough
diff --git a/emcee/utils.py b/src/emcee/utils.py
similarity index 78%
rename from emcee/utils.py
rename to src/emcee/utils.py
index e5bd559e..42bbcfd4 100644
--- a/emcee/utils.py
+++ b/src/emcee/utils.py
@@ -1,14 +1,12 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import division, print_function
-
-__all__ = ["sample_ball", "deprecated", "deprecation_warning"]
-
 import warnings
 from functools import wraps
 
 import numpy as np
 
+__all__ = ["sample_ball", "deprecated", "deprecation_warning"]
+
 
 def deprecation_warning(msg):
     warnings.warn(msg, category=DeprecationWarning, stacklevel=2)
@@ -24,6 +22,7 @@ def wrapper(func, alternate=alternate):
         def f(*args, **kwargs):
             deprecation_warning(msg)
             return func(*args, **kwargs)
+
         return f
 
     return wrapper
@@ -39,9 +38,10 @@ def sample_ball(p0, std, size=1):
     :param size: The number of samples to produce.
 
     """
-    assert(len(p0) == len(std))
-    return np.vstack([p0 + std * np.random.normal(size=len(p0))
-                      for i in range(size)])
+    assert len(p0) == len(std)
+    return np.vstack(
+        [p0 + std * np.random.normal(size=len(p0)) for i in range(size)]
+    )
 
 
 @deprecated(None)
@@ -57,6 +57,6 @@ def sample_ellipsoid(p0, covmat, size=1):
     :param size: The number of samples to produce.
 
     """
-    return np.random.multivariate_normal(np.atleast_1d(p0),
-                                         np.atleast_2d(covmat),
-                                         size=size)
+    return np.random.multivariate_normal(
+        np.atleast_1d(p0), np.atleast_2d(covmat), size=size
+    )
diff --git a/tox.ini b/tox.ini
index d1bc4886..6c09ae37 100644
--- a/tox.ini
+++ b/tox.ini
@@ -8,4 +8,4 @@ deps =
     numpy
     scipy
     h5py
-commands = pytest -v emcee/tests
+commands = pytest -v src/emcee/tests