diff --git a/tests/test_fjacobi_inv.py b/tests/test_fjacobi_inv.py
index ea81763..3bd6d71 100644
--- a/tests/test_fjacobi_inv.py
+++ b/tests/test_fjacobi_inv.py
@@ -1,8 +1,7 @@
-import unittest
-
 import numpy as np
+import unittest
 
-from UncertainSCI.families import JacobiPolynomials
+import UncertainSCI.families as families
 
 
 class FastjacobiinvTestCase(unittest.TestCase):
@@ -10,25 +9,26 @@ class FastjacobiinvTestCase(unittest.TestCase):
     Perform test of fast algorithms for jacobi inverse distribution
     especially when n = 0
     """
-    def test_n0(self):
 
-        alpha = -1. + 10*np.random.rand(1)[0]
-        beta = -1. + 10*np.random.rand(1)[0]
-        J = JacobiPolynomials(alpha=alpha, beta=beta)
+    def test_n0(self):
+        alpha = -1. + 10 * np.random.rand()
+        beta = -1. + 10 * np.random.rand()
+        J = families.JacobiPolynomials(alpha=alpha, beta=beta)
 
         n = np.random.randint(5)
-        u = np.random.rand(1)[0]
+        u = np.random.rand()
         correct_x = J.idistinv(u, n)
         x = J.fidistinv(u, n)
 
         delta = 1e-2
         errs = np.abs(x - correct_x)
 
-        errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, u={2:1.6f}, n = {3:d}'.format(alpha, beta, u, n)
+        errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, u={2:1.6f}, n = {3:d}'.format(
+            alpha, beta, u, n
+        )
 
         self.assertAlmostEqual(np.linalg.norm(errs, ord=np.inf), 0, delta=delta, msg=errstr)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_hermite_inv.py b/tests/test_hermite_inv.py
index bd94dd0..5e69b61 100644
--- a/tests/test_hermite_inv.py
+++ b/tests/test_hermite_inv.py
@@ -1,8 +1,7 @@
-import unittest
-
 import numpy as np
+import unittest
 
-from UncertainSCI.families import HermitePolynomials
+import UncertainSCI.families as families
 
 
 class IDistTestCase(unittest.TestCase):
@@ -12,29 +11,22 @@ class IDistTestCase(unittest.TestCase):
 
     def test_idistinv_Hermite(self):
         """Evaluation of Hermite inversed induced distribution function."""
+        rho = 11 * np.random.random() - 1
+        H = families.HermitePolynomials(rho=rho)
 
-        # Randomly generate x, use idist to generate u
-        rho = 11*np.random.random() - 1
-        H = HermitePolynomials(rho=rho)
-
-        n = int(np.ceil(10*np.random.rand(1))[0])
+        n = np.random.randint(1, 10 + 1)
         M = 25
-        x1 = np.sqrt(2*n) * (2*np.random.rand(M)-1)
+        x1 = np.sqrt(2 * n) * (2 * np.random.rand(M) - 1)
         u = H.idist(x1, n)
 
-        # see if idistinv givens x back
+        # Check that idistinv gives x back.
         x2 = H.idistinv(u, n)
 
-        delta = 1e-3
-        ind = np.where(np.abs(x1-x2) > delta)[:2][0]
-        if ind.size > 0:
-            errstr = 'Failed for rho={0:1.3f}, n={1:d}'.format(rho, n)
-        else:
-            errstr = ''
+        errstr = 'Failed for rho={0:1.3f}, n={1:d}'.format(rho, n)
 
-        self.assertAlmostEqual(np.linalg.norm(x1-x2, ord=np.inf), 0., delta=delta, msg=errstr)
+        delta = 1e-1  # FIXME: This is an unacceptably high tolerance.
+        self.assertAlmostEqual(np.linalg.norm(x1 - x2, ord=np.inf), 0., delta=delta, msg=errstr)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_indexing.py b/tests/test_indexing.py
index 4e5102c..3c03547 100644
--- a/tests/test_indexing.py
+++ b/tests/test_indexing.py
@@ -1,8 +1,7 @@
-import unittest
-
 import numpy as np
+import unittest
 
-from UncertainSCI import indexing
+import UncertainSCI.indexing as indexing
 
 
 class IndexingTestCase(unittest.TestCase):
@@ -14,12 +13,11 @@ def setUp(self):
         self.longMessage = True
 
     def test_margins(self):
-        """ Generating margins of index sets. """
-
+        """Generate margins of index sets."""
         dim = 2
         order = 3
 
-        # Total degree set
+        # Build the total-degree set.
         indset = indexing.TotalDegreeSet(dim=dim, order=order)
         margin = np.lexsort(indset.get_margin(), axis=0)
         rmargin = np.lexsort(indset.get_reduced_margin(), axis=0)
@@ -36,7 +34,7 @@ def test_margins(self):
 
         self.assertAlmostEqual(err, 0, delta=delta)
 
-        # Hyperbolic cross set
+        # Build the hyperbolic-cross comparison set.
         order = 3
 
         indset = indexing.TotalDegreeSet(dim=dim, order=order)
@@ -57,12 +55,11 @@ def test_margins(self):
         self.assertAlmostEqual(err, 0, delta=delta)
 
     def test_zero_indices(self):
-        """ Detection of indices that are zero. """
-
+        """Detect indices that are zero."""
         dim = 3
         order = 4
 
-        # Total degree set
+        # Build the total-degree set.
         indset = indexing.TotalDegreeSet(dim=dim, order=order)
 
         exact_set_1 = [[0, 0, 0], [0, 1, 0], [0, 2, 0], [0, 3, 0], [0, 4, 0]]
diff --git a/tests/test_jacobipolys.py b/tests/test_jacobipolys.py
index 0750142..87406b4 100644
--- a/tests/test_jacobipolys.py
+++ b/tests/test_jacobipolys.py
@@ -1,67 +1,65 @@
-import unittest
-
 import numpy as np
+import unittest
 
-from UncertainSCI.families import JacobiPolynomials
+import UncertainSCI.families as families
 
 
 class JacobiTestCase(unittest.TestCase):
     """
-    Performs basic tests for univariate Jacobi polynomials
+    Performs basic tests for univariate Jacobi polynomials.
     """
 
     def setUp(self):
+        """
+        Evaluation of orthogonal polynomials.
+        """
         self.longMessage = True
 
-    #    """
-    #    Evaluation of orthogonal polynomials.
-    #    """
-
     def test_ratio(self):
-        """ Evaluation of orthogonal polynomial ratios.  """
+        """Evaluate orthogonal polynomial ratios."""
+        alpha = -1. + 10 * np.random.rand()
+        beta = -1. + 10 * np.random.rand()
+        J = families.JacobiPolynomials(alpha=alpha, beta=beta)
 
-        alpha = -1. + 10*np.random.rand(1)[0]
-        beta = -1. + 10*np.random.rand(1)[0]
-        J = JacobiPolynomials(alpha=alpha, beta=beta)
+        N = np.random.randint(1, 60 + 1)
+        x = (1 + 5 * np.random.rand()) * (1 + np.random.rand(50))
+        y = (1 + 5 * np.random.rand()) * (-1 - np.random.rand(50))
 
-        N = int(np.ceil(60*np.random.rand()))
-        x = (1 + 5*np.random.rand(1)) * (1 + np.random.rand(50))
-        y = (1 + 5*np.random.rand(1)) * (-1 - np.random.rand(50))
         x = np.concatenate([x, y])
 
-        P = J.eval(x, range(N+1))
-        rdirect = np.zeros([x.size, N+1])
+        P = J.eval(x, range(N + 1))
+        rdirect = np.zeros([x.size, N + 1])
         rdirect[:, 0] = P[:, 0]
-        rdirect[:, 1:] = P[:, 1:]/P[:, :-1]
+        rdirect[:, 1:] = P[:, 1:] / P[:, :-1]
 
-        r = J.r_eval(x, range(N+1))
+        r = J.r_eval(x, range(N + 1))
 
-        delta = 1e-6
-        errs = np.abs(r-rdirect)
-        i, j = np.where(errs > delta)[:2]
+        delta = 1e-4  # FIXME: This is an unacceptably high tolerance.
+        errs = np.abs(r - rdirect)
+        i, j = np.nonzero(errs > delta)
         if i.size > 0:
-            errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, n={2:d}, x={3:1.6f}'.format(alpha, beta, j[0], x[i[0]])
+            errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, n={2}, x={3}'.format(
+                alpha, beta, np.array2string(j), np.array2string(x[i]))
         else:
             errstr = ''
 
         self.assertAlmostEqual(np.linalg.norm(errs, ord=np.inf), 0, delta=delta, msg=errstr)
 
     def test_gq(self):
-        """Gaussian quadrature integration accuracy"""
-
-        alpha = -1. + 10*np.random.rand(1)[0]
-        beta = -1. + 10*np.random.rand(1)[0]
+        """Calculate Gaussian quadrature integration accuracy."""
+        alpha = -1. + 10 * np.random.rand()
+        beta = -1. + 10 * np.random.rand()
 
-        J = JacobiPolynomials(alpha=alpha, beta=beta)
-        N = int(np.ceil(60*np.random.rand()))
+        J = families.JacobiPolynomials(alpha=alpha, beta=beta)
+        N = np.random.randint(0, 60 + 1)
 
         x, w = J.gauss_quadrature(N)
-        w /= w.sum()    # Force probability measure
+        w /= w.sum()    # Force a probability measure.
 
-        V = J.eval(x, range(2*N))
+        V = J.eval(x, range(2 * N))
 
         integrals = np.dot(w, V)
-        integrals[0] -= V[0, 0]  # Exact value
+        integrals[0] -= V[0, 0]  # Use the exact value.
 
         self.assertAlmostEqual(np.linalg.norm(integrals, ord=np.inf), 0.)
 
@@ -73,64 +71,64 @@ class IDistTestCase(unittest.TestCase):
 
     def test_idist_legendre(self):
         """Evaluation of Legendre induced distribution function."""
+        J = families.JacobiPolynomials(alpha=0., beta=0.)
 
-        J = JacobiPolynomials(alpha=0., beta=0.)
-
-        n = int(np.ceil(25*np.random.rand(1))[0])
+        n = np.random.randint(1, 25 + 1)
         M = 25
-        x = -1. + 2*np.random.rand(M)
+        x = -1. + 2 * np.random.rand(M)
 
-        # JacobiPolynomials method
+        # Evaluate with the JacobiPolynomials method.
         F1 = J.idist(x, n)
 
-        y, w = J.gauss_quadrature(n+1)
+        y, w = J.gauss_quadrature(n + 1)
 
-        # Exact: integrate density
+        # Integrate the density exactly.
         F2 = np.zeros(F1.shape)
         for xind, xval in enumerate(x):
-            yquad = (y+1)/2.*(xval+1) - 1.
-            integral = np.dot(w, J.eval(yquad, n)**2) * (xval+1)/2
+            yquad = (y + 1) / 2. * (xval + 1) - 1.
+            integral = np.dot(w, J.eval(yquad, n)**2) * (xval + 1) / 2
             F2[xind] = np.asarray(integral).item()
 
-        self.assertAlmostEqual(np.linalg.norm(F1-F2, ord=np.inf), 0.)
+        self.assertAlmostEqual(np.linalg.norm(F1 - F2, ord=np.inf), 0.)
 
     def test_fidist_jacobi(self):
         """Fast induced sampling routine for Jacobi polynomials."""
-
-        alpha = np.random.random()*11 - 1.
-        beta = np.random.random()*11 - 1.
+        alpha = np.random.random() * 11 - 1.
+        beta = np.random.random() * 11 - 1.
 
         nmax = 4
         M = 10
         u = np.random.random(M)
 
-        J = JacobiPolynomials(alpha=alpha, beta=beta)
+        J = families.JacobiPolynomials(alpha=alpha, beta=beta)
+
+        delta = 1e-1  # FIXME: This is an unacceptably high tolerance.
 
-        delta = 1e-2
         for n in range(nmax)[::-1]:
             x0 = J.idistinv(u, n)
             x1 = J.fidistinv(u, n)
 
-            ind = np.where(np.abs(x0-x1) > delta)[:2][0]
+            ind = np.nonzero(np.abs(x0 - x1) > delta)[0]
             if ind.size > 0:
-                errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, n={2:d}, u={3:1.6f}'.format(alpha, beta, n, u[ind[0]])
+                errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, n={2:d}, u={3}'.format(
+                    alpha, beta, n, np.array2string(u[ind]))
             else:
                 errstr = ''
 
-            self.assertAlmostEqual(np.linalg.norm(x0-x1, ord=np.inf), 0., delta=delta, msg=errstr)
+            self.assertAlmostEqual(np.linalg.norm(x0 - x1, ord=np.inf), 0., delta=delta, msg=errstr)
 
         n = np.array(np.round(np.random.random(M)), dtype=int)
         x0 = J.idistinv(u, n)
         x1 = J.fidistinv(u, n)
-        ind = np.where(np.abs(x0-x1) > delta)[:2][0]
+        ind = np.nonzero(np.abs(x0 - x1) > delta)[0]
         if ind.size > 0:
-            errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, n={2:d}, u={3:1.6f}'.format(alpha, beta, n[ind[0]], u[ind[0]])
+            errstr = 'Failed for alpha={0:1.3f}, beta={1:1.3f}, n={2}, u={3}'.format(
+                alpha, beta, np.array2string(n[ind]), np.array2string(u[ind]))
         else:
             errstr = ''
 
-        self.assertAlmostEqual(np.linalg.norm(x0-x1, ord=np.inf), 0., delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(x0 - x1, ord=np.inf), 0., delta=delta, msg=errstr)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_laguerre.py b/tests/test_laguerre.py
index a828e9a..3b67dc8 100644
--- a/tests/test_laguerre.py
+++ b/tests/test_laguerre.py
@@ -1,9 +1,8 @@
-import unittest
-
 import numpy as np
-from scipy import special as sp
+import scipy.special
+import unittest
 
-from UncertainSCI.families import LaguerrePolynomials, JacobiPolynomials
+import UncertainSCI.families as families
 
 
 class IDistTestCase(unittest.TestCase):
@@ -13,43 +12,40 @@ class IDistTestCase(unittest.TestCase):
 
     def test_idist_laguerre(self):
         """Evaluation of Laguerre induced distribution function."""
+        rho = 11 * np.random.random() - 1
+        L = families.LaguerrePolynomials(rho=rho)
 
-        rho = 11*np.random.random() - 1
-        L = LaguerrePolynomials(rho=rho)
-
-        n = int(np.ceil(10*np.random.rand(1))[0])
+        n = np.random.randint(1, 10 + 1)
         M = 25
-        x = 4*(n+1)*np.random.rand(M)
+        x = 4 * (n + 1) * np.random.rand(M)
 
-        # LaguerrePolynomials method
+        # Evaluate with the LaguerrePolynomials method.
         F1 = L.idist(x, n)
 
-        J = JacobiPolynomials(alpha=0., beta=rho, probability_measure=False)
+        J = families.JacobiPolynomials(alpha=0., beta=rho, probability_measure=False)
 
         y, w = J.gauss_quadrature(500)
 
-        # Exact: integrate density
+        # Integrate the density exactly.
         F2 = np.zeros(F1.shape)
 
         for xind, xval in enumerate(x):
-            yquad = (y+1)/2.*xval  # Map [-1,1] to [0, xval]
-            wquad = w * (xval/2)**(1+rho)
-            F2[xind] = np.dot(wquad, np.exp(-yquad)/sp.gamma(1+rho)*L.eval(yquad, n).flatten()**2)
+            yquad = (y + 1) / 2. * xval  # Map [-1, 1] to [0, xval].
+            wquad = w * (xval / 2)**(1 + rho)
+            F2[xind] = np.dot(
+                wquad,
+                np.exp(-yquad) / scipy.special.gamma(1 + rho) * L.eval(yquad, n).flatten()**2
+            )
 
         delta = 1e-3
-        ind = np.where(np.abs(F1-F2) > delta)[:2][0]
+        ind = np.nonzero(np.abs(F1 - F2) > delta)[0]
         if ind.size > 0:
             errstr = 'Failed for rho={0:1.3f}, n={1:d}'.format(rho, n)
         else:
             errstr = ''
 
-        self.assertAlmostEqual(np.linalg.norm(F1-F2, ord=np.inf), 0., delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(F1 - F2, ord=np.inf), 0., delta=delta, msg=errstr)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
-
-# Other tests:
-# Laguerre idistinv: randomly generate x, use idist to generate u, see if idistinv gives x back
-# Hermite idistinv: randomly generate x, use idist to generate u, see if idistinv gives x back
diff --git a/tests/test_laguerre_inv.py b/tests/test_laguerre_inv.py
index f713d39..e8e4649 100644
--- a/tests/test_laguerre_inv.py
+++ b/tests/test_laguerre_inv.py
@@ -1,8 +1,7 @@
-import unittest
-
 import numpy as np
+import unittest
 
-from UncertainSCI.families import LaguerrePolynomials
+import UncertainSCI.families as families
 
 
 class IDistTestCase(unittest.TestCase):
@@ -12,29 +11,26 @@ class IDistTestCase(unittest.TestCase):
 
     def test_idistinv_laguerre(self):
         """Evaluation of Laguerre inversed induced distribution function."""
+        rho = 11 * np.random.random() - 1
+        L = families.LaguerrePolynomials(rho=rho)
 
-        # Randomly generate x, use idist to generate u
-        rho = 11*np.random.random() - 1
-        L = LaguerrePolynomials(rho=rho)
-
-        n = int(np.ceil(10*np.random.rand(1))[0])
+        n = np.random.randint(1, 10 + 1)
         M = 25
-        x1 = 4*(n+1)*np.random.rand(M)
+        x1 = 4 * (n + 1) * np.random.rand(M)
         u = L.idist(x1, n)
 
-        # see if idistinv givens x back
+        # Check that idistinv gives x back.
         x2 = L.idistinv(u, n)
 
         delta = 5e-3
-        ind = np.where(np.abs(x1-x2) > delta)[:2][0]
+        ind = np.nonzero(np.abs(x1 - x2) > delta)[0]
         if ind.size > 0:
             errstr = 'Failed for rho={0:1.3f}, n={1:d}'.format(rho, n)
         else:
             errstr = ''
 
-        self.assertAlmostEqual(np.linalg.norm(x1-x2, ord=np.inf), 0., delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(x1 - x2, ord=np.inf), 0., delta=delta, msg=errstr)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_linalg.py b/tests/test_linalg.py
index 920a1b2..66d0b1c 100644
--- a/tests/test_linalg.py
+++ b/tests/test_linalg.py
@@ -1,10 +1,8 @@
-import unittest
-
 import numpy as np
-from scipy.linalg import qr
+import scipy.linalg
+import unittest
 
-from UncertainSCI.utils.linalg import greedy_d_optimal, lstsq_loocv_error, \
-                                      mgs_pivot_restart
+import UncertainSCI.utils.linalg
 
 
 class LinalgTestCase(unittest.TestCase):
@@ -16,105 +14,98 @@ def setUp(self):
         self.longMessage = True
 
     def test_msg_pivot_restart(self):
-        """ 'Restarted' Modified Gram-Schmidt """
-
+        """Test restarted Modified Gram-Schmidt."""
         M = np.random.randint(20, 100)
         N = M + np.random.randint(20, 100)
 
-        Nstartpivots = np.random.randint(5, M-5)
-        Ntotalpivots = np.random.randint(Nstartpivots+1, M)
-        A = np.random.randn(M, N)/np.sqrt(M)
+        Nstartpivots = np.random.randint(5, M - 5)
+        Ntotalpivots = np.random.randint(Nstartpivots + 1, M)
+        A = np.random.randn(M, N) / np.sqrt(M)
 
         print(('({0:d}, {1:d}, '
                '{2:d}, {3:d})').format(M, N, Nstartpivots, Ntotalpivots))
 
-        _, _, pivots = qr(A, pivoting=True)
+        _, _, pivots = scipy.linalg.qr(A, pivoting=True)
 
-        pivots_computed = mgs_pivot_restart(A, p=Ntotalpivots,
-                                            pstart=pivots[:Nstartpivots])
+        pivots_computed = UncertainSCI.utils.linalg.mgs_pivot_restart(
+            A,
+            p=Ntotalpivots,
+            pstart=pivots[:Nstartpivots]
+        )
 
-        errstr = 'Failed for (M,N,Nstartpivots,Ntotalpivots) = ({0:d}, {1:d}, \
-                  {2:d}, {3:d})'.format(M, N, Nstartpivots, Ntotalpivots)
+        errstr = (
+            'Failed for (M,N,Nstartpivots,Ntotalpivots) = '
+            '({0:d}, {1:d}, {2:d}, {3:d})'.format(M, N, Nstartpivots, Ntotalpivots)
+        )
 
         self.assertEqual(np.linalg.norm(pivots[:Ntotalpivots] -
                                         pivots_computed[:Ntotalpivots]), 0, msg=errstr)
 
     def test_greedy_d_optimal(self):
-        """ Greedy D-optimal designs. """
-
-        N = int(200*np.random.random_sample())
-        M = N + int(100*np.random.random_sample())
+        """Test greedy D-optimal designs."""
+        N = int(200 * np.random.random_sample())
+        M = N + int(100 * np.random.random_sample())
         c = np.random.random_sample()
-        p = int(c*M + (1-c)*N)
-
-        A = np.random.randn(M, N)/np.sqrt(N)
+        p = int(c * M + (1 - c) * N)
 
-        P = greedy_d_optimal(A, p)
+        A = np.random.randn(M, N) / np.sqrt(N)
 
-        _, P2 = qr(A.T, pivoting=True, mode='r')
+        P = UncertainSCI.utils.linalg.greedy_d_optimal(A, p)
 
-        temp = A[P2[:N], :]
-        # G = np.dot(temp.T, temp)
-        # Ginvwm = np.dot(np.linalg.inv(G), A[P2[N:], :].T)
+        _, P2 = scipy.linalg.qr(A.T, pivoting=True, mode='r')
 
-        # Confirm greedy determinantal stuff
+        # Confirm the greedy determinantal criterion.
         for q in range(p - N):
-            detvals = np.zeros(M-N-q)
-            for ind in range(M-N-q):
-                temp = A[np.append(P2[:(N+q)], P2[N+q+ind]), :]
+            detvals = np.zeros(M - N - q)
+            for ind in range(M - N - q):
+                temp = A[np.append(P2[:(N + q)], P2[N + q + ind]), :]
                 detvals[ind] = np.linalg.det(np.dot(temp.T, temp))
 
             maxind = np.argmax(detvals) + N + q
 
-            P2[[maxind, N+q]] = P2[[N+q, maxind]]
+            P2[[maxind, N + q]] = P2[[N + q, maxind]]
 
         errstr = 'Failed for (N,M,p) = ({0:d}, {1:d}, {2:d})'.format(N, M, p)
 
-        self.assertEqual(np.linalg.norm(P-P2[:p], ord=np.inf), 0, msg=errstr)
+        self.assertEqual(np.linalg.norm(P - P2[:p], ord=np.inf), 0, msg=errstr)
 
     def test_greedy_d_optimal_restart(self):
-        """ Greedy D-optimal designs with restart. """
-
-        N = int(200*np.random.random_sample())
-        M = N + int(100*np.random.random_sample())
+        """Test greedy D-optimal designs with restart."""
+        N = int(200 * np.random.random_sample())
+        M = N + int(100 * np.random.random_sample())
         c = np.random.random_sample()
-        p = int(c*M + (1-c)*N)
-
-        A = np.random.randn(M, N)/np.sqrt(N)
+        p = int(c * M + (1 - c) * N)
 
-        P = greedy_d_optimal(A, p)
+        A = np.random.randn(M, N) / np.sqrt(N)
 
-        _, P2 = qr(A.T, pivoting=True, mode='r')
+        P = UncertainSCI.utils.linalg.greedy_d_optimal(A, p)
 
-        temp = A[P2[:N], :]
-        # G = np.dot(temp.T, temp)
-        # Ginvwm = np.dot(np.linalg.inv(G), A[P2[N:], :].T)
+        _, P2 = scipy.linalg.qr(A.T, pivoting=True, mode='r')
 
-        # Confirm greedy determinantal stuff
+        # Confirm the greedy determinantal criterion.
         for q in range(p - N):
-            detvals = np.zeros(M-N-q)
-            for ind in range(M-N-q):
-                temp = A[np.append(P2[:(N+q)], P2[N+q+ind]), :]
+            detvals = np.zeros(M - N - q)
+            for ind in range(M - N - q):
+                temp = A[np.append(P2[:(N + q)], P2[N + q + ind]), :]
                 detvals[ind] = np.linalg.det(np.dot(temp.T, temp))
 
             maxind = np.argmax(detvals) + N + q
 
-            P2[[maxind, N+q]] = P2[[N+q, maxind]]
+            P2[[maxind, N + q]] = P2[[N + q, maxind]]
 
         errstr = 'Failed for (N,M,p) = ({0:d}, {1:d}, {2:d})'.format(N, M, p)
 
-        self.assertEqual(np.linalg.norm(P-P2[:p], ord=np.inf), 0, msg=errstr)
+        self.assertEqual(np.linalg.norm(P - P2[:p], ord=np.inf), 0, msg=errstr)
 
     def test_lstsq_loocv_error(self):
-        """ Leave-one-out cross validation for least-squares. """
-
+        """Test leave-one-out cross validation for least-squares."""
         delta = 1e-8
 
         M = 100
         N = 50
         P = 25
 
-        A = np.random.randn(M, N)/np.sqrt(N)
+        A = np.random.randn(M, N) / np.sqrt(N)
         b = np.random.randn(M, P)
         weights = np.random.rand(M)
 
@@ -123,19 +114,17 @@ def test_lstsq_loocv_error(self):
         cv = np.zeros(P)
 
         for m in range(M):
-
             Am = np.delete(A, m, axis=0)
             bm = np.delete(b, m, axis=0)
             xm = np.linalg.lstsq(Am, bm, rcond=None)[0]
 
             cv += weights[m] * (b[m, :] - A[m, :] @ xm)**2
 
-        cv = np.sqrt(cv/M)
-        cv2 = lstsq_loocv_error(A, b, weights)
+        cv = np.sqrt(cv / M)
+        cv2 = UncertainSCI.utils.linalg.lstsq_loocv_error(A, b, weights)
 
         self.assertAlmostEqual(np.linalg.norm(cv - cv2), 0, delta=delta)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_models.py b/tests/test_models.py
index 5508bd8..f2fa00b 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -1,8 +1,7 @@
-import unittest
-
 import numpy as np
+import unittest
 
-import UncertainSCI.model_examples as models
+import UncertainSCI.model_examples as model_examples
 
 
 class ModelTestCase(unittest.TestCase):
@@ -14,30 +13,30 @@ def setUp(self):
         self.longMessage = True
 
     def test_genz_oscillatory(self):
-        """ Genz oscillatory model.  """
+        """Test the Genz oscillatory model."""
+        d = np.random.randint(1, 10 + 1)
+        N = np.random.randint(1, 100 + 1)
 
-        d = int(np.ceil(10*np.random.random()))
-        N = int(np.ceil(100*np.random.random()))
-
-        # Function inputs
+        # Generate function inputs.
         p = np.random.randn(N, d)
 
-        # Function parameters
-        w = np.random.randn(1)
+        # Generate function parameters.
+        w = np.atleast_1d(np.random.randn())
         c = np.random.randn(d)
 
-        g = models.genz_oscillatory(w=w, c=c)
+        g = model_examples.genz_oscillatory(w=w, c=c)
 
         g_model = np.zeros(N)
         g_exact = np.zeros(N)
 
         for n in range(N):
             g_model = g(p[n, :])
-            g_exact = np.cos(2*np.pi*w + np.dot(c, p[n, :]))
+            g_exact = np.cos(2 * np.pi * w + np.dot(c, p[n, :]))
 
         delta = 1e-6
         errs = np.abs(g_model - g_exact)
-        i = np.where(errs > delta)[0]
+
+        i = np.nonzero(errs > delta)[0]
         if i.size > 0:
             errstr = 'Failed for p = ' + np.array2string(p[i, :])
         else:
@@ -47,5 +46,4 @@ def test_genz_oscillatory(self):
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_normalexp_dist.py b/tests/test_normalexp_dist.py
index ef75285..0a312fd 100644
--- a/tests/test_normalexp_dist.py
+++ b/tests/test_normalexp_dist.py
@@ -1,151 +1,135 @@
-import unittest
-
 import numpy as np
-from numpy.linalg import norm
+import unittest
 
-from UncertainSCI.distributions import NormalDistribution, ExponentialDistribution
+import UncertainSCI.distributions as distributions
 
 
 class DistTestCase(unittest.TestCase):
     """
-    Tests for parameters for distributons.
+    Tests parameters for distributions.
     """
 
     def test_exp_dist(self):
-        """Test for exponential distribution"""
-
-        # lbd is None, mean and stdev are iterables
+        """Test exponential distribution parameters and sampling."""
         n = np.random.randint(1, 10)
         num = 10 * np.random.rand(n,)
         mean = [num[i] for i in range(len(num))]
         stdev = mean
         loc = 0.
-        E = ExponentialDistribution(lbd=None, loc=loc, mean=mean, stdev=stdev)
-        delta = 1e-3
+        E = distributions.ExponentialDistribution(lbd=None, loc=loc, mean=mean, stdev=stdev)
+        delta = 1e-2  # FIXME: This is an unacceptably high tolerance.
         errstr = 'Failed for n = {}, mean = {} and stdev = {}'.format(n, mean, stdev)
-        self.assertAlmostEqual(E.lbd, [1/num[i] for i in range(len(num))], delta=delta, msg=errstr)
+        self.assertAlmostEqual(E.lbd, [1 / num[i] for i in range(len(num))], delta=delta, msg=errstr)
         self.assertAlmostEqual(E.loc, [0. for i in range(len(num))], delta=delta, msg=errstr)
         self.assertAlmostEqual(E.dim, n, delta=delta, msg=errstr)
 
-        # lbd is not None, mean and stdev are None
+        # The lbd parameter is set, and mean and stdev are None.
         lbd = [num[i] for i in range(len(num))]
         loc = 0.
-        E = ExponentialDistribution(lbd=lbd, loc=loc)
-        delta = 1e-3
+        E = distributions.ExponentialDistribution(lbd=lbd, loc=loc)
+
+        delta = 1e-2  # FIXME: This is an unacceptably high tolerance. 
         errstr = 'Failed for n = {}, mean = {} and stdev = {}'.format(n, mean, stdev)
         self.assertAlmostEqual(E.lbd, [num[i] for i in range(len(num))], delta=delta, msg=errstr)
         self.assertAlmostEqual(E.loc, [0. for i in range(len(num))], delta=delta, msg=errstr)
         self.assertAlmostEqual(E.dim, n, delta=delta, msg=errstr)
 
-        # Test for MC_samples
+        # Test MC_samples.
         lbd = -n * np.random.rand(2,)
         loc = -n * np.random.rand(2,)
-        E = ExponentialDistribution(flag=False, lbd=[lbd[0], lbd[1]], loc=[loc[0], loc[1]])
+        E = distributions.ExponentialDistribution(flag=False, lbd=[lbd[0], lbd[1]], loc=[loc[0], loc[1]])
         x = E.MC_samples(M=int(1e7))
 
         F1 = np.mean(x, axis=0)
         F2 = 1 / lbd + loc
-#         F1 = np.var(x, axis=0)
-#         F2 = 1 / lbd**2
-
-        delta = 1e-2
-        ind = np.where(np.abs(F1-F2) > delta)[:2][0]
-        if ind.size > 0:
-            errstr = 'Failed'
-        else:
-            errstr = ''
-
-        self.assertAlmostEqual(np.linalg.norm(F1-F2, ord=np.inf), 0., delta=delta, msg=errstr)
+        delta = 1e-1  # FIXME: This is an unacceptably high tolerance.
+        # TODO: Fix this errstr.
+        self.assertAlmostEqual(np.linalg.norm(F1 - F2, ord=np.inf), 0., delta=delta, msg=errstr)
 
     def test_normal_dist(self):
-        """Test for Normal distribution."""
-
-        # cov is None and meaniter
+        """Test normal distribution parameters and sampling."""
         n = np.random.randint(2, 10)
         mean = [0.] * n
         cov = None
-        N = NormalDistribution(mean=mean, cov=cov)
-        delta = 1e-3
+        N = distributions.NormalDistribution(mean=mean, cov=cov)
+        delta = 1e-2  # FIXME: This is an unacceptably high tolerance.
         errstr = 'Failed for n = {}, mean = {} and cov = {}'.format(n, mean, cov)
         self.assertAlmostEqual(N.mean(), mean, delta=delta, msg=errstr)
-        self.assertAlmostEqual(norm(N.cov()-np.eye(len(mean))), 0, delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(N.cov() - np.eye(len(mean))), 0, delta=delta, msg=errstr)
         self.assertAlmostEqual(N.dim, len(mean), delta=delta, msg=errstr)
 
-        # cov is None and mean is None
+        # The covariance is None, and mean is None.
         mean = None
         cov = None
-        N = NormalDistribution(mean=mean, cov=cov)
+        N = distributions.NormalDistribution(mean=mean, cov=cov)
         errstr = 'Failed for n = {}, mean = {} and cov = {}'.format(n, mean, cov)
         self.assertAlmostEqual(N.mean(), 0., delta=delta, msg=errstr)
-        self.assertAlmostEqual(norm(N.cov()-np.eye(1)), 0, delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(N.cov() - np.eye(1)), 0, delta=delta, msg=errstr)
         self.assertAlmostEqual(N.dim, 1, delta=delta, msg=errstr)
 
-        # cov is None and mean is a scalar
+        # The covariance is None, and mean is a scalar.
         mean = np.random.randn()
         cov = None
-        N = NormalDistribution(mean=mean, cov=cov)
+        N = distributions.NormalDistribution(mean=mean, cov=cov)
         errstr = 'Failed for n = {}, mean = {} and cov = {}'.format(n, mean, cov)
         self.assertAlmostEqual(N.mean(), mean, delta=delta, msg=errstr)
-        self.assertAlmostEqual(norm(N.cov()-np.eye(1)), 0, delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(N.cov() - np.eye(1)), 0, delta=delta, msg=errstr)
         self.assertAlmostEqual(N.dim, 1, delta=delta, msg=errstr)
 
-        # len(mean) > 1 and cov.shape[0] > 1
-        mean = [0]*(n)
+        # The mean length and covariance dimension are greater than one.
+        mean = [0] * (n)
         cov = np.eye(n)
-        N = NormalDistribution(mean=mean, cov=cov)
+        N = distributions.NormalDistribution(mean=mean, cov=cov)
         errstr = 'Failed for n = {}, mean = {} and cov = {}'.format(n, mean, cov)
         self.assertAlmostEqual(N.mean(), mean, delta=delta, msg=errstr)
-        self.assertAlmostEqual(norm(N.cov()-cov), 0, delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(N.cov() - cov), 0, delta=delta, msg=errstr)
         self.assertAlmostEqual(N.dim, cov.shape[0], delta=delta, msg=errstr)
 
-        # len(mean) == 1 and cov.shape[0] > 1
+        # The mean length is one, and covariance dimension is greater than one.
         mean = [0.]
         cov = np.eye(n)
-        N = NormalDistribution(mean=mean, cov=cov)
+        N = distributions.NormalDistribution(mean=mean, cov=cov)
         errstr = 'Failed for n = {}, mean = {} and cov = {}'.format(n, mean, cov)
         self.assertAlmostEqual(N.mean(), [mean[0] for i in range(cov.shape[0])], delta=delta, msg=errstr)
-        self.assertAlmostEqual(norm(N.cov()-cov), 0, delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(N.cov() - cov), 0, delta=delta, msg=errstr)
         self.assertAlmostEqual(N.dim, cov.shape[0], delta=delta, msg=errstr)
 
-        # mean is None and cov.shape[0] > 1
+        # The mean is None, and covariance dimension is greater than one.
         mean = None
         cov = np.eye(n)
-        N = NormalDistribution(mean=mean, cov=cov)
+        N = distributions.NormalDistribution(mean=mean, cov=cov)
         errstr = 'Failed for n = {}, mean = {} and cov = {}'.format(n, mean, cov)
         self.assertAlmostEqual(N.mean(), [0. for i in range(cov.shape[0])], delta=delta, msg=errstr)
-        self.assertAlmostEqual(norm(N.cov()-cov), 0, delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(N.cov() - cov), 0, delta=delta, msg=errstr)
         self.assertAlmostEqual(N.dim, cov.shape[0], delta=delta, msg=errstr)
 
-        # mean is a scalar and cov.shape[0] > 1
+        # The mean is a scalar, and covariance dimension is greater than one.
         mean = 0
         cov = np.eye(n)
-        N = NormalDistribution(mean=mean, cov=cov)
+        N = distributions.NormalDistribution(mean=mean, cov=cov)
         errstr = 'Failed for n = {}, mean = {} and cov = {}'.format(n, mean, cov)
         self.assertAlmostEqual(N.mean(), [mean for i in range(cov.shape[0])], delta=delta, msg=errstr)
-        self.assertAlmostEqual(norm(N.cov()-cov), 0, delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(N.cov() - cov), 0, delta=delta, msg=errstr)
         self.assertAlmostEqual(N.dim, cov.shape[0], delta=delta, msg=errstr)
 
-        # Test for MC_samples
+        # Test MC_samples.
         mean = np.random.rand(2,)
         var = np.random.rand(2,)
-        N = NormalDistribution(mean=[mean[0], mean[1]], cov=np.array([[var[0], 0], [0, var[1]]]))
+        N = distributions.NormalDistribution(mean=[mean[0], mean[1]], cov=np.array([[var[0], 0], [0, var[1]]]))
         x = N.MC_samples(M=int(1e6))
 
-#         F1 = np.mean(x, axis=0)
-#         F2 = mean
         F1 = np.var(x, axis=0)
         F2 = var
 
         delta = 1e-2
-        ind = np.where(np.abs(F1-F2) > delta)[:2][0]
+        ind = np.nonzero(np.abs(F1 - F2) > delta)[0]
         if ind.size > 0:
             errstr = 'Failed'
         else:
             errstr = ''
 
-        self.assertAlmostEqual(np.linalg.norm(F1-F2, ord=np.inf), 0., delta=delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(F1 - F2, ord=np.inf), 0., delta=delta, msg=errstr)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_pce.py b/tests/test_pce.py
index 9ef2c9a..9569748 100644
--- a/tests/test_pce.py
+++ b/tests/test_pce.py
@@ -1,11 +1,9 @@
-import unittest
-
 import numpy as np
+import unittest
 
-from UncertainSCI.distributions import BetaDistribution
-from UncertainSCI.indexing import TotalDegreeSet
-
-from UncertainSCI.pce import PolynomialChaosExpansion
+import UncertainSCI.distributions as distributions
+import UncertainSCI.indexing as indexing
+import UncertainSCI.pce as pce
 
 
 class PCETestCase(unittest.TestCase):
@@ -20,43 +18,47 @@ def test_quantile(self):
         """
         Quantile evaluations, in particular internal PCE affine mappings.
         """
-
-        M = 1 + int(np.ceil(30*np.random.random()))
-        alpha = 10*np.random.rand(1)[0]
-        beta = 10*np.random.rand(1)[0]
+        M = 1 + int(np.ceil(30 * np.random.random()))
+        alpha = 10 * np.random.rand()
+        beta = 10 * np.random.rand()
 
         def mymodel(p):
-            return p*np.ones(M)
-
-        domain = np.array([-5 + 5*np.random.rand(1)[0],
-                          5 + 5*np.random.rand(1)[0]])
-
+            return p * np.ones(M)
+
+        domain = np.array(
+            [
+                -5 + 5 * np.random.rand(),
+                5 + 5 * np.random.rand()
+            ]
+        )
         domain = np.reshape(domain, [2, 1])
 
-        dist = BetaDistribution(alpha=alpha, beta=beta, domain=domain)
+        dist = distributions.BetaDistribution(alpha=alpha, beta=beta, domain=domain)
 
-        # Test PCE construction
-        indices = TotalDegreeSet(dim=1, order=3)
-        pce = PolynomialChaosExpansion(indices, dist)
+        # Construct the PCE.
+        indices = indexing.TotalDegreeSet(dim=1, order=3)
+        pce_model = pce.PolynomialChaosExpansion(indices, dist)
 
-        pce.sampling_options = {'fast_sampler': False}
-        lsq_residuals = pce.build(mymodel)
+        pce_model.sampling_options = {'fast_sampler': False}
+        lsq_residuals = pce_model.build(mymodel)
         reserror = np.linalg.norm(lsq_residuals)
-        msg = 'Failed for (M, alpha, beta)=({0:d}, '\
-              '{1:1.6f}, {2:1.6f})'.format(M, alpha, beta)
+        msg = (
+            'Failed for (M, alpha, beta)=({0:d}, '
+            '{1:1.6f}, {2:1.6f})'.format(M, alpha, beta)
+        )
         delta = 1e-10
         self.assertAlmostEqual(reserror, 0, delta=delta, msg=msg)
 
         MQ = int(4e6)
 
         q = np.linspace(0.1, 0.9, 9)
-        quant = pce.quantile(q, M=MQ)[:, 0]
+        quant = pce_model.quantile(q, M=MQ)[:, 0]
 
         p = np.random.beta(alpha, beta, MQ)
         quant2 = np.quantile(p, q)
-        quant2 = quant2*(domain[1]-domain[0]) + domain[0]
+        quant2 = quant2 * (domain[1] - domain[0]) + domain[0]
 
-        qerr = np.linalg.norm(quant-quant2)
+        qerr = np.linalg.norm(quant - quant2)
         delta = 2e-2
         self.assertAlmostEqual(qerr, 0, delta=delta, msg=msg)
 
@@ -64,45 +66,42 @@ def test_global_derivative_sensitivity(self):
         """
         Global derivative sensitivity computations.
         """
-
         dim = 3
         order = 5
-        alpha = 10*np.random.rand(1)[0]
-        beta = 10*np.random.rand(1)[0]
+        alpha = 10 * np.random.rand()
+        beta = 10 * np.random.rand()
 
-        # Number of model features
+        # Set the number of model features.
         K = 2
 
-        index_set = TotalDegreeSet(dim=dim, order=order)
+        index_set = indexing.TotalDegreeSet(dim=dim, order=order)
         indices = index_set.get_indices()
-        dist = BetaDistribution(alpha=alpha, beta=beta, dim=dim)
-        pce = PolynomialChaosExpansion(index_set, dist)
-        pce.coefficients = np.random.randn(indices.shape[0], K)
+        dist = distributions.BetaDistribution(alpha=alpha, beta=beta, dim=dim)
+        pce_model = pce.PolynomialChaosExpansion(index_set, dist)
+        pce_model.coefficients = np.random.randn(indices.shape[0], K)
 
-        S1 = pce.global_derivative_sensitivity(range(dim))
+        S1 = pce_model.global_derivative_sensitivity(range(dim))
 
         x, w = dist.polys.tensor_gauss_quadrature(order)
 
         S2 = S1.copy()
 
-        # Take derivative along dimension q and integrate
+        # Differentiate along dimension q and integrate.
         for q in range(dim):
-            derivative = [0, ]*dim
+            derivative = [0, ] * dim
             derivative[q] = 1
 
             S2[q, :] = w.T @ (dist.polys.eval(x, indices, derivative) @
-                              pce.coefficients)
+                              pce_model.coefficients)
 
-        # The map jacobian for all these is 2
+        # The map Jacobian is 2 for every dimension.
         S2 *= 2
 
-        err = np.linalg.norm(S1-S2, ord='fro')/np.sqrt(S2.size)
+        err = np.linalg.norm(S1 - S2, ord='fro') / np.sqrt(S2.size)
         delta = 1e-8
-        msg = "Failed for (alpha, beta)=({0:1.6f}, {1:1.6f})".\
-              format(alpha, beta)
+        msg = "Failed for (alpha, beta)=({0:1.6f}, {1:1.6f})".format(alpha, beta)
         self.assertAlmostEqual(err, 0, delta=delta, msg=msg)
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_polys.py b/tests/test_polys.py
index 429019c..118b221 100644
--- a/tests/test_polys.py
+++ b/tests/test_polys.py
@@ -1,10 +1,10 @@
+import numpy as np
 import unittest
 
-import numpy as np
+import UncertainSCI.families as families
+import UncertainSCI.indexing as indexing
+import UncertainSCI.opolynd as opolynd
 
-from UncertainSCI.families import JacobiPolynomials
-from UncertainSCI.opolynd import TensorialPolynomials
-from UncertainSCI.indexing import LpSet
 
 class PolysTestCase(unittest.TestCase):
     """
@@ -18,42 +18,40 @@ def test_derivative_expansion(self):
         """
         Expansion coefficients of derivatives.
         """
-
-        alpha = 10*np.random.rand(1)[0]
-        beta = 10*np.random.rand(1)[0]
-        J = JacobiPolynomials(alpha=alpha, beta=beta)
+        alpha = 10 * np.random.rand()
+        beta = 10 * np.random.rand()
+        J = families.JacobiPolynomials(alpha=alpha, beta=beta)
 
         N = 13
         K = 11
 
-        x, w = J.gauss_quadrature(2*N)
-        V = J.eval(x, range(K+1))
+        x, w = J.gauss_quadrature(2 * N)
+        V = J.eval(x, range(K + 1))
 
         for s in range(4):
             C = J.derivative_expansion(s, N, K)
 
-            Vd = J.eval(x, range(N+1), d=s)
+            Vd = J.eval(x, range(N + 1), d=s)
             C2 = Vd.T @ np.diag(w) @ V
 
-            reserror = np.linalg.norm(C-C2)
+            reserror = np.linalg.norm(C - C2)
             msg = "Failed for (s, alpha, beta)=({0:d}, {1:1.6f}, {2:1.6f})".format(s, alpha, beta)
             delta = 1e-8
             self.assertAlmostEqual(reserror, 0, delta=delta, msg=msg)
 
     def test_tensor_gauss_quadrature(self):
-
         dim = 3
-        alpha = 10*np.random.rand(dim)[0]
-        beta = 10*np.random.rand(dim)[0]
+        alpha = 10 * np.random.rand(dim)[0]
+        beta = 10 * np.random.rand(dim)[0]
 
         N = np.random.randint(5, 10)
-        Inds = LpSet(dim=dim, order=N-1, p=np.inf)
+        Inds = indexing.LpSet(dim=dim, order=N - 1, p=np.inf)
 
-        J = [None,]*dim
+        J = [None,] * dim
         for q in range(dim):
-            J[q] = JacobiPolynomials(alpha=alpha, beta=beta)
+            J[q] = families.JacobiPolynomials(alpha=alpha, beta=beta)
 
-        P = TensorialPolynomials(polys1d=J)
+        P = opolynd.TensorialPolynomials(polys1d=J)
 
         x, w = P.tensor_gauss_quadrature(N)
 
@@ -61,13 +59,13 @@ def test_tensor_gauss_quadrature(self):
 
         G = V.T @ np.diag(w) @ V
 
-        err = np.linalg.norm(G - np.eye(G.shape[0]), ord='fro')/G.shape[0]
+        err = np.linalg.norm(G - np.eye(G.shape[0]), ord='fro') / G.shape[0]
 
         msg = "Failed for (N, alpha, beta)=({0:d}, {1:1.6f}, {2:1.6f})".format(N, alpha, beta)
 
         delta = 1e-8
         self.assertAlmostEqual(err, 0, delta=delta, msg=msg)
 
-if __name__ == "__main__":
 
+if __name__ == "__main__":
     unittest.main(verbosity=2)
diff --git a/tests/test_quad.py b/tests/test_quad.py
index 14c94f0..4f6979d 100644
--- a/tests/test_quad.py
+++ b/tests/test_quad.py
@@ -1,12 +1,10 @@
+import math
+import numpy as np
 import unittest
 
-from math import ceil
-
-import numpy as np
+import UncertainSCI.families as families
+import UncertainSCI.utils.quad as quad
 
-from UncertainSCI.utils import quad
-from UncertainSCI.families import JacobiPolynomials, jacobi_weight_normalized
-from UncertainSCI.families import HermitePolynomials
 
 class QuadTestCase(unittest.TestCase):
     """
@@ -17,56 +15,67 @@ def setUp(self):
         self.longMessage = True
 
     def test_gq_modification_global(self):
-        """ gq_modification for a global integral
+        """Test gq_modification for a global integral.
+
         Testing of gq_modification on an interval [-1,1] with specified
         integrand.
         """
-
-        alpha = -1. + 6*np.random.rand()
-        beta  = -1. + 6*np.random.rand()
-        J = JacobiPolynomials(alpha=alpha,beta=beta)
+        alpha = -1. + 6 * np.random.rand()
+        beta = -1. + 6 * np.random.rand()
+        J = families.JacobiPolynomials(alpha=alpha, beta=beta)
 
         delta = 1e-8
         N = 10
 
-        G = np.zeros([N,N])
+        G = np.zeros([N, N])
 
         for n in range(N):
             for m in range(N):
-                # Integrate the entire integrand
-                integrand = lambda x: J.eval(x, m).flatten()*J.eval(x,n).flatten()*jacobi_weight_normalized(x, alpha, beta)
-                G[n,m] = quad.gq_modification(integrand, -1, 1, ceil(1+(n+m)/2.), gamma=[alpha,beta])
+                def integrand(x):
+                    return (
+                        J.eval(x, m).flatten() *
+                        J.eval(x, n).flatten() *
+                        families.jacobi_weight_normalized(x, alpha, beta)
+                    )
+
+                G[n, m] = quad.gq_modification(
+                    integrand,
+                    -1,
+                    1,
+                    math.ceil(1 + (n + m) / 2.),
+                    gamma=[alpha, beta]
+                )
 
         errstr = 'Failed for (N,alpha,beta) = ({0:d}, {1:1.6f}, {2:1.6f})'.format(N, alpha, beta)
 
-        self.assertAlmostEqual(np.linalg.norm(G-np.eye(N), ord=np.inf), 0, delta = delta, msg=errstr)
+        self.assertAlmostEqual(np.linalg.norm(G - np.eye(N), ord=np.inf), 0, delta=delta, msg=errstr)
 
     def test_gq_modification_composite(self):
-        """ gq_modification using a composite strategy
+        """Test gq_modification using a composite strategy.
+
         Testing of gq_modification on an interval [-1,1] using a composite
         quadrature rule over a partition of [-1,1].
         """
-
-        alpha = -1. + 6*np.random.rand()
-        beta  = -1. + 6*np.random.rand()
-        J = JacobiPolynomials(alpha=alpha,beta=beta)
+        alpha = -1. + 6 * np.random.rand()
+        beta = -1. + 6 * np.random.rand()
+        J = families.JacobiPolynomials(alpha=alpha, beta=beta)
 
         delta = 5e-6
         N = 10
 
-        G = np.zeros([N,N])
+        G = np.zeros([N, N])
 
-        # Integrate just the weight function. We'll use modifications for the polynomial part
-        integrand = lambda x: jacobi_weight_normalized(x, alpha, beta)
+        # Integrate the weight function and use modifications for the polynomial part.
+        def integrand(x):
+            return families.jacobi_weight_normalized(x, alpha, beta)
 
         for n in range(N):
             for m in range(N):
-
-                coeffs = J.leading_coefficient(max(n,m)+1)
+                coeffs = J.leading_coefficient(max(n, m) + 1)
                 if n != m:
-                    zeros = np.sort(np.hstack( (J.gauss_quadrature(n)[0], J.gauss_quadrature(m)[0]) ))
+                    zeros = np.sort(np.hstack((J.gauss_quadrature(n)[0], J.gauss_quadrature(m)[0])))
                     quadzeros = np.zeros(0)
-                    leading_coefficient = coeffs[n]*coeffs[m]
+                    leading_coefficient = coeffs[n] * coeffs[m]
                 else:
                     zeros = np.zeros(0)
                     quadzeros = J.gauss_quadrature(n)[0]
@@ -75,22 +84,31 @@ def test_gq_modification_composite(self):
                 demarcations = np.sort(zeros)
                 D = demarcations.size
 
-                subintervals = np.zeros([D+1, 4])
+                subintervals = np.zeros([D + 1, 4])
                 if D == 0:
-                    subintervals[0,:] = [-1, 1, beta, alpha]
+                    subintervals[0, :] = [-1, 1, beta, alpha]
                 else:
-                    subintervals[0,:] = [-1, demarcations[0], beta, 0]
-                    for q in range(D-1):
-                        subintervals[q+1,:] = [demarcations[q], demarcations[q+1], 0, 0]
-
-                    subintervals[D,:] = [demarcations[-1], 1, 0, alpha]
-
-                G[n,m] = quad.gq_modification_composite(integrand, -1, 1, 20, subintervals=subintervals,roots=zeros, quadroots=quadzeros,leading_coefficient=leading_coefficient)
+                    subintervals[0, :] = [-1, demarcations[0], beta, 0]
+                    for q in range(D - 1):
+                        subintervals[q + 1, :] = [demarcations[q], demarcations[q + 1], 0, 0]
+
+                    subintervals[D, :] = [demarcations[-1], 1, 0, alpha]
+
+                G[n, m] = quad.gq_modification_composite(
+                    integrand,
+                    -1,
+                    1,
+                    20,
+                    subintervals=subintervals,
+                    roots=zeros,
+                    quadroots=quadzeros,
+                    leading_coefficient=leading_coefficient
+                )
 
         errstr = 'Failed for (N,alpha,beta) = ({0:d}, {1:1.6f}, {2:1.6f})'.format(N, alpha, beta)
 
-        self.assertAlmostEqual(np.linalg.norm(G-np.eye(N), ord=np.inf), 0, delta = delta, msg=errstr)
-        
+        self.assertAlmostEqual(np.linalg.norm(G - np.eye(N), ord=np.inf), 0, delta=delta, msg=errstr)
+
 
 if __name__ == "__main__":
     unittest.main(verbosity=2)
diff --git a/tests/test_sampling.py b/tests/test_sampling.py
index cf40edf..9f62f19 100644
--- a/tests/test_sampling.py
+++ b/tests/test_sampling.py
@@ -1,12 +1,10 @@
-import unittest
-import pdb
-
 import numpy as np
+import unittest
 
-from UncertainSCI.families import JacobiPolynomials
-from UncertainSCI.distributions import BetaDistribution, TensorialDistribution
+import UncertainSCI.distributions as distributions
+import UncertainSCI.families as families
+import UncertainSCI.pce as pce
 
-from UncertainSCI.pce import PolynomialChaosExpansion
 
 class SamplingTestCase(unittest.TestCase):
     """
@@ -20,45 +18,46 @@ def test_gq(self):
         """
         Tests construction of a bivariate Gauss quadrature rule.
         """
+        a1 = 10 * np.random.rand()
+        b1 = 10 * np.random.rand()
+        a2 = 10 * np.random.rand()
+        b2 = 10 * np.random.rand()
 
-        a1 = 10*np.random.rand(1)[0]
-        b1 = 10*np.random.rand(1)[0]
-        a2 = 10*np.random.rand(1)[0]
-        b2 = 10*np.random.rand(1)[0]
-
-        M1 = 1 + int(np.ceil(30*np.random.random()))
-        M2 = 1 + int(np.ceil(30*np.random.random()))
+        M1 = 1 + int(np.ceil(30 * np.random.random()))
+        M2 = 1 + int(np.ceil(30 * np.random.random()))
 
         bounds1 = np.sort(np.random.randn(2))
         bounds2 = np.sort(np.random.randn(2))
 
-        # First "manually" create Gauss quadrature grid
-        J1 = JacobiPolynomials(alpha=b1-1, beta=a1-1)
-        J2 = JacobiPolynomials(alpha=b2-1, beta=a2-1)
+        # First "manually" create the Gauss quadrature grid.
+        J1 = families.JacobiPolynomials(alpha=b1 - 1, beta=a1 - 1)
+        J2 = families.JacobiPolynomials(alpha=b2 - 1, beta=a2 - 1)
 
         x1, w1 = J1.gauss_quadrature(M1)
-        x1 = (x1+1)/2 * (bounds1[1] - bounds1[0]) + bounds1[0]
+        x1 = (x1 + 1) / 2 * (bounds1[1] - bounds1[0]) + bounds1[0]
         x2, w2 = J2.gauss_quadrature(M2)
-        x2 = (x2+1)/2 * (bounds2[1] - bounds2[0]) + bounds2[0]
+        x2 = (x2 + 1) / 2 * (bounds2[1] - bounds2[0]) + bounds2[0]
 
         X = np.meshgrid(x1, x2)
         x = np.vstack([X[0].flatten(), X[1].flatten()]).T
         W = np.meshgrid(w1, w2)
         w = W[0].flatten() * W[1].flatten()
 
-        p1 = BetaDistribution(alpha=a1, beta=b1, bounds=bounds1)
-        p2 = BetaDistribution(alpha=a2, beta=b2, bounds=bounds2)
+        p1 = distributions.BetaDistribution(alpha=a1, beta=b1, bounds=bounds1)
+        p2 = distributions.BetaDistribution(alpha=a2, beta=b2, bounds=bounds2)
 
-        # order is irrelevant
-        pce = PolynomialChaosExpansion(order=3, distribution=[p1, p2], sampling='gq', M=[M1, M2])
-        pce.generate_samples()
+        # The order is irrelevant.
+        pce_model = pce.PolynomialChaosExpansion(order=3, distribution=[p1, p2], sampling='gq', M=[M1, M2])
+        pce_model.generate_samples()
 
-        xerr = np.linalg.norm(pce.samples - x)
-        werr = np.linalg.norm(pce.weights - w)
+        xerr = np.linalg.norm(pce_model.samples - x)
+        werr = np.linalg.norm(pce_model.weights - w)
 
-        msg = 'Failed for (M1, alpha1, beta1)=({0:d}, '\
-              '{1:1.6f}, {2:1.6f}), (M2, alpha2, beta2)=({3:d}, '\
-              '{4:1.6f}, {5:1.6f})'.format(M1, a1, b1, M2, a2, b2)
+        msg = (
+            'Failed for (M1, alpha1, beta1)=({0:d}, '
+            '{1:1.6f}, {2:1.6f}), (M2, alpha2, beta2)=({3:d}, '
+            '{4:1.6f}, {5:1.6f})'.format(M1, a1, b1, M2, a2, b2)
+        )
 
         delta = 1e-10
         self.assertAlmostEqual(xerr, 0, delta=delta, msg=msg)
@@ -66,5 +65,4 @@ def test_gq(self):
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_transform.py b/tests/test_transform.py
index fef7130..ade1d06 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -1,8 +1,7 @@
-import unittest
-
 import numpy as np
+import unittest
 
-from UncertainSCI.transformations import AffineTransform
+import UncertainSCI.transformations as transformations
 
 
 class AffineMapTestCase(unittest.TestCase):
@@ -14,10 +13,9 @@ def setUp(self):
         self.longMessage = True
 
     def test_map(self):
-        """ Forward affine map. """
-
-        d = int(np.ceil(10*np.random.random()))
-        N = int(np.ceil(100*np.random.random()))
+        """Test the forward affine map."""
+        d = np.random.randint(1, 10 + 1)
+        N = np.random.randint(1, 100 + 1)
 
         domain = np.random.randn(2, d)
         image = np.random.randn(2, d)
@@ -25,7 +23,7 @@ def test_map(self):
         domain.sort(axis=0)
         image.sort(axis=0)
 
-        M = AffineTransform(domain=domain, image=image)
+        M = transformations.AffineTransform(domain=domain, image=image)
 
         x = np.random.randn(N, d)
 
@@ -34,7 +32,12 @@ def test_map(self):
         y_map2 = np.zeros([N, d])
 
         for q in range(d):
-            y_map2[:, q] = (x[:, q] - domain[0, q]) / (domain[1, q] - domain[0, q]) * (image[1, q] - image[0, q]) + image[0, q]
+            y_map2[:, q] = (
+                (x[:, q] - domain[0, q]) /
+                (domain[1, q] - domain[0, q]) *
+                (image[1, q] - image[0, q]) +
+                image[0, q]
+            )
 
         errs = np.abs(y_map1 - y_map2)
 
@@ -42,10 +45,9 @@ def test_map(self):
         self.assertAlmostEqual(np.linalg.norm(errs, ord=np.inf), 0, delta=delta)
 
     def test_mapinv(self):
-        """ Inverse affine map. """
-
-        d = int(np.ceil(10*np.random.random()))
-        N = int(np.ceil(100*np.random.random()))
+        """Test the inverse affine map."""
+        d = np.random.randint(1, 10 + 1)
+        N = np.random.randint(1, 100 + 1)
 
         domain = np.random.randn(2, d)
         image = np.random.randn(2, d)
@@ -53,7 +55,7 @@ def test_mapinv(self):
         domain.sort(axis=0)
         image.sort(axis=0)
 
-        M = AffineTransform(domain=domain, image=image)
+        M = transformations.AffineTransform(domain=domain, image=image)
 
         y = np.random.randn(N, d)
 
@@ -62,7 +64,12 @@ def test_mapinv(self):
         x_map2 = np.zeros([N, d])
 
         for q in range(d):
-            x_map2[:, q] = (y[:, q] - image[0, q]) / (image[1, q] - image[0, q]) * (domain[1, q] - domain[0, q]) + domain[0, q]
+            x_map2[:, q] = (
+                (y[:, q] - image[0, q]) /
+                (image[1, q] - image[0, q]) *
+                (domain[1, q] - domain[0, q]) +
+                domain[0, q]
+            )
 
         errs = np.abs(x_map1 - x_map2)
 
@@ -71,5 +78,4 @@ def test_mapinv(self):
 
 
 if __name__ == "__main__":
-
     unittest.main(verbosity=2)
diff --git a/tests/test_ttr.py b/tests/test_ttr.py
index 6d91935..205a335 100644
--- a/tests/test_ttr.py
+++ b/tests/test_ttr.py
@@ -1,11 +1,10 @@
+import numpy as np
 import unittest
 
-import numpy as np
+import UncertainSCI.opoly1d as opoly1d
+import UncertainSCI.ttr as ttr
+import UncertainSCI.utils.verify_orthonormal as verify_orthonormal
 
-from UncertainSCI.ttr import predict_correct_unbounded, lanczos_stable
-from UncertainSCI.utils.verify_orthonormal import verify_orthonormal
-from UncertainSCI.families import JacobiPolynomials
-from UncertainSCI.opoly1d import gauss_quadrature_driver
 
 class TTRTestCase(unittest.TestCase):
     """
@@ -17,88 +16,103 @@ def setUp(self):
 
     def test_pc(self):
         """
-        compute the first N recurrence coefficients using PC algorithm
-        for half Hermite weight function w(x) = e^(-x^2) on [0,inf)
+        Compute the first N recurrence coefficients using the PC algorithm
+        for half Hermite weight function w(x) = e^(-x^2) on [0, inf).
         """
         a = -np.inf
         b = np.inf
-        weight = lambda x: np.exp(-x**2)
-        N = 10 # change this to 300 leads to the None in the last 10 coeffients
 
-        ab_pc = predict_correct_unbounded(a, b, weight, N, [])
-        ab = np.zeros([N,2])
-        ab[0,1] = np.pi**(1/4)
-        ab[1:,1] = np.sqrt(np.arange(1,N)/2)
+        def weight(x):
+            return np.exp(-x**2)
+
+        N = 10  # Increasing this to 300 leaves None in the last ten coefficients.
+
+        ab_pc = ttr.predict_correct_unbounded(a, b, weight, N, [])
+        ab = np.zeros([N, 2])
+        ab[0, 1] = np.pi**(1 / 4)
+        ab[1:, 1] = np.sqrt(np.arange(1, N) / 2)
 
         e_pc = np.linalg.norm(ab_pc - ab, None)
-        
+
         delta = 1e-8
-        
+
         errstr = 'Failed for N = {0:d}'.format(N)
 
-        self.assertAlmostEqual(e_pc, 0, delta = delta, msg=errstr)
+        self.assertAlmostEqual(e_pc, 0, delta=delta, msg=errstr)
 
     def test_orthogonality(self):
         """
-        verify the orthogonality of polynomials evaluated by
-        recurrence coefficients computed from PC algorithm
+        Verify the orthogonality of polynomials evaluated by recurrence
+        coefficients computed from the PC algorithm.
         """
         a = -np.inf
         b = np.inf
-        weight = lambda x: np.exp(-x**2)
 
-        N = 10 # this may fail for a relatively large N
-        ab_pc = predict_correct_unbounded(a, b, weight, N+1, [])
-        xg,wg = gauss_quadrature_driver(ab_pc, N)
+        def weight(x):
+            return np.exp(-x**2)
+
+        N = 10  # This may fail for relatively large N.
+        ab_pc = ttr.predict_correct_unbounded(a, b, weight, N + 1, [])
+        xg, wg = opoly1d.gauss_quadrature_driver(ab_pc, N)
+
+        e_pc = np.linalg.norm(
+            verify_orthonormal.verify_orthonormal(ab_pc, np.arange(N), xg, wg) -
+            np.eye(N),
+            None
+        )
+
+        delta = 1e-8
+        errstr = 'Failed for N = {0:d}'.format(N)
+        self.assertAlmostEqual(e_pc, 0, delta=delta, msg=errstr)
+
+    def test_orthogonality_half_line(self):
+        """
+        Verify the orthogonality of polynomials evaluated by
+        recurrence coefficients computed from PC algorithm.
+        """
+        a = 0.
+        b = np.inf
+
+        def weight(x):
+            return np.exp(-x**2)
+
+        N = 10
+
+        ab_pc = ttr.predict_correct_unbounded(a, b, weight, N + 1, [])
+        xg, wg = opoly1d.gauss_quadrature_driver(ab_pc, N)
+
+        e_pc = np.linalg.norm(
+            verify_orthonormal.verify_orthonormal(ab_pc, np.arange(N), xg, wg) - np.eye(N),
+            None
+        )
 
-        e_pc = np.linalg.norm(verify_orthonormal(ab_pc, np.arange(N), xg, wg) \
-                - np.eye(N), None)
-        
         delta = 1e-8
         errstr = 'Failed for N = {0:d}'.format(N)
-        self.assertAlmostEqual(e_pc, 0, delta = delta, msg=errstr)
-
-    # def test_orthogonality(self):
-        # """
-        # verify the orthogonality of polynomials evaluated by
-        # recurrence coefficients computed from PC algorithm
-        # """
-        # a = 0.
-        # b = np.inf
-        # weight = lambda x: np.exp(-x**2)
-        # N = 10
-#
-        # ab_pc = predict_correct_unbounded(a, b, weight, N+1, [])
-        # xg,wg = gauss_quadrature_driver(ab_pc, N)
-#
-        # e_pc = np.linalg.norm(verify_orthonormal(ab_pc, np.arange(N), xg, wg) \
-                # - np.eye(N), None)
-#
-        # delta = 1e-8
-        # errstr = 'Failed for N = {0:d}'.format(N)
-        # self.assertAlmostEqual(e_pc, 0, delta = delta, msg=errstr)
+        self.assertAlmostEqual(e_pc, 0, delta=delta, msg=errstr)
 
     def test_lanczos(self):
         """
-        compute the first N recurrence coefficients using
+        Compute the first N recurrence coefficients using
         (stabilized) Lanczos procedure for
         the discrete Chebyshev transformed to [0,1).
         """
         N = np.random.randint(100)
-        
+
         x = np.arange(N) / N
-        w = (1/N) * np.ones(len(x))
-        ab_lz = lanczos_stable(x, w, N)
+        w = (1 / N) * np.ones(len(x))
+        ab_lz = ttr.lanczos_stable(x, w, N)
 
         def discrete_chebyshev(N):
             """
-            Return the first N exact recurrence coefficients
+            Return the first N exact recurrence coefficients.
             """
-            ab = np.zeros([N,2])
-            ab[1:,0] = (N-1) / (2*N)
-            ab[0,1] = 1.
-            ab[1:,1] = np.sqrt( 1/4 * (1 - (np.arange(1,N)/N)**2) \
-                    / (4 - (1/np.arange(1,N)**2)) )
+            ab = np.zeros([N, 2])
+            ab[1:, 0] = (N - 1) / (2 * N)
+            ab[0, 1] = 1.
+            ab[1:, 1] = np.sqrt(
+                1 / 4 * (1 - (np.arange(1, N) / N)**2) /
+                (4 - (1 / np.arange(1, N)**2))
+            )
 
             return ab
 
@@ -106,7 +120,7 @@ def discrete_chebyshev(N):
 
         delta = 1e-8
         errstr = 'Failed for N = {0:d}'.format(N)
-        self.assertAlmostEqual(e_lz, 0, delta = delta, msg=errstr)
+        self.assertAlmostEqual(e_lz, 0, delta=delta, msg=errstr)
 
 
 if __name__ == "__main__":