Add causality and K-best example
This commit is contained in:
Daniel
@@ -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.
|
||||
@@ -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
|
||||
)
|
||||
)
|
||||
|
||||
@@ -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()
|
||||
@@ -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}")
|
||||
@@ -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)
|
||||
Reference in New Issue
Block a user