diff --git a/ILP/Kaffee/caffeinesynthesis.py b/ILP/Kaffee/caffeinesynthesis.py new file mode 100644 index 0000000..6019e0d --- /dev/null +++ b/ILP/Kaffee/caffeinesynthesis.py @@ -0,0 +1,127 @@ +import gurobipy as gp +from gurobipy import GRB, Model, quicksum + +HYPEREDGES = { + 1: ([], ['Xanthine']), + 2: (['Xanthine'], ['p_{0,0}']), + 3: (['Xanthine'], ['p_{0,1}']), + 4: (['Xanthine'], ['p_{0,2}']), + 5: (['p_{0,0}'], ['p_{0,3}']), + 6: (['p_{0,0}'], ['p_{0,4}']), + 7: (['p_{0,1}'], ['p_{0,3}']), + 8: (['p_{0,1}'], ['p_{0,5}']), + 9: (['p_{0,2}'], ['p_{0,4}']), + 10: (['p_{0,2}'], ['p_{0,5}']), + 11: (['p_{0,3}'], ['Caffeine']), + 12: (['p_{0,4}'], ['Caffeine']), + 13: (['p_{0,5}'], ['Caffeine']), + 14: (['Caffeine'], []), +} + +VERTICES = ['Xanthine', 'p_{0,0}', 'p_{0,1}', 'p_{0,2}', 'p_{0,3}', 'p_{0,4}', 'p_{0,5}', 'Caffeine'] + +#Change to only have one likelihood, where with manual order association +NMRLIKELYHOODS = [0.0, 0.52, 0.43, 0.62, 0.53, 0.56, 0.68, 0.0] +#Results for comparison with 7-Methylxanthine +NMRLIKELYHOODS1 = [0.54, 0.52, 0.43, 0.62, 0.52, 0.56, 0.48, 0.5] +#Results for comparison with 3,7-Methylxanthine +NMRLIKELYHOODS2 = [0.52, 0.52, 0.52, 0.61, 0.52, 0.56, 0.61, 0.59] +#Results for comparison with 3,7-Methylxanthine with Methyl C of spectrum 2 bei +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): + 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") + + + for v, nmr in zip(vertices, nmrlikelihoods): + n[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}") + + for e_id, value in FIXED_FLOWS.items(): + 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}") + + reaction_path = {} + + + if excluded_support: + 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.ModelSense = GRB.MAXIMIZE + + + #Multiply node value with infow or outflow + 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", + ) + + model.setObjectiveN( + quicksum(-1 * x[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): + 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): + print(f"\n{title}:") + for e_id in sorted(flow_solution): + flow = flow_solution[e_id] + tails, heads = hyperedges[e_id] + print(f"Hyperedge {e_id}: Flow = {flow}, Tails = {tails}, Heads = {heads}") + print("\nBinary Variables:") + for e_id in sorted(binary_solution): + print(f"Binary Variable b_{e_id} = {binary_solution[e_id]}") + + print(f"\nTotal flow: {sum(flow_solution.values())}") + print(f"Number of used hyperedges: {len(binary_solution)}") + +def main(): + model, x, b = build_model("HypergraphFlow", HYPEREDGES, VERTICES, NMRLIKELYHOODS) + model.optimize() + if model.status != GRB.Status.OPTIMAL: + print("No optimal solution found for the first model.") + return + optimal_solution = positive_entries(x) + optimal_binary_solution = positive_entries(b) + print_solution("Optimal Solution", optimal_solution, optimal_binary_solution, HYPEREDGES) + + """ excluded_support = list(optimal_binary_solution.keys()) + second_model, x2, b2 = build_model("SecondBestHypergraphFlow", HYPEREDGES, excluded_support=excluded_support,) + + second_model.optimize() + if second_model.status == GRB.Status.OPTIMAL: + second_solution = positive_entries(x2) + second_binary_solution = positive_entries(b2) + print_solution("Second Best Solution", second_solution, second_binary_solution, HYPEREDGES, VERTICES) + else: + print("No optimal solution found for the second best model.") + """ +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/ILP/Kaffee/nmrSimilarityCaffeine.py b/ILP/Kaffee/nmrSimilarityCaffeine.py new file mode 100644 index 0000000..3506d8f --- /dev/null +++ b/ILP/Kaffee/nmrSimilarityCaffeine.py @@ -0,0 +1,319 @@ +#Binning mostly for broader peaks? +#0.66 für H und 8.4 für C bei anderen TMS Werten +#Gute 13C Ergebnisse für alles +11 ppm: Im Vergeich mit Coffein haben alle disubstituierten bei wenigen Hohen Werten falsche Zuordnung, bei mono und nicht substituierten sogar keine Falsche zuornung (0.1 bis 5 mit 0.1 Schritten) +#+11 nicht universell, aber 9 bis 13 bei allen sweet spot + +import math +import numpy as np + +#Xanthine +HXANTHINE = { + 1: ([7.96], [1]), + 2: ([9.45], [1]), + 3: ([7.725], [1]), + 4: ([7.625], [1]), +} +CXANTHINE = { + 1: ([159.40], [1]), + 2: ([164.01], [1]), + 3: ([120.94], [1]), + 4: ([161.24], [1]), + 5: ([146.98], [1]), +} + +#1-Methylxanthine +H1XANTHINE = { + 1: ([7.93], [1]), + 2: ([9.45], [1]), + 3: ([4.05], [3]), + 4: ([7.91], [1]), +} +C1XANTHINE = { + 1: ([161.50], [1]), + 2: ([166.28], [1]), + 3: ([120.65], [1]), + 4: ([158.74], [1]), + 5: ([146.25], [1]), + 6: ([38.55], [1]), +} + +#3-Methylxanthine +H3XANTHINE = { + 1: ([4.15], [3]), + 2: ([7.73], [1]), + 3: ([7.99], [1]), + 4: ([9.49], [1]), +} +C3XANTHINE = { + 1: ([161.83], [1]), + 2: ([163.37], [1]), + 3: ([121.24], [1]), + 4: ([163.15], [1]), + 5: ([146.49], [1]), + 6: ([39.71], [1]), +} + +#7-Methylxanthine +H7XANTHINE = { + 1: ([7.55], [1]), + 2: ([4.47], [3]), + 3: ([7.72], [1]), + 4: ([7.655], [1]), +} +C7XANTHINE= { + 1: ([159.50], [1]), + 2: ([165.47], [1]), + 3: ([122.15], [1]), + 4: ([162.31], [1]), + 5: ([151.55], [1]), + 6: ([45.06], [1]), +} + +#Theophylline +H13XANTHINE = { + 1: ([4.03], [3]), + 2: ([7.98], [1]), + 3: ([4.19], [3]), + 4: ([9.49], [1]), +} +C13XANTHINE = { + 1: ([163.77], [1]), + 2: ([165.26], [1]), + 3: ([120.73], [1]), + 4: ([160.99], [1]), + 5: ([145.80], [1]), + 6: ([40.42], [1]), + 7: ([37.60], [1]), +} + +#Paraxanthine +H17XANTHINE = { + 1: ([4.50], [3]), + 2: ([7.70], [1]), + 3: ([3.98], [3]), + 4: ([7.82], [1]), +} +C17XANTHINE = { + 1: ([161.41], [1]), + 2: ([167.17], [1]), + 3: ([121.81], [1]), + 4: ([160.18], [1]), + 5: ([151.09], [1]), + 6: ([45.17], [1]), + 7: ([36.96], [1]), +} + +CPARAXANTHINE = { + 1: ([26.7], [1]), + 2: ([32.9], [1]), + 3: ([151.1], [1]), + 4: ([106.5], [1]), + 5: ([147.4], [1]), + 6: ([155.3], [1]), + 7: ([143.0], [1]), +} + + +#Theobromine +H37XANTHINE = { + 1: ([4.49], [3]), + 2: ([7.75], [1]), + 3: ([4.11], [3]), + 4: ([7.65], [1]), +} +C37XANTHINE = { + 1: ([161.76], [1]), + 2: ([164.86], [1]), + 3: ([122.51], [1]), + 4: ([164.29], [1]), + 5: ([151.10], [1]), + 6: ([39.33], [1]), + 7: ([45.07], [1]), +} + +#Caffeine +H137XANTHINE = { + 1: ([7.73], [1]), + 2: ([4.15], [3]), + 3: ([4.52], [3]), + 4: ([4.01], [3]), +} +C137XANTHINE = { + 1: ([163.66], [1]), + 2: ([166.66], [1]), + 3: ([122.03], [1]), + 4: ([162.23], [1]), + 5: ([150.50], [1]), + 6: ([40.09], [1]), + 7: ([45.23], [1]), + 8: ([37.17], [1]), +} + +CCAFFEINE = { + 1: ([155.7], [1]), #166 + 2: ([148.8], [1]), #159 + 3: ([107.7], [1]), #118 + 4: ([152.2], [1]), #163 + 5: ([143.0], [1]), #154 + 6: ([27.2], [1]), #38 + 7: ([29.1], [1]), #40 + 8: ([32.9], [1]), #44 +} + +CCAFFEINEADJUSTED = { + 1: ([166.7], [1]), #166 + 2: ([159.8], [1]), #159 + 3: ([118.7], [1]), #118 + 4: ([163.2], [1]), #163 + 5: ([154.0], [1]), #154 + 6: ([38.2], [1]), #38 + 7: ([40.1], [1]), #40 + 8: ([44.9], [1]), #44 +} +CCAFFEINE2 = { + 1: ([27.7], [1]), + 2: ([29.3], [1]), + 3: ([33.1], [1]), + 4: ([151.0], [1]), + 5: ([148.1], [1]), + 6: ([106.6], [1]), + 7: ([154.5], [1]), + 8: ([142.8], [1]), +} + +C137XANTHINEADJUSTED = { + 1: ([155.62], [1]), + 2: ([158.29], [1]), + 3: ([113.66], [1]), + 4: ([153.86], [1]), + 5: ([142.13], [1]), + 6: ([31.72], [1]), + 7: ([36.86], [1]), + 8: ([28.8], [1]), +} + +#Experimental 7-Methylxanthine nmr +#Secundary source 11.52, 3.81 +HNMR1= { + 1: ([10.85], [1]), + 2: ([11.50], [1]), + 3: ([3.82], [3]), + 4: ([7.88], [1]), +} + +CNMR1= { + 1: ([155.85], [1]), + 2: ([151.35], [1]), + 3: ([149.30], [1]), + 4: ([143.01], [1]), + 5: ([106.90], [1]), + 6: ([33.03], [1]), +} + + + +#Experimental Theobromine nmr +HNMR2 = { + 1: ([11.10], [1]), + 2: ([3.33], [3]), + 3: ([3.84], [3]), + 4: ([7.97], [1]), +} + +CNMR2 = { + 1: ([154.9], [1]), + 2: ([149.8], [1]), + 3: ([107.1], [1]), + 4: ([151.0], [1]), + 5: ([142.8], [1]), + 6: ([29.3], [1]), + 7: ([33.9], [1]), +} + +#Combination of methyl group and base purine rings from two papers +CNMR3 = { + 1: ([154.9], [1]), + 2: ([150.0], [1]), + 3: ([108.1], [1]), + 4: ([153.1], [1]), + 5: ([142.6], [1]), + 6: ([29.3], [1]), + 7: ([33.9], [1]), +} + +def overlap(listref, listnew): + twoleft = np.sum(np.multiply(np.concatenate((listref, [0, 0])), np.concatenate(([0, 0], listnew)))) + oneleft = np.sum(np.multiply(np.concatenate((listref, [0])), np.concatenate(([0], listnew)))) + neutral = np.sum(np.multiply(listref,listnew)) + oneright = np.sum(np.multiply(np.concatenate(([0], listref)), np.concatenate((listnew, [0])))) + tworight = np.sum(np.multiply(np.concatenate(([0, 0], listref)), np.concatenate((listnew, [0, 0])))) + overlap = (oneleft + oneright)* 0.5 + neutral + return overlap + +def bin_array(spectra, highest_ppm, lowest_ppm, bin_width): + binnumber = math.ceil((highest_ppm - lowest_ppm)/bin_width) + bin = [0] * binnumber + for peak in spectra: + (shift, height) = spectra[peak] + binindex = math.floor((shift[0] - lowest_ppm) / bin_width) + bin[binindex] += height[0] + normalizedbin = np.divide(bin, np.sum(bin)) + return normalizedbin + +def define_border_values(spectraref, spectranew, bin_width): + shifts = [] + for _,(shift,_) in spectraref.items(): + shifts.append(shift[0]) + for _,(shift,_) in spectranew.items(): + shifts.append(shift[0]) + highest_ppm = math.ceil(max(shifts)) + bin_width + lowest_ppm = math.floor(min(shifts)) - bin_width + #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): + #Maximize likelihood or minimize Deviation + #Values for two spectra and optimize largest for both different? + #Spectra in Nodes to allow maximize overlapp with both spectra or one spectra. + #5.4.2 Eliminating X–H signals from 1H NMR spectra + lowest_ppm, highest_ppm = define_border_values(spectraref, spectranew, bin_width) + binref = bin_array(spectraref, highest_ppm, lowest_ppm, bin_width) + binnew = bin_array(spectranew, highest_ppm, lowest_ppm, bin_width) + crosscorr = overlap(binref, binnew) + refselfcorr = overlap(binref, binref) + newselfcorr = overlap(binnew, binnew) + simidx = crosscorr / math.sqrt(refselfcorr * newselfcorr) + return(simidx) + +def correction(spectra, corretionppm): + newspectra = {} + for id, (shift, height) in spectra.items(): + shiftvalue = shift[0] + adjustedshift = shiftvalue + corretionppm + newspectra[id] = ([adjustedshift], height) + return newspectra + +def main(): + spectrumref = CNMR3 + #1H-NMR Spectra ignoriert, da meiste H sauer, da an N gebunden + #spectra = [HXANTHINE, H1XANTHINE, H3XANTHINE, H7XANTHINE, H1XANTHINE, H17XANTHINE, H37XANTHINE, H137XANTHINE] + spectra = [CXANTHINE, C1XANTHINE, C3XANTHINE, C7XANTHINE, C1XANTHINE, C17XANTHINE, C37XANTHINE, C137XANTHINE] + spectranames = ["XANTHINE", "1XANTHINE", "3XANTHINE", "7XANTHINE", "1XANTHINE", "17XANTHINE", "37XANTHINE", "137XANTHINE"] + likelihood = [] + for spectrumtrue in spectra: + similaritybycorrection = [] + #Paper Chemical reviews Carbons bound to Heavy atoms (TMS) to high -> this could be reason for too high values. + correctionvalues = [11] #np.arange(8.37, 9.88, 0.01) #for C tested np.arange(8.0, 16.1, 0.1) range(8, 12) 8.4, 8.37, 11, 9.4 (for CNMR3), 9.87 (true for all ref, 8.37 + 1.5 for the precision), 9.37 (good for first, ok for second, third because only 7 better/equal but for first much higher) np.arange(6.13, 9.86, 0.01) (only for first), for H 0.66 (not good), np.arange(0.4, 1.0, 0.01), 0.6 for first, second never first either 17 or caf higher np.arange(0.51, 0.82, 0.01) good measure + 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 + similaritylist = [] + binwidthlist = np.arange(0.1, 3.9, 0.1) + for i in binwidthlist: + similaritylist.append(similarity_nmr(spectrumtrue, spectrumrefcorrected, i)) + similaritymean = sum(similaritylist) / len(similaritylist) + similaritybycorrection.append(similaritymean) + likelihood.append(round(sum(similaritybycorrection)/len(similaritybycorrection), 2)) + print(likelihood) + +if __name__ == "__main__": + main() \ No newline at end of file