jjhelmus / nmrglue

A module for working with NMR data in Python
BSD 3-Clause "New" or "Revised" License
209 stars 86 forks source link

Peak Fitting of Grossly Overlapped peaks in SSNMR #74

Open witheve17 opened 6 years ago

witheve17 commented 6 years ago

Hello,

I am trying to use the the ng.linesh functionality to fit peaks of SSNMR spectra that are sharp narrow Lorenztian with in a broad peak. the peak piking function does not result in a reasonable list of peaks because bot the narrow peak and the broad peak have the same center value. There because I know I want to fit to just two peaks I attempted to bypass the find peak step to specify the input in to the fit function. However it is not very clear from the example the correct formatting that would allow for this to happen. This is because the example is dependent on passing information from another embedded function.

I tried the defining the following for a simple two gauss fit

lineshapes=['g', 'g'] params=[[intens1, mean1,fwhm1],[intens2,mean2,fwhm2] amps=[amp1,amp2] (because all value must be positive) bounds=[[[(0,None),(0,None),(0,None)],[(0,None),(0,None),(0,None)]]] ampbounds=[amp1max,amp2max] centers=[center1, center2] error_flag = True ng.linesh.fit_NDregion=(data, lineshapes, params, amps, bounds, ampbounds, centers, error_flag)

I get the following error


ValueError Traceback (most recent call last)

in () 7 error_flag = True 8 ----> 9 params_best, amp_best, iers = ng.linesh.fit_NDregion(data, lineshapes, params, amps, bounds, ampbounds, centers, error_flag) /usr/local/lib/python2.7/site-packages/nmrglue/analysis/linesh.pyc in fit_NDregion(region, lineshapes, params, amps, bounds, ampbounds, wmask, error_flag, **kw) 423 if len(guess) != ndim: 424 err = "Incorrect number of params for peak %i" --> 425 raise ValueError(err % (i)) 426 427 for j, dim_guess in enumerate(guess): # dimension loop ValueError: Incorrect number of params for peak 0 Which I assume I am not comprehended the gaussian function employed by this fitting function. I looked at the documentation for the gaussian function being employed and it does just need an guess intensity, a mean, and a fwhm. I do not know how to proceed because the peaks are grossly overlapping. Any suggestions. I was also wondering if there were plans to run Simpson through python and not just convert the data. Best, Velencia
jjhelmus commented 6 years ago

@witheve17 The fit_NDregion function can be confusing. It is very particular about the layout of the parameters and the error messages are not very helpful. There looks to be a few issues with the code snippet your provided.

Perhaps the following example would be helpful. Here a 1D spectrum with two overlapping peaks is simulated with sim_NDregion and them fit using fit_NDregion. There are a number of comments in the code that try to explain the layout of the various arguments to fit_NDregion.

import nmrglue as ng

# simulate a 1D spectrum with overlapped peaks
sim_lineshapes = ['g']
sim_params = [[(50, 22.0)], [(62, 2.1)]]
sim_amps = [100, 150]
region = ng.linesh.sim_NDregion((128, ), sim_lineshapes, sim_params, sim_amps)

# plot the simulated spectrum
import matplotlib.pyplot as plt
plt.plot(region)
plt.show()

# one dimensional data, all peaks are fit to a Gaussian lineshape
lineshapes = ['g']

# two, one-dimensional peaks each with two parameters; peak center and linewidth
center1 = 49.5
lw1 = 21.0

center2 = 61.5
lw2 = 2.2

#            peak 1,           peak 2
params = [ [(center1, lw1)], [(center2, lw2)] ]

# estimated amplitudes for the two peaks
amps = [99.5, 149.5]

# bounds (min, max) for each parameter for each peak
# both the peak center and linewidth must be greater than 0
#         peak 1: center      linewidth  peak 2: center, linewidth
bounds = [      [((0, None), (0, None))],   [((0, None), (0, None))] ]

# amplitude bounds (min, max) for each of the two peaks
#            peak 1   , peak 2
ampbounds = [(0, None), (0, None)]
error_flag = False

params_best, amps_best, ier = ng.linesh.fit_NDregion(
    region, lineshapes, params, amps, bounds=bounds, ampbounds=ampbounds,
    error_flag=error_flag)
print(params_best)
print(amps_best)
print(ier)

results

[[[49.999999999999993, 21.999999999999996]], [[62.0, 2.0999999999999983]]]
[ 100.  150.]
2

figure_1