diff --git a/ILP/caffeinesynthesis.py b/ILP/caffeinesynthesis.py index 38f2294..eefced8 100644 --- a/ILP/caffeinesynthesis.py +++ b/ILP/caffeinesynthesis.py @@ -33,6 +33,10 @@ VERTICES1 = { 8: ([0, 0], ['Caffeine']), } +#25014CX3 + +#03 + VERTICES2 = { 1: ([0], ['Xanthine']), 2: ([0], ['p_{0,0}']), @@ -47,69 +51,107 @@ VERTICES2 = { VERTICES = ['Xanthine', 'p_{0,0}', 'p_{0,1}', 'p_{0,2}', 'p_{0,3}', 'p_{0,4}', 'p_{0,5}', 'Caffeine'] NMRLIKELYHOODS = [0.0, 0.2, 0.7, 0.1, 0.1, 0.2, 0.7, 0.0] -NMRLIKELYHOODS1 = [0.75, 0.59, 0.36, 0.89, 0.59, 0.86, 0.78, 0.65] -NMRLIKELYHOODS2 = [0.56, 0.75, 0.45, 0.79, 0.75, 0.76, 0.84, 0.74] -NMRLIKELYHOODS3 = [0.54, 0.53, 0.7, 0.77, 0.53, 0.85, 0.94, 0.79] +NMRLIKELYHOODS1 = [0.54, 0.52, 0.43, 0.62, 0.52, 0.56, 0.48, 0.5] +NMRLIKELYHOODS2 = [0.52, 0.52, 0.52, 0.61, 0.52, 0.56, 0.61, 0.59] +NMRLIKELYHOODS3 = [0.5, 0.53, 0.59, 0.59, 0.53, 0.56, 0.68, 0.62] FIXED_FLOWS = { 1: 1, 14: 1, } -def build_model(name, hyperedges, vertices, nmrlikelihoods, excluded_support=None): +def build_model(name, hyperedges, vertices, nmrlikelihoods1, nmrlikelihoods2, nmrlikelihoods3, 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} - n = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb=0.0, ub=1.0, name="nmr") + 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 = 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") + c1 = model.addVars(vertices, vtype = GRB.BINARY, name = "n1_choice") + c2 = model.addVars(vertices, vtype = GRB.BINARY, name = "n2_choice") + c3 = model.addVars(vertices, vtype = GRB.BINARY, name = "n3_choice") + - for v, nmr in zip(vertices, nmrlikelihoods): - n[v] = nmr - print(f'Vertice: {v}, Similarity: {n[v]}') + for v, nmr in zip(vertices, nmrlikelihoods1): + n1[v] = nmr + #print(f'Vertice: {v}, Similarity: {n1[v]}') + for v, nmr in zip(vertices, nmrlikelihoods2): + n2[v] = nmr + + for v, nmr in zip(vertices, nmrlikelihoods3): + n3[v] = nmr vertices = set(v for tails, heads in hyperedges.values() for v in tails + heads) for v in vertices: inflow = quicksum(x[e_id] for e_id, (_, heads) in hyperedges.items() if v in heads) outflow = quicksum(x[e_id] for e_id, (tails, _) in hyperedges.items() if v in tails) - model.addConstr(inflow == outflow, name=f"flow_conservation_{v}") + 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}") + 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}") + 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}") reaction_path = {} - #for v_id in vertices: - #model.addConstr() if excluded_support: - model.addConstr(quicksum(b[e_id] for e_id in excluded_support) <= len(excluded_support) - 1, name="different_hyperedges",) + model.addConstr(quicksum(b[e_id] for e_id in excluded_support) <= len(excluded_support) - 1, name = "different_hyperedges",) #Multiplizier den node Wert mit infow + outflow - - #model.setObjective(quicksum(n[v_id] for v_id in vertices), GRB.MINIMIZE) #Maximize Similarity of Nodes model.ModelSense = GRB.MAXIMIZE + #Multiply with a list of binarys where only one value is one and the rest zero and not both the same position + ''' model.setObjectiveN( - quicksum(n[t_id[0]] * x[e_id] for e_id, (_, t_id) in hyperedges.items() if t_id != []), - index=0, - priority=2, - name="maximize_nmr_similarity", + quicksum(n1[v_id] * c1[v_id] for v_id in vertices if v_id != []) + + quicksum(n3[v_id] * c3[v_id] for v_id in vertices if v_id != []), + index = 0, + priority = 2, + name = "maximize_nmr_similarity", ) + ''' + #Multiply node value with infow or outflow + ''' + model.setObjectiveN( + quicksum(n1[t_id[0]] * c1[t_id[0]] * x[e_id] for e_id, (_, t_id) in hyperedges.items() if t_id != []) + +quicksum(n3[t_id[0]] * c3[t_id[0]] * x[e_id] for e_id, (_, t_id) in hyperedges.items() if t_id != []), + index = 0, + priority = 2, + name = "maximize_nmr_similarity", + ) + ''' + + #Hybrid of both above + model.setObjective( + quicksum(n1[t_id[0]] * c1[t_id[0]] * x[e_id] for e_id, (_, t_id) in hyperedges.items() if t_id != []) + +quicksum(n3[t_id[0]] * c3[t_id[0]] * x[e_id] for e_id, (_, t_id) in hyperedges.items() if t_id != []), + sense = GRB.MAXIMIZE, + ) + + + model.addConstr(gp.quicksum(c1[v_id] for v_id in vertices) == 1) + model.addConstr(gp.quicksum(c3[v_id] for v_id in vertices) == 1) + for v_id in vertices: + model.addConstr(c1[v_id] + c3[v_id] <= 1) + + model.addConstr(gp.quicksum(x[e_id] for e_id, (_, t_id) in hyperedges.items()) == 5) + '''' model.setObjectiveN( quicksum(-1 *b[e_id] for e_id in hyperedges), index=1, priority=1, name="minimize_used_hyperedges", ) + ''' return model, x, b -def positive_entries(variable_dict, threshold=0.5): +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 print_solution(title, flow_solution, binary_solution, hyperedges): @@ -126,7 +168,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", HYPEREDGES, VERTICES, NMRLIKELYHOODS) + model, x, b = build_model("HypergraphFlow", HYPEREDGES, VERTICES, NMRLIKELYHOODS1, NMRLIKELYHOODS2, NMRLIKELYHOODS3) model.optimize() if model.status != GRB.Status.OPTIMAL: print("No optimal solution found for the first model.") diff --git a/ILP/nmrSimilarity.py b/ILP/nmrSimilarity.py index d6873ef..8a52f5e 100644 --- a/ILP/nmrSimilarity.py +++ b/ILP/nmrSimilarity.py @@ -232,10 +232,10 @@ CNMR2 = { #Combination of methyl group and base purine rings from two papers CNMR3 = { - 1: ([153.1], [1]), + 1: ([154.9], [1]), 2: ([150.0], [1]), 3: ([108.1], [1]), - 4: ([154.9], [1]), + 4: ([153.1], [1]), 5: ([142.6], [1]), 6: ([29.3], [1]), 7: ([33.9], [1]), @@ -268,7 +268,7 @@ def define_border_values(spectraref, spectranew, bin_width): shifts.append(shift[0]) highest_ppm = math.ceil(max(shifts)) lowest_ppm = math.floor(min(shifts)) - lowest_ppm = min(shifts) - bin_width/2 #Worse result. None of the previously wrong (except 0.6) become right + #lowest_ppm = min(shifts) - bin_width/2 #Worse result. None of the previously wrong (except 0.6) become right return (lowest_ppm, highest_ppm) def similarity_nmr(spectraref, spectranew, bin_width): @@ -294,21 +294,25 @@ def correction(spectra, corretionppm): return newspectra def main(): - spectrumref = CNMR3 - spectrafalse = [CXANTHINE, C1XANTHINE, C3XANTHINE, C7XANTHINE, C1XANTHINE, C17XANTHINE, C37XANTHINE, C137XANTHINE] + spectrumref = CCAFFEINE2 + spectra = [CXANTHINE, C1XANTHINE, C3XANTHINE, C7XANTHINE, C1XANTHINE, C17XANTHINE, C37XANTHINE, C137XANTHINE] + spectranames = ["CXANTHINE", "C1XANTHINE", "C3XANTHINE", "C7XANTHINE", "C1XANTHINE", "C17XANTHINE", "C37XANTHINE", "C137XANTHINE"] likelihood = [] - for spectrumtrue in spectrafalse: + for spectrumtrue in spectra: #errorlist = {} - errorlist = [] - for correctionvalue in range(8, 12): #also tested 8.4, 8.37, 11 + #errorlist = [] + similaritybycorrection = [] + correctionvalues = range(8, 16) #also tested 8.4, 8.37, 11, 9.4 (for CNMR3) + for correctionvalue in correctionvalues: spectrumrefcorrected = correction(spectrumref, correctionvalue) #CCAFFEINE 11 (klappt hier sehr gut) CCAFFEINE2 12 CPARAXANTHINE 10 CNMR1 9, 10 o 11 (sehr gut) CNMR2 10 o 11 - error = 0 + #Likelihood by number of higher similarity than all others. + '''error = 0 total = 0 - for spectrumfalse in spectrafalse: + for spectrumfalse in spectra: positive = 0 negative = 0 bad_binwidth = [] - for i in np.arange(0.1, 3.9, 0.1): #Irgendwie hat 0.99999999999 einen Fehler den 1.0 nicht hat + for i in np.arange(0.1, 2.6, 0.1): #successfull at max 3.9, but max 1.7 is lowest where nmr3 correctly classified, 1.6 increases likelihood of 3,7 over 1,7 even with 8.4 correction truesimilarity = similarity_nmr(spectrumtrue, spectrumrefcorrected, i) falsesimilarity = similarity_nmr(spectrumfalse, spectrumrefcorrected, i) #print(truesimilarity) @@ -324,8 +328,19 @@ def main(): error += negative #errorlist[correctionvalue] = error errorlist.append(error) - #print(min(range(len(errorlist)), key=errorlist.__getitem__)) - likelihood.append(round((total - min(errorlist))/total, 2)) + print(min(range(len(errorlist)), key=errorlist.__getitem__)) + likelihood.append(round((total - min(errorlist))/total, 2))''' + #Likelihood by mean similarity + #This method demonstrates the same problems as the other likelihood method + similaritylist = [] + for i in np.arange(0.1, 3.9, 0.1): + similaritylist.append(similarity_nmr(spectrumtrue, spectrumrefcorrected, i)) + similaritymean = sum(similaritylist) / len(similaritylist) + similaritybycorrection.append(similaritymean) + name = spectranames[spectra.index(spectrumtrue)] + correctionindex = max(range(len(similaritybycorrection)), key=similaritybycorrection.__getitem__) + print(f'{name}: {correctionindex} = {correctionvalues[correctionindex]}') + likelihood.append(round(max(similaritybycorrection), 2)) print(likelihood) '''for i in np.arange(0.01, 0.07, 0.01): print(f'Increment i: {i}')