Problem with line fit

[n-salvi]

Hi all,

I am getting the following error when trying to fit a spectrum:

Traceback (most recent call last): File "/Users/nicola/GoogleDrive/Python/NMR/peak_picking.py", line 48, in ng.analysis.linesh.fit_spectrum(data, lineshapes, params, amps, bounds, ampbounds, centers, rIDs, box_width, error_flag=True) File "/sw/lib/python2.7/site-packages/nmrglue/analysis/linesh.py", line 277, in fit_spectrum campbounds, wmask, error_flag, _kw) File "/sw/lib/python2.7/site-packages/nmrglue/analysis/linesh.py", line 521, in fit_NDregion r = f_NDregion(region, ls_classes, p0, p_bounds, n_peaks, wmask, _kw) File "/sw/lib/python2.7/site-packages/nmrglue/analysis/linesh.py", line 807, in f_NDregion p_best = leastsqbound(err_NDregion, p0, bounds=p_bounds, args=args, *_kw) File "/sw/lib/python2.7/site-packages/nmrglue/analysis/leastsqbound.py", line 253, in leastsqbound m = _check_func('leastsq', 'func', func, x0, args, n)[0] File "/sw/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 20, in _check_func res = atleast1d(thefunc(((x0[:numinputs],) + args))) File "/sw/lib/python2.7/site-packages/nmrglue/analysis/linesh.py", line 771, in err_NDregion sim_region = s_NDregion(list(p), shape, ls_classes, n_peaks) File "/sw/lib/python2.7/site-packages/nmrglue/analysis/linesh.py", line 725, in s_NDregion r = s_single_NDregion([A] + curr_p, shape, ls_classes) File "/sw/lib/python2.7/site-packages/nmrglue/analysis/linesh.py", line 763, in s_single_NDregion r = np.kron(r, ls) # vector direct product flattened File "/sw/lib/python2.7/site-packages/numpy/lib/shape_base.py", line 782, in kron result = concatenate(result, axis=axis) ValueError: need at least one array to concatenate

Here is the script I am using:

import nmrglue as ng import numpy as np import sys

dic, data = ng.pipe.read(sys.argv[1]) noise = data.std()

peaks = ng.analysis.peakpick.pick(data, pthres=3*noise, est_params=True, cluster=True)

lineshapes = ['pv', 'pv']

params_x = [(x, lw, 0.5) for x, lw in zip(peaks['X_AXIS'], peaks['X_LW'])] params_y = [(y, lw, 0.5) for y, lw in zip(peaks['Y_AXIS'], peaks['Y_LW'])] params = [[par_x, par_y] for par_x, par_y in zip(params_x, params_y)]

amps = peaks['VOL'] bounds = [[((None, None), (0.0, None), (0.0, 1.0)), ((None, None), (0.0, None), (0.0, 1.0))] for peak in peaks] ampbounds = [(None, None) for peak in peaks] centers = [(x, y) for x, y in zip(peaks['X_AXIS'], peaks['Y_AXIS'])] rIDs = peaks['cID']

box_width = (np.mean(peaks['X_LW']), np.mean(peaks['Y_LW']))

ng.analysis.linesh.fit_spectrum(data, lineshapes, params, amps, bounds, ampbounds, centers, rIDs, box_width, error_flag=True)

Am I missing something?

Thanks, Nicola

[jjhelmus]

Withing running this with actually data I have to guess at the cause. My bet is that one of the regions being fit has one or more dimensions which have zero size, check that the box_width is not zero in either dimension and that the are not any peaks very near the edge of the spectrum.

[jjhelmus]

@n-salvi Can you provide a copy of the data you are working on or at least a description of the number and sizes of the dimensions.

[n-salvi]

Thank you @jjhelmus for your prompt reply. Would a copy of one processed spectrum work?

https://dl.dropboxusercontent.com/u/23799377/test.ft2

By the way, would your recommend any procedure for fitting lineshapes in pseudo-3d experiments, e.g. for measuring relaxation? Ideally, it should be possible to vary only the signal intensity among the planes while keeping the same peak center and line widths for all of them.

Thanks again, Nicola

[jjhelmus]

You need to swap the order of the axes when setting up the parameters for ng.linesh.fit_spectrum. This still causes some issue when estimating the errors but the following works with no changes:

import nmrglue as ng
import numpy as np
import sys

dic, data = ng.pipe.read('test.ft2')
noise = data.std()

peaks = ng.analysis.peakpick.pick(
    data, pthres=3*noise, est_params=True, cluster=True)

lineshapes = ['pv', 'pv']

params_x = [(x, lw, 0.5) for x, lw in zip(peaks['X_AXIS'], peaks['X_LW'])]
params_y = [(y, lw, 0.5) for y, lw in zip(peaks['Y_AXIS'], peaks['Y_LW'])]
params = [[par_y, par_x] for par_x, par_y in zip(params_x, params_y)]

amps = peaks['VOL']
bounds = [[
    ((None, None), (0.0, None), (0.0, 1.0)),
    ((None, None), (0.0, None), (0.0, 1.0))] for peak in peaks]
ampbounds = [(None, None) for peak in peaks]
centers = [(x, y) for x, y in zip(peaks['Y_AXIS'], peaks['X_AXIS'])]
rIDs = peaks['cID']
box_width = (np.mean(peaks['Y_LW']), np.mean(peaks['X_LW']))
print box_width

ng.analysis.linesh.fit_spectrum(
    data, lineshapes, params, amps, bounds, ampbounds, centers, rIDs,
    box_width, error_flag=False)

The issue stems from the fact that the peakpick treats the leading dimension as Y and the last dimension as X, so the parameters and more importantly the box centers must be specified as Y, X. Switching these two results in empty regions being selected out and then trying to fit peaks to an empty array. A better error should probably be reported when this is detected.

Even with the swap if error_flag is True there will still be issues as some of the fitting fails to converge. I'll be adding a fix so that None is returned for the estimated parameter errors when this occurs in a few moments.

[jjhelmus]

As for fitting pseudo-3D experiments with nmrglue, the scale lineshape can be used to model the pseudo/relaxation dimension. This will fit peak intensities along with the peak location and linewidth in all spectra. This lineshape can be specified with the 's' or 'scale' abbreviation.

In my own work I did not find that this technique gave good results. I found better results by summing the intensities of all the point in a box around each peak in the spectrum and repeating this for each spectra in the experiment. Then the summations for each peak could be treated as trajectory for further analysis (for example fitting to an decaying exponential in the case of relaxation experiments). The fitting t1 data example provides code which can be used for this. When peaks have significant overlap this technique has problems.

In addition the section of the nmrglue paper on the "Analysis of NMR relaxation data" provides details on this method, the example presented in the paper are also available in the Relaxation trajectory analysis example.

[jjhelmus]

@n-salvi Is this issue still causing you problems or can I close this issue?