Dateien nach "ILP/Vanilla" hochladen

ILP/Vanilla/nmrSimilarityVanillien.py

@@ -0,0 +1,358 @@
+#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 = CNMR2
+    #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:
+        #errorlist = {}
+        #errorlist = []
+        similaritybycorrection = []
+        #Paper Chemical reviews Carbons bound to Heavy atoms (TMS) to high -> this could be reason for to high values.
+        correctionvalues =  [9.37] #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
+            #Likelihood by number of higher similarity than all others.
+            '''error = 0
+            total = 0
+            for spectrumfalse in spectra:
+                positive = 0
+                negative = 0
+                bad_binwidth = []
+                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)
+                    #print(falsesimilarity)
+                    if(truesimilarity - falsesimilarity < 0 or truesimilarity == 0):
+                        negative += 1
+                        bad_binwidth.append(i)
+                    else:
+                        positive += 1
+                    total += 1
+                #print(f'Wrong similarity result: {negative} and Right similarity result: {positive}')
+                #print(bad_binwidth)
+                error += negative
+            #errorlist[correctionvalue] = error
+            errorlist.append(error)
+        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 = []
+            binwidthlist = np.arange(0.1, 3.9, 0.1) #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)
+        name = spectranames[spectra.index(spectrumtrue)]
+        correctionindex = max(range(len(similaritybycorrection)), key=similaritybycorrection.__getitem__)
+        print(f'{name}: {correctionindex} = {correctionvalues[correctionindex]}')
+        #Maybe not the best but a mean instead?
+        likelihood.append(round(sum(similaritybycorrection)/len(similaritybycorrection), 2))
+        #likelihood.append(round(max(similaritybycorrection), 2))
+    print(likelihood)
+    '''for i in np.arange(0.01, 0.07, 0.01):
+        print(f'Increment i: {i}')
+        print(similarity_nmr(HNMR1, HNMR2, i))
+        print(similarity_nmr(H1XANTHINE, HNMR1, i))
+        print(similarity_nmr(H3XANTHINE, HNMR1, i))
+        print(similarity_nmr(H7XANTHINE, HNMR1, i))
+    '''
+if __name__ == "__main__":
+    main()

ILP/Vanilla/vanillasynthesis.py

@@ -0,0 +1,128 @@
+import gurobipy as gp
+from gurobipy import GRB, Model, quicksum
+
+HYPEREDGES = {
+    1: ([], ['Cinnamicacid']),
+    2: (['Cinnamicacid'], ['p_{0,0}']),
+    3: (['Cinnamicacid'], ['p_{0,1}']),
+    4: (['Cinnamicacid'], ['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}'], ['34Dihydroxybenzaldehyd']),
+    12: (['p_{0,4}'], ['34Dihydroxybenzaldehyd']),
+    13: (['p_{0,5}'], ['34Dihydroxybenzaldehyd']),
+    14: (['34Dihydroxybenzaldehyd'], []),
+}
+
+#Similyrity of Nodes NMR to Measured NMR
+#Can how to add Values along a Path?
+#Can you count position on path?
+#Can you save spectra values to add hypergraphposition to spectrainformation?
+#Compositespectra vs oredered set of spectra
+
+VERTICES = ['Cinnamicacid', 'p_{0,0}', 'p_{0,1}', 'p_{0,2}', 'p_{0,3}', 'p_{0,4}', 'p_{0,5}',  '34Dihydroxybenzaldehyd']
+
+#Change to only have one likelihood, where with manual order association
+NMRLIKELYHOODS = [0.0, 0.2, 0.7, 0.1, 0.1, 0.2, 0.7, 0.0]
+
+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()