Skip to content

Commit 8dac0f6

Browse files
gjbexgjbex
committed
Add illustration of using numpy arrays directly in Cython
1 parent 3b03a76 commit 8dac0f6

File tree

9 files changed

+379
-0
lines changed

9 files changed

+379
-0
lines changed

‎source-code/cython/Numpy/Convolution/Makefile‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
VERSION = cpython-311-x86_64-linux-gnu
2+
CONVOLUTION_LIB = convolution_cython.$(VERSION).so
3+
CONVOLUTION_INDEXED_LIB = convolution_cython_indexed.$(VERSION).so
4+
CONVOLUTION_NO_CHECKS_LIB = convolution_cython_no_checks.$(VERSION).so
5+
6+
all: $(CONVOLUTION_LIB) $(CONVOLUTION_INDEXED_LIB) $(CONVOLUTION_NO_CHECKS_LIB)
7+
8+
$(CONVOLUTION_LIB): convolution.pyx
9+
python setup.py build_ext --inplace
10+
11+
$(CONVOLUTION_LIB): convolution_indexed.pyx
12+
python setup.py build_ext --inplace
13+
14+
$(CONVOLUTION_LIB): convolution_no_checks.pyx
15+
python setup.py build_ext --inplace
16+
17+
clean:
18+
python setup.py clean
19+
$(RM) $(wildcard *.c) $(wildcard *.so)
20+
$(RM) -r build

‎source-code/cython/Numpy/Convolution/README.md‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
# Convolution
2+
3+
Illustration of how to improve perfomance when using numpy arrays in Cython
4+
by doing indexing right. The algorithm used is a naive implementation of
5+
convolution.
6+
7+
8+
## What is it?
9+
10+
1. `convolution.py`: Python implementation of the algorithm (uses numpy).
11+
1. `convolution.pyx`: Cython implementation that adds types.
12+
1. `convolution_indexed.pyx`: Cython implementation that adds
13+
numpy array element types and indexing information.
14+
1. `convolution_no_checks.pyx`: Cython implementation that adds
15+
function decorator to avoid index checks for numpy arrays.
16+
1. `setup.py`: script to build the extensions.
17+
1. `Makefile`: make file to conveniently build the extensions.
18+
1. `driver.py`: Python scripts to run benchmark tests on the various
19+
implementations.

‎source-code/cython/Numpy/Convolution/convolution.py‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
import numpy as np
2+
3+
4+
def convolve(input, kernel):
5+
'''Naive convolution of matrices input and kernel
6+
7+
Parameters
8+
----------
9+
input : numpy.ndarray
10+
input matrix
11+
kernel : numpy.ndarray
12+
filter kernel matrix, dimensions must be odd
13+
14+
Returns
15+
-------
16+
output : numpy.ndarray
17+
output matrix, it is not cropped
18+
19+
Raises
20+
------
21+
ValueError
22+
if dimensions of kernel are not odd
23+
'''
24+
if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
25+
raise ValueError("Only odd dimensions on filter supported")
26+
# s_mid and t_mid are number of pixels between the center pixel
27+
# and the edge, ie for a 5x5 filter they will be 2.
28+
#
29+
# The output size is calculated by adding s_mid, t_mid to each
30+
# side of the dimensions of the input image.
31+
s_mid = kernel.shape[0] // 2
32+
t_mid = kernel.shape[1] // 2
33+
x_max = input.shape[0] + 2 * s_mid
34+
y_max = input.shape[1] + 2 * t_mid
35+
# Allocate result image.
36+
output = np.zeros([x_max, y_max], dtype=input.dtype)
37+
# Do convolution
38+
for x in range(x_max):
39+
for y in range(y_max):
40+
# Calculate pixel value for h at (x,y). Sum one component
41+
# for each pixel (s, t) of the filter kernel.
42+
s_from = max(s_mid - x, -s_mid)
43+
s_to = min((x_max - x) - s_mid, s_mid + 1)
44+
t_from = max(t_mid - y, -t_mid)
45+
t_to = min((y_max - y) - t_mid, t_mid + 1)
46+
value = 0
47+
for s in range(s_from, s_to):
48+
for t in range(t_from, t_to):
49+
v = x - s_mid + s
50+
w = y - t_mid + t
51+
value += kernel[s_mid - s, t_mid - t]*input[v, w]
52+
output[x, y] = value
53+
return output

‎source-code/cython/Numpy/Convolution/convolution.pyx‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
import numpy as np
2+
cimport numpy as cnp
3+
cnp.import_array()
4+
5+
DTYPE = np.float64
6+
ctypedef cnp.float64_t DTYPE_t
7+
8+
def convolve(cnp.ndarray input, cnp.ndarray kernel):
9+
'''Naive convolution of matrices input and kernel
10+
11+
Parameters
12+
----------
13+
input : numpy.ndarray
14+
input matrix
15+
kernel : numpy.ndarray
16+
filter kernel matrix, dimensions must be odd
17+
18+
Returns
19+
-------
20+
output : numpy.ndarray
21+
output matrix, it is not cropped
22+
23+
Raises
24+
------
25+
ValueError
26+
if dimensions of kernel are not odd
27+
'''
28+
if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
29+
raise ValueError("Only odd dimensions on filter supported")
30+
assert input.dtype == DTYPE and kernel.dtype == DTYPE
31+
# s_mid and t_mid are number of pixels between the center pixel
32+
# and the edge, ie for a 5x5 filter they will be 2.
33+
#
34+
# The output size is calculated by adding s_mid, t_mid to each
35+
# side of the dimensions of the input image.
36+
cdef int s_mid = kernel.shape[0] // 2
37+
cdef int t_mid = kernel.shape[1] // 2
38+
cdef int x_max = input.shape[0] + 2 * s_mid
39+
cdef int y_max = input.shape[1] + 2 * t_mid
40+
# Allocate result image.
41+
cdef cnp.ndarray output = np.zeros([x_max, y_max], dtype=input.dtype)
42+
# Do convolution
43+
cdef int x, y, s, t, v, w
44+
cdef int s_from, s_to, t_from, t_to
45+
cdef DTYPE_t value
46+
for x in range(x_max):
47+
for y in range(y_max):
48+
# Calculate pixel value for h at (x,y). Sum one component
49+
# for each pixel (s, t) of the filter kernel.
50+
s_from = max(s_mid - x, -s_mid)
51+
s_to = min((x_max - x) - s_mid, s_mid + 1)
52+
t_from = max(t_mid - y, -t_mid)
53+
t_to = min((y_max - y) - t_mid, t_mid + 1)
54+
value = 0
55+
for s in range(s_from, s_to):
56+
for t in range(t_from, t_to):
57+
v = x - s_mid + s
58+
w = y - t_mid + t
59+
value += kernel[s_mid - s, t_mid - t]*input[v, w]
60+
output[x, y] = value
61+
return output

‎source-code/cython/Numpy/Convolution/convolution_indexed.pyx‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
import numpy as np
2+
cimport numpy as cnp
3+
cnp.import_array()
4+
5+
DTYPE = np.float64
6+
ctypedef cnp.float64_t DTYPE_t
7+
8+
def convolve(cnp.ndarray[DTYPE_t, ndim=2] input, cnp.ndarray[DTYPE_t, ndim=2] kernel):
9+
'''Naive convolution of matrices input and kernel
10+
11+
Parameters
12+
----------
13+
input : numpy.ndarray
14+
input matrix
15+
kernel : numpy.ndarray
16+
filter kernel matrix, dimensions must be odd
17+
18+
Returns
19+
-------
20+
output : numpy.ndarray
21+
output matrix, it is not cropped
22+
23+
Raises
24+
------
25+
ValueError
26+
if dimensions of kernel are not odd
27+
'''
28+
if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
29+
raise ValueError("Only odd dimensions on filter supported")
30+
assert input.dtype == DTYPE and kernel.dtype == DTYPE
31+
# s_mid and t_mid are number of pixels between the center pixel
32+
# and the edge, ie for a 5x5 filter they will be 2.
33+
#
34+
# The output size is calculated by adding s_mid, t_mid to each
35+
# side of the dimensions of the input image.
36+
cdef int s_mid = kernel.shape[0] // 2
37+
cdef int t_mid = kernel.shape[1] // 2
38+
cdef int x_max = input.shape[0] + 2 * s_mid
39+
cdef int y_max = input.shape[1] + 2 * t_mid
40+
# Allocate result image.
41+
cdef cnp.ndarray[DTYPE_t, ndim=2] output = np.zeros([x_max, y_max], dtype=input.dtype)
42+
# Do convolution
43+
cdef int x, y, s, t, v, w
44+
cdef int s_from, s_to, t_from, t_to
45+
cdef DTYPE_t value
46+
for x in range(x_max):
47+
for y in range(y_max):
48+
# Calculate pixel value for h at (x,y). Sum one component
49+
# for each pixel (s, t) of the filter kernel.
50+
s_from = max(s_mid - x, -s_mid)
51+
s_to = min((x_max - x) - s_mid, s_mid + 1)
52+
t_from = max(t_mid - y, -t_mid)
53+
t_to = min((y_max - y) - t_mid, t_mid + 1)
54+
value = 0
55+
for s in range(s_from, s_to):
56+
for t in range(t_from, t_to):
57+
v = x - s_mid + s
58+
w = y - t_mid + t
59+
value += kernel[s_mid - s, t_mid - t]*input[v, w]
60+
output[x, y] = value
61+
return output

‎source-code/cython/Numpy/Convolution/convolution_no_checks.pyx‎

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
cimport cython
2+
import numpy as np
3+
cimport numpy as cnp
4+
cnp.import_array()
5+
6+
DTYPE = np.float64
7+
ctypedef cnp.float64_t DTYPE_t
8+
9+
@cython.boundscheck(False)
10+
@cython.wraparound(False)
11+
def convolve(cnp.ndarray[DTYPE_t, ndim=2] input, cnp.ndarray[DTYPE_t, ndim=2] kernel):
12+
'''Naive convolution of matrices input and kernel
13+
14+
Parameters
15+
----------
16+
input : numpy.ndarray
17+
input matrix
18+
kernel : numpy.ndarray
19+
filter kernel matrix, dimensions must be odd
20+
21+
Returns
22+
-------
23+
output : numpy.ndarray
24+
output matrix, it is not cropped
25+
26+
Raises
27+
------
28+
ValueError
29+
if dimensions of kernel are not odd
30+
'''
31+
if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
32+
raise ValueError("Only odd dimensions on filter supported")
33+
assert input.dtype == DTYPE and kernel.dtype == DTYPE
34+
# s_mid and t_mid are number of pixels between the center pixel
35+
# and the edge, ie for a 5x5 filter they will be 2.
36+
#
37+
# The output size is calculated by adding s_mid, t_mid to each
38+
# side of the dimensions of the input image.
39+
cdef int s_mid = kernel.shape[0] // 2
40+
cdef int t_mid = kernel.shape[1] // 2
41+
cdef int x_max = input.shape[0] + 2 * s_mid
42+
cdef int y_max = input.shape[1] + 2 * t_mid
43+
# Allocate result image.
44+
cdef cnp.ndarray[DTYPE_t, ndim=2] output = np.zeros([x_max, y_max], dtype=input.dtype)
45+
# Do convolution
46+
cdef int x, y, s, t, v, w
47+
cdef int s_from, s_to, t_from, t_to
48+
cdef DTYPE_t value
49+
for x in range(x_max):
50+
for y in range(y_max):
51+
# Calculate pixel value for h at (x,y). Sum one component
52+
# for each pixel (s, t) of the filter kernel.
53+
s_from = max(s_mid - x, -s_mid)
54+
s_to = min((x_max - x) - s_mid, s_mid + 1)
55+
t_from = max(t_mid - y, -t_mid)
56+
t_to = min((y_max - y) - t_mid, t_mid + 1)
57+
value = 0
58+
for s in range(s_from, s_to):
59+
for t in range(t_from, t_to):
60+
v = x - s_mid + s
61+
w = y - t_mid + t
62+
value += kernel[s_mid - s, t_mid - t]*input[v, w]
63+
output[x, y] = value
64+
return output

‎source-code/cython/Numpy/Convolution/driver.py‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
#!/usr/bin/env python
2+
3+
# script to run various implementations of the convolution filter
4+
# to compare their performance.
5+
# The script has the following input parameters:
6+
# --input_x: width of the input matrix
7+
# --input_y: height of the input matrix
8+
# --kernel_x: width of the kernel matrix
9+
# --kernel_y: height of the kernel matrix
10+
# --iterations: number of iterations to run
11+
# --implementation: which implementation to run
12+
13+
import argparse
14+
import numpy as np
15+
import timeit
16+
import sys
17+
18+
def main():
19+
parser = argparse.ArgumentParser(description='Run convolution filter')
20+
parser.add_argument('--input_x', type=int, default=100,
21+
help='width of the input matrix')
22+
parser.add_argument('--input_y', type=int, default=100,
23+
help='height of the input matrix')
24+
parser.add_argument('--kernel_x', type=int, default=9,
25+
help='width of the kernel matrix')
26+
parser.add_argument('--kernel_y', type=int, default=9,
27+
help='height of the kernel matrix')
28+
parser.add_argument('--iterations', type=int, default=1,
29+
help='number of iterations to run')
30+
parser.add_argument('--implementation', type=str,
31+
choices=['python', 'cython', 'cython_indexed', 'cython_no_checks'],
32+
default='python',
33+
help='which implementation to run')
34+
parser.add_argument('--check', action='store_true',
35+
help='check the results')
36+
args = parser.parse_args()
37+
38+
# Create input and kernel matrices
39+
input = np.random.rand(args.input_x, args.input_y)
40+
kernel = np.random.rand(args.kernel_x, args.kernel_y)
41+
42+
# Run the requested implementation
43+
if args.implementation == 'python':
44+
from convolution import convolve
45+
elif args.implementation == 'cython':
46+
from convolution_cython import convolve
47+
elif args.implementation == 'cython_indexed':
48+
from convolution_cython_indexed import convolve
49+
elif args.implementation == 'cython_no_checks':
50+
from convolution_cython_no_checks import convolve
51+
52+
# Measure the execution time
53+
time = timeit.timeit(lambda: convolve(input, kernel), number=args.iterations)
54+
print(f'Time: {time/args.iterations} s per iteration, {args.iterations} iterations')
55+
56+
# Check the results
57+
if args.check:
58+
from convolution import convolve as convolve_python
59+
result_python = convolve_python(input, kernel)
60+
result = convolve(input, kernel)
61+
if np.allclose(result_python, result):
62+
print('Results are correct')
63+
else:
64+
print('Results are incorrect')
65+
sys.exit(1)
66+
67+
if __name__ == '__main__':
68+
main()

‎source-code/cython/Numpy/Convolution/setup.py‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
#!/usr/bin/env python
2+
3+
from distutils.core import setup
4+
from Cython.Build import cythonize
5+
from pathlib import Path
6+
from distutils.extension import Extension
7+
8+
home_dir = str(Path.home())
9+
10+
extensions = [
11+
Extension(
12+
"convolution_cython",
13+
["convolution.pyx"],
14+
include_dirs=[f'{home_dir}/mambaforge/envs/python_for_hpc/lib/python3.11/site-packages/numpy/core/include/'],
15+
),
16+
Extension(
17+
"convolution_cython_indexed",
18+
["convolution_indexed.pyx"],
19+
include_dirs=[f'{home_dir}/mambaforge/envs/python_for_hpc/lib/python3.11/site-packages/numpy/core/include/'],
20+
),
21+
Extension(
22+
"convolution_cython_no_checks",
23+
["convolution_no_checks.pyx"],
24+
include_dirs=[f'{home_dir}/mambaforge/envs/python_for_hpc/lib/python3.11/site-packages/numpy/core/include/'],
25+
),
26+
]
27+
28+
setup(
29+
ext_modules=cythonize(extensions,
30+
language_level='3str')
31+
)

‎source-code/cython/Numpy/README.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,3 +24,5 @@ to `np.sum`.
2424
the pure Python Cython implementation of the array sum.
2525
1. `Makefile`: make file for building the cython extension and profiling
2626
the three implementations.
27+
1. `Convolution`: illustration of using numpy arrays directly in Cython
28+
code.

0 commit comments

Comments
 (0)