Merge branch 'development' into feature/update_cython_slides

Author: gjbex

source-code/cython/Exceptions/.gitignore

@@ -0,0 +1,2 @@
+average.c
+average_pure.c

source-code/cython/Exceptions/Makefile

@@ -1,11 +1,16 @@
-VERSION = cpython-34m
+VERSION = cpython-311-x86_64-linux-gnu
 AVERAGE_LIB = average.$(VERSION).so
+AVERAGE_PURE_LIB = average_pure.$(VERSION).so
 
-all: $(AVERAGE_LIB)
+all: $(AVERAGE_LIB) $(AVERAGE_PURE_LIB)
 
 $(AVERAGE_LIB): average.pyx
 	python setup.py build_ext --inplace
 
+$(AVERAGE_PURE_LIB): average_pure.py
+	python setup.py build_ext --inplace
+
 clean:
 	python setup.py clean
-	rm -f average.c $(AVERAGE_LIB)
+	$(RM) average.c average_pure.c $(AVERAGE_LIB) $(AVERAGE_PURE_LIB)
+	$(RM) -r build

source-code/cython/Exceptions/README.md

@@ -2,9 +2,12 @@
 Error handling in Cython code.
 
 ## What is it?
-1. `average.pyx`: code to be compiled using Cython, implements two
+1. `average.pyx`: Cython code to be compiled using Cython, implements two
     functions to compute the average of an array slice, one with, the
     other without error handling.
+1. `average_pure.py`: pure Python code to be compiled using Cython,
+   implements two functions to compute the average of an array slice,
+   one with, the other without error handling.
 1. `setup.py`: Python build script.
 1. `Makefile`: make file to build the extension.
 1. `compute_average.py`: script to load the compiled module and call the

source-code/cython/Exceptions/average_pure.py

@@ -0,0 +1,29 @@
+import cython
+
+
+def average(data, m=0, n=None):
+    if n is None:
+        n = len(data)
+    return _average(memoryview(data), m, n)
+
+def average_no_except(data, m=0, n=None):
+    if n is None:
+        n = len(data)
+    return _average_no_except(memoryview(data), m, n)
+
+@cython.cfunc
+@cython.exceptval(-1.0, check=True)
+def _average(data, m: cython.int=0, n: cython.int=-1) -> cython.double:
+    i: cython.int
+    mean: cython.double = 0.0
+    for i in range(m, n):
+        mean += data[i]
+    return mean/(n - m + 1)
+
+@cython.cfunc
+def _average_no_except(data, m: cython.int=0, n: cython.int=-1) -> cython.double:
+    i: cython.int
+    mean: cython.double = 0.0
+    for i in range(m, n):
+        mean += data[i]
+    return mean/(n - m + 1)

source-code/cython/Exceptions/compute_average.py

@@ -2,24 +2,31 @@
 
 from argparse import ArgumentParser
 import array
-import average
 
-
-size = 10
+DEFAULT_SIZE = 10
 
 arg_parser = ArgumentParser(description='test Cython errors')
 arg_parser.add_argument('--m', type=int, default=0, help='lower bound')
-arg_parser.add_argument('--n', type=int, default=size, help='upper bound')
+arg_parser.add_argument('--n', type=int, default=DEFAULT_SIZE,
+                        help='upper bound')
+arg_parser.add_argument('--size', type=int, default=DEFAULT_SIZE,
+                        help='array size')
+arg_parser.add_argument('--implementation', choices=['pure', 'cython'],
+                        default='cython', help='implementation to use')
 options = arg_parser.parse_args()
-data = array.array('d', list(range(size)))
+if options.implementation == 'cython':
+    from average import average, average_no_except
+elif options.implementation == 'pure':
+    from average_pure import average, average_no_except
+data = array.array('d', list(range(options.size)))
 print('with except:')
 try:
-    print(average.average(data, options.m, options.n))
+    print(average(data, options.m, options.n))
 except Exception as e:
     print('caught exception {0}: {1}'.format(str(e.__class__), str(e)))
 print('without except:')
 try:
-    print(average.average_no_except(data, options.m, options.n))
+    print(average_no_except(data, options.m, options.n))
     print('no exception caught')
 except Exception as e:
     print('caught exception {0}: {1}'.format(e.__class__, str(e)))

source-code/cython/Exceptions/setup.py

@@ -4,5 +4,6 @@
 from Cython.Build import cythonize
 
 setup(
-    ext_modules=cythonize('average.pyx')
+    ext_modules=cythonize(['average.pyx', 'average_pure.py'],
+                          language_level='3str')
 )

source-code/cython/Numpy/Convolution/Makefile

@@ -0,0 +1,20 @@
+VERSION = cpython-311-x86_64-linux-gnu
+CONVOLUTION_LIB = convolution_cython.$(VERSION).so
+CONVOLUTION_INDEXED_LIB = convolution_cython_indexed.$(VERSION).so
+CONVOLUTION_NO_CHECKS_LIB = convolution_cython_no_checks.$(VERSION).so
+
+all: $(CONVOLUTION_LIB) $(CONVOLUTION_INDEXED_LIB) $(CONVOLUTION_NO_CHECKS_LIB)
+
+$(CONVOLUTION_LIB): convolution.pyx
+	python setup.py build_ext --inplace
+
+$(CONVOLUTION_LIB): convolution_indexed.pyx
+	python setup.py build_ext --inplace
+
+$(CONVOLUTION_LIB): convolution_no_checks.pyx
+	python setup.py build_ext --inplace
+
+clean:
+	python setup.py clean
+	$(RM) $(wildcard *.c) $(wildcard *.so)
+	$(RM) -r build

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

@@ -0,0 +1,19 @@
+# Convolution
+
+Illustration of how to improve perfomance when using numpy arrays in Cython
+by doing indexing right.  The algorithm used is a naive implementation of
+convolution.
+
+
+## What is it?
+
+1. `convolution.py`: Python implementation of the algorithm (uses numpy).
+1. `convolution.pyx`: Cython implementation that adds types.
+1. `convolution_indexed.pyx`: Cython implementation that adds
+   numpy array element types and indexing information.
+1. `convolution_no_checks.pyx`: Cython implementation that adds
+   function decorator to avoid index checks for numpy arrays.
+1. `setup.py`: script to build the extensions.
+1. `Makefile`: make file to conveniently build the extensions.
+1. `driver.py`: Python scripts to run benchmark tests on the various
+   implementations.

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

@@ -0,0 +1,53 @@
+import numpy as np
+
+
+def convolve(input, kernel):
+    '''Naive convolution of matrices input and kernel
+
+    Parameters
+    ----------
+    input : numpy.ndarray
+        input matrix
+    kernel : numpy.ndarray
+        filter kernel matrix, dimensions must be odd
+
+    Returns
+    -------
+    output : numpy.ndarray
+        output matrix, it is not cropped
+
+    Raises
+    ------
+    ValueError
+        if dimensions of kernel are not odd
+    '''
+    if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
+        raise ValueError("Only odd dimensions on filter supported")
+    # s_mid and t_mid are number of pixels between the center pixel
+    # and the edge, ie for a 5x5 filter they will be 2.
+    #
+    # The output size is calculated by adding s_mid, t_mid to each
+    # side of the dimensions of the input image.
+    s_mid = kernel.shape[0] // 2
+    t_mid = kernel.shape[1] // 2
+    x_max = input.shape[0] + 2 * s_mid
+    y_max = input.shape[1] + 2 * t_mid
+    # Allocate result image.
+    output = np.zeros([x_max, y_max], dtype=input.dtype)
+    # Do convolution
+    for x in range(x_max):
+        for y in range(y_max):
+            # Calculate pixel value for h at (x,y). Sum one component
+            # for each pixel (s, t) of the filter kernel.
+            s_from = max(s_mid - x, -s_mid)
+            s_to = min((x_max - x) - s_mid, s_mid + 1)
+            t_from = max(t_mid - y, -t_mid)
+            t_to = min((y_max - y) - t_mid, t_mid + 1)
+            value = 0
+            for s in range(s_from, s_to):
+                for t in range(t_from, t_to):
+                    v = x - s_mid + s
+                    w = y - t_mid + t
+                    value += kernel[s_mid - s, t_mid - t]*input[v, w]
+            output[x, y] = value
+    return output

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

@@ -0,0 +1,61 @@
+import numpy as np
+cimport numpy as cnp
+cnp.import_array()
+
+DTYPE = np.float64
+ctypedef cnp.float64_t DTYPE_t
+
+def convolve(cnp.ndarray input, cnp.ndarray kernel):
+    '''Naive convolution of matrices input and kernel
+
+    Parameters
+    ----------
+    input : numpy.ndarray
+        input matrix
+    kernel : numpy.ndarray
+        filter kernel matrix, dimensions must be odd
+
+    Returns
+    -------
+    output : numpy.ndarray
+        output matrix, it is not cropped
+
+    Raises
+    ------
+    ValueError
+        if dimensions of kernel are not odd
+    '''
+    if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
+        raise ValueError("Only odd dimensions on filter supported")
+    assert input.dtype == DTYPE and kernel.dtype == DTYPE
+    # s_mid and t_mid are number of pixels between the center pixel
+    # and the edge, ie for a 5x5 filter they will be 2.
+    #
+    # The output size is calculated by adding s_mid, t_mid to each
+    # side of the dimensions of the input image.
+    cdef int s_mid = kernel.shape[0] // 2
+    cdef int t_mid = kernel.shape[1] // 2
+    cdef int x_max = input.shape[0] + 2 * s_mid
+    cdef int y_max = input.shape[1] + 2 * t_mid
+    # Allocate result image.
+    cdef cnp.ndarray output = np.zeros([x_max, y_max], dtype=input.dtype)
+    # Do convolution
+    cdef int x, y, s, t, v, w
+    cdef int s_from, s_to, t_from, t_to
+    cdef DTYPE_t value
+    for x in range(x_max):
+        for y in range(y_max):
+            # Calculate pixel value for h at (x,y). Sum one component
+            # for each pixel (s, t) of the filter kernel.
+            s_from = max(s_mid - x, -s_mid)
+            s_to = min((x_max - x) - s_mid, s_mid + 1)
+            t_from = max(t_mid - y, -t_mid)
+            t_to = min((y_max - y) - t_mid, t_mid + 1)
+            value = 0
+            for s in range(s_from, s_to):
+                for t in range(t_from, t_to):
+                    v = x - s_mid + s
+                    w = y - t_mid + t
+                    value += kernel[s_mid - s, t_mid - t]*input[v, w]
+            output[x, y] = value
+    return output

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

@@ -0,0 +1,61 @@
+import numpy as np
+cimport numpy as cnp
+cnp.import_array()
+
+DTYPE = np.float64
+ctypedef cnp.float64_t DTYPE_t
+
+def convolve(cnp.ndarray[DTYPE_t, ndim=2] input, cnp.ndarray[DTYPE_t, ndim=2] kernel):
+    '''Naive convolution of matrices input and kernel
+
+    Parameters
+    ----------
+    input : numpy.ndarray
+        input matrix
+    kernel : numpy.ndarray
+        filter kernel matrix, dimensions must be odd
+
+    Returns
+    -------
+    output : numpy.ndarray
+        output matrix, it is not cropped
+
+    Raises
+    ------
+    ValueError
+        if dimensions of kernel are not odd
+    '''
+    if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
+        raise ValueError("Only odd dimensions on filter supported")
+    assert input.dtype == DTYPE and kernel.dtype == DTYPE
+    # s_mid and t_mid are number of pixels between the center pixel
+    # and the edge, ie for a 5x5 filter they will be 2.
+    #
+    # The output size is calculated by adding s_mid, t_mid to each
+    # side of the dimensions of the input image.
+    cdef int s_mid = kernel.shape[0] // 2
+    cdef int t_mid = kernel.shape[1] // 2
+    cdef int x_max = input.shape[0] + 2 * s_mid
+    cdef int y_max = input.shape[1] + 2 * t_mid
+    # Allocate result image.
+    cdef cnp.ndarray[DTYPE_t, ndim=2] output = np.zeros([x_max, y_max], dtype=input.dtype)
+    # Do convolution
+    cdef int x, y, s, t, v, w
+    cdef int s_from, s_to, t_from, t_to
+    cdef DTYPE_t value
+    for x in range(x_max):
+        for y in range(y_max):
+            # Calculate pixel value for h at (x,y). Sum one component
+            # for each pixel (s, t) of the filter kernel.
+            s_from = max(s_mid - x, -s_mid)
+            s_to = min((x_max - x) - s_mid, s_mid + 1)
+            t_from = max(t_mid - y, -t_mid)
+            t_to = min((y_max - y) - t_mid, t_mid + 1)
+            value = 0
+            for s in range(s_from, s_to):
+                for t in range(t_from, t_to):
+                    v = x - s_mid + s
+                    w = y - t_mid + t
+                    value += kernel[s_mid - s, t_mid - t]*input[v, w]
+            output[x, y] = value
+    return output

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

@@ -0,0 +1,64 @@
+cimport cython
+import numpy as np
+cimport numpy as cnp
+cnp.import_array()
+
+DTYPE = np.float64
+ctypedef cnp.float64_t DTYPE_t
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def convolve(cnp.ndarray[DTYPE_t, ndim=2] input, cnp.ndarray[DTYPE_t, ndim=2] kernel):
+    '''Naive convolution of matrices input and kernel
+
+    Parameters
+    ----------
+    input : numpy.ndarray
+        input matrix
+    kernel : numpy.ndarray
+        filter kernel matrix, dimensions must be odd
+
+    Returns
+    -------
+    output : numpy.ndarray
+        output matrix, it is not cropped
+
+    Raises
+    ------
+    ValueError
+        if dimensions of kernel are not odd
+    '''
+    if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
+        raise ValueError("Only odd dimensions on filter supported")
+    assert input.dtype == DTYPE and kernel.dtype == DTYPE
+    # s_mid and t_mid are number of pixels between the center pixel
+    # and the edge, ie for a 5x5 filter they will be 2.
+    #
+    # The output size is calculated by adding s_mid, t_mid to each
+    # side of the dimensions of the input image.
+    cdef int s_mid = kernel.shape[0] // 2
+    cdef int t_mid = kernel.shape[1] // 2
+    cdef int x_max = input.shape[0] + 2 * s_mid
+    cdef int y_max = input.shape[1] + 2 * t_mid
+    # Allocate result image.
+    cdef cnp.ndarray[DTYPE_t, ndim=2] output = np.zeros([x_max, y_max], dtype=input.dtype)
+    # Do convolution
+    cdef int x, y, s, t, v, w
+    cdef int s_from, s_to, t_from, t_to
+    cdef DTYPE_t value
+    for x in range(x_max):
+        for y in range(y_max):
+            # Calculate pixel value for h at (x,y). Sum one component
+            # for each pixel (s, t) of the filter kernel.
+            s_from = max(s_mid - x, -s_mid)
+            s_to = min((x_max - x) - s_mid, s_mid + 1)
+            t_from = max(t_mid - y, -t_mid)
+            t_to = min((y_max - y) - t_mid, t_mid + 1)
+            value = 0
+            for s in range(s_from, s_to):
+                for t in range(t_from, t_to):
+                    v = x - s_mid + s
+                    w = y - t_mid + t
+                    value += kernel[s_mid - s, t_mid - t]*input[v, w]
+            output[x, y] = value
+    return output