Re: [Numpy-discussion] PRNGs and multi-threading

[Sturla Molden]

Xavier Saint-Mleux skrev:
> Of course, the mathematically correct way would be to use a correct
> jumpahead function, but all the implementations that I know of are GPL. 
> A recent article about this is:
>
> www.iro.umontreal.ca/~lecuyer/myftp/papers/jumpmt.pdf
>
>   
I know of no efficient "jumpahead" function for MT. Several seconds for 
1000 jumps ahead is not impressive -- just generating the deviates is 
faster!

With DCMT it is easy to create "independent" MTs with smaller periods. 
Independence here means that the "characteristic polynomials are 
relatively prime to each other". A "small" period of e.g. 2**521 - 1 
means that if we produce 1 billion deviates per minute, it would still 
take the MT about 10**143 years to cycle. Chances are we will not be 
around to see that happen. It also seems that nvidia has endorsed this 
method:

http://developer.download.nvidia.com/compute/cuda/sdk/website/projects/MersenneTwister/doc/MersenneTwister.pdf



S.M.

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

Re: [Numpy-discussion] PRNGs and multi-threading

[Xavier Saint-Mleux]

Robert Kern wrote:
> On Fri, Aug 21, 2009 at 20:50, Sturla Molden wrote:
>   
>> Sturla Molden skrev:
>> 
>>> It seems there is a special version of the Mersenne Twister for this.
>>> The code is LGPL (annoying for SciPy but ok for me).
>>>   
>> http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dgene.pdf
>> 
>>
>> Basically it encodes the thread-ids in the characteristic polynomial of
>> the MT, producing multiple small-period, independent MTs. That solves it
>> then. Too bad this is LGPL. It would be a very useful enhancement to
>> RandomState.
>> 
>
> I agree. It might be possible to re-implement it from the original
> papers, but it's a chunk of work.
>
>   
I use the following PRNG class, derived from RandomState, which allows a
PRNG to create multiple different sub-PRNGs in a deterministic way:

http://bazaar.launchpad.net/~piaget-dev/piaget/dev/annotate/head%3A/piaget/math/prng.py

It is written in Python and has an Apache license (BSD-like).  There is
no mathematical proof that different PRNGs will have states "far enough"
from each other, but it works well in practice (I've had bad surprises
using Python's random.jumpahead, which is not a real jumpahead).

Of course, the mathematically correct way would be to use a correct
jumpahead function, but all the implementations that I know of are GPL. 
A recent article about this is:

www.iro.umontreal.ca/~lecuyer/myftp/papers/jumpmt.pdf


Xavier Saint-Mleux



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

Re: [Numpy-discussion] PRNGs and multi-threading

[Robert Kern]

On Fri, Aug 21, 2009 at 20:50, Sturla Molden wrote:
> Sturla Molden skrev:
>> It seems there is a special version of the Mersenne Twister for this.
>> The code is LGPL (annoying for SciPy but ok for me).
>
> http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dgene.pdf
> 
>
> Basically it encodes the thread-ids in the characteristic polynomial of
> the MT, producing multiple small-period, independent MTs. That solves it
> then. Too bad this is LGPL. It would be a very useful enhancement to
> RandomState.

I agree. It might be possible to re-implement it from the original
papers, but it's a chunk of work.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] PRNGs and multi-threading

[Sturla Molden]

Sturla Molden skrev:
> It seems there is a special version of the Mersenne Twister for this. 
> The code is LGPL (annoying for SciPy but ok for me).

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dgene.pdf 


Basically it encodes the thread-ids in the characteristic polynomial of 
the MT, producing multiple small-period, independent MTs. That solves it 
then. Too bad this is LGPL. It would be a very useful enhancement to 
RandomState.

Sturla Molden



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

Re: [Numpy-discussion] PRNGs and multi-threading

[Sturla Molden]

Robert Kern skrev:
> As long as you use different seeds, I believe this is fine. The state
> size of MT is so enormous that almost any reasonable use will not find
> overlaps.
>   

It seems there is a special version of the Mersenne Twister for this. 
The code is LGPL (annoying for SciPy but ok for me).

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/DC/dc.html 



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

Re: [Numpy-discussion] PRNGs and multi-threading

[Sturla Molden]

Robert Kern skrev:
> Although you don't really have re-entrancy issues, you will usually
> want one PRNG per thread for determinism.
>   
I see... numpy.random.rand does not have re-entrancy issues because of 
the GIL, but I get indeterminism from the OS scheduling the threads. 
RandomState might not release the GIL either, but preserves determinism 
in presence of multiple threads. Thanks.  :-)
 

Regards,
Sturla Molden







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

Re: [Numpy-discussion] PRNGs and multi-threading

[Robert Kern]

On Fri, Aug 21, 2009 at 17:48, Sturla Molden wrote:
>
> I am not sure if this is the right place to discuss this issue. However,
> I problem I keep running into in Monte Carlo simulations is generating
> pseudorandom numbers with multiple threads. PRNGs such as the Mersenne
> Twister keep an internal state, which prevents the PRNG from being
> re-entrant and thread-safe.

C extension function calls are always atomic unless if they release
the GIL. numpy.random does not. You have de facto locks thanks to the
GIL.

> Possible solutions:
>
>
> 1. Use multiple instances of PRNG states (numpy.random.RandomState), one
> for each thread. This should give no contention, but is this
> mathematically acceptable? I don't know. At least I have not seen any
> proof that it is.

As long as you use different seeds, I believe this is fine. The state
size of MT is so enormous that almost any reasonable use will not find
overlaps.

Although you don't really have re-entrancy issues, you will usually
want one PRNG per thread for determinism.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] PRNGs and multi-threading

[Sturla Molden]

I am not sure if this is the right place to discuss this issue. However, 
I problem I keep running into in Monte Carlo simulations is generating 
pseudorandom numbers with multiple threads. PRNGs such as the Mersenne 
Twister keep an internal state, which prevents the PRNG from being 
re-entrant and thread-safe.


Possible solutions:


1. Use multiple instances of PRNG states (numpy.random.RandomState), one 
for each thread. This should give no contention, but is this 
mathematically acceptable? I don't know. At least I have not seen any 
proof that it is.



2. Protect the PRNG internally with a spinlock. In Windows lingo, that is:

#include 
static volatile long spinlock = 0;
#define ACQUIRE_SPINLOCK while(InterlockedExchangeAcquire(&spinlock, 1));
#define RELEASE_SPINLOCK InterlockedExchangeAcquire(&spinlock, 0);

Problem: possible contention between thread, idle work if threads spins 
a lot.



3. Use a conventional mutex object to protect the PRNG (e.g. 
threading.Lock in Python or CRITICAL_SECTION in Windows).  Problem: 
contention, context shifting, and mutexes tend to be slow. Possibly the 
worst solution.



4. Put the PRNG in a dedicated thread, fill up rather big arrays with 
pseudo-random numbers, and write them to a queue. Problem: Context 
shifting unless a CPU is dedicated to this task. Unless producing random 
numbers consitutes a major portion of the simulation, this should not 
lead to much contention.

import threading
import Queue
import numpy
from numpy.random import rand

class PRNG(threading.Thread):

''' a thread that generate arrays with random numbers
and dumps them to a queue '''

def __init__(self, nthreads, shape):
self.shape = shape
self.count = numpy.prod(shape)
self.queue = Queue.Queue(4 * nthreads) # magic number
self.stop_evt = threading.Event()
   
def generate(self, block=True, timeout=None):
return self.queue.get(block,timeout)
   
def join(self, timeout=None):
self.stop_evt.set()
super(PRNG,self).join(timeout)
   
def run(self):
tmp = rand(count).reshape(shape)
while 1:
try:
self.queue.put(tmp, block=True, timeout=2.0)
except Queue.Full:
if self.stop_evt.isSet(): break
else:
tmp = rand(count).reshape(shape)





Do you have any view on this? Is there any way of creating multiple 
independent random states that will work correctly? I know of SPRNG 
(Scalable PRNG), but it is made to work with MPI (which I don't use).


Regards,
Sturla Molden












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