threshold.py 30 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981
  1. """
  2. Threshold Graphs - Creation, manipulation and identification.
  3. """
  4. from math import sqrt
  5. import networkx as nx
  6. from networkx.utils import py_random_state
  7. __all__ = ["is_threshold_graph", "find_threshold_graph"]
  8. @nx._dispatchable
  9. def is_threshold_graph(G):
  10. """
  11. Returns `True` if `G` is a threshold graph.
  12. Parameters
  13. ----------
  14. G : NetworkX graph instance
  15. An instance of `Graph`, `DiGraph`, `MultiGraph` or `MultiDiGraph`
  16. Returns
  17. -------
  18. bool
  19. `True` if `G` is a threshold graph, `False` otherwise.
  20. Examples
  21. --------
  22. >>> from networkx.algorithms.threshold import is_threshold_graph
  23. >>> G = nx.path_graph(3)
  24. >>> is_threshold_graph(G)
  25. True
  26. >>> G = nx.barbell_graph(3, 3)
  27. >>> is_threshold_graph(G)
  28. False
  29. References
  30. ----------
  31. .. [1] Threshold graphs: https://en.wikipedia.org/wiki/Threshold_graph
  32. """
  33. return is_threshold_sequence([d for n, d in G.degree()])
  34. def is_threshold_sequence(degree_sequence):
  35. """
  36. Returns True if the sequence is a threshold degree sequence.
  37. Uses the property that a threshold graph must be constructed by
  38. adding either dominating or isolated nodes. Thus, it can be
  39. deconstructed iteratively by removing a node of degree zero or a
  40. node that connects to the remaining nodes. If this deconstruction
  41. fails then the sequence is not a threshold sequence.
  42. """
  43. ds = degree_sequence[:] # get a copy so we don't destroy original
  44. ds.sort()
  45. while ds:
  46. if ds[0] == 0: # if isolated node
  47. ds.pop(0) # remove it
  48. continue
  49. if ds[-1] != len(ds) - 1: # is the largest degree node dominating?
  50. return False # no, not a threshold degree sequence
  51. ds.pop() # yes, largest is the dominating node
  52. ds = [d - 1 for d in ds] # remove it and decrement all degrees
  53. return True
  54. def creation_sequence(degree_sequence, with_labels=False, compact=False):
  55. """
  56. Determines the creation sequence for the given threshold degree sequence.
  57. The creation sequence is a list of single characters 'd'
  58. or 'i': 'd' for dominating or 'i' for isolated vertices.
  59. Dominating vertices are connected to all vertices present when it
  60. is added. The first node added is by convention 'd'.
  61. This list can be converted to a string if desired using "".join(cs)
  62. If with_labels==True:
  63. Returns a list of 2-tuples containing the vertex number
  64. and a character 'd' or 'i' which describes the type of vertex.
  65. If compact==True:
  66. Returns the creation sequence in a compact form that is the number
  67. of 'i's and 'd's alternating.
  68. Examples:
  69. [1,2,2,3] represents d,i,i,d,d,i,i,i
  70. [3,1,2] represents d,d,d,i,d,d
  71. Notice that the first number is the first vertex to be used for
  72. construction and so is always 'd'.
  73. with_labels and compact cannot both be True.
  74. Returns None if the sequence is not a threshold sequence
  75. """
  76. if with_labels and compact:
  77. raise ValueError("compact sequences cannot be labeled")
  78. # make an indexed copy
  79. if isinstance(degree_sequence, dict): # labeled degree sequence
  80. ds = [[degree, label] for (label, degree) in degree_sequence.items()]
  81. else:
  82. ds = [[d, i] for i, d in enumerate(degree_sequence)]
  83. ds.sort()
  84. cs = [] # creation sequence
  85. while ds:
  86. if ds[0][0] == 0: # isolated node
  87. (d, v) = ds.pop(0)
  88. if len(ds) > 0: # make sure we start with a d
  89. cs.insert(0, (v, "i"))
  90. else:
  91. cs.insert(0, (v, "d"))
  92. continue
  93. if ds[-1][0] != len(ds) - 1: # Not dominating node
  94. return None # not a threshold degree sequence
  95. (d, v) = ds.pop()
  96. cs.insert(0, (v, "d"))
  97. ds = [[d[0] - 1, d[1]] for d in ds] # decrement due to removing node
  98. if with_labels:
  99. return cs
  100. if compact:
  101. return make_compact(cs)
  102. return [v[1] for v in cs] # not labeled
  103. def make_compact(creation_sequence):
  104. """
  105. Returns the creation sequence in a compact form
  106. that is the number of 'i's and 'd's alternating.
  107. Examples
  108. --------
  109. >>> from networkx.algorithms.threshold import make_compact
  110. >>> make_compact(["d", "i", "i", "d", "d", "i", "i", "i"])
  111. [1, 2, 2, 3]
  112. >>> make_compact(["d", "d", "d", "i", "d", "d"])
  113. [3, 1, 2]
  114. Notice that the first number is the first vertex
  115. to be used for construction and so is always 'd'.
  116. Labeled creation sequences lose their labels in the
  117. compact representation.
  118. >>> make_compact([3, 1, 2])
  119. [3, 1, 2]
  120. """
  121. first = creation_sequence[0]
  122. if isinstance(first, str): # creation sequence
  123. cs = creation_sequence[:]
  124. elif isinstance(first, tuple): # labeled creation sequence
  125. cs = [s[1] for s in creation_sequence]
  126. elif isinstance(first, int): # compact creation sequence
  127. return creation_sequence
  128. else:
  129. raise TypeError("Not a valid creation sequence type")
  130. ccs = []
  131. count = 1 # count the run lengths of d's or i's.
  132. for i in range(1, len(cs)):
  133. if cs[i] == cs[i - 1]:
  134. count += 1
  135. else:
  136. ccs.append(count)
  137. count = 1
  138. ccs.append(count) # don't forget the last one
  139. return ccs
  140. def uncompact(creation_sequence):
  141. """
  142. Converts a compact creation sequence for a threshold
  143. graph to a standard creation sequence (unlabeled).
  144. If the creation_sequence is already standard, return it.
  145. See creation_sequence.
  146. """
  147. first = creation_sequence[0]
  148. if isinstance(first, str): # creation sequence
  149. return creation_sequence
  150. elif isinstance(first, tuple): # labeled creation sequence
  151. return creation_sequence
  152. elif isinstance(first, int): # compact creation sequence
  153. ccscopy = creation_sequence[:]
  154. else:
  155. raise TypeError("Not a valid creation sequence type")
  156. cs = []
  157. while ccscopy:
  158. cs.extend(ccscopy.pop(0) * ["d"])
  159. if ccscopy:
  160. cs.extend(ccscopy.pop(0) * ["i"])
  161. return cs
  162. def creation_sequence_to_weights(creation_sequence):
  163. """
  164. Returns a list of node weights which create the threshold
  165. graph designated by the creation sequence. The weights
  166. are scaled so that the threshold is 1.0. The order of the
  167. nodes is the same as that in the creation sequence.
  168. """
  169. # Turn input sequence into a labeled creation sequence
  170. first = creation_sequence[0]
  171. if isinstance(first, str): # creation sequence
  172. if isinstance(creation_sequence, list):
  173. wseq = creation_sequence[:]
  174. else:
  175. wseq = list(creation_sequence) # string like 'ddidid'
  176. elif isinstance(first, tuple): # labeled creation sequence
  177. wseq = [v[1] for v in creation_sequence]
  178. elif isinstance(first, int): # compact creation sequence
  179. wseq = uncompact(creation_sequence)
  180. else:
  181. raise TypeError("Not a valid creation sequence type")
  182. # pass through twice--first backwards
  183. wseq.reverse()
  184. w = 0
  185. prev = "i"
  186. for j, s in enumerate(wseq):
  187. if s == "i":
  188. wseq[j] = w
  189. prev = s
  190. elif prev == "i":
  191. prev = s
  192. w += 1
  193. wseq.reverse() # now pass through forwards
  194. for j, s in enumerate(wseq):
  195. if s == "d":
  196. wseq[j] = w
  197. prev = s
  198. elif prev == "d":
  199. prev = s
  200. w += 1
  201. # Now scale weights
  202. if prev == "d":
  203. w += 1
  204. wscale = 1 / w
  205. return [ww * wscale for ww in wseq]
  206. # return wseq
  207. def weights_to_creation_sequence(
  208. weights, threshold=1, with_labels=False, compact=False
  209. ):
  210. """
  211. Returns a creation sequence for a threshold graph
  212. determined by the weights and threshold given as input.
  213. If the sum of two node weights is greater than the
  214. threshold value, an edge is created between these nodes.
  215. The creation sequence is a list of single characters 'd'
  216. or 'i': 'd' for dominating or 'i' for isolated vertices.
  217. Dominating vertices are connected to all vertices present
  218. when it is added. The first node added is by convention 'd'.
  219. If with_labels==True:
  220. Returns a list of 2-tuples containing the vertex number
  221. and a character 'd' or 'i' which describes the type of vertex.
  222. If compact==True:
  223. Returns the creation sequence in a compact form that is the number
  224. of 'i's and 'd's alternating.
  225. Examples:
  226. [1,2,2,3] represents d,i,i,d,d,i,i,i
  227. [3,1,2] represents d,d,d,i,d,d
  228. Notice that the first number is the first vertex to be used for
  229. construction and so is always 'd'.
  230. with_labels and compact cannot both be True.
  231. """
  232. if with_labels and compact:
  233. raise ValueError("compact sequences cannot be labeled")
  234. # make an indexed copy
  235. if isinstance(weights, dict): # labeled weights
  236. wseq = [[w, label] for (label, w) in weights.items()]
  237. else:
  238. wseq = [[w, i] for i, w in enumerate(weights)]
  239. wseq.sort()
  240. cs = [] # creation sequence
  241. cutoff = threshold - wseq[-1][0]
  242. while wseq:
  243. if wseq[0][0] < cutoff: # isolated node
  244. (w, label) = wseq.pop(0)
  245. cs.append((label, "i"))
  246. else:
  247. (w, label) = wseq.pop()
  248. cs.append((label, "d"))
  249. cutoff = threshold - wseq[-1][0]
  250. if len(wseq) == 1: # make sure we start with a d
  251. (w, label) = wseq.pop()
  252. cs.append((label, "d"))
  253. # put in correct order
  254. cs.reverse()
  255. if with_labels:
  256. return cs
  257. if compact:
  258. return make_compact(cs)
  259. return [v[1] for v in cs] # not labeled
  260. # Manipulating NetworkX.Graphs in context of threshold graphs
  261. @nx._dispatchable(graphs=None, returns_graph=True)
  262. def threshold_graph(creation_sequence, create_using=None):
  263. """
  264. Create a threshold graph from the creation sequence or compact
  265. creation_sequence.
  266. The input sequence can be a
  267. creation sequence (e.g. ['d','i','d','d','d','i'])
  268. labeled creation sequence (e.g. [(0,'d'),(2,'d'),(1,'i')])
  269. compact creation sequence (e.g. [2,1,1,2,0])
  270. Use cs=creation_sequence(degree_sequence,labeled=True)
  271. to convert a degree sequence to a creation sequence.
  272. Returns None if the sequence is not valid
  273. """
  274. # Turn input sequence into a labeled creation sequence
  275. first = creation_sequence[0]
  276. if isinstance(first, str): # creation sequence
  277. ci = list(enumerate(creation_sequence))
  278. elif isinstance(first, tuple): # labeled creation sequence
  279. ci = creation_sequence[:]
  280. elif isinstance(first, int): # compact creation sequence
  281. cs = uncompact(creation_sequence)
  282. ci = list(enumerate(cs))
  283. else:
  284. raise ValueError("not a valid creation sequence")
  285. G = nx.empty_graph(0, create_using)
  286. if G.is_directed():
  287. raise nx.NetworkXError("Directed Graph not supported")
  288. G.name = "Threshold Graph"
  289. # add nodes and edges
  290. # if type is 'i' just add nodea
  291. # if type is a d connect to everything previous
  292. while ci:
  293. (v, node_type) = ci.pop(0)
  294. if node_type == "d": # dominating type, connect to all existing nodes
  295. # We use `for u in list(G):` instead of
  296. # `for u in G:` because we edit the graph `G` in
  297. # the loop. Hence using an iterator will result in
  298. # `RuntimeError: dictionary changed size during iteration`
  299. for u in list(G):
  300. G.add_edge(v, u)
  301. G.add_node(v)
  302. return G
  303. @nx._dispatchable
  304. def find_alternating_4_cycle(G):
  305. """
  306. Returns False if there aren't any alternating 4 cycles.
  307. Otherwise returns the cycle as [a,b,c,d] where (a,b)
  308. and (c,d) are edges and (a,c) and (b,d) are not.
  309. """
  310. for u, v in G.edges():
  311. for w in G.nodes():
  312. if not G.has_edge(u, w) and u != w:
  313. for x in G.neighbors(w):
  314. if not G.has_edge(v, x) and v != x:
  315. return [u, v, w, x]
  316. return False
  317. @nx._dispatchable(returns_graph=True)
  318. def find_threshold_graph(G, create_using=None):
  319. """
  320. Returns a threshold subgraph that is close to largest in `G`.
  321. The threshold graph will contain the largest degree node in G.
  322. Parameters
  323. ----------
  324. G : NetworkX graph instance
  325. An instance of `Graph`, or `MultiDiGraph`
  326. create_using : NetworkX graph class or `None` (default), optional
  327. Type of graph to use when constructing the threshold graph.
  328. If `None`, infer the appropriate graph type from the input.
  329. Returns
  330. -------
  331. graph :
  332. A graph instance representing the threshold graph
  333. Examples
  334. --------
  335. >>> from networkx.algorithms.threshold import find_threshold_graph
  336. >>> G = nx.barbell_graph(3, 3)
  337. >>> T = find_threshold_graph(G)
  338. >>> T.nodes # may vary
  339. NodeView((7, 8, 5, 6))
  340. References
  341. ----------
  342. .. [1] Threshold graphs: https://en.wikipedia.org/wiki/Threshold_graph
  343. """
  344. return threshold_graph(find_creation_sequence(G), create_using)
  345. @nx._dispatchable
  346. def find_creation_sequence(G):
  347. """
  348. Find a threshold subgraph that is close to largest in G.
  349. Returns the labeled creation sequence of that threshold graph.
  350. """
  351. cs = []
  352. # get a local pointer to the working part of the graph
  353. H = G
  354. while H.order() > 0:
  355. # get new degree sequence on subgraph
  356. dsdict = dict(H.degree())
  357. ds = [(d, v) for v, d in dsdict.items()]
  358. ds.sort()
  359. # Update threshold graph nodes
  360. if ds[-1][0] == 0: # all are isolated
  361. cs.extend(zip(dsdict, ["i"] * (len(ds) - 1) + ["d"]))
  362. break # Done!
  363. # pull off isolated nodes
  364. while ds[0][0] == 0:
  365. (d, iso) = ds.pop(0)
  366. cs.append((iso, "i"))
  367. # find new biggest node
  368. (d, bigv) = ds.pop()
  369. # add edges of star to t_g
  370. cs.append((bigv, "d"))
  371. # form subgraph of neighbors of big node
  372. H = H.subgraph(H.neighbors(bigv))
  373. cs.reverse()
  374. return cs
  375. # Properties of Threshold Graphs
  376. def triangles(creation_sequence):
  377. """
  378. Compute number of triangles in the threshold graph with the
  379. given creation sequence.
  380. """
  381. # shortcut algorithm that doesn't require computing number
  382. # of triangles at each node.
  383. cs = creation_sequence # alias
  384. dr = cs.count("d") # number of d's in sequence
  385. ntri = dr * (dr - 1) * (dr - 2) / 6 # number of triangles in clique of nd d's
  386. # now add dr choose 2 triangles for every 'i' in sequence where
  387. # dr is the number of d's to the right of the current i
  388. for i, typ in enumerate(cs):
  389. if typ == "i":
  390. ntri += dr * (dr - 1) / 2
  391. else:
  392. dr -= 1
  393. return ntri
  394. def triangle_sequence(creation_sequence):
  395. """
  396. Return triangle sequence for the given threshold graph creation sequence.
  397. """
  398. cs = creation_sequence
  399. seq = []
  400. dr = cs.count("d") # number of d's to the right of the current pos
  401. dcur = (dr - 1) * (dr - 2) // 2 # number of triangles through a node of clique dr
  402. irun = 0 # number of i's in the last run
  403. drun = 0 # number of d's in the last run
  404. for i, sym in enumerate(cs):
  405. if sym == "d":
  406. drun += 1
  407. tri = dcur + (dr - 1) * irun # new triangles at this d
  408. else: # cs[i]="i":
  409. if prevsym == "d": # new string of i's
  410. dcur += (dr - 1) * irun # accumulate shared shortest paths
  411. irun = 0 # reset i run counter
  412. dr -= drun # reduce number of d's to right
  413. drun = 0 # reset d run counter
  414. irun += 1
  415. tri = dr * (dr - 1) // 2 # new triangles at this i
  416. seq.append(tri)
  417. prevsym = sym
  418. return seq
  419. def cluster_sequence(creation_sequence):
  420. """
  421. Return cluster sequence for the given threshold graph creation sequence.
  422. """
  423. triseq = triangle_sequence(creation_sequence)
  424. degseq = degree_sequence(creation_sequence)
  425. cseq = []
  426. for i, deg in enumerate(degseq):
  427. tri = triseq[i]
  428. if deg <= 1: # isolated vertex or single pair gets cc 0
  429. cseq.append(0)
  430. continue
  431. max_size = (deg * (deg - 1)) // 2
  432. cseq.append(tri / max_size)
  433. return cseq
  434. def degree_sequence(creation_sequence):
  435. """
  436. Return degree sequence for the threshold graph with the given
  437. creation sequence
  438. """
  439. cs = creation_sequence # alias
  440. seq = []
  441. rd = cs.count("d") # number of d to the right
  442. for i, sym in enumerate(cs):
  443. if sym == "d":
  444. rd -= 1
  445. seq.append(rd + i)
  446. else:
  447. seq.append(rd)
  448. return seq
  449. def density(creation_sequence):
  450. """
  451. Return the density of the graph with this creation_sequence.
  452. The density is the fraction of possible edges present.
  453. """
  454. N = len(creation_sequence)
  455. two_size = sum(degree_sequence(creation_sequence))
  456. two_possible = N * (N - 1)
  457. den = two_size / two_possible
  458. return den
  459. def degree_correlation(creation_sequence):
  460. """
  461. Return the degree-degree correlation over all edges.
  462. """
  463. cs = creation_sequence
  464. s1 = 0 # deg_i*deg_j
  465. s2 = 0 # deg_i^2+deg_j^2
  466. s3 = 0 # deg_i+deg_j
  467. m = 0 # number of edges
  468. rd = cs.count("d") # number of d nodes to the right
  469. rdi = [i for i, sym in enumerate(cs) if sym == "d"] # index of "d"s
  470. ds = degree_sequence(cs)
  471. for i, sym in enumerate(cs):
  472. if sym == "d":
  473. rdi.pop(0)
  474. degi = ds[i]
  475. for dj in rdi:
  476. degj = ds[dj]
  477. s1 += degj * degi
  478. s2 += degi**2 + degj**2
  479. s3 += degi + degj
  480. m += 1
  481. denom = 2 * m * s2 - s3 * s3
  482. numer = 4 * m * s1 - s3 * s3
  483. if denom == 0:
  484. if numer == 0:
  485. return 1
  486. raise ValueError(f"Zero Denominator but Numerator is {numer}")
  487. return numer / denom
  488. def shortest_path(creation_sequence, u, v):
  489. """
  490. Find the shortest path between u and v in a
  491. threshold graph G with the given creation_sequence.
  492. For an unlabeled creation_sequence, the vertices
  493. u and v must be integers in (0,len(sequence)) referring
  494. to the position of the desired vertices in the sequence.
  495. For a labeled creation_sequence, u and v are labels of vertices.
  496. Use cs=creation_sequence(degree_sequence,with_labels=True)
  497. to convert a degree sequence to a creation sequence.
  498. Returns a list of vertices from u to v.
  499. Example: if they are neighbors, it returns [u,v]
  500. """
  501. # Turn input sequence into a labeled creation sequence
  502. first = creation_sequence[0]
  503. if isinstance(first, str): # creation sequence
  504. cs = [(i, creation_sequence[i]) for i in range(len(creation_sequence))]
  505. elif isinstance(first, tuple): # labeled creation sequence
  506. cs = creation_sequence[:]
  507. elif isinstance(first, int): # compact creation sequence
  508. ci = uncompact(creation_sequence)
  509. cs = [(i, ci[i]) for i in range(len(ci))]
  510. else:
  511. raise TypeError("Not a valid creation sequence type")
  512. verts = [s[0] for s in cs]
  513. if v not in verts:
  514. raise ValueError(f"Vertex {v} not in graph from creation_sequence")
  515. if u not in verts:
  516. raise ValueError(f"Vertex {u} not in graph from creation_sequence")
  517. # Done checking
  518. if u == v:
  519. return [u]
  520. uindex = verts.index(u)
  521. vindex = verts.index(v)
  522. bigind = max(uindex, vindex)
  523. if cs[bigind][1] == "d":
  524. return [u, v]
  525. # must be that cs[bigind][1]=='i'
  526. cs = cs[bigind:]
  527. while cs:
  528. vert = cs.pop()
  529. if vert[1] == "d":
  530. return [u, vert[0], v]
  531. # All after u are type 'i' so no connection
  532. return -1
  533. def shortest_path_length(creation_sequence, i):
  534. """
  535. Return the shortest path length from indicated node to
  536. every other node for the threshold graph with the given
  537. creation sequence.
  538. Node is indicated by index i in creation_sequence unless
  539. creation_sequence is labeled in which case, i is taken to
  540. be the label of the node.
  541. Paths lengths in threshold graphs are at most 2.
  542. Length to unreachable nodes is set to -1.
  543. """
  544. # Turn input sequence into a labeled creation sequence
  545. first = creation_sequence[0]
  546. if isinstance(first, str): # creation sequence
  547. if isinstance(creation_sequence, list):
  548. cs = creation_sequence[:]
  549. else:
  550. cs = list(creation_sequence)
  551. elif isinstance(first, tuple): # labeled creation sequence
  552. cs = [v[1] for v in creation_sequence]
  553. i = [v[0] for v in creation_sequence].index(i)
  554. elif isinstance(first, int): # compact creation sequence
  555. cs = uncompact(creation_sequence)
  556. else:
  557. raise TypeError("Not a valid creation sequence type")
  558. # Compute
  559. N = len(cs)
  560. spl = [2] * N # length 2 to every node
  561. spl[i] = 0 # except self which is 0
  562. # 1 for all d's to the right
  563. for j in range(i + 1, N):
  564. if cs[j] == "d":
  565. spl[j] = 1
  566. if cs[i] == "d": # 1 for all nodes to the left
  567. for j in range(i):
  568. spl[j] = 1
  569. # and -1 for any trailing i to indicate unreachable
  570. for j in range(N - 1, 0, -1):
  571. if cs[j] == "d":
  572. break
  573. spl[j] = -1
  574. return spl
  575. def betweenness_sequence(creation_sequence, normalized=True):
  576. """
  577. Return betweenness for the threshold graph with the given creation
  578. sequence. The result is unscaled. To scale the values
  579. to the interval [0,1] divide by (n-1)*(n-2).
  580. """
  581. cs = creation_sequence
  582. seq = [] # betweenness
  583. lastchar = "d" # first node is always a 'd'
  584. dr = float(cs.count("d")) # number of d's to the right of current pos
  585. irun = 0 # number of i's in the last run
  586. drun = 0 # number of d's in the last run
  587. dlast = 0.0 # betweenness of last d
  588. for i, c in enumerate(cs):
  589. if c == "d": # cs[i]=="d":
  590. # betweenness = amt shared with earlier d's and i's
  591. # + new isolated nodes covered
  592. # + new paths to all previous nodes
  593. b = dlast + (irun - 1) * irun / dr + 2 * irun * (i - drun - irun) / dr
  594. drun += 1 # update counter
  595. else: # cs[i]="i":
  596. if lastchar == "d": # if this is a new run of i's
  597. dlast = b # accumulate betweenness
  598. dr -= drun # update number of d's to the right
  599. drun = 0 # reset d counter
  600. irun = 0 # reset i counter
  601. b = 0 # isolated nodes have zero betweenness
  602. irun += 1 # add another i to the run
  603. seq.append(float(b))
  604. lastchar = c
  605. # normalize by the number of possible shortest paths
  606. if normalized:
  607. order = len(cs)
  608. scale = 1.0 / ((order - 1) * (order - 2))
  609. seq = [s * scale for s in seq]
  610. return seq
  611. def eigenvectors(creation_sequence):
  612. """
  613. Return a 2-tuple of Laplacian eigenvalues and eigenvectors
  614. for the threshold network with creation_sequence.
  615. The first value is a list of eigenvalues.
  616. The second value is a list of eigenvectors.
  617. The lists are in the same order so corresponding eigenvectors
  618. and eigenvalues are in the same position in the two lists.
  619. Notice that the order of the eigenvalues returned by eigenvalues(cs)
  620. may not correspond to the order of these eigenvectors.
  621. """
  622. ccs = make_compact(creation_sequence)
  623. N = sum(ccs)
  624. vec = [0] * N
  625. val = vec[:]
  626. # get number of type d nodes to the right (all for first node)
  627. dr = sum(ccs[::2])
  628. nn = ccs[0]
  629. vec[0] = [1.0 / sqrt(N)] * N
  630. val[0] = 0
  631. e = dr
  632. dr -= nn
  633. type_d = True
  634. i = 1
  635. dd = 1
  636. while dd < nn:
  637. scale = 1.0 / sqrt(dd * dd + i)
  638. vec[i] = i * [-scale] + [dd * scale] + [0] * (N - i - 1)
  639. val[i] = e
  640. i += 1
  641. dd += 1
  642. if len(ccs) == 1:
  643. return (val, vec)
  644. for nn in ccs[1:]:
  645. scale = 1.0 / sqrt(nn * i * (i + nn))
  646. vec[i] = i * [-nn * scale] + nn * [i * scale] + [0] * (N - i - nn)
  647. # find eigenvalue
  648. type_d = not type_d
  649. if type_d:
  650. e = i + dr
  651. dr -= nn
  652. else:
  653. e = dr
  654. val[i] = e
  655. st = i
  656. i += 1
  657. dd = 1
  658. while dd < nn:
  659. scale = 1.0 / sqrt(i - st + dd * dd)
  660. vec[i] = [0] * st + (i - st) * [-scale] + [dd * scale] + [0] * (N - i - 1)
  661. val[i] = e
  662. i += 1
  663. dd += 1
  664. return (val, vec)
  665. def spectral_projection(u, eigenpairs):
  666. """
  667. Returns the coefficients of each eigenvector
  668. in a projection of the vector u onto the normalized
  669. eigenvectors which are contained in eigenpairs.
  670. eigenpairs should be a list of two objects. The
  671. first is a list of eigenvalues and the second a list
  672. of eigenvectors. The eigenvectors should be lists.
  673. There's not a lot of error checking on lengths of
  674. arrays, etc. so be careful.
  675. """
  676. coeff = []
  677. evect = eigenpairs[1]
  678. for ev in evect:
  679. c = sum(evv * uv for (evv, uv) in zip(ev, u))
  680. coeff.append(c)
  681. return coeff
  682. def eigenvalues(creation_sequence):
  683. """
  684. Return sequence of eigenvalues of the Laplacian of the threshold
  685. graph for the given creation_sequence.
  686. Based on the Ferrer's diagram method. The spectrum is integral
  687. and is the conjugate of the degree sequence.
  688. See::
  689. @Article{degree-merris-1994,
  690. author = {Russel Merris},
  691. title = {Degree maximal graphs are Laplacian integral},
  692. journal = {Linear Algebra Appl.},
  693. year = {1994},
  694. volume = {199},
  695. pages = {381--389},
  696. }
  697. """
  698. degseq = degree_sequence(creation_sequence)
  699. degseq.sort()
  700. eiglist = [] # zero is always one eigenvalue
  701. eig = 0
  702. row = len(degseq)
  703. bigdeg = degseq.pop()
  704. while row:
  705. if bigdeg < row:
  706. eiglist.append(eig)
  707. row -= 1
  708. else:
  709. eig += 1
  710. if degseq:
  711. bigdeg = degseq.pop()
  712. else:
  713. bigdeg = 0
  714. return eiglist
  715. # Threshold graph creation routines
  716. @py_random_state(2)
  717. def random_threshold_sequence(n, p, seed=None):
  718. """
  719. Create a random threshold sequence of size n.
  720. A creation sequence is built by randomly choosing d's with
  721. probability p and i's with probability 1-p.
  722. s=nx.random_threshold_sequence(10,0.5)
  723. returns a threshold sequence of length 10 with equal
  724. probably of an i or a d at each position.
  725. A "random" threshold graph can be built with
  726. G=nx.threshold_graph(s)
  727. seed : integer, random_state, or None (default)
  728. Indicator of random number generation state.
  729. See :ref:`Randomness<randomness>`.
  730. """
  731. if not (0 <= p <= 1):
  732. raise ValueError("p must be in [0,1]")
  733. cs = ["d"] # threshold sequences always start with a d
  734. for i in range(1, n):
  735. if seed.random() < p:
  736. cs.append("d")
  737. else:
  738. cs.append("i")
  739. return cs
  740. # maybe *_d_threshold_sequence routines should
  741. # be (or be called from) a single routine with a more descriptive name
  742. # and a keyword parameter?
  743. def right_d_threshold_sequence(n, m):
  744. """
  745. Returns a "right-dominated" threshold sequence with `n` vertices and `m` edges.
  746. Each vertex in the sequence is either dominant or isolated.
  747. In the "right-dominated" version, once the basic sequence is formed,
  748. isolated vertices may be flipped to dominant from the right in order
  749. to reach the target number of edges.
  750. Parameters
  751. ----------
  752. n : int
  753. Number of vertices.
  754. m : int
  755. Number of edges.
  756. Returns
  757. -------
  758. A list of 'd' (dominant) and 'i' (isolated) forming a right-dominated threshold sequence.
  759. Raises
  760. ------
  761. ValueError
  762. If `m` exceeds the maximum number of edges.
  763. Examples
  764. --------
  765. >>> from networkx.algorithms.threshold import right_d_threshold_sequence
  766. >>> right_d_threshold_sequence(5, 3)
  767. ['d', 'i', 'i', 'd', 'i']
  768. """
  769. cs = ["d"] + ["i"] * (n - 1) # create sequence with n insolated nodes
  770. # m <n : not enough edges, make disconnected
  771. if m < n:
  772. cs[m] = "d"
  773. return cs
  774. # too many edges
  775. if m > n * (n - 1) / 2:
  776. raise ValueError("Too many edges for this many nodes.")
  777. # connected case m >n-1
  778. ind = n - 1
  779. sum = n - 1
  780. while sum < m:
  781. cs[ind] = "d"
  782. ind -= 1
  783. sum += ind
  784. ind = m - (sum - ind)
  785. cs[ind] = "d"
  786. return cs
  787. def left_d_threshold_sequence(n, m):
  788. """
  789. Returns a "left-dominated" threshold sequence with `n` vertices and `m` edges.
  790. Each vertex in the sequence is either dominant or isolated.
  791. In the "left-dominated" version, once the basic sequence is formed,
  792. isolated vertices may be flipped to dominant from the left in order
  793. to reach the target number of edges.
  794. Parameters
  795. ----------
  796. n : int
  797. Number of vertices.
  798. m : int
  799. Number of edges.
  800. Returns
  801. -------
  802. A list of 'd' (dominant) and 'i' (isolated) forming a left-dominated threshold sequence.
  803. Raises
  804. ------
  805. ValueError
  806. If `m` exceeds the maximum number of edges.
  807. Examples
  808. --------
  809. For certain small cases, both left and right dominated versions produce
  810. the same sequence. However, for larger values of `m`, the difference in
  811. flipping order becomes evident. For instance, compare the sequences for
  812. ``n=6, m=8``:
  813. >>> from networkx.algorithms.threshold import left_d_threshold_sequence
  814. >>> seq = left_d_threshold_sequence(6, 8)
  815. >>> seq
  816. ['d', 'd', 'd', 'i', 'i', 'd']
  817. In contrast, the right-dominated version yields:
  818. >>> from networkx.algorithms.threshold import right_d_threshold_sequence
  819. >>> right_seq = right_d_threshold_sequence(6, 8)
  820. >>> right_seq
  821. ['d', 'i', 'i', 'd', 'i', 'd']
  822. """
  823. cs = ["d"] + ["i"] * (n - 1) # create sequence with n insolated nodes
  824. # m <n : not enough edges, make disconnected
  825. if m < n:
  826. cs[m] = "d"
  827. return cs
  828. # too many edges
  829. if m > n * (n - 1) / 2:
  830. raise ValueError("Too many edges for this many nodes.")
  831. # Connected case when M>N-1
  832. cs[n - 1] = "d"
  833. sum = n - 1
  834. ind = 1
  835. while sum < m:
  836. cs[ind] = "d"
  837. sum += ind
  838. ind += 1
  839. if sum > m: # be sure not to change the first vertex
  840. cs[sum - m] = "i"
  841. return cs