From dda6cc62ee82841c37738331666f959d5d8caa82 Mon Sep 17 00:00:00 2001 From: justinehansen Date: Tue, 13 Jul 2021 08:30:21 -0700 Subject: [PATCH 1/4] [ENH] Added betzel null --- netneurotools/networks.py | 178 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 178 insertions(+) diff --git a/netneurotools/networks.py b/netneurotools/networks.py index 6a36077..7531ff9 100644 --- a/netneurotools/networks.py +++ b/netneurotools/networks.py @@ -376,3 +376,181 @@ def threshold_network(network, retain=10): mst = np.array((mst + mst.T) != 0, dtype=int) return mst + + +def math_length_degree_distribution(W, D, nbins=10, nswap=1000, + replacement=False, weighted=True, + seed=None): + """ + Generates degree- and edge length-preserving surrogate connectomes. + + Parameters + ---------- + W : (N, N) array-like + weighted or binary symmetric connectivity matrix. + D : (N, N) array-like + symmetric distance matrix. + nbins : int + number of distance bins (edge length matrix is performed by swapping + connections in the same bin). Default = 10. + nswap : int + total number of edge swaps to perform. Recommended = nnodes * 20 + Default = 1000. + replacement : bool, optional + if True all the edges are available for swapping. Default= False + weighted : bool, optional + Whether to return weighted rewired connectivity matrix. Default = True + seed : float, optional + Random seed. Default = None + + Returns + ------- + newB : (N, N) array-like + binary rewired matrix + newW: (N, N) array-like + weighted rewired matrix. Returns matrix of zeros of weighted=False. + nr : int + number of successful rewires + + Notes + ----- + Takes a weighted, symmetric connectivity matrix `data` and Euclidean/fiber + length matrix `distance` and generates a randomized network with: + 1. exactly the same degree sequence + 2. approximately the same edge length distribution + 3. exactly the same edge weight distribution + 4. approximately the same weight-length relationship + + Reference + --------- + Betzel, R. F., Bassett, D. S. (2018) Specificity and robustness of + long-distance connections in weighted, interareal connectomes. PNAS. + + """ + + rs = check_random_state(seed) + N = len(W) + # divide the distances by lengths + bins = np.linspace(D[D.nonzero()].min(), D[D.nonzero()].max(), nbins + 1) + bins[-1] += 1 + L = np.zeros((N, N)) + for n in range(nbins): + i, j = np.where(np.logical_and(bins[n] <= D, D < bins[n + 1])) + L[i, j] = n + 1 + + # binarized connectivity + B = (W != 0).astype(np.int_) + + # existing edges (only upper triangular cause it's symmetric) + cn_x, cn_y = np.where(np.triu((B != 0) * B, k=1)) + + tries = 0 + nr = 0 + newB = np.copy(B) + + while((len(cn_x) >= 2) & (nr < nswap)): + # choose randomly the edge to be rewired + r = rs.randint(len(cn_x)) + n_x, n_y = cn_x[r], cn_y[r] + tries += 1 + + # options to rewire with + # connected nodes that doesn't involve (n_x, n_y) + index = (cn_x != n_x) & (cn_y != n_y) & (cn_y != n_x) & (cn_x != n_y) + if(len(np.where(index)[0]) == 0): + cn_x = np.delete(cn_x, r) + cn_y = np.delete(cn_y, r) + + else: + ops1_x, ops1_y = cn_x[index], cn_y[index] + # options that will preserve the distances + # (ops1_x, ops1_y) such that + # L(n_x,n_y) = L(n_x, ops1_x) & L(ops1_x,ops1_y) = L(n_y, ops1_y) + index = (L[n_x, n_y] == L[n_x, ops1_x]) & ( + L[ops1_x, ops1_y] == L[n_y, ops1_y]) + if(len(np.where(index)[0]) == 0): + cn_x = np.delete(cn_x, r) + cn_y = np.delete(cn_y, r) + + else: + ops2_x, ops2_y = ops1_x[index], ops1_y[index] + # options of edges that didn't exist before + index = [(newB[min(n_x, ops2_x[i])][max(n_x, ops2_x[i])] == 0) + & (newB[min(n_y, ops2_y[i])][max(n_y, + ops2_y[i])] == 0) + for i in range(len(ops2_x))] + if(len(np.where(index)[0]) == 0): + cn_x = np.delete(cn_x, r) + cn_y = np.delete(cn_y, r) + + else: + ops3_x, ops3_y = ops2_x[index], ops2_y[index] + + # choose randomly one edge from the final options + r1 = rs.randint(len(ops3_x)) + nn_x, nn_y = ops3_x[r1], ops3_y[r1] + + # Disconnect the existing edges + newB[n_x, n_y] = 0 + newB[nn_x, nn_y] = 0 + # Connect the new edges + newB[min(n_x, nn_x), max(n_x, nn_x)] = 1 + newB[min(n_y, nn_y), max(n_y, nn_y)] = 1 + # one successfull rewire! + nr += 1 + + # rewire with replacement + if replacement: + cn_x[r], cn_y[r] = min(n_x, nn_x), max(n_x, nn_x) + index = np.where((cn_x == nn_x) & (cn_y == nn_y))[0] + cn_x[index], cn_y[index] = min( + n_y, nn_y), max(n_y, nn_y) + # rewire without replacement + else: + cn_x = np.delete(cn_x, r) + cn_y = np.delete(cn_y, r) + index = np.where((cn_x == nn_x) & (cn_y == nn_y))[0] + cn_x = np.delete(cn_x, index) + cn_y = np.delete(cn_y, index) + + if(nr < nswap): + print(f"I didn't finish, out of rewirable edges: {len(cn_x)}") + + i, j = np.triu_indices(N, k=1) + # Make the connectivity matrix symmetric + newB[j, i] = newB[i, j] + + # check the number of edges is preserved + if(len(np.where(B != 0)[0]) != len(np.where(newB != 0)[0])): + print( + f"ERROR --- number of edges changed, \ + B:{len(np.where(B!=0)[0])}, newB:{len(np.where(newB!=0)[0])}") + # check that the degree of the nodes it's the same + for i in range(N): + if(np.sum(B[i]) != np.sum(newB[i])): + print( + f"ERROR --- node {i} changed k by: \ + {np.sum(B[i]) - np.sum(newB[i])}") + + newW = np.zeros((N, N)) + if(weighted): + # Reassign the weights + mask = np.triu(B != 0, k=1) + inids = D[mask] + iniws = W[mask] + inids_index = np.argsort(inids) + # Weights from the shortest to largest edges + iniws = iniws[inids_index] + mask = np.triu(newB != 0, k=1) + finds = D[mask] + i, j = np.where(mask) + # Sort the new edges from the shortest to the largest + finds_index = np.argsort(finds) + i_sort = i[finds_index] + j_sort = j[finds_index] + # Assign the initial sorted weights + newW[i_sort, j_sort] = iniws + # Make it symmetrical + newW[j_sort, i_sort] = iniws + + return newB, newW, nr From 6160a0f6e83816658e6c46462822c8370ad57a5d Mon Sep 17 00:00:00 2001 From: justinehansen Date: Wed, 4 Aug 2021 13:43:42 -0700 Subject: [PATCH 2/4] [FIX] Fix typo in description --- netneurotools/networks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/netneurotools/networks.py b/netneurotools/networks.py index 7531ff9..e602c3b 100644 --- a/netneurotools/networks.py +++ b/netneurotools/networks.py @@ -408,7 +408,7 @@ def math_length_degree_distribution(W, D, nbins=10, nswap=1000, newB : (N, N) array-like binary rewired matrix newW: (N, N) array-like - weighted rewired matrix. Returns matrix of zeros of weighted=False. + weighted rewired matrix. Returns matrix of zeros if weighted=False. nr : int number of successful rewires From 723477823178dc168ee9a0f2953f58ca455dae27 Mon Sep 17 00:00:00 2001 From: Zhen-Qi Liu Date: Tue, 4 Oct 2022 10:33:53 -0400 Subject: [PATCH 3/4] Update networks.py Fix a typo --- netneurotools/networks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/netneurotools/networks.py b/netneurotools/networks.py index e602c3b..751c59a 100644 --- a/netneurotools/networks.py +++ b/netneurotools/networks.py @@ -378,7 +378,7 @@ def threshold_network(network, retain=10): return mst -def math_length_degree_distribution(W, D, nbins=10, nswap=1000, +def match_length_degree_distribution(W, D, nbins=10, nswap=1000, replacement=False, weighted=True, seed=None): """ From 79a58779014fa11293e9778f0296d79b8f35b313 Mon Sep 17 00:00:00 2001 From: Zhen-Qi Liu Date: Tue, 4 Oct 2022 11:03:21 -0400 Subject: [PATCH 4/4] [FIX] Fix style --- netneurotools/networks.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/netneurotools/networks.py b/netneurotools/networks.py index 751c59a..62c7dea 100644 --- a/netneurotools/networks.py +++ b/netneurotools/networks.py @@ -379,8 +379,8 @@ def threshold_network(network, retain=10): def match_length_degree_distribution(W, D, nbins=10, nswap=1000, - replacement=False, weighted=True, - seed=None): + replacement=False, weighted=True, + seed=None): """ Generates degree- and edge length-preserving surrogate connectomes. @@ -448,7 +448,7 @@ def match_length_degree_distribution(W, D, nbins=10, nswap=1000, nr = 0 newB = np.copy(B) - while((len(cn_x) >= 2) & (nr < nswap)): + while ((len(cn_x) >= 2) & (nr < nswap)): # choose randomly the edge to be rewired r = rs.randint(len(cn_x)) n_x, n_y = cn_x[r], cn_y[r] @@ -457,7 +457,7 @@ def match_length_degree_distribution(W, D, nbins=10, nswap=1000, # options to rewire with # connected nodes that doesn't involve (n_x, n_y) index = (cn_x != n_x) & (cn_y != n_y) & (cn_y != n_x) & (cn_x != n_y) - if(len(np.where(index)[0]) == 0): + if len(np.where(index)[0]) == 0: cn_x = np.delete(cn_x, r) cn_y = np.delete(cn_y, r) @@ -468,7 +468,7 @@ def match_length_degree_distribution(W, D, nbins=10, nswap=1000, # L(n_x,n_y) = L(n_x, ops1_x) & L(ops1_x,ops1_y) = L(n_y, ops1_y) index = (L[n_x, n_y] == L[n_x, ops1_x]) & ( L[ops1_x, ops1_y] == L[n_y, ops1_y]) - if(len(np.where(index)[0]) == 0): + if len(np.where(index)[0]) == 0: cn_x = np.delete(cn_x, r) cn_y = np.delete(cn_y, r) @@ -479,7 +479,7 @@ def match_length_degree_distribution(W, D, nbins=10, nswap=1000, & (newB[min(n_y, ops2_y[i])][max(n_y, ops2_y[i])] == 0) for i in range(len(ops2_x))] - if(len(np.where(index)[0]) == 0): + if (len(np.where(index)[0]) == 0): cn_x = np.delete(cn_x, r) cn_y = np.delete(cn_y, r) @@ -513,7 +513,7 @@ def match_length_degree_distribution(W, D, nbins=10, nswap=1000, cn_x = np.delete(cn_x, index) cn_y = np.delete(cn_y, index) - if(nr < nswap): + if nr < nswap: print(f"I didn't finish, out of rewirable edges: {len(cn_x)}") i, j = np.triu_indices(N, k=1) @@ -521,19 +521,19 @@ def match_length_degree_distribution(W, D, nbins=10, nswap=1000, newB[j, i] = newB[i, j] # check the number of edges is preserved - if(len(np.where(B != 0)[0]) != len(np.where(newB != 0)[0])): + if len(np.where(B != 0)[0]) != len(np.where(newB != 0)[0]): print( f"ERROR --- number of edges changed, \ B:{len(np.where(B!=0)[0])}, newB:{len(np.where(newB!=0)[0])}") # check that the degree of the nodes it's the same for i in range(N): - if(np.sum(B[i]) != np.sum(newB[i])): + if np.sum(B[i]) != np.sum(newB[i]): print( f"ERROR --- node {i} changed k by: \ {np.sum(B[i]) - np.sum(newB[i])}") newW = np.zeros((N, N)) - if(weighted): + if weighted: # Reassign the weights mask = np.triu(B != 0, k=1) inids = D[mask]