Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-07 Thread Sturla Molden 
On 06.02.2012 22:27, Sturla Molden wrote:



 # Make a 4D view of this data, such that b[i,j]
 # is a 2D block with shape (4,4) (e.g. b[0,0] is
 # the same as a[:4, :4]).
 b = as_strided(a, shape=(a.shape[0]/4, a.shape[1]/4, 4, 4),
 strides=(4*a.strides[0], 4*a.strides[1], a.strides[0], 
 a.strides[1]))


 Yes :-) Being used to Fortran (and also MATLAB) this is the kind of mapping 
 It never occurs for me to think about. What else but NumPy is flexible enough 
 to do this? :-)

Actually, using as_strided is not needed. We can just reshape like this:

(m,n) --- (m//4, 4, n//4, 4)

and then use np.any along the two length-4 dimensions.

   m,n = data.shape
   cond = lamda x : (x = t1)  (x = t2)
   x = cond(data).reshape((m//4, 4, n//4, 4))
   found = np.any(np.any(x, axis=1), axis=2)


Sturla
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-07 Thread eat 
Hi

This is elegant and very fast as well!
On Tue, Feb 7, 2012 at 2:57 PM, Sturla Molden stu...@molden.no wrote:

 On 06.02.2012 22:27, Sturla Molden wrote:
 
 
 
  # Make a 4D view of this data, such that b[i,j]
  # is a 2D block with shape (4,4) (e.g. b[0,0] is
  # the same as a[:4, :4]).
  b = as_strided(a, shape=(a.shape[0]/4, a.shape[1]/4, 4, 4),
  strides=(4*a.strides[0], 4*a.strides[1], a.strides[0],
 a.strides[1]))
 
 
  Yes :-) Being used to Fortran (and also MATLAB) this is the kind of
 mapping It never occurs for me to think about. What else but NumPy is
 flexible enough to do this? :-)

 Actually, using as_strided is not needed. We can just reshape like this:

(m,n) --- (m//4, 4, n//4, 4)

 and then use np.any along the two length-4 dimensions.

   m,n = data.shape
   cond = lamda x : (x = t1)  (x = t2)

I guess you meant here cond= lambda x: (x= t1) (x= t2)

   x = cond(data).reshape((m//4, 4, n//4, 4))
   found = np.any(np.any(x, axis=1), axis=2)

Regards,
eat



 Sturla
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-07 Thread Sturla Molden 
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

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-07 Thread Sturla Molden 

   for i in range(m):
   for j in range(n):
   found[i//4,j//4] = cond(x[i,j])


Blah, that should be

 found[i//4,j//4] |= cond(x[i,j])


Sturla


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-06 Thread Sturla Molden 
Short answer: Create 16 view arrays, each with a stride of 4 in both 
dimensions. Test them against the conditions and combine the tests with an |= 
operator. Thus you replace the nested loop with one that has only 16 
iterations. Or reshape to 3 dimensions, the last with length 4, and you can do 
the same with only four view arrays.

Sturla

Sendt fra min iPad

Den 6. feb. 2012 kl. 20:16 skrev Moroney, Catherine M (388D) 
catherine.m.moro...@jpl.nasa.gov:

 Hello,
 
 I have to write a code to downsample an array in a specific way, and I am 
 hoping that
 somebody can tell me how to do this without the nested do-loops.  Here is the 
 problem
 statement:  Segment a (MXN) array into 4x4 squares and set a flag if any of 
 the pixels
 in that 4x4 square meet a certain condition.
 
 Here is the code that I want to rewrite avoiding loops:
 
 shape_out = (data_in.shape[0]/4, data_in.shape[1]/4)
 found = numpy.zeros(shape_out).astype(numpy.bool)
 
 for i in xrange(0, shape_out[0]):
for j in xrange(0, shape_out[1]):
 
excerpt = data_in[i*4:(i+1)*4, j*4:(j+1)*4]
mask = numpy.where( (excerpt = t1)  (excerpt = t2), True, False)
if (numpy.any(mask)):
found[i,j] = True
 
 Thank you for any hints and education!
 
 Catherine
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-06 Thread Warren Weckesser 
On Mon, Feb 6, 2012 at 2:57 PM, Sturla Molden stu...@molden.no wrote:

 Short answer: Create 16 view arrays, each with a stride of 4 in both
 dimensions. Test them against the conditions and combine the tests with an
 |= operator. Thus you replace the nested loop with one that has only 16
 iterations. Or reshape to 3 dimensions, the last with length 4, and you can
 do the same with only four view arrays.

 Sturla



Or use 'as_strided' from the stride_tricks module:

import numpy as np
from numpy.lib.stride_tricks import as_strided

# Make some demonstration data.
a = np.random.random_integers(0, 99, size=(12,16))

# Find 4x4 blocks whose values are between these limits.
low = 5
high = 94

# Make a 4D view of this data, such that b[i,j]
# is a 2D block with shape (4,4) (e.g. b[0,0] is
# the same as a[:4, :4]).
b = as_strided(a, shape=(a.shape[0]/4, a.shape[1]/4, 4, 4),
   strides=(4*a.strides[0], 4*a.strides[1], a.strides[0],
a.strides[1]))

# Reduce this to a 2D array of boolean values.  The value is
# True if the corresponding 4x4 block contains only values
# between low and high (inclusive).
result = np.all(np.all((b = low)  (b = high), axis=-1), axis=-1)

print a
print result


Output (from ipython)t:

In [111]: run 2Dto4by4blocks.py
[[ 5 50 62 43 52 21 67 70 55 12 25 21  1 95 73 31]
 [44  4 60 27 93 54 25 87 15 22 37 31 46 84 10 46]
 [35 26 11 76 91 79 58 92 57 62 27 27 17 19 39 86]
 [94 96 21 36 90 18 80 62 91 68 39 22 68 31 71 48]
 [89 34 52 80 93 73 54 13 25 28 57 32 55 42  3 13]
 [65 68 41 82 55 81 64 59 73  4 44 46 91  1 86 52]
 [99 87 21 70 26 64  2 10 62 82 52 67 85 88 45 53]
 [33 10  6  3 46 71 17 58 20 56 30 18 19 17 60 76]
 [18 22 62 53 45 21 83 86 69 35 32 36 33 74 81 70]
 [24 39 93 12 37  4  4 16 45 59 46  4 90 24  1 13]
 [26 37 11 11 24 58  6 44 43 44 94 55 22  8  7 85]
 [26 91 31 75 72 25 23 89  3 30 45 93 62 72 96 39]]
[[False  True  True False]
 [False False False False]
 [ True False False False]]


Warren



 Sendt fra min iPad

 Den 6. feb. 2012 kl. 20:16 skrev Moroney, Catherine M (388D) 
 catherine.m.moro...@jpl.nasa.gov:

  Hello,
 
  I have to write a code to downsample an array in a specific way, and I
 am hoping that
  somebody can tell me how to do this without the nested do-loops.  Here
 is the problem
  statement:  Segment a (MXN) array into 4x4 squares and set a flag if any
 of the pixels
  in that 4x4 square meet a certain condition.
 
  Here is the code that I want to rewrite avoiding loops:
 
  shape_out = (data_in.shape[0]/4, data_in.shape[1]/4)
  found = numpy.zeros(shape_out).astype(numpy.bool)
 
  for i in xrange(0, shape_out[0]):
 for j in xrange(0, shape_out[1]):
 
 excerpt = data_in[i*4:(i+1)*4, j*4:(j+1)*4]
 mask = numpy.where( (excerpt = t1)  (excerpt = t2), True,
 False)
 if (numpy.any(mask)):
 found[i,j] = True
 
  Thank you for any hints and education!
 
  Catherine
  ___
  NumPy-Discussion mailing list
  NumPy-Discussion@scipy.org
  http://mail.scipy.org/mailman/listinfo/numpy-discussion
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-06 Thread Sturla Molden 

Something like this:

m,n = data.shape
x = data.reshape((m,n//4,4))
z  = (x[0::4,...] = t1)  (x[0::4,...] = t1)
z |= (x[1::4,...] = t1)  (x[1::4,...] = t1)
z |= (x[2::4,...] = t1)  (x[2::4,...] = t1)
z |= (x[3::4,...] = t1)  (x[3::4,...] = t1)
found = np.any(z, axis=2)

Sturla

Sendt fra min iPad

Den 6. feb. 2012 kl. 21:57 skrev Sturla Molden stu...@molden.no:

 Short answer: Create 16 view arrays, each with a stride of 4 in both 
 dimensions. Test them against the conditions and combine the tests with an |= 
 operator. Thus you replace the nested loop with one that has only 16 
 iterations. Or reshape to 3 dimensions, the last with length 4, and you can 
 do the same with only four view arrays.
 
 Sturla
 
 Sendt fra min iPad
 
 Den 6. feb. 2012 kl. 20:16 skrev Moroney, Catherine M (388D) 
 catherine.m.moro...@jpl.nasa.gov:
 
 Hello,
 
 I have to write a code to downsample an array in a specific way, and I am 
 hoping that
 somebody can tell me how to do this without the nested do-loops.  Here is 
 the problem
 statement:  Segment a (MXN) array into 4x4 squares and set a flag if any of 
 the pixels
 in that 4x4 square meet a certain condition.
 
 Here is the code that I want to rewrite avoiding loops:
 
 shape_out = (data_in.shape[0]/4, data_in.shape[1]/4)
 found = numpy.zeros(shape_out).astype(numpy.bool)
 
 for i in xrange(0, shape_out[0]):
   for j in xrange(0, shape_out[1]):
 
   excerpt = data_in[i*4:(i+1)*4, j*4:(j+1)*4]
   mask = numpy.where( (excerpt = t1)  (excerpt = t2), True, False)
   if (numpy.any(mask)):
   found[i,j] = True
 
 Thank you for any hints and education!
 
 Catherine
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-06 Thread Sturla Molden 
The last t1 on each lineis of course t2. Sorry for the typo. Hard to code on an 
ipad ;-)

Sturla

Sendt fra min iPad

Den 6. feb. 2012 kl. 22:12 skrev Sturla Molden stu...@molden.no:

 
 Something like this:
 
 m,n = data.shape
 x = data.reshape((m,n//4,4))
 z  = (x[0::4,...] = t1)  (x[0::4,...] = t1)
 z |= (x[1::4,...] = t1)  (x[1::4,...] = t1)
 z |= (x[2::4,...] = t1)  (x[2::4,...] = t1)
 z |= (x[3::4,...] = t1)  (x[3::4,...] = t1)
 found = np.any(z, axis=2)
 
 Sturla
 
 Sendt fra min iPad
 
 Den 6. feb. 2012 kl. 21:57 skrev Sturla Molden stu...@molden.no:
 
 Short answer: Create 16 view arrays, each with a stride of 4 in both 
 dimensions. Test them against the conditions and combine the tests with an 
 |= operator. Thus you replace the nested loop with one that has only 16 
 iterations. Or reshape to 3 dimensions, the last with length 4, and you can 
 do the same with only four view arrays.
 
 Sturla
 
 Sendt fra min iPad
 
 Den 6. feb. 2012 kl. 20:16 skrev Moroney, Catherine M (388D) 
 catherine.m.moro...@jpl.nasa.gov:
 
 Hello,
 
 I have to write a code to downsample an array in a specific way, and I am 
 hoping that
 somebody can tell me how to do this without the nested do-loops.  Here is 
 the problem
 statement:  Segment a (MXN) array into 4x4 squares and set a flag if any of 
 the pixels
 in that 4x4 square meet a certain condition.
 
 Here is the code that I want to rewrite avoiding loops:
 
 shape_out = (data_in.shape[0]/4, data_in.shape[1]/4)
 found = numpy.zeros(shape_out).astype(numpy.bool)
 
 for i in xrange(0, shape_out[0]):
  for j in xrange(0, shape_out[1]):
 
  excerpt = data_in[i*4:(i+1)*4, j*4:(j+1)*4]
  mask = numpy.where( (excerpt = t1)  (excerpt = t2), True, False)
  if (numpy.any(mask)):
  found[i,j] = True
 
 Thank you for any hints and education!
 
 Catherine
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-06 Thread Sturla Molden 


 
 # Make a 4D view of this data, such that b[i,j]
 # is a 2D block with shape (4,4) (e.g. b[0,0] is
 # the same as a[:4, :4]).
 b = as_strided(a, shape=(a.shape[0]/4, a.shape[1]/4, 4, 4),
strides=(4*a.strides[0], 4*a.strides[1], a.strides[0], 
 a.strides[1]))
 

Yes :-) Being used to Fortran (and also MATLAB) this is the kind of mapping It 
never occurs for me to think about. What else but NumPy is flexible enough to 
do this? :-)

Sturla
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-06 Thread eat 
Hi,

On Mon, Feb 6, 2012 at 9:16 PM, Moroney, Catherine M (388D) 
catherine.m.moro...@jpl.nasa.gov wrote:

 Hello,

 I have to write a code to downsample an array in a specific way, and I am
 hoping that
 somebody can tell me how to do this without the nested do-loops.  Here is
 the problem
 statement:  Segment a (MXN) array into 4x4 squares and set a flag if any
 of the pixels
 in that 4x4 square meet a certain condition.

 Here is the code that I want to rewrite avoiding loops:

 shape_out = (data_in.shape[0]/4, data_in.shape[1]/4)
 found = numpy.zeros(shape_out).astype(numpy.bool)

 for i in xrange(0, shape_out[0]):
for j in xrange(0, shape_out[1]):

excerpt = data_in[i*4:(i+1)*4, j*4:(j+1)*4]
mask = numpy.where( (excerpt = t1)  (excerpt = t2),
 True, False)
if (numpy.any(mask)):
found[i,j] = True

 Thank you for any hints and education!

 Catherine

Following Warrens answer a slight demonstration of code like this:

import numpy as np

def ds_0(data_in, t1= 1, t2= 4):

shape_out= (data_in.shape[0]/ 4, data_in.shape[1]/ 4)

found= np.zeros(shape_out).astype(np.bool)

for i in xrange(0, shape_out[0]):

for j in xrange(0, shape_out[1]):

excerpt= data_in[i* 4: (i+ 1)* 4, j* 4: (j+ 1)* 4]

mask= np.where((excerpt= t1) (excerpt= t2), True, False)

if (np.any(mask)):

found[i, j]= True

return found

# with stride_tricks you may cook up something like this:

from numpy.lib.stride_tricks import as_strided as ast

def _ss(dt, ds, s):

return {'shape': (ds[0]/ s[0], ds[1]/ s[1])+ s,

'strides': (s[0]* dt[0], s[1]* dt[1])+ dt}

def _view(D, shape= (4, 4)):

return ast(D, **_ss(D.strides, D.shape, shape))

def ds_1(data_in, t1= 1, t2= 4):

# return _view(data_in)

excerpt= _view(data_in)

mask= np.where((excerpt= t1) (excerpt= t2), True, False)

return mask.sum(2).sum(2).astype(np.bool)

 if __name__ == '__main__':

from numpy.random import randint

r= randint(777, size= (64, 288)); print r

print np.allclose(ds_0(r), ds_1(r))




My 2 cents,
eat

 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] avoiding loops when downsampling arrays

2012-02-06 Thread eat 
Hi,

Sorry for my latest post, hands way too quick ;(

On Mon, Feb 6, 2012 at 9:16 PM, Moroney, Catherine M (388D) 
catherine.m.moro...@jpl.nasa.gov wrote:

 Hello,

 I have to write a code to downsample an array in a specific way, and I am
 hoping that
 somebody can tell me how to do this without the nested do-loops.  Here is
 the problem
 statement:  Segment a (MXN) array into 4x4 squares and set a flag if any
 of the pixels
 in that 4x4 square meet a certain condition.

 Here is the code that I want to rewrite avoiding loops:

 shape_out = (data_in.shape[0]/4, data_in.shape[1]/4)
 found = numpy.zeros(shape_out).astype(numpy.bool)

 for i in xrange(0, shape_out[0]):
for j in xrange(0, shape_out[1]):

excerpt = data_in[i*4:(i+1)*4, j*4:(j+1)*4]
mask = numpy.where( (excerpt = t1)  (excerpt = t2),
 True, False)
if (numpy.any(mask)):
found[i,j] = True

 Thank you for any hints and education!

Following closely with Warrens answer a slight demonstration of code like
this:
import numpy as np

 def ds_0(data_in, t1= 1, t2= 4):

shape_out= (data_in.shape[0]/ 4, data_in.shape[1]/ 4)

found= np.zeros(shape_out).astype(np.bool)

for i in xrange(0, shape_out[0]):

for j in xrange(0, shape_out[1]):

excerpt= data_in[i* 4: (i+ 1)* 4, j* 4: (j+ 1)* 4]

mask= np.where((excerpt= t1) (excerpt= t2), True, False)

if (np.any(mask)):

found[i, j]= True

return found


 # with stride_tricks you may cook up something like this:

from numpy.lib.stride_tricks import as_strided as ast


 def _ss(dt, ds, s):

return {'shape': (ds[0]/ s[0], ds[1]/ s[1])+ s,

'strides': (s[0]* dt[0], s[1]* dt[1])+ dt}


 def _view(D, shape= (4, 4)):

return ast(D, **_ss(D.strides, D.shape, shape))


 def ds_1(data_in, t1= 1, t2= 4):

excerpt= _view(data_in)

mask= np.where((excerpt= t1) (excerpt= t2), True, False)

return mask.sum(2).sum(2).astype(np.bool)


 if __name__ == '__main__':

from numpy.random import randint

r= randint(777, size= (64, 288)); print r

print np.allclose(ds_0(r), ds_1(r))


and when run, it will yield like:

In []: run dsa

[[ 60 470 521 ..., 147 435 295]

 [246 127 662 ..., 718 525 256]

 [354 384 205 ..., 225 364 239]

 ...,

 [277 428 201 ..., 460 282 433]

 [ 27 407 130 ..., 245 346 309]

 [649 157 153 ..., 316 613 570]]

True

and compared in performance wise:
In []: %timeit ds_0(r)
10 loops, best of 3: 56.3 ms per loop

In []: %timeit ds_1(r)
100 loops, best of 3: 2.17 ms per loop


My 2 cents,

eat



 Catherine
 ___
 NumPy-Discussion mailing list
 NumPy-Discussion@scipy.org
 http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion