diff --git a/mod/butadien/butadienCausality.py b/mod/butadien/butadienCausality.py
new file mode 100644
index 0000000..51dae99
--- /dev/null
+++ b/mod/butadien/butadienCausality.py
@@ -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}")