<br><br><div class="gmail_quote">On Tue, Mar 29, 2011 at 11:59 AM, Sturla Molden <span dir="ltr"><<a href="mailto:sturla@molden.no">sturla@molden.no</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">

Den 29.03.2011 16:49, skrev Sturla Molden:<br>
><br>
> This will not work. A boolean array is not compactly stored, but an<br>
> array of bytes. Only the first bit 0 is 1 with probability p, bits 1 to<br>
> 7 bits are 1 with probability 0. We thus have to do this 8 times for<br>
> each byte, shift left by range(8), and combine them with binary or.<br>
> Also the main use of random bits is crypto, which requires the use of<br>
> /dev/urandom or CrypGenRandom instead of Mersenne Twister (np.random.rand).<br>
<br>
<br>
Here's a cleaned-up one for those who might be interested :-)<br></blockquote><div><br><br><br>How about adding it to <a href="http://www.scipy.org/Cookbook">http://www.scipy.org/Cookbook</a>?<br><br>Warren<br><br> <br>

</div><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<br>
Sturla<br>
<br>
<br>
<br>
import numpy as np<br>
import os<br>
<br>
<br>
def randombits(n, p, intention='numerical'):<br>
     """<br>
     Returns an array with packed bits drawn from n Bernoulli<br>
     trials with successrate p.<br>
     """<br>
     assert (intention in ('numerical','crypto'))<br>
     # number of bytes<br>
     m = (n >> 3) + 1 if n % 8 else n >> 3<br>
     if intention == 'numerical':<br>
         # Mersenne Twister<br>
         rflt = np.random.rand(m*8)<br>
     else:<br>
         # /dev/urandom on Linux, Apple, et al.,<br>
         # CryptGenRandom on Windows<br>
         rflt = np.frombuffer(os.urandom(m*64),dtype=np.uint64)<br>
         rflt = rflt / float(2**64)<br>
     b = (rflt.reshape((m,8))<p)<br>
     # pack the bits<br>
     b = (b<<range(8)).sum(axis=1).astype(np.uint8)<br>
     # zero the trailing m*8 - n bits<br>
     b[-1] &= (0xFF >> (m*8 - n))<br>
     return b<br>
<br>
<br>
def bitcount(a, bytewise=False):<br>
     """<br>
     Count the number of set bits in an array of np.uint8.<br>
     """<br>
     assert(a.dtype == np.uint8)<br>
     c = a[:,np.newaxis].repeat(8, axis=1) >> range(8) & 0x01<br>
     return (c.sum(axis=1) if bytewise else c.sum())<br>
<br>
<br>
if __name__ == '__main__':<br>
<br>
     b = randombits(int(1e6), .1, intent='numerical')<br>
     print bitcount(b) # should be close to 100000<br>
     b = randombits(int(1e6), .1, intent='crypto')<br>
     print bitcount(b) # should be close to 100000<br>
<br>
<br>
<br>
<br>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
</blockquote></div><br>