Coverage for src/driada/network/net_base.py: 53.30%
379 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-25 15:40 +0300
1from .graph_utils import *
2from .randomization import *
3from .drawing import *
4from .spectral import *
5from .quantum import *
7from scipy import linalg as la
8from scipy.sparse.linalg import eigs
9from sklearn.neighbors import NearestNeighbors
10from scipy.stats import entropy
12UNDIR_MATRIX_TYPES = ['adj', 'trans', 'lap', 'nlap', 'rwlap']
13DIR_MATRIX_TYPES = ['adj', 'lap_out', 'lap_in']
14MATRIX_TYPES = UNDIR_MATRIX_TYPES + DIR_MATRIX_TYPES
15SUPPORTED_GRAPH_TYPES = [nx.Graph, nx.DiGraph]
18def check_matrix_type(mode, is_directed):
19 if mode not in MATRIX_TYPES:
20 raise ValueError(f'Matrix type {mode} is not in allowed matrix types: {MATRIX_TYPES}')
22 if is_directed and mode not in DIR_MATRIX_TYPES:
23 raise ValueError(f'Matrix type {mode} is not allowed for directed networks.'
24 f'Supported options are: {DIR_MATRIX_TYPES}')
26 if not is_directed and mode not in UNDIR_MATRIX_TYPES:
27 raise ValueError(f'Matrix type {mode} is not allowed for undirected networks.'
28 f'Supported options are: {UNDIR_MATRIX_TYPES}')
31def check_adjacency(a):
32 if a.shape[0] != a.shape[1]:
33 raise Exception('Non-square adjacency matrix!')
36def check_directed(directed, real_world):
37 if real_world:
38 if int(directed) not in [0, 1]:
39 raise Exception('Fractional direction is not valid for a real network')
40 elif directed < 0 or directed > 1:
41 raise Exception('Wrong "directed" parameter value:', directed)
44def check_weights_and_directions(a, weighted, directed):
45 is_directed = not np.allclose(a.data, a.T.data)
46 is_weighted = not np.allclose(a.data, a.astype(bool).astype(int).data)
48 symm_text = 'asymmetric' if is_directed else 'symmetric'
49 if is_directed != bool(directed):
50 raise Exception(f'Error in network construction: "directed" set to {directed},'
51 f' but the adjacency matrix is {symm_text}')
53 w_text = 'weighted' if is_weighted else 'not weighted'
54 if is_weighted != bool(weighted):
55 raise Exception(f'Error in network construction: "weighted" set to {weighted},'
56 f' but the adjacency matrix {w_text}')
59def select_construction_pipeline(a, graph):
60 # TODO: add partial directions
61 if a is None:
62 if graph is None:
63 raise ValueError('Either "adj" or "graph" argument must be non-empty')
64 else:
65 if not np.any([isinstance(graph, gtype) for gtype in SUPPORTED_GRAPH_TYPES]):
66 raise TypeError(f'graph should have one of supported graph types: {SUPPORTED_GRAPH_TYPES}')
67 else:
68 pipeline = 'graph'
70 else:
71 if graph is None:
72 pipeline = 'adj'
73 else:
74 raise ValueError('Either "adj" or "graph" should be given, not both')
76 return pipeline
79class Network:
80 """
81 An object for network analysis with the focus on spectral graph theory
82 """
84 def __init__(self,
85 adj=None,
86 graph=None,
87 preprocessing='giant_cc',
88 name='',
89 pos=None,
90 verbose=False,
91 create_nx_graph=True,
92 node_attrs=None,
93 **network_args):
95 self.name = name
96 self.verbose = verbose
97 self.network_params = network_args
98 self.create_nx_graph = create_nx_graph
100 self.init_method = select_construction_pipeline(adj, graph)
102 self.directed = network_args.get('directed')
103 if self.directed is None:
104 if self.init_method == 'adj':
105 self.directed = not np.allclose(adj.data, adj.T.data)
106 elif self.init_method == 'graph':
107 self.directed = nx.is_directed(graph)
109 self.weighted = network_args.get('weighted')
110 if self.weighted is None:
111 if self.init_method == 'adj':
112 self.weighted = not np.allclose(adj.data, adj.astype(bool).astype(int).data)
113 elif self.init_method == 'graph':
114 self.weighted = nx.is_weighted(graph)
116 self.real_world = network_args.get('real_world')
117 if self.real_world is None:
118 self.real_world = True
120 check_directed(self.directed, self.real_world)
122 # set empty attributes for different matrix and data types
123 valid_mtypes = DIR_MATRIX_TYPES if self.directed else UNDIR_MATRIX_TYPES
124 for mt in valid_mtypes:
125 setattr(self, mt, None)
126 setattr(self, mt + '_spectrum', None)
127 setattr(self, mt + '_eigenvectors', None)
128 setattr(self, mt + '_zvalues', None)
129 setattr(self, mt + '_ipr', None)
131 # initialize adjacency matrix and (probably) associated graph
132 if self.init_method == 'adj':
133 # initialize Network object from sparse matrix
134 self._preprocess_adj_and_data(a=adj,
135 pos=pos,
136 node_attrs=node_attrs,
137 preprocessing=preprocessing,
138 create_graph=create_nx_graph)
140 if self.init_method == 'graph':
141 # initialize Network object from NetworkX graph or digraph
142 self._preprocess_graph_and_data(graph=graph,
143 pos=pos,
144 node_attrs=node_attrs,
145 preprocessing=preprocessing)
147 # each network object has out- and in-degree sequences from its initialization
148 self.get_node_degrees()
150 # TODO: remove lem_emb
151 self.lem_emb = None
153 def _preprocess_graph_and_data(self,
154 graph=None,
155 preprocessing=None,
156 pos=None,
157 node_attrs=None):
159 if preprocessing is None:
160 if self.verbose:
161 print('No preprocessing specified, this may lead to unexpected errors in graph connectivity!')
162 fgraph = remove_selfloops_from_graph(graph)
164 elif preprocessing == 'remove_isolates':
165 fgraph = remove_isolates_and_selfloops_from_graph(graph)
167 elif preprocessing == 'giant_cc':
168 g_ = remove_selfloops_from_graph(graph)
169 fgraph = get_giant_cc_from_graph(g_)
170 self.n_cc = 1
172 elif preprocessing == 'giant_scc':
173 g_ = remove_selfloops_from_graph(graph)
174 fgraph = get_giant_scc_from_graph(g_)
175 self.n_cc = 1
176 self.n_scc = 1
178 else:
179 raise ValueError('Wrong preprocessing type!')
181 lost_nodes = graph.number_of_nodes() - fgraph.number_of_nodes()
182 lost_edges = graph.number_of_edges() - fgraph.number_of_edges()
183 if lost_nodes + lost_edges != 0 and self.verbose:
184 print(f'{lost_nodes} nodes and {lost_edges} edges removed ')
186 # add node positions if provided
187 if pos is not None:
188 self.pos = {node: pos[node] for node in graph.nodes() if node in fgraph.nodes()}
189 else:
190 self.pos = None
192 if node_attrs is not None:
193 self.node_attrs = {node: node_attrs[node] for node in graph.nodes() if node in fgraph.nodes()}
194 else:
195 self.node_attrs = None
197 self.graph = fgraph
198 self.adj = nx.adjacency_matrix(fgraph)
199 self.n = nx.number_of_nodes(self.graph)
201 def _preprocess_adj_and_data(self,
202 a=None,
203 preprocessing=None,
204 pos=None,
205 node_attrs=None,
206 create_graph=True):
208 # if NetworkX graph should be created, we revert to graph-based initialization for simplicity
209 if create_graph:
210 gtype = nx.DiGraph if self.directed else nx.Graph
211 graph = nx.from_scipy_sparse_array(a, create_using=gtype)
212 self._preprocess_graph_and_data(graph=graph,
213 pos=pos,
214 node_attrs=node_attrs,
215 preprocessing=preprocessing)
217 return
219 if preprocessing is None:
220 if self.verbose:
221 print('No preprocessing specified, this may lead to unexpected errors in graph connectivity!')
222 fadj = remove_selfloops_from_adj(a)
223 nodes_range = range(fadj.shape[0])
224 node_mapping = dict(zip(nodes_range, nodes_range)) # no nodes have been deleted
226 elif preprocessing == 'remove_isolates':
227 a_ = remove_selfloops_from_adj(a)
228 fadj, node_mapping = remove_isolates_from_adj(a_)
230 elif preprocessing == 'giant_cc':
231 a_ = remove_selfloops_from_adj(a)
232 fadj, node_mapping = get_giant_cc_from_adj(a_)
233 self.n_cc = 1
235 elif preprocessing == 'giant_scc':
236 a_ = remove_selfloops_from_adj(a)
237 fadj, node_mapping = get_giant_scc_from_adj(a_)
238 self.n_cc = 1
239 self.n_scc = 1
241 else:
242 raise ValueError('Wrong preprocessing type!')
244 lost_nodes = a.shape[0] - fadj.shape[0]
245 lost_edges = a.nnz - fadj.nnz
246 if not self.directed:
247 lost_edges = lost_edges // 2
249 if lost_nodes + lost_edges != 0 and self.verbose:
250 print(f'{lost_nodes} nodes and {lost_edges} edges removed')
252 # add node positions if provided
253 if pos is not None:
254 self.pos = {node: pos[node] for node in range(a.shape[0]) if node in node_mapping}
255 else:
256 self.pos = None
258 if node_attrs is not None:
259 self.node_attrs = {node: node_attrs[node] for node in range(a.shape[0]) if node in node_mapping}
260 else:
261 self.node_attrs = None
263 self.graph = None
264 self.adj = fadj
265 self._init_to_final_node_mapping = node_mapping
266 self.n = self.adj.shape[0]
268 def is_connected(self):
269 ccs = list(get_ccs_from_adj(self.adj))
270 return len(ccs) == 1
272 def randomize(self, rmode='shuffle'):
273 # TODO: update routines
274 if rmode == 'graph_iom':
275 if self.directed:
276 g = nx.DiGraph(self.graph)
277 else:
278 g = nx.Graph(self.graph)
280 new_graph = random_rewiring_IOM_preserving(g, r=10)
281 rand_adj = nx.adjacency_matrix(new_graph)
283 elif rmode == 'adj_iom':
284 rand_adj = adj_random_rewiring_iom_preserving(self.adj,
285 is_weighted=self.weighted,
286 r=10)
288 elif rmode == 'complete':
289 rand_adj = random_rewiring_complete_graph(self.adj)
291 elif rmode == 'shuffle':
292 rand_adj = random_rewiring_dense_graph(self.adj)
294 else:
295 raise ValueError('Unknown randomization method')
297 rand_net = Network(sp.csr_matrix(rand_adj),
298 self.network_params,
299 name=self.name + f' {rmode} rand',
300 pos=self.pos,
301 real_world=False,
302 verbose=False)
304 return rand_net
306 def get_node_degrees(self):
307 # convert sparse matrix to 0-1 format and sum over specific axis
308 self.outdeg = np.array(self.adj.astype(bool).astype(int).sum(axis=0)).ravel()
309 self.indeg = np.array(self.adj.astype(bool).astype(int).sum(axis=1)).ravel()
310 self.deg = np.array((self.adj + self.adj.T).astype(bool).astype(int).sum(axis=1)).ravel()
312 min_out = min(self.outdeg)
313 min_in = min(self.indeg)
314 max_out = max(self.outdeg)
315 max_in = max(self.indeg)
317 if max_out != min_out:
318 self.scaled_outdeg = np.array([1.0 * (deg - min_out) / (max_out - min_out) for deg in self.outdeg])
319 else:
320 self.scaled_outdeg = np.ones(len(self.outdeg))
322 if min_in != max_in:
323 self.scaled_indeg = np.array([1.0 * (deg - min_in) / (max_in - min_in) for deg in self.indeg])
324 else:
325 self.scaled_outdeg = np.ones(len(self.indeg))
327 def get_degree_distr(self, mode='all'):
328 if mode == 'all':
329 deg = self.deg
330 elif mode == 'out':
331 deg = self.outdeg
332 elif mode == 'in':
333 deg = self.indeg
334 else:
335 raise ValueError('Wrong mode for degree distribution.')
337 hist, bins = np.histogram(deg, bins=max(deg) - min(deg), density=True)
338 return hist
340 def get_matrix(self, mode):
341 check_matrix_type(mode, self.directed)
342 matrix = getattr(self, mode)
343 if matrix is None:
344 if mode == 'lap' or mode == 'lap_out':
345 matrix = get_laplacian(self.adj)
346 elif mode == 'nlap':
347 matrix = get_norm_laplacian(self.adj)
348 elif mode == 'rwlap':
349 matrix = get_rw_laplacian(self.adj)
350 elif mode == 'trans':
351 matrix = get_trans_matrix(self.adj)
352 elif mode == 'adj':
353 matrix = self.adj
354 else:
355 raise Exception(f'Wrong matrix type: {mode}')
357 setattr(self, mode, matrix)
359 matrix = getattr(self, mode)
360 return matrix
362 def get_spectrum(self, mode):
363 check_matrix_type(mode, self.directed)
364 spectrum = getattr(self, mode + '_spectrum')
365 if spectrum is None:
366 self.diagonalize(mode=mode)
367 spectrum = getattr(self, mode + '_spectrum')
369 return spectrum
371 def get_eigenvectors(self, mode):
372 check_matrix_type(mode, self.directed)
373 eigenvectors = getattr(self, mode + '_eigenvectors')
374 if eigenvectors is None:
375 self.diagonalize(mode=mode)
376 eigenvectors = getattr(self, mode + '_eigenvectors')
378 return eigenvectors
380 def get_ipr(self, mode):
381 check_matrix_type(mode, self.directed)
382 ipr = getattr(self, mode + '_ipr')
383 if ipr is None:
384 self.calculate_ipr(mode=mode)
385 ipr = getattr(self, mode + '_ipr')
387 return ipr
389 def get_z_values(self, mode):
390 check_matrix_type(mode, self.directed)
391 zvals = getattr(self, mode + '_zvalues')
392 if zvals is None:
393 self.calculate_z_values(mode=mode)
394 zvals = getattr(self, mode + '_zvalues')
396 return zvals
398 def partial_diagonalize(self, spectrum_params):
399 '''
400 noise = self.spectrum_params['noise']
401 neigs = self.spectrum_params['neigs']
402 if noise == 0:
403 R = np.zeros((n, n))
404 else:
405 R = np.multiply(matrix.todense(), np.random.normal(loc=0.0, scale=noise, size=(n, n)))
406 '''
407 raise Exception('this method is under construction')
409 def diagonalize(self, mode='lap', verbose=None):
410 if verbose is None:
411 verbose = self.verbose
413 check_matrix_type(mode, self.directed)
414 if verbose:
415 print('Preparing for diagonalizing...')
417 outdeg = self.outdeg
418 indeg = self.indeg
419 deg = self.deg
421 A = self.adj.astype(float)
422 n = self.n
424 if n != np.count_nonzero(outdeg) and verbose:
425 print(n - np.count_nonzero(outdeg), 'nodes without out-edges')
426 if n != np.count_nonzero(indeg) and verbose:
427 print(n - np.count_nonzero(indeg), 'nodes without in-edges')
429 nz = np.count_nonzero(deg)
430 if nz != n and self.n_cc == 1:
431 print('Graph has', str(n - nz), 'isolated nodes!')
432 raise Exception('Graph is not connected!')
434 if not self.weighted and not self.directed and not np.allclose(outdeg, indeg):
435 raise Exception('out- and in- degrees do not coincide in boolean')
437 matrix = self.get_matrix(mode)
439 if verbose:
440 print('Performing diagonalization...')
442 matrix_is_symmetric = np.allclose(matrix.data, matrix.T.data)
443 if matrix_is_symmetric:
444 raw_eigs, right_eigvecs = la.eigh(matrix.todense())
445 else:
446 raw_eigs, right_eigvecs = la.eig(matrix.todense(), right=True)
448 raw_eigs = np.around(raw_eigs, decimals=12)
449 sorted_eigs = np.sort(raw_eigs)
451 if 'lap' in mode:
452 n_comp = len(raw_eigs[np.abs(raw_eigs) == 0])
453 if n_comp != 1 and not self.weighted and self.n_cc == 1:
454 print('eigenvalues:', sorted_eigs)
455 raise Exception('Graph has %d components!' % n_comp)
457 setattr(self, mode, matrix)
459 if np.allclose(np.imag(sorted_eigs), np.zeros(len(sorted_eigs)), atol=1e-12):
460 sorted_eigs = np.real(sorted_eigs)
461 else:
462 if not self.directed:
463 raise ValueError('Complex eigenvalues found in non-directed network!')
465 setattr(self, mode + '_spectrum', sorted_eigs)
467 sorted_eigenvectors = right_eigvecs[np.ix_(range(len(sorted_eigs)), np.argsort(raw_eigs))]
468 if np.allclose(np.imag(sorted_eigenvectors), np.zeros(sorted_eigenvectors.shape), atol=1e-8):
469 sorted_eigenvectors = np.real(sorted_eigenvectors)
470 else:
471 if not self.directed:
472 raise ValueError('Complex eigenvectors found in non-directed network!')
474 setattr(self, mode + '_eigenvectors', sorted_eigenvectors)
476 if verbose:
477 print('Diagonalizing finished')
479 # TODO: add Gromov hyperbolicity
481 def calculate_z_values(self, mode='lap'):
482 spectrum = self.get_spectrum(mode)
483 seigs = sorted(list(set(spectrum)), key=np.abs)
484 if len(seigs) != len(spectrum) and self.verbose:
485 print('WARNING:', len(spectrum) - len(seigs), 'repeated eigenvalues discarded')
487 if self.verbose:
488 print('Computing nearest neighbours...')
490 X = np.array([[np.real(x), np.imag(x)] for x in seigs])
491 nbrs = NearestNeighbors(n_neighbors=3, algorithm='ball_tree').fit(X)
492 distances, indices = nbrs.kneighbors(X)
493 nnbs = [seigs[x] for x in indices[:, 1]]
494 nnnbs = [seigs[x] for x in indices[:, 2]]
495 '''
496 nndist = np.array([(nnbs[i] - eigs[i]) for i in range(len(eigs))])
497 nnndist = np.array([(nnnbs[i] - eigs[i]) for i in range(len(eigs))])
498 '''
499 zlist = np.array([(nnbs[i] - seigs[i]) / (nnnbs[i] - seigs[i]) for i in range(len(seigs))])
500 zdict = dict(zip(seigs, zlist))
502 setattr(self, mode + '_zvalues', zdict)
504 def calculate_ipr(self, mode='adj'):
505 eigenvectors = self.get_eigenvectors(mode)
506 nvecs = eigenvectors.shape[1]
507 ipr = np.zeros(nvecs)
508 eig_entropy = np.zeros(nvecs)
510 for i in range(nvecs):
511 ipr[i] = sum([np.abs(v) ** 4 for v in eigenvectors[:, i]])
512 # entropy[i] = -np.log(ipr[i]) # erdos entropy (deprecated)
513 eig_entropy[i] = entropy(np.array([np.abs(v) ** 2 for v in eigenvectors[:, i]]))
515 setattr(self, mode + '_ipr', ipr)
516 # self.eigenvector_entropy = eig_entropy / np.log(self.n)
518 def _get_lap_spectrum(self, norm=False):
519 if not self.directed:
520 if norm:
521 spectrum = self.get_spectrum('nlap') # could be rwlap as well
522 else:
523 spectrum = self.get_spectrum('lap')
524 else:
525 if norm:
526 raise NotImplementedError('Normalized Laplacian not implemented for directed networks')
527 else:
528 spectrum = self.get_spectrum('lap_out')
530 return spectrum
532 def calculate_thermodynamic_entropy(self, tlist, verbose=False, norm=False):
533 spectrum = self._get_lap_spectrum(norm=norm)
534 res = [spectral_entropy(spectrum, t, verbose=verbose) for t in tlist]
535 return res
537 def calculate_free_entropy(self, tlist, norm=False):
538 spectrum = self._get_lap_spectrum(norm=norm)
539 res = [free_entropy(spectrum, t) for t in tlist]
540 return res
542 def calculate_q_entropy(self, q, tlist, norm=False):
543 spectrum = self._get_lap_spectrum(norm=norm)
544 res = [q_entropy(spectrum, t, q=q) for t in tlist]
545 return res
547 def calculate_estrada_communicability(self):
548 adj_spectrum = self.get_spectrum('adj')
549 self.estrada_communicability = sum([np.exp(e) for e in adj_spectrum])
550 return self.estrada_communicability
552 def get_estrada_bipartivity_index(self):
553 adj_spectrum = self.get_spectrum('adj')
554 esi1 = sum([np.exp(-e) for e in adj_spectrum])
555 esi2 = sum([np.exp(e) for e in adj_spectrum])
556 self.estrada_bipartivity = esi1 / esi2
557 return self.estrada_bipartivity
559 def localization_signatures(self, mode='lap'):
560 zvals = self.get_z_values(mode)
562 mean_cos_phi = np.mean(np.array([np.cos(np.angle(x)) for x in zvals]))
563 rvals = [1. / (np.abs(z)) ** 2 for z in zvals]
564 mean_inv_r_sq = np.mean(np.array(rvals))
566 if self.verbose:
567 print('mean cos phi complex:', mean_cos_phi)
568 print('mean 1/r^2 real:', mean_inv_r_sq)
570 return mean_inv_r_sq, mean_cos_phi
572 def construct_lem_embedding(self, dim):
573 if self.directed:
574 raise Exception('LEM embedding is not implemented for directed graphs')
576 A = self.adj
577 A = A.asfptype()
578 print('Performing spectral decomposition...')
579 K = A.shape[0]
580 NL = get_norm_laplacian(A)
581 DH = get_inv_sqrt_diag_matrix(A)
583 start_v = np.ones(K)
584 eigvals, eigvecs = eigs(NL, k=dim + 1, which='LR', v0=start_v, maxiter=K * 1000)
585 eigvals = np.asarray([np.round(np.real(x), 6) for x in eigvals])
587 if np.count_nonzero(eigvals == 1.0) > 1:
588 raise Exception('Error while LEM embedding construction: graph is not connected!')
589 else:
590 vecs = eigvecs.T[1:]
591 vec_norms = np.array([np.real(sum([x * x for x in v])) for v in vecs])
592 vecs = vecs / vec_norms[:, np.newaxis]
593 # explanation: https://jlmelville.github.io/smallvis/spectral.html
594 vecs = DH.dot(sp.csr_array(vecs, dtype=float))
596 self.lem_emb = vecs