Files
nmrproject/ILP/butadien/nmrSimilarityButadien.py

295 lines
7.7 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
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 XH 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()