Dateien nach "ILP/butadien" hochladen

2026-08-28 12:10:40 +02:00
parent bded95d692
commit 7db935744c
2 changed files with 293 additions and 3 deletions
@@ -2,6 +2,7 @@ import gurobipy as gp
 from gurobipy import GRB, Model, quicksum

 HYPEREDGES = {
+    1: ([], ['Butadien']),
    4: (['Butadien', 'Butadien'], []), 
    5: (['Butadien', 'Butadien'], ['Butadien', 'Butadien']), 
    7: (['Butadien', 'Butadien'], []), 
@@ -180,7 +181,12 @@ HYPEREDGES = {
 }

 VERTICES = ['Butadien', 'p_{0,0}', 'p_{0,1}', 'p_{0,2}', 'p_{0,3}', 'p_{0,4}', 'p_{0,5}', 'p_{0,6}', 'p_{0,7}', 'p_{0,8}', 'p_{0,9}', 'p_{0,10}', 'p_{0,11}', 'p_{0,12}', 'p_{0,13}', 'p_{0,14}', 'p_{0,15}', 'p_{0,16}', 'p_{0,17}', 'p_{0,18}']
-
+#Vergleich mit dem NMR von Ethylen und Hexatrien
+NMR1 = [0.32, 0.5, 0.87, 0.0, 0.11, 0.58, 0.06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.09, 0.0, 0.0, 0.17, 0.0, 0.09]
+#Vergleich mit dem NMR von Ethylen und Octrien
+NMR2 = [0.33, 0.43, 0.6, 0.0, 0.13, 0.9, 0.07, 0.01, 0.0, 0.0, 0.01, 0.0, 0.0, 0.01, 0.1, 0.0, 0.0, 0.22, 0.0, 0.13]
+#Vergleich mit dem NMR von Ethylen und Benzol
+NMR3 = [0.0, 0.32, 0.0, 0.0, 0.0, 0.0, 0.05, 0.0, 0.0, 0.02, 0.0, 0.29, 0.95, 0.17, 0.0, 0.0, 0.0, 0.0, 0.0, 0.05]
 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']

 FIXED_FLOWS = {
@@ -192,8 +198,9 @@ def build_model(name, hyperedges, vertices, nmrlikelihoods, excluded_support=Non

    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")
-    
+    n1 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 2.0, name = "nmr1")
+    n2 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 2.0, name = "nmr2")
+    n3 = model.addVars(vertices, vtype=GRB.CONTINUOUS, lb = 0.0, ub = 2.0, name = "nmr3")

    for v, nmr in zip(vertices, nmrlikelihoods):
        n[v] = nmr
@@ -0,0 +1,283 @@
+import math
+import numpy as np
+
+#C=CC=C Butadien
+CBUTADIEN = {
+    1: ([122.26], [2]),
+    2: ([145.21], [2]),
+}
+
+#C=C Ethylen
+CP0 = {
+    1: ([128.33], [2]),
+}
+
+#C=C/C=C/C=C Hexatrien
+CP1 = {
+    1: ([122.47], [2]),
+    2: ([144.32], [2]),
+    3: ([140.50], [2]),
+}
+
+#C1CCC(C=C)CC=1
+CP2 = {
+    1: ([135.43], [1]),
+    2: ([31.80], [1]),
+    3: ([29.20], [1]),
+    4: ([42.60], [1]),
+    5: ([152.20], [1]),
+    6: ([114.57], [1]),
+    7: ([38.56], [1]),
+    8: ([134.81], [1]),
+}
+
+#C(CCC(C=C)CC=C)=C
+CP3 = {
+    1: ([147.47], [1]),
+    2: ([37.13], [1]),
+    3: ([33.61], [1]),
+    4: ([52.35], [1]),
+    5: ([151.07], [1]),
+    6: ([119.60], [1]),
+    7: ([49.64], [1]),
+    8: ([143.68], [1]),
+    9: ([121.81], [1]),
+    10: ([115.67], [1]),
+}
+
+#C=CC=CC=CC=C Octatetraen
+CP4 = {
+    1: ([122.10], [2]),
+    2: ([144.40], [2]),
+    3: ([140.76], [2]),
+    4: ([139.87], [2]),
+}
+
+#C1CCC(C=CC=C)CC=1
+CP5 = {
+    1: ([135.39], [1]),
+    2: ([31.78], [1]),
+    3: ([30.01], [1]),
+    4: ([42.13], [1]),
+    5: ([147.99], [1]),
+    6: ([133.91], [1]),
+    7: ([144.86], [1]),
+    8: ([118.59], [1]),
+    9: ([38.62], [1]),
+    10: ([134.73], [1]),
+}
+
+#C=CC1C=CCCC1
+CP6 = {
+    1: ([117.64], [1]),
+    2: ([152.55], [1]),
+    3: ([50.21], [1]),
+    4: ([138.81], [1]),
+    5: ([135.87], [1]),
+    6: ([31.05], [1]),
+    7: ([27.64], [1]),
+    8: ([35.53], [1]),
+}
+
+#C1CCCCC=1
+CP7 = {
+    1: ([135.79], [2]),
+    2: ([31.39], [2]),
+    3: ([28.19], [2]),
+}
+
+#C=CC1CC=CCC1C=C
+CP8 = {
+    1: ([119.30], [2]),
+    2: ([151.13], [2]),
+    3: ([46.76], [2]),
+    4: ([39.99], [2]),
+    5: ([134.37], [2]),
+}
+
+#C1CCC(C=C)C(C=C)C=1
+CP9 = {
+    1: ([135.24], [1]),
+    2: ([31.15], [1]),
+    3: ([34.06], [1]),
+    4: ([52.33], [1]),
+    5: ([151.22], [1]),
+    6: ([118.67], [1]),
+    7: ([54.94], [1]),
+    8: ([150.28], [1]),
+    9: ([119.75], [1]),
+    10: ([138.06], [1]),
+}
+
+#C(C1CC(C=C)C=CC1)=C
+CP10 = {
+    1: ([152.90], [1]),
+    2: ([47.12], [1]),
+    3: ([36.63], [1]),
+    4: ([46.21], [1]),
+    5: ([150.05], [1]),
+    6: ([115.14], [1]),
+    7: ([137.57], [1]),
+    8: ([132.93], [1]),
+    9: ([38.14], [1]),
+    10: ([117.15], [1]),
+}
+
+#C1C=CC=CC=1 Benzol
+CP11 = {
+    1: ([132.96], [6])
+}
+
+#C(CCC1C=CC=CC1)=C
+CP12 = {
+    1: ([148.24], [1]),
+    2: ([39.58], [1]),
+    3: ([40.28], [1]),
+    4: ([39.83], [1]),
+    5: ([138.15], [1]),
+    6: ([131.12], [1]),
+    7: ([131.68], [1]),
+    8: ([133.29], [1]),
+    9: ([35.91], [1]),
+    10: ([118.42], [1]),
+}
+
+#C1C(C=CC=C)CCCC=1
+CP13 = {
+    1: ([137.84], [1]),
+    2: ([45.55], [1]),
+    3: ([147.44], [1]),
+    4: ([135.54], [1]),
+    5: ([144.68], [1]),
+    6: ([118.72], [1]),
+    7: ([36.87], [1]),
+    8: ([28.62], [1]),
+    9: ([31.29], [1]),
+    10: ([137.16], [1]),
+}
+
+#C=CC(C=C)CCCC=C
+CP14 = {
+    1: ([118.15], [2]),
+    2: ([149.96], [2]),
+    3: ([51.02], [1]),
+    4: ([29.41], [1]),
+    5: ([32.52], [1]),
+    6: ([36.00], [1]),
+    7: ([146.13], [1]),
+    8: ([119.77], [1]),
+}
+
+#C=CCCCCC=C
+CP15 = {
+    1: ([118.87], [2]),
+    2: ([147.62], [2]),
+    3: ([42.24], [2]),
+    4: ([36.95], [2]),
+}
+
+#C=CC1C=CC(C=C)CC1
+CP16 = {
+    1: ([116.08], [2]),
+    2: ([150.56], [2]),
+    3: ([44.33], [2]),
+    4: ([140.69], [2]),
+    5: ([26.35], [2]),
+}
+
+#C1CC2CCCCC2CC=1
+CP17 = {
+    1: ([135.03], [2]),
+    2: ([37.74], [2]),
+    3: ([34.32], [2]),
+    4: ([30.76], [2]),
+    5: ([24.84], [2]),
+}
+
+#C1CC2C=CC=CC2CC=1
+CP18 = {
+    1: ([139.14], [2]),
+    2: ([36.29], [2]),
+    3: ([40.97], [2]),
+    4: ([137.28], [2]),
+    5: ([127.67], [2]),
+}
+
+
+
+
+
+
+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 addspectra(spectrum1, spectrum2):
+    spectrum = spectrum1.copy()
+    for _, ([ppm2], [height2]) in spectrum2.items():
+        for peak, ([ppm1], [height1]) in spectrum1.items():
+            if ppm1 == ppm2:
+                spectrum[peak] == ([ppm1], [height1 + height2])
+                continue
+        spectrum[len(spectrum) + 1] = ([ppm2], [height2])
+    return(spectrum)
+
+
+def main():
+    spectrumrefs = [addspectra(CP0, CP1), addspectra(CP0, CP4), addspectra(CP0, CP11)]
+    #1H-NMR Spectra ignoriert, da meiste H sauer, da an N gebunden
+    #spectra = [HCINNAMICACID, HPCOUMARICACID, HMCOUMARICACID, HBENZALDEHYD, HCAFFEICACID, H3HYDROXYBENZALDEHYD, H4HYDROXYBENZALDEHYD, H34DIHYDROXYBENZALDEHYD]
+    spectra = [CBUTADIEN, CP0, CP1, CP2, CP3, CP4, CP5, CP6, CP7, CP8, CP9, CP10, CP11, CP12, CP13, CP14, CP15, CP16, CP17, CP18]
+    spectranames = ['C=CC=C or Butadien', 'C=C or Ethyen or CP0', 'C=CC=CC=C or Hexatrien or CP1', 'C1CCC(C=C)CC=1 or CP2', 'C(CCC(C=C)CC=C)=C or CP3', 'C=CC=CC=CC=C or Octatetraen or CP4', 'C1CCC(C=CC=C)CC=1 or CP5', 'C=CC1C=CCCC1 or CP6', 'C1CCCCC=1 or Cyclohexen or CP7', 'C=CC1CC=CCC1C=C or CP8', 'C1CCC(C=C)C(C=C)C=1 or CP9', 'C(C1CC(C=C)C=CC1)=C or CP10', 'C1C=CC=CC=1 or CP11', 'C(CCC1C=CC=CC1)=C or CP12', 'C1C(C=CC=C)CCCC=1 or CP13', 'C=CC(C=C)CCCC=C or CP14', 'C=CCCCCC=C or CP15', 'C=CC1C=CC(C=C)CC1 or CP16', 'C1CC2CCCCC2CC=1 or CP17', 'C1CC2C=CC=CC2CC=1 or CP18']
+    for spectrumref in spectrumrefs:
+        likelihood = []
+        for spectrumtrue in spectra:
+            similaritylist = []
+            binwidthlist = np.arange(0.1, 1.1, 0.1)
+            for i in binwidthlist:
+                similaritylist.append(similarity_nmr(spectrumtrue, spectrumref, i))
+            similaritymean = round(sum(similaritylist) / len(similaritylist), 2)
+            likelihood.append(similaritymean)
+        print(likelihood)
+if __name__ == "__main__":
+    main()