import random import sys from gurobipy import GRB, Model, quicksum def set_random_seed(seed=None): if seed is not None: random.seed(seed) print(f"Random seed set to: {seed}") else: random.seed() print("Random seed set to system default.") def build_model(name, hyperedges, r_values, excluded_support=None): model = Model(name) x = {e_id: model.addVar(vtype=GRB.INTEGER, lb=0, name=f"x_{e_id}") for e_id in hyperedges} b = {e_id: model.addVar(vtype=GRB.BINARY, name=f"b_{e_id}") for e_id in hyperedges} z = model.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name="z") vertices = set(v for tails, heads in hyperedges.values() for v in tails + heads) for v in vertices: inflow = quicksum(x[e_id] for e_id, (_, heads) in hyperedges.items() if v in heads) outflow = quicksum(x[e_id] for e_id, (tails, _) in hyperedges.items() if v in tails) model.addConstr(inflow == outflow, name=f"flow_conservation_{v}") fixed_flows = { 41: 1, 42: 6, 43: 5, 44: 1, } for e_id, value in fixed_flows.items(): model.addConstr(x[e_id] == value, name=f"fixed_flow_{e_id}") for e_id in hyperedges: # If a hyperedge is not selected, its flow must be zero. model.addGenConstrIndicator(b[e_id], 0, x[e_id] == 0, name=f"unused_implies_zero_{e_id}") # If a hyperedge is selected, it must carry at least one unit of flow. model.addConstr(x[e_id] >= b[e_id], name=f"used_implies_positive_flow_{e_id}") # z upper-bounds the barrier of every used hyperedge. model.addConstr(z >= r_values[e_id] * b[e_id], name=f"barrier_link_{e_id}") if excluded_support: model.addConstr( quicksum(b[e_id] for e_id in excluded_support) <= len(excluded_support) - 1, name="different_hyperedges", ) # Primary objective: minimize the maximal barrier. model.ModelSense = GRB.MINIMIZE model.setObjectiveN(z, index=0, priority=2, name="minimize_max_barrier") # Secondary objective: among equally good barriers, minimize the number of used hyperedges. model.setObjectiveN( quicksum(b[e_id] for e_id in hyperedges), index=1, priority=1, name="minimize_used_hyperedges", ) return model, x, b, z def positive_entries(variable_dict, threshold=0.5): return {e_id: var.X for e_id, var in variable_dict.items() if var.X > threshold} def print_solution(title, solution, binary_solution, hyperedges, r_values, z_value, binary_prefix="b"): used_hyperedges = sorted(binary_solution) actual_max_barrier = max((r_values[e_id] for e_id in used_hyperedges), default=0.0) print(f"\n{title}:") for e_id in sorted(solution): flow = solution[e_id] tails, heads = hyperedges[e_id] print( f"Hyperedge {e_id}: Flow = {flow}, " f"Tails = {tails}, Heads = {heads}, r_value = {r_values[e_id]}" ) print("\nBinary Variables:") for e_id in sorted(binary_solution): print(f"Binary Variable {binary_prefix}_{e_id} = {binary_solution[e_id]}") print(f"\nObjective variable z: {z_value}") print(f"Max r_{{e_id}} value among used hyperedges: {actual_max_barrier}") print(f"Number of used hyperedges: {len(binary_solution)}") def main(seed=None): if seed is not None: set_random_seed(seed) else: set_random_seed(42) hyperedges = { 3: (["Ribulose-5-Phosphate"], ["p_{0,0}"]), 6: (["Ribulose-5-Phosphate", "p_{0,0}"], ["p_{0,1}", "p_{0,2}"]), 9: (["p_{0,1}", "p_{0,2}"], ["Fructose-6-Phosphate", "p_{0,3}"]), 11: (["p_{0,0}", "p_{0,1}"], ["p_{0,3}", "p_{0,4}"]), 13: (["Ribulose-5-Phosphate", "p_{0,2}"], ["Fructose-6-Phosphate", "p_{0,5}"]), 14: (["Fructose-6-Phosphate", "p_{0,3}"], ["Fructose-6-Phosphate", "p_{0,3}"]), 16: (["Fructose-6-Phosphate", "p_{0,5}"], ["p_{0,3}", "p_{0,6}"]), 17: (["Fructose-6-Phosphate", "p_{0,2}"], ["Ribulose-5-Phosphate", "p_{0,3}"]), 18: (["Fructose-6-Phosphate", "p_{0,0}"], ["p_{0,1}", "p_{0,3}"]), 20: (["p_{0,3}", "p_{0,4}"], ["Fructose-6-Phosphate", "p_{0,7}"]), 21: (["Ribulose-5-Phosphate", "p_{0,3}"], ["Fructose-6-Phosphate", "p_{0,2}"]), 22: (["p_{0,1}", "p_{0,3}"], ["Fructose-6-Phosphate", "p_{0,0}"]), 23: (["p_{0,4}", "p_{0,5}"], ["p_{0,6}", "p_{0,7}"]), 24: (["p_{0,2}", "p_{0,4}"], ["Ribulose-5-Phosphate", "p_{0,7}"]), 25: (["p_{0,0}", "p_{0,4}"], ["p_{0,1}", "p_{0,7}"]), 26: (["Ribulose-5-Phosphate", "p_{0,5}"], ["p_{0,2}", "p_{0,6}"]), 27: (["p_{0,1}", "p_{0,5}"], ["p_{0,0}", "p_{0,6}"]), 29: (["p_{0,2}"], ["p_{0,8}"]), 31: (["p_{0,0}", "p_{0,8}"], ["p_{0,9}"]), 33: (["p_{0,2}", "p_{0,8}"], ["p_{0,10}"]), 36: (["H2O", "p_{0,9}"], ["p_{0,11}", "p_{0,12}"]), 37: (["H2O", "p_{0,9}"], ["p_{0,4}", "p_{0,11}"]), 39: (["H2O", "p_{0,10}"], ["p_{0,11}", "p_{0,13}"]), 40: (["H2O", "p_{0,10}"], ["Fructose-6-Phosphate", "p_{0,11}"]), 41: ([], ["H2O"]), 42: ([], ["Ribulose-5-Phosphate"]), 43: (["Fructose-6-Phosphate"], []), 44: (["p_{0,11}"], []), } r_values = { e_id: random.random() if e_id not in [41, 42, 43, 44] else 0.0 for e_id in hyperedges } model, x, b, z = build_model("HypergraphFlow", hyperedges, r_values) model.optimize() optimal_solution = positive_entries(x) optimal_binary_solution = positive_entries(b) print_solution( "Optimal Solution", optimal_solution, optimal_binary_solution, hyperedges, r_values, z.X, binary_prefix="b", ) excluded_support = list(optimal_binary_solution.keys()) second_model, x2, b2, z2 = build_model( "SecondBestHypergraphFlow", hyperedges, r_values, excluded_support=excluded_support, ) second_model.optimize() if second_model.status == GRB.Status.OPTIMAL: second_solution = positive_entries(x2) second_binary_solution = positive_entries(b2) print_solution( "Second Best Solution", second_solution, second_binary_solution, hyperedges, r_values, z2.X, binary_prefix="b2", ) else: print("No optimal solution found for the second best model.") if __name__ == "__main__": cli_seed = int(sys.argv[1]) if len(sys.argv) > 1 else None main(cli_seed)