Translating a 1d convolution from Python to Julia

General Usage
, ,
1

Hi all,

I am trying to port some code written in Python. It performs a 1d convolution on a 3d array. The Python and Julia code output arrays with different dimensions and values. Does anyone know the equivalent code on Julia?

Julia

using DSP
x = [0 0 1 2 3 5 0 0; 0 0 4 5 6 5 0 0 ;;; 0 0 2 3 4 5 0 0 ; 0 0 5 6 7 5 0 0]
z = [5,6,7,10]
    conv(x, z)

output:

5×8×2 Array{Int64, 3}:
[:, :, 1] =
 0  0   5  10  15  25  0  0
 0  0  26  37  48  55  0  0
 0  0  31  44  57  65  0  0
 0  0  38  55  72  85  0  0
 0  0  40  50  60  50  0  0

[:, :, 2] =
 0  0  10  15  20  25  0  0
 0  0  37  48  59  55  0  0
 0  0  44  57  70  65  0  0
 0  0  55  72  89  85  0  0
 0  0  50  60  70  50  0  0

Python

import scipy
x1 = np.array([[0,0,1,2,3,5,0,0],[0,0,4,5,6,5,0,0]])
x2 = np.array([[0,0,2,3,4,5,0,0],[0,0,5,6,7,5,0,0]])
x = np.dstack((x1, x2))
z = [5,6,7,10]
scipy.ndimage.convolve1d(x, z, axis=1, mode='constant', cval=0)

output:

array([[[  5,  10],
        [ 16,  27],
        [ 34,  52],
        [ 67,  90],
        [ 71,  88],
        [ 65,  75],
        [ 50,  50],
        [  0,   0]],

       [[ 20,  25],
        [ 49,  60],
        [ 88, 106],
        [136, 159],
        [122, 139],
        [ 95, 105],
        [ 50,  50],
        [  0,   0]]])
2

I suggest determining how the convolve1d is defined in numpy.
Do the same for the conv operation in DSP.
Then generate the correct calls in Julia.

In practice, I’ve found the simplest is to apply the operations to simple
arrays and see what happens by inspection. You can do that in both
scipy and julia REPLs.

Things to look for

  • Expanding dimensions versus circular correlation or implicit padding
  • Dimension order and index offsets
  • Differences in print or display representations
  • Applying a 1-D operation along a specific dimension of a multi-dimensional data array
  • You julia and numpy versions (I think repeated ;;; inside [] changed in recent versions of julia

N.B. I haven’t used scipy/numpy and not used the DSP package in Julia.

1 Like
3

Thank you. This is good advice. I read through the documentation and code to the best of my ability and tried to discern patterns with simpler examples. The problem is that the behavior of the functions become more divergent as the number of dimensions increases.

The scipy documentation provides the following example:

convolve1d([2, 8, 0, 4, 1, 9, 9, 0], [1, 3], cval=0)
array([14, 24,  4, 13, 12, 36, 27,  0])

In Julia, the code produces:

conv([2, 8, 0, 4, 1, 9, 9, 0], [1, 3])
9-element Vector{Int64}:
  2
 14
 24
  4
 13
 12
 36
 27
  0

The results are similar except for an additional 2 in the Julia output. In my first example involving a 3d array, the arrays differ in ways that I have yet to understand.

4

The following produces the same result:

using DSP
x = [0 0 1 2 3 5 0 0; 0 0 4 5 6 5 0 0 ;;; 0 0 2 3 4 5 0 0 ; 0 0 5 6 7 5 0 0]
z = [5,6,7,10]

n, m = length(z)÷2, size(x,2)
r0 = conv.(eachslice(x,dims=1), (z,))
r = [view(u,n+1:n+m,:) for u in r0]

julia> r[1]
8×2 Matrix{Int64}:
  5  10
 16  27
 34  52
 67  90
 71  88
 65  75
 50  50
  0   0

julia> r[2]
8×2 Matrix{Int64}:
  20   25
  49   60
  88  106
 136  159
 122  139
  95  105
  50   50
   0    0

The result of the DSP convolution must be cut off at the beginning and end, so that its output has the same length as the input.

1 Like
5

Thank you for your help! I will study this in more detail tomorrow.

6

Adding view() inside the comprehension above provides better performance (NB: edited the post), but another way to trim the convolution output is:

r = selectdim.(r0, 1, (n+1:n+m,))
1 Like
7

Thanks for your help! I will mark your answer as a solution. I will provide an update if I am clever enough to find a more efficient solution (I suspect that convolve1d in scipy is more efficient because it is not manipulating intermediate arrays). But in either case, having a correct solution is very helpful!

1 Like