amachine.am_create.am_structured_composition_optimal
1from collections import defaultdict 2import copy 3import random 4import warnings 5 6import numpy as np 7import networkx as nx 8 9from ..am_hmm import HMM 10from ..am_causal_state import CausalState 11from ..am_transition import Transition 12 13from ..am_random import exp_uniform_blend, resolve_rng 14from ..am_vocabulary import Vocabulary 15 16def optimal_node_permutation( 17 G: nx.DiGraph, 18 cohesion_matrix: np.ndarray, 19 instances_per_component: int, 20 max_search_time: int, 21) -> list[int] | None: 22 23 from ortools.sat.python import cp_model 24 25 nodes = list(G.nodes()) 26 N = len(nodes) 27 n_components = cohesion_matrix.shape[0] 28 node_to_idx = {node: i for i, node in enumerate(nodes)} 29 30 SCALE = 10_000 31 cohesion_int = np.round(cohesion_matrix * SCALE).astype(int) 32 flat = cohesion_int.flatten().tolist() 33 score_lo, score_hi = int(cohesion_int.min()), int(cohesion_int.max()) 34 35 model = cp_model.CpModel() 36 37 perm = [model.NewIntVar(0, N - 1, f"p{i}") for i in range(N)] 38 model.AddAllDifferent(perm) 39 40 edge_scores = [] 41 for u, v in G.edges(): 42 i, j = node_to_idx[u], node_to_idx[v] 43 44 # Map instance index → component index (equiv. of x // instances_per_component) 45 comp_i = model.NewIntVar(0, n_components - 1, f"ci_{i}_{j}") 46 comp_j = model.NewIntVar(0, n_components - 1, f"cj_{i}_{j}") 47 model.AddDivisionEquality(comp_i, perm[i], instances_per_component) 48 model.AddDivisionEquality(comp_j, perm[j], instances_per_component) 49 50 flat_idx = model.NewIntVar(0, n_components * n_components - 1, f"idx_{i}_{j}") 51 model.Add(flat_idx == n_components * comp_i + comp_j) 52 53 score = model.NewIntVar(score_lo, score_hi, f"s_{i}_{j}") 54 model.AddElement(flat_idx, flat, score) 55 edge_scores.append(score) 56 57 model.Maximize(cp_model.LinearExpr.Sum(edge_scores)) 58 59 solver = cp_model.CpSolver() 60 solver.parameters.num_search_workers = 3 61 solver.parameters.max_time_in_seconds = max_search_time 62 solver.parameters.linearization_level = 2 63 64 status = solver.Solve(model) 65 66 if status in (cp_model.OPTIMAL, cp_model.FEASIBLE): 67 return [solver.Value(perm[i]) for i in range(N)] 68 69 raise Exception("Solver couldn't find a feasible solution.") 70 71def structured_composition_optimal( 72 composition_id : str, 73 core_alphabet : list[str], 74 composition_rigidity : float, 75 instances_per_component : int, 76 component_groups : dict[ str, list[HMM] ], 77 component_residency_factor : float, 78 composition_connectivity_factor : float, 79 repetition_penalty : float = 0.75, 80 end_of_component_symbols : set[str] | None = None, 81 symbol_cohesion : np.ndarray | None = None, 82 max_search_time : int = 60, 83 random_seed : int | None = None ) -> HMM : 84 85 if end_of_component_symbols is not None and len( end_of_component_symbols ) == 0 : 86 raise ValueError( "end_of_component_symbols is empty" ) 87 88 if ( end_of_component_symbols is not None 89 and len( end_of_component_symbols ) == 1 90 and composition_connectivity_factor != 0.0 ) : 91 warnings.warn( "Only 1 EOC symbol, composition_connectivity_factor will be unused. " ) 92 93 if end_of_component_symbols is None and symbol_cohesion is None : 94 raise ValueError( "If using standard symbols instead of EOC symbols, must pass cohesion matrix." ) 95 96 if random_seed is not None: 97 py_rng = random.Random(random_seed) 98 np_rng = np.random.default_rng(random_seed) 99 else: 100 py_rng = random 101 np_rng = resolve_rng(None) 102 103 component_pool = [] 104 for groupd_id, components in component_groups.items() : 105 for c in components : 106 component_pool.append( ( c, groupd_id ) ) 107 108 n_components = len( component_pool ) 109 110 # This is the cohesion score for component i flowing into component j 111 component_cohesion = np.zeros( ( n_components, n_components ) ) 112 for i in range( n_components ) : 113 component_cohesion[ i, : ] = exp_uniform_blend( 114 n=n_components, 115 alpha=(1.0-composition_rigidity), 116 np_rng=np_rng ) 117 118 component_cohesion[ i, i ] *= ( 1.0 - repetition_penalty ) 119 120 total_subgraphs = n_components*instances_per_component 121 122 alphabet = core_alphabet.copy() 123 symbol_idx_map = { s : i for i, s in enumerate( alphabet ) } 124 125 symbol_set = set( alphabet ) 126 127 if end_of_component_symbols is not None : 128 # it is possible EOC symbols is already in the alphabet 129 # if not add them 130 for s in end_of_component_symbols : 131 if s in symbol_set : 132 continue 133 symbol_idx_map[ s ] = len( alphabet ) 134 alphabet.append( s ) 135 136 ############################################################## 137 # add the states and transitions 138 139 cstate_name_to_global_idxs = { 140 group_id : defaultdict(list) 141 for group_id in component_groups.keys() 142 } 143 144 cstate_instance_to_global_idx = [ {} for _ in range( total_subgraphs ) ] 145 146 global_idx_to_component_state = {} 147 global_idx_to_group = {} 148 g_sidx = 0 149 instance_idx = 0 150 151 #--------------------------------------------------------------------------- 152 153 states = [] 154 transitions = [] 155 156 for m_idx, ( m, groupd_id ) in enumerate( component_pool ) : 157 158 for _ in range( instances_per_component ) : 159 160 for s_idx, state in enumerate( m.states ) : 161 162 # should be assumed isoclasses are global properties 163 # in current usual case it will be the same as the group id 164 # but in general that might not be true 165 # what is true is iscoclasses are scoped within a group 166 # real problem is m_idx, g_idx may exist as a class name from a previous composition 167 # hierarchical composition requires hierarchical group and machine class names 168 # simple fix is use global group_id's and machine_id's 169 170 m_classes = { 171 f"g_{groupd_id}", 172 f"cm_{composition_id}_{m_idx}", 173 f"cmi_{composition_id}_{instance_idx}", 174 f"c_{composition_id}", 175 f"isoclass_{m.isoclass}" 176 } 177 178 # Note, later we will have to adjust/add the isomorphs since the names have changed. 179 180 states.append( CausalState( 181 name=f"{g_sidx}/{state.name}", 182 classes=( state.classes | m_classes ) 183 ) ) 184 185 global_idx_to_group[ g_sidx ] = groupd_id 186 global_idx_to_component_state[ g_sidx ] = state 187 cstate_name_to_global_idxs[ groupd_id ][ state.name ].append( g_sidx ) 188 189 cstate_instance_to_global_idx[ instance_idx ][ s_idx ] = g_sidx 190 191 g_sidx += 1 192 193 for tr in m.transitions : 194 origin = cstate_instance_to_global_idx[ instance_idx ][ tr.origin_state_idx ] 195 target = cstate_instance_to_global_idx[ instance_idx ][ tr.target_state_idx ] 196 transitions.append( 197 Transition( 198 origin_state_idx=origin, 199 target_state_idx=target, 200 prob=tr.prob, 201 symbol_idx=tr.symbol_idx, 202 pq=None 203 ) ) 204 205 instance_idx += 1 206 #--------------------------------------------------------------------------- 207 208 # Remap and add the isomorphs 209 for g_idx, s in enumerate( states ) : 210 groupd_id = global_idx_to_group[ g_idx ] 211 c_state = global_idx_to_component_state[ g_idx ] 212 for c_iso in c_state.isomorphs : 213 # Compnents are duplicated, which means isomorphs are too 214 g_isos = cstate_name_to_global_idxs[ groupd_id ][ c_iso ] 215 for g_iso in g_isos : 216 states[ g_idx ].add_isomorph( states[ g_iso ].name ) 217 218 ############################################################################################### 219 # Construct a cohesive component graph skeleton before adding the actual compononent connections 220 221 G = nx.DiGraph() 222 G.add_nodes_from( [ i for i in range( total_subgraphs ) ] ) 223 224 # minus self edges 225 all_edges = [ 226 (u, v) 227 for u in range(total_subgraphs) 228 for v in range(total_subgraphs) 229 if u != v 230 ] 231 232 py_rng.shuffle( all_edges ) 233 234 n_edges_from = defaultdict(int) 235 n_edges_to = defaultdict(int) 236 edges_added = 0 237 238 cycle_edges = [] 239 extra_edges = [] 240 241 max_edges = ( 1.0 + composition_connectivity_factor ) * total_subgraphs 242 243 # start with hamiltonian cycle 244 for u, v in all_edges: 245 246 if n_edges_from[u] >= 1 : 247 continue 248 249 if n_edges_to[v] >= 1: 250 continue 251 252 if nx.has_path( G, v, u ) and edges_added < total_subgraphs - 1: 253 continue 254 255 G.add_edge(u, v) 256 edges_added += 1 257 258 n_edges_from[ u ] += 1 259 n_edges_to[v] += 1 260 261 cycle_edges.append( ( u,v ) ) 262 263 if edges_added >= total_subgraphs : 264 if nx.is_strongly_connected(G): 265 break 266 267 # Get some candidate extra edges 268 for u, v in all_edges : 269 270 if G.has_edge( u , v ) : 271 continue 272 273 if end_of_component_symbols is not None and n_edges_from[u] >= len( end_of_component_symbols ) : 274 continue 275 276 edges_added += 1 277 278 n_edges_from[ u ] += 1 279 n_edges_to[ v ] += 1 280 281 extra_edges.append( (u,v) ) 282 283 if edges_added >= max_edges : 284 break 285 286 ############################################################## 287 # calculate best entry and exit nodes for each subgraph 288 289 # lower is better 290 def score_as_from( m, s_idx ) : 291 n_outgoing = 0 292 has_self_loop = False 293 for tr in m.transitions : 294 if tr.origin_state_idx == s_idx : 295 n_outgoing += 1 296 if tr.origin_state_idx == s_idx and tr.target_state_idx == s_idx : 297 has_self_loop = True 298 return ( n_outgoing, not has_self_loop ) 299 300 # higher is better 301 def score_as_to( m, s_idx, path_length ) : 302 303 n_incoming = 0 304 has_self_loop = False 305 for tr in m.transitions : 306 if tr.target_state_idx == s_idx : 307 n_incoming += 1 308 if tr.origin_state_idx == s_idx and tr.target_state_idx == s_idx : 309 has_self_loop = True 310 311 return ( path_length, 1.0/n_incoming, has_self_loop ) 312 313 best_entry_and_exits = [] 314 315 for cm, _ in component_pool: 316 317 n_states = len(cm.states) 318 cm_dg = cm.as_digraph() 319 320 best_from = min(range(n_states), key=lambda s: score_as_from(cm, s)) 321 322 lengths_to_from = nx.single_source_shortest_path_length( 323 cm_dg.reverse(), best_from) 324 325 best_to = max( 326 range(n_states), 327 key=lambda s: score_as_to(cm, s, lengths_to_from.get(s, 0))) 328 329 best_entry_and_exits.append( ( best_to, best_from ) ) 330 331 ############################################################## 332 333 node_permutation = optimal_node_permutation( 334 G=G, 335 cohesion_matrix=component_cohesion, 336 instances_per_component=instances_per_component, 337 max_search_time=60 338 ) 339 340 nodes = list( G.nodes() ) 341 mapping = { i : node_permutation[i] for i in range(total_subgraphs) } 342 cycle_edges = [ ( mapping[u], mapping[v] ) for u, v in cycle_edges ] 343 extra_edges = [ ( mapping[u], mapping[v] ) for u, v in extra_edges ] 344 345 ############################################################## 346 # Connect the components 347 348 # In theory could score to and from nodes also based on what kind of cohesion scores are possible 349 # Given what free symbols they would have, but this is a bit complicated, considering adding it later. 350 351 extra_edges_added = 0 352 for edge_idx, ( u, v ) in enumerate( cycle_edges + extra_edges ) : 353 354 is_cycle_edge = edge_idx < len( cycle_edges ) 355 356 m_u_idx = u // instances_per_component 357 m_v_idx = v // instances_per_component 358 359 m_u, _ = component_pool[ m_u_idx ] 360 m_v, _ = component_pool[ m_v_idx ] 361 362 best_from = best_entry_and_exits[ m_u_idx ][ 1 ] 363 best_to = best_entry_and_exits[ m_v_idx ][ 0 ] 364 365 origin = cstate_instance_to_global_idx[ u ][ best_from ] 366 target = cstate_instance_to_global_idx[ v ][ best_to ] 367 368 # if end_of_component_symbols is None we use the whole alphabet, otherwise eoc only 369 free_eoc_symbols = end_of_component_symbols.copy() if end_of_component_symbols is not None else set() 370 free_normal_symbols = set( alphabet ) if end_of_component_symbols is None else set() 371 372 # reweight existing transitions and figure out which eoc symbols are available 373 for tr in transitions : 374 if tr.origin_state_idx == origin : 375 tr.prob = tr.prob * component_residency_factor 376 symbol = alphabet[ tr.symbol_idx ] 377 if symbol in free_eoc_symbols : 378 free_eoc_symbols.remove( symbol ) 379 if symbol in free_normal_symbols : 380 free_normal_symbols.remove( symbol ) 381 382 # Shouldn't happen on a cycle edge 383 if len( free_eoc_symbols ) == 0 and end_of_component_symbols is not None : 384 if is_cycle_edge : 385 raise ValueError( "Ran out of eoc symbols attempting to connect components." ) 386 else : 387 continue 388 389 # This could happen on a cycle edge 390 if len( free_normal_symbols ) == 0 and end_of_component_symbols is None : 391 if is_cycle_edge : 392 raise ValueError( "Ran out of symbols attempting to connect components." ) 393 else : 394 continue 395 396 if not is_cycle_edge : 397 extra_edges_added += 1 398 399 if len( free_eoc_symbols ) > 0 : 400 symbol = py_rng.choice( sorted( free_eoc_symbols ) ) 401 402 elif len( free_normal_symbols ) > 0 : 403 avail = list( free_normal_symbols ) 404 cohesion_scores = [ 0 for _ in avail ] 405 for i, s in enumerate( avail ) : 406 for tr in transitions : 407 if tr.target_state_idx == origin : 408 cohesion_scores[ i ] += symbol_cohesion[ tr.symbol_idx, symbol_idx_map[ s ] ] 409 elif tr.origin_state_idx == target : 410 cohesion_scores[ i ] += symbol_cohesion[ symbol_idx_map[ s ], tr.symbol_idx ] 411 best_i = np.argmax( np.array( cohesion_scores ) ) 412 symbol = avail[ best_i ] 413 414 transitions.append( 415 Transition( 416 origin_state_idx=origin, 417 target_state_idx=target, 418 prob=( 1.0 - component_residency_factor ), 419 symbol_idx=symbol_idx_map[ symbol ], 420 pq=None 421 ) ) 422 423 print( f"Was able to add {extra_edges_added}/{len(extra_edges)} extra cross-component transitions beyond a Hamiltonian cycle" ) 424 425 426 return HMM( 427 states=states, 428 transitions=transitions, 429 start_state=0, 430 alphabet=alphabet 431 )
def
optimal_node_permutation( G: networkx.classes.digraph.DiGraph, cohesion_matrix: numpy.ndarray, instances_per_component: int, max_search_time: int) -> list[int] | None:
17def optimal_node_permutation( 18 G: nx.DiGraph, 19 cohesion_matrix: np.ndarray, 20 instances_per_component: int, 21 max_search_time: int, 22) -> list[int] | None: 23 24 from ortools.sat.python import cp_model 25 26 nodes = list(G.nodes()) 27 N = len(nodes) 28 n_components = cohesion_matrix.shape[0] 29 node_to_idx = {node: i for i, node in enumerate(nodes)} 30 31 SCALE = 10_000 32 cohesion_int = np.round(cohesion_matrix * SCALE).astype(int) 33 flat = cohesion_int.flatten().tolist() 34 score_lo, score_hi = int(cohesion_int.min()), int(cohesion_int.max()) 35 36 model = cp_model.CpModel() 37 38 perm = [model.NewIntVar(0, N - 1, f"p{i}") for i in range(N)] 39 model.AddAllDifferent(perm) 40 41 edge_scores = [] 42 for u, v in G.edges(): 43 i, j = node_to_idx[u], node_to_idx[v] 44 45 # Map instance index → component index (equiv. of x // instances_per_component) 46 comp_i = model.NewIntVar(0, n_components - 1, f"ci_{i}_{j}") 47 comp_j = model.NewIntVar(0, n_components - 1, f"cj_{i}_{j}") 48 model.AddDivisionEquality(comp_i, perm[i], instances_per_component) 49 model.AddDivisionEquality(comp_j, perm[j], instances_per_component) 50 51 flat_idx = model.NewIntVar(0, n_components * n_components - 1, f"idx_{i}_{j}") 52 model.Add(flat_idx == n_components * comp_i + comp_j) 53 54 score = model.NewIntVar(score_lo, score_hi, f"s_{i}_{j}") 55 model.AddElement(flat_idx, flat, score) 56 edge_scores.append(score) 57 58 model.Maximize(cp_model.LinearExpr.Sum(edge_scores)) 59 60 solver = cp_model.CpSolver() 61 solver.parameters.num_search_workers = 3 62 solver.parameters.max_time_in_seconds = max_search_time 63 solver.parameters.linearization_level = 2 64 65 status = solver.Solve(model) 66 67 if status in (cp_model.OPTIMAL, cp_model.FEASIBLE): 68 return [solver.Value(perm[i]) for i in range(N)] 69 70 raise Exception("Solver couldn't find a feasible solution.")
def
structured_composition_optimal( composition_id: str, core_alphabet: list[str], composition_rigidity: float, instances_per_component: int, component_groups: dict[str, list[amachine.am_hmm.HMM]], component_residency_factor: float, composition_connectivity_factor: float, repetition_penalty: float = 0.75, end_of_component_symbols: set[str] | None = None, symbol_cohesion: numpy.ndarray | None = None, max_search_time: int = 60, random_seed: int | None = None) -> amachine.am_hmm.HMM:
72def structured_composition_optimal( 73 composition_id : str, 74 core_alphabet : list[str], 75 composition_rigidity : float, 76 instances_per_component : int, 77 component_groups : dict[ str, list[HMM] ], 78 component_residency_factor : float, 79 composition_connectivity_factor : float, 80 repetition_penalty : float = 0.75, 81 end_of_component_symbols : set[str] | None = None, 82 symbol_cohesion : np.ndarray | None = None, 83 max_search_time : int = 60, 84 random_seed : int | None = None ) -> HMM : 85 86 if end_of_component_symbols is not None and len( end_of_component_symbols ) == 0 : 87 raise ValueError( "end_of_component_symbols is empty" ) 88 89 if ( end_of_component_symbols is not None 90 and len( end_of_component_symbols ) == 1 91 and composition_connectivity_factor != 0.0 ) : 92 warnings.warn( "Only 1 EOC symbol, composition_connectivity_factor will be unused. " ) 93 94 if end_of_component_symbols is None and symbol_cohesion is None : 95 raise ValueError( "If using standard symbols instead of EOC symbols, must pass cohesion matrix." ) 96 97 if random_seed is not None: 98 py_rng = random.Random(random_seed) 99 np_rng = np.random.default_rng(random_seed) 100 else: 101 py_rng = random 102 np_rng = resolve_rng(None) 103 104 component_pool = [] 105 for groupd_id, components in component_groups.items() : 106 for c in components : 107 component_pool.append( ( c, groupd_id ) ) 108 109 n_components = len( component_pool ) 110 111 # This is the cohesion score for component i flowing into component j 112 component_cohesion = np.zeros( ( n_components, n_components ) ) 113 for i in range( n_components ) : 114 component_cohesion[ i, : ] = exp_uniform_blend( 115 n=n_components, 116 alpha=(1.0-composition_rigidity), 117 np_rng=np_rng ) 118 119 component_cohesion[ i, i ] *= ( 1.0 - repetition_penalty ) 120 121 total_subgraphs = n_components*instances_per_component 122 123 alphabet = core_alphabet.copy() 124 symbol_idx_map = { s : i for i, s in enumerate( alphabet ) } 125 126 symbol_set = set( alphabet ) 127 128 if end_of_component_symbols is not None : 129 # it is possible EOC symbols is already in the alphabet 130 # if not add them 131 for s in end_of_component_symbols : 132 if s in symbol_set : 133 continue 134 symbol_idx_map[ s ] = len( alphabet ) 135 alphabet.append( s ) 136 137 ############################################################## 138 # add the states and transitions 139 140 cstate_name_to_global_idxs = { 141 group_id : defaultdict(list) 142 for group_id in component_groups.keys() 143 } 144 145 cstate_instance_to_global_idx = [ {} for _ in range( total_subgraphs ) ] 146 147 global_idx_to_component_state = {} 148 global_idx_to_group = {} 149 g_sidx = 0 150 instance_idx = 0 151 152 #--------------------------------------------------------------------------- 153 154 states = [] 155 transitions = [] 156 157 for m_idx, ( m, groupd_id ) in enumerate( component_pool ) : 158 159 for _ in range( instances_per_component ) : 160 161 for s_idx, state in enumerate( m.states ) : 162 163 # should be assumed isoclasses are global properties 164 # in current usual case it will be the same as the group id 165 # but in general that might not be true 166 # what is true is iscoclasses are scoped within a group 167 # real problem is m_idx, g_idx may exist as a class name from a previous composition 168 # hierarchical composition requires hierarchical group and machine class names 169 # simple fix is use global group_id's and machine_id's 170 171 m_classes = { 172 f"g_{groupd_id}", 173 f"cm_{composition_id}_{m_idx}", 174 f"cmi_{composition_id}_{instance_idx}", 175 f"c_{composition_id}", 176 f"isoclass_{m.isoclass}" 177 } 178 179 # Note, later we will have to adjust/add the isomorphs since the names have changed. 180 181 states.append( CausalState( 182 name=f"{g_sidx}/{state.name}", 183 classes=( state.classes | m_classes ) 184 ) ) 185 186 global_idx_to_group[ g_sidx ] = groupd_id 187 global_idx_to_component_state[ g_sidx ] = state 188 cstate_name_to_global_idxs[ groupd_id ][ state.name ].append( g_sidx ) 189 190 cstate_instance_to_global_idx[ instance_idx ][ s_idx ] = g_sidx 191 192 g_sidx += 1 193 194 for tr in m.transitions : 195 origin = cstate_instance_to_global_idx[ instance_idx ][ tr.origin_state_idx ] 196 target = cstate_instance_to_global_idx[ instance_idx ][ tr.target_state_idx ] 197 transitions.append( 198 Transition( 199 origin_state_idx=origin, 200 target_state_idx=target, 201 prob=tr.prob, 202 symbol_idx=tr.symbol_idx, 203 pq=None 204 ) ) 205 206 instance_idx += 1 207 #--------------------------------------------------------------------------- 208 209 # Remap and add the isomorphs 210 for g_idx, s in enumerate( states ) : 211 groupd_id = global_idx_to_group[ g_idx ] 212 c_state = global_idx_to_component_state[ g_idx ] 213 for c_iso in c_state.isomorphs : 214 # Compnents are duplicated, which means isomorphs are too 215 g_isos = cstate_name_to_global_idxs[ groupd_id ][ c_iso ] 216 for g_iso in g_isos : 217 states[ g_idx ].add_isomorph( states[ g_iso ].name ) 218 219 ############################################################################################### 220 # Construct a cohesive component graph skeleton before adding the actual compononent connections 221 222 G = nx.DiGraph() 223 G.add_nodes_from( [ i for i in range( total_subgraphs ) ] ) 224 225 # minus self edges 226 all_edges = [ 227 (u, v) 228 for u in range(total_subgraphs) 229 for v in range(total_subgraphs) 230 if u != v 231 ] 232 233 py_rng.shuffle( all_edges ) 234 235 n_edges_from = defaultdict(int) 236 n_edges_to = defaultdict(int) 237 edges_added = 0 238 239 cycle_edges = [] 240 extra_edges = [] 241 242 max_edges = ( 1.0 + composition_connectivity_factor ) * total_subgraphs 243 244 # start with hamiltonian cycle 245 for u, v in all_edges: 246 247 if n_edges_from[u] >= 1 : 248 continue 249 250 if n_edges_to[v] >= 1: 251 continue 252 253 if nx.has_path( G, v, u ) and edges_added < total_subgraphs - 1: 254 continue 255 256 G.add_edge(u, v) 257 edges_added += 1 258 259 n_edges_from[ u ] += 1 260 n_edges_to[v] += 1 261 262 cycle_edges.append( ( u,v ) ) 263 264 if edges_added >= total_subgraphs : 265 if nx.is_strongly_connected(G): 266 break 267 268 # Get some candidate extra edges 269 for u, v in all_edges : 270 271 if G.has_edge( u , v ) : 272 continue 273 274 if end_of_component_symbols is not None and n_edges_from[u] >= len( end_of_component_symbols ) : 275 continue 276 277 edges_added += 1 278 279 n_edges_from[ u ] += 1 280 n_edges_to[ v ] += 1 281 282 extra_edges.append( (u,v) ) 283 284 if edges_added >= max_edges : 285 break 286 287 ############################################################## 288 # calculate best entry and exit nodes for each subgraph 289 290 # lower is better 291 def score_as_from( m, s_idx ) : 292 n_outgoing = 0 293 has_self_loop = False 294 for tr in m.transitions : 295 if tr.origin_state_idx == s_idx : 296 n_outgoing += 1 297 if tr.origin_state_idx == s_idx and tr.target_state_idx == s_idx : 298 has_self_loop = True 299 return ( n_outgoing, not has_self_loop ) 300 301 # higher is better 302 def score_as_to( m, s_idx, path_length ) : 303 304 n_incoming = 0 305 has_self_loop = False 306 for tr in m.transitions : 307 if tr.target_state_idx == s_idx : 308 n_incoming += 1 309 if tr.origin_state_idx == s_idx and tr.target_state_idx == s_idx : 310 has_self_loop = True 311 312 return ( path_length, 1.0/n_incoming, has_self_loop ) 313 314 best_entry_and_exits = [] 315 316 for cm, _ in component_pool: 317 318 n_states = len(cm.states) 319 cm_dg = cm.as_digraph() 320 321 best_from = min(range(n_states), key=lambda s: score_as_from(cm, s)) 322 323 lengths_to_from = nx.single_source_shortest_path_length( 324 cm_dg.reverse(), best_from) 325 326 best_to = max( 327 range(n_states), 328 key=lambda s: score_as_to(cm, s, lengths_to_from.get(s, 0))) 329 330 best_entry_and_exits.append( ( best_to, best_from ) ) 331 332 ############################################################## 333 334 node_permutation = optimal_node_permutation( 335 G=G, 336 cohesion_matrix=component_cohesion, 337 instances_per_component=instances_per_component, 338 max_search_time=60 339 ) 340 341 nodes = list( G.nodes() ) 342 mapping = { i : node_permutation[i] for i in range(total_subgraphs) } 343 cycle_edges = [ ( mapping[u], mapping[v] ) for u, v in cycle_edges ] 344 extra_edges = [ ( mapping[u], mapping[v] ) for u, v in extra_edges ] 345 346 ############################################################## 347 # Connect the components 348 349 # In theory could score to and from nodes also based on what kind of cohesion scores are possible 350 # Given what free symbols they would have, but this is a bit complicated, considering adding it later. 351 352 extra_edges_added = 0 353 for edge_idx, ( u, v ) in enumerate( cycle_edges + extra_edges ) : 354 355 is_cycle_edge = edge_idx < len( cycle_edges ) 356 357 m_u_idx = u // instances_per_component 358 m_v_idx = v // instances_per_component 359 360 m_u, _ = component_pool[ m_u_idx ] 361 m_v, _ = component_pool[ m_v_idx ] 362 363 best_from = best_entry_and_exits[ m_u_idx ][ 1 ] 364 best_to = best_entry_and_exits[ m_v_idx ][ 0 ] 365 366 origin = cstate_instance_to_global_idx[ u ][ best_from ] 367 target = cstate_instance_to_global_idx[ v ][ best_to ] 368 369 # if end_of_component_symbols is None we use the whole alphabet, otherwise eoc only 370 free_eoc_symbols = end_of_component_symbols.copy() if end_of_component_symbols is not None else set() 371 free_normal_symbols = set( alphabet ) if end_of_component_symbols is None else set() 372 373 # reweight existing transitions and figure out which eoc symbols are available 374 for tr in transitions : 375 if tr.origin_state_idx == origin : 376 tr.prob = tr.prob * component_residency_factor 377 symbol = alphabet[ tr.symbol_idx ] 378 if symbol in free_eoc_symbols : 379 free_eoc_symbols.remove( symbol ) 380 if symbol in free_normal_symbols : 381 free_normal_symbols.remove( symbol ) 382 383 # Shouldn't happen on a cycle edge 384 if len( free_eoc_symbols ) == 0 and end_of_component_symbols is not None : 385 if is_cycle_edge : 386 raise ValueError( "Ran out of eoc symbols attempting to connect components." ) 387 else : 388 continue 389 390 # This could happen on a cycle edge 391 if len( free_normal_symbols ) == 0 and end_of_component_symbols is None : 392 if is_cycle_edge : 393 raise ValueError( "Ran out of symbols attempting to connect components." ) 394 else : 395 continue 396 397 if not is_cycle_edge : 398 extra_edges_added += 1 399 400 if len( free_eoc_symbols ) > 0 : 401 symbol = py_rng.choice( sorted( free_eoc_symbols ) ) 402 403 elif len( free_normal_symbols ) > 0 : 404 avail = list( free_normal_symbols ) 405 cohesion_scores = [ 0 for _ in avail ] 406 for i, s in enumerate( avail ) : 407 for tr in transitions : 408 if tr.target_state_idx == origin : 409 cohesion_scores[ i ] += symbol_cohesion[ tr.symbol_idx, symbol_idx_map[ s ] ] 410 elif tr.origin_state_idx == target : 411 cohesion_scores[ i ] += symbol_cohesion[ symbol_idx_map[ s ], tr.symbol_idx ] 412 best_i = np.argmax( np.array( cohesion_scores ) ) 413 symbol = avail[ best_i ] 414 415 transitions.append( 416 Transition( 417 origin_state_idx=origin, 418 target_state_idx=target, 419 prob=( 1.0 - component_residency_factor ), 420 symbol_idx=symbol_idx_map[ symbol ], 421 pq=None 422 ) ) 423 424 print( f"Was able to add {extra_edges_added}/{len(extra_edges)} extra cross-component transitions beyond a Hamiltonian cycle" ) 425 426 427 return HMM( 428 states=states, 429 transitions=transitions, 430 start_state=0, 431 alphabet=alphabet 432 )