diff --git a/causality-and-k-best/README.md b/causality-and-k-best/README.md new file mode 100644 index 0000000..4622f8d --- /dev/null +++ b/causality-and-k-best/README.md @@ -0,0 +1,42 @@ +# Causality and K-best + +This directory is a self-contained example based on the current Butadiene code. +Nothing in `ILP/butadien` or `mod/butadien` is changed. + +## 1. Enumerate K best supports with Gurobi + +Set `K` near the top of `butadiensynthesis_k_best.py` and run: + +```bash +python3 butadiensynthesis_k_best.py +``` + +The script is the existing Butadiene ILP with one small generalisation: instead of +computing only the best and second-best support, it repeatedly excludes every +previously found support. It writes all solutions to `gurobi_flows.json`. + +Thus "K best" here means K best distinct reaction supports, matching the meaning +of the existing second-best code. + +## 2. Check the Gurobi solutions for causality + +Activate the MØD environment and run this directory with MØD: + +```bash +mod -f check_causality.py +``` + +The checker reads `gurobi_flows.json`, maps each reaction to the current MØD +derivation graph by educt/product SMILES, fixes exactly that hyperflow, and runs +`causality.RealisabilityQuery`. + +Typical output is: + +```text +solution 1: causally realisable: YES +solution 2: causally realisable: NO +... +``` + +`butadienNetwork.py` and `graph.py` are local helper copies for constructing +the same Butadiene MØD network. diff --git a/causality-and-k-best/butadienNetwork.py b/causality-and-k-best/butadienNetwork.py new file mode 100644 index 0000000..ee86401 --- /dev/null +++ b/causality-and-k-best/butadienNetwork.py @@ -0,0 +1,232 @@ +config.common.ignoreDeprecation = False + +config.ilp.solver="CPLEX" +#import mod +#from mod import * +import networkx as nx + +# import Graph class from graph.py +from graph import GraphObj + +#ethylen = Graph.fromGMLString( +"""graph [ + node [ id 0 label "H" ] + node [ id 1 label "H" ] + node [ id 2 label "C" ] + node [ id 3 label "C" ] + node [ id 4 label "H" ] + node [ id 5 label "H" ] + edge [ source 0 target 2 label "-" ] + edge [ source 1 target 2 label "-" ] + edge [ source 2 target 3 label "=" ] + edge [ source 3 target 4 label "-" ] + edge [ source 3 target 5 label "-" ] +]""" +#, name="Ethylen") + +butadien = Graph.fromGMLString( +"""graph [ + node [ id 0 label "C" ] + node [ id 1 label "C" ] + node [ id 2 label "C" ] + node [ id 3 label "C" ] + node [ id 4 label "H" ] + node [ id 5 label "H" ] + node [ id 6 label "H" ] + node [ id 7 label "H" ] + node [ id 8 label "H" ] + node [ id 9 label "H" ] + edge [ source 0 target 1 label "=" ] + edge [ source 1 target 2 label "-" ] + edge [ source 2 target 3 label "=" ] + edge [ source 0 target 4 label "-" ] + edge [ source 0 target 5 label "-" ] + edge [ source 1 target 6 label "-" ] + edge [ source 2 target 7 label "-" ] + edge [ source 3 target 8 label "-" ] + edge [ source 3 target 9 label "-" ] +]""" +, name="Butadien") + +#benzene = Graph.fromSMILES('c1=cc=cc=c1') + +#benzene = Graph.fromGMLString( +"""graph [ + node [ id 0 label "C" ] + node [ id 1 label "H" ] + node [ id 2 label "C" ] + node [ id 3 label "H" ] + node [ id 4 label "C" ] + node [ id 5 label "H" ] + node [ id 6 label "C" ] + node [ id 7 label "H" ] + node [ id 8 label "C" ] + node [ id 9 label "H" ] + node [ id 10 label "C" ] + node [ id 11 label "H" ] + edge [ source 0 target 1 label "-" ] + edge [ source 0 target 2 label "=" ] + edge [ source 2 target 3 label "-" ] + edge [ source 2 target 4 label "-" ] + edge [ source 4 target 5 label "-" ] + edge [ source 4 target 6 label "=" ] + edge [ source 6 target 7 label "-" ] + edge [ source 6 target 8 label "-" ] + edge [ source 8 target 9 label "-" ] + edge [ source 8 target 10 label "=" ] + edge [ source 10 target 11 label "-" ] + edge [ source 10 target 0 label "-" ] +]""" +#, name="Benzene") + + +#pentadien = Graph.fromGMLString( +"""graph +[ + node [ id 0 label "C" ] + node [ id 1 label "C" ] + node [ id 2 label "C" ] + node [ id 3 label "C" ] + node [ id 4 label "H" ] + node [ id 5 label "H" ] + node [ id 6 label "H" ] + node [ id 7 label "H" ] + node [ id 8 label "H" ] + node [ id 9 label "C" ] + node [ id 10 label "H" ] + node [ id 11 label "H" ] + node [ id 12 label "H" ] + edge [ source 0 target 1 label "=" ] + edge [ source 1 target 2 label "-" ] + edge [ source 2 target 3 label "=" ] + edge [ source 0 target 4 label "-" ] + edge [ source 0 target 5 label "-" ] + edge [ source 1 target 6 label "-" ] + edge [ source 2 target 7 label "-" ] + edge [ source 3 target 8 label "-" ] + edge [ source 3 target 9 label "-" ] + edge [ source 9 target 10 label "-" ] + edge [ source 9 target 11 label "-" ] + edge [ source 9 target 12 label "-" ] +] +""" +#, name="Pentadien") + +restswap = Rule.fromGMLString( +"""rule [ + left [ + edge [ source 1 target 2 label "=" ] + edge [ source 3 target 4 label "=" ] + ] + context [ + node [ id 1 label "C" ] + node [ id 2 label "C"] + node [ id 3 label "C"] + node [ id 4 label "C"] + ] + right [ + edge [ source 1 target 3 label "=" ] + edge [ source 2 target 4 label "=" ] + ] +]""" +) + +dielsalder = Rule.fromGMLString( +"""rule [ + left [ + edge [ source 1 target 2 label "=" ] + edge [ source 2 target 3 label "-" ] + edge [ source 3 target 4 label "=" ] + edge [ source 5 target 6 label "=" ] + + ] + context [ + node [ id 1 label "C" ] + node [ id 2 label "C"] + node [ id 3 label "C"] + node [ id 4 label "C"] + node [ id 5 label "C"] + node [ id 6 label "C"] + ] + right [ + edge [ source 1 target 2 label "-" ] + edge [ source 2 target 3 label "=" ] + edge [ source 3 target 4 label "-" ] + edge [ source 4 target 5 label "-" ] + edge [ source 5 target 6 label "-" ] + edge [ source 6 target 1 label "-" ] + ] +]""" +) + +def cyclesizes(g): + nxGraph = GraphObj(g).nx_graph + + #Restriction to E/Z? Benzen can't be created or consumed by the reaction due to steric difficulties + + #Chain restrictions: (No steric reason) + for nodestart in nxGraph.nodes: + for nodeend in nxGraph.nodes: + if nodeend != nodestart: + paths = sorted(list(nx.all_shortest_paths(nxGraph, nodestart, nodeend))) + for path in paths: + if len(path) > 10: #First node included in length + return True + + #Cycle restrictions: + cycles = sorted(list(nx.chordless_cycles(nxGraph))) + cycle_lengths = [len(x) for x in cycles] + if cycles == []: + return False + #One Atom can't be in four different cycles + for vertice in nxGraph.nodes: + counter = 0 + for cycle in cycles: + if vertice in cycle: + counter += 1 + if counter >= 4: + return True + #One Cycle can't overlapp with another on over 2 Connection points + for cycle in cycles: + for cycleref in cycles: + if cycle != cycleref and len(list(set(cycle) & set(cycleref))) >= 3: + return True + #Only cordless cycles of length 5,6 and 7 are acceptable + if min(cycle_lengths) >= 6 and max(cycle_lengths) <=6: + return False + return True + +#Disallows Allenes (Two Doublebonds on same carbon) +def doubledoublebond(g): + for vertice in g.vertices: + if (vertice.stringLabel == "C" and vertice.degree == 2): # and vertice.edges.label == ["=", "="] + for edge in vertice.incidentEdges: + print(type(edge.bondType)) + bonds = [type(edge.bondType) for edge in vertice.incidentEdges] + if (bonds[0] == bonds[1] and len(bonds) == 2): + return True + return False + + +def restriction(dg): + for a in dg.right: + if a.vLabelCount("C") > 10: + return False + if doubledoublebond(a): + return False + if cyclesizes(a): + return False + return True + +ethylen = Graph.fromSMILES('C=C') +goal = Graph.fromSMILES('C1C=CC=CC=1') +dg = DG(graphDatabase=inputGraphs) +dg.build().execute( + addSubset(inputGraphs) + >> rightPredicate[ + restriction + ]( + repeat(revive(inputRules)) #Revive not necessary + ) +) + diff --git a/causality-and-k-best/butadiensynthesis_k_best.py b/causality-and-k-best/butadiensynthesis_k_best.py new file mode 100644 index 0000000..0bb015d --- /dev/null +++ b/causality-and-k-best/butadiensynthesis_k_best.py @@ -0,0 +1,367 @@ +import json +from pathlib import Path + +import numpy as np +import gurobipy as gp +from gurobipy import GRB, Model, quicksum + +HYPERGRAPH = { + 0: (['Butadien'], []), + 2: (['p_{0,0}'], []), + 3: (['p_{0,1}'], []), + 4: (['Butadien', 'Butadien'], ['p_{0,0}', 'p_{0,1}']), + 6: (['p_{0,2}'], []), + 7: (['Butadien', 'Butadien'], ['p_{0,2}']), + 10: (['p_{0,0}', 'p_{0,1}'], ['Butadien', 'Butadien']), + 12: (['p_{0,3}'], []), + 13: (['p_{0,0}', 'p_{0,2}'], ['p_{0,3}']), + 17: (['p_{0,4}'], []), + 18: (['p_{0,1}', 'p_{0,1}'], ['Butadien', 'p_{0,4}']), + 19: (['p_{0,5}'], []), + 20: (['p_{0,1}', 'p_{0,2}'], ['Butadien', 'p_{0,5}']), + 23: (['Butadien', 'p_{0,1}'], ['p_{0,0}', 'p_{0,4}']), + 26: (['Butadien', 'p_{0,2}'], ['p_{0,0}', 'p_{0,5}']), + 27: (['p_{0,6}'], []), + 28: (['p_{0,0}', 'p_{0,1}'], ['p_{0,6}']), + 29: (['p_{0,7}'], []), + 30: (['Butadien', 'p_{0,0}'], ['p_{0,7}']), + 31: (['p_{0,8}'], []), + 32: (['Butadien', 'p_{0,1}'], ['p_{0,8}']), + 33: (['Butadien', 'p_{0,1}'], ['p_{0,5}']), + 34: (['p_{0,9}'], []), + 35: (['Butadien', 'p_{0,1}'], ['p_{0,9}']), + 36: (['p_{0,10}'], []), + 37: (['Butadien', 'p_{0,1}'], ['p_{0,10}']), + 38: (['p_{0,3}'], ['p_{0,0}', 'p_{0,2}']), + 41: (['p_{0,11}'], []), + 42: (['p_{0,4}'], ['p_{0,0}', 'p_{0,11}']), + 43: (['p_{0,12}'], []), + 44: (['p_{0,5}'], ['p_{0,12}']), + 63: (['p_{0,13}'], []), + 64: (['p_{0,4}', 'p_{0,6}'], ['p_{0,1}', 'p_{0,13}']), + 65: (['p_{0,2}', 'p_{0,4}'], ['p_{0,1}', 'p_{0,5}']), + 66: (['p_{0,0}', 'p_{0,4}'], ['Butadien', 'p_{0,1}']), + 67: (['Butadien', 'p_{0,4}'], ['p_{0,1}', 'p_{0,1}']), + 77: (['p_{0,5}', 'p_{0,6}'], ['p_{0,2}', 'p_{0,13}']), + 79: (['p_{0,0}', 'p_{0,5}'], ['Butadien', 'p_{0,2}']), + 80: (['Butadien', 'p_{0,5}'], ['p_{0,1}', 'p_{0,2}']), + 83: (['p_{0,1}', 'p_{0,5}'], ['p_{0,2}', 'p_{0,4}']), + 96: (['Butadien', 'p_{0,6}'], ['p_{0,0}', 'p_{0,13}']), + 97: (['p_{0,1}', 'p_{0,6}'], ['Butadien', 'p_{0,13}']), + 99: (['p_{0,14}'], []), + 100: (['p_{0,0}', 'p_{0,6}'], ['p_{0,14}']), + 101: (['p_{0,15}'], []), + 102: (['p_{0,0}', 'p_{0,7}'], ['p_{0,15}']), + 121: (['p_{0,0}', 'p_{0,4}'], ['p_{0,13}']), + 122: (['p_{0,16}'], []), + 123: (['p_{0,0}', 'p_{0,4}'], ['p_{0,16}']), + 124: (['p_{0,17}'], []), + 125: (['Butadien', 'p_{0,7}'], ['p_{0,17}']), + 127: (['p_{0,12}'], ['p_{0,5}']), + 130: (['p_{0,14}'], ['p_{0,0}', 'p_{0,6}']), + 132: (['p_{0,15}'], ['p_{0,0}', 'p_{0,7}']), + 135: (['p_{0,0}', 'p_{0,11}'], ['p_{0,4}']), + 153: (['Butadien', 'p_{0,13}'], ['p_{0,1}', 'p_{0,6}']), + 155: (['p_{0,0}', 'p_{0,13}'], ['Butadien', 'p_{0,6}']), + 157: (['p_{0,1}', 'p_{0,13}'], ['p_{0,4}', 'p_{0,6}']), + 158: (['p_{0,2}', 'p_{0,13}'], ['p_{0,5}', 'p_{0,6}']), + 210: (['p_{0,18}'], []), + 211: (['Butadien', 'p_{0,11}'], ['p_{0,18}']), + 213: ([], ['Butadien']), +} + +HYPERGRAPH2 = { + 2: (['p_{0,0}'], []), + 3: (['p_{0,1}'], []), + 4: (['Butadien', 'Butadien'], ['p_{0,0}', 'p_{0,1}']), + 6: (['p_{0,2}'], []), + 7: (['Butadien', 'Butadien'], ['p_{0,2}']), + 17: (['p_{0,4}'], []), + 23: (['Butadien', 'p_{0,1}'], ['p_{0,0}', 'p_{0,4}']), + 41: (['p_{0,11}'], []), + 42: (['p_{0,4}'], ['p_{0,0}', 'p_{0,11}']), + 213: ([], ['Butadien']), +} + +HYPERGRAPH3 = { + 0: (['Butadien'], []), + 213: ([], ['Butadien']), + 4: (['Butadien', 'Butadien'], ['p_{0,0}', 'Butadien']), + 2: (['p_{0,0}'], []), +} + + +FIXED_FLOWS = { + 213: 3, + #2: 3, + 4: 1, + #23: 1, + #41: 1, + #42: 1, +} + +K = 10 + +def build_model(name, hyperedges, elmax, el1, el2, el3, excluded_supports=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} + n1 = {e_id: model.addVar(vtype=GRB.BINARY, name = f"n1_{e_id}") for e_id in hyperedges} + n2 = {e_id: model.addVar(vtype=GRB.BINARY, name = f"n2_{e_id}") for e_id in hyperedges} + n3 = {e_id: model.addVar(vtype=GRB.BINARY, name = f"n3_{e_id}") for e_id in hyperedges} + + + vertices = set(v for tails, heads in hyperedges.values() for v in tails + heads) + + #Every item created has to be consumed: + ''' + for v in vertices: + count[v] = 0 + count = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, name = "count") + for mol in vertices: + for edge, (tails, heads) in hyperedges.items(): + count[mol] += x[edge] * (tails.count(mol) - heads.count(mol)) + model.addConstr(count[mol] == 0) + ''' + + for v in vertices: + inflow = quicksum(x[e_id] * heads.count(v) for e_id, (_, heads) in hyperedges.items()) + outflow = quicksum(x[e_id] * tails.count(v) for e_id, (tails, _) in hyperedges.items()) + model.addConstr(inflow == outflow, name = f"flow_conservation_{v}") + + 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: + model.addGenConstrIndicator(b[e_id], 0, x[e_id] == 0, name = f"unused_implies_zero_{e_id}") + model.addConstr(x[e_id] >= b[e_id], name = f"used_implies_positive_flow_{e_id}") + #Only if an edge has flow, can it contribute to the likelihood + model.addConstr(b[e_id] >= n1[e_id], name = f"only_used_contribute_to_nmr1_{e_id}") + model.addConstr(b[e_id] >= n2[e_id], name = f"only_used_contribute_to_nmr2_{e_id}") + model.addConstr(b[e_id] >= n3[e_id], name = f"only_used_contribute_to_nmr3_{e_id}") + #One reaction can only contribute once + model.addConstr(n1[e_id] + n2[e_id] + n3[e_id] <= 1) + + #Only one edge per nmr can contribute + model.addConstr(quicksum(n1[e_id] for e_id in hyperedges) == 1) + model.addConstr(quicksum(n2[e_id] for e_id in hyperedges) == 1) + model.addConstr(quicksum(n3[e_id] for e_id in hyperedges) == 1) + + reaction_path = {} + + + if excluded_supports: + for i, support in enumerate(excluded_supports): + model.addConstr( + quicksum(b[e_id] for e_id in support) <= len(support) - 1, + name=f"different_hyperedges_{i}", + ) + + model.setObjective(quicksum(1000 * (el1[e_id] * n1[e_id] + el2[e_id] * n2[e_id] + el3[e_id] * n3[e_id]) - (1 / elmax[e_id]) * x[e_id] for e_id in hyperedges if elmax[e_id] != 0),GRB.MAXIMIZE) + + + #model.setObjective(quicksum(1000 * elmax[e_id] * b[e_id] - (1 - elmax[e_id]) * x[e_id] for e_id in hyperedges),GRB.MAXIMIZE) + #model.setObjective(quicksum(1000 * elmax[e_id] * b[e_id] - (1 / elmax[e_id]) * x[e_id] for e_id in hyperedges if elmax[e_id] != 0),GRB.MAXIMIZE) + #model.setObjective(quicksum(1000 * elmax[e_id] * b[e_id] + np.log(elmax[e_id]) * x[e_id] for e_id in hyperedges if elmax[e_id] != 0),GRB.MAXIMIZE) + #Excluding creation and destruction only three reactions for three nmr + #model.addConstr(quicksum(b[e_id] for e_id, (tails, heads) in hyperedges.items() if tails != [] and heads != []) == 3) + + #2 Butadien create first different molecule and it has to be created first: + startmolecule = ["Butadien"] + + #Including it or not changes first and second solution (sometimes flipped) + #model.addConstr(quicksum(b[e_id] for e_id, (tails, _) in hyperedges.items() if list(set(tails)) == startmolecule) == 1) + #model.addConstr(b[4] + b[7] == 1) + + model.addConstr(quicksum(b[e_id] for e_id, (tails, heads) in hyperedges.items() if tails == [] and heads == startmolecule)== 1) + #model.addConstr(b[213] == 1) + + #model.addConstr(b[41] == 1) + + + + #No cyclic reaction pairs: + for e_id1, (tails1, heads1) in hyperedges.items(): + for e_id2, (tails2, heads2) in hyperedges.items(): + if (heads1, tails1) == (tails2, heads2): + #print(e_id1, e_id2) + model.addConstr((b[e_id1] + b[e_id2]) <= 1) + + return model, x, b + + +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 edgelikelihoods(variable_dict, edgelikelihood, threshold = 0.5): + return {e_id: edgelikelihood[e_id] for e_id, var in variable_dict.items() if var.X > threshold} + +def flow_to_json(flow_solution, smiles): + inputs = {} + outputs = {} + reactions = [] + + for e_id, value in flow_solution.items(): + value = int(round(value)) + tails, heads = HYPERGRAPH[e_id] + + if tails == []: + for molecule in heads: + inputs[smiles[molecule]] = inputs.get(smiles[molecule], 0) + value + elif heads == []: + for molecule in tails: + outputs[smiles[molecule]] = outputs.get(smiles[molecule], 0) + value + else: + reactions.append({ + "flow": value, + "educts": [smiles[molecule] for molecule in tails], + "products": [smiles[molecule] for molecule in heads], + }) + + return { + "inputs": inputs, + "outputs": outputs, + "reactions": reactions, + } + +def write_flows_json(flows): + path = Path(__file__).with_name("gurobi_flows.json") + path.write_text(json.dumps(flows, indent=2) + "\n") + print(f"\nWrote {path}") + +def print_solution(title, flow_solution, binary_solution, edge_likelihoods, hyperedges): + print(f"\n{title}:") + for e_id in sorted(flow_solution): + flow = flow_solution[e_id] + tails, heads = hyperedges[e_id] + likelihood = edge_likelihoods[e_id] + print(f"Hyperedge {e_id}: Flow = {flow}, Tails = {tails}, Heads = {heads}, With likelihood = {likelihood}") + print("\nBinary Variables:") + for e_id in sorted(binary_solution): + print(f"Binary Variable b_{e_id} = {binary_solution[e_id]}") + + print(f"\nTotal flow: {sum(flow_solution.values())}") + print(f"Number of used hyperedges: {len(binary_solution)}") + +def main(): + VERTICES = ['Butadien', 'p_{0,0}', 'p_{0,1}', 'p_{0,2}', 'p_{0,3}', 'p_{0,4}', 'p_{0,5}', 'p_{0,6}', 'p_{0,7}', 'p_{0,8}', 'p_{0,9}', 'p_{0,10}', 'p_{0,11}', 'p_{0,12}', 'p_{0,13}', 'p_{0,14}', 'p_{0,15}', 'p_{0,16}', 'p_{0,17}', 'p_{0,18}'] + #Vergleich mit dem NMR von Ethylen und Hexatrien + NMR1 = [0.32, 0.5, 0.87, 0.0, 0.11, 0.58, 0.06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.09, 0.0, 0.0, 0.17, 0.0, 0.09] + #Vergleich mit dem NMR von Ethylen und Octrien + NMR2 = [0.33, 0.43, 0.6, 0.0, 0.13, 0.9, 0.07, 0.01, 0.0, 0.0, 0.01, 0.0, 0.0, 0.01, 0.1, 0.0, 0.0, 0.22, 0.0, 0.13] + #Vergleich mit dem NMR von Ethylen und Benzol + NMR3 = [0.0, 0.32, 0.0, 0.0, 0.0, 0.0, 0.05, 0.0, 0.0, 0.02, 0.0, 0.29, 0.95, 0.17, 0.0, 0.0, 0.0, 0.0, 0.0, 0.05] + VERTICESSMILES = ['C=CC=C', 'C=C', 'C=CC=CC=C', 'C1CCC(C=C)CC=1', 'C(CCC(C=C)CC=C)=C', 'C=CC=CC=CC=C', 'C1CCC(C=CC=C)CC=1', 'C=CC1C=CCCC1', 'C1CCCCC=1', 'C=CC1CC=CCC1C=C', 'C1CCC(C=C)C(C=C)C=1', 'C(C1CC(C=C)C=CC1)=C', 'C1C=CC=CC=1', 'C(CCC1C=CC=CC1)=C', 'C1C(C=CC=C)CCCC=1', 'C=CC(C=C)CCCC=C', 'C=CCCCCC=C', 'C=CC1C=CC(C=C)CC1', 'C1CC2CCCCC2CC=1', 'C1CC2C=CC=CC2CC=1'] + SMILES = dict(zip(VERTICES, VERTICESSMILES)) + + #Chosable parameters + modes = ["Product", "Average"] + mode = modes[1] + normalize = True + + if normalize: + NMR1 = [round(l/sum(NMR1), 2) for l in NMR1] + NMR2 = [round(l/sum(NMR2), 2) for l in NMR2] + NMR3 = [round(l/sum(NMR3), 2) for l in NMR3] + #print(NMR1) + #print(NMR2) + #print(NMR3) + + + #Kombiniert Molekül it Wahrscheinlichkeit für NMR1: + VERTICE1 = {} + for vertice, likelihood in zip(VERTICES, NMR1): + VERTICE1[vertice] = likelihood + #Kombiniert Molekül it Wahrscheinlichkeit für NMR2: + VERTICE2 = {} + for vertice, likelihood in zip(VERTICES, NMR2): + VERTICE2[vertice] = likelihood + #Kombiniert Molekül it Wahrscheinlichkeit für NMR3: + VERTICE3 = {} + for vertice, likelihood in zip(VERTICES, NMR3): + VERTICE3[vertice] = likelihood + + #Kantenwahrscheinlichkeiten für NMR1: + EDGE1 = {} + #Kantenwahrscheinlichkeiten für NMR2: + EDGE2 = {} + #Kantenwahrscheinlichkeiten für NMR3: + EDGE3 = {} + #Maximum der verschiedenen Kantenwahrscheinlichkeiten: + EDGEMAX = {} + + #Ist es besser, die Edukte mit zu berücksichtigen? Dann hin rück gleich wahrscheinlich, aber nur 1 kann gewählt werden + #Ordnet die Kantenwahrscheinlichkeit basierend auf dem Mittel oder Produkt der Wahrscheinlichkeiten der Produkte + for edge, (tails, heads) in HYPERGRAPH.items(): + if heads == [] or tails == []: + EDGE1[edge] = 0.0 + EDGE2[edge] = 0.0 + EDGE3[edge] = 0.0 + elif mode == "Product": + #print([VERTICE1[head] for head in heads]) + EDGE1[edge] = round(np.prod([VERTICE1[head] for head in heads]),2) + EDGE2[edge] = round(np.prod([VERTICE2[head] for head in heads]),2) + EDGE3[edge] = round(np.prod([VERTICE3[head] for head in heads]),2) + elif mode == "Average": + EDGE1[edge] = round(sum(VERTICE1[head] for head in heads)/len(heads),2) + EDGE2[edge] = round(sum(VERTICE2[head] for head in heads)/len(heads),2) + EDGE3[edge] = round(sum(VERTICE3[head] for head in heads)/len(heads),2) + + EDGEMAX[edge] = max([EDGE1[edge], EDGE2[edge], EDGE3[edge]]) + #print(EDGEMAX[edge], HYPERGRAPH[edge]) + + ''' + #Vertices from example: + verticelist = [0,1,2,3,4,5,11,12] + print([VERTICES[a] for a in verticelist]) + print([NMR1[a] for a in verticelist]) + print([NMR2[a] for a in verticelist]) + print([NMR3[a] for a in verticelist]) + #Edges from example + print(EDGEMAX[4], HYPERGRAPH[4]) + print(EDGEMAX[7], HYPERGRAPH[7]) + print(EDGEMAX[13], HYPERGRAPH[13]) + print(EDGEMAX[23], HYPERGRAPH[23]) + print(EDGEMAX[37], HYPERGRAPH[37]) + print(EDGEMAX[42], HYPERGRAPH[42]) + print(EDGEMAX[67], HYPERGRAPH[67]) + ''' + + excluded_supports = [] + flows = [] + + for rank in range(1, K + 1): + model, x, b = build_model( + f"HypergraphFlow_{rank}", + HYPERGRAPH, + EDGEMAX, + EDGE1, + EDGE2, + EDGE3, + excluded_supports=excluded_supports, + ) + model.optimize() + + if model.status != GRB.Status.OPTIMAL: + break + + solution = positive_entries(x) + binary_solution = positive_entries(b) + likelihoods = edgelikelihoods(b, EDGEMAX) + + print_solution( + f"Solution {rank}", + solution, + binary_solution, + likelihoods, + HYPERGRAPH, + ) + + flows.append(flow_to_json(solution, SMILES)) + excluded_supports.append(list(binary_solution.keys())) + + write_flows_json(flows) + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/causality-and-k-best/check_causality.py b/causality-and-k-best/check_causality.py new file mode 100644 index 0000000..6d45713 --- /dev/null +++ b/causality-and-k-best/check_causality.py @@ -0,0 +1,114 @@ +import json +from collections import Counter + +include("butadienNetwork.py") + +with open("gurobi_flows.json") as f: + gurobi_flows = json.load(f) + + +def canonical_smiles(smiles): + return str(Graph.fromSMILES(smiles).smiles) + + +def vertex_key(v): + return str(v.graph.smiles) + + +vertices_by_smiles = {} +for v in dg.vertices: + vertices_by_smiles.setdefault(vertex_key(v), []).append(v) + + +def current_vertex(smiles): + matches = vertices_by_smiles.get(canonical_smiles(smiles), []) + if len(matches) != 1: + raise RuntimeError( + f"{smiles}: expected one MØD vertex, found {len(matches)}" + ) + return matches[0] + + +def current_reaction(educts, products): + wanted_educts = Counter(canonical_smiles(x) for x in educts) + wanted_products = Counter(canonical_smiles(x) for x in products) + + matches = [] + for e in dg.edges: + current_educts = Counter(vertex_key(v) for v in e.sources) + current_products = Counter(vertex_key(v) for v in e.targets) + if current_educts == wanted_educts and current_products == wanted_products: + matches.append(e) + + if len(matches) != 1: + raise RuntimeError( + f"{educts} -> {products}: expected one MØD edge, found {len(matches)}" + ) + return matches[0] + + +def check_causality(gurobi_flow, rank): + source_flow = { + current_vertex(smiles): value + for smiles, value in gurobi_flow["inputs"].items() + } + output_flow = { + current_vertex(smiles): value + for smiles, value in gurobi_flow["outputs"].items() + } + + reaction_flow = {} + for reaction in gurobi_flow["reactions"]: + edge = current_reaction(reaction["educts"], reaction["products"]) + reaction_flow[edge] = reaction_flow.get(edge, 0) + reaction["flow"] + + flow = hyperflow.Model(dg) + + for v in source_flow: + flow.addSource(v.graph) + for v in output_flow: + flow.addSink(v.graph) + + for e in dg.edges: + flow.addConstraint(edgeFlow[e] == reaction_flow.get(e, 0)) + + for v in dg.vertices: + flow.addConstraint(inFlow[v] == source_flow.get(v, 0)) + flow.addConstraint(outFlow[v] == output_flow.get(v, 0)) + + flow.findSolutions(maxNumSolutions=1) + + try: + solution = next(iter(flow.solutions)) + except StopIteration: + print(f"solution {rank}: not a valid MØD hyperflow") + return "INVALID" + + query = causality.RealisabilityQuery(flow.dg) + dag = query.findDAG(solution) + + if not dag: + print(f"solution {rank}: causally realisable: NO") + return "NO" + + print(f"solution {rank}: causally realisable: YES") + + printer = DGPrinter() + printer.graphvizPrefix = 'layout = "dot";' + + def vertex_visible(v): + return solution.eval(vertexFlow[v]) != 0 + + printer.pushVertexVisible(vertex_visible) + dg.print(printer=printer, data=dag.getPrintData(False)) + + return "YES" + + +results = [] +for rank, gurobi_flow in enumerate(gurobi_flows, start=1): + results.append(check_causality(gurobi_flow, rank)) + +print("\n=== causality summary ===") +for rank, result in enumerate(results, start=1): + print(f"solution {rank}: {result}") diff --git a/causality-and-k-best/graph.py b/causality-and-k-best/graph.py new file mode 100644 index 0000000..c432191 --- /dev/null +++ b/causality-and-k-best/graph.py @@ -0,0 +1,136 @@ +import mod +import networkx as nx +import time + +# from equilibrator_assets.generate_compound import create_compounds + +# each vertex has attribute "label" which stores string label +# each edge has attribute "bond" which stores bond order (-, =, etc.) + +class GraphObj: + def __init__(self, graph, verbose=False): + self._verbose = verbose + + mod_lg_class = mod.libpymod.Rule.LeftGraph + mod_rg_class = mod.libpymod.Rule.RightGraph + mod_g_class = mod.libpymod.Graph + if isinstance(graph, mod_lg_class) or isinstance(graph, mod_rg_class) or isinstance(graph, mod_g_class): + self._nx_graph = self.mod_to_nx_graph(graph) + self._gml_string = self.nx_graph_to_GML_string(self._nx_graph) + elif isinstance(graph, nx.Graph): + self._nx_graph = graph + self._gml_string = self.nx_graph_to_GML_string(self._nx_graph) + elif isinstance(graph, str): + self._gml_string = graph + self._nx_graph = self.GML_to_nx_graph(self._gml_string) + else: + print(f"ERROR: Graph class cannot identify graph type in constructor: {type(graph)}") + + # up to the caller to check the number of components created + ccs = [self._nx_graph.subgraph(c).copy() for c in nx.connected_components(self._nx_graph)] + ccs = sorted(ccs, key=len, reverse=True) + self._num_components = len(ccs) + if self._num_components > 1: + self._components = [Graph(c) for c in ccs] + else: + self._components = [self] + + if self._num_components == 1: + self._mod_graph = self.nx_graph_to_mod(self._nx_graph) + + # self._equ_compounds = self.equ_compounds() + + def equ_compounds(self): + smiles = [c.mod_graph.smiles for c in self._components] + equ_comps = create_compounds(smiles, mol_format="smiles", bypass_chemaxon=True, save_empty_compounds=True) + return equ_comps + + def mod_to_nx_graph(self, modGraph: mod.Graph): + g = nx.Graph() + + for v in modGraph.vertices: + g.add_node(int(v.id), label=str(v.stringLabel), modID=int(v.id)) + + for e in modGraph.edges: + g.add_edge(int(e.source.id), int(e.target.id), bond=str(e.bondType)) + + return g + + def nx_graph_to_mod(self, nxGraph): + try: + return mod.Graph.fromGMLString(self.nx_graph_to_GML_string(nxGraph)) + except mod.libpymod.InputError: # graph is not connected probably + if self._verbose: + print("Error converting nxGraph to mod graph. Likely graph is not connected. This will not affect rule generation.") + return None + + def GML_to_nx_graph(self, gml_string): + g = nx.Graph() + lines = gml_string.split("\n") + for line in lines: + tokens = line.split() + if len(tokens) < 2: + continue + if tokens[0] == "node": + (_,_,_, mid, _, l, _) = tokens + g.add_node(int(mid), label=l[1:-1], modID=int(mid)) + elif tokens[0] == "edge": + (_, _, _, u, _, v, _, bondOrder, _ ) = tokens + g.add_edge(int(u), int(v), bond=bondOrder[1:-1]) + return g + + def nx_graph_to_GML_string(self, nxGraph): + out = [] + out.append("graph [") + + out.extend([f"\t\tnode [ id {nxGraph.nodes[node]['modID']} label \"{nxGraph.nodes[node]['label']}\" ]" for node in nxGraph.nodes]) + out.extend([f"\t\tedge [ source {u} target {v} label \"{nxGraph[u][v]['bond']}\" ]" + for (u,v) in nxGraph.edges]) + + out.append("]") + + return "\n".join(out) + + @property + def nx_graph(self): + return self._nx_graph + + @property + def mod_graph(self): + if self.num_components == 1: + return self._mod_graph + else: + return mod.Graph.fromGMLString(self.gml_string) + + @property + def gml_string(self): + return self._gml_string + + @property + def gml(self): + return self._gml_string + + @property + def edges(self): + return self._nx_graph.edges + + @property + def nodes(self): + return self._nx_graph.nodes + + @property + def num_components(self): + return self._num_components + + @property + def connected_components(self): + return self._components + + def mod_print(self): + self._mod_graph.print() + + def __str__(self) -> str: + return self._gml_string + + def __hash__(self) -> int: + return hash(self._gml_string)