Generate all heritable regulatory architectures (HRAs) with n nodes.
This version emphasizes correctness while optimizing where possible.
11def generate_hras_accurate(n, verbose=True):
12 """
13 Generate all heritable regulatory architectures (HRAs) with n nodes.
14 This version emphasizes correctness while optimizing where possible.
15 """
16 start_time = time.time()
17 all_nodes = list(range(1, n + 1))
18 all_edges = list(itertools.permutations(all_nodes, r=2))
19
20 if verbose:
21 print(f"Finding weakly connected heritable topologies for n={n}...")
22
23
24
25 weakly_connected_heritable_non_iso_digraphs = []
26
27
28
29 min_edges = n
30
31
32
33
34
35
36 valid_edge_sets = []
37 for k in range(min_edges, len(all_edges) + 1):
38 if verbose and k == min_edges:
39 print(f"Generating edge sets with at least {min_edges} edges...")
40 edge_combinations = itertools.combinations(all_edges, r=k)
41 valid_edge_sets.extend(edge_combinations)
42
43 if verbose:
44 print(f"Checking {len(valid_edge_sets)} potential edge configurations...")
45
46
47 for idx, edge_set in enumerate(valid_edge_sets):
48 if verbose and idx % 10000 == 0 and idx > 0:
49 print(f"Processed {idx}/{len(valid_edge_sets)} edge sets...")
50
51 g = nx.DiGraph()
52 g.add_nodes_from(all_nodes)
53 g.add_edges_from(edge_set)
54
55
56 if not nx.is_weakly_connected(g):
57 continue
58
59 heritable = all(g.in_degree(i) > 0 for i in all_nodes)
60 if not heritable:
61 continue
62
63
64 found_isomorphic = False
65 for existing in weakly_connected_heritable_non_iso_digraphs:
66 if nx.is_isomorphic(g, existing):
67 found_isomorphic = True
68 break
69
70 if not found_isomorphic:
71 weakly_connected_heritable_non_iso_digraphs.append(g)
72
73 topology_time = time.time()
74 if verbose:
75 print(
76 f"Found {len(weakly_connected_heritable_non_iso_digraphs)} heritable topologies in {topology_time - start_time:.2f} seconds"
77 )
78
79
80 non_iso_regulatory_architectures = []
81
82 if verbose:
83 print("Generating regulatory architectures...")
84
85
86 for g_idx, g in enumerate(weakly_connected_heritable_non_iso_digraphs):
87 if verbose and g_idx % 10 == 0:
88 print(
89 f"Processing topology {g_idx+1}/{len(weakly_connected_heritable_non_iso_digraphs)}..."
90 )
91
92 edge_set = list(g.edges())
93 edge_count = len(edge_set)
94
95
96
97 for regulation_pattern in itertools.product(
98 ["gray", "black"], repeat=edge_count
99 ):
100
101 regulated_g = nx.DiGraph()
102 regulated_g.add_nodes_from(all_nodes)
103
104
105 regulation_to_edges = {}
106 for idx, (u, v) in enumerate(edge_set):
107 regulated_g.add_edge(u, v)
108 regulation_to_edges[(u, v)] = {"regulation": regulation_pattern[idx]}
109
110 nx.set_edge_attributes(regulated_g, regulation_to_edges)
111
112
113 found_isomorphic = False
114 for existing in non_iso_regulatory_architectures:
115 matcher = isomorphism.DiGraphMatcher(
116 regulated_g,
117 existing,
118 edge_match=lambda e1, e2: e1.get("regulation")
119 == e2.get("regulation"),
120 )
121 if matcher.is_isomorphic():
122 found_isomorphic = True
123 break
124
125 if not found_isomorphic:
126 non_iso_regulatory_architectures.append(regulated_g)
127
128 regulatory_time = time.time()
129 if verbose:
130 print(
131 f"Found {len(non_iso_regulatory_architectures)} non-isomorphic regulatory architectures in {regulatory_time - topology_time:.2f} seconds"
132 )
133
134
135
136 heritable_non_iso_regulatory_architectures = []
137
138 if verbose:
139 print("Filtering for heritable regulatory architectures...")
140
141 for regulated_g in non_iso_regulatory_architectures:
142 heritable = True
143 for node in all_nodes:
144 has_gray_incoming = False
145 for pred in regulated_g.predecessors(node):
146 if regulated_g[pred][node].get("regulation") == "gray":
147 has_gray_incoming = True
148 break
149
150 if not has_gray_incoming:
151 heritable = False
152 break
153
154 if heritable:
155 heritable_non_iso_regulatory_architectures.append(regulated_g)
156
157 end_time = time.time()
158 if verbose:
159 print(
160 f"Found {len(heritable_non_iso_regulatory_architectures)} heritable regulatory architectures in {end_time - start_time:.2f} seconds"
161 )
162 print(f"Total time: {end_time - start_time:.2f} seconds")
163
164 return {
165 "topologies": weakly_connected_heritable_non_iso_digraphs,
166 "regulatory": non_iso_regulatory_architectures,
167 "heritable": heritable_non_iso_regulatory_architectures,
168 }
169
170