All-Python implementation of n-dimensional convolution

Question

To better understand how discrete finite convolution works (read educational purposes) I wrote an all-python implementation of the convolution function. The idea was for it to give the same output as numpy.convolve, including the mode options. I generalized the code so that it functions for n-dimensional convolutions rather than just for 1-dimensional convolutions. This means it more closely matches the behavior of scipy.signal.convolve.

def py_nd_convolve(s, k, mode='full'):
 # All python implementation of n-dimensional scipy.signal.convolve for educational purposes
 # First ensure that one of the arrays fits inside of the other and by convention call this one
 # the kernel "k" and call the larger array the signal "s".
 if all(np.less_equal(k.shape, s.shape)):
 pass
 elif all(np.greater(k.shape, s.shape)):
 s, k = k.astype(float), s.astype(float)
 else:
 raise ValueError(f'Arrays have mismatched shapes: {k.shape} and {s.shape}')
 dim_list = np.arange(s.ndim)
 n_s = np.array(s.shape, dtype=int)
 n_k = np.array(k.shape, dtype=int)
 n_full = n_s + n_k - 1
 out = np.zeros(n_full, dtype=float)
 zero_arr = np.zeros(s.ndim, dtype=int)
 for ind, val in np.ndenumerate(out):
 ind_arr = np.array(ind) # convert index tuple to numpy array so it can be used for arithmetic
 # i_min and i_max specify the allowed range which indices can take on in each dimension
 # to perform the convolution summation. The values for these are non-trivial near the
 # boundaries of the output array.
 i_min = np.maximum(zero_arr, ind_arr - n_k + 1)
 i_max = np.minimum(ind_arr, n_s - 1)
 # create a d-dimensional tuple of slices to slice both the kernel and signal to perform the convolution.
 s_slicer = tuple([slice(i_min[d], i_max[d] + 1) for d in dim_list])
 k_slicer = tuple([slice(ind[d]-i_max[d], ind[d]-i_min[d]+1) for d in dim_list])
 # element-wise product and summation of subset of signal and reversed kernel
 out[tuple(ind)] = np.multiply(s[s_slicer], np.flip(k[k_slicer])).sum()
 # This section of code clips the output data as per user specification to manage boundary effects
 n_min = np.zeros_like(n_full)
 n_max = n_full
 if mode == 'full': # take all possible output points, maximal boundary effects
 pass
 elif mode == 'valid': # eliminate any part of the output which is effect by boundary effects
 n_min = n_k-1
 n_max = n_s
 elif mode == 'same': # clip the output so that it has the same shape as s
 n_min = (n_k-1)/2
 n_max = n_s + (n_k-1)/2
 slicer = [slice(int(n_min[d]), int(n_max[d])) for d in dim_list] # d-dimensional slice tuple for clipping
 return out[tuple(slicer)]

Since this code was for educational purposes I did not desire it to run fast, I namely desired for it to be very readable. As written the code is very very slow. While scipy.signal.convolve can convolve a (200, 200) array with a (50, 50) array in a few milliseconds this code takes a few seconds.

I am not interested in making this code as fast as possible because I know that for that I could simply use the built in functions. I also know that I could possibly perform this operation faster by using an fft approach.

Question

Rather, I am curious if there are any glaring performance issues with the code that can be easily modified so that I can maintain a similar level of readability and logical flow but a higher level of performance.

For example, am I making an egregious numpy mistakes such as copying arrays unnecessarily or misusing memory in a very bad way?

• I've noticed that I could flip the k array outside of the loop but it doesn't look like this will save much time.
• It looks like a big chunk of time is spent on defining s_slicer and k_slicer.

• Jagerber48: Originally asked on Stackoverflow but I think it should have been here originally: stackoverflow.com/questions/57245494/…

Answer 1

Specific suggestions:

1. All t abbr mks th c hard 2 rd. Plz repl abbrs w/full wrds. Think of variable names as comments – they don't matter at all to the computer, but they do matter to anybody reading the code. A good rule of thumb is that if any of your comments explain what is coming up rather than why it was done like that, you can probably replace the comment with more elegant code. I've been working over a year with some colleagues on a medium-size new system, and we've taken this approach to heart: about 0.5% of the lines in Python files have comments on them.
2. Conversions often take a lot of time, especially on a lot of items, and even more especially in a loop. Working with only lists or only tuples might speed up things a bit. In terms of readability, tuples are generally used when there is a fixed number of items, often of different types. Neither of these seem to apply to this code.
3. The contents of the function can probably be split into at least three more for readability:
   1. validation,
   2. mode handling, and
   3. the loop contents.

4. The modes could be related as Enum values:

   class Mode(Enum):
    FULL = "full"
    VALID = "valid"
    SAME = "same"
   def ...(..., mode: Mode = Mode.FULL):
   

General suggestions:

1. black can automatically format your code to be more idiomatic.
2. isort can group and sort your imports automatically.

3. flake8 with a strict complexity limit will give you more hints to write idiomatic Python:

   [flake8]
   max-complexity = 4
   ignore = W503,E203
   

4. I would then recommend adding type hints everywhere and validating them using a strict mypy configuration:

   [mypy]
   check_untyped_defs = true
   disallow_untyped_defs = true
   ignore_missing_imports = true
   no_implicit_optional = true
   warn_redundant_casts = true
   warn_return_any = true
   warn_unused_ignores = true
   

• Jagerber48: Thank you for your recommendations. I rearranged the loop a bit and I do less type conversions. That seems to have helped a lot. Is it standard for me to post the updated code? Two questions. 1) I've been looking up about Enum values but I haven't found any example where they are used to help with function options. Could you create or point me to a basic example? 2) Regarding abbreviations I've kept some variables short because I think of them as variables in equations. I have a science background so lots of math. Do you think things like n_s should be renamed to signal_shape?
• l0b0: You can post another question with the new code - that's welcome and fairly common.
• l0b0: n_s is almost completely meaningless to someone unfamiliar with the code, so signal_shape would almost certainly be better (if that is really what it is). It would help if I first knew what s and k are - to me they look interchangeable, like they be named first and second or target and other, or something.