On 07.02.2012 15:27, eat wrote:
> This is elegant and very fast as well!
Just be aware that it depends on C ordered input. So:
m,n = data.shape
cond = lamda x : (x >= t1) & (x <= t2)
x = cond(np.ascontiguousarray(data)).reshape((m//4, 4, n//4, 4))
found = np.any(np.any(x, axis=1), axis=2)
With Fortran ordered data, we must reshape in the opposite direction:
(m,n) ---> (4, m//4, 4, n//4)
That is:
m,n = data.shape
cond = lamda x : (x >= t1) & (x <= t2)
x = cond(np.asfortranarray(data)).reshape((4, m//4, 4, n//4))
found = np.any(np.any(x, axis=0), axis=1)
On the other hand, using as_strided instead of reshape should work with
any ordering. Think of as_strided as a generalization of reshape.
A problem with both these approaches is that it implicitely loops over
discontiguous memory. The only solution to that is to use C, Fortran or
Cython.
So in Cython:
import numpy as np
cimport numpy as np
cimport cython
cdef inline np.npy_bool cond(np.float64_t x,
np.float64_t t1, np.float64_t t2):
return 1 if (x >= t1 and x <= t2) else 0
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
def find(object data, np.float64_t t1, np.float64_t t2):
cdef np.ndarray[np.float64_t, ndim=2, mode='c'] x
cdef np.ndarray[np.npy_bool, ndim=2, mode='c'] found
cdef Py_ssize_t m, n, i, j
x = np.ascontiguousarray(data)
m,n = x.shape[0], x.shape[1]
found = np.zeros((m//4,n//4),dtype=bool)
for i in range(m):
for j in range(n):
found[i//4,j//4] = cond(x[i,j])
return found
It might be that the C compiler has difficulties optimizing the loop if
it think x and found cound be aliased, in which case the next logical
step would be C99 or Fortran...
Sturla
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
