On Tue, Mar 29, 2011 at 11:59 AM, Sturla Molden <stu...@molden.no> wrote:
> 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 :-) > How about adding it to http://www.scipy.org/Cookbook? Warren > > 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 >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion