nmrproject/ILP/nmrSimilarity.py
kilian 20824b6400 Dateien nach "ILP" hochladen
Now ILP solver picks optimal combination of two different similarities and similarity can also be calcuated using the mean of similarity instead of the comparitive value
2026-06-24 11:29:14 +02:00

353 lines
10 KiB
Python
Raw 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.

#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))
lowest_ppm = math.floor(min(shifts))
#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 main():
spectrumref = CCAFFEINE2
spectra = [CXANTHINE, C1XANTHINE, C3XANTHINE, C7XANTHINE, C1XANTHINE, C17XANTHINE, C37XANTHINE, C137XANTHINE]
spectranames = ["CXANTHINE", "C1XANTHINE", "C3XANTHINE", "C7XANTHINE", "C1XANTHINE", "C17XANTHINE", "C37XANTHINE", "C137XANTHINE"]
likelihood = []
for spectrumtrue in spectra:
#errorlist = {}
#errorlist = []
similaritybycorrection = []
correctionvalues = range(8, 16) #also tested 8.4, 8.37, 11, 9.4 (for CNMR3)
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 = []
for i in np.arange(0.1, 3.9, 0.1):
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]}')
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()