Merge pull request #40 from neuro-ml/dev

Dev

.github/workflows/tests.yml

@@ -61,7 +61,7 @@ jobs:
           pytest tests/test_backend.py -m nonumba
       - name: Generate coverage report
         run: |
-          coverage xml -o reports/coverage-${{ matrix.python-version }}.xml
+          coverage xml -o reports/coverage-${{ matrix.python-version }}.xml --omit=imops/src/_numba_zoom.py
           sed -i -e "s|$MODULE_PARENT/||g" reports/coverage-${{ matrix.python-version }}.xml
           sed -i -e "s|$(echo $MODULE_PARENT/ | tr "/" .)||g" reports/coverage-${{ matrix.python-version }}.xml
 

benchmarks/benchmark_numeric.py

@@ -29,14 +29,14 @@ def time_add_value(self, backend, num_threads):
 
     @discard_arg(2)
     def time_copy(self, backend, num_threads):
-        copy(self.nums1_3d, num_threads, backend)
+        copy(self.nums1_3d, num_threads=num_threads, backend=backend)
 
     def time_full(self, backend, dtype, num_threads):
-        full(self.shape, 42, dtype, num_threads, backend)
+        full(self.shape, 42, dtype, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
     def time_fill_(self, backend, num_threads):
-        fill_(self.empty_3d, 42, num_threads, backend)
+        fill_(self.empty_3d, 42, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
     def peakmem_add_array(self, backend, num_threads):
@@ -48,11 +48,11 @@ def peakmem_add_value(self, backend, num_threads):
 
     @discard_arg(2)
     def peakmem_copy(self, backend, num_threads):
-        copy(self.nums1_3d, num_threads, backend)
+        copy(self.nums1_3d, num_threads=num_threads, backend=backend)
 
     def peakmem_full(self, backend, dtype, num_threads):
-        full(self.shape, 42, dtype, num_threads, backend)
+        full(self.shape, 42, dtype, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
     def peakmem_fill_(self, backend, num_threads):
-        fill_(self.empty_3d, 42, num_threads, backend)
+        fill_(self.empty_3d, 42, num_threads=num_threads, backend=backend)

imops/__version__.py

@@ -1 +1 @@
-__version__ = '0.8.3'
+__version__ = '0.8.4'

imops/box.py

@@ -5,6 +5,9 @@
 
 import numpy as np
 
+# for backward compatibility
+from .utils import build_slices  # noqa: F401
+
 
 # Immutable numpy array
 Box = np.ndarray

imops/morphology.py

@@ -3,7 +3,12 @@
 
 import numpy as np
 from scipy.ndimage import generate_binary_structure
-from skimage.morphology import binary_dilation as scipy_binary_dilation, binary_erosion as scipy_binary_erosion
+from skimage.morphology import (
+    binary_closing as scipy_binary_closing,
+    binary_dilation as scipy_binary_dilation,
+    binary_erosion as scipy_binary_erosion,
+    binary_opening as scipy_binary_opening,
+)
 
 from .backend import BackendLike, Cython, Scipy, resolve_backend
 from .box import add_margin, box_to_shape, mask_to_box, shape_to_box
@@ -14,7 +19,7 @@
     _binary_erosion as cython_fast_binary_erosion,
 )
 from .src._morphology import _binary_dilation as cython_binary_dilation, _binary_erosion as cython_binary_erosion
-from .utils import composition_args, morphology_composition_args, normalize_num_threads
+from .utils import morphology_composition_args, normalize_num_threads
 
 
 def morphology_op_wrapper(
@@ -42,8 +47,12 @@ def wrapped(
 
         if output is None:
             output = np.empty_like(image, dtype=bool)
+        elif boxed:
+            raise ValueError('`boxed==True` is incompatible with provided `output`')
         elif output.shape != image.shape:
             raise ValueError('Input image and output image shapes must be the same.')
+        elif output.dtype != bool:
+            raise ValueError(f'Output image must have `bool` dtype, got {output.dtype}.')
         elif not output.data.c_contiguous:
             # TODO: Implement morphology for `output` of arbitrary layout
             raise ValueError('`output` must be a C-contiguous array.')
@@ -53,7 +62,7 @@ def wrapped(
         if backend.name == 'Scipy':
             if boxed:
                 raise ValueError('`boxed==True` is incompatible with "Scipy" backend.')
-            output = src_op(image, footprint)
+            src_op(image, footprint, out=output)
 
             return output
 
@@ -63,7 +72,7 @@ def wrapped(
                 "Falling back to scipy's implementation.",
                 stacklevel=3,
             )
-            output = backend2src_op[Scipy()](image, footprint)
+            backend2src_op[Scipy()](image, footprint, out=output)
 
             return output
 
@@ -95,7 +104,7 @@ def wrapped(
         if n_dummy:
             output = output[(0,) * n_dummy]
 
-        return output.astype(bool, copy=False)
+        return output
 
     return wrapped
 
@@ -244,7 +253,7 @@ def binary_erosion(
 _binary_closing = morphology_op_wrapper(
     'binary_closing',
     {
-        Scipy(): composition_args(scipy_binary_erosion, scipy_binary_dilation),
+        Scipy(): scipy_binary_closing,
         Cython(fast=False): morphology_composition_args(cython_binary_erosion, cython_binary_dilation),
         Cython(fast=True): morphology_composition_args(cython_fast_binary_erosion, cython_fast_binary_dilation),
     },
@@ -297,7 +306,7 @@ def binary_closing(
 _binary_opening = morphology_op_wrapper(
     'binary_opening',
     {
-        Scipy(): composition_args(scipy_binary_dilation, scipy_binary_erosion),
+        Scipy(): scipy_binary_opening,
         Cython(fast=False): morphology_composition_args(cython_binary_dilation, cython_binary_erosion),
         Cython(fast=True): morphology_composition_args(cython_fast_binary_dilation, cython_fast_binary_erosion),
     },

imops/numeric.py

@@ -258,7 +258,7 @@ def full(
     backend: BackendLike = None,
 ) -> np.ndarray:
     """
-    Return a new array of given shape and type, filled with `fill_value`.
+    Return a new array of given shape and dtype, filled with `fill_value`.
 
     Uses a fast parallelizable implementation for fp16-32-64 and int16-32-64 inputs and ndim <= 4.
 
@@ -284,10 +284,10 @@ def full(
     >>> x = full((2, 3, 4), 1.5, dtype=int)  # same as np.ones((2, 3, 4), dtype=int)
     >>> x = full((2, 3, 4), 1, dtype='uint16')  # will fail because of unsupported uint16 dtype
     """
-    nums = np.empty(shape, dtype=dtype, order=order)
+    dtype = dtype or np.array(fill_value).dtype
 
-    if dtype is not None:
-        fill_value = nums.dtype.type(fill_value)
+    nums = np.empty(shape, dtype=dtype, order=order)
+    fill_value = nums.dtype.type(fill_value)
 
     fill_(nums, fill_value, num_threads, backend)
 

imops/testing.py

@@ -12,13 +12,12 @@ def sk_iradon(xs):
         return np.stack([iradon_(x) for x in xs])
 
 
-def sk_radon(xs, strict=True):
+def sk_radon(xs):
     with warnings.catch_warnings():
         warnings.filterwarnings('ignore', module='numpy')
         warnings.simplefilter('ignore', DeprecationWarning)
         warnings.simplefilter('ignore', np.VisibleDeprecationWarning)
-        if strict:
-            warnings.filterwarnings('error', '.*image must be zero.*', module='skimage')
+        warnings.filterwarnings('error', '.*image must be zero.*', module='skimage')
         return np.stack([radon_(x) for x in xs])
 
 

imops/utils.py

@@ -151,14 +151,6 @@ def broadcast_to_axis(axis: AxesLike, *arrays: AxesParams):
     return tuple(np.repeat(x, len(axis) // len(x), 0) for x in arrays)
 
 
-# TODO: come up with a better name
-def composition_args(f: Callable, g: Callable) -> Callable:
-    def inner(*args):
-        return f(g(*args), *args[1:])
-
-    return inner
-
-
 def morphology_composition_args(f, g) -> Callable:
     def wrapper(
         image: np.ndarray,

tests/test_interp1d.py

@@ -95,6 +95,34 @@ def test_length_inequality_exception(backend):
         interp1d(x, y, axis=0, backend=backend)
 
 
+def test_nans(backend):
+    if backend.name == 'Scipy':
+        return
+
+    x = np.array([0, 1, 2])
+    y = np.array([np.inf, -np.inf, np.inf])
+
+    with pytest.raises(RuntimeError):
+        interp1d(x, y, axis=0, fill_value=0, backend=backend)(x / 2)
+
+    x = np.array([0, 1, 2, 3, 4, 5])
+    y = np.array([np.inf, 0, 1, 2, -np.inf, np.inf])
+
+    with pytest.raises(RuntimeError):
+        interp1d(x, y, axis=0, fill_value=0, backend=backend)(x)
+
+    y = np.array([np.inf, 0, 1, np.inf, -np.inf, np.inf])
+
+    allclose(
+        interp1d(x, y, axis=0, fill_value=0, backend=backend)(x / 2),
+        np.array([np.inf, np.inf, np.inf, 0.5, 1, np.inf]),
+    )
+    allclose(
+        interp1d(x, -y, axis=0, fill_value=0, backend=backend)(x / 2),
+        np.array([-np.inf, -np.inf, -np.inf, -0.5, -1, -np.inf]),
+    )
+
+
 def test_extrapolation(backend):
     for i in range(n_samples):
         shape = np.random.randint(16, 64, size=np.random.randint(1, 4))
@@ -152,7 +180,15 @@ def test_stress(backend):
         old_locations = np.random.randn(shape[axis])
         new_locations = np.random.randn(np.random.randint(shape[axis] // 2, shape[axis] * 2))
 
-        out = interp1d(old_locations, inp, axis=axis, bounds_error=False, fill_value=0, backend=backend)(new_locations)
+        out = interp1d(
+            old_locations,
+            inp,
+            axis=axis,
+            copy=np.random.binomial(1, 0.5),
+            bounds_error=False,
+            fill_value=0,
+            backend=backend,
+        )(new_locations)
         desired_out = scipy_interp1d(old_locations, inp, axis=axis, bounds_error=False, fill_value=0)(new_locations)
 
         allclose(out, desired_out, rtol=1e-6, err_msg=f'{i, shape}')

tests/test_measure.py

@@ -338,7 +338,7 @@ def test_labeled_center_of_mass(backend, dtype, label_dtype):
             else np.random.choice(np.array([[False], [True], [False, True], [True, False]], dtype=object))
         )
 
-        out = center_of_mass(inp, labels, index, backend=backend)
+        out = center_of_mass(inp, labels, index, num_threads=1, backend=backend)
         desired_out = scipy_center_of_mass(inp, labels, index)
 
         for x, y in zip(out, desired_out):

tests/test_morphology.py

@@ -150,12 +150,12 @@ def take_by_coords(array, coords):
             box_pos = np.asarray([np.random.randint(0, s - bs + 1) for bs, s in zip(box_size, shape)])
             box_coord = np.array([box_pos, box_pos + box_size])
             inp = np.random.binomial(1, 0.7, box_size)
-            inp = restore_crop(inp, box_coord, shape, 0)
+            inp = restore_crop(inp, box_coord, shape, 0).astype(bool)
         else:
-            inp = np.random.binomial(1, 0.5, shape)
+            inp = np.random.binomial(1, 0.5, shape).astype(bool)
 
         footprint_shape = footprint_shape_modifier(np.random.randint(1, 4, size=inp.ndim))
-        footprint = np.random.binomial(1, 0.5, footprint_shape) if np.random.binomial(1, 0.5, 1) else None
+        footprint = np.random.binomial(1, 0.5, footprint_shape) if np.random.binomial(1, 0.5) else None
 
         if backend == Scipy() and boxed:
             with pytest.raises(ValueError):
@@ -173,9 +173,15 @@ def take_by_coords(array, coords):
             return
 
         desired_out = sk_op(inp, footprint)
+        output = np.empty_like(inp)
+
+        if np.random.binomial(1, 0.5) or boxed:
+            output = imops_op(inp, footprint, backend=backend, boxed=boxed)
+        else:
+            imops_op(inp, footprint, output=output, backend=backend, boxed=boxed)
 
         assert_eq(
-            imops_op(inp, footprint, backend=backend, boxed=boxed),
+            output,
             desired_out,
             err_msg=f'{i, shape, footprint, box_coord if boxed else None}',
         )

tests/test_numeric.py

@@ -272,8 +272,10 @@ def sample_value(dtype):
         shape = np.random.randint(32, 64, size=np.random.randint(1, 5))
         fill_value = sample_value(np.zeros(1, dtype=dtype).dtype)
 
-        nums = full(shape, fill_value, dtype, num_threads=num_threads, backend=backend)
-        desired_nums = np.full(shape, fill_value, dtype if np.random.binomial(1, 0.5) else None)
+        dtype_or_none = dtype if np.random.binomial(1, 0.5) else None
+
+        nums = full(shape, fill_value, dtype_or_none, num_threads=num_threads, backend=backend)
+        desired_nums = np.full(shape, fill_value, dtype_or_none)
 
         if dtype in ('int16', 'int32', 'int64'):
             assert_eq(nums, desired_nums)

tests/test_utils.py

@@ -1,10 +1,24 @@
 import os
 from unittest import mock
 
+import numpy as np
 import pytest
 
 from imops.backend import Cython
-from imops.utils import build_slices, check_len, imops_num_threads, normalize_num_threads, set_num_threads
+from imops.utils import (
+    broadcast_axis,
+    broadcast_to_axis,
+    build_slices,
+    check_len,
+    imops_num_threads,
+    normalize_num_threads,
+    set_num_threads,
+)
+
+
+assert_eq = np.testing.assert_array_equal
+
+MANY_THREADS = 42069
 
 
 def test_check_len():
@@ -48,3 +62,51 @@ def test_imops_num_threads():
         assert normalize_num_threads(-1, Cython()) == min(10, len(os.sched_getaffinity(0)))
 
     assert normalize_num_threads(-1, Cython()) == len(os.sched_getaffinity(0))
+
+
+@mock.patch.dict(os.environ, {}, clear=True)
+def test_many_threads_warning_os():
+    with pytest.warns(UserWarning):
+        normalize_num_threads(MANY_THREADS, Cython())
+
+
+@mock.patch.dict(os.environ, {'OMP_NUM_THREADS': '2'}, clear=True)
+def test_many_threads_warning_omp():
+    with pytest.warns(UserWarning):
+        normalize_num_threads(MANY_THREADS, Cython())
+
+
+@mock.patch.dict(os.environ, {}, clear=True)
+def test_many_threads_warning_imops():
+    with imops_num_threads(10):
+        with pytest.warns(UserWarning):
+            normalize_num_threads(MANY_THREADS, Cython())
+
+
+def test_broadcast_to_axis():
+    arrays = np.ones((1, 2)), np.ones((3, 4, 5)), np.ones(1), 1
+    axis = [0, 0, 0]
+
+    for x, out in zip((np.ones((3, 2)), np.ones((3, 4, 5)), np.ones(3), np.ones(3)), broadcast_to_axis(axis, *arrays)):
+        assert_eq(x, out)
+
+    with pytest.raises(ValueError):
+        broadcast_to_axis(axis)
+
+    with pytest.raises(ValueError):
+        broadcast_to_axis(None, *arrays)
+
+    with pytest.raises(ValueError):
+        broadcast_to_axis([0, 0], *arrays)
+
+
+def test_broadcast_axis():
+    arrays = np.ones((1, 3)), np.ones((2, 3))
+
+    for out in broadcast_axis([0, 1], 2, *arrays)[1:]:
+        assert_eq(out, np.ones((2, 3)))
+
+    arrays = np.ones((3, 1)), np.ones((2, 3))
+
+    with pytest.raises(ValueError):
+        broadcast_axis([0, 1], 2, *arrays)