diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 19986753..efd34178 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,13 +1,18 @@ name: nmrglue tests -on: [push] +on: + push: + branches: master + pull_request: + branches: master jobs: build: 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: - uses: actions/checkout@v3 diff --git a/CHANGELOG.md b/CHANGELOG.md index 8969d2a6..f1d1a7e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,22 @@ -0.12 (Unreleased) +0.13 (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/examples/fitting_data/t1_measurements/extract_trajs.py b/examples/fitting_data/t1_measurements/extract_trajs.py index f96628b1..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.recfromtxt("boxes.in", encoding="UTF-8") -spectra_list = np.recfromtxt("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 eaa967a7..2bccc8c4 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", dtype=None) # 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..96705870 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", 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 4645f623..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.recfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in", dtype=None) # 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..71af4bc8 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", dtype=None) # 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..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.recfromtxt("boxes.in", names=True) -spectra_list = np.recfromtxt("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 a364b938..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.recfromtxt("boxes.in", names=True) -spectra_list = np.recfromtxt("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): @@ -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..3aeef6c4 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', dtype=None) 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..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.recfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in", dtype=None) # 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..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.recfromtxt("limits.in") +peak_list = np.genfromtxt("limits.in", dtype=None) # 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..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.recfromtxt('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 diff --git a/nmrglue/__init__.py b/nmrglue/__init__.py index 88c7a986..8f141ce7 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.13-dev' diff --git a/nmrglue/fileio/bruker.py b/nmrglue/fileio/bruker.py index 9de2f52a..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: @@ -426,22 +430,29 @@ 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 + 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, estimated_dims=estimated_dims) + + return dic, data @@ -768,9 +779,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: @@ -1052,19 +1063,25 @@ def guess_shape(dic): except KeyError: dtypa = 0 # default value, int32 data + if dtypa == 0: + bytesize = 4 + elif dtypa == 2: + bytesize = 8 + else: + 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: - 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 @@ -1521,7 +1538,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. @@ -1570,7 +1587,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 @@ -2155,7 +2183,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 +2201,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): """ @@ -2706,4 +2751,49 @@ 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 + +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 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. diff --git a/nmrglue/fileio/jcampdx.py b/nmrglue/fileio/jcampdx.py index 1d41e2b5..570c360a 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 = [] @@ -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/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): 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/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..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 @@ -24,8 +25,8 @@ TNTMAGIC_RE = re.compile(br"^TNT1\.\d\d\d$") -TNTMAGIC = np.dtype('a8') -TNTTLV = np.dtype([('tag', 'a4'), ('bool', ' 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_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""" 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 diff --git a/nmrglue/process/pipe_proc.py b/nmrglue/process/pipe_proc.py index 85583543..400e1daf 100644 --- a/nmrglue/process/pipe_proc.py +++ b/nmrglue/process/pipe_proc.py @@ -2033,13 +2033,13 @@ 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) 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 diff --git a/nmrglue/process/proc_autophase.py b/nmrglue/process/proc_autophase.py index 7d5479cc..84b40afd 100644 --- a/nmrglue/process/proc_autophase.py +++ b/nmrglue/process/proc_autophase.py @@ -12,8 +12,20 @@ from .proc_base import ps - -def autops(data, fn, p0=0.0, p1=0.0, return_phases=False, peak_width=100, **kwargs): +# fmin keyword arguments and their scipy.optimize.minimize(Nelder-Mead) +# option equivalents, used when autops is called with p1_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, + p1_bounds=None, **kwargs): """ Automatic linear phase correction @@ -29,29 +41,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. + 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 *p1_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 *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``, + ``ftol``, ``maxiter``, ``maxfun``, ``disp`` and ``initial_simplex`` + are also accepted and translated to the matching ``options`` entries + (explicit ``options`` take precedence). Returns ------- ndata : ndarray Phased NMR data. + opt : ndarray, 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') + + Restrict the first-order phase to +-1800 degrees: + + >>> phased, phases = ng.proc_autophase.autops( + ... data, 'acme', p1_bounds=(-1800, 1800), return_phases=True) + + Zero-order only — fix p1 to zero: + + >>> phased = ng.proc_autophase.autops(data, 'acme', p1_bounds=(0, 0)) """ - arguments = [data,] + arguments = [data] if not callable(fn): if fn == 'peak_minima': @@ -63,11 +111,37 @@ 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 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( + 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 p1_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=[(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) phasedspc = ps(data, p0=opt[0], p1=opt[1]) diff --git a/nmrglue/process/proc_base.py b/nmrglue/process/proc_base.py index d2bcfd6b..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 @@ -324,7 +336,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 +357,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 +384,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 +414,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 @@ -2644,7 +2659,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) diff --git a/nmrglue/process/proc_lp.py b/nmrglue/process/proc_lp.py index 4a8ef948..e309db7d 100644 --- a/nmrglue/process/proc_lp.py +++ b/nmrglue/process/proc_lp.py @@ -528,8 +528,9 @@ 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]) + # first K columns of V are first K rows of Vh + Vlh = np.asmatrix(Vh[:K, :]) sl = s[:K] if min_var: # adjust singular values using minimum variance method @@ -537,7 +538,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 @@ -656,7 +657,7 @@ def lp_model(trace, slice=slice(None), order=8, mode="f", mirror=None, # build the B matrix (a Vandermonde matrix) and solve for the coefficients poles = np.array(poles) - B = np.row_stack([poles ** (i) for i in range(len(x))]) + B = np.vstack([poles ** (i) for i in range(len(x))]) z, resid, rank, s = np.linalg.lstsq(B, np.array(x)) # Now the z_n = amp_n*exp(phase_n*i), use this to determine the amplitudes @@ -803,10 +804,11 @@ 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) + # the next 3 lines and the pinv_diagsvd function were adapted from the # scipy.linalg.pinv2 function - jjh eps = np.finfo('float').eps @@ -833,16 +835,16 @@ 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 + # 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): @@ -852,9 +854,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 +1020,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 diff --git a/setup.py b/setup.py index 5666a22b..e657b8f9 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.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', @@ -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', diff --git a/tests/test_autophase.py b/tests/test_autophase.py new file mode 100644 index 00000000..61eaf7df --- /dev/null +++ b/tests/test_autophase.py @@ -0,0 +1,281 @@ +""" +tests/test_autophase.py — Tests for nmrglue.process.proc_autophase + +Covers autops() with and without p1_bounds, both built-in scoring functions, +callable scoring functions, return_phases flag, and error handling. + +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 +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 + +# 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.01, + 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) + # 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) + + +# ============================================================================= +# 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) + # 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]}" + + +# ============================================================================= +# p1-bounded path — scipy.optimize.minimize(Nelder-Mead) +# ============================================================================= + +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', p1_bounds=(-1800, 1800)) + assert isinstance(result, np.ndarray) + assert result.shape == data.shape + + +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', p1_bounds=(-1800, 1800), + return_phases=True) + assert isinstance(phased, np.ndarray) + assert len(opt) == 2 + + +def test_autops_p1_fixed_to_zero(): + """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}" + + +@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_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 abs(p0_true + opt[0]) < PHASE_TOL, f"p0={opt[0]}" + assert abs(p1_true + opt[1]) < PHASE_TOL, f"p1={opt[1]}" + + +@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 abs(_wrapped_error(p0_true, opt[0])) < PHASE_TOL, f"p0={opt[0]}" + assert abs(opt[1]) < PHASE_TOL, f"p1={opt[1]}" + + +@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 -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_p1_bounded_peak_minima(): + """p1-bounded autops works with 'peak_minima' scoring function.""" + data = _lorentzian_spectrum(p0=40.0) + result = autops(data, 'peak_minima', p1_bounds=(-1800, 1800)) + assert isinstance(result, np.ndarray) + assert result.shape == data.shape + + +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, p1_bounds=(-1800, 1800), + return_phases=True) + assert isinstance(phased, np.ndarray) + assert len(opt) == 2 + + +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', p1_bounds=(-1800, 1800)) + + +@pytest.mark.parametrize('p1_bounds', [ + (0,), + (-1800, 0, 1800), + [(-180, 180), (-1800, 1800)], +]) +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='p1_bounds must be a'): + autops(data, 'acme', p1_bounds=p1_bounds) + + +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', p1_bounds=(-1800, 1800), + return_phases=True, + options={'xatol': 1e-2, 'fatol': 1e-2, 'disp': False}) + assert isinstance(phased, np.ndarray) + + +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', 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_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', p1_bounds=(-1800, 1800), + return_phases=True) + assert_allclose(opt_b, opt_ub, atol=1.0) + + +# ============================================================================= +# 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']) 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)''' 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) diff --git a/tests/test_tecmag.py b/tests/test_tecmag.py index 864f74b6..8acf1b87 100644 --- a/tests/test_tecmag.py +++ b/tests/test_tecmag.py @@ -1,21 +1,349 @@ +""" +tests/test_tecmag.py — Tests for nmrglue.fileio.tecmag + +Covers two layers: + +1. Integration tests — read a real .tnt file from DATA_DIR and verify + the parsed data against a known reference. + +2. Unit tests — synthesise minimal valid .tnt byte sequences to exercise + read(), guess_tables(), and the internal helpers _read_pascal_pair() + and _convert_si() without requiring external test data. +""" import os +import struct import numpy as np -from numpy.testing import assert_array_almost_equal +from numpy.testing import assert_allclose, assert_array_almost_equal +import pytest import nmrglue as ng +from nmrglue.fileio.tecmag import ( + TNTTMAG, TNTTMG2, guess_tables, _read_pascal_pair, _convert_si, +) from setup import DATA_DIR -def test_tecmag_load_time_domain(): - ref1, data = ng.tecmag.read(os.path.join( - DATA_DIR, 'tecmag', 'LiCl_ref1.tnt')) +# ============================================================================= +# Helpers — synthesise minimal valid .tnt byte sequences +# ============================================================================= + +def _pascal_string(s): + """Build a length-prefixed Pascal string (4-byte LE length + bytes).""" + 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]) + + +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(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(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(tmp_path, [ + ('rfamp', '1.0 5.0 10.0'), + ('tp', '50u'), + ('deadtime', '20u'), + ]) + tables = guess_tables(path) + assert 'rfamp' in tables + assert 'tp' not in tables + assert 'deadtime' not in tables + + +def test_guess_tables_name_filter(tmp_path): + """The `names` argument restricts results to the requested subset.""" + path = _write_tnt(tmp_path, [ + ('rfamp', '1.0 2.0 3.0'), + ('d1', '1m 2m 3m'), + ('d2', '1s 2s 3s'), + ]) + tables = guess_tables(path, names=['rfamp', 'd1']) + assert set(tables.keys()) == {'rfamp', 'd1'} + + +def test_guess_tables_convert_si_disabled(tmp_path): + """When convert_si=False, values are returned as raw string token lists.""" + 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(tmp_path): + """A file with no PSEQ section returns an empty dict without raising.""" + 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(tmp_path): + """A name with a mismatched length prefix is silently skipped.""" + bogus = struct.pack('