diff --git a/ILP/butadien/butadiensynthesis.py b/ILP/butadien/butadiensynthesis.py index 54e658c..ebf65a1 100644 --- a/ILP/butadien/butadiensynthesis.py +++ b/ILP/butadien/butadiensynthesis.py @@ -103,7 +103,46 @@ NMR3 = [0.0, 0.17, 0.0, 0.0, 0.0, 0.0, 0.03, 0.0, 0.0, 0.01, 0.0, 0.16, 0.51, 0. 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'] NMRMAX = [max(n1,n2,n3) for n1,n2,n3 in zip(NMR1,NMR2,NMR3)] -#Implement Edgelikelihood before Model +#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 = {} + +modes = ["Product", "Average"] +mode = modes[0] +#Ordnet die Kantenwahrscheinlichkeit basierend auf dem Mittel oder Produkt der Wahrscheinlichkeiten der Produkte +for edge, (heads, tails) in HYPERGRAPH.items(): + if heads != [] and tails != [] and mode == "Product": + #Produkt klappt nur, wenn normalisierte Likelihoods + EDGE1[edge] = math.prod(VERTICE1[head] for head in heads) + EDGE2[edge] = math.prod(VERTICE2[head] for head in heads) + EDGE2[edge] = math.prod(VERTICE2[head] for head in heads) + if heads != [] and tails != [] and mode == "Average": + EDGE1[edge] = quicksum(VERTICE1[head] for head in heads)/len(heads) + EDGE2[edge] = quicksum(VERTICE2[head] for head in heads)/len(heads) + EDGE3[edge] = quicksum(VERTICE3[head] for head in heads)/len(heads) + else: + EDGE1[edge] = 0.0 + EDGE2[edge] = 0.0 + EDGE3[edge] = 0.0 + EDGEMAX[edge] = max([EDGE1[edge], EDGE2[edge], EDGE2[edge]]) FIXED_FLOWS = { @@ -115,64 +154,37 @@ FIXED_FLOWS = { #42: 1, } -def build_model(name, hyperedges, vertices, nmrlikelihoodsmax, nmrlikelihoods1, nmrlikelihoods2, nmrlikelihoods3, excluded_support=None): +def build_model(name, hyperedges, edgelikelihoodsmax, edgelikelihoods1, edgelikelihoods2, edgelikelihoods3, excluded_support=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} - count = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, name = "count") - nmax = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "nmrmax") - n1 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "nmr1") - n2 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "nmr2") - n3 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "nmr3") - enmax = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmrmax") - en1 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmr1") - en2 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmr2") - en3 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmr3") - - #Assigns every Molecule the likelihood compared to the three different reference spectra - for v, nmrmax, nmr1, nmr2, nmr3 in zip(vertices, nmrlikelihoodsmax, nmrlikelihoods1, nmrlikelihoods2, nmrlikelihoods3): - nmax[v] = nmrmax - n1[v] = nmr1 - n2[v] = nmr2 - n3[v] = nmr3 - count[v] = 0 - - #Assigns the edges a likelihood based on the products - for e, (tails, heads) in hyperedges.items(): - if heads != [] and tails != []: #Multiplication better? - #Produkt klappt nur, wenn normalisierte Likelihoods - en1[e] = math.prod(n1[head] for head in heads) - en2[e] = math.prod(n2[head] for head in heads) - en3[e] = math.prod(n3[head] for head in heads) - enmax[e] = max([en1[e], en2[e], en3[e]]) - else: - enmax[e] = 0.0 - en1[e] = 0.0 - en2[e] = 0.0 - en3[e] = 0.0 - ''' - if heads != [] and tails != []: - enmax[e] = quicksum(nmax[head] for head in heads)/len(heads) - en1[e] = quicksum(n1[head] for head in heads)/len(heads) - en2[e] = quicksum(n2[head] for head in heads)/len(heads) - en3[e] = quicksum(n3[head] for head in heads)/len(heads) - else: - enmax[e] = 0.0 - en1[e] = 0.0 - en2[e] = 0.0 - en3[e] = 0.0 - ''' + enmax = model.addVars(hyperedges, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmrmax") + en1 = model.addVars(hyperedges, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmr1") + en2 = model.addVars(hyperedges, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmr2") + en3 = model.addVars(hyperedges, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 1.0, name = "edgenmr3") + + for edge, elmax, el1, el2, el3 in zip(hyperedges, edgelikelihoodsmax, edgelikelihoods1, edgelikelihoods2, edgelikelihoods3): + enmax[edge] = elmax + en1[edge] = el1 + en2[edge] = el2 + en3[edge] = el3 - - print(enmax[2]) - print(enmax[3]) - print(enmax[17]) - print(enmax[41]) 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, (heads, tails) in hyperedges.items(): + count[mol] += x[edge] * (heads.count(mol) - tails.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()) @@ -195,6 +207,7 @@ def build_model(name, hyperedges, vertices, nmrlikelihoodsmax, nmrlikelihoods1, model.ModelSense = GRB.MAXIMIZE #Multiply edgelikelihood with the edge use boolean + #Adapt to have model.setObjectiveN( quicksum(enmax[e_id] * b[e_id] for e_id, (_, _) in hyperedges.items()), index = 0, @@ -219,13 +232,7 @@ def build_model(name, hyperedges, vertices, nmrlikelihoodsmax, nmrlikelihoods1, model.addConstr(b[213] == 1) model.addConstr(b[41] == 1) - #Every item created has to be consumed: - ''' - for mol in vertices: - for edge, (heads, tails) in hyperedges.items(): - count[mol] += x[edge] * (heads.count(mol) - tails.count(mol)) - model.addConstr(count[mol] == 0) - ''' + #No cyclic reaction pairs: for e_id1, (heads1, tails1) in hyperedges.items(): @@ -254,7 +261,7 @@ def print_solution(title, flow_solution, binary_solution, hyperedges): print(f"Number of used hyperedges: {len(binary_solution)}") def main(): - model, x, b = build_model("HypergraphFlow", HYPERGRAPH, VERTICES, NMRMAX, NMR1, NMR2, NMR3) + model, x, b = build_model("HypergraphFlow", HYPERGRAPH, EDGEMAX, EDGE1, EDGE2, EDGE3) model.optimize() if model.status != GRB.Status.OPTIMAL: print("No optimal solution found for the first model.") @@ -264,7 +271,7 @@ def main(): print_solution("Optimal Solution", optimal_solution, optimal_binary_solution, HYPERGRAPH) excluded_support = list(optimal_binary_solution.keys()) - second_model, x2, b2 = build_model("SecondBestHypergraphFlow", HYPERGRAPH, VERTICES, NMRMAX, NMR1, NMR2, NMR3, excluded_support=excluded_support,) + second_model, x2, b2 = build_model("SecondBestHypergraphFlow", HYPERGRAPH, EDGEMAX, EDGE1, EDGE2, EDGE3, excluded_support=excluded_support,) second_model.optimize() if second_model.status == GRB.Status.OPTIMAL: