Pure Python
You might want to test your code, it doesn't seem to do what you expect. Please run this script, test your code against mine and inspect the output:
import numpy as np
def find_first(a, index, value):
while index<a.size and a[index]!=value:
index += 1
return index
def find_end(a, index, value):
while index<a.size and a[index]==value:
index += 1
return index
def replace_run(a, begin, end, threshold, replace):
if end-begin+1 > threshold:
a[begin:end] = replace
def process_row(a, value, threshold, replace):
first = 0
while first < a.size:
if a[first]==value:
end = find_end(a, first, value)
replace_run(a, first, end, threshold, replace)
first = end
else:
first = find_first(a, first, value)
def replace_py(a, value, length, replace):
mat = a.copy()
for row in mat:
process_row(row, value, length, replace)
return mat
################################################################################
# Your code as posted in the question:
def replace_runs(a, search, run_length, replace = 2):
a_copy = a.copy() # Don't modify original
for i, row in enumerate(a):
runs = []
current_run = []
for j, val in enumerate(row):
if val == search:
current_run.append(j)
else:
if len(current_run) >= run_length or j == len(row) -1:
runs.append(current_run)
current_run = []
if len(current_run) >= run_length or j == len(row) -1:
runs.append(current_run)
for run in runs:
for col in run:
a_copy[i][col] = replace
return a_copy
# End of your code
################################################################################
def print_mismatch(a, b):
print 'Elementwise equals'
mat_equals = a==b
print mat_equals
print 'Reduced to rows'
for i, outcome in enumerate(np.logical_and.reduce(mat_equals, 1)):
print i, outcome
if __name__=='__main__':
np.random.seed(31)
shape = (20, 10)
mat = np.asarray(a=np.random.binomial(1, p=0.5, size=shape), dtype=np.int32)
mat.reshape(shape)
runs = replace_runs(mat, 1, 3, 2)
py = replace_py(mat, 1, 3, 2)
print 'Original'
print mat
print 'replace_runs()'
print runs
print 'replace_py()'
print py
print 'Mismatch between replace_runs() and replace_py()'
print_mismatch(runs, py)
Benchmarking doesn't make sense until your code is not fixed. So I will use mine, the replace_py() function for benchmarking.
The replace_py() implementation, which I think does what you want, is not pythonic, it has many anti-patterns. Nevertheless, it seems to be correct.
Timing:
np.random.seed(31)
shape = (100000, 10)
mat = np.asarray(a=np.random.binomial(1, p=0.5, size=shape), dtype=np.int32)
mat.reshape(shape)
%timeit replace_py(mat, 1, 3, 2)
1 loops, best of 3: 9.49 s per loop
Cython
I don't think your problem can be easily rewritten to use Numpy and vectorization. Perhaps a Numpy guru can do that but I am afraid the code will be either really obscure or slow (or both). To quote one of the Numpy developers:
[...] when it either requires a PhD in NumPy-ology to vectorize the
solution or it results in too much memory overhead, you can reach for
Cython [...]
So I re-wrote the replace_py() and the functions it calls in Cython, using typed memoryviews:
# cython: infer_types=True
# cython: boundscheck=False
# cython: wraparound=False
import numpy as np
cimport numpy as np
cdef inline int find_first(int[:] a, int index, int n, int value) nogil:
while index<n and a[index]!=value:
index += 1
return index
cdef inline int find_end(int[:] a, int index, int n, int value) nogil:
while index<n and a[index]==value:
index += 1
return index
cdef inline void replace_run(int[:] a, int begin, int end, int threshold, int replace) nogil:
if end-begin+1 > threshold:
for i in xrange(begin, end):
a[i] = replace
cdef inline void process_row(int[:] a, int value, int threshold, int replace) nogil:
cdef int first, end, n
first = 0
n = a.shape[0]
while first < n:
if a[first]==value:
end = find_end(a, first, n, value)
replace_run(a, first, end, threshold, replace)
first = end
else:
first = find_first(a, first, n, value)
def replace_cy(np.ndarray[np.int32_t, ndim=2] a, int value, int length, int replace):
cdef int[:, ::1] vmat
cdef int i, n
mat = a.copy()
vmat = mat
n = vmat.shape[0]
for i in xrange(n):
process_row(vmat[i], value, length, replace)
return mat
It required a little massaging and the code is more cluttered than the corresponding Python code given above. But it wasn't too much work and it was quite straightforward.
Timing:
np.random.seed(31)
shape = (100000, 10)
mat = np.asarray(a=np.random.binomial(1, p=0.5, size=shape), dtype=np.int32)
mat.reshape(shape)
%timeit replace_cy(mat, 1, 3, 2)
100 loops, best of 3: 8.16 ms per loop
That's a 1163x speed-up!
I have received help on Github and now the Numba version is also working; I just added @autojit to the pure Python code, except a[begin:end] = replace, see the discussion on Github where I got this workaround.
import numpy as np
from numba import autojit
@autojit
def find_first(a, index, value):
while index<a.size and a[index]!=value:
index += 1
return index
@autojit
def find_end(a, index, value):
while index<a.size and a[index]==value:
index += 1
return index
@autojit
def replace_run(a, begin, end, threshold, replace):
if end-begin+1 > threshold:
for i in xrange(begin, end):
a[i] = replace
@autojit
def process_row(a, value, threshold, replace):
first = 0
while first < a.size:
if a[first]==value:
end = find_end(a, first, value)
replace_run(a, first, end, threshold, replace)
first = end
else:
first = find_first(a, first, value)
@autojit
def replace_numba(a, value, length, replace):
mat = a.copy()
for row in mat:
process_row(row, value, length, replace)
return mat
Timing (with the usual input as above, code omitted):
1 loops, best of 3: 86.5 ms per loop
That's a 110x speed-up compared to the pure Python code for basically free!!! The Numba version is still a 10x slower than Cython, most likely due to not inlining the tiny functions, but I think it's amazing to get this speed-up basically for free, without messing up our Python code!
