kilian/nmrproject
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Dateien nach "ILP/Kaffee" hochladen

Browse Source
This commit is contained in:
kilian
2026-08-12 12:35:44 +02:00
parent 9950985030
commit 2a54d96d85
2 changed files with 446 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
ILP/Kaffee/caffeinesynthesis.py
+127
View File
@@ -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()
ILP/Kaffee/nmrSimilarityCaffeine.py
+319
View File
@@ -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()