Dateien nach "mod/butadien" hochladen
This commit is contained in:
@@ -0,0 +1,330 @@
|
||||
config.ilp.solver="CPLEX"
|
||||
#import mod
|
||||
#from mod import *
|
||||
import networkx as nx
|
||||
|
||||
# import Graph class from graph.py
|
||||
from graph import GraphObj
|
||||
|
||||
# Search enough flows to include the known rare causal solution.
|
||||
MAX_SOLUTIONS = 10
|
||||
PRINT_FULL_PPP_OUTPUT = False
|
||||
|
||||
# Important: the default analysis checks causal realisability from the
|
||||
# explicitly declared source molecules only. We do not add catalysts/seeds
|
||||
# to rescue an otherwise non-realisable flow.
|
||||
USE_CATALYSTS = False
|
||||
|
||||
#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
|
||||
|
||||
flowPrinter = FlowPrinter()
|
||||
flowPrinter.printUnfiltered = False
|
||||
|
||||
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
|
||||
)
|
||||
)
|
||||
|
||||
|
||||
# Deliberately no automatic DG/product printing here. Use
|
||||
# PRINT_FULL_PPP_OUTPUT below if the complete PPP-style report is wanted.
|
||||
flow = Flow(dg)
|
||||
flow.addSource(butadien)
|
||||
|
||||
flow.addSink(ethylen)
|
||||
flow.addSink(goal)
|
||||
|
||||
flow.addConstraint(inFlow[butadien] >= 2)
|
||||
flow.addConstraint(outFlow[goal] == 1)
|
||||
|
||||
def printStuff(flowToPrint):
|
||||
# Same output structure as in the PPP example.
|
||||
post.summarySection("Loaded Graphs")
|
||||
for a in inputGraphs:
|
||||
a.print()
|
||||
|
||||
post.summarySection("Loaded Rules")
|
||||
for a in inputRules:
|
||||
a.print()
|
||||
|
||||
dg.print()
|
||||
|
||||
post.summarySection("Product Graphs")
|
||||
for v in dg.vertices:
|
||||
v.graph.print()
|
||||
|
||||
# Important: flowToPrint contains only the one selected causal solution.
|
||||
flowToPrint.solutions.print()
|
||||
|
||||
|
||||
def makeSingleSolutionFlow(solution):
|
||||
"""
|
||||
Copy the original hyperflow model and fix all core flow variables to the
|
||||
values of one already-computed solution. The copied model can then be
|
||||
printed with the normal PPP printStuff(), but contains only this solution.
|
||||
"""
|
||||
selectedFlow = hyperflow.Model(flow)
|
||||
|
||||
# Fix every reaction/hyperedge multiplicity.
|
||||
for e in dg.edges:
|
||||
selectedFlow.addConstraint(
|
||||
edgeFlow[e] == solution.eval(edgeFlow[e])
|
||||
)
|
||||
|
||||
# Also fix all external input/output multiplicities.
|
||||
for v in dg.vertices:
|
||||
selectedFlow.addConstraint(
|
||||
inFlow[v] == solution.eval(inFlow[v])
|
||||
)
|
||||
selectedFlow.addConstraint(
|
||||
outFlow[v] == solution.eval(outFlow[v])
|
||||
)
|
||||
|
||||
selectedFlow.findSolutions(maxNumSolutions=1)
|
||||
return selectedFlow
|
||||
|
||||
|
||||
# Keep the test small.
|
||||
flow.findSolutions(maxNumSolutions=MAX_SOLUTIONS)
|
||||
|
||||
numSolutions = 0
|
||||
numRealisable = 0
|
||||
|
||||
for s in flow.solutions:
|
||||
numSolutions += 1
|
||||
|
||||
# Exact PPP causality analysis / rendering.
|
||||
query = causality.RealisabilityQuery(flow.dg)
|
||||
printer = DGPrinter()
|
||||
printer.graphvizPrefix = 'layout = "dot";'
|
||||
def vertexVisible(v):
|
||||
return s.eval(vertexFlow[v]) != 0
|
||||
|
||||
printer.pushVertexVisible(vertexVisible)
|
||||
|
||||
dag = query.findDAG(s)
|
||||
|
||||
if dag:
|
||||
numRealisable += 1
|
||||
data = dag.getPrintData(False)
|
||||
dg.print(printer=printer, data=data)
|
||||
|
||||
print("\n=== causality summary ===")
|
||||
print(f"hyperflow solutions checked: {numSolutions}")
|
||||
print(f"causally realisable: {numRealisable}")
|
||||
print(f"not causally realisable: {numSolutions - numRealisable}")
|
||||
Reference in New Issue
Block a user