Changeset 306 for trunk/src/random.py

Timestamp:

04/24/08 15:38:32 (5 years ago)

Author:

apdavison

Message:

The NumpyRNG class now has some extra options for running parallel simulations. If parallel_safe is True (and if the MPI rank and number of processes are correctly supplied), then the sequence of random numbers seen by the simulation script should be independent of the number of MPI processes (although not yet the number of threads).

Files:

• trunk/src/random.py (modified) (3 diffs)

trunk/src/random.py

@@ -56 +56 @@
     """Wrapper for the numpy.random.RandomState class (Mersenne Twister PRNG)."""
     
-    def __init__(self, seed=None):
+    def __init__(self, seed=None, rank=0, num_processes=1, parallel_safe=False):
         AbstractRNG.__init__(self, seed)
         self.rng = numpy.random.RandomState()
@@ -63 +63 @@
         else:
             self.rng.seed()
+        self.rank = rank # MPI rank
+        self.num_processes = num_processes # total number of MPI processes
+        self.parallel_safe = parallel_safe
             
     def __getattr__(self, name):
@@ -72 +75 @@
         
         If n is 1, return a float, if n > 1, return a numpy array,
-        if n < 0, raise an Exception."""
+        if n < 0, raise an Exception."""      
         if n == 0:
-            return numpy.random.rand(0) # We return an empty array
-        if n > 0:
-            return getattr(self.rng, distribution)(size=n, *parameters)
+            rarr = numpy.random.rand(0) # We return an empty array
+        elif n > 0:
+            if self.num_processes > 1 and not self.parallel_safe:
+                # n is the number for the whole model, so if we do not care about
+                # having exactly the same random numbers independent of the
+                # number of processors (m), we only need generate n/m+1 per node
+                # (assuming round-robin distribution of cells between processors)
+                n = n/self.num_processes + 1 
+            rarr = getattr(self.rng, distribution)(size=n, *parameters)
         #elif n == 1:
-        #    return getattr(self.rng, distribution)(size=1, *parameters)[0]
+        #    rarr = getattr(self.rng, distribution)(size=1, *parameters)[0]
         else:
             raise ValueError, "The sample number must be positive"
-
+        if self.parallel_safe and self.num_processes > 1:
+            # strip out the random numbers that should be used on other processors
+            rarr = rarr[numpy.arange(self.rank, len(rarr), self.num_processes)]
+        return rarr
 
 class GSLRNG(AbstractRNG):