Commit c10d130cbb9 by Pierre Augier

skimage benchmarks: first conclusions

parent 25415654a4e
......@@ -11,7 +11,6 @@ See the `ROADMAP.rst file
- Cython backend: less bugs, better support for fused types, nonecheck,
cdivision, ...
- Fix default parameters for Pythran
0.4.0 (2019-09-22)
......
CLANG=clang-6.0
CXXLANG=clang++-6.0
all:
python setup_pyx.py
python setup.py
all_clang:
CXX=$(CXXLANG) CC=$(CLANG) python setup_pyx.py
CXX=$(CXXLANG) CC=$(CLANG) python setup.py
clean_all:
rm -f pyx/*.so pyx/*.html pyx/*.c
rm -rf future/__pythran__ future/__cython__ future/__numba__
\ No newline at end of file
rm -rf future/__pythran__ future/__cython__ future/__numba__
......@@ -36,18 +36,30 @@ python setup.py`.
## What we need now
Pieces of code to create realistic input parameters for all boosted functions (for benchmarking):
- _colormixer.py (add, multiply, brightness, sigmoid_gamma, gamma,
py_hsv_2_rgb, py_rgb_2_hsv, hsv_add, hsv_multiply)
- _convex_hull.py (possible_hull)
- _greyreconstruct.py (reconstruction_loop)
- _hessian_det_appx.py (_hessian_matrix_det)
- _moments_cy.py (moments_hu)
- _radon_transform.py (sart_projection_update)
- _unwrap_1d.py (unwrap_1d)
- brief_cy.py (_brief_loop)
- cmorph.py (_dilate, _erode)
Pieces of code to create realistic input parameters for all boosted functions
(for benchmarking, see the directory `setup_codes`):
- [ ] _colormixer.py
- [ ] add
- [ ] multiply
- [ ] brightness
- [ ] sigmoid_gamma
- [ ] gamma
- [ ] py_hsv_2_rgb
- [ ] py_rgb_2_hsv
- [ ] hsv_add
- [ ] hsv_multiply
- [ ] _convex_hull.py (possible_hull)
- [x] _greyreconstruct.py (reconstruction_loop)
- [ ] _hessian_det_appx.py (_hessian_matrix_det)
- [ ] _moments_cy.py (moments_hu)
- [ ] _radon_transform.py (sart_projection_update)
- [ ] _unwrap_1d.py (unwrap_1d)
- [ ] brief_cy.py (_brief_loop)
- [ ] cmorph.py
- [x] _dilate
- [ ] _erode
## Next tasks
......@@ -67,11 +79,73 @@ Pieces of code to create realistic input parameters for all boosted functions (f
And more complicated things!
- [ ] `with nogil:` (for example for _convex_hull.py)
- [ ] `from libc.math cimport exp, pow` (how?)
- [ ] casting: `<np.uint8_t> op_result`, which in Cython is different from
`np.uint8(op_result)` :-(
`np.uint8(op_result)` :-(
- [ ] `cdef void foo(int[:] a) nogil:` (bugs Cython)
- [ ] `with nogil:`
A question: do we have to support C array creation (`cdef cnp.uint8_t
lut[256]` or `cdef float HSV[3]`)? Is it *really* more efficient than standard Numpy array?
Questions:
- Do we have to support C array creation (`cdef cnp.uint8_t lut[256]` or
`cdef float HSV[3]`)? Is it *really* more efficient than standard Numpy
array?
- Would it be ok to have all function defined with `cpdef` ?
## About the Transonic code
- [ ] _colormixer.py
It is a difficult case, with advanced Cython code (`with nogil:`, `cdef`, C
arrays, `<cnp.int16_t>`, `from libc.math cimport exp, pow`).
- [ ] add
- [ ] multiply
- [ ] brightness
- [ ] sigmoid_gamma
- [ ] gamma
- [ ] py_hsv_2_rgb
- [ ] py_rgb_2_hsv
- [ ] hsv_add
- [ ] hsv_multiply
- [x] _convex_hull.py (possible_hull)
Done (`cython -a` white) except one (important) `with nogil:`.
- [x] _greyreconstruct.py (reconstruction_loop)
Done (`cython -a` white) except one (important) `with nogil:`. The 3
Transonic backends produce codes faster than the old .pyx (for Cython, it
should be related to the missing `with nogil:`).
- [ ] _hessian_det_appx.py (_hessian_matrix_det)
Not an easy one (`cdef ... nogil`, `with nogil`).
- [x] _moments_cy.py (moments_hu)
This one is easy!
- [ ] _radon_transform.py (sart_projection_update)
Needs cast, `with nogil` and `from libc.math cimport ...`.
- [ ] _unwrap_1d.py (unwrap_1d)
Needs cast and `from libc.math cimport ...`.
- [ ] brief_cy.py (_brief_loop)
Needs `with nogil`.
- [x] cmorph.py
- [x] _dilate
- [x] _erode
`from libc.stdlib cimport malloc, free` is used in the .pyx file! Thus the
Transonic code and the generated Cython code is simpler (using `np.empty`).
The generated Cython seems to be ~10% slower but the generated Pythran is
faster than the .pyx from skimage.
from bench_util import bench_one, statements
for mod_name, func_name in statements.keys():
bench_one(mod_name)
import argparse
from bench_util import bench_one
parser = argparse.ArgumentParser(description="Run one benchmark")
parser.add_argument("module", help="Module name")
args = parser.parse_args()
bench_one(args.module)
from subprocess import getoutput
from pathlib import Path
from transonic.util import timeit
statements = {
("cmorph", "_dilate"): "_dilate(image, selem, out, shift_x, shift_y)",
(
......@@ -9,54 +11,71 @@ statements = {
): "reconstruction_loop(ranks, prev, next_, strides, current_idx, image_stride)",
}
name_module = "cmorph"
name_function = "_dilate"
import_from_skimage = {
(
"_greyreconstruct",
"reconstruction_loop",
): "from skimage.morphology._greyreconstruct import reconstruction_loop"
}
name_module = "_greyreconstruct"
name_function = "reconstruction_loop"
print(f"module: {name_module}")
stmt = statements[(name_module, name_function)]
print(stmt)
def bench_one(name_module="cmorph", func=None):
path_setup = Path("setup_codes") / f"{name_module}_{name_function}.py"
if func is not None:
raise NotImplementedError
if not path_setup.exists():
raise RuntimeError
functions = [
(mod, func_) for (mod, func_) in statements.keys() if mod == name_module
]
with open(path_setup) as file:
setup = file.read()
if not functions:
raise ValueError(f"bad name_module: {name_module}")
setup_pyx = setup.replace(
f"from future.{name_module} import", f"from pyx.{name_module} import"
)
name_function = functions[0][1]
code_pyx = f"""
from transonic.util import timeit
setup = '''{setup_pyx}'''
stmt = '''{stmt}'''
print(timeit(stmt, setup, total_duration=2))
"""
print(f"module: {name_module}")
stmt = statements[(name_module, name_function)]
print(stmt)
path_setup = Path("setup_codes") / f"{name_module}_{name_function}.py"
if not path_setup.exists():
raise RuntimeError
with open(path_setup) as file:
setup = file.read()
code = f"""
setup_pyx = setup.replace(
f"from future.{name_module} import", f"from pyx.{name_module} import"
)
code = f"""
from transonic.util import timeit
setup = '''{setup}'''
stmt = '''{stmt}'''
print(timeit(stmt, setup, total_duration=2))
"""
"""
time_old = timeit(stmt, setup_pyx, total_duration=2)
print(f'cython pyx skimage {time_old:.2e} s (= norm)')
with open("tmp.py", "w") as file:
file.write(code_pyx)
with open("tmp.py", "w") as file:
file.write(code)
time_old = float(getoutput("python tmp.py"))
for backend in ("cython", "pythran", "numba"):
time = float(getoutput(f"TRANSONIC_BACKEND='{backend}' python tmp.py"))
print(f"{backend:18s} {time:.2e} s (= {time/time_old:.4f} * norm)")
print(f'cython "skimage" {time_old:.2e} s (= norm)')
# print(getoutput("TRANSONIC_NO_REPLACE=1 python tmp.py"))
with open("tmp.py", "w") as file:
file.write(code.format(package="future"))
if (name_module, name_function) not in import_from_skimage:
return
for backend in ("cython", "pythran", "numba"):
time = float(getoutput(f"TRANSONIC_BACKEND='{backend}' python tmp.py"))
print(f"{backend:16s} {time:.2e} s (= {time/time_old:.4f} * norm)")
setup_from_skimage = setup.replace(
f"from future.{name_module} import {name_function}",
import_from_skimage[(name_module, name_function)],
)
time = timeit(stmt, setup_from_skimage, total_duration=2)
# print(getoutput("TRANSONIC_NO_REPLACE=1 python tmp.py"))
print(f"{'from skimage':18s} {time:.2e} s (= {time/time_old:.4f} * norm)")
import numpy as np
from transonic import boost
from transonic import boost, Array
@boost
@boost(wraparound=False, boundscheck=False, cdivision=True, nonecheck=False)
def possible_hull(img: "uint8[:,:]"):
"""Return positions of pixels that possibly belong to the convex hull.
......@@ -19,8 +19,10 @@ def possible_hull(img: "uint8[:,:]"):
the convex hull.
"""
rows = img.shape[0]
cols = img.shape[1]
r: np.intp
c: np.intp
rows: np.intp = img.shape[0]
cols: np.intp = img.shape[1]
# Output: rows storage slots for left boundary pixels
# cols storage slots for top boundary pixels
......@@ -29,8 +31,10 @@ def possible_hull(img: "uint8[:,:]"):
coords = np.ones((2 * (rows + cols), 2), dtype=np.intp)
coords *= -1
rows_cols = rows + cols
rows_2_cols = 2 * rows + cols
nonzero: Array[np.intp, "2d", "C", "memview"] = coords
rows_cols: np.intp = rows + cols
rows_2_cols: np.intp = 2 * rows + cols
for r in range(rows):
......@@ -44,20 +48,23 @@ def possible_hull(img: "uint8[:,:]"):
rows_2_cols_c = rows_2_cols + c
# Left check
if coords[r, 1] == -1:
coords[r, 0] = r
coords[r, 1] = c
if nonzero[r, 1] == -1:
nonzero[r, 0] = r
nonzero[r, 1] = c
# Right check
elif coords[rows_cols_r, 1] < c:
coords[rows_cols_r] = r, c
elif nonzero[rows_cols_r, 1] < c:
nonzero[rows_cols_r, 0] = r
nonzero[rows_cols_r, 1] = c
# Top check
if coords[rows_c, 1] == -1:
coords[rows_c] = r, c
if nonzero[rows_c, 1] == -1:
nonzero[rows_c, 0] = r
nonzero[rows_c, 1] = c
# Bottom check
elif coords[rows_2_cols_c, 0] < r:
coords[rows_2_cols_c] = r, c
elif nonzero[rows_2_cols_c, 0] < r:
nonzero[rows_2_cols_c, 0] = r
nonzero[rows_2_cols_c, 1] = c
return coords[coords[:, 0] != -1]
import numpy as np
from transonic import boost
from transonic import boost, Array
@boost
@boost(wraparound=False, cdivision=True, nonecheck=False)
def moments_hu(nu: "float64[:,:]"):
hu = np.zeros((7, ), dtype=np.double)
t0 = nu[3, 0] + nu[1, 2]
t1 = nu[2, 1] + nu[0, 3]
q0 = t0 * t0
q1 = t1 * t1
n4 = 4 * nu[1, 1]
s = nu[2, 0] + nu[0, 2]
d = nu[2, 0] - nu[0, 2]
hu: Array[np.float64, "1d", "C"] = np.zeros((7,), dtype=np.float64)
t0: np.float64 = nu[3, 0] + nu[1, 2]
t1: np.float64 = nu[2, 1] + nu[0, 3]
q0: np.float64 = t0 * t0
q1: np.float64 = t1 * t1
n4: np.float64 = 4 * nu[1, 1]
s: np.float64 = nu[2, 0] + nu[0, 2]
d: np.float64 = nu[2, 0] - nu[0, 2]
hu[0] = s
hu[1] = d * d + n4 * nu[1, 1]
hu[3] = q0 + q1
......
......@@ -4,6 +4,8 @@ from transonic import boost, Optional, Array
A = Array[np.uint8, "2d", "memview"]
A1d = Array[np.intp, "1d", "memview"]
@boost(wraparound=False, boundscheck=False, cdivision=True, nonecheck=False)
def _dilate(
......@@ -44,16 +46,16 @@ def _dilate(
srows: np.intp = selem.shape[0]
scols: np.intp = selem.shape[1]
centre_r: np.intp = int(selem.shape[0] / 2) - shift_y
centre_c: np.intp = int(selem.shape[1] / 2) - shift_x
centre_r: np.intp = int(srows / 2) - shift_y
centre_c: np.intp = int(scols / 2) - shift_x
image = np.ascontiguousarray(image)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint8)
selem_num: np.intp = np.sum(np.asarray(selem) != 0)
sr: Array[np.intp, "1d", "memview"] = np.empty(selem_num, dtype=np.intp)
sc: Array[np.intp, "1d", "memview"] = np.empty(selem_num, dtype=np.intp)
sr: A1d = np.empty(selem_num, dtype=np.intp)
sc: A1d = np.empty(selem_num, dtype=np.intp)
s: np.intp = 0
r: np.intp
......@@ -87,7 +89,7 @@ def _dilate(
return np.asarray(out)
@boost
@boost(wraparound=False, boundscheck=False, cdivision=True, nonecheck=False)
def _erode(
image: A,
selem: A,
......@@ -120,23 +122,26 @@ def _erode(
The result of the morphological erosion.
"""
rows = image.shape[0]
cols = image.shape[1]
srows = selem.shape[0]
scols = selem.shape[1]
rows: np.intp = image.shape[0]
cols: np.intp = image.shape[1]
srows: np.intp = selem.shape[0]
scols: np.intp = selem.shape[1]
centre_r = int(selem.shape[0] / 2) - shift_y
centre_c = int(selem.shape[1] / 2) - shift_x
centre_r: np.intp = int(srows / 2) - shift_y
centre_c: np.intp = int(scols / 2) - shift_x
image = np.ascontiguousarray(image)
if out is None:
out = np.zeros((rows, cols), dtype=np.uint8)
selem_num = np.sum(np.asarray(selem) != 0)
sr = np.empty(selem_num, dtype=np.intp)
sc = np.empty(selem_num, dtype=np.intp)
selem_num: np.intp = np.sum(np.asarray(selem) != 0)
sr: A1d = np.empty(selem_num, dtype=np.intp)
sc: A1d = np.empty(selem_num, dtype=np.intp)
s: np.intp = 0
r: np.intp
c: np.intp
s = 0
for r in range(srows):
for c in range(scols):
if selem[r, c] != 0:
......@@ -144,6 +149,11 @@ def _erode(
sc[s] = c - centre_c
s += 1
local_max: np.uint8
value: np.uint8
rr: np.intp
cc: np.intp
for r in range(rows):
for c in range(cols):
local_min = 255
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment