[Numpy-discussion] random number genration

[Sturla Molden]

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