Repository: jjhelmus/nmrglue

Issue:
Issue with ng.peakpick.pick for "thres" and "thres-fast" algorithms with Numpy ^v1.22
**Problem:**
When using `ng.peakpick.pick` with the algorithms "thres" and "thres-fast" and the Numpy version ^1.22, the following error occurs:

```
Error:
peakpick.py:369 # find the rectangular region around the segment
peakpick.py:370 region = data[seg_slice]
peakpick.py:371 edge = [s.start for s in seg_slice]
peakpick.py:372 rlocation = [l - s.start for l, s in zip(location, seg_slice)]

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`), and integer or boolean arrays are valid indices.
```

It appears to be related to the function call `find_all_thres_fast(data, pthres, msep, True)` since it returns a list for `seg_slice`.

**Potential Fix:**
Modifying line 370 in peakpick.py to:

```python
peakpick.py:370 region = data[seg_slice[0]]
```
or adjusting the return of the find_all_thres_fast function accordingly.

Somehow this is not needed for numpy versions below 1.22 but for above versions it becomes a problem.

Discussion:

Patch:
diff --git a/examples/coadd/coadd_2d_pipe/coadd_2d.py b/examples/coadd/coadd_2d_pipe/coadd_2d.py
index aa6ce54..27a35a8 100755
--- a/examples/coadd/coadd_2d_pipe/coadd_2d.py
+++ b/examples/coadd/coadd_2d_pipe/coadd_2d.py
@@ -15,10 +15,10 @@ coadd_dic = dict(dic)
 
 # loop over the files, adding them to the coadded array
 for f in flist:
-    print("Coadding file:", f)
+    print "Coadding file:", f
     dic, data = ng.pipe.read(f)
     coadd_data += data / len(flist)
 
 # write out the file
-print("Writing out coadded data")
+print "Writing out coadded data"
 ng.pipe.write("coadded.fid", coadd_dic, coadd_data, overwrite=True)
diff --git a/nmrglue/analysis/peakpick.py b/nmrglue/analysis/peakpick.py
index 8eeda0d..17e46f5 100644
--- a/nmrglue/analysis/peakpick.py
+++ b/nmrglue/analysis/peakpick.py
@@ -535,7 +535,7 @@ def find_pseg_slice(data, location, thres):
             stop = stop + 1
             al[dim] = stop
         seg_slice.append(slice(start + 1, stop))
-    return seg_slice
+    return tuple(seg_slice)
 
 
 def find_nseg_slice(data, location, thres):
@@ -558,4 +558,4 @@ def find_nseg_slice(data, location, thres):
             stop = stop + 1
             al[dim] = stop
         seg_slice.append(slice(start + 1, stop))
-    return seg_slice
+    return tuple(seg_slice)
diff --git a/nmrglue/process/proc_base.py b/nmrglue/process/proc_base.py
index d2bcfd6..43b79e6 100644
--- a/nmrglue/process/proc_base.py
+++ b/nmrglue/process/proc_base.py
@@ -2667,18 +2667,18 @@ def expand_nus(data, shape, nuslist, aqorder=None, quadrature_order=None, allow_
     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 possibilities compatible with the dimension are 
+        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 default, this uses
+        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 generated.
+        flag that allows datasets higher than 4D to be genereted.
         By default, this is not allowed to prevent generation of large datasets
         by mistake
 
@@ -2693,11 +2693,11 @@ def expand_nus(data, shape, nuslist, aqorder=None, quadrature_order=None, allow_
     ------
     ValueError
         if shape of the final data is not compatible
-        with the expected shape from nuslist
+        with the expeted shape from nuslist
     ValueError
         if dataset dimension is less than 2
     ValueError
-        if dataset dimension is greater than 4 and
+        if dataset dimension is greated than 4 and
         the flag to allow this is False
 
     """
@@ -2709,7 +2709,7 @@ def expand_nus(data, shape, nuslist, aqorder=None, quadrature_order=None, allow_
         raise ValueError(f"expected {ndim}D data but got {ndim}D nuslist")
 
     if ndim < 2:
-        raise ValueError(f"Needs to be be at least a 2D dataset (not {ndim}D)")
+        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):
diff --git a/nmrglue/process/proc_bl.py b/nmrglue/process/proc_bl.py
index 1f2dc90..ab6aff9 100644
--- a/nmrglue/process/proc_bl.py
+++ b/nmrglue/process/proc_bl.py
@@ -237,7 +237,7 @@ def baseline_corrector(data, wd=20):
     sigma = _find_noise_sd(sd_set, 0.999)
     sn_vector = _is_signal(sigma, sd_set, 3)
     s_start = _get_signal_start(sn_vector)
-    s_end = np.sort(len(data) - 1 - _get_signal_start(sn_vector[::-1]))
+    s_end = np.sort(len(data) - 1 - _get_signal_start((sn_vector[::-1])))
     r = _get_temp_baseline(data, s_start, s_end, 7)
     baseline = _smooth(r, 60)
     return data - baseline


Test patch:
diff --git a/nmrglue/analysis/tests/test_analysis_integration.py b/nmrglue/analysis/tests/test_analysis_integration.py
index 4660e1a..81e1301 100644
--- a/nmrglue/analysis/tests/test_analysis_integration.py
+++ b/nmrglue/analysis/tests/test_analysis_integration.py
@@ -58,7 +58,7 @@ def test_1d_integrate():
     assert_array_almost_equal(values_2,  [0.5 * 3, 3.])
 
 
-def test_1d_integrate_with_noise():
+def test_1d_integrate_withnoise():
     """ Test integration of synthetic 1D data with two peaks and noise"""
 
     # seed random number
@@ -88,11 +88,11 @@ def test_1d_integrate_with_noise():
     assert abs(integrate(data, uc, limits) - 1) <= max_error
 
     # Test renormalization of norms
-    results = integrate(data, uc, ((4, 6), (7, 9)), noise_limits=(1, 2),
+    resutls = integrate(data, uc, ((4, 6), (7, 9)), noise_limits=(1, 2),
                         norm_to_range=1)
 
     # Test renormalization of values.
-    assert abs(results[0, 0] - 0.5) <= max_error
+    assert abs(resutls[0, 0] - 0.5) <= max_error
 
 
 def test_1d_ndintegrate():
diff --git a/tests/test_peakpick.py b/tests/test_peakpick.py
new file mode 100644
index 0000000..5895715
--- /dev/null
+++ b/tests/test_peakpick.py
@@ -0,0 +1,245 @@
+import nmrglue as ng
+import numpy as np
+import pytest
+
+
+_ONE_D_PEAKS = {
+    "shape": (1024,),
+    "positions": [100, 200, 300, 400, 500],
+    "lws": [10, 20, 10, 20, 10],
+    "amps": [100, 200, 300, 300, 150],
+    "vparams": [0.2, 0.4, 0.2, 0.5, 1.0],
+}
+
+_TWO_D_PEAKS = {
+    "shape": (512, 128),
+    "positions": [(100, 100), (200, 53), (300, 110), (400, 75)],
+    "lws": [(10, 5), (20, 5), (10, 5), (20, 5)],
+    "amps": [100, 200, 300, 300],
+    "vparams": [(0.2, 0.2), (0.4, 0.4), (0.2, 0.2), (0.5, 0.4)],
+}
+
+
+_THREE_D_PEAKS = {
+    "shape": (256, 64, 64),
+    "positions": [(150, 24, 22), (200, 10, 50), (210, 50, 10), (77, 30, 15)],
+    "lws": [(5, 5, 5), (3, 8, 5), (7, 5, 5), (7, 6, 5)],
+    "amps": [100, 200, 300, 300],
+    "vparams": [(0.2, 0.2, 0.3), (0.5, 0.4, 0.4), (0.1, 0.2, 0.2), (0.3, 0.5, 0.4)],
+}
+
+
+def _generate_1d_data(shape, positions, lws, amps, vparams):
+    """
+    Generates a test 1d dataset with multiple peaks
+
+    Parameters
+    ----------
+    shape : Iterable[int]
+        shape of the numpy array to be created
+    positions : Iterable[int]
+        a list of the positions of the peaks
+    lws : Iterable[float|int]
+        a list of linewidths for each peak
+    amps : Iterable[float|int]
+        a list of amplitudes for each peak
+    vparams : Iterable[float|int]
+        a list of list containing the eta parameter
+        for the pseudo voigt lineshape
+
+    Returns
+    -------
+    numpy.ndarray
+        simulated one-d dataset
+
+    """
+    data = ng.linesh.sim_NDregion(
+        shape=shape,
+        lineshapes=["pv"],
+        params=[[(pos, lws, vp)] for pos, lws, vp in zip(positions, lws, vparams)],
+        amps=amps,
+    )
+    return data
+
+
+def _generate_2d_data(dataset):
+    """
+    Generates a test 2d dataset with multiple peaks
+
+    Parameters
+    ----------
+    shape : Iterable[Iterable[int, int]]
+        shape of the numpy array to be created
+    positions : Iterable[Iterable[int, int]]
+        a list of list of two positions for each peak
+    lws : Iterable[Iterable[float, float]]
+        a list of list of two linewidths for each peak
+    amps : Iterable[Iterable[float, float]]
+        a list of list of two amplitudes for each peak
+    vparams : Iterable[Iterable[float, float]]
+        a list of list containing 2 values for the
+        eta parameter for the pseud-voigt lineshape
+
+    Returns
+    -------
+    numpy.ndarray
+        simulated two-d dataset
+
+    """
+    params = []
+    for i in range(len(dataset["positions"])):
+        d = [[], []]
+        for j in range(2):
+            d[j] = (
+                dataset["positions"][i][j],
+                dataset["lws"][i][j],
+                dataset["vparams"][i][j],
+            )
+
+        params.append(d)
+
+    data = ng.linesh.sim_NDregion(
+        shape=(512, 128), lineshapes=["pv", "pv"], params=params, amps=dataset["amps"]
+    )
+
+    return data
+
+
+def _generate_3d_data(dataset):
+    """
+    Generates a test 3d dataset with multiple peaks
+
+    Parameters
+    ----------
+    shape : Iterable[Iterable[int, int, int]]
+        shape of the numpy array to be created
+    positions : Iterable[Iterable[int, int, int]]
+        a list of list of three positions for each peak
+    lws : Iterable[Iterable[float, float, float]]
+        a list of list of three linewidths for each peak
+    amps : Iterable[Iterable[float, float, float]]
+        a list of list of three amplitudes for each peak
+    vparams : Iterable[Iterable[float, float, float]]
+        a list of list containing three values for the
+        eta parameter for the pseud-voigt lineshape
+
+    Returns
+    -------
+    numpy.ndarray
+        simulated three-d dataset
+
+    """
+    params = []
+    for i in range(len(dataset["positions"])):
+        d = [[], [], []]
+        for j in range(3):
+            d[j] = (
+                dataset["positions"][i][j],
+                dataset["lws"][i][j],
+                dataset["vparams"][i][j],
+            )
+
+        params.append(d)
+
+    data = ng.linesh.sim_NDregion(
+        shape=(256, 64, 64),
+        lineshapes=["pv", "pv", "pv"],
+        params=params,
+        amps=dataset["amps"],
+    )
+
+    return data
+
+
+def _test_1d(dataset, algorithm, msep=None, rtol=1):
+    """test 1d peak picking"""
+
+    data = _generate_1d_data(**dataset)
+    peaks = ng.peakpick.pick(data, pthres=50, algorithm=algorithm, msep=msep)
+    assert np.allclose(peaks.X_AXIS, dataset["positions"], rtol=1)
+
+
+def _test_2d(dataset, algorithm, msep=None, rtol=1):
+    """test 2d peak picking"""
+
+    data = _generate_2d_data(dataset)
+    peaks = ng.peakpick.pick(data, pthres=50, algorithm=algorithm, msep=msep)
+    assert np.allclose(peaks.X_AXIS, [i[1] for i in dataset["positions"]], rtol=1)
+    assert np.allclose(peaks.Y_AXIS, [i[0] for i in dataset["positions"]], rtol=1)
+
+
+def _test_3d(dataset, algorithm, msep=None, rtol=1):
+    """test 3d peak picking"""
+
+    data = _generate_3d_data(dataset)
+    peaks = ng.peakpick.pick(data, pthres=50, algorithm=algorithm, msep=msep)
+
+    assert np.allclose(
+        sorted(peaks.X_AXIS), sorted([i[2] for i in dataset["positions"]]), rtol=rtol
+    )
+    assert np.allclose(
+        sorted(peaks.Y_AXIS), sorted([i[1] for i in dataset["positions"]]), rtol=rtol
+    )
+    assert np.allclose(
+        sorted(peaks.Z_AXIS), sorted([i[0] for i in dataset["positions"]]), rtol=rtol
+    )
+
+
+@pytest.mark.fast
+def test_1d_connected():
+    _test_1d(dataset=_ONE_D_PEAKS, algorithm="connected")
+
+
+@pytest.mark.fast
+def test_1d_downward():
+    _test_1d(dataset=_ONE_D_PEAKS, algorithm="downward")
+
+
+@pytest.mark.fast
+def test_1d_thres():
+    _test_1d(dataset=_ONE_D_PEAKS, algorithm="thres", msep=(5,))
+
+
+@pytest.mark.fast
+def test_1d_thres_fast():
+    _test_1d(dataset=_ONE_D_PEAKS, algorithm="thres-fast", msep=(5,))
+
+
+@pytest.mark.fast
+def test_2d_connected():
+    _test_2d(dataset=_TWO_D_PEAKS, algorithm="connected")
+
+
+@pytest.mark.fast
+def test_2d_downward():
+    _test_2d(dataset=_TWO_D_PEAKS, algorithm="downward")
+
+
+@pytest.mark.fast
+def test_2d_thres():
+    _test_2d(dataset=_TWO_D_PEAKS, algorithm="thres", msep=(5, 5))
+
+
+@pytest.mark.fast
+def test_2d_thres_fast():
+    _test_2d(dataset=_TWO_D_PEAKS, algorithm="thres-fast", msep=(5, 5))
+
+
+@pytest.mark.fast
+def test_3d_connected():
+    _test_3d(dataset=_THREE_D_PEAKS, algorithm="connected")
+
+
+@pytest.mark.fast
+def test_3d_downward():
+    _test_3d(dataset=_THREE_D_PEAKS, algorithm="downward")
+
+
+@pytest.mark.fast
+def test_3d_thres():
+    _test_3d(dataset=_THREE_D_PEAKS, algorithm="thres", msep=(2, 2, 2))
+
+
+@pytest.mark.fast
+def test_3d_thres_fast():
+    _test_3d(dataset=_THREE_D_PEAKS, algorithm="thres-fast", msep=(2, 2, 2))