diff --git a/ILP/Vanilla/vanillasynthesis.py b/ILP/Vanilla/vanillasynthesis.py
index 6cf9da4..614af99 100644
--- a/ILP/Vanilla/vanillasynthesis.py
+++ b/ILP/Vanilla/vanillasynthesis.py
@@ -20,8 +20,8 @@ HYPERGRAPH = {
 }
 
 FIXED_FLOWS = {
-    #1: 1,
-    #14: 1,
+    1: 1,
+    14: 1,
 }
 
 def build_model(name, hyperedges, vertices, ele, el1, el2, excluded_support=None):
@@ -29,7 +29,9 @@ def build_model(name, hyperedges, vertices, ele, el1, el2, excluded_support=None
 
     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}
+
     vertices = set(v for tails, heads in hyperedges.values() for v in tails + heads)
 
     for v in vertices:
@@ -43,6 +45,15 @@ def build_model(name, hyperedges, vertices, ele, el1, el2, excluded_support=None
     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}")
+        #One reaction can only contribute once
+        model.addConstr(n1[e_id] +  n2[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)
 
     reaction_path = {}
 
@@ -51,9 +62,12 @@ def build_model(name, hyperedges, vertices, ele, el1, el2, excluded_support=None
         model.addConstr(quicksum(b[e_id] for e_id in excluded_support) <= len(excluded_support) - 1, name = "different_hyperedges",)
     
     #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)
+    #model.addConstr(quicksum(b[e_id] for e_id, (tails, heads) in hyperedges.items() if tails != [] and heads != []) == 3)
 
-    model.setObjective(quicksum(1000 * ele[e_id] * b[e_id] - x[e_id] for e_id in hyperedges),GRB.MAXIMIZE)
+    model.setObjective(quicksum(1000 * (ele[e_id] * b[e_id]) - (1 / ele[e_id]) * x[e_id] for e_id in hyperedges if ele[e_id] != 0),GRB.MAXIMIZE)
+
+    #ILP picks best pair
+    #model.setObjective(quicksum(1000 * (el1[e_id] * n1[e_id] + el2[e_id] * n2[e_id]) - (1 / ele[e_id]) * x[e_id] for e_id in hyperedges if ele[e_id] != 0),GRB.MAXIMIZE)
 
     return model, x, b
 
@@ -89,7 +103,7 @@ def main():
     #Chosable parameters
     modes = ["Product", "Average"]
     mode = modes[0] 
-    normalize = False
+    normalize = True
     
     if normalize:
         NMR1 = [round(l/sum(NMR1), 2) for l in NMR1]