Den 29.03.2011 16:49, skrev Sturla Molden:
>
> This will not work. A boolean array is not compactly stored, but an
> array of bytes. Only the first bit 0 is 1 with probability p, bits 1 to
> 7 bits are 1 with probability 0. We thus have to do this 8 times for
> each byte, shift left by range(8), and combine them with binary or.
> Also the main use of random bits is crypto, which requires the use of
> /dev/urandom or CrypGenRandom instead of Mersenne Twister (np.random.rand).
Here's a cleaned-up one for those who might be interested :-)
Sturla
import numpy as np
import os
def randombits(n, p, intention='numerical'):
"""
Returns an array with packed bits drawn from n Bernoulli
trials with successrate p.
"""
assert (intention in ('numerical','crypto'))
# number of bytes
m = (n >> 3) + 1 if n % 8 else n >> 3
if intention == 'numerical':
# Mersenne Twister
rflt = np.random.rand(m*8)
else:
# /dev/urandom on Linux, Apple, et al.,
# CryptGenRandom on Windows
rflt = np.frombuffer(os.urandom(m*64),dtype=np.uint64)
rflt = rflt / float(2**64)
b = (rflt.reshape((m,8))<p)
# pack the bits
b = (b<<range(8)).sum(axis=1).astype(np.uint8)
# zero the trailing m*8 - n bits
b[-1] &= (0xFF >> (m*8 - n))
return b
def bitcount(a, bytewise=False):
"""
Count the number of set bits in an array of np.uint8.
"""
assert(a.dtype == np.uint8)
c = a[:,np.newaxis].repeat(8, axis=1) >> range(8) & 0x01
return (c.sum(axis=1) if bytewise else c.sum())
if __name__ == '__main__':
b = randombits(int(1e6), .1, intent='numerical')
print bitcount(b) # should be close to 100000
b = randombits(int(1e6), .1, intent='crypto')
print bitcount(b) # should be close to 100000
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
