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 correction(spectra, corretionppm): newspectra = {} for id, (shift, height) in spectra.items(): shiftvalue = shift[0] adjustedshift = shiftvalue + corretionppm newspectra[id] = ([adjustedshift], height) return newspectra 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(): normalize = False spectrumrefs = [addspectra(CP0, CP1), addspectra(CP0, CP4), addspectra(CP0, CP11)] 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: likelihoods = [] 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 = sum(similaritylist) / len(similaritylist) likelihoods.append(similaritymean) if normalize: normalizedlikelihood = [round(likelihood/np.sum(likelihoods), 2) for likelihood in likelihoods] print(normalizedlikelihood) if not normalize: notnormalizedlikelihood = [round(likelihood, 2) for likelihood in likelihoods] print(notnormalizedlikelihood) if __name__ == "__main__": main()