Merge pull request #189 from kaustubhmote/read_nus_bruker

Expanding NUS data

doc/source/reference/bruker.rst

@@ -35,8 +35,8 @@ These are functions which are targeted for users of nmrglue.
     write_pprog
     guess_udic
     create_dic
-
-
+    read_nuslist
+    
 
 Developer Information
 ---------------------

doc/source/reference/proc_base.rst

@@ -186,3 +186,12 @@ are included here for completeness.
     std_flt
     sum_flt
     largest_power_of_2
+
+
+Other processing functions
+--------------------------
+
+.. autosummary::
+    :toctree: generated/
+
+    expand_nus

nmrglue/fileio/bruker.py

@@ -428,8 +428,20 @@ def read(dir=".", bin_file=None, acqus_files=None, pprog_file=None, shape=None,
 
     # read the binary file
     f = os.path.join(dir, bin_file)
-    null, data = read_binary(f, shape=shape, cplex=cplex, big=big,
+    _, 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")
+    except KeyError:
+        # old datasets do not have the FnTYPE parameter in acqus files. 
+        # also fails silently when acqus file is absent.
+        pass
+
     return dic, data
 
 
@@ -537,8 +549,20 @@ def read_lowmem(dir=".", bin_file=None, acqus_files=None, pprog_file=None,
 
     # read the binary file
     f = os.path.join(dir, bin_file)
-    null, data = read_binary_lowmem(f, shape=shape, cplex=cplex, big=big,
+    _, data = read_binary_lowmem(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")
+    except KeyError:
+        # old datasets do not have the FnTYPE parameter in acqus files. 
+        # also fails silently when acqus file is absent.
+        pass
+
     return dic, data
 
 
@@ -2581,3 +2605,49 @@ def _merge_dict(a, b):
     c = a.copy()
     c.update(b)
     return c
+
+
+def read_nuslist(dirc=".", fname="nuslist"):
+    """
+    Reads nuslist in bruker format
+
+    Parameters
+    ----------
+    dirc : str, optional
+        directory for the data, by default "."
+    fname : str, optional
+        name of the file that has the nuslist, by default 'nuslist'
+
+    Returns
+    -------
+    converted_nuslist: list of n-tuples
+        nuslist
+
+    Raises
+    ------
+    OSError
+        if directory is absent
+
+    FileNotFoundError
+        if file is absent
+
+    """
+    if not os.path.isdir(dirc):
+        raise OSError(f"directory {dirc} does not exist")
+
+    if fname is None:
+        fname = "nuslist"
+
+    try:
+        with open(os.path.join(dirc, fname)) as f:
+            nuslist = f.read().splitlines()
+    except FileNotFoundError:
+        raise FileNotFoundError(f"nuslist file ({fname}) not found in directory {dirc}")
+
+    converted_nuslist = []
+    for line in nuslist:
+        numbers = tuple(int(i) for i in line.split())
+        converted_nuslist.append(numbers)
+
+    return converted_nuslist
+

nmrglue/process/proc_base.py

@@ -7,6 +7,7 @@
 
 # TODO determine which of these work on N-dimension and which assume 2D
 
+from itertools import product
 import numpy as np
 import scipy.signal
 import scipy.linalg
@@ -2641,3 +2642,101 @@ 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)
     return zd(data, window, x0=x0, slope=slope)
+
+
+def expand_nus(data, shape, nuslist, aqorder=None, quadrature_order=None, allow_high_dims=False):
+    """
+    Converts a non-uniformaly sampled dataset to a fully sampled
+    dataset. FIDs are sorted according the nuslist, and missing
+    fids replaced by zeros. All dimensions are supported, but
+    only 2D, 3D, and 4D are allowed by default to avoid creation
+    of large datasets by mistake.
+
+    Parameters
+    ----------
+    data : np.ndarray
+        non-uniformly sampled dataset
+    shape : tuple
+        shape of the fully sampled data
+    nuslist : list of n-tuples
+        nuslist
+    aqorder : list | None
+        the order in which indirect dimensions are acquired.
+        defaults to [1, 0] for 3D and [2, 1, 0] for 4D datasets. 
+        All other possibilites compatible with the dimension are 
+        allowed, for eg, [0, 1] for 3D, and [0, 1, 2], [0, 2, 1], 
+        [2, 0, 1], [1, 0, 2], [1, 2, 0] for 4D data.  
+    quadrature_order : list | None
+        ordering of quadrature points. by defualt, this uses
+        the order from itertools.product, which seems to be 
+        fine most of the common acquistions. For example, for a 2D dataset,
+        this will be [(0,), (1,)], for 3D, this will be [(0, 0), (0, 1), (1, 0), (1, 1)]
+        and for 4D [(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1),
+        (1, 1, 0), (1, 1, 1)], where 0=real and 1=imag point
+    allow_high_dims : bool
+        flag that allows datasets higher than 4D to be genereted.
+        By default, this is not allowed to prevent generation of large datasets
+        by mistake
+
+    Returns
+    -------
+    full_data : np.ndarray
+        data from the nus dataset sorted according to
+        the nuslist, and with missing fids replaced by
+        zeros
+
+    Raises
+    ------
+    ValueError
+        if shape of the final data is not compatible
+        with the expeted shape from nuslist
+    ValueError
+        if dataset dimension is less than 2
+    ValueError
+        if dataset dimension is greated than 4 and
+        the flag to allow this is False
+
+    """
+
+    idims = len(nuslist[0])
+    ndim = idims + 1
+
+    if len(shape) != ndim:
+        raise ValueError(f"expected {ndim}D data but got {ndim}D nuslist")
+
+    if ndim < 2:
+        raise ValueError(f"Needs to be be atleast a 2D dataset (not {ndim}D)")
+
+    # protection against unintended generation of large datasets
+    if (ndim > 4) and (allow_high_dims == False):
+        raise ValueError(
+            f"converting {ndim}D data to a fully sampled data not allowed by default. "
+            f"Use the flag 'allow_high_dims=True' if you really want to do this."
+        )
+
+    # make sure that nuslist is correctly ordered and typed
+    if aqorder is None:
+        aqorder = list(range(idims))[::-1]
+
+    reordered_nuslist = []
+    for item in nuslist:
+        reordered_item = []
+        for pos in aqorder:
+            reordered_item.append(item[pos])
+        reordered_nuslist.append(tuple(reordered_item))
+
+    # generate a full dataset
+    full_data = np.zeros(shape, dtype=data.dtype)
+
+    # quadrature points in all dimensions
+    if quadrature_order is None:
+        quadrature_order = list(product((0, 1), repeat=idims))
+
+    # add fids at appropriate positions
+    for i, nus_pos in enumerate(reordered_nuslist):
+        for j, quad_pos in enumerate(quadrature_order):
+            fid_pos = i * (2 ** idims) + j
+            final_pos = tuple(p * 2 + q for p, q in zip(nus_pos, quad_pos))
+            full_data[final_pos] = data[fid_pos]
+
+    return full_data

tests/test_bruker.py

@@ -210,3 +210,16 @@ def test_write_pdata_2d():
                                      'pdata', '1'), read_acqus=False)
     write_readback_pdata(dic=dic, data=data)
     write_readback_pdata(dic=dic, data=data, pdata_folder=90)
+
+
+def test_read_nuslist():
+    """ reading nuslist """
+    with open("tmp_nuslist", "w") as f:
+        f.write("""0 0\n10 20\n50 21\n9 8\n7 8\n20 20""")
+
+    nuslist = ng.bruker.read_nuslist(fname="tmp_nuslist")
+    assert nuslist == [(0, 0), (10, 20), (50, 21), (9, 8), (7, 8), (20, 20)]
+
+    os.remove("tmp_nuslist")
+
+    pass

tests/test_proc_base.py

@@ -0,0 +1,147 @@
+import numpy as np
+import nmrglue as ng
+
+
+def test_reorder_nus_data_2d():
+    """reorder 2d nus data"""
+    nus_data = np.linspace(0, 1, 128 * 4 * 2).reshape(-1, 128)
+    nuslist = [(0,), (2,), (5,), (7,)]
+    full_shape = (16, 128)
+    full_data = ng.proc_base.expand_nus(
+        data=nus_data, shape=full_shape, nuslist=nuslist
+    )
+
+    assert full_data.shape == full_shape
+    assert np.allclose(full_data[0], nus_data[0], atol=1e-7)
+    assert np.allclose(full_data[1], nus_data[1], atol=1e-7)
+    assert np.allclose(full_data[4], nus_data[2], atol=1e-7)
+    assert np.allclose(full_data[5], nus_data[3], atol=1e-7)
+    assert np.allclose(full_data[10], nus_data[4], atol=1e-7)
+    assert np.allclose(full_data[11], nus_data[5], atol=1e-7)
+    assert np.allclose(full_data[14], nus_data[6], atol=1e-7)
+    assert np.allclose(full_data[15], nus_data[7], atol=1e-7)
+
+
+def test_reorder_nus_data_3d():
+    """reorder 3d nus data"""
+    nus_data = np.linspace(0, 1, 128 * 4 * 2 * 2).reshape(-1, 128)
+    nuslist = [(0, 0), (2, 3), (5, 3), (6, 7)]
+    full_shape = (16, 16, 128)
+    full_data = ng.proc_base.expand_nus(
+        data=nus_data, shape=full_shape, nuslist=nuslist
+    )
+
+    assert full_data.shape == full_shape
+    assert np.allclose(full_data[0, 0], nus_data[0], atol=1e-7)
+    assert np.allclose(full_data[0, 1], nus_data[1], atol=1e-7)
+    assert np.allclose(full_data[1, 0], nus_data[2], atol=1e-7)
+    assert np.allclose(full_data[1, 1], nus_data[3], atol=1e-7)
+
+    assert np.allclose(full_data[6, 4], nus_data[4], atol=1e-7)
+    assert np.allclose(full_data[6, 5], nus_data[5], atol=1e-7)
+    assert np.allclose(full_data[7, 4], nus_data[6], atol=1e-7)
+    assert np.allclose(full_data[7, 5], nus_data[7], atol=1e-7)
+
+    assert np.allclose(full_data[6, 10], nus_data[8], atol=1e-7)
+    assert np.allclose(full_data[6, 11], nus_data[9], atol=1e-7)
+    assert np.allclose(full_data[7, 10], nus_data[10], atol=1e-7)
+    assert np.allclose(full_data[7, 11], nus_data[11], atol=1e-7)
+
+    assert np.allclose(full_data[14, 12], nus_data[12], atol=1e-7)
+    assert np.allclose(full_data[14, 13], nus_data[13], atol=1e-7)
+    assert np.allclose(full_data[15, 12], nus_data[14], atol=1e-7)
+    assert np.allclose(full_data[15, 13], nus_data[15], atol=1e-7)
+
+
+def test_reorder_nus_data_4d():
+    """reorder 4d nus data"""
+    nus_data = np.linspace(0, 1, 32 * 2 * 2 * 2 * 2).reshape(-1, 32)
+    nuslist = [(0, 0, 0), (2, 3, 1)]
+    full_shape = (8, 8, 8, 32)
+    full_data = ng.proc_base.expand_nus(
+        data=nus_data, shape=full_shape, nuslist=nuslist
+    )
+
+    assert full_data.shape == full_shape
+    assert np.allclose(full_data[0, 0, 0], nus_data[0], atol=1e-7)
+    assert np.allclose(full_data[0, 0, 1], nus_data[1], atol=1e-7)
+    assert np.allclose(full_data[0, 1, 0], nus_data[2], atol=1e-7)
+    assert np.allclose(full_data[0, 1, 1], nus_data[3], atol=1e-7)
+    assert np.allclose(full_data[1, 0, 0], nus_data[4], atol=1e-7)
+    assert np.allclose(full_data[1, 0, 1], nus_data[5], atol=1e-7)
+    assert np.allclose(full_data[1, 1, 0], nus_data[6], atol=1e-7)
+    assert np.allclose(full_data[1, 1, 1], nus_data[7], atol=1e-7)
+
+    assert np.allclose(full_data[2, 6, 4], nus_data[8], atol=1e-7)
+    assert np.allclose(full_data[2, 6, 5], nus_data[9], atol=1e-7)
+    assert np.allclose(full_data[2, 7, 4], nus_data[10], atol=1e-7)
+    assert np.allclose(full_data[2, 7, 5], nus_data[11], atol=1e-7)
+    assert np.allclose(full_data[3, 6, 4], nus_data[12], atol=1e-7)
+    assert np.allclose(full_data[3, 6, 5], nus_data[13], atol=1e-7)
+    assert np.allclose(full_data[3, 7, 4], nus_data[14], atol=1e-7)
+    assert np.allclose(full_data[3, 7, 5], nus_data[15], atol=1e-7)
+
+
+def test_reorder_nus_rev_aqorder():
+    """reorder 3d nus data"""
+    nus_data = np.linspace(0, 1, 128 * 4 * 2 * 2).reshape(-1, 128)
+    nuslist = [(0, 0), (2, 3), (5, 3), (6, 7)]
+    full_shape = (16, 16, 128)
+    full_data = ng.proc_base.expand_nus(
+        data=nus_data, shape=full_shape, nuslist=nuslist, aqorder=[0, 1]
+    )
+
+    assert full_data.shape == full_shape
+    assert np.allclose(full_data[0, 0], nus_data[0], atol=1e-7)
+    assert np.allclose(full_data[0, 1], nus_data[1], atol=1e-7)
+    assert np.allclose(full_data[1, 0], nus_data[2], atol=1e-7)
+    assert np.allclose(full_data[1, 1], nus_data[3], atol=1e-7)
+
+    assert np.allclose(full_data[4, 6], nus_data[4], atol=1e-7)
+    assert np.allclose(full_data[4, 7], nus_data[5], atol=1e-7)
+    assert np.allclose(full_data[5, 6], nus_data[6], atol=1e-7)
+    assert np.allclose(full_data[5, 7], nus_data[7], atol=1e-7)
+
+    assert np.allclose(full_data[10, 6], nus_data[8], atol=1e-7)
+    assert np.allclose(full_data[10, 7], nus_data[9], atol=1e-7)
+    assert np.allclose(full_data[11, 6], nus_data[10], atol=1e-7)
+    assert np.allclose(full_data[11, 7], nus_data[11], atol=1e-7)
+
+    assert np.allclose(full_data[12, 14], nus_data[12], atol=1e-7)
+    assert np.allclose(full_data[12, 15], nus_data[13], atol=1e-7)
+    assert np.allclose(full_data[13, 14], nus_data[14], atol=1e-7)
+    assert np.allclose(full_data[13, 15], nus_data[15], atol=1e-7)
+
+
+def test_reorder_nus_3d_quadorder():
+    """reorder 3d nus data"""
+    nus_data = np.linspace(0, 1, 128 * 4 * 2 * 2).reshape(-1, 128)
+    nuslist = [(0, 0), (2, 3), (5, 3), (6, 7)]
+    full_shape = (16, 16, 128)
+    full_data = ng.proc_base.expand_nus(
+        data=nus_data,
+        shape=full_shape,
+        nuslist=nuslist,
+        quadrature_order=[(0, 0), (1, 0), (0, 1), (1, 1)],
+    )
+
+    assert full_data.shape == full_shape
+    assert np.allclose(full_data[0, 0], nus_data[0], atol=1e-7)
+    assert np.allclose(full_data[1, 0], nus_data[1], atol=1e-7)
+    assert np.allclose(full_data[0, 1], nus_data[2], atol=1e-7)
+    assert np.allclose(full_data[1, 1], nus_data[3], atol=1e-7)
+
+    assert np.allclose(full_data[6, 4], nus_data[4], atol=1e-7)
+    assert np.allclose(full_data[7, 4], nus_data[5], atol=1e-7)
+    assert np.allclose(full_data[6, 5], nus_data[6], atol=1e-7)
+    assert np.allclose(full_data[7, 5], nus_data[7], atol=1e-7)
+
+    assert np.allclose(full_data[6, 10], nus_data[8], atol=1e-7)
+    assert np.allclose(full_data[7, 10], nus_data[9], atol=1e-7)
+    assert np.allclose(full_data[6, 11], nus_data[10], atol=1e-7)
+    assert np.allclose(full_data[7, 11], nus_data[11], atol=1e-7)
+
+    assert np.allclose(full_data[14, 12], nus_data[12], atol=1e-7)
+    assert np.allclose(full_data[15, 12], nus_data[13], atol=1e-7)
+    assert np.allclose(full_data[14, 13], nus_data[14], atol=1e-7)
+    assert np.allclose(full_data[15, 13], nus_data[15], atol=1e-7)