From d90931e939a31de06dd869fad5f1b5ab9b5685a6 Mon Sep 17 00:00:00 2001 From: Juuso Lehtivarjo Date: Wed, 13 Nov 2024 09:36:38 +0200 Subject: [PATCH 01/38] Function rename for NumPy >2.0 --- nmrglue/process/proc_lp.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/nmrglue/process/proc_lp.py b/nmrglue/process/proc_lp.py index 4a8ef948..0269dc91 100644 --- a/nmrglue/process/proc_lp.py +++ b/nmrglue/process/proc_lp.py @@ -528,8 +528,8 @@ def cadzow_single(x, M, K, min_var=False): U, s, Vh = scipy.linalg.svd(X) # correct the singular values and truncate the rank K - Ul = np.mat(U[:, :K]) - Vlh = np.mat(Vh[:K, :]) # first K columns of V are first K rows of Vh + Ul = np.asmatrix(U[:, :K]) + Vlh = np.asmatrix(Vh[:K, :]) # first K columns of V are first K rows of Vh sl = s[:K] if min_var: # adjust singular values using minimum variance method @@ -537,7 +537,7 @@ def cadzow_single(x, M, K, min_var=False): s2 = (1. / (M - K)) * np.power(s[K:], 2).sum() sl = np.array([l - s2 / l for l in sl]) - Sl = np.mat(np.diag(sl)) + Sl = np.asmatrix(np.diag(sl)) # compute enhanced data vector for rank-reduced data matrix Xp = Ul * Sl * Vlh @@ -803,7 +803,7 @@ def find_lpc_svd(D, d): L = D.shape[0] m = D.shape[1] U, s, Vh = scipy.linalg.svd(D) # SVD decomposition - U, Vh = np.mat(U), np.mat(Vh) # make U and Vh matrices + U, Vh = np.asmatrix(U), np.asmatrix(Vh) # make U and Vh matrices Si = pinv_diagsvd(s, m, L) # construct the pseudo-inverse sigma matrix return np.array(Vh.H * Si * U.H * d) @@ -833,7 +833,7 @@ def find_lpc_qr(D, d): Find linear prediction filter using QR decomposition. """ q, r = scipy.linalg.qr(D) - q, r = np.mat(q), np.mat(r) + q, r = np.asmatrix(q), np.asmatrix(r) # SPEED # the next line is slow and the use of pinv2 should be avoided as @@ -852,9 +852,9 @@ def find_lpc_cholesky(D, d): # form the normal equation (D.H*D)*a = D.H*d # SPEED # this can be improved by using the Hankel nature of D - D = np.mat(D) - DhD = np.mat(np.dot(D.H, D)) - Dhd = np.mat(np.dot(D.H, d)) + D = np.asmatrix(D) + DhD = np.asmatrix(np.dot(D.H, D)) + Dhd = np.asmatrix(np.dot(D.H, d)) c, lower = scipy.linalg.cho_factor(DhD) # Compute Cholesky decomp. return scipy.linalg.cho_solve((c, lower), Dhd) # solve normal equation @@ -1018,7 +1018,7 @@ def find_lproots_hsvd(x, M, K, mode, zmethod='sm'): # SVD of data matrix and truncation of U to form Uk U, s, Vh = scipy.linalg.svd(X) - Uk = np.mat(U[:, :K]) # truncated U matrix of rank K + Uk = np.asmatrix(U[:, :K]) # truncated U matrix of rank K Ub = Uk[:-1] # Uk with bottom row removed Ut = Uk[1:] # Uk with top row removed From d18cfb7c2debac2a13d544d1ccfff98aba836dcd Mon Sep 17 00:00:00 2001 From: Juuso Lehtivarjo Date: Wed, 13 Nov 2024 10:18:23 +0200 Subject: [PATCH 02/38] More numpy 2.0 compatibility changes with pylint --- .../t1_measurements/extract_trajs.py | 4 +- .../t1_measurements/fit_exp_leastsq.py | 8 ++- examples/fitting_data/t1_measurements/pt.py | 5 +- .../integration/integrate_1d/integrate_1d.py | 4 +- .../integration/integrate_2d/integrate_2d.py | 6 +- .../extract_traj.py | 4 +- .../s12-s15_relaxation_analysis/plot_boxes.py | 16 ++--- .../plot_trajectories.py | 11 ++-- examples/plotting/plot_2d/plot_assignments.py | 4 +- examples/plotting/plot_2d/plot_boxes.py | 12 ++-- .../simulation/simulate_ucsf/make_ucsf.py | 2 +- nmrglue/fileio/pipe.py | 23 ++++--- nmrglue/fileio/table.py | 2 +- nmrglue/fileio/tecmag.py | 62 +++++++++---------- 14 files changed, 89 insertions(+), 74 deletions(-) diff --git a/examples/fitting_data/t1_measurements/extract_trajs.py b/examples/fitting_data/t1_measurements/extract_trajs.py index f96628b1..dbffd751 100755 --- a/examples/fitting_data/t1_measurements/extract_trajs.py +++ b/examples/fitting_data/t1_measurements/extract_trajs.py @@ -5,8 +5,8 @@ import numpy as np # read in the integration limits and list of spectra -peak_list = np.recfromtxt("boxes.in", encoding="UTF-8") -spectra_list = np.recfromtxt("spectra.in", encoding="UTF-8") +peak_list = np.genfromtxt("boxes.in", encoding="UTF-8") +spectra_list = np.genfromtxt("spectra.in", encoding="UTF-8") # prepare the trajs records array num_spec = spectra_list.size diff --git a/examples/fitting_data/t1_measurements/fit_exp_leastsq.py b/examples/fitting_data/t1_measurements/fit_exp_leastsq.py index eaa967a7..efeef81b 100755 --- a/examples/fitting_data/t1_measurements/fit_exp_leastsq.py +++ b/examples/fitting_data/t1_measurements/fit_exp_leastsq.py @@ -7,9 +7,11 @@ # read in the trajectories and times trajs = np.load("traj.npy") -t1 = np.recfromtxt("time.dat") +t1 = np.genfromtxt("time.dat") # fitting function and residual calculation + + def fit_func(p, x): A, R2 = p # bound A between 0.98 and 1.02 (although fits do not reflect this) @@ -20,11 +22,13 @@ def fit_func(p, x): return A * np.exp(-1.0 * np.array(x) * R2 / 1.0e6) + def residuals(p, y, x): err = y - fit_func(p, x) return err -p0 = [1.0, 0.05] # initial guess + +p0 = [1.0, 0.05] # initial guess fits = {} # loop over the peak trajectories diff --git a/examples/fitting_data/t1_measurements/pt.py b/examples/fitting_data/t1_measurements/pt.py index fa2710c6..3dec254c 100755 --- a/examples/fitting_data/t1_measurements/pt.py +++ b/examples/fitting_data/t1_measurements/pt.py @@ -6,6 +6,8 @@ import numpy as np # the same fit_func as in fit_exp_leastsq.py + + def fit_func(p, x): A, R2 = p @@ -18,10 +20,11 @@ def fit_func(p, x): return A * np.exp(-1.0 * np.array(x) * R2 / 1.0e6) + # read in the trajectories, fitting results, and times fits = pickle.load(open("fits.pickle", "rb")) trajs = np.load("traj.npy") -times = np.recfromtxt("time.dat") +times = np.genfromtxt("time.dat") sim_times = np.linspace(times[0], times[-1], 2000) diff --git a/examples/integration/integrate_1d/integrate_1d.py b/examples/integration/integrate_1d/integrate_1d.py index 4645f623..c47bb239 100755 --- a/examples/integration/integrate_1d/integrate_1d.py +++ b/examples/integration/integrate_1d/integrate_1d.py @@ -10,7 +10,7 @@ length = data.shape[0] # read in the integration limits -peak_list = np.recfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in") # determine the ppm scale uc = ng.pipe.make_uc(dic, data) @@ -40,7 +40,7 @@ ax.plot(peak_scale, peak.cumsum() / 100. + peak.max(), 'g-') ax.plot(peak_scale, [0] * len(peak_scale), 'r-') ax.text(peak_scale[0], 0.5 * peak.sum() / 100. + peak.max(), name, - fontsize=8) + fontsize=8) # write out the integration info tup = (name, peak_scale[0], peak_scale[-1], peak.sum()) diff --git a/examples/integration/integrate_2d/integrate_2d.py b/examples/integration/integrate_2d/integrate_2d.py index 46a6271f..13f81a88 100755 --- a/examples/integration/integrate_2d/integrate_2d.py +++ b/examples/integration/integrate_2d/integrate_2d.py @@ -8,10 +8,10 @@ dic, data = ng.pipe.read("nmrpipe_2d/test.ft2") # read in the integration limits -peak_list = np.recfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in") # prepare the output file -f = open("volumes.out",'w') +f = open("volumes.out", 'w') f.write("# Name\tVolume\n") # loop over the integration limits @@ -23,7 +23,7 @@ y0, y1 = y1, y0 vol = data[y0:y1 + 1, x0:x1 + 1].sum() - f.write("%s\t%.3f\n"%(name, vol)) + f.write("%s\t%.3f\n" % (name, vol)) # close the output file f.close() diff --git a/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py b/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py index 0d837103..65110227 100644 --- a/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py +++ b/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py @@ -2,8 +2,8 @@ import numpy as np # read the integration limits and list of spectra -peak_list = np.recfromtxt("boxes.in", names=True) -spectra_list = np.recfromtxt("spectra.in") +peak_list = np.genfromtxt("boxes.in", names=True) +spectra_list = np.genfromtxt("spectra.in") # create an array to hold the trajectories trajectories = np.empty((peak_list.size, spectra_list.size), dtype='float') diff --git a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py index a364b938..939b79d5 100644 --- a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py +++ b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py @@ -12,8 +12,8 @@ cl = 30000 * 1.20 ** np.arange(20) # read in the box limits and list of spectra -peak_list = np.recfromtxt("boxes.in", names=True) -spectra_list = np.recfromtxt("spectra.in") +peak_list = np.genfromtxt("boxes.in", names=True) +spectra_list = np.genfromtxt("spectra.in") # loop over the spectra for spec_number, spectra in enumerate(spectra_list): @@ -30,16 +30,16 @@ y0, y1 = y1, y0 # slice the data around the peak - slice = data[y0 - ypad:y1 + 1 + ypad, x0 - xpad:x1 + 1 + xpad] + slice_ = data[y0 - ypad:y1 + 1 + ypad, x0 - xpad:x1 + 1 + xpad] # create the figure fig = plt.figure() ax = fig.add_subplot(111) # plot the contours - print "Plotting:", peak, spec_number + print("Plotting:", peak, spec_number) extent = (x0 - xpad + 1, x1 + xpad - 1, y0 - ypad + 1, y1 + ypad - 1) - ax.contour(slice, cl, cmap=cmap, extent=extent) + ax.contour(slice_, cl, cmap=cmap, extent=extent) # draw a box around the peak ax.plot([x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0], 'k--') @@ -51,6 +51,6 @@ [y0 + 1, y0 + 1, y1 - 1, y1 - 1, y0 + 1], 'k--', alpha=0.35) # set the title, save the figure - ax.set_title('Peak: %s Spectrum: %i'%(peak, spec_number)) - fig.savefig('peak_%s_spectrum_%i'%(peak, spec_number)) - del(fig) + ax.set_title('Peak: %s Spectrum: %i' % (peak, spec_number)) + fig.savefig('peak_%s_spectrum_%i' % (peak, spec_number)) + del fig diff --git a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py index 7137c964..d67d2dc7 100644 --- a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py +++ b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py @@ -2,18 +2,21 @@ import matplotlib.pyplot as plt # exponential function used to fit the data -def fit_func(p,x): + + +def fit_func(p, x): A, R2 = p return A * np.exp(-1.0 * np.array(x) * R2 / 1.0e6) -fitting_results = np.recfromtxt('fits.txt') + +fitting_results = np.genfromtxt('fits.txt') experimental_relaxation_times = np.loadtxt("relaxation_times.in") -simulated_relaxation_times = np.linspace(0,4000000,2000) +simulated_relaxation_times = np.linspace(0, 4000000, 2000) # loop over the fitting results for peak, A, R2, ier in fitting_results: - print "Plotting:", peak + print("Plotting:", peak) # load the experimental and simulated relaxation trajectories experimental_trajectory = np.loadtxt(peak + '.dat') diff --git a/examples/plotting/plot_2d/plot_assignments.py b/examples/plotting/plot_2d/plot_assignments.py index 8165258c..9abf5ffa 100755 --- a/examples/plotting/plot_2d/plot_assignments.py +++ b/examples/plotting/plot_2d/plot_assignments.py @@ -20,7 +20,7 @@ dic, data = ng.pipe.read("nmrpipe_2d/test.ft2") # read in the integration limits -peak_list = np.recfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in") # create the figure fig = plt.figure() @@ -28,7 +28,7 @@ # plot the contours ax.contour(data, cl, cmap=cmap, - extent=(0, data.shape[1] - 1, 0, data.shape[0] - 1)) + extent=(0, data.shape[1] - 1, 0, data.shape[0] - 1)) # loop over the peaks for name, x0, y0, x1, y1 in peak_list: diff --git a/examples/plotting/plot_2d/plot_boxes.py b/examples/plotting/plot_2d/plot_boxes.py index 7c30606e..ebd4e6bf 100755 --- a/examples/plotting/plot_2d/plot_boxes.py +++ b/examples/plotting/plot_2d/plot_boxes.py @@ -21,7 +21,7 @@ dic, data = ng.pipe.read("nmrpipe_2d/test.ft2") # read in the integration limits -peak_list = np.recfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in") # loop over the peaks for name, x0, y0, x1, y1 in peak_list: @@ -32,16 +32,16 @@ y0, y1 = y1, y0 # slice the data around the peak - slice = data[y0 - ypad:y1 + 1 + ypad, x0 - xpad:x1 + 1 + xpad] + slice_ = data[y0 - ypad:y1 + 1 + ypad, x0 - xpad:x1 + 1 + xpad] # create the figure fig = plt.figure() ax = fig.add_subplot(111) # plot the contours - print "Plotting:", name + print("Plotting:", name) etup = (x0 - xpad + 1, x1 + xpad - 1, y0 - ypad + 1, y1 + ypad - 1) - ax.contour(slice, cl, cmap=cmap, extent=etup) + ax.contour(slice_, cl, cmap=cmap, extent=etup) # draw a box around the peak ax.plot([x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0], 'k--') @@ -50,11 +50,11 @@ ax.plot([x0 - 1, x1 + 1, x1 + 1, x0 - 1, x0 - 1], [y0 - 1, y0 - 1, y1 + 1, y1 + 1, y0 - 1], 'k--', alpha=0.35) ax.plot([x0 + 1, x1 - 1, x1 - 1, x0 + 1, x0 + 1], - [y0 + 1, y0 + 1, y1 - 1, y1 - 1, y0 + 1], 'k--',alpha=0.35) + [y0 + 1, y0 + 1, y1 - 1, y1 - 1, y0 + 1], 'k--', alpha=0.35) # set the title ax.set_title(name) # save the figure fig.savefig(f"{name}.png") - del(fig) + del fig diff --git a/examples/simulation/simulate_ucsf/make_ucsf.py b/examples/simulation/simulate_ucsf/make_ucsf.py index bc7b5126..892b8954 100644 --- a/examples/simulation/simulate_ucsf/make_ucsf.py +++ b/examples/simulation/simulate_ucsf/make_ucsf.py @@ -32,7 +32,7 @@ data = np.empty((512, 1024), dtype='float32') # read in the peak list -peak_list = np.recfromtxt('peaks.txt', names=True) +peak_list = np.genfromtxt('peaks.txt', names=True) npeaks = len(peak_list) # convert the peak list from PPM to points diff --git a/nmrglue/fileio/pipe.py b/nmrglue/fileio/pipe.py index f4f8ef60..9964c2f8 100644 --- a/nmrglue/fileio/pipe.py +++ b/nmrglue/fileio/pipe.py @@ -81,7 +81,7 @@ def read_table(filename): # print(dtd['names'],dtd['formats']) s = [l.encode('utf-8') for l in dl] - rec = np.recfromtxt(s, dtype=dtd, comments='XXXXXXXXXXX') + rec = np.genfromtxt(s, dtype=dtd, comments='XXXXXXXXXXX') return cl, pformat, np.atleast_1d(rec) @@ -192,6 +192,7 @@ def make_uc(dic, data, dim=-1): # dictionary/data creation # ############################ + fd2dphase_dic = {"magnitude": 0, "tppi": 1, "states": 2, "image": 3} @@ -318,7 +319,7 @@ def create_dic(udic, datetimeobj=datetime.datetime.now()): if ((dic["FDF1QUADFLAG"] == dic["FDF2QUADFLAG"] == dic["FDF3QUADFLAG"]) and (dic["FDF1QUADFLAG"] == dic["FDF4QUADFLAG"] == 1)): - dic["FDQUADFLAG"] = 1.0 + dic["FDQUADFLAG"] = 1.0 return dic @@ -1111,7 +1112,7 @@ def write_slice_3D(filemask, dic, data, shape, slices): # file doesn't exist, create a empty one ndata = np.zeros((dy, dx), dtype=data.dtype) write_single(f, dic, data, False) - del(ndata) + del ndata # mmap the [new] file mdata = np.memmap(f, dtype='float32', offset=512 * 4, mode='r+') @@ -1132,11 +1133,11 @@ def write_slice_3D(filemask, dic, data, shape, slices): if data.dtype == 'complex64': idata[sy, sx] = data.imag[i] idata.flush() - del(idata) + del idata # clean up - del(rdata) - del(mdata) + del rdata + del mdata # iter3D tools (xyz2pipe and pipe2xyz replacements) # Notes for iter3D implementation @@ -1260,6 +1261,7 @@ class iter3D: ziter.write("ft/test%03d.ft3",XZplane,dic) """ + def __init__(self, filemask, in_lead="x", out_lead="DEFAULT"): """ Create a iter3D object @@ -1806,9 +1808,9 @@ def __init__(self, filemask, order=(0, 1, 2), fcheck=False): f3 = "FDF" + str(int(dic["FDDIMORDER3"])) quadrature_factor = [2, 1][int(dic[f3 + 'QUADFLAG'])] - #Checking whether "nmrPipe -fn EXT ..." has been applied to z-dim or not. - #If EXT has been applied, FDF*XN is not zero. - #If z-dim is in time-domain, data-size given by FDF*X1 and FDF*XN has to be doubled. + # Checking whether "nmrPipe -fn EXT ..." has been applied to z-dim or not. + # If EXT has been applied, FDF*XN is not zero. + # If z-dim is in time-domain, data-size given by FDF*X1 and FDF*XN has to be doubled. if dic[f3 + 'FTFLAG']: if int(dic[f3 + 'XN']) == 0: @@ -1900,6 +1902,7 @@ class pipestream_3d(fileiobase.data_nd): Ordering of axes against file. """ + def __init__(self, filename, order=(0, 1, 2)): """ Create and set up object @@ -1995,6 +1998,7 @@ class pipe_4d(fileiobase.data_nd): set exist. Raises a IOError if files are missing. Default is False. """ + def __init__(self, filemask, order=(0, 1, 2, 3), fcheck=False): """ Create and set up object, check that files exist if fcheck is True @@ -2195,6 +2199,7 @@ def __fgetitem__(self, slices): f.close() return out + # data, see fdata.h fdata_nums = { 'FDMAGIC': '0', diff --git a/nmrglue/fileio/table.py b/nmrglue/fileio/table.py index 16cb9ece..16c1726b 100644 --- a/nmrglue/fileio/table.py +++ b/nmrglue/fileio/table.py @@ -144,7 +144,7 @@ def read(filename): dtd['formats'] = comments.pop(dl[0])[9:].split() # return the data as a records array - return comments, np.atleast_1d(np.recfromtxt(filename, dtype=dtd)) + return comments, np.atleast_1d(np.genfromtxt(filename, dtype=dtd)) def write(filename, comments, rec, overwrite=False): diff --git a/nmrglue/fileio/tecmag.py b/nmrglue/fileio/tecmag.py index a5931e3a..beab90c4 100644 --- a/nmrglue/fileio/tecmag.py +++ b/nmrglue/fileio/tecmag.py @@ -24,8 +24,8 @@ TNTMAGIC_RE = re.compile(br"^TNT1\.\d\d\d$") -TNTMAGIC = np.dtype('a8') -TNTTLV = np.dtype([('tag', 'a4'), ('bool', ' Date: Wed, 13 Nov 2024 11:04:05 +0200 Subject: [PATCH 03/38] One more numpy 2.0 fix... --- nmrglue/fileio/tecmag.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nmrglue/fileio/tecmag.py b/nmrglue/fileio/tecmag.py index beab90c4..c22e2f1a 100644 --- a/nmrglue/fileio/tecmag.py +++ b/nmrglue/fileio/tecmag.py @@ -93,7 +93,7 @@ ('DF_DN', ' Date: Wed, 27 Nov 2024 10:03:30 +0200 Subject: [PATCH 04/38] ...and more NumPy 2.0 fixes --- nmrglue/fileio/tecmag.py | 2 +- nmrglue/process/proc_lp.py | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/nmrglue/fileio/tecmag.py b/nmrglue/fileio/tecmag.py index c22e2f1a..34d7a5ba 100644 --- a/nmrglue/fileio/tecmag.py +++ b/nmrglue/fileio/tecmag.py @@ -97,7 +97,7 @@ ('LatchLP', ' Date: Wed, 22 Jan 2025 14:43:37 +0200 Subject: [PATCH 05/38] Added dtype=None to np.genfromtxt calls --- examples/fitting_data/t1_measurements/extract_trajs.py | 4 ++-- examples/fitting_data/t1_measurements/fit_exp_leastsq.py | 2 +- examples/fitting_data/t1_measurements/pt.py | 2 +- examples/integration/integrate_1d/integrate_1d.py | 2 +- examples/integration/integrate_2d/integrate_2d.py | 2 +- .../s12-s15_relaxation_analysis/extract_traj.py | 4 ++-- .../jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py | 4 ++-- .../s12-s15_relaxation_analysis/plot_trajectories.py | 2 +- examples/plotting/plot_2d/plot_assignments.py | 2 +- examples/plotting/plot_2d/plot_boxes.py | 2 +- examples/simulation/simulate_ucsf/make_ucsf.py | 2 +- 11 files changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/fitting_data/t1_measurements/extract_trajs.py b/examples/fitting_data/t1_measurements/extract_trajs.py index dbffd751..298e1872 100755 --- a/examples/fitting_data/t1_measurements/extract_trajs.py +++ b/examples/fitting_data/t1_measurements/extract_trajs.py @@ -5,8 +5,8 @@ import numpy as np # read in the integration limits and list of spectra -peak_list = np.genfromtxt("boxes.in", encoding="UTF-8") -spectra_list = np.genfromtxt("spectra.in", encoding="UTF-8") +peak_list = np.genfromtxt("boxes.in", dtype=None, encoding="UTF-8") +spectra_list = np.genfromtxt("spectra.in", dtype=None, encoding="UTF-8") # prepare the trajs records array num_spec = spectra_list.size diff --git a/examples/fitting_data/t1_measurements/fit_exp_leastsq.py b/examples/fitting_data/t1_measurements/fit_exp_leastsq.py index efeef81b..2bccc8c4 100755 --- a/examples/fitting_data/t1_measurements/fit_exp_leastsq.py +++ b/examples/fitting_data/t1_measurements/fit_exp_leastsq.py @@ -7,7 +7,7 @@ # read in the trajectories and times trajs = np.load("traj.npy") -t1 = np.genfromtxt("time.dat") +t1 = np.genfromtxt("time.dat", dtype=None) # fitting function and residual calculation diff --git a/examples/fitting_data/t1_measurements/pt.py b/examples/fitting_data/t1_measurements/pt.py index 3dec254c..96705870 100755 --- a/examples/fitting_data/t1_measurements/pt.py +++ b/examples/fitting_data/t1_measurements/pt.py @@ -24,7 +24,7 @@ def fit_func(p, x): # read in the trajectories, fitting results, and times fits = pickle.load(open("fits.pickle", "rb")) trajs = np.load("traj.npy") -times = np.genfromtxt("time.dat") +times = np.genfromtxt("time.dat", dtype=None) sim_times = np.linspace(times[0], times[-1], 2000) diff --git a/examples/integration/integrate_1d/integrate_1d.py b/examples/integration/integrate_1d/integrate_1d.py index c47bb239..5141b6d5 100755 --- a/examples/integration/integrate_1d/integrate_1d.py +++ b/examples/integration/integrate_1d/integrate_1d.py @@ -10,7 +10,7 @@ length = data.shape[0] # read in the integration limits -peak_list = np.genfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in", dtype=None) # determine the ppm scale uc = ng.pipe.make_uc(dic, data) diff --git a/examples/integration/integrate_2d/integrate_2d.py b/examples/integration/integrate_2d/integrate_2d.py index 13f81a88..71af4bc8 100755 --- a/examples/integration/integrate_2d/integrate_2d.py +++ b/examples/integration/integrate_2d/integrate_2d.py @@ -8,7 +8,7 @@ dic, data = ng.pipe.read("nmrpipe_2d/test.ft2") # read in the integration limits -peak_list = np.genfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in", dtype=None) # prepare the output file f = open("volumes.out", 'w') diff --git a/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py b/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py index 65110227..d3301bcb 100644 --- a/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py +++ b/examples/jbnmr_examples/s12-s15_relaxation_analysis/extract_traj.py @@ -2,8 +2,8 @@ import numpy as np # read the integration limits and list of spectra -peak_list = np.genfromtxt("boxes.in", names=True) -spectra_list = np.genfromtxt("spectra.in") +peak_list = np.genfromtxt("boxes.in", dtype=None, names=True) +spectra_list = np.genfromtxt("spectra.in", dtype=None) # create an array to hold the trajectories trajectories = np.empty((peak_list.size, spectra_list.size), dtype='float') diff --git a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py index 939b79d5..4a923f96 100644 --- a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py +++ b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_boxes.py @@ -12,8 +12,8 @@ cl = 30000 * 1.20 ** np.arange(20) # read in the box limits and list of spectra -peak_list = np.genfromtxt("boxes.in", names=True) -spectra_list = np.genfromtxt("spectra.in") +peak_list = np.genfromtxt("boxes.in", dtype=None, names=True) +spectra_list = np.genfromtxt("spectra.in", dtype=None) # loop over the spectra for spec_number, spectra in enumerate(spectra_list): diff --git a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py index d67d2dc7..3aeef6c4 100644 --- a/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py +++ b/examples/jbnmr_examples/s12-s15_relaxation_analysis/plot_trajectories.py @@ -9,7 +9,7 @@ def fit_func(p, x): return A * np.exp(-1.0 * np.array(x) * R2 / 1.0e6) -fitting_results = np.genfromtxt('fits.txt') +fitting_results = np.genfromtxt('fits.txt', dtype=None) experimental_relaxation_times = np.loadtxt("relaxation_times.in") simulated_relaxation_times = np.linspace(0, 4000000, 2000) diff --git a/examples/plotting/plot_2d/plot_assignments.py b/examples/plotting/plot_2d/plot_assignments.py index 9abf5ffa..4f20a039 100755 --- a/examples/plotting/plot_2d/plot_assignments.py +++ b/examples/plotting/plot_2d/plot_assignments.py @@ -20,7 +20,7 @@ dic, data = ng.pipe.read("nmrpipe_2d/test.ft2") # read in the integration limits -peak_list = np.genfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in", dtype=None) # create the figure fig = plt.figure() diff --git a/examples/plotting/plot_2d/plot_boxes.py b/examples/plotting/plot_2d/plot_boxes.py index ebd4e6bf..fa4da3da 100755 --- a/examples/plotting/plot_2d/plot_boxes.py +++ b/examples/plotting/plot_2d/plot_boxes.py @@ -21,7 +21,7 @@ dic, data = ng.pipe.read("nmrpipe_2d/test.ft2") # read in the integration limits -peak_list = np.genfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in", dtype=None) # loop over the peaks for name, x0, y0, x1, y1 in peak_list: diff --git a/examples/simulation/simulate_ucsf/make_ucsf.py b/examples/simulation/simulate_ucsf/make_ucsf.py index 892b8954..a783ab40 100644 --- a/examples/simulation/simulate_ucsf/make_ucsf.py +++ b/examples/simulation/simulate_ucsf/make_ucsf.py @@ -32,7 +32,7 @@ data = np.empty((512, 1024), dtype='float32') # read in the peak list -peak_list = np.genfromtxt('peaks.txt', names=True) +peak_list = np.genfromtxt('peaks.txt', dtype=None, names=True) npeaks = len(peak_list) # convert the peak list from PPM to points From 439b0e641b37767b8bd56a4df6898925c490a91c Mon Sep 17 00:00:00 2001 From: Maria Hecker Date: Thu, 26 Jun 2025 16:01:58 +0000 Subject: [PATCH 06/38] Fix deprecated scipy calls --- nmrglue/process/proc_base.py | 2 +- nmrglue/process/proc_bl.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/nmrglue/process/proc_base.py b/nmrglue/process/proc_base.py index d2bcfd6b..675ecda7 100644 --- a/nmrglue/process/proc_base.py +++ b/nmrglue/process/proc_base.py @@ -2644,7 +2644,7 @@ def zd_gaussian(data, wide=1.0, x0=0.0, slope=1.0, g=1): """ tln2 = np.sqrt(2 * np.log(2)) - window = 1 - scipy.signal.gaussian(2 * wide + 1, g / tln2) + window = 1 - scipy.signal.windows.gaussian(2 * wide + 1, g / tln2) return zd(data, window, x0=x0, slope=slope) diff --git a/nmrglue/process/proc_bl.py b/nmrglue/process/proc_bl.py index 1f2dc908..4d6f1622 100644 --- a/nmrglue/process/proc_bl.py +++ b/nmrglue/process/proc_bl.py @@ -206,7 +206,7 @@ def calc_bl_med(x, mw, sf, sigma): # described algorithm but is MUCH faster # convolve with a gaussian - g = scipy.signal.gaussian(sf, sigma) + g = scipy.signal.windows.gaussian(sf, sigma) g = g / g.sum() return scipy.signal.convolve(m, g, mode='same') @@ -360,7 +360,7 @@ def sol_gaussian(data, w=16, mode='same'): """ Solvent filter with square gaussian filter. See :py:func:`sol_general`. """ - filter = scipy.signal.gaussian(w, w / 2.) + filter = scipy.signal.windows.gaussian(w, w / 2.) return sol_general(data, filter, w=w, mode=mode) From 6ccb7a7d5a8623c1551a2411a69ac11cf7db4d02 Mon Sep 17 00:00:00 2001 From: Angela Maennel Date: Thu, 3 Jul 2025 14:01:34 +0200 Subject: [PATCH 07/38] try prefered encoding first but also other --- nmrglue/fileio/bruker.py | 65 +++++++++++++++++++++++++--------------- 1 file changed, 41 insertions(+), 24 deletions(-) diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index 9de2f52a..e97b79f0 100644 --- a/nmrglue/fileio/bruker.py +++ b/nmrglue/fileio/bruker.py @@ -2155,7 +2155,7 @@ def read_jcamp(filename, encoding=locale.getpreferredencoding()): filename : str Filename of Bruker JCAMP-DX file. encoding : str - Encoding of Bruker JCAMP-DX file. Defaults to the system default locale + Encoding of Bruker JCAMP-DX file. Defaults to the system default locale. Returns ------- @@ -2173,32 +2173,49 @@ def read_jcamp(filename, encoding=locale.getpreferredencoding()): """ dic = {"_coreheader": [], "_comments": []} # create empty dictionary - - with open(filename, 'r', encoding=encoding) as f: - while True: # loop until end of file is found - - line = f.readline().rstrip() # read a line - if line == '': # end of file found - break - - if line[:6] == "##END=": - # print("End of file") - break - elif line[:2] == "$$": - dic["_comments"].append(line) - elif line[:2] == "##" and line[2] != "$": - dic["_coreheader"].append(line) - elif line[:3] == "##$": - try: - key, value = parse_jcamp_line(line, f) - dic[key] = value - except: - warn(f"Unable to correctly parse line: {line}") - else: - warn(f"Extraneous line: {line}") + try: + with open(filename, 'r', encoding=encoding) as f: + dic=parse_jcamp_file(f,dic) + except: + if encoding == "utf-8": + with open(filename, 'r', encoding="cp1252") as f: + dic=parse_jcamp_file(f,dic) + else: + with open(filename, 'r', encoding="utf-8") as f: + dic=parse_jcamp_file(f,dic) return dic +def parse_jcamp_file(f,dic): + """ + Parse a Bruker JCAMP-DX file into the given dictionary. + + The dictionary shoul contain the keys _coreheader and _comments intialized + with empty lists. Bruker parameter "$FOO" are extracted into strings, + floats or lists and assigned to dic["FOO"] + """ + while True: # loop until end of file is found + + line = f.readline().rstrip() # read a line + if line == '': # end of file found + break + + if line[:6] == "##END=": + # print("End of file") + break + elif line[:2] == "$$": + dic["_comments"].append(line) + elif line[:2] == "##" and line[2] != "$": + dic["_coreheader"].append(line) + elif line[:3] == "##$": + try: + key, value = parse_jcamp_line(line, f) + dic[key] = value + except: + warn(f"Unable to correctly parse line: {line}") + else: + warn(f"Extraneous line: {line}") + return dic def parse_jcamp_line(line, f): """ From c47f249bf6fd4feefb4676a556dfbfba82acd288 Mon Sep 17 00:00:00 2001 From: Juuso Lehtivarjo Date: Mon, 14 Jul 2025 13:54:02 +0300 Subject: [PATCH 08/38] Fix new Numpy deprecation (binary mode of np.fromstring) --- nmrglue/fileio/nmrml.py | 2 +- nmrglue/fileio/tecmag.py | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/nmrglue/fileio/nmrml.py b/nmrglue/fileio/nmrml.py index 28dfb450..3956301f 100644 --- a/nmrglue/fileio/nmrml.py +++ b/nmrglue/fileio/nmrml.py @@ -76,6 +76,6 @@ def _get_nmrml_data(fid_dict, data_dtype): # convert to ndarray if data_dtype is None: data_dtype = BYTE_FORMAT_TO_DTYPE[byte_format] - data = np.fromstring(uncompressed_bytes, dtype=data_dtype) + data = np.frombuffer(uncompressed_bytes, dtype=data_dtype) return data diff --git a/nmrglue/fileio/tecmag.py b/nmrglue/fileio/tecmag.py index 34d7a5ba..d2d9b56c 100644 --- a/nmrglue/fileio/tecmag.py +++ b/nmrglue/fileio/tecmag.py @@ -253,7 +253,7 @@ def read(filename): with open(filename, 'rb') as tntfile: - tntmagic = np.fromstring(tntfile.read(TNTMAGIC.itemsize), + tntmagic = np.frombuffer(tntfile.read(TNTMAGIC.itemsize), TNTMAGIC, count=1)[0] if not TNTMAGIC_RE.match(tntmagic): @@ -263,21 +263,21 @@ def read(filename): # Read in the section headers tnthdrbytes = tntfile.read(TNTTLV.itemsize) - while (TNTTLV.itemsize == len(tnthdrbytes)): - tlv = np.fromstring(tnthdrbytes, TNTTLV)[0] + while TNTTLV.itemsize == len(tnthdrbytes): + tlv = np.frombuffer(tnthdrbytes, TNTTLV)[0] data_length = tlv['length'] hdrdict = {'offset': tntfile.tell(), 'length': data_length, 'bool': bool(tlv['bool'])} if data_length <= 4096: hdrdict['data'] = tntfile.read(data_length) - assert (len(hdrdict['data']) == data_length) + assert len(hdrdict['data']) == data_length else: tntfile.seek(data_length, os.SEEK_CUR) tnt_sections[tlv['tag'].decode()] = hdrdict tnthdrbytes = tntfile.read(TNTTLV.itemsize) - assert (tnt_sections['TMAG']['length'] == TNTTMAG.itemsize) + assert tnt_sections['TMAG']['length'] == TNTTMAG.itemsize tmag = np.fromstring(tnt_sections['TMAG']['data'], TNTTMAG, count=1)[0] assert (tnt_sections['DATA']['length'] == @@ -291,7 +291,7 @@ def read(filename): shape=tmag['actual_npts'].prod()) data = np.reshape(data, tmag['actual_npts'], order='F') - assert (tnt_sections['TMG2']['length'] == TNTTMG2.itemsize) + assert tnt_sections['TMG2']['length'] == TNTTMG2.itemsize tmg2 = np.fromstring(tnt_sections['TMG2']['data'], TNTTMG2, count=1)[0] dic = dict() From 725a1a10f99d5c4fc1b502d94792f108dc7f05f3 Mon Sep 17 00:00:00 2001 From: Jiwoo Seo Date: Fri, 19 Sep 2025 18:48:11 -0400 Subject: [PATCH 09/38] bugfix: Fix zero-shift bug in `proc_base.ls` - Fixed a bug where `proc_base.ls()` returned zeroed out data when left shifting by zero (pts=0), by adding a guard clause. - Minor changes: Changed the default parameter from `pts=0.0` to `pts=0` for clarity in a few functions, especially with the docstring describing `pts` as int --- nmrglue/process/proc_base.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/nmrglue/process/proc_base.py b/nmrglue/process/proc_base.py index d2bcfd6b..ff8bb3d6 100644 --- a/nmrglue/process/proc_base.py +++ b/nmrglue/process/proc_base.py @@ -324,7 +324,7 @@ def tri(data, loc="auto", lHi=0.0, rHi=0.0, inv=False, rev=False): ################### -def rs(data, pts=0.0): +def rs(data, pts=0): """ Right shift and zero fill. @@ -345,12 +345,13 @@ def rs(data, pts=0.0): roll : shift without zero filling. """ + data = np.roll(data, int(pts), axis=-1) data[..., :int(pts)] = 0 return data -def ls(data, pts=0.0): +def ls(data, pts=0): """ Left shift and fill with zero @@ -371,12 +372,14 @@ def ls(data, pts=0.0): roll : shift without zero filling. """ + if int(pts) == 0: + return data data = np.roll(data, -int(pts), axis=-1) data[..., -int(pts):] = 0 return data -def cs(data, pts=0.0, neg=False): +def cs(data, pts=0, neg=False): """ Circular shift @@ -399,7 +402,7 @@ def cs(data, pts=0.0, neg=False): return roll(data, pts, neg) -def roll(data, pts=0.0, neg=False): +def roll(data, pts=0, neg=False): """ Roll axis From 4d3599b5ed0e43c683ac4f00c9ac2330e3c81317 Mon Sep 17 00:00:00 2001 From: Petr Padrta Date: Sat, 29 Nov 2025 15:10:14 +0100 Subject: [PATCH 10/38] write() was ignoring overwrite flag for proc* files --- nmrglue/fileio/bruker.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index e97b79f0..ef6c3a5e 100644 --- a/nmrglue/fileio/bruker.py +++ b/nmrglue/fileio/bruker.py @@ -768,9 +768,9 @@ def write(dir, dic, data, bin_file=None, acqus_files=None, procs_files=None, pdata_path = dir for f in procs_files: - write_jcamp(dic[f], os.path.join(pdata_path, f)) + write_jcamp(dic[f], os.path.join(pdata_path, f), overwrite=overwrite) for f in proc_files: - write_jcamp(dic[f+'s'], os.path.join(pdata_path, f)) + write_jcamp(dic[f+'s'], os.path.join(pdata_path, f), overwrite=overwrite) # write out the pulse program if write_prog: @@ -2723,4 +2723,4 @@ def read_vdlist(dirc, fname='vdlist'): # Convert to floats vdlist = [float(i) for i in vdlist] - return vdlist \ No newline at end of file + return vdlist From 833e14a204bcedd15877a00b81e7e87033bd3c94 Mon Sep 17 00:00:00 2001 From: Jonathan Helmus Date: Tue, 23 Dec 2025 14:37:33 -0600 Subject: [PATCH 11/38] fix unpacking in jeol --- nmrglue/fileio/jeol.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nmrglue/fileio/jeol.py b/nmrglue/fileio/jeol.py index 53f60d62..f6788689 100644 --- a/nmrglue/fileio/jeol.py +++ b/nmrglue/fileio/jeol.py @@ -216,7 +216,7 @@ def truncate_data(dic, data): slices = slices[::-1] - return data[*slices] + return data[tuple(slices)] def get_data_shape(dic): From 41d764f4304773f0b7ae5bf186e9befe53831d7c Mon Sep 17 00:00:00 2001 From: Jonathan Helmus Date: Tue, 23 Dec 2025 14:53:00 -0600 Subject: [PATCH 12/38] enable ci workflow on PRs --- .github/workflows/ci.yml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 19986753..5439f643 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,6 +1,10 @@ name: nmrglue tests -on: [push] +on: + push: + branches: master + pull_request: + branches: master jobs: build: From 028a7861aaf26d33b38bd2af3e4a84464b2e5313 Mon Sep 17 00:00:00 2001 From: Jonathan Helmus Date: Tue, 23 Dec 2025 14:56:23 -0600 Subject: [PATCH 13/38] do not fail fast --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 5439f643..de029f44 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,6 +12,7 @@ jobs: strategy: matrix: python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + fail-fast: false steps: - uses: actions/checkout@v3 From 2d0c507697c33d5467de158daf1f84c1604fa1f7 Mon Sep 17 00:00:00 2001 From: "Kaustubh R. Mote" Date: Thu, 25 Dec 2025 09:17:00 +0530 Subject: [PATCH 14/38] fixed deprecated fromstring call --- nmrglue/fileio/tecmag.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nmrglue/fileio/tecmag.py b/nmrglue/fileio/tecmag.py index d2d9b56c..0d2e547a 100644 --- a/nmrglue/fileio/tecmag.py +++ b/nmrglue/fileio/tecmag.py @@ -278,7 +278,7 @@ def read(filename): tnthdrbytes = tntfile.read(TNTTLV.itemsize) assert tnt_sections['TMAG']['length'] == TNTTMAG.itemsize - tmag = np.fromstring(tnt_sections['TMAG']['data'], TNTTMAG, count=1)[0] + tmag = np.frombuffer(tnt_sections['TMAG']['data'], TNTTMAG, count=1)[0] assert (tnt_sections['DATA']['length'] == tmag['actual_npts'].prod() * 8) @@ -292,7 +292,7 @@ def read(filename): data = np.reshape(data, tmag['actual_npts'], order='F') assert tnt_sections['TMG2']['length'] == TNTTMG2.itemsize - tmg2 = np.fromstring(tnt_sections['TMG2']['data'], TNTTMG2, count=1)[0] + tmg2 = np.frombuffer(tnt_sections['TMG2']['data'], TNTTMG2, count=1)[0] dic = dict() for name in TNTTMAG.names: From 295c95399b1d4c86988554a2c9478b1a0efe65ee Mon Sep 17 00:00:00 2001 From: "Kaustubh R. Mote" Date: Fri, 25 Apr 2025 08:05:49 +0530 Subject: [PATCH 15/38] better parsing of topspin version --- nmrglue/fileio/bruker.py | 45 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index ef6c3a5e..ea93f66b 100644 --- a/nmrglue/fileio/bruker.py +++ b/nmrglue/fileio/bruker.py @@ -2724,3 +2724,48 @@ def read_vdlist(dirc, fname='vdlist'): vdlist = [float(i) for i in vdlist] return vdlist + +def guess_topspin_version(dic): + """ + Guess the version of topspin on which the data was acquired + + Parameters + ---------- + dic : dict + dictionary associated with the data + + Returns + ------- + tuple + (version, tag, instrument) + + """ + version = dic["acqus"]["_coreheader"][0] + version = version.split("##TITLE= Parameter file, ")[1] + version = ( + version.replace(" ", "") + .replace("\t", "") + .replace("Version", "") + .replace("TopSpin", "topspin.") + .replace("TOPSPIN", "topspin.") + ) + if "XWIN-NMR" in version: + version = version.replace("XWIN-NMR", "xwin-nmr.").split(".") + else: + version = version.replace(" ", "").replace("pl", ".pl.").split(".") + + parsed_version = [] + for i in version: + try: + parsed_version.append(int(i)) + except (TypeError, ValueError): + parsed_version.append(i) + + version = f"{parsed_version[1]}.{parsed_version[2]}" + instrument = parsed_version[0] + try: + tag = ''.join(str(i) for i in parsed_version[3:]) + except IndexError: + tag = '' + + return version, tag, instrument From fe210eebfca18687e6c1f4bc71fc60da43b54e22 Mon Sep 17 00:00:00 2001 From: "Kaustubh R. Mote" Date: Thu, 24 Apr 2025 19:16:09 +0530 Subject: [PATCH 16/38] correctly use DTYPA for determining the data shape --- nmrglue/fileio/bruker.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index ea93f66b..04190d33 100644 --- a/nmrglue/fileio/bruker.py +++ b/nmrglue/fileio/bruker.py @@ -1052,6 +1052,12 @@ def guess_shape(dic): except KeyError: dtypa = 0 # default value, int32 data + version = guess_topspin_version(dic)[1] + if (version >= 4) and (dtypa == 2): + bytesize = 8 + else: + bytesize = 4 + # last (direct) dimension is given by "TD" parameter in acqus file # rounded up to nearest (1024/(bytes per point)) # next-to-last dimension may be given by "TD" in acqu2s. In 3D+ data @@ -1063,8 +1069,8 @@ def guess_shape(dic): # additional dimension given by data size if shape[2] != 0 and shape[3] != 0: - shape[1] = fsize // (shape[3] * shape[2] * 4) - shape[0] = fsize // (shape[3] * shape[2] * shape[1] * 4) + shape[1] = fsize // (shape[3] * shape[2] * bytesize) + shape[0] = fsize // (shape[3] * shape[2] * shape[1] * bytesize) # if there in no pulse program parameters in dictionary return current # shape after removing zeros From 363724c9f34705621ac3bdceff30c78374ea642a Mon Sep 17 00:00:00 2001 From: "Kaustubh R. Mote" Date: Tue, 30 Dec 2025 12:30:27 +0530 Subject: [PATCH 17/38] data reshaping now handles NUS datasets correctly --- nmrglue/fileio/bruker.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index 04190d33..cd175fef 100644 --- a/nmrglue/fileio/bruker.py +++ b/nmrglue/fileio/bruker.py @@ -426,22 +426,24 @@ def read(dir=".", bin_file=None, acqus_files=None, pprog_file=None, shape=None, else: isfloat = False - # read the binary file - f = os.path.join(dir, bin_file) - _, data = read_binary(f, shape=shape, cplex=cplex, big=big, - isfloat=isfloat) - try: if dic['acqus']['FnTYPE'] == 2: # non-uniformly sampled data try: dic['nuslist'] = read_nuslist(dir) except FileNotFoundError: warn("NUS data detected, but nuslist was not found") + shape = (-1, shape[-1]) except KeyError: # old datasets do not have the FnTYPE parameter in acqus files. # also fails silently when acqus file is absent. pass + # read the binary file + f = os.path.join(dir, bin_file) + _, data = read_binary(f, shape=shape, cplex=cplex, big=big, + isfloat=isfloat) + + return dic, data @@ -1052,20 +1054,20 @@ def guess_shape(dic): except KeyError: dtypa = 0 # default value, int32 data - version = guess_topspin_version(dic)[1] - if (version >= 4) and (dtypa == 2): + if dtypa == 0: + bytesize = 4 + elif dtypa == 2: bytesize = 8 else: - bytesize = 4 + raise ValueError(f'DTYPA ({dtypa}) is inconsistent with expected values of 0 or 2') # last (direct) dimension is given by "TD" parameter in acqus file # rounded up to nearest (1024/(bytes per point)) # next-to-last dimension may be given by "TD" in acqu2s. In 3D+ data # this is often the sum of the indirect dimensions - if dtypa == 2: - shape = [0, 0, td2, int(np.ceil(td0 / 128.) * 128.)] - else: - shape = [0, 0, td2, int(np.ceil(td0 / 256.) * 256.)] + + pointsize = 1024.0 / bytesize + shape = [0, 0, td2, int(np.ceil(td0 / pointsize) * pointsize)] # additional dimension given by data size if shape[2] != 0 and shape[3] != 0: From 6571048a3655fcff7528e13c01596af7a8fdb2c0 Mon Sep 17 00:00:00 2001 From: "Kaustubh R. Mote" Date: Wed, 31 Dec 2025 15:52:50 +0530 Subject: [PATCH 18/38] added a function to update the unit conversion object --- nmrglue/fileio/fileiobase.py | 42 ++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/nmrglue/fileio/fileiobase.py b/nmrglue/fileio/fileiobase.py index d53ca4ef..e83262f2 100644 --- a/nmrglue/fileio/fileiobase.py +++ b/nmrglue/fileio/fileiobase.py @@ -330,6 +330,48 @@ def us_scale(self): __call__ = i # calling the object x is the same as x.i +def update_uc(uc, **kwargs): + """ + Updates the unit_conversion object with new parameters + and returns a new unit_conversion object + + Parameters + ---------- + uc : unit_conversion + unit_conversion object to update + + Returns + ------- + unit_conversion + Updated unit_conversion object + + Raises + ------ + TypeError + If the provided uc is not a unit conversion object + TypeError + If the parameter to update is not a valid parameter. + + """ + if not isinstance(uc, unit_conversion): + raise TypeError(f'{uc} is not a unit_conversion object.') + + if not kwargs: + return uc + + params = ('size', 'cplx', 'sw', 'obs', 'car') + new_params = {p: uc.__getattribute__(f'_{p}') for p in params} + + for k, v in kwargs.items(): + if k in params: + new_params[k] = v + else: + raise TypeError( + f'{k} cannot be used to construct a unit conversion object. Acceptable keywords are {(params)}.') + + return unit_conversion(*new_params.values()) + + def uc_from_udic(udic, dim=-1): """ Create a unit conversion object from a Universal dictionary. From 986e98549415322d68fd7c8a21f929430e42095d Mon Sep 17 00:00:00 2001 From: "Kaustubh R. Mote" Date: Wed, 31 Dec 2025 16:00:16 +0530 Subject: [PATCH 19/38] added a test for update_uc --- nmrglue/fileio/tests/test_fileiobase.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/nmrglue/fileio/tests/test_fileiobase.py b/nmrglue/fileio/tests/test_fileiobase.py index 0206ca65..0de21894 100644 --- a/nmrglue/fileio/tests/test_fileiobase.py +++ b/nmrglue/fileio/tests/test_fileiobase.py @@ -46,3 +46,14 @@ def test_uc_with_float_size(): size=64.0, cplx=False, sw=1.0, obs=1.0, car=1.0) scale = uc.ppm_scale() assert len(scale) == 64 + + +def test_update_uc(): + uc = ng.fileiobase.unit_conversion( + size=64.0, cplx=False, sw=1.0, obs=1.0, car=1.0) + uc2 = ng.fileiobase.update_uc(uc, size=10, cplx=True, sw=2.0, car=5.2, obs=3.0) + assert uc2._size == 10 + assert uc2._cplx is True + assert abs(uc2._sw - 2.0) < 1e-5 + assert abs(uc2._car - 5.2) < 1e-5 + assert abs(uc2._obs - 3.0) < 1e-5 \ No newline at end of file From 8e21c021440f3f8b6699b714435b0e5d62a69ff0 Mon Sep 17 00:00:00 2001 From: Simon Dollinger Date: Mon, 19 Jan 2026 16:28:47 +0100 Subject: [PATCH 20/38] Ignore BOM signature when parsing jdx files --- nmrglue/fileio/jcampdx.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nmrglue/fileio/jcampdx.py b/nmrglue/fileio/jcampdx.py index 1d41e2b5..6582fa07 100644 --- a/nmrglue/fileio/jcampdx.py +++ b/nmrglue/fileio/jcampdx.py @@ -51,7 +51,7 @@ def _parsejcampdx(filename): # when encountering ##END, push the ready dict to another list readyblocklist = [] - filein = open(filename, 'r', encoding="utf-8", errors="replace") + filein = open(filename, 'r', encoding="utf-8-sig", errors="replace") currentkey = None currentvaluestrings = [] From ab04a8a2af32cbd9621188c81dad5d58ac2f1e3d Mon Sep 17 00:00:00 2001 From: "Kaustubh R. Mote" Date: Mon, 9 Feb 2026 09:29:17 +0530 Subject: [PATCH 21/38] handles partially acquired datasets in bruker --- nmrglue/fileio/bruker.py | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index cd175fef..c2de4573 100644 --- a/nmrglue/fileio/bruker.py +++ b/nmrglue/fileio/bruker.py @@ -438,10 +438,15 @@ def read(dir=".", bin_file=None, acqus_files=None, pprog_file=None, shape=None, # also fails silently when acqus file is absent. pass + estimated_dims = 1 + for i in (2, 3, 4): + if f'acqu{i}s' in dic.keys(): + estimated_dims += 1 + # read the binary file f = os.path.join(dir, bin_file) _, data = read_binary(f, shape=shape, cplex=cplex, big=big, - isfloat=isfloat) + isfloat=isfloat, estimated_dims=estimated_dims) return dic, data @@ -1529,7 +1534,7 @@ def reorder_submatrix(data, shape, submatrix_shape, reverse=False): # Bruker binary (fid/ser) reading and writing -def read_binary(filename, shape=(1), cplex=True, big=True, isfloat=False): +def read_binary(filename, shape=(1), cplex=True, big=True, isfloat=False, estimated_dims=None): """ Read Bruker binary data from file and return dic,data pair. @@ -1578,7 +1583,18 @@ def read_binary(filename, shape=(1), cplex=True, big=True, isfloat=False): return dic, data.reshape(shape) except ValueError: - warn(f"{data.shape} cannot be shaped into {shape}") + try: + data = data.reshape(-1, shape[-1]) + if estimated_dims and (estimated_dims > 2): + warn( + "Data is inconsistent with acquistion parameters. " + "This usually happens with partially acquired datasets with dims >2. " + "A 2D dataset will be returned. This will requires further reshaping before processing." + ) + return dic, data + except ValueError: + warn(f"{data.shape} cannot be shaped into {shape} or a consistent 2D array. A 1D array will be returned.") + return dic, data From fad974a2c6cc6527123ea94d2cf3e08fbce40235 Mon Sep 17 00:00:00 2001 From: Daniel Pizetta Date: Thu, 28 May 2026 21:11:19 -0300 Subject: [PATCH 22/38] Add guess_tables() to tecmag for reading sequence tables from PSEQ Closes #123. The PSEQ section of Tecmag .tnt files stores sequence tables (rfamp, delay tables, etc.) as pairs of length-prefixed Pascal strings. This adds guess_tables() to expose them as a dict. The pytnt library (chatcannon/pytnt, GPL-3.0) reads delay tables by matching the default TNMR naming pattern `de[0-9]+:[0-9]`. As noted by robbyblum in #123, most users assign descriptive names to their tables (e.g. "rfamp", "tp"), which pytnt's regex does not catch. This implementation generalizes the parser to find any valid Pascal-string pair in PSEQ, with strict length-prefix validation to reject coincidental substring matches. The function follows the API suggested by kaustubhmote in #123: a standalone guess_tables() function with optional name filtering. Single-value (scalar) entries are excluded since they typically represent sequence variable defaults rather than swept tables. SI-prefix conversion (e.g. '100u' -> 1e-4) mirrors pytnt's convert_si() behavior so unit-suffixed tokens are returned in base SI units. Tests cover the numeric parsing, SI conversion, scalar exclusion, name filtering, the convert_si=False mode, files without a PSEQ section, rejection of invalid length prefixes, rejection of non-ASCII names, and mixed whitespace separators. --- nmrglue/fileio/tecmag.py | 194 +++++++++++++++- nmrglue/fileio/tests/test_tecmag.py | 331 ++++++++++++++++++++++++++++ 2 files changed, 524 insertions(+), 1 deletion(-) create mode 100644 nmrglue/fileio/tests/test_tecmag.py diff --git a/nmrglue/fileio/tecmag.py b/nmrglue/fileio/tecmag.py index 0d2e547a..11eb0dcf 100644 --- a/nmrglue/fileio/tecmag.py +++ b/nmrglue/fileio/tecmag.py @@ -15,6 +15,7 @@ import io import os +import struct import re import numpy as np @@ -38,7 +39,7 @@ ('repeat_times', ' 1e-4. +# Same set of prefixes used by pytnt (chatcannon/pytnt, GPL-3.0). +_SI_PREFIX = { + 'y': 1e-24, 'z': 1e-21, 'a': 1e-18, 'f': 1e-15, 'p': 1e-12, + 'n': 1e-9, 'u': 1e-6, 'm': 1e-3, 'c': 1e-2, 'd': 1e-1, + 's': 1.0, 'k': 1e3, 'M': 1e6, 'G': 1e9, 'T': 1e12, + 'P': 1e15, 'E': 1e18, 'Z': 1e21, 'Y': 1e24, +} + +# Maximum allowed length of a single table name (characters) +_NAME_MAX_LEN = 64 +# Maximum allowed length of a single table's content (bytes) +_CONTENT_MAX_LEN = 1_000_000 + + +def _convert_si(tokens): + """Convert a list of strings to a numpy float64 array, handling SI suffixes. + + Each token is parsed as a plain float; if that fails and the last character + is an SI prefix (y, z, ..., Y), the rest is parsed as float and multiplied + by the prefix value. Raises ValueError on unparseable tokens. + """ + out = np.empty(len(tokens), dtype=np.float64) + for i, t in enumerate(tokens): + try: + out[i] = float(t) + except ValueError: + if t and t[-1] in _SI_PREFIX: + out[i] = _SI_PREFIX[t[-1]] * float(t[:-1]) + else: + raise ValueError( + "could not convert {0!r} to a number".format(t)) + return out + + +def _read_pascal_pair(buf, pos): + """Try to read two consecutive Pascal strings at byte offset ``pos``. + + A Pascal string is a 4-byte little-endian length followed by that many + bytes of ASCII content. Returns ``(name, content, next_pos)`` if both + strings validate, else ``None``. + + Validation rules: + - Name length must be 1..64 inclusive. + - Name bytes must be printable ASCII (0x20..0x7e). + - Name must not begin with a digit (table names start with a letter + or underscore in TNMR). + - Content length must be 1..1_000_000 inclusive. + - Content bytes must be printable ASCII or whitespace. + """ + if pos < 0 or pos + 4 > len(buf): + return None + name_len = struct.unpack_from(' len(buf): + return None + content_len = struct.unpack_from('>> import nmrglue as ng + >>> tables = ng.tecmag.guess_tables('data.tnt') + >>> sorted(tables.keys()) + ['d1', 'rfamp'] + >>> tables['rfamp'] + array([ 1. , 4.22, 7.44, 10.67, 13.89, 17.11, 20.33, 23.56, 26.78, 30. ]) + >>> ng.tecmag.guess_tables('data.tnt', names=['rfamp']) + {'rfamp': array([1. , 4.22, 7.44, ...])} + """ + with open(filename, 'rb') as f: + data = f.read() + + # Locate PSEQ section payload (skip 4-byte tag + 4-byte bool + 4-byte length) + pseq_pos = data.find(b'PSEQ') + if pseq_pos < 0: + return {} + region_start = pseq_pos + 4 + 4 + 4 + + name_filter = set(names) if names is not None else None + tables = {} + + pos = region_start + while pos < len(data) - 8: + pair = _read_pascal_pair(data, pos) + if pair is None: + pos += 1 + continue + name, content, next_pos = pair + + # Skip if user requested a name filter and this isn't on it + if name_filter is not None and name not in name_filter: + pos = next_pos + continue + # Skip duplicate names; keep first occurrence + if name in tables: + pos = next_pos + continue + # Skip scalars: whitespace-split content must have >1 token + tokens = content.split() + if len(tokens) <= 1: + pos = next_pos + continue + + if convert_si: + try: + tables[name] = _convert_si(tokens) + except ValueError: + # Non-numeric content; keep as token list + tables[name] = tokens + else: + tables[name] = tokens + pos = next_pos + + return tables diff --git a/nmrglue/fileio/tests/test_tecmag.py b/nmrglue/fileio/tests/test_tecmag.py new file mode 100644 index 00000000..6e63087a --- /dev/null +++ b/nmrglue/fileio/tests/test_tecmag.py @@ -0,0 +1,331 @@ +"""Unit tests for nmrglue/fileio/tecmag.py module. + +These tests use hand-crafted byte sequences to exercise the parser +without requiring external .tnt test data. The synthetic files mimic +the minimum structure needed: the TNT1.NNN magic header, TLV section +records for TMAG/PSEQ/DATA/TMG2, and valid Pascal-string table pairs +inside PSEQ. +""" +import os +import struct +import tempfile + +import numpy as np +from numpy.testing import assert_allclose, assert_array_equal +import pytest + +import nmrglue as ng +from nmrglue.fileio.tecmag import ( + TNTTMAG, TNTTMG2, guess_tables, _read_pascal_pair, _convert_si +) + + +# ============================================================================= +# Helpers — synthesize minimal valid .tnt byte sequences +# ============================================================================= + +def _pascal_string(s): + """Build a length-prefixed Pascal string.""" + if isinstance(s, str): + s = s.encode('ascii') + return struct.pack(' Date: Fri, 19 Jun 2026 12:32:23 -0300 Subject: [PATCH 23/38] Move guess_table tests to the main test file --- nmrglue/fileio/tests/test_tecmag.py | 331 --------------------- tests/test_tecmag.py | 440 +++++++++++++++++++++++++++- 2 files changed, 433 insertions(+), 338 deletions(-) delete mode 100644 nmrglue/fileio/tests/test_tecmag.py diff --git a/nmrglue/fileio/tests/test_tecmag.py b/nmrglue/fileio/tests/test_tecmag.py deleted file mode 100644 index 6e63087a..00000000 --- a/nmrglue/fileio/tests/test_tecmag.py +++ /dev/null @@ -1,331 +0,0 @@ -"""Unit tests for nmrglue/fileio/tecmag.py module. - -These tests use hand-crafted byte sequences to exercise the parser -without requiring external .tnt test data. The synthetic files mimic -the minimum structure needed: the TNT1.NNN magic header, TLV section -records for TMAG/PSEQ/DATA/TMG2, and valid Pascal-string table pairs -inside PSEQ. -""" -import os -import struct -import tempfile - -import numpy as np -from numpy.testing import assert_allclose, assert_array_equal -import pytest - -import nmrglue as ng -from nmrglue.fileio.tecmag import ( - TNTTMAG, TNTTMG2, guess_tables, _read_pascal_pair, _convert_si -) - - -# ============================================================================= -# Helpers — synthesize minimal valid .tnt byte sequences -# ============================================================================= - -def _pascal_string(s): - """Build a length-prefixed Pascal string.""" - if isinstance(s, str): - s = s.encode('ascii') - return struct.pack(' 0) and np.any(got < 0) + assert_allclose(got, [-0.5, -0.1, 0.0, 0.5, 1.0, 0.5, 0.0, -0.1, -0.5]) + finally: + os.unlink(path) + + +def test_guess_tables_high_precision_floats(): + """High-precision floats (9 decimal places) round-trip within float64 tolerance.""" + values = [0.0, 0.123456789, 0.344343467, 0.577350269, + 0.707106781, 0.866025404, 0.951056516, 1.0] + content = " ".join(f"{v:.9f}" for v in values) + path = _write_tnt([('rfamp_hp', content)]) + try: + assert_allclose( + guess_tables(path)['rfamp_hp'], + values, + atol=1e-7, + ) + finally: + os.unlink(path) + + +def test_guess_tables_signed_hz_offsets(): + """Plain signed float Hz values (no SI suffix) are parsed correctly.""" + values = [0.0, 1234.567890123, -987.654321098, 500.123456789, -333.444555666] + content = " ".join(f"{v:.9f}" for v in values) + path = _write_tnt([('freq_offset', content)]) + try: + assert_allclose( + guess_tables(path)['freq_offset'], + values, + atol=1e-7, + ) + finally: + os.unlink(path) + + +def test_guess_tables_skips_scalars(): + """Single-token (scalar) entries are excluded from the result.""" + path = _write_tnt([ + ('rfamp', '1.0 5.0 10.0'), + ('tp', '50u'), + ('deadtime', '20u'), + ]) + try: + tables = guess_tables(path) + assert 'rfamp' in tables + assert 'tp' not in tables + assert 'deadtime' not in tables + finally: + os.unlink(path) + + +def test_guess_tables_name_filter(): + """The `names` argument restricts results to the requested subset.""" + path = _write_tnt([ + ('rfamp', '1.0 2.0 3.0'), + ('d1', '1m 2m 3m'), + ('d2', '1s 2s 3s'), + ]) + try: + tables = guess_tables(path, names=['rfamp', 'd1']) + assert set(tables.keys()) == {'rfamp', 'd1'} + finally: + os.unlink(path) + + +def test_guess_tables_convert_si_disabled(): + """When convert_si=False, values are returned as raw string token lists.""" + path = _write_tnt([('d1', '1m 2m 3m')]) + try: + tables = guess_tables(path, convert_si=False) + assert tables['d1'] == ['1m', '2m', '3m'] + finally: + os.unlink(path) + + +def test_guess_tables_no_pseq_section(): + """A file with no PSEQ section returns an empty dict without raising.""" + fd, path = tempfile.mkstemp(suffix='.tnt') + try: + with os.fdopen(fd, 'wb') as f: + f.write(b'TNT1.008' + b'\x00' * 256) + assert guess_tables(path) == {} + finally: + os.unlink(path) + + +def test_guess_tables_rejects_invalid_length_prefix(): + """A name with a mismatched length prefix is silently skipped.""" + fd, path = tempfile.mkstemp(suffix='.tnt') + try: + # Length says 99, only 5 bytes of name follow → invalid + bogus = struct.pack(' Date: Fri, 19 Jun 2026 12:48:17 -0300 Subject: [PATCH 24/38] Use temp path and skip test if data dir is not available --- tests/test_tecmag.py | 376 ++++++++++++++++--------------------------- 1 file changed, 139 insertions(+), 237 deletions(-) diff --git a/tests/test_tecmag.py b/tests/test_tecmag.py index 35589c74..8acf1b87 100644 --- a/tests/test_tecmag.py +++ b/tests/test_tecmag.py @@ -10,10 +10,8 @@ read(), guess_tables(), and the internal helpers _read_pascal_pair() and _convert_si() without requiring external test data. """ - import os import struct -import tempfile import numpy as np from numpy.testing import assert_allclose, assert_array_almost_equal @@ -39,7 +37,7 @@ def _pascal_string(s): def _tlv_section(tag, payload, flag=0): - """Build a TLV section record: 4-byte tag + 4-byte flag + 4-byte length + payload.""" + """Build a TLV record: 4-byte tag + 4-byte flag + 4-byte length + payload.""" tag_bytes = tag.encode('ascii') if isinstance(tag, str) else tag assert len(tag_bytes) == 4, "Tag must be exactly 4 bytes" return tag_bytes + struct.pack(' 0) and np.any(got < 0) - assert_allclose(got, [-0.5, -0.1, 0.0, 0.5, 1.0, 0.5, 0.0, -0.1, -0.5]) - finally: - os.unlink(path) + path = _write_tnt(tmp_path, [('shape_bip', '-0.5 -0.1 0.0 0.5 1.0 0.5 0.0 -0.1 -0.5')]) + got = guess_tables(path)['shape_bip'] + assert np.any(got > 0) and np.any(got < 0) + assert_allclose(got, [-0.5, -0.1, 0.0, 0.5, 1.0, 0.5, 0.0, -0.1, -0.5]) -def test_guess_tables_high_precision_floats(): +def test_guess_tables_high_precision_floats(tmp_path): """High-precision floats (9 decimal places) round-trip within float64 tolerance.""" values = [0.0, 0.123456789, 0.344343467, 0.577350269, 0.707106781, 0.866025404, 0.951056516, 1.0] content = " ".join(f"{v:.9f}" for v in values) - path = _write_tnt([('rfamp_hp', content)]) - try: - assert_allclose( - guess_tables(path)['rfamp_hp'], - values, - atol=1e-7, - ) - finally: - os.unlink(path) - - -def test_guess_tables_signed_hz_offsets(): + path = _write_tnt(tmp_path, [('rfamp_hp', content)]) + assert_allclose(guess_tables(path)['rfamp_hp'], values, atol=1e-7) + + +def test_guess_tables_signed_hz_offsets(tmp_path): """Plain signed float Hz values (no SI suffix) are parsed correctly.""" values = [0.0, 1234.567890123, -987.654321098, 500.123456789, -333.444555666] content = " ".join(f"{v:.9f}" for v in values) - path = _write_tnt([('freq_offset', content)]) - try: - assert_allclose( - guess_tables(path)['freq_offset'], - values, - atol=1e-7, - ) - finally: - os.unlink(path) - - -def test_guess_tables_skips_scalars(): + path = _write_tnt(tmp_path, [('freq_offset', content)]) + assert_allclose(guess_tables(path)['freq_offset'], values, atol=1e-7) + + +def test_guess_tables_skips_scalars(tmp_path): """Single-token (scalar) entries are excluded from the result.""" - path = _write_tnt([ + path = _write_tnt(tmp_path, [ ('rfamp', '1.0 5.0 10.0'), ('tp', '50u'), ('deadtime', '20u'), ]) - try: - tables = guess_tables(path) - assert 'rfamp' in tables - assert 'tp' not in tables - assert 'deadtime' not in tables - finally: - os.unlink(path) + tables = guess_tables(path) + assert 'rfamp' in tables + assert 'tp' not in tables + assert 'deadtime' not in tables -def test_guess_tables_name_filter(): +def test_guess_tables_name_filter(tmp_path): """The `names` argument restricts results to the requested subset.""" - path = _write_tnt([ + path = _write_tnt(tmp_path, [ ('rfamp', '1.0 2.0 3.0'), ('d1', '1m 2m 3m'), ('d2', '1s 2s 3s'), ]) - try: - tables = guess_tables(path, names=['rfamp', 'd1']) - assert set(tables.keys()) == {'rfamp', 'd1'} - finally: - os.unlink(path) + tables = guess_tables(path, names=['rfamp', 'd1']) + assert set(tables.keys()) == {'rfamp', 'd1'} -def test_guess_tables_convert_si_disabled(): +def test_guess_tables_convert_si_disabled(tmp_path): """When convert_si=False, values are returned as raw string token lists.""" - path = _write_tnt([('d1', '1m 2m 3m')]) - try: - tables = guess_tables(path, convert_si=False) - assert tables['d1'] == ['1m', '2m', '3m'] - finally: - os.unlink(path) + path = _write_tnt(tmp_path, [('d1', '1m 2m 3m')]) + assert guess_tables(path, convert_si=False)['d1'] == ['1m', '2m', '3m'] -def test_guess_tables_no_pseq_section(): +def test_guess_tables_no_pseq_section(tmp_path): """A file with no PSEQ section returns an empty dict without raising.""" - fd, path = tempfile.mkstemp(suffix='.tnt') - try: - with os.fdopen(fd, 'wb') as f: - f.write(b'TNT1.008' + b'\x00' * 256) - assert guess_tables(path) == {} - finally: - os.unlink(path) + path = tmp_path / "no_pseq.tnt" + path.write_bytes(b'TNT1.008' + b'\x00' * 256) + assert guess_tables(str(path)) == {} -def test_guess_tables_rejects_invalid_length_prefix(): +def test_guess_tables_rejects_invalid_length_prefix(tmp_path): """A name with a mismatched length prefix is silently skipped.""" - fd, path = tempfile.mkstemp(suffix='.tnt') - try: - # Length says 99, only 5 bytes of name follow → invalid - bogus = struct.pack(' Date: Sun, 16 Aug 2026 17:48:03 -0500 Subject: [PATCH 25/38] REL: 0.12 --- CHANGELOG.md | 13 ++++++++++++- nmrglue/__init__.py | 2 +- setup.py | 2 +- 3 files changed, 14 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8969d2a6..7684e327 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ -0.12 (Unreleased) +0.12 (2026-08-16) ================= +* Bruker to csdm conversion returns the correct dimension units (#227) +* When reading bruker files both utf-8 and cp1252 encodings are tried (#239) +* Correctly read Bruker data in floating point format (#249) +* Support reading partially acquired Bruker data (#252) +* Add support for reading JEOL jdf format (#228) +* Add support for reading RS2D/QUAD format (#230) +* Add the update_uc function to update a unit_conversion object (#250) +* Add a guess_tables function for use when reading Tecmat sequence tables (#257) +* The read function in nmrml can optionally return processed data (#232) +* Add support for Numpy 2.0+ (#224, #248) +* Various minor fixes (#233, #236, #242, #238, #246, #245, #251) 0.11 (2024-10-29) diff --git a/nmrglue/__init__.py b/nmrglue/__init__.py index 88c7a986..07c23d06 100644 --- a/nmrglue/__init__.py +++ b/nmrglue/__init__.py @@ -3,4 +3,4 @@ from .util import * from .analysis import * -__version__ = '0.12-dev' +__version__ = '0.12' diff --git a/setup.py b/setup.py index 5666a22b..9161a70e 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ setup( name='nmrglue', - version='0.12-dev', # change this in nmrglue/__init__.py also + version='0.12', # change this in nmrglue/__init__.py also description='A module for working with NMR data in Python', long_description=long_description, url='http://www.nmrglue.com', From 7af35af2387726a4974b69330fddbb41bffdbe96 Mon Sep 17 00:00:00 2001 From: Jonathan Helmus Date: Sun, 16 Aug 2026 17:56:30 -0500 Subject: [PATCH 26/38] start v0.13-dev --- CHANGELOG.md | 4 ++++ nmrglue/__init__.py | 2 +- setup.py | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7684e327..f1d1a7e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,7 @@ +0.13 (Unreleased) +================= + + 0.12 (2026-08-16) ================= * Bruker to csdm conversion returns the correct dimension units (#227) diff --git a/nmrglue/__init__.py b/nmrglue/__init__.py index 07c23d06..8f141ce7 100644 --- a/nmrglue/__init__.py +++ b/nmrglue/__init__.py @@ -3,4 +3,4 @@ from .util import * from .analysis import * -__version__ = '0.12' +__version__ = '0.13-dev' diff --git a/setup.py b/setup.py index 9161a70e..955abe55 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ setup( name='nmrglue', - version='0.12', # change this in nmrglue/__init__.py also + version='0.13-dev', # change this in nmrglue/__init__.py also description='A module for working with NMR data in Python', long_description=long_description, url='http://www.nmrglue.com', From 65e83e4039250d9c6a74a4f5f606c7e5e82bdac5 Mon Sep 17 00:00:00 2001 From: Jonathan Helmus Date: Sun, 16 Aug 2026 17:58:42 -0500 Subject: [PATCH 27/38] drop support for Python 3.8,9 add 3.13,14 --- .github/workflows/ci.yml | 2 +- setup.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index de029f44..efd34178 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,7 +11,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] fail-fast: false steps: diff --git a/setup.py b/setup.py index 955abe55..e657b8f9 100644 --- a/setup.py +++ b/setup.py @@ -26,11 +26,11 @@ 'Intended Audience :: Developers', 'License :: OSI Approved :: BSD License', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.8', - 'Programming Language :: Python :: 3.9', 'Programming Language :: Python :: 3.10', 'Programming Language :: Python :: 3.11', 'Programming Language :: Python :: 3.12', + 'Programming Language :: Python :: 3.13', + 'Programming Language :: Python :: 3.14', 'Topic :: Scientific/Engineering', 'Operating System :: MacOS :: MacOS X', 'Operating System :: Microsoft :: Windows', From e6e47f4bb859021e8f7d8728314f401155ecebca Mon Sep 17 00:00:00 2001 From: Christian Fernandez Date: Mon, 21 Sep 2026 19:05:22 +0200 Subject: [PATCH 28/38] Document GM Gaussian constant --- nmrglue/process/proc_base.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/nmrglue/process/proc_base.py b/nmrglue/process/proc_base.py index d353a207..e55a4320 100644 --- a/nmrglue/process/proc_base.py +++ b/nmrglue/process/proc_base.py @@ -93,6 +93,18 @@ def gm(data, g1=0.0, g2=0.0, g3=0.0, inv=False, rev=False): ndata : ndarray Array of NMR data with apodization applied. + Notes + ----- + The factor ``0.6`` in the Gaussian term matches the NMRPipe GM function. + It approximates + + .. math:: + \\frac{1}{2\\sqrt{\\ln(2)}} \\approx 0.600561, + + the coefficient for which the Fourier transform of the Gaussian + contribution has a full width at half maximum of ``g2``. The rounded + value is retained for numerical compatibility with NMRPipe. + """ size = data.shape[-1] e = pi * np.arange(size) * g1 From 3505954195d5380ab5cdba6663c7dcc892a57e2d Mon Sep 17 00:00:00 2001 From: Christian Fernandez Date: Mon, 21 Sep 2026 19:34:32 +0200 Subject: [PATCH 29/38] Fix Bruker magnitude quadrature metadata --- nmrglue/fileio/bruker.py | 4 ++ nmrglue/fileio/tests/test_bruker.py | 67 ++++++++++++++++++++++++++++- 2 files changed, 69 insertions(+), 2 deletions(-) diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index c2de4573..f778412a 100644 --- a/nmrglue/fileio/bruker.py +++ b/nmrglue/fileio/bruker.py @@ -182,8 +182,10 @@ def add_axis_to_udic(udic, dic, udim, strip_fake): udic[udim]["encoding"] = "undefined" elif aq_mod == 1: udic[udim]["encoding"] = "magnitude" # qf + udic[udim]["complex"] = False elif aq_mod == 2: udic[udim]["encoding"] = "magnitude" # qsec + udic[udim]["complex"] = False elif aq_mod == 3: udic[udim]["encoding"] = "tppi" elif aq_mod == 4: @@ -203,8 +205,10 @@ def add_axis_to_udic(udic, dic, udim, strip_fake): aq_mod = dic[pro_file]["MC2"] if aq_mod == 0: udic[udim]["encoding"] = "magnitude" # qf + udic[udim]["complex"] = False elif aq_mod == 1: udic[udim]["encoding"] = "magnitude" # qsec + udic[udim]["complex"] = False elif aq_mod == 2: udic[udim]["encoding"] = "tppi" elif aq_mod == 3: diff --git a/nmrglue/fileio/tests/test_bruker.py b/nmrglue/fileio/tests/test_bruker.py index 55a04f5c..b29f97bb 100644 --- a/nmrglue/fileio/tests/test_bruker.py +++ b/nmrglue/fileio/tests/test_bruker.py @@ -1,15 +1,78 @@ """ Unit tests for nmrglue/fileio/bruker.py module """ import os +import shutil +import tempfile import numpy as np import nmrglue as ng -import shutil -import tempfile +import pytest # Test data. DATA_DIR = os.path.join(os.path.dirname(__file__), 'bruker_test_data') + +def _make_2d_dic(indirect_dic): + """Return minimal Bruker parameters for a 2D dataset.""" + return { + 'acqus': { + 'AQ_mod': 3, + 'NUC1': '1H', + 'O1': 0.0, + 'SFO1': 400.0, + 'SW_h': 10000.0, + }, + 'procs': {'SF': 400.0}, + **indirect_dic, + } + + +@pytest.mark.parametrize('fnmode', [1, 2]) +def test_guess_udic_magnitude_indirect_dimension(fnmode): + """QF and QSEC indirect dimensions are real, not complex.""" + dic = _make_2d_dic({ + 'acqu2s': { + 'FnMODE': fnmode, + 'NUC1': '1H', + 'SFO1': 400.0, + 'SW': 10.0, + }, + 'proc2s': {'SF': 400.0}, + }) + data = np.zeros((40, 128), dtype=np.complex128) + + udic = ng.bruker.guess_udic(dic, data) + converter = ng.convert.converter() + converter.from_bruker(dic, data) + pipe_dic, pipe_data = converter.to_pipe() + _, transposed = ng.pipe_proc.tp(pipe_dic.copy(), pipe_data) + + assert udic[0]['encoding'] == 'magnitude' + assert udic[0]['complex'] is False + assert pipe_dic['FDF2QUADFLAG'] == 0.0 + assert pipe_dic['FDF1QUADFLAG'] == 1.0 + assert transposed.shape == (128, 40) + + +@pytest.mark.parametrize('mc2', [0, 1]) +def test_guess_udic_magnitude_indirect_dimension_from_proc_params(mc2): + """Processed QF and QSEC indirect dimensions are real, not complex.""" + dic = _make_2d_dic({ + 'proc2s': { + 'AXNUC': '1H', + 'MC2': mc2, + 'OFFSET': 0.0, + 'SF': 400.0, + 'SW_p': 4000.0, + }, + }) + data = np.zeros((40, 128), dtype=np.complex128) + + udic = ng.bruker.guess_udic(dic, data) + + assert udic[0]['encoding'] == 'magnitude' + assert udic[0]['complex'] is False + def test_read_pdata(): """Reading processed bruker 1D data""" From 1a1d2dc88f9a9b407959df5ef4e93c41b0b16d38 Mon Sep 17 00:00:00 2001 From: Christian Fernandez Date: Thu, 1 Oct 2026 00:08:29 +0200 Subject: [PATCH 30/38] Replace removed pinv2 in linear prediction --- nmrglue/process/proc_lp.py | 4 ++-- tests/test_proc_lp.py | 21 +++++++++++++++++++++ 2 files changed, 23 insertions(+), 2 deletions(-) create mode 100644 tests/test_proc_lp.py diff --git a/nmrglue/process/proc_lp.py b/nmrglue/process/proc_lp.py index 05012eaa..e309db7d 100644 --- a/nmrglue/process/proc_lp.py +++ b/nmrglue/process/proc_lp.py @@ -838,13 +838,13 @@ def find_lpc_qr(D, d): q, r = np.asmatrix(q), np.asmatrix(r) # SPEED - # the next line is slow and the use of pinv2 should be avoided as + # the next line is slow and the use of a pseudo-inverse should be avoided as # pseudo inversion of r involves a computationally expensive SVD # decomposition which is not needed. Rather r*x = q.H*d should be # solved for x using LAPACK's ZTRTRS function (or similar function with # different prefix). This is not currently available in scipy/numpy and # therefore is not used here. - return scipy.linalg.pinv2(r) * q.H * d + return scipy.linalg.pinv(r) * q.H * d def find_lpc_cholesky(D, d): diff --git a/tests/test_proc_lp.py b/tests/test_proc_lp.py new file mode 100644 index 00000000..b6ea0a56 --- /dev/null +++ b/tests/test_proc_lp.py @@ -0,0 +1,21 @@ +import numpy as np + +import nmrglue as ng + + +def test_find_lpc_qr(): + D = np.array( + [ + [1.0, 0.0], + [0.0, 1.0], + [1.0, 1.0], + ], + dtype=np.float64, + ) + d = np.array([[1.0], [2.0], [3.0]], dtype=np.float64) + + result = np.asarray(ng.proc_lp.find_lpc_qr(D, d)) + + assert result.shape == (2, 1) + assert result.dtype == np.dtype("float64") + np.testing.assert_allclose(result, [[1.0], [2.0]], rtol=1e-14, atol=1e-14) From 61a45e6639522c887643d9545cabe40b490dc32a Mon Sep 17 00:00:00 2001 From: "Daniel C. Pizetta" Date: Fri, 19 Jun 2026 14:27:26 -0300 Subject: [PATCH 31/38] Add optional bounds parameters to autophase --- nmrglue/process/proc_autophase.py | 83 ++++++++-- tests/test_autophase.py | 248 ++++++++++++++++++++++++++++++ 2 files changed, 314 insertions(+), 17 deletions(-) create mode 100644 tests/test_autophase.py diff --git a/nmrglue/process/proc_autophase.py b/nmrglue/process/proc_autophase.py index 7d5479cc..9f7b5957 100644 --- a/nmrglue/process/proc_autophase.py +++ b/nmrglue/process/proc_autophase.py @@ -13,7 +13,8 @@ from .proc_base import ps -def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, **kwargs): +def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, + bounds=None, **kwargs): """ Automatic linear phase correction @@ -29,29 +30,65 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, **kwar p1 : float Initial first order phase in degrees. peak_width : int - Width of the ROI for peak_minima optimization (number of surrounding indexes) - return_phases : Bool - returns a list of optimized values of phases [p0, p1] in addition - to the phased data - kwargs : additional key-word arguments to be passed to the solver - The are documented at the following link: + Width of the ROI for peak_minima optimization (number of surrounding + indexes). Only used when fn='peak_minima'. + return_phases : bool + If True, also return the optimised ``[p0, p1]`` list in addition to + the phased data. + bounds : sequence of two (min, max) pairs or None, optional + Bounds on the zero- and first-order phases, respectively, in degrees. + For example ``bounds=[(-180, 180), (-1800, 1800)]`` restricts p0 to + +-180 and p1 to +-1800 degrees. When *bounds* is provided the + optimisation uses :func:`scipy.optimize.minimize` with the Nelder-Mead + method (same simplex algorithm as :func:`scipy.optimize.fmin`) so the + solver behaviour is equivalent but constrained. When *bounds* is + ``None`` (default) the original unconstrained + :func:`scipy.optimize.fmin` is used and the function is fully + backward-compatible. + + Bounds are useful to find only the zero-order phase (set p1 bounds to + ``(0, 0)``), or to prevent the optimiser from wandering into physically + unreasonable regions when the spectrum has poor initial phase. + kwargs : additional key-word arguments + Passed directly to the underlying solver. + + When *bounds* is ``None`` (fmin): see https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin.html - Some of the more useful ones for this use case: + Useful options: ``disp`` (bool), ``ftol`` (float). - * disp : Bool - Turns on or off the printing of convergence messages - By default, this is set to True. - * ftol : float - Absolute error in fn between iterations that is acceptable for - convergence. The default is 1e-4. + When *bounds* is provided (minimize): see + https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html + The method is fixed to 'Nelder-Mead'; pass e.g. + ``options={'xatol': 1e-4}``. Returns ------- ndata : ndarray Phased NMR data. + opt : list of float, optional + ``[p0, p1]`` optimised phase values in degrees. Only returned when + *return_phases* is ``True``. + + Examples + -------- + Unconstrained autophase (original behaviour): + + >>> phased = ng.proc_autophase.autops(data, 'acme') + + Bounded autophase — restrict search to +-180 in p0, +-1800 in p1: + + >>> phased, phases = ng.proc_autophase.autops( + ... data, 'acme', p0=0.0, p1=0.0, + ... bounds=[(-180, 180), (-1800, 1800)], + ... return_phases=True) + + Zero-order only — fix p1 to zero: + + >>> phased = ng.proc_autophase.autops( + ... data, 'acme', bounds=[(-180, 180), (0, 0)]) """ - arguments = [data,] + arguments = [data] if not callable(fn): if fn == 'peak_minima': @@ -63,11 +100,23 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, **kwar 'acme': _ps_acme_score, }[fn] except KeyError: - raise KeyError(f'Unable to find algorithm "{fn}". Use "acme", "peak_minima" or pass a custom function.') + raise KeyError( + f'Unable to find algorithm "{fn}". ' + 'Use "acme", "peak_minima" or pass a custom function.' + ) opt = [p0, p1] - opt = scipy.optimize.fmin(fn, x0=opt, args=tuple(arguments), **kwargs) + if bounds is not None: + result = scipy.optimize.minimize( + fn, x0=opt, args=tuple(arguments), + method='Nelder-Mead', + bounds=bounds, + **kwargs, + ) + opt = result.x + else: + opt = scipy.optimize.fmin(fn, x0=opt, args=tuple(arguments), **kwargs) phasedspc = ps(data, p0=opt[0], p1=opt[1]) diff --git a/tests/test_autophase.py b/tests/test_autophase.py new file mode 100644 index 00000000..cd991d3a --- /dev/null +++ b/tests/test_autophase.py @@ -0,0 +1,248 @@ +""" +tests/test_autophase.py — Tests for nmrglue.process.proc_autophase + +Covers autops() with and without bounds, both built-in scoring functions, +callable scoring functions, return_phases flag, and error handling. + +The bounded path uses scipy.optimize.minimize (Nelder-Mead); the unbounded +path uses scipy.optimize.fmin. Both are exercised here. +""" +import numpy as np +from numpy.testing import assert_allclose +import pytest + +import nmrglue as ng +from nmrglue.process.proc_autophase import ( + autops, + _ps_acme_score, + _ps_peak_minima_score, +) +from nmrglue.process.proc_base import ps + + +# ============================================================================= +# Helpers +# ============================================================================= + +def _lorentzian_spectrum(n=512, sw=10000.0, f0=1000.0, t2=0.1, p0=0.0, p1=0.0): + """Return a phase-shifted Lorentzian spectrum (complex ndarray). + + Parameters + ---------- + n : int — number of points + sw : float — spectral width (Hz) + f0 : float — peak frequency (Hz) + t2 : float — transverse relaxation time (s) + p0 : float — zero-order phase error (degrees) + p1 : float — first-order phase error (degrees) + """ + t = np.arange(n) / sw + fid = np.exp(1j * 2 * np.pi * f0 * t) * np.exp(-t / t2) + spec = np.fft.fftshift(np.fft.fft(fid)) + return ps(spec, p0=p0, p1=p1) + + +# ============================================================================= +# Backward-compatibility — unbounded (original fmin path) +# ============================================================================= + +def test_autops_acme_returns_ndarray(): + """autops with 'acme' returns an ndarray of the same shape.""" + data = _lorentzian_spectrum(p0=30.0) + result = autops(data, 'acme', disp=False) + assert isinstance(result, np.ndarray) + assert result.shape == data.shape + + +def test_autops_peak_minima_returns_ndarray(): + """autops with 'peak_minima' returns an ndarray of the same shape.""" + data = _lorentzian_spectrum(p0=30.0) + result = autops(data, 'peak_minima', disp=False) + assert isinstance(result, np.ndarray) + assert result.shape == data.shape + + +def test_autops_return_phases_unbounded(): + """return_phases=True returns (ndarray, array-like) with two elements.""" + data = _lorentzian_spectrum(p0=45.0) + result = autops(data, 'acme', return_phases=True, disp=False) + assert isinstance(result, tuple) and len(result) == 2 + phased, opt = result + assert isinstance(phased, np.ndarray) + assert len(opt) == 2 + + +def test_autops_callable_fn_unbounded(): + """A callable scoring function is accepted in place of a string alias.""" + data = _lorentzian_spectrum(p0=20.0) + result = autops(data, _ps_acme_score, disp=False) + assert isinstance(result, np.ndarray) + + +def test_autops_unknown_fn_raises(): + """An unrecognised string fn raises KeyError with a helpful message.""" + data = _lorentzian_spectrum() + with pytest.raises(KeyError, match='Unable to find algorithm'): + autops(data, 'not_a_real_algorithm') + + +def test_autops_acme_reduces_phase_error(): + """autops 'acme' reduces a known phase error towards zero.""" + p0_true = 60.0 + data = _lorentzian_spectrum(p0=p0_true) + phased, opt = autops(data, 'acme', return_phases=True, disp=False) + # After correction the residual phase should be much smaller than the + # original error (we don't require exact recovery — ACME is heuristic). + assert abs(opt[0]) < abs(p0_true) + + +# ============================================================================= +# Bounded path — scipy.optimize.minimize(Nelder-Mead) +# ============================================================================= + +def test_autops_bounded_returns_ndarray(): + """autops with bounds= returns an ndarray of the same shape.""" + data = _lorentzian_spectrum(p0=30.0) + result = autops(data, 'acme', bounds=[(-180, 180), (-1800, 1800)]) + assert isinstance(result, np.ndarray) + assert result.shape == data.shape + + +def test_autops_bounded_return_phases(): + """return_phases=True with bounds returns (ndarray, array-like).""" + data = _lorentzian_spectrum(p0=45.0) + phased, opt = autops(data, 'acme', + bounds=[(-180, 180), (-1800, 1800)], + return_phases=True) + assert isinstance(phased, np.ndarray) + assert len(opt) == 2 + + +def test_autops_p1_fixed_to_zero(): + """Bounds (0, 0) on p1 keep the first-order phase exactly zero.""" + data = _lorentzian_spectrum(p0=50.0, p1=0.0) + _, opt = autops(data, 'acme', + bounds=[(-180, 180), (0, 0)], + return_phases=True) + assert abs(opt[1]) < 1e-6, f"p1 should be 0.0, got {opt[1]}" + + +def test_autops_p0_fixed_to_zero(): + """Bounds (0, 0) on p0 keep the zero-order phase exactly zero.""" + data = _lorentzian_spectrum(p0=0.0, p1=200.0) + _, opt = autops(data, 'acme', + bounds=[(0, 0), (-1800, 1800)], + return_phases=True) + assert abs(opt[0]) < 1e-6, f"p0 should be 0.0, got {opt[0]}" + + +def test_autops_bounded_respects_p0_range(): + """Optimised p0 stays within the specified bounds.""" + data = _lorentzian_spectrum(p0=30.0) + b0 = (-90, 90) + _, opt = autops(data, 'acme', + bounds=[b0, (-1800, 1800)], + return_phases=True) + assert b0[0] <= opt[0] <= b0[1], f"p0={opt[0]} outside {b0}" + + +def test_autops_bounded_respects_p1_range(): + """Optimised p1 stays within the specified bounds.""" + data = _lorentzian_spectrum(p0=0.0, p1=300.0) + b1 = (-500, 500) + _, opt = autops(data, 'acme', + bounds=[(-180, 180), b1], + return_phases=True) + assert b1[0] <= opt[1] <= b1[1], f"p1={opt[1]} outside {b1}" + + +def test_autops_bounded_acme_reduces_phase_error(): + """Bounded autops 'acme' reduces a known phase error towards zero.""" + p0_true = 60.0 + data = _lorentzian_spectrum(p0=p0_true) + _, opt = autops(data, 'acme', + bounds=[(-180, 180), (-1800, 1800)], + return_phases=True) + assert abs(opt[0]) < abs(p0_true) + + +def test_autops_bounded_peak_minima(): + """Bounded autops works with 'peak_minima' scoring function.""" + data = _lorentzian_spectrum(p0=40.0) + result = autops(data, 'peak_minima', + bounds=[(-180, 180), (-1800, 1800)]) + assert isinstance(result, np.ndarray) + assert result.shape == data.shape + + +def test_autops_bounded_callable_fn(): + """Bounded autops accepts a callable scoring function.""" + data = _lorentzian_spectrum(p0=25.0) + phased, opt = autops(data, _ps_acme_score, + bounds=[(-180, 180), (-1800, 1800)], + return_phases=True) + assert isinstance(phased, np.ndarray) + assert len(opt) == 2 + + +def test_autops_bounded_unknown_fn_raises(): + """An unrecognised string fn raises KeyError even when bounds are set.""" + data = _lorentzian_spectrum() + with pytest.raises(KeyError, match='Unable to find algorithm'): + autops(data, 'not_a_real_algorithm', + bounds=[(-180, 180), (-1800, 1800)]) + + +def test_autops_bounded_kwargs_passed_to_minimize(): + """Extra kwargs are forwarded to scipy.optimize.minimize.""" + data = _lorentzian_spectrum(p0=30.0) + # options= is a minimize-specific kwarg; should not raise + phased, opt = autops(data, 'acme', + bounds=[(-180, 180), (-1800, 1800)], + return_phases=True, + options={'xatol': 1e-2, 'fatol': 1e-2, 'disp': False}) + assert isinstance(phased, np.ndarray) + + +def test_autops_bounded_vs_unbounded_agreement(): + """Bounded and unbounded results are consistent on a well-conditioned spectrum.""" + data = _lorentzian_spectrum(p0=30.0, p1=0.0) + _, opt_ub = autops(data, 'acme', return_phases=True, disp=False) + _, opt_b = autops(data, 'acme', + bounds=[(-180, 180), (-1800, 1800)], + return_phases=True) + # Both should find similar p0 (within solver tolerance + ACME heuristic noise) + assert abs(opt_ub[0] - opt_b[0]) < 10.0, ( + f"Unbounded p0={opt_ub[0]:.2f} and bounded p0={opt_b[0]:.2f} disagree") + + +# ============================================================================= +# Score function unit tests +# ============================================================================= + +def test_acme_score_returns_scalar(): + """_ps_acme_score returns a finite scalar.""" + data = _lorentzian_spectrum() + score = _ps_acme_score((0.0, 0.0), data) + assert np.isfinite(score) + assert np.isscalar(score) or score.ndim == 0 + + +def test_acme_score_lower_for_phased(): + """_ps_acme_score is lower for correctly phased than misphased data.""" + data = _lorentzian_spectrum() + score_good = _ps_acme_score((0.0, 0.0), data) + score_bad = _ps_acme_score((90.0, 0.0), data) + assert score_good < score_bad + + +def test_peak_minima_score_returns_scalar(): + """_ps_peak_minima_score returns a finite scalar.""" + data = _lorentzian_spectrum() + score = _ps_peak_minima_score((0.0, 0.0), data, peak_width=50) + assert np.isfinite(score) + + +if __name__ == '__main__': + import pytest as _pytest + _pytest.main([__file__, '-v']) From c5aa62661ee538312a273fc00841c9f71f7a6bbe Mon Sep 17 00:00:00 2001 From: "Daniel C. Pizetta" Date: Wed, 7 Oct 2026 08:33:37 -0300 Subject: [PATCH 32/38] Address review: stricter autophase tests, accept fmin kwargs with bounds - Use t2=0.01 and halve the first FID point in the synthetic spectrum so it is not truncated and needs no intrinsic phase correction. - Check that the recovered phases undo the applied error (opt[0] ~ -p0, p1 ~ 0) instead of only comparing magnitudes. - Translate fmin-style kwargs (xtol, ftol, maxiter, maxfun, disp, initial_simplex) into minimize options when bounds are given. Co-Authored-By: Claude Opus 5.5 --- nmrglue/process/proc_autophase.py | 24 +++++++++++++++++++++++- tests/test_autophase.py | 27 ++++++++++++++++++++++----- 2 files changed, 45 insertions(+), 6 deletions(-) diff --git a/nmrglue/process/proc_autophase.py b/nmrglue/process/proc_autophase.py index 9f7b5957..0085e19e 100644 --- a/nmrglue/process/proc_autophase.py +++ b/nmrglue/process/proc_autophase.py @@ -12,6 +12,17 @@ from .proc_base import ps +# fmin keyword arguments and their scipy.optimize.minimize(Nelder-Mead) +# option equivalents, used when autops is called with bounds. +_FMIN_TO_MINIMIZE_OPTIONS = { + 'xtol': 'xatol', + 'ftol': 'fatol', + 'maxiter': 'maxiter', + 'maxfun': 'maxfev', + 'disp': 'disp', + 'initial_simplex': 'initial_simplex', +} + def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, bounds=None, **kwargs): @@ -59,7 +70,10 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, When *bounds* is provided (minimize): see https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html The method is fixed to 'Nelder-Mead'; pass e.g. - ``options={'xatol': 1e-4}``. + ``options={'xatol': 1e-4}``. The fmin-style arguments ``xtol``, + ``ftol``, ``maxiter``, ``maxfun``, ``disp`` and ``initial_simplex`` + are also accepted and translated to the matching ``options`` entries + (explicit ``options`` take precedence). Returns ------- @@ -108,10 +122,18 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, opt = [p0, p1] if bounds is not None: + # Translate fmin-style keyword arguments into minimize options so the + # same call works with and without bounds. + options = {} + for fmin_key, option_key in _FMIN_TO_MINIMIZE_OPTIONS.items(): + if fmin_key in kwargs: + options[option_key] = kwargs.pop(fmin_key) + options.update(kwargs.pop('options', {})) result = scipy.optimize.minimize( fn, x0=opt, args=tuple(arguments), method='Nelder-Mead', bounds=bounds, + options=options, **kwargs, ) opt = result.x diff --git a/tests/test_autophase.py b/tests/test_autophase.py index cd991d3a..0ab63888 100644 --- a/tests/test_autophase.py +++ b/tests/test_autophase.py @@ -19,12 +19,15 @@ ) from nmrglue.process.proc_base import ps +# Tolerance (degrees) when comparing recovered phases with the applied error. +PHASE_TOL = 5.0 + # ============================================================================= # Helpers # ============================================================================= -def _lorentzian_spectrum(n=512, sw=10000.0, f0=1000.0, t2=0.1, p0=0.0, p1=0.0): +def _lorentzian_spectrum(n=512, sw=10000.0, f0=1000.0, t2=0.01, p0=0.0, p1=0.0): """Return a phase-shifted Lorentzian spectrum (complex ndarray). Parameters @@ -38,6 +41,9 @@ def _lorentzian_spectrum(n=512, sw=10000.0, f0=1000.0, t2=0.1, p0=0.0, p1=0.0): """ t = np.arange(n) / sw fid = np.exp(1j * 2 * np.pi * f0 * t) * np.exp(-t / t2) + # Halve the first point so the DFT does not add a baseline offset that + # ACME would compensate with a spurious p1. + fid[0] *= 0.5 spec = np.fft.fftshift(np.fft.fft(fid)) return ps(spec, p0=p0, p1=p1) @@ -91,9 +97,9 @@ def test_autops_acme_reduces_phase_error(): p0_true = 60.0 data = _lorentzian_spectrum(p0=p0_true) phased, opt = autops(data, 'acme', return_phases=True, disp=False) - # After correction the residual phase should be much smaller than the - # original error (we don't require exact recovery — ACME is heuristic). - assert abs(opt[0]) < abs(p0_true) + # The correction must undo the applied error: opt[0] ~ -p0_true, p1 ~ 0. + assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" + assert abs(opt[1]) < PHASE_TOL, f"p1={opt[1]}" # ============================================================================= @@ -163,7 +169,8 @@ def test_autops_bounded_acme_reduces_phase_error(): _, opt = autops(data, 'acme', bounds=[(-180, 180), (-1800, 1800)], return_phases=True) - assert abs(opt[0]) < abs(p0_true) + assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" + assert abs(opt[1]) < PHASE_TOL, f"p1={opt[1]}" def test_autops_bounded_peak_minima(): @@ -204,6 +211,16 @@ def test_autops_bounded_kwargs_passed_to_minimize(): assert isinstance(phased, np.ndarray) +def test_autops_bounded_accepts_fmin_kwargs(): + """fmin-style kwargs (disp, xtol, ftol, ...) also work with bounds.""" + data = _lorentzian_spectrum(p0=30.0) + _, opt = autops(data, 'acme', + bounds=[(-180, 180), (-1800, 1800)], + return_phases=True, + disp=False, xtol=1e-4, ftol=1e-4, maxiter=2000) + assert abs(30.0 + opt[0]) < PHASE_TOL, f"p0={opt[0]}" + + def test_autops_bounded_vs_unbounded_agreement(): """Bounded and unbounded results are consistent on a well-conditioned spectrum.""" data = _lorentzian_spectrum(p0=30.0, p1=0.0) From b0514c44870cb0c1aef19efcaca7c6461570dccc Mon Sep 17 00:00:00 2001 From: "Daniel C. Pizetta" Date: Wed, 7 Oct 2026 08:35:00 -0300 Subject: [PATCH 33/38] Validate autops bounds shape and fix documented return type scipy silently broadcasts a single (min, max) pair to both phases, so bounds=[(-180, 180)] also constrained p1. Require exactly two pairs. Co-Authored-By: Claude Opus 5.5 --- nmrglue/process/proc_autophase.py | 7 ++++++- tests/test_autophase.py | 12 ++++++++++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/nmrglue/process/proc_autophase.py b/nmrglue/process/proc_autophase.py index 0085e19e..3662a9f5 100644 --- a/nmrglue/process/proc_autophase.py +++ b/nmrglue/process/proc_autophase.py @@ -79,7 +79,7 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, ------- ndata : ndarray Phased NMR data. - opt : list of float, optional + opt : ndarray, optional ``[p0, p1]`` optimised phase values in degrees. Only returned when *return_phases* is ``True``. @@ -122,6 +122,11 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, opt = [p0, p1] if bounds is not None: + if len(bounds) != 2 or any(len(b) != 2 for b in bounds): + raise ValueError( + 'bounds must contain two (min, max) pairs, one for p0 and ' + f'one for p1, got {bounds!r}' + ) # Translate fmin-style keyword arguments into minimize options so the # same call works with and without bounds. options = {} diff --git a/tests/test_autophase.py b/tests/test_autophase.py index 0ab63888..3393ad72 100644 --- a/tests/test_autophase.py +++ b/tests/test_autophase.py @@ -200,6 +200,18 @@ def test_autops_bounded_unknown_fn_raises(): bounds=[(-180, 180), (-1800, 1800)]) +@pytest.mark.parametrize('bounds', [ + [(-180, 180)], + [(-180, 180), (-1800, 1800), (0, 0)], + [(-180, 180), (0,)], +]) +def test_autops_invalid_bounds_raises(bounds): + """bounds must be exactly two (min, max) pairs.""" + data = _lorentzian_spectrum(p0=30.0) + with pytest.raises(ValueError, match='two \\(min, max\\) pairs'): + autops(data, 'acme', bounds=bounds) + + def test_autops_bounded_kwargs_passed_to_minimize(): """Extra kwargs are forwarded to scipy.optimize.minimize.""" data = _lorentzian_spectrum(p0=30.0) From 9c6803dd6973e006e087141bde318c99593bffbb Mon Sep 17 00:00:00 2001 From: "Daniel C. Pizetta" Date: Wed, 7 Oct 2026 08:41:53 -0300 Subject: [PATCH 34/38] Document single-phase bounds and the p0 limit caveat in autops - Explain (None, None) to leave one phase unbounded, with examples. - Note that hard p0 limits can trap the optimiser near +-180 since the phase is periodic; the zero-order example now leaves p0 unbounded. - Test single-phase bounds and p0 recovery near 180 with p0 unbounded. Co-Authored-By: Claude Opus 5.5 --- nmrglue/process/proc_autophase.py | 29 ++++++++++++++++++++++++----- tests/test_autophase.py | 24 ++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 5 deletions(-) diff --git a/nmrglue/process/proc_autophase.py b/nmrglue/process/proc_autophase.py index 3662a9f5..70496126 100644 --- a/nmrglue/process/proc_autophase.py +++ b/nmrglue/process/proc_autophase.py @@ -57,9 +57,23 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, :func:`scipy.optimize.fmin` is used and the function is fully backward-compatible. - Bounds are useful to find only the zero-order phase (set p1 bounds to - ``(0, 0)``), or to prevent the optimiser from wandering into physically - unreasonable regions when the spectrum has poor initial phase. + Use ``(None, None)`` to leave one of the phases unbounded, e.g. + ``bounds=[(None, None), (-90, 90)]`` bounds only p1 and + ``bounds=[(-180, 180), (None, None)]`` bounds only p0. A pair with + equal limits fixes that phase, e.g. ``(0, 0)`` for p1 to find only + the zero-order phase. + + Bounds are useful to find only the zero-order phase, or to prevent + the optimiser from wandering into physically unreasonable regions + when the spectrum has poor initial phase. + + .. note:: The bounds are hard limits, while the zero-order phase is + periodic (360 degrees). If the required p0 correction lies close + to a p0 limit, the optimiser can stop at that limit and compensate + with a wrong p1 instead of reaching the correct value. When only + p1 needs to be restricted, leave p0 unbounded with + ``(None, None)``; otherwise keep a margin between the p0 limits + and the expected correction. kwargs : additional key-word arguments Passed directly to the underlying solver. @@ -96,10 +110,15 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, ... bounds=[(-180, 180), (-1800, 1800)], ... return_phases=True) - Zero-order only — fix p1 to zero: + Zero-order only — fix p1 to zero and leave p0 unbounded: >>> phased = ng.proc_autophase.autops( - ... data, 'acme', bounds=[(-180, 180), (0, 0)]) + ... data, 'acme', bounds=[(None, None), (0, 0)]) + + Bound only the first-order phase: + + >>> phased = ng.proc_autophase.autops( + ... data, 'acme', bounds=[(None, None), (-90, 90)]) """ arguments = [data] diff --git a/tests/test_autophase.py b/tests/test_autophase.py index 3393ad72..3c834095 100644 --- a/tests/test_autophase.py +++ b/tests/test_autophase.py @@ -200,6 +200,30 @@ def test_autops_bounded_unknown_fn_raises(): bounds=[(-180, 180), (-1800, 1800)]) +@pytest.mark.parametrize('bounds', [ + [(-180, 180), (None, None)], + [(None, None), (-90, 90)], + [(-np.inf, np.inf), (-1800, 1800)], +]) +def test_autops_bounds_single_phase(bounds): + """(None, None) or infinite limits leave one phase unbounded.""" + p0_true, p1_true = 60.0, 30.0 + data = _lorentzian_spectrum(p0=p0_true, p1=p1_true) + _, opt = autops(data, 'acme', bounds=bounds, return_phases=True) + assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" + assert abs(p1_true + opt[1]) < PHASE_TOL, f"p1={opt[1]}" + + +def test_autops_p0_only_unbounded_p0(): + """Fixing p1 with p0 unbounded recovers p0 near the +-180 limit.""" + p0_true = 170.0 + data = _lorentzian_spectrum(p0=p0_true) + _, opt = autops(data, 'acme', bounds=[(None, None), (0, 0)], + return_phases=True) + assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" + assert opt[1] == 0.0 + + @pytest.mark.parametrize('bounds', [ [(-180, 180)], [(-180, 180), (-1800, 1800), (0, 0)], From 2743c012d8cf64175c9d9e360029ca8be798582e Mon Sep 17 00:00:00 2001 From: "Daniel C. Pizetta" Date: Wed, 7 Oct 2026 10:13:38 -0300 Subject: [PATCH 35/38] Replace autops bounds with p1_bounds and leave p0 unbounded The zero-order phase is periodic, so hard limits on it can trap the optimiser at a limit while p1 absorbs the remaining error (e.g. a 170 degree p0 error gave p0=-180, p1=22). p0 is now always optimised without limits and wrapped into [-180, 180) on return; only p1 can be bounded, via p1_bounds=(min, max). Co-Authored-By: Claude Opus 5.5 --- nmrglue/process/proc_autophase.py | 79 +++++------- tests/test_autophase.py | 192 +++++++++++++----------------- 2 files changed, 115 insertions(+), 156 deletions(-) diff --git a/nmrglue/process/proc_autophase.py b/nmrglue/process/proc_autophase.py index 70496126..84b40afd 100644 --- a/nmrglue/process/proc_autophase.py +++ b/nmrglue/process/proc_autophase.py @@ -13,7 +13,7 @@ from .proc_base import ps # fmin keyword arguments and their scipy.optimize.minimize(Nelder-Mead) -# option equivalents, used when autops is called with bounds. +# option equivalents, used when autops is called with p1_bounds. _FMIN_TO_MINIMIZE_OPTIONS = { 'xtol': 'xatol', 'ftol': 'fatol', @@ -25,7 +25,7 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, - bounds=None, **kwargs): + p1_bounds=None, **kwargs): """ Automatic linear phase correction @@ -46,42 +46,28 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, return_phases : bool If True, also return the optimised ``[p0, p1]`` list in addition to the phased data. - bounds : sequence of two (min, max) pairs or None, optional - Bounds on the zero- and first-order phases, respectively, in degrees. - For example ``bounds=[(-180, 180), (-1800, 1800)]`` restricts p0 to - +-180 and p1 to +-1800 degrees. When *bounds* is provided the - optimisation uses :func:`scipy.optimize.minimize` with the Nelder-Mead - method (same simplex algorithm as :func:`scipy.optimize.fmin`) so the - solver behaviour is equivalent but constrained. When *bounds* is - ``None`` (default) the original unconstrained - :func:`scipy.optimize.fmin` is used and the function is fully - backward-compatible. - - Use ``(None, None)`` to leave one of the phases unbounded, e.g. - ``bounds=[(None, None), (-90, 90)]`` bounds only p1 and - ``bounds=[(-180, 180), (None, None)]`` bounds only p0. A pair with - equal limits fixes that phase, e.g. ``(0, 0)`` for p1 to find only - the zero-order phase. - - Bounds are useful to find only the zero-order phase, or to prevent - the optimiser from wandering into physically unreasonable regions - when the spectrum has poor initial phase. - - .. note:: The bounds are hard limits, while the zero-order phase is - periodic (360 degrees). If the required p0 correction lies close - to a p0 limit, the optimiser can stop at that limit and compensate - with a wrong p1 instead of reaching the correct value. When only - p1 needs to be restricted, leave p0 unbounded with - ``(None, None)``; otherwise keep a margin between the p0 limits - and the expected correction. + p1_bounds : (min, max) or None, optional + Bounds on the first-order phase in degrees, e.g. ``(-1800, 1800)``. + Either limit may be ``None`` for no bound on that side, and equal + limits fix p1, e.g. ``(0, 0)`` to find only the zero-order phase. + When *p1_bounds* is provided the optimisation uses + :func:`scipy.optimize.minimize` with the Nelder-Mead method (same + simplex algorithm as :func:`scipy.optimize.fmin`) and the returned p0 + is wrapped into [-180, 180) degrees. When *p1_bounds* is ``None`` + (default) the original unconstrained :func:`scipy.optimize.fmin` is + used and the function is fully backward-compatible. + + The zero-order phase is never bounded: it is periodic, and hard + limits on it can trap the optimiser at a limit while p1 compensates + for the remaining error. kwargs : additional key-word arguments Passed directly to the underlying solver. - When *bounds* is ``None`` (fmin): see + When *p1_bounds* is ``None`` (fmin): see https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin.html Useful options: ``disp`` (bool), ``ftol`` (float). - When *bounds* is provided (minimize): see + When *p1_bounds* is provided (minimize): see https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html The method is fixed to 'Nelder-Mead'; pass e.g. ``options={'xatol': 1e-4}``. The fmin-style arguments ``xtol``, @@ -103,22 +89,14 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, >>> phased = ng.proc_autophase.autops(data, 'acme') - Bounded autophase — restrict search to +-180 in p0, +-1800 in p1: + Restrict the first-order phase to +-1800 degrees: >>> phased, phases = ng.proc_autophase.autops( - ... data, 'acme', p0=0.0, p1=0.0, - ... bounds=[(-180, 180), (-1800, 1800)], - ... return_phases=True) + ... data, 'acme', p1_bounds=(-1800, 1800), return_phases=True) - Zero-order only — fix p1 to zero and leave p0 unbounded: + Zero-order only — fix p1 to zero: - >>> phased = ng.proc_autophase.autops( - ... data, 'acme', bounds=[(None, None), (0, 0)]) - - Bound only the first-order phase: - - >>> phased = ng.proc_autophase.autops( - ... data, 'acme', bounds=[(None, None), (-90, 90)]) + >>> phased = ng.proc_autophase.autops(data, 'acme', p1_bounds=(0, 0)) """ arguments = [data] @@ -140,14 +118,14 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, opt = [p0, p1] - if bounds is not None: - if len(bounds) != 2 or any(len(b) != 2 for b in bounds): + if p1_bounds is not None: + if len(p1_bounds) != 2 or not all( + b is None or np.isscalar(b) for b in p1_bounds): raise ValueError( - 'bounds must contain two (min, max) pairs, one for p0 and ' - f'one for p1, got {bounds!r}' + f'p1_bounds must be a (min, max) pair, got {p1_bounds!r}' ) # Translate fmin-style keyword arguments into minimize options so the - # same call works with and without bounds. + # same call works with and without p1_bounds. options = {} for fmin_key, option_key in _FMIN_TO_MINIMIZE_OPTIONS.items(): if fmin_key in kwargs: @@ -156,11 +134,12 @@ def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, result = scipy.optimize.minimize( fn, x0=opt, args=tuple(arguments), method='Nelder-Mead', - bounds=bounds, + bounds=[(None, None), tuple(p1_bounds)], options=options, **kwargs, ) opt = result.x + opt[0] = (opt[0] + 180) % 360 - 180 else: opt = scipy.optimize.fmin(fn, x0=opt, args=tuple(arguments), **kwargs) diff --git a/tests/test_autophase.py b/tests/test_autophase.py index 3c834095..61eaf7df 100644 --- a/tests/test_autophase.py +++ b/tests/test_autophase.py @@ -1,10 +1,10 @@ """ tests/test_autophase.py — Tests for nmrglue.process.proc_autophase -Covers autops() with and without bounds, both built-in scoring functions, +Covers autops() with and without p1_bounds, both built-in scoring functions, callable scoring functions, return_phases flag, and error handling. -The bounded path uses scipy.optimize.minimize (Nelder-Mead); the unbounded +The p1-bounded path uses scipy.optimize.minimize (Nelder-Mead); the unbounded path uses scipy.optimize.fmin. Both are exercised here. """ import numpy as np @@ -27,7 +27,8 @@ # Helpers # ============================================================================= -def _lorentzian_spectrum(n=512, sw=10000.0, f0=1000.0, t2=0.01, p0=0.0, p1=0.0): +def _lorentzian_spectrum(n=512, sw=10000.0, f0=1000.0, t2=0.01, + p0=0.0, p1=0.0): """Return a phase-shifted Lorentzian spectrum (complex ndarray). Parameters @@ -103,170 +104,149 @@ def test_autops_acme_reduces_phase_error(): # ============================================================================= -# Bounded path — scipy.optimize.minimize(Nelder-Mead) +# p1-bounded path — scipy.optimize.minimize(Nelder-Mead) # ============================================================================= -def test_autops_bounded_returns_ndarray(): - """autops with bounds= returns an ndarray of the same shape.""" +def _wrapped_error(p0_true, p0_found): + """Residual zero-order phase error wrapped into [-180, 180).""" + return (p0_true + p0_found + 180) % 360 - 180 + + +def test_autops_p1_bounded_returns_ndarray(): + """autops with p1_bounds= returns an ndarray of the same shape.""" data = _lorentzian_spectrum(p0=30.0) - result = autops(data, 'acme', bounds=[(-180, 180), (-1800, 1800)]) + result = autops(data, 'acme', p1_bounds=(-1800, 1800)) assert isinstance(result, np.ndarray) assert result.shape == data.shape -def test_autops_bounded_return_phases(): - """return_phases=True with bounds returns (ndarray, array-like).""" +def test_autops_p1_bounded_return_phases(): + """return_phases=True with p1_bounds returns (ndarray, array-like).""" data = _lorentzian_spectrum(p0=45.0) - phased, opt = autops(data, 'acme', - bounds=[(-180, 180), (-1800, 1800)], + phased, opt = autops(data, 'acme', p1_bounds=(-1800, 1800), return_phases=True) assert isinstance(phased, np.ndarray) assert len(opt) == 2 def test_autops_p1_fixed_to_zero(): - """Bounds (0, 0) on p1 keep the first-order phase exactly zero.""" - data = _lorentzian_spectrum(p0=50.0, p1=0.0) - _, opt = autops(data, 'acme', - bounds=[(-180, 180), (0, 0)], - return_phases=True) - assert abs(opt[1]) < 1e-6, f"p1 should be 0.0, got {opt[1]}" + """p1_bounds=(0, 0) keeps the first-order phase exactly zero.""" + data = _lorentzian_spectrum(p0=50.0) + _, opt = autops(data, 'acme', p1_bounds=(0, 0), return_phases=True) + assert opt[1] == 0.0, f"p1 should be 0.0, got {opt[1]}" + assert abs(_wrapped_error(50.0, opt[0])) < PHASE_TOL, f"p0={opt[0]}" + + +def test_autops_p1_bounded_respects_range(): + """Optimised p1 stays within p1_bounds.""" + data = _lorentzian_spectrum(p1=300.0) + b1 = (-200, 200) + _, opt = autops(data, 'acme', p1_bounds=b1, return_phases=True) + assert b1[0] <= opt[1] <= b1[1], f"p1={opt[1]} outside {b1}" -def test_autops_p0_fixed_to_zero(): - """Bounds (0, 0) on p0 keep the zero-order phase exactly zero.""" - data = _lorentzian_spectrum(p0=0.0, p1=200.0) - _, opt = autops(data, 'acme', - bounds=[(0, 0), (-1800, 1800)], - return_phases=True) - assert abs(opt[0]) < 1e-6, f"p0 should be 0.0, got {opt[0]}" +@pytest.mark.parametrize('p1_bounds', [(None, 1800), (-1800, None), + (None, None)]) +def test_autops_p1_bounds_open_side(p1_bounds): + """None leaves that side of the p1 range unbounded.""" + data = _lorentzian_spectrum(p0=60.0, p1=30.0) + _, opt = autops(data, 'acme', p1_bounds=p1_bounds, return_phases=True) + assert abs(_wrapped_error(60.0, opt[0])) < PHASE_TOL, f"p0={opt[0]}" + assert abs(30.0 + opt[1]) < PHASE_TOL, f"p1={opt[1]}" -def test_autops_bounded_respects_p0_range(): - """Optimised p0 stays within the specified bounds.""" - data = _lorentzian_spectrum(p0=30.0) - b0 = (-90, 90) - _, opt = autops(data, 'acme', - bounds=[b0, (-1800, 1800)], +def test_autops_p1_bounded_acme_reduces_phase_error(): + """p1-bounded autops 'acme' undoes a known phase error.""" + p0_true, p1_true = 60.0, 30.0 + data = _lorentzian_spectrum(p0=p0_true, p1=p1_true) + _, opt = autops(data, 'acme', p1_bounds=(-1800, 1800), return_phases=True) - assert b0[0] <= opt[0] <= b0[1], f"p0={opt[0]} outside {b0}" + assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" + assert abs(p1_true + opt[1]) < PHASE_TOL, f"p1={opt[1]}" -def test_autops_bounded_respects_p1_range(): - """Optimised p1 stays within the specified bounds.""" - data = _lorentzian_spectrum(p0=0.0, p1=300.0) - b1 = (-500, 500) - _, opt = autops(data, 'acme', - bounds=[(-180, 180), b1], +@pytest.mark.parametrize('p0_true', range(-175, 180, 5)) +def test_autops_p1_bounded_recovers_any_p0(p0_true): + """p0 is recovered over the whole circle, including near +-180.""" + data = _lorentzian_spectrum(p0=p0_true) + _, opt = autops(data, 'acme', p1_bounds=(-1800, 1800), return_phases=True) - assert b1[0] <= opt[1] <= b1[1], f"p1={opt[1]} outside {b1}" + assert abs(_wrapped_error(p0_true, opt[0])) < PHASE_TOL, f"p0={opt[0]}" + assert abs(opt[1]) < PHASE_TOL, f"p1={opt[1]}" -def test_autops_bounded_acme_reduces_phase_error(): - """Bounded autops 'acme' reduces a known phase error towards zero.""" - p0_true = 60.0 - data = _lorentzian_spectrum(p0=p0_true) - _, opt = autops(data, 'acme', - bounds=[(-180, 180), (-1800, 1800)], +@pytest.mark.parametrize('p0_start', [190.0, -530.0]) +def test_autops_p1_bounded_p0_wrapped(p0_start): + """The returned p0 is wrapped into [-180, 180).""" + data = _lorentzian_spectrum(p0=170.0) + _, opt = autops(data, 'acme', p0=p0_start, p1_bounds=(-1800, 1800), return_phases=True) - assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" - assert abs(opt[1]) < PHASE_TOL, f"p1={opt[1]}" + assert -180 <= opt[0] < 180, f"p0={opt[0]}" + assert abs(_wrapped_error(170.0, opt[0])) < PHASE_TOL, f"p0={opt[0]}" -def test_autops_bounded_peak_minima(): - """Bounded autops works with 'peak_minima' scoring function.""" +def test_autops_p1_bounded_peak_minima(): + """p1-bounded autops works with 'peak_minima' scoring function.""" data = _lorentzian_spectrum(p0=40.0) - result = autops(data, 'peak_minima', - bounds=[(-180, 180), (-1800, 1800)]) + result = autops(data, 'peak_minima', p1_bounds=(-1800, 1800)) assert isinstance(result, np.ndarray) assert result.shape == data.shape -def test_autops_bounded_callable_fn(): - """Bounded autops accepts a callable scoring function.""" +def test_autops_p1_bounded_callable_fn(): + """p1-bounded autops accepts a callable scoring function.""" data = _lorentzian_spectrum(p0=25.0) - phased, opt = autops(data, _ps_acme_score, - bounds=[(-180, 180), (-1800, 1800)], + phased, opt = autops(data, _ps_acme_score, p1_bounds=(-1800, 1800), return_phases=True) assert isinstance(phased, np.ndarray) assert len(opt) == 2 -def test_autops_bounded_unknown_fn_raises(): - """An unrecognised string fn raises KeyError even when bounds are set.""" +def test_autops_p1_bounded_unknown_fn_raises(): + """An unrecognised string fn raises KeyError even with p1_bounds.""" data = _lorentzian_spectrum() with pytest.raises(KeyError, match='Unable to find algorithm'): - autops(data, 'not_a_real_algorithm', - bounds=[(-180, 180), (-1800, 1800)]) - - -@pytest.mark.parametrize('bounds', [ - [(-180, 180), (None, None)], - [(None, None), (-90, 90)], - [(-np.inf, np.inf), (-1800, 1800)], -]) -def test_autops_bounds_single_phase(bounds): - """(None, None) or infinite limits leave one phase unbounded.""" - p0_true, p1_true = 60.0, 30.0 - data = _lorentzian_spectrum(p0=p0_true, p1=p1_true) - _, opt = autops(data, 'acme', bounds=bounds, return_phases=True) - assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" - assert abs(p1_true + opt[1]) < PHASE_TOL, f"p1={opt[1]}" - - -def test_autops_p0_only_unbounded_p0(): - """Fixing p1 with p0 unbounded recovers p0 near the +-180 limit.""" - p0_true = 170.0 - data = _lorentzian_spectrum(p0=p0_true) - _, opt = autops(data, 'acme', bounds=[(None, None), (0, 0)], - return_phases=True) - assert abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" - assert opt[1] == 0.0 + autops(data, 'not_a_real_algorithm', p1_bounds=(-1800, 1800)) -@pytest.mark.parametrize('bounds', [ - [(-180, 180)], - [(-180, 180), (-1800, 1800), (0, 0)], - [(-180, 180), (0,)], +@pytest.mark.parametrize('p1_bounds', [ + (0,), + (-1800, 0, 1800), + [(-180, 180), (-1800, 1800)], ]) -def test_autops_invalid_bounds_raises(bounds): - """bounds must be exactly two (min, max) pairs.""" +def test_autops_invalid_p1_bounds_raises(p1_bounds): + """p1_bounds must be a single (min, max) pair.""" data = _lorentzian_spectrum(p0=30.0) - with pytest.raises(ValueError, match='two \\(min, max\\) pairs'): - autops(data, 'acme', bounds=bounds) + with pytest.raises(ValueError, match='p1_bounds must be a'): + autops(data, 'acme', p1_bounds=p1_bounds) -def test_autops_bounded_kwargs_passed_to_minimize(): +def test_autops_p1_bounded_kwargs_passed_to_minimize(): """Extra kwargs are forwarded to scipy.optimize.minimize.""" data = _lorentzian_spectrum(p0=30.0) # options= is a minimize-specific kwarg; should not raise - phased, opt = autops(data, 'acme', - bounds=[(-180, 180), (-1800, 1800)], + phased, opt = autops(data, 'acme', p1_bounds=(-1800, 1800), return_phases=True, options={'xatol': 1e-2, 'fatol': 1e-2, 'disp': False}) assert isinstance(phased, np.ndarray) -def test_autops_bounded_accepts_fmin_kwargs(): - """fmin-style kwargs (disp, xtol, ftol, ...) also work with bounds.""" +def test_autops_p1_bounded_accepts_fmin_kwargs(): + """fmin-style kwargs (disp, xtol, ftol, ...) also work with p1_bounds.""" data = _lorentzian_spectrum(p0=30.0) - _, opt = autops(data, 'acme', - bounds=[(-180, 180), (-1800, 1800)], + _, opt = autops(data, 'acme', p1_bounds=(-1800, 1800), return_phases=True, disp=False, xtol=1e-4, ftol=1e-4, maxiter=2000) assert abs(30.0 + opt[0]) < PHASE_TOL, f"p0={opt[0]}" -def test_autops_bounded_vs_unbounded_agreement(): - """Bounded and unbounded results are consistent on a well-conditioned spectrum.""" - data = _lorentzian_spectrum(p0=30.0, p1=0.0) +def test_autops_p1_bounded_vs_unbounded_agreement(): + """p1-bounded and unbounded results agree on a clean spectrum.""" + data = _lorentzian_spectrum(p0=30.0, p1=20.0) _, opt_ub = autops(data, 'acme', return_phases=True, disp=False) - _, opt_b = autops(data, 'acme', - bounds=[(-180, 180), (-1800, 1800)], - return_phases=True) - # Both should find similar p0 (within solver tolerance + ACME heuristic noise) - assert abs(opt_ub[0] - opt_b[0]) < 10.0, ( - f"Unbounded p0={opt_ub[0]:.2f} and bounded p0={opt_b[0]:.2f} disagree") + _, opt_b = autops(data, 'acme', p1_bounds=(-1800, 1800), + return_phases=True) + assert_allclose(opt_b, opt_ub, atol=1.0) # ============================================================================= From 7851bef1fcb52636fc0335c2b5a92f5b92ea1247 Mon Sep 17 00:00:00 2001 From: reneeotten Date: Wed, 7 Oct 2026 09:37:44 -0400 Subject: [PATCH 36/38] FIX: update FDSPECNUM properly when using pipe_proc.ext on 2D data set Closes: #276 --- nmrglue/process/pipe_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nmrglue/process/pipe_proc.py b/nmrglue/process/pipe_proc.py index 85583543..2f6802de 100644 --- a/nmrglue/process/pipe_proc.py +++ b/nmrglue/process/pipe_proc.py @@ -2039,7 +2039,7 @@ def ext(dic, data, x1="default", xn="default", y1="default", yn="default", # print("xmin:",x_min,"xmax:",x_max) data = data[int(y_min):int(y_max), int(x_min):int(x_max)] - if y_min != 1 and y_max != data.shape[0]: # only update when sliced + if y_max - y_min != old_y: # only update when sliced dic["FDSLICECOUNT"] = y_max - y_min dic["FDSPECNUM"] = y_max - y_min dic["FDSIZE"] = x_max - x_min From b72a9595151eab60b8f8deb2ce93ab54f76c1e12 Mon Sep 17 00:00:00 2001 From: reneeotten Date: Wed, 7 Oct 2026 23:08:27 -0400 Subject: [PATCH 37/38] FIX: correct type, define y_max instead of y_min --- nmrglue/process/pipe_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nmrglue/process/pipe_proc.py b/nmrglue/process/pipe_proc.py index 2f6802de..400e1daf 100644 --- a/nmrglue/process/pipe_proc.py +++ b/nmrglue/process/pipe_proc.py @@ -2033,7 +2033,7 @@ def ext(dic, data, x1="default", xn="default", y1="default", yn="default", if y_max > data.shape[0]: y_min = y_min - (y_max - data.shape[0]) - y_min = data.shape[0] + y_max = data.shape[0] # print("ymin:",y_min,"ymax:",y_max) # print("xmin:",x_min,"xmax:",x_max) From cd0c2c9d1cac4f1125a89dbe6ecdcd8ca194b59d Mon Sep 17 00:00:00 2001 From: Juuso Lehtivarjo Date: Thu, 8 Oct 2026 16:12:16 +0300 Subject: [PATCH 38/38] Fixes to JCAMP-DX FID handling (#231) * Fixes to JCAMP-DX udic & FID handling --------- Co-authored-by: Juuso Lehtivarjo --- nmrglue/fileio/jcampdx.py | 92 +++++++++++++++++++++++++++++---------- tests/test_jcampdx.py | 37 ++++++++++++---- 2 files changed, 99 insertions(+), 30 deletions(-) diff --git a/nmrglue/fileio/jcampdx.py b/nmrglue/fileio/jcampdx.py index 6582fa07..570c360a 100644 --- a/nmrglue/fileio/jcampdx.py +++ b/nmrglue/fileio/jcampdx.py @@ -136,7 +136,6 @@ def _parsejcampdx(filename): if currentkey is None: warn(f"JCAMP-DX data line without associated key: {line}") continue - currentvaluestrings.append(commentsplit[0]) # push possible non-closed blocks @@ -375,7 +374,8 @@ def _parse_pseudo(datalines): valuechar = _DUP_DIGITS[char] newmode = 3 except KeyError: - warn(f"Unknown pseudo-digit: {char} at line: {dataline}") + warn( + f"Unknown pseudo-digit: {char} at line: {dataline}") return None # finish previous number @@ -746,6 +746,23 @@ def _find_firstx_lastx(dic): return firstx, lastx, isppm +def get_complex_array(data): + """ + By default JCAMP FID data are read as two separate arrays + for real and imaginary parts. This function returns + them combined for single complex array. + """ + + if len(data) != 2: + warn("data is not list of arrays [real, imag]") + return None + + complexdata = np.empty((len(data[0]), ), dtype='complex128') + complexdata.real = data[0][:] + complexdata.imag = data[1][:] + return complexdata + + def guess_udic(dic, data): """ Guess parameters of universal dictionary from dic, data pair. @@ -790,32 +807,63 @@ def guess_udic(dic, data): pass # "size" - if isinstance(data, list): - data = data[0] # if list [R,I] + npoints = None if data is not None: - udic[0]["size"] = len(data) + if isinstance(data, list): + npoints = len(data[0]) # if list [R,I] + else: + npoints = len(data) + udic[0]["size"] = npoints else: warn('No data, cannot set udic size') - # "sw" - # get firstx, lastx and unit - firstx, lastx, isppm = _find_firstx_lastx(dic) - - # ppm data: convert to Hz - if isppm: - if obs_freq: - firstx = firstx * obs_freq - lastx = lastx * obs_freq - else: - firstx, lastx = (None, None) - warn('Data is in ppm but have no frequency, cannot set udic sweep') + # check if this is fid or processed: + is_processed = True # by default, expect processed data + try: + datatype = dic["DATATYPE"][0] + if datatype.strip().upper().replace(" ", "") == "NMRFID": + is_processed = False + except KeyError: + pass + # "sw" and "car" + # get firstx, lastx and unit + firstx, lastx, is_ppm = _find_firstx_lastx(dic) if firstx is not None and lastx is not None: - udic[0]["sw"] = abs(lastx - firstx) + if is_processed: + # ppm data: convert to Hz + if is_ppm: + if obs_freq: + firstx = firstx * obs_freq + lastx = lastx * obs_freq + else: + firstx, lastx = (None, None) + warn('Data is in ppm but base frequency is unknown, \ + cannot set udic spectral width') + if firstx is not None and lastx is not None: + udic[0]["sw"] = abs(lastx - firstx) + udic[0]["car"] = (lastx + firstx) / 2 + else: + # FID: + if npoints: + aqtime = lastx - firstx + sw = (npoints - 1) / aqtime + udic[0]["sw"] = sw + # note: in FIDs "car" is left to default as there is no required + # standard tag in JCAMP for it. Quite often manufacturers store + # it under their own tags though. else: - warn('Cannot set udic sweep') - - # keys not found in standard&required JCAMP-DX keys and thus left default: - # car, complex, encoding + warn('No data ranges found from JCAMP, cannot set udic sw & car') + + # "time" & "freq" + udic[0]["freq"] = is_processed + udic[0]["time"] = not is_processed + + # "complex" always false in JCAMP (R&I in separate arrays) + udic[0]["complex"] = False + # ...unless combined by nmrglue user with get_complex_array + if not isinstance(data, list): + if np.iscomplexobj(data): + udic[0]["complex"] = True return udic diff --git a/tests/test_jcampdx.py b/tests/test_jcampdx.py index ce9173a0..cb2c0270 100644 --- a/tests/test_jcampdx.py +++ b/tests/test_jcampdx.py @@ -109,10 +109,13 @@ def test_jcampdx1(): # check udic: udic = ng.jcampdx.guess_udic(dic, data) - assert np.abs(udic[0]["obs"]-100.4) < epsilon_e - assert np.abs(udic[0]["sw"]-24038.5) < epsilon_e + assert np.abs(udic[0]["obs"] - 100.4) < epsilon_e + assert np.abs(udic[0]["sw"] - 24038.5) < epsilon_e assert udic[0]["size"] == npoints_target assert udic[0]["label"] == "13C" + assert udic[0]["time"] is False + assert udic[0]["freq"] is True + assert udic[0]["complex"] is False def test_jcampdx2(): @@ -122,7 +125,8 @@ def test_jcampdx2(): # npoints, first, last, freq, sweep # note: first and last are raw values from datalines for convenience, # i.e. not scaled with YFACTORS - cases.append(("TESTFID.DX", 16384, 573, -11584, 100.4, 0.6815317)) + cases.append(("TESTFID.DX", 16384, 573, - + 11584, 100.4, (16384-1)/0.6815317)) cases.append(("bruker1.dx", 16384, -5, -51, 200.13, 4098.3606557377)) cases.append(("bruker2.dx", 16384, 42, 422, 300.13336767, 6024.096385479)) cases.append(("bruker3.dx", 16384, 22, -313, 300.13336729, 6024.096385479)) @@ -143,9 +147,11 @@ def test_jcampdx2(): print(case[0]) # read casepath = os.path.join(DATA_DIR, "jcampdx", case[0]) - dic, data = ng.jcampdx.read(casepath) - if isinstance(data, list): - data = data[0] # for data with both R&I, check only R + dic, rawdata = ng.jcampdx.read(casepath) + if isinstance(rawdata, list): + data = rawdata[0] # for data with both R&I, check only R + else: + data = rawdata # since first and last are raw values, do yfactor # back-scaling here @@ -164,8 +170,22 @@ def test_jcampdx2(): # check udic udic = ng.jcampdx.guess_udic(dic, data) - assert np.abs(udic[0]["obs"]-case[4]) < epsilon - assert np.abs(udic[0]["sw"]-case[5]) < epsilon + assert np.abs(udic[0]["obs"] - case[4]) < epsilon + assert np.abs(udic[0]["sw"] - case[5]) < epsilon + if dic["DATATYPE"][0] == "NMR FID": + assert udic[0]["time"] is True + assert udic[0]["freq"] is False + assert udic[0]["complex"] is False + assert len(rawdata) == 2 # [R,I] + # FID: merge raw arrays and assert changes + data = ng.jcampdx.get_complex_array(rawdata) + assert len(data) == case[1] + udic = ng.jcampdx.guess_udic(dic, data) + assert udic[0]["complex"] is True + else: + assert udic[0]["time"] is False + assert udic[0]["freq"] is True + assert udic[0]["complex"] is False def test_jcampdx_dicstructure(): @@ -195,6 +215,7 @@ def test_jcampdx_dicstructure(): assert dic["_datatype_NA"][0]["_comments"][0] == "comment line" assert dic["_datatype_NA"][0]["_comments"][1] == "another comment" + def test_jcampdx_dicstructure2(): '''JCAMP-DX read: ensure correct dic structure (nested LINKs, multiple spectra)'''