fast, non-Cython implementation for correlate_sparse

[grlee77]

Description

This PR replaces the recently added correlate_sparse Cython version with a simple slicing-based variant. This was done so that this function can be used with other array types implementing the NumPy API, such as CuPy. There should be no change in behavior or API.

A secondary benefit is that this also profiled modestly faster on 2d and 3d cases. See the threshold_sauvola benchmarks for v0.18.1 vs. this PR below. I tried out varying sparsity levels as well and wasn't able to find a case where it is worse than the current implementation.

Benchmark Results

v0.18.1
· Discovering benchmarks
· Running 2 total benchmarks (1 commits * 1 environments * 2 benchmarks)
[ 0.00%] ·· Benchmarking existing-py_home_lee8rx_miniconda3_envs_skimage_dev_bin_python3.7
[ 25.00%] ··· benchmark_filters.ThresholdSauvolaSuite.time_sauvola 252±0ms
[ 50.00%] ··· benchmark_filters.ThresholdSauvolaSuite.time_sauvola_3d 458±0ms

This PR
· Discovering benchmarks
· Running 2 total benchmarks (1 commits * 1 environments * 2 benchmarks)
[ 0.00%] ·· Benchmarking existing-py_home_lee8rx_miniconda3_envs_skimage_dev_bin_python3.7
[ 25.00%] ··· benchmark_filters.ThresholdSauvolaSuite.time_sauvola 188±0ms
[ 50.00%] ··· benchmark_filters.ThresholdSauvolaSuite.time_sauvola_3d 387±0ms

Checklist

• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Benchmark in ./benchmarks, if your changes aren't covered by an
  existing benchmark
• Unit tests
• Clean style in the spirit of PEP8

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.

[rfezzani]

Thank you @grlee77, I like the idea of reducing code base 😉

[rfezzani]

Suggested change

	if (0, ) * kernel.ndim not in indices:
	if kernel.reshape(-1)[0] == 0:

[grlee77]

I don't think this is equivalent? What I want is if the sparse kernel does not have an entry at the (0,) * ndim corner location, then an entry should be added so that there will not be a shift in the output of _correlate_sparse. (Internal to _correlate_sparse, the output array is initialized from this first point, which needs to be the corner location in order to avoid an unexpected spatial shift in the output).

See if e4cff50 makes things any clearer

[rfezzani]

Isn't (0,) * ndim the index of the first entry of kernel? kernel.reshape(-1)[0] is the first entry of kernel. If it is 0 then (0,) * ndim is not in indices.

[grlee77]

Yes, I agree now that what you suggest is equivalent and would be fine with changing to that.

[grlee77]

On closer look, I would prefer to leave this check as-is so that it avoids a device->host transfer if using a CuPy backend. (For CuPy, kernel.reshape(-1)[0] will give a 0-dimensional cupy.ndarray rather than a host scalar and comparing to 0 requires overhead of transferring that value from the GPU back to the host to do the comparison)

[pep8speaks]

Hello @grlee77! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2021-01-07 13:19:11 UTC

[emmanuelle]

Thanks @grlee77 ! I had the same comment as @rfezzani with respect to assertion vs error but it's not in a function of the public API and the public-API function calling the helper function will not cause problems, so this should not block the PR. Also, I noted that correlate_sparse is in the public API but it's not use in any gallery example, I'm not sure it is very discoverable for users (but once again this should not block the PR).

[jni]

Just some minor comments from me. Thanks @grlee77, very nice!

[jni]

Given that one of the goals of this PR is to work with alternative array implementations, wouldn't it be safer to just assume that we need a copy here? ie we should aim to target the Array API rather than just NumPy and CuPy?

[grlee77]

I am fine with that and don't think it will substantially reduce performance, but given the reliance on scipy.ndimage elsewhere in the library for many things related to filtering/morphology/labeling (as well as usage of FFTs or sparse matrices in various places), we will already not be able to realistically do many things purely with the Array API. I think NumPy+SciPy, Dask+dask-image and CuPy are the only implementations that provide most of the needed features currently.

[jni]

Well, dask arrays don't have .flags either...

In [1]: arr = da.random.random((5,))
In [2]: arr.flags
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-2-c469d1119f64> in <module>
----> 1 arr.flags

AttributeError: 'Array' object has no attribute 'flags'

[jni]

Is this faster than -1 * v? I presume yes — then maybe add a comment above the start of the if/else to say we optimise to avoid multiplication/copy for these values?

[grlee77]

I suspected it may be faster, but hadn't timed it. I tried just now and it was faster on very small arrays, but the same on large arrays, so probably not important to have it as a special case

import numpy as np
a = np.ones((4096, 4096))

%timeit (-a)
22.6 ms ± 130 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit (-1 * a)
22.6 ms ± 33.4 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

a = np.ones((16, 16))

%timeit (-a)
400 ns ± 2.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%timeit (-1 * a)
913 ns ± 3.76 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

[grlee77]

The same relative performance pattern occurs for CuPy as reported for NumPy above (it is necessary to synchronize to measure the actual time since CUDA kernels run asynchronously, so you will only measure the CPU overhead without that)

import cupy as cp
a = cp.ones((4096, 4096))
d = cp.cuda.Device()

%timeit (-a); d.synchronize()
701 µs ± 621 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit (-1 * a); d.synchronize()
705 µs ± 713 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

a = cp.ones((16, 16))

%timeit (-a); d.synchronize()
8.89 µs ± 86.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit (-1 * a); d.synchronize()
12.8 µs ± 228 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

This also demonstrates how CuPy is much faster when working with large arrays, but NumPy is much faster on tiny arrays because it doesn't the kernel launch overhead.

[jni]

Coolio, let's remove that special case then?

[jni]

As a minor point, from an API niceness standpoint I'd rather kernel_indices and kernel_values be two separate arrays that are zipped inside the function. Not critical though.

[grlee77]

The CI failures seen here should be resolved by #5175

[stefanv]

Thanks, Greg!