diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
new file mode 100644
index 0000000..1957d94
--- /dev/null
+++ b/.github/workflows/ci.yaml
@@ -0,0 +1,28 @@
+name: CI
+on:
+  push:
+  pull_request:
+
+jobs:
+  test:
+    name: test ${{ matrix.python_version }}
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        python_version:
+          - "3.10"
+          - "3.9"
+          - "3.8"
+          - "3.7"
+    steps:
+      - uses: actions/setup-python@v2
+        with:
+          python-version: ${{ matrix.python_version }}
+      - uses: actions/checkout@v2
+      - name: Install tox-gh
+        run: python -m pip install tox-gh
+      - name: Setup test suite
+        run: tox -vv --notest
+      - name: Run test suite
+        run: tox --skip-pkg-install
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index 0a33029..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,16 +0,0 @@
-language: python
-
-matrix:
-  include:
-    - python: 2.7
-    - python: 3.5
-    - python: 3.6
-    - python: 3.7
-      dist: xenial
-      sudo: true
-
-# command to install dependencies
-install: pip install tox-travis
-
-# command to run tests
-script: tox
diff --git a/README.md b/README.md
index c811e81..ebb0aa1 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 Author: Roan LaPlante <rlaplant@nmr.mgh.harvard.edu>
 
-Tested against python 2.7 and 3.9.
+Tested against python 3.7+.
 
 ## Copyright information
 
diff --git a/bct/algorithms/clustering.py b/bct/algorithms/clustering.py
index 484cc23..5732b3d 100644
--- a/bct/algorithms/clustering.py
+++ b/bct/algorithms/clustering.py
@@ -250,7 +250,7 @@ def clustering_coef_wu_sign(W, coef_type='default'):
     '''
     Returns the weighted clustering coefficient generalized or separated
     for positive and negative weights.
-  
+
     Three Algorithms are supported; herefore referred to as default, zhang,
     and costantini.
 
@@ -384,9 +384,8 @@ def consensus_und(D, tau, reps=1000, seed=None):
     reps : int
         number of times the clustering algorithm is reapplied. default value
         is 1000.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -485,7 +484,7 @@ def get_components(A, no_depend=False):
     if not np.all(A == A.T):  # ensure matrix is undirected
         raise BCTParamError('get_components can only be computed for undirected'
                             ' matrices.  If your matrix is noisy, correct it with np.around')
-    
+
     A = binarize(A, copy=True)
     n = len(A)
     np.fill_diagonal(A, 1)
@@ -503,7 +502,7 @@ def get_components(A, no_depend=False):
         temp.append(item)
         union_sets = temp
 
-    comps = np.array([i+1 for v in range(n) for i in 
+    comps = np.array([i+1 for v in range(n) for i in
         range(len(union_sets)) if v in union_sets[i]])
     comp_sizes = np.array([len(s) for s in union_sets])
 
diff --git a/bct/algorithms/core.py b/bct/algorithms/core.py
index 2bba5a0..fb0dcef 100644
--- a/bct/algorithms/core.py
+++ b/bct/algorithms/core.py
@@ -141,8 +141,8 @@ def assortativity_wei(CIJ, flag=0):
 
 
 def core_periphery_dir(W, gamma=1, C0=None, seed=None):
-    ''' 
-    The optimal core/periphery subdivision is a partition of the network 
+    '''
+    The optimal core/periphery subdivision is a partition of the network
     into two nonoverlapping groups of nodes, a core group and a periphery
     group. The number of core-group edges is maximized, and the number of
     within periphery edges is minimized.
@@ -166,16 +166,15 @@ def core_periphery_dir(W, gamma=1, C0=None, seed=None):
         0 < gamma < 1 detects large core, small periphery
     C0 : NxN np.ndarray
         Initial core structure
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
     '''
     rng = get_rng(seed)
     n = len(W)
     np.fill_diagonal(W, 0)
 
     if C0 == None:
-        C = rng.randint(2, size=(n,))
+        C = rng.integers(2, size=(n,))
     else:
         C = C0.copy()
 
@@ -195,20 +194,20 @@ def core_periphery_dir(W, gamma=1, C0=None, seed=None):
     flag = True
     it = 0
     while flag:
-        it += 1  
+        it += 1
         if it > 100:
             raise BCTParamError('Infinite Loop aborted')
 
         flag = False
         #initial node indices
-        ixes = np.arange(n)    
+        ixes = np.arange(n)
 
         Ct = C.copy()
         while len(ixes) > 0:
             Qt = np.zeros((n,))
             ctix, = np.where(Ct)
             nctix, = np.where(np.logical_not(Ct))
-            q0 = (np.sum(B[np.ix_(ctix, ctix)]) - 
+            q0 = (np.sum(B[np.ix_(ctix, ctix)]) -
                   np.sum(B[np.ix_(nctix, nctix)]))
             Qt[ctix] = q0 - 2 * np.sum(B[ctix, :], axis=1)
             Qt[nctix] = q0 + 2 * np.sum(B[nctix, :], axis=1)
@@ -216,18 +215,18 @@ def core_periphery_dir(W, gamma=1, C0=None, seed=None):
             max_Qt = np.max(Qt[ixes])
             u, = np.where(np.abs(Qt[ixes]-max_Qt) < 1e-10)
             #tunourn
-            u = u[rng.randint(len(u))]
+            u = u[rng.integers(len(u))]
             Ct[ixes[u]] = np.logical_not(Ct[ixes[u]])
             #casga
 
             ixes = np.delete(ixes, u)
-            
+
             if max_Qt - q > 1e-10:
                 flag = True
                 C = Ct.copy()
                 cix, = np.where(C)
                 ncix, = np.where(np.logical_not(C))
-                q = (np.sum(B[np.ix_(cix, cix)]) - 
+                q = (np.sum(B[np.ix_(cix, cix)]) -
                      np.sum(B[np.ix_(ncix, ncix)]))
 
     cix, = np.where(C)
@@ -384,7 +383,7 @@ def local_assortativity_wu_sign(W):
     ----------
     W : NxN np.ndarray
         undirected connection matrix with positive and negative weights
-    
+
     Returns
     -------
     loc_assort_pos : Nx1 np.ndarray
@@ -405,13 +404,13 @@ def local_assortativity_wu_sign(W):
 
     for curr_node in range(n):
         jp = np.where(W[curr_node, :] > 0)
-        loc_assort_pos[curr_node] = np.sum(np.abs(str_pos[jp] - 
+        loc_assort_pos[curr_node] = np.sum(np.abs(str_pos[jp] -
             str_pos[curr_node])) / str_pos[curr_node]
         jn = np.where(W[curr_node, :] < 0)
         loc_assort_neg[curr_node] = np.sum(np.abs(str_neg[jn] -
             str_neg[curr_node])) / str_neg[curr_node]
 
-    loc_assort_pos = ((r_pos + 1) / n - 
+    loc_assort_pos = ((r_pos + 1) / n -
         loc_assort_pos / np.sum(loc_assort_pos))
     loc_assort_neg = ((r_neg + 1) / n -
         loc_assort_neg / np.sum(loc_assort_neg))
diff --git a/bct/algorithms/generative.py b/bct/algorithms/generative.py
index fe5c805..f5b46ed 100644
--- a/bct/algorithms/generative.py
+++ b/bct/algorithms/generative.py
@@ -10,7 +10,7 @@
 
 
 @due.dcite(BibTeX(BETZEL2016), description="Generative models")
-def generative_model(A, D, m, eta, gamma=None, model_type='matching', 
+def generative_model(A, D, m, eta, gamma=None, model_type='matching',
     model_var='powerlaw', epsilon=1e-6, copy=True, seed=None):
     '''
     Generates synthetic networks using the models described in
@@ -32,10 +32,10 @@ def generative_model(A, D, m, eta, gamma=None, model_type='matching',
     D : np.ndarray
         Matrix of euclidean distances or other distances between nodes
     m : int
-        Number of connections that should be present in the final synthetic 
+        Number of connections that should be present in the final synthetic
         network
     eta : np.ndarray
-        A vector describing a range of values to estimate for eta, the 
+        A vector describing a range of values to estimate for eta, the
         hyperparameter describing exponential weighting of the euclidean
         distance.
     gamma : np.ndarray
@@ -44,7 +44,7 @@ def generative_model(A, D, m, eta, gamma=None, model_type='matching',
         algorithm. If model_type='euclidean' or another distance metric,
         this can be None.
     model_type : Enum(str)
-        euclidean : Uses only euclidean distances to generate connection 
+        euclidean : Uses only euclidean distances to generate connection
             probabilities
         neighbors : count of common neighbors
         matching : matching index, the normalized overlap in neighborhoods
@@ -67,16 +67,15 @@ def generative_model(A, D, m, eta, gamma=None, model_type='matching',
     copy : bool
         Some algorithms add edges directly to the input matrix. Set this flag
         to make a copy of the input matrix instead. Defaults to True.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
     '''
     rng = get_rng(seed)
     if copy:
         A = A.copy()
 
     n = len(D)
-    
+
     #These parameters don't do any of the voronoi narrowing.
     #Its a list of eta values paired with gamma values.
     #To try 3 eta and 3 gamma pairs, should use 9 list values.
@@ -93,7 +92,7 @@ def k_avg(K):
             epsilon)
 
     def k_diff(K):
-        return np.abs(np.tile(K, (n, 1)) - 
+        return np.abs(np.tile(K, (n, 1)) -
                       np.transpose(np.tile(K, (n, 1)))) + epsilon
 
     def k_max(K):
@@ -117,7 +116,7 @@ def s_diff(K, sc):
 
     def s_min(K, sc):
         return np.where(K < sc, K + epsilon, sc + epsilon)
-    
+
     def s_max(K, sc):
         #return np.max((K, sc.T), axis=0)
         return np.where(K > sc, K + epsilon, sc + epsilon)
@@ -169,12 +168,12 @@ def clu_gen(A, K, D, m, eta, gamma, model_var, x_fun):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         c = clustering_coef_bu(A)
         k = np.sum(A, axis=1)
@@ -204,7 +203,7 @@ def clu_gen(A, K, D, m, eta, gamma, model_var, x_fun):
             c[k<=1] = 0
             bth[uu] = 1
             bth[vv] = 1
-    
+
             k_result = x_fun(c, bth)
 
             #print(np.shape(k_result))
@@ -239,12 +238,12 @@ def deg_gen(A, K, D, m, eta, gamma, model_var, s_fun):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         P = Fd * Fk * np.logical_not(A)
         u,v = np.where(np.triu(np.ones((n,n)), 1))
@@ -258,7 +257,7 @@ def deg_gen(A, K, D, m, eta, gamma, model_var, s_fun):
 #        print(np.shape(np.where(A[u,v])), 'sqishy')
 #        print(np.shape(P), 'squnnaq')
 
-        #b[:mseed] = np.where(A[np.ix_(u,v)]) 
+        #b[:mseed] = np.where(A[np.ix_(u,v)])
         b[:mseed] = np.squeeze(np.where(A[u,v]))
         #print(mseed, m)
         for i in range(mseed, m):
@@ -318,16 +317,16 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         Ff = Fd * Fk * np.logical_not(A)
         u,v = np.where(np.triu(np.ones((n,n)), 1))
-    
+
         for ii in range(mseed, m):
             C = np.append(0, np.cumsum(Ff[u,v]))
             r = np.sum(rng.random_sample()*C[-1] >= C)
@@ -343,7 +342,7 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
             for i in range(len(updateuu)):
                 j = updateuu[i]
                 c2 = np.append(A[:,j], A[j,:])
-    
+
                 use = np.logical_or(c1, c2)
                 use[uu] = use[uu+n] = use[j] = use[j+n] = 0
                 ncon = np.sum(c1[use]) + np.sum(c2[use])
@@ -356,12 +355,12 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
             updatevv, = np.where(np.inner(A, A[:,vv]))
             np.delete(updatevv, np.where(updatevv == uu))
             np.delete(updatevv, np.where(updatevv == vv))
-        
+
             c1 = np.append(A[:,vv], A[vv,:])
             for i in range(len(updatevv)):
                 j = updatevv[i]
                 c2 = np.append(A[:,j], A[j,:])
-    
+
                 use = np.logical_or(c1, c2)
                 use[vv] = use[vv+n] = use[j] = use[j+n] = 0
                 ncon = np.sum(c1[use]) + np.sum(c2[use])
@@ -374,7 +373,7 @@ def matching_gen(A, K, D, m, eta, gamma, model_var):
             Ff = Fd * Fk * np.logical_not(A)
 
         return A
-    
+
     def neighbors_gen(A, K, D, m, eta, gamma, model_var):
         K += epsilon
 
@@ -388,16 +387,16 @@ def neighbors_gen(A, K, D, m, eta, gamma, model_var):
         if mv1 in ('powerlaw', 'power_law'):
             Fd = D**eta
         elif mv1 in ('exponential',):
-            Fd = np.exp(eta*D) 
+            Fd = np.exp(eta*D)
 
         if mv2 in ('powerlaw', 'power_law'):
             Fk = K**gamma
         elif mv2 in ('exponential',):
-            Fk = np.exp(gamma*K) 
+            Fk = np.exp(gamma*K)
 
         Ff = Fd * Fk * np.logical_not(A)
         u,v = np.where(np.triu(np.ones((n,n)), 1))
-    
+
         for ii in range(mseed, m):
             C = np.append(0, np.cumsum(Ff[u,v]))
             r = np.sum(rng.random_sample()*C[-1] >= C)
@@ -407,7 +406,7 @@ def neighbors_gen(A, K, D, m, eta, gamma, model_var):
 
             x = A[uu, :].astype(int)
             y = A[:, vv].astype(int)
-    
+
             K[uu, y] += 1
             K[y, uu] += 1
             K[vv, x] += 1
@@ -416,7 +415,7 @@ def neighbors_gen(A, K, D, m, eta, gamma, model_var):
             if mv2 in ('powerlaw', 'power_law'):
                 Fk = K**gamma
             elif mv2 in ('exponential',):
-                Fk = np.exp(gamma*K) 
+                Fk = np.exp(gamma*K)
 
             if mv2 in ('powerlaw', 'power_law'):
                 Ff[uu, y] = Ff[y, uu] = Fd[uu, y] * (K[uu, y] ** gamma)
@@ -474,12 +473,12 @@ def euclidean_gen(A, D, m, eta, model_var):
     elif model_type in ('clu-max', 'clu_max'):
         Kseed = k_max(clustering_coef_bu(A))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
-            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_max) 
+            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_max)
 
     elif model_type in ('clu-min', 'clu_min'):
         Kseed = k_min(clustering_coef_bu(A))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
-            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_min) 
+            B[:,:,j] = clu_gen(A, Kseed, D, m, ep, gp, model_var, x_min)
 
     elif model_type in ('clu-prod', 'clu_prod'):
         Kseed = k_prod(clustering_coef_bu(A))
@@ -495,7 +494,7 @@ def euclidean_gen(A, D, m, eta, model_var):
         Kseed = k_diff(np.sum(A, axis=1))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
             B[:,:,j] = deg_gen(A, Kseed, D, m, ep, gp, model_var, s_diff)
-    
+
     elif model_type in ('deg-max', 'deg_max'):
         Kseed = k_max(np.sum(A, axis=1))
         for j, (ep, gp) in enumerate(zip(eta, gamma)):
@@ -525,11 +524,11 @@ def euclidean_gen(A, D, m, eta, model_var):
 
     elif model_type in ('spatial', 'geometric', 'euclidean'):
         for j, ep in enumerate(eta):
-            B[:,:,j] = euclidean_gen(A, D, m, ep, model_var) 
+            B[:,:,j] = euclidean_gen(A, D, m, ep, model_var)
 
     return np.squeeze(B)
 
-def evaluate_generative_model(A, Atgt, D, eta, gamma=None, 
+def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
     model_type='matching', model_var='powerlaw', epsilon=1e-6, seed=None):
     '''
     Generates synthetic networks with parameters provided and evaluates their
@@ -537,7 +536,7 @@ def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
     Basically it takes the Kolmogorov-Smirnov statistics of 4 network
     measures; comparing the degree distributions, clustering coefficients,
     betweenness centrality, and Euclidean distances between connected regions.
-    
+
     The energy is globally low if the synthetic network matches the target.
     Energy is defined as the maximum difference across the four statistics.
     '''
@@ -548,11 +547,11 @@ def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
     xb = betweenness_bin(Atgt)
     xe = D[np.triu(Atgt, 1) > 0]
 
-    B = generative_model(A, D, m, eta, gamma, model_type=model_type, 
+    B = generative_model(A, D, m, eta, gamma, model_type=model_type,
                          model_var=model_var, epsilon=epsilon, copy=True, seed=seed)
 
     #if eta != gamma then an error is thrown within generative model
-    
+
     nB = len(eta)
 
     if nB == 1:
@@ -562,7 +561,7 @@ def evaluate_generative_model(A, Atgt, D, eta, gamma=None,
 
     def kstats(x, y):
         bin_edges = np.concatenate([[-np.inf],
-                                    np.sort(np.concatenate((x, y))), 
+                                    np.sort(np.concatenate((x, y))),
                                     [np.inf]])
 
         bin_x,_ = np.histogram(x, bin_edges)
diff --git a/bct/algorithms/modularity.py b/bct/algorithms/modularity.py
index 90bb8c4..efaec38 100644
--- a/bct/algorithms/modularity.py
+++ b/bct/algorithms/modularity.py
@@ -92,15 +92,14 @@ def community_louvain(W, gamma=1, ci=None, B='modularity', seed=None):
         initial community affiliation vector. default value=None
     B : str | NxN np.arraylike
         string describing objective function type, or provides a custom
-        NxN objective-function matrix. builtin values 
+        NxN objective-function matrix. builtin values
             'modularity' uses Q-metric as objective function
             'potts' uses Potts model Hamiltonian.
             'negative_sym' symmetric treatment of negative weights
             'negative_asym' asymmetric treatment of negative weights
         The default value is to use the Q-metric
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -167,7 +166,7 @@ def community_louvain(W, gamma=1, ci=None, B='modularity', seed=None):
             print ('Warning: objective function matrix not symmetric, '
                    'symmetrizing')
             B = (B + B.T) / 2
-    
+
     Hnm = np.zeros((n, n))
     for m in range(1, n + 1):
         Hnm[:, m - 1] = np.sum(B[:, ci == m], axis=1)  # node to module degree
@@ -236,7 +235,7 @@ def community_louvain(W, gamma=1, ci=None, B='modularity', seed=None):
         q0 = q
 
         q = np.trace(B)  # compute modularity
-    
+
     # Workaround to normalize
     if not renormalize:
         return ci, q/s
@@ -607,9 +606,8 @@ def modularity_finetune_dir(W, ci=None, gamma=1, seed=None):
     gamma : float
         resolution parameter. default value=1. Values 0 <= gamma < 1 detect
         larger modules while gamma > 1 detects smaller modules.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -711,9 +709,8 @@ def modularity_finetune_und(W, ci=None, gamma=1, seed=None):
     gamma : float
         resolution parameter. default value=1. Values 0 <= gamma < 1 detect
         larger modules while gamma > 1 detects smaller modules.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -815,9 +812,8 @@ def modularity_finetune_und_sign(W, qtype='sta', gamma=1, ci=None, seed=None):
         larger modules while gamma > 1 detects smaller modules.
     ci : Nx1 np.ndarray | None
         initial community affiliation vector
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -949,9 +945,8 @@ def modularity_louvain_dir(W, gamma=1, hierarchy=False, seed=None):
         larger modules while gamma > 1 detects smaller modules.
     hierarchy : bool
         Enables hierarchical output. Defalut value=False
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1082,9 +1077,8 @@ def modularity_louvain_und(W, gamma=1, hierarchy=False, seed=None):
         larger modules while gamma > 1 detects smaller modules.
     hierarchy : bool
         Enables hierarchical output. Defalut value=False
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1222,9 +1216,8 @@ def modularity_louvain_und_sign(W, gamma=1, qtype='sta', seed=None):
     gamma : float
         resolution parameter. default value=1. Values 0 <= gamma < 1 detect
         larger modules while gamma > 1 detects smaller modules.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1392,9 +1385,8 @@ def modularity_probtune_und_sign(W, qtype='sta', gamma=1, ci=None, p=.45,
         initial community affiliation vector
     p : float
         probability of random node moves. Default value = 0.45
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1460,9 +1452,9 @@ def modularity_probtune_und_sign(W, qtype='sta', gamma=1, ci=None, p=.45,
 
     for u in rng.permutation(n):  # loop over nodes in random order
         ma = ci[u] - 1  # current module
-        r = rng.random_sample() < p
+        r = rng.random() < p
         if r:
-            mb = rng.randint(n)  # select new module randomly
+            mb = rng.integers(n)  # select new module randomly
         else:
             dq0 = ((Knm0[u, :] + W0[u, u] - Knm0[u, ma]) -
                    gamma * Kn0[u] * (Km0 + Kn0[u] - Km0[ma]) / s0)
diff --git a/bct/algorithms/physical_connectivity.py b/bct/algorithms/physical_connectivity.py
index 55cf7cf..aecaf27 100644
--- a/bct/algorithms/physical_connectivity.py
+++ b/bct/algorithms/physical_connectivity.py
@@ -93,9 +93,8 @@ def rentian_scaling(A, xyz, n, seed=None):
     n : int
         Number of partitions to compute. Each partition is a data point; you
         want a large enough number to adequately compute Rent's exponent.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -141,7 +140,7 @@ def rentian_scaling(A, xyz, n, seed=None):
     # and the number of edges traversing the boundary of the partition (e)
     while count < n:
         # define cube endpoints
-        randx = np.sort((1 + nmax - nmin) * rng.random_sample((2,)))
+        randx = np.sort((1 + nmax - nmin) * rng.random((2,)))
 
         # find nodes in cube
         l1 = xyzn[:, 0] > randx[0]
diff --git a/bct/algorithms/reference.py b/bct/algorithms/reference.py
index a259e55..77876ac 100644
--- a/bct/algorithms/reference.py
+++ b/bct/algorithms/reference.py
@@ -27,9 +27,8 @@ def latmio_dir_connected(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -149,9 +148,8 @@ def latmio_dir(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -196,10 +194,10 @@ def latmio_dir(R, itr, D=None, seed=None):
         att = 0
         while att <= max_attempts:  # while not rewired
             while True:
-                e1 = rng.randint(k)
-                e2 = rng.randint(k)
+                e1 = rng.integers(k)
+                e2 = rng.integers(k)
                 while e1 == e2:
-                    e2 = rng.randint(k)
+                    e2 = rng.integers(k)
                 a = i[e1]
                 b = j[e1]
                 c = i[e2]
@@ -249,9 +247,8 @@ def latmio_und_connected(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -302,10 +299,10 @@ def latmio_und_connected(R, itr, D=None, seed=None):
         while att <= max_attempts:
             rewire = True
             while True:
-                e1 = rng.randint(k)
-                e2 = rng.randint(k)
+                e1 = rng.integers(k)
+                e2 = rng.integers(k)
                 while e1 == e2:
-                    e2 = rng.randint(k)
+                    e2 = rng.integers(k)
                 a = i[e1]
                 b = j[e1]
                 c = i[e2]
@@ -314,7 +311,7 @@ def latmio_und_connected(R, itr, D=None, seed=None):
                 if a != c and a != d and b != c and b != d:
                     break
 
-            if rng.random_sample() > .5:
+            if rng.random() > .5:
                 i.setflags(write=True)
                 j.setflags(write=True)
                 i[e2] = d
@@ -384,9 +381,8 @@ def latmio_und(R, itr, D=None, seed=None):
     D : np.ndarray | None
         distance-to-diagonal matrix. Defaults to the actual distance matrix
         if not specified.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -430,10 +426,10 @@ def latmio_und(R, itr, D=None, seed=None):
         att = 0
         while att <= max_attempts:
             while True:
-                e1 = rng.randint(k)
-                e2 = rng.randint(k)
+                e1 = rng.integers(k)
+                e2 = rng.integers(k)
                 while e1 == e2:
-                    e2 = rng.randint(k)
+                    e2 = rng.integers(k)
                 a = i[e1]
                 b = j[e1]
                 c = i[e2]
@@ -442,7 +438,7 @@ def latmio_und(R, itr, D=None, seed=None):
                 if a != c and a != d and b != c and b != d:
                     break
 
-            if rng.random_sample() > .5:
+            if rng.random() > .5:
                 i.setflags(write=True)
                 j.setflags(write=True)
                 i[e2] = d
@@ -488,9 +484,8 @@ def makeevenCIJ(n, k, sz_cl, seed=None):
         number of edges
     sz_cl : int
         size of clusters (must be power of 2)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -566,9 +561,8 @@ def makefractalCIJ(mx_lvl, E, sz_cl, seed=None):
         connection density fall off per level
     sz_cl : int
         size of clusters (must be power of 2)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -603,7 +597,7 @@ def makefractalCIJ(mx_lvl, E, sz_cl, seed=None):
     ee = mx_lvl - CIJ - sz_cl
     ee = (ee > 0) * ee
     prob = (1 / E**ee) * (np.ones((s, s)) - np.eye(s))
-    CIJ = (prob > rng.random_sample((n, n)))
+    CIJ = (prob > rng.random((n, n)))
 
     # count connections
     k = np.sum(CIJ)
@@ -622,9 +616,8 @@ def makerandCIJdegreesfixed(inv, outv, seed=None):
         in-degree vector
     outv : Nx1 np.ndarray
         out-degree vector
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -674,9 +667,9 @@ def makerandCIJdegreesfixed(inv, outv, seed=None):
                 if len(tried) == k:
                     raise BCTParamError('Could not resolve the given '
                                         'in and out vectors')
-                switch = rng.randint(k)
+                switch = rng.integers(k)
                 while switch in tried:
-                    switch = rng.randint(k)
+                    switch = rng.integers(k)
                 if not (CIJ[edges[0, i], edges[1, switch]] or
                         CIJ[edges[0, switch], edges[1, i]]):
                     CIJ[edges[0, switch], edges[1, switch]] = 0
@@ -706,9 +699,8 @@ def makerandCIJ_dir(n, k, seed=None):
         number of vertices
     K : int
         number of edges
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -738,9 +730,8 @@ def makerandCIJ_und(n, k, seed=None):
         number of vertices
     K : int
         number of edges
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -771,9 +762,8 @@ def makeringlatticeCIJ(n, k, seed=None):
         number of vertices
     K : int
         number of edges
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -829,9 +819,8 @@ def maketoeplitzCIJ(n, k, s, seed=None):
         number of edges
     s : float
         standard deviation of toeplitz
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -851,7 +840,7 @@ def maketoeplitzCIJ(n, k, s, seed=None):
     CIJ = np.zeros((n, n))
     itr = 0
     while np.sum(CIJ) != k:
-        CIJ = (rng.random_sample((n, n)) < template)
+        CIJ = (rng.random((n, n)) < template)
         itr += 1
         if itr > 10000:
             raise BCTParamError('Infinite loop was caught generating toeplitz '
@@ -881,9 +870,8 @@ def null_model_dir_sign(W, bin_swaps=5, wei_freq=.1, seed=None):
         wei_freq == 0.1 implies that weights sorted each 10th step (faster,
             default value)
         wei_freq == 0 implies no sorting of weights (not recommended)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1005,9 +993,8 @@ def null_model_und_sign(W, bin_swaps=5, wei_freq=.1, seed=None):
         wei_freq == 0.1 implies that weights sorted each 10th step (faster,
             default value)
         wei_freq == 0 implies no sorting of weights (not recommended)
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1131,9 +1118,8 @@ def randmio_dir_connected(R, itr, seed=None):
         directed binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1157,10 +1143,10 @@ def randmio_dir_connected(R, itr, seed=None):
         while att <= max_attempts:  # while not rewired
             rewire = True
             while True:
-                e1 = rng.randint(k)
-                e2 = rng.randint(k)
+                e1 = rng.integers(k)
+                e2 = rng.integers(k)
                 while e1 == e2:
-                    e2 = rng.randint(k)
+                    e2 = rng.integers(k)
                 a = i[e1]
                 b = j[e1]
                 c = i[e2]
@@ -1223,9 +1209,8 @@ def randmio_dir(R, itr, seed=None):
         directed binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1248,10 +1233,10 @@ def randmio_dir(R, itr, seed=None):
         att = 0
         while att <= max_attempts:  # while not rewired
             while True:
-                e1 = rng.randint(k)
-                e2 = rng.randint(k)
+                e1 = rng.integers(k)
+                e2 = rng.integers(k)
                 while e1 == e2:
-                    e2 = rng.randint(k)
+                    e2 = rng.integers(k)
                 a = i[e1]
                 b = j[e1]
                 c = i[e2]
@@ -1300,9 +1285,8 @@ def randmio_und_connected(R, itr, seed=None):
         undirected binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1334,10 +1318,10 @@ def randmio_und_connected(R, itr, seed=None):
         while att <= max_attempts:  # while not rewired
             rewire = True
             while True:
-                e1 = rng.randint(k)
-                e2 = rng.randint(k)
+                e1 = rng.integers(k)
+                e2 = rng.integers(k)
                 while e1 == e2:
-                    e2 = rng.randint(k)
+                    e2 = rng.integers(k)
                 a = i[e1]
                 b = j[e1]
                 c = i[e2]
@@ -1346,7 +1330,7 @@ def randmio_und_connected(R, itr, seed=None):
                 if a != c and a != d and b != c and b != d:
                     break  # all 4 vertices must be different
 
-            if rng.random_sample() > .5:
+            if rng.random() > .5:
 
                 i.setflags(write=True)
                 j.setflags(write=True)
@@ -1412,9 +1396,8 @@ def randmio_dir_signed(R, itr, seed=None):
         directed binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1488,9 +1471,8 @@ def randmio_und(R, itr, seed=None):
         undirected binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1517,9 +1499,9 @@ def randmio_und(R, itr, seed=None):
         att = 0
         while att <= max_attempts:  # while not rewired
             while True:
-                e1, e2 = rng.randint(k, size=(2,))
+                e1, e2 = rng.integers(k, size=(2,))
                 while e1 == e2:
-                    e2 = rng.randint(k)
+                    e2 = rng.integers(k)
                 a = i[e1]
                 b = j[e1]
                 c = i[e2]
@@ -1528,7 +1510,7 @@ def randmio_und(R, itr, seed=None):
                 if a != c and a != d and b != c and b != d:
                     break  # all 4 vertices must be different
 
-            if rng.random_sample() > .5:
+            if rng.random() > .5:
                 i.setflags(write=True)
                 j.setflags(write=True)
                 i[e2] = d
@@ -1571,9 +1553,8 @@ def randmio_und_signed(R, itr, seed=None):
         undirected binary/weighted connection matrix
     itr : int
         rewiring parameter. Each edge is rewired approximately itr times.
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1636,9 +1617,8 @@ def randomize_graph_partial_und(A, B, maxswap, seed=None):
         mask; edges to avoid
     maxswap : int
         number of rewirings
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1663,9 +1643,9 @@ def randomize_graph_partial_und(A, B, maxswap, seed=None):
     nswap = 0
     while nswap < maxswap:
         while True:
-            e1, e2 = rng.randint(m, size=(2,))
+            e1, e2 = rng.integers(m, size=(2,))
             while e1 == e2:
-                e2 = rng.randint(m)
+                e2 = rng.integers(m)
             a = i[e1]
             b = j[e1]
             c = i[e2]
@@ -1674,7 +1654,7 @@ def randomize_graph_partial_und(A, B, maxswap, seed=None):
             if a != c and a != d and b != c and b != d:
                 break  # all 4 vertices must be different
 
-        if rng.random_sample() > .5:
+        if rng.random() > .5:
             i[e2] = d
             j[e2] = c  # flip edge c-d with 50% probability
             c = i[e2]
@@ -1711,9 +1691,8 @@ def randomizer_bin_und(R, alpha, seed=None):
         binary undirected connection matrix
     alpha : float
         fraction of edges to rewire
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -1761,7 +1740,7 @@ def randomizer_bin_und(R, alpha, seed=None):
         raise BCTParamError("No possible randomization")
 
     for it in range(k):
-        if rng.random_sample() > alpha:
+        if rng.random() > alpha:
             continue  # rewire alpha% of edges
 
         a = i[it]
@@ -1779,10 +1758,10 @@ def randomizer_bin_und(R, alpha, seed=None):
         if np.size(ii):
             # choose one randomly
             nummates = np.size(ii)
-            mate = rng.randint(nummates)
+            mate = rng.integers(nummates)
 
             # randomly orient the second edge
-            if rng.random_sample() > .5:
+            if rng.random() > .5:
                 c = i_intersect[ii[mate]]
                 d = i_intersect[jj[mate]]
             else:
diff --git a/bct/nbs.py b/bct/nbs.py
index ce5006c..be397fd 100644
--- a/bct/nbs.py
+++ b/bct/nbs.py
@@ -27,7 +27,7 @@ def nbs_bct(x, y, thresh, k=1000, tail='both', paired=False, verbose=False, seed
     thresh : float
         minimum t-value used as threshold
     k : int
-        number of permutations used to estimate the empirical null 
+        number of permutations used to estimate the empirical null
         distribution
     tail : {'left', 'right', 'both'}
         enables specification of particular alternative hypothesis
@@ -39,9 +39,8 @@ def nbs_bct(x, y, thresh, k=1000, tail='both', paired=False, verbose=False, seed
         subject populations to have equal N. default value = False
     verbose : bool
         print some extra information each iteration. defaults value = False
-    seed : hashable, optional
-        If None (default), use the np.random's global random state to generate random numbers.
-        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+    seed : None, int, or numpy.random.Generator
+        Seed (or RNG itself) used to generate random numbers.
 
     Returns
     -------
@@ -55,42 +54,42 @@ def nbs_bct(x, y, thresh, k=1000, tail='both', paired=False, verbose=False, seed
         an adjacency matrix identifying the edges comprising each component.
         edges are assigned indexed values.
     null : Kx1 np.ndarray
-        A vector of K sampled from the null distribution of maximal component 
+        A vector of K sampled from the null distribution of maximal component
         size.
 
     Notes
     -----
-    ALGORITHM DESCRIPTION 
-    The NBS is a nonparametric statistical test used to isolate the 
-    components of an N x N undirected connectivity matrix that differ 
-    significantly between two distinct populations. Each element of the 
-    connectivity matrix stores a connectivity value and each member of 
-    the two populations possesses a distinct connectivity matrix. A 
-    component of a connectivity matrix is defined as a set of 
-    interconnected edges. 
-
-    The NBS is essentially a procedure to control the family-wise error 
-    rate, in the weak sense, when the null hypothesis is tested 
+    ALGORITHM DESCRIPTION
+    The NBS is a nonparametric statistical test used to isolate the
+    components of an N x N undirected connectivity matrix that differ
+    significantly between two distinct populations. Each element of the
+    connectivity matrix stores a connectivity value and each member of
+    the two populations possesses a distinct connectivity matrix. A
+    component of a connectivity matrix is defined as a set of
+    interconnected edges.
+
+    The NBS is essentially a procedure to control the family-wise error
+    rate, in the weak sense, when the null hypothesis is tested
     independently at each of the N(N-1)/2 edges comprising the undirected
-    connectivity matrix. The NBS can provide greater statistical power 
-    than conventional procedures for controlling the family-wise error 
+    connectivity matrix. The NBS can provide greater statistical power
+    than conventional procedures for controlling the family-wise error
     rate, such as the false discovery rate, if the set of edges at which
     the null hypothesis is rejected constitues a large component or
     components.
     The NBS comprises fours steps:
     1. Perform a two-sample T-test at each edge indepedently to test the
        hypothesis that the value of connectivity between the two
-       populations come from distributions with equal means. 
+       populations come from distributions with equal means.
     2. Threshold the T-statistic available at each edge to form a set of
-       suprathreshold edges. 
+       suprathreshold edges.
     3. Identify any components in the adjacency matrix defined by the set
-       of suprathreshold edges. These are referred to as observed 
-       components. Compute the size of each observed component 
-       identified; that is, the number of edges it comprises. 
+       of suprathreshold edges. These are referred to as observed
+       components. Compute the size of each observed component
+       identified; that is, the number of edges it comprises.
     4. Repeat K times steps 1-3, each time randomly permuting members of
-       the two populations and storing the size of the largest component 
+       the two populations and storing the size of the largest component
        identified for each permuation. This yields an empirical estimate
-       of the null distribution of maximal component size. A corrected 
+       of the null distribution of maximal component size. A corrected
        p-value for each observed component is then calculated using this
        null distribution.
 
diff --git a/bct/utils/miscellaneous_utilities.py b/bct/utils/miscellaneous_utilities.py
index b4a9f36..9f7f86a 100644
--- a/bct/utils/miscellaneous_utilities.py
+++ b/bct/utils/miscellaneous_utilities.py
@@ -1,5 +1,4 @@
 from __future__ import division, print_function
-import random
 import numpy as np
 
 
@@ -26,7 +25,7 @@ def pick_four_unique_nodes_quickly(n, seed=None):
     clever but still substantially slower.
     '''
     rng = get_rng(seed)
-    k = rng.randint(n**4)
+    k = rng.integers(n**4)
     a = k % n
     b = k // n % n
     c = k // n ** 2 % n
@@ -87,28 +86,18 @@ def dummyvar(cis, return_sparse=False):
 
 
 def get_rng(seed=None):
-    """
-    By default, or if `seed` is np.random, return the global RandomState
-    instance used by np.random.
-    If `seed` is a RandomState instance, return it unchanged.
-    Otherwise, use the passed (hashable) argument to seed a new instance
-    of RandomState and return it.
+    """Return a random number generator.
+
+    This method is a simple wrapper over ``numpy.random.default_rng``,
+    used for legacy compatibility purposes;
+    see that function's docs for more information.
 
     Parameters
     ----------
-    seed : hashable or np.random.RandomState or np.random, optional
+    seed : None, int, array_like[ints], SeedSequence, BitGenerator, Generator
 
     Returns
     -------
-    np.random.RandomState
+    np.random.Generator
     """
-    if seed is None or seed == np.random:
-        return np.random.mtrand._rand
-    elif isinstance(seed, np.random.RandomState):
-        return seed
-    try:
-        rstate =  np.random.RandomState(seed)
-    except ValueError:
-        rstate = np.random.RandomState(random.Random(seed).randint(0, 2**32-1))
-    return rstate
-
+    return np.random.default_rng(seed)
diff --git a/test/conftest.py b/test/conftest.py
index c9bf6cd..c59f1c9 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -1,4 +1,5 @@
 import pytest
+import numpy as np
 
 
 def pytest_configure(config):
@@ -20,3 +21,15 @@ def pytest_collection_modifyitems(config, items):
     for item in items:
         if "long" in item.keywords:
             item.add_marker(skip)
+
+
+@pytest.fixture
+def stable_rng():
+    """
+    To make use of algorithmic improvements in numpy's random number generators,
+    bctpy uses numpy.random.default_rng() when creating RNGs.
+    As such, results are not guaranteed to be deterministic between numpy versions.
+    Therefore, for testing purposes we explicitly set the BitGenerator class.
+    """
+    seed = 1991
+    return np.random.Generator(np.random.PCG64(seed))
diff --git a/test/modularity_test.py b/test/modularity_test.py
index 44a788e..12b18e4 100644
--- a/test/modularity_test.py
+++ b/test/modularity_test.py
@@ -15,12 +15,11 @@ def test_modularity_und():
     # package numerical instability of eigendecompositions
 
 
-def test_modularity_louvain_und():
+def test_modularity_louvain_und(stable_rng):
     x = load_sample(thres=.4)
 
-    seed = 38429004
-    _, q = bct.modularity_louvain_und(x, seed=seed)
-    assert np.allclose(q, 0.25892588)
+    _, q = bct.modularity_louvain_und(x, seed=stable_rng)
+    assert np.allclose(q, 0.253_148_271)
 
     fails = 0
     for i in range(100):
@@ -33,17 +32,15 @@ def test_modularity_louvain_und():
             else:
                 fails += 1
 
-    seed = 94885236
-    _, q = bct.modularity_finetune_und(x, seed=seed)
-    assert np.allclose(q, .25879794)
+    _, q = bct.modularity_finetune_und(x, seed=stable_rng)
+    assert np.allclose(q, 0.242_632_732)
 
 
-def test_modularity_finetune_und():
+def test_modularity_finetune_und(stable_rng):
     x = load_sample(thres=.4)
 
-    seed = 94885236
-    _, q = bct.modularity_finetune_und(x, seed=seed)
-    assert np.allclose(q, .25879794)
+    _, q = bct.modularity_finetune_und(x, seed=stable_rng)
+    assert np.allclose(q, 0.253_148_271)
 
     fails = 0
     for i in range(100):
@@ -56,9 +53,8 @@ def test_modularity_finetune_und():
             else:
                 fails += 1
 
-    seed = 71040925
-    ci, oq = bct.modularity_louvain_und(x, seed=seed)
-    _, q = bct.modularity_finetune_und(x, ci=ci, seed=seed)
+    ci, oq = bct.modularity_louvain_und(x, seed=stable_rng)
+    _, q = bct.modularity_finetune_und(x, ci=ci, seed=stable_rng)
     print(q, oq)
     # assert np.allclose(q, .25892588)
     assert np.allclose(q, .25856714)
@@ -86,48 +82,44 @@ def test_modularity_finetune_und():
     #(i.e. it is unstable both when the modular structure is identical and not)
 
 
-def test_modularity_louvain_und_sign_seed():
+def test_modularity_louvain_und_sign_seed(stable_rng):
     # performance is same as matlab if randomness is quashed
     x = load_signed_sample()
-    seed = 90772777
-    _, q = bct.modularity_louvain_und_sign(x, seed=seed)
+    _, q = bct.modularity_louvain_und_sign(x, seed=stable_rng)
     print(q)
-    assert np.allclose(q, .46605515)
+    assert np.allclose(q, 0.453_375_647)
 
 
-def test_modularity_finetune_und_sign_actually_finetune():
+def test_modularity_finetune_und_sign_actually_finetune(stable_rng):
     x = load_signed_sample()
-    seed = 34908314
-    ci, oq = bct.modularity_louvain_und_sign(x, seed=seed)
-    _, q = bct.modularity_finetune_und_sign(x, seed=seed, ci=ci)
+    ci, oq = bct.modularity_louvain_und_sign(x, seed=stable_rng)
+    _, q = bct.modularity_finetune_und_sign(x, seed=stable_rng, ci=ci)
     print(q)
-    assert np.allclose(q, .47282924)
+    assert np.allclose(q, 0.462_292_161)
     assert q >= oq
 
-    seed = 88215881
-    np.random.seed(seed)
-    randomized_sample = np.random.random_sample(size=(len(x), len(x)))
+    randomized_sample = stable_rng.random(size=(len(x), len(x)))
     randomized_sample = randomized_sample + randomized_sample.T
     x[np.where(bct.threshold_proportional(randomized_sample, .2))] = 0
 
-    ci, oq = bct.modularity_louvain_und_sign(x, seed=seed)
+    ci, oq = bct.modularity_louvain_und_sign(x, seed=stable_rng)
     print(oq)
-    assert np.allclose(oq, .45254522)
+    assert np.allclose(oq, 0.458_063_807)
     for i in range(100):
         _, q = bct.modularity_finetune_und_sign(x, ci=ci)
         assert q >= oq
 
 
-def test_modularity_probtune_und_sign():
+def test_modularity_probtune_und_sign(stable_rng):
     x = load_signed_sample()
-    seed = 59468096
-    ci, q = bct.modularity_probtune_und_sign(x, seed=seed)
+    ci, q = bct.modularity_probtune_und_sign(x, seed=stable_rng)
     print(q)
-    assert np.allclose(q, .07885327)
+    # N.B. this result is quite different
+    # assert np.allclose(q, .07885327)  # legacy numpy.random.RandomState
+    assert np.allclose(q, 0.114_732_792)  # stable_rng
 
-    seed = 1742447
-    ci, _ = bct.modularity_louvain_und_sign(x, seed=seed)
-    _, oq = bct.modularity_finetune_und_sign(x, seed=seed, ci=ci)
+    ci, _ = bct.modularity_louvain_und_sign(x, seed=stable_rng)
+    _, oq = bct.modularity_finetune_und_sign(x, seed=stable_rng, ci=ci)
 
     for i in np.arange(.05, .5, .02):
         fails = 0
@@ -148,11 +140,10 @@ def test_modularity_dir_low_modularity():
     assert np.allclose(q, .06450290)
 
 
-def test_modularity_louvain_dir_low_modularity():
+def test_modularity_louvain_dir_low_modularity(stable_rng):
     x = load_directed_low_modularity_sample(thres=.67)
-    seed = 28917147
-    _, q = bct.modularity_louvain_dir(x, seed=seed)
-    assert np.allclose(q, .06934894)
+    _, q = bct.modularity_louvain_dir(x, seed=stable_rng)
+    assert np.allclose(q, 0.069_334_607)
 
 # def test_modularity_finetune_dir_low_modularity():
 #	x = load_directed_low_modularity_sample(thres=.67)
@@ -173,11 +164,11 @@ def test_modularity_dir():
     assert np.allclose(q, .32742787)
 
 
-def test_modularity_louvain_dir():
+def test_modularity_louvain_dir(stable_rng):
     x = load_directed_sample()
-    seed = 43938304
-    _, q = bct.modularity_louvain_dir(x, seed=seed)
-    assert np.allclose(q, .32697921)
+    _, q = bct.modularity_louvain_dir(x, seed=stable_rng)
+    # assert np.allclose(q, .32697921)  # legacy np.random.RandomState
+    assert np.allclose(q, 0.373_475_890)  # stable_rng
 
 # def test_modularity_finetune_dir():
 #	x = load_directed_sample()
@@ -194,10 +185,9 @@ def test_modularity_louvain_dir():
     # what is wrong
 
 
-def test_community_louvain():
+def test_community_louvain(stable_rng):
     x = load_sample(thres=0.4)
-    seed = 39185
-    ci, q = bct.community_louvain(x, seed=seed)
+    ci, q = bct.community_louvain(x, seed=stable_rng)
     print(q)
     assert np.allclose(q, 0.2583, atol=0.015)
 
diff --git a/tox.ini b/tox.ini
index 79932ad..92bd572 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,5 +1,5 @@
 [tox]
-envlist = py{27,34,35,36,37}
+envlist = py{37,38,39,310}
 
 [testenv]
 deps = -rrequirements.txt
@@ -7,3 +7,10 @@ commands =
     pip install -U pip
     pip install .
     pytest --skiplong -v
+
+[gh]
+python =
+    3.7 = py37
+    3.8 = py38
+    3.9 = py39
+    3.10 = py310