Merge pull request #1736 from pymor/global_random_state

New approach to handling randomness in pyMOR

.codecov.yml

@@ -16,8 +16,11 @@ coverage:
     changes: false
     project:
       pymor:
+        # since we have no policy on this might as well turn this off entirely
+        target: 0%
         flags:
           - gitlab_ci
           - github_actions
     patch:
       pymor:
+        target: 0%

docs/source/substitutions.py

@@ -164,6 +164,8 @@
 
 .. |SymplecticBasis| replace:: :class:`~pymor.algorithms.symplectic.SymplecticBasis`
 .. |CanonicalSymplecticFormOperator| replace:: :class:`~pymor.operators.symplectic.CanonicalSymplecticFormOperator`
+
+.. |RNG| replace:: :class:`random number generator <pymor.tools.random.RNG>`
 '''
 
 substitutions = interfaces + common

src/pymor/algorithms/eigs.py

@@ -9,10 +9,11 @@
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator, InverseOperator
 from pymor.operators.interface import Operator
+from pymor.tools.random import new_rng
 
 
 def eigs(A, E=None, k=3, sigma=None, which='LM', b=None, l=None, maxiter=1000, tol=1e-13,
-         imag_tol=1e-12, complex_pair_tol=1e-12, complex_evp=False, left_evp=False, seed=0):
+         imag_tol=1e-12, complex_pair_tol=1e-12, complex_evp=False, left_evp=False):
     """Approximate a few eigenvalues of a linear |Operator|.
 
     Computes `k` eigenvalues `w` with corresponding eigenvectors `v` which solve
@@ -67,9 +68,6 @@ def eigs(A, E=None, k=3, sigma=None, which='LM', b=None, l=None, maxiter=1000, t
         are real setting this argument to `False` will increase stability and performance.
     left_evp
         If set to `True` compute left eigenvectors else compute right eigenvectors.
-    seed
-        Random seed which is used for computing the initial vector for the Arnoldi
-        iteration.
 
     Returns
     -------
@@ -93,7 +91,8 @@ def eigs(A, E=None, k=3, sigma=None, which='LM', b=None, l=None, maxiter=1000, t
         assert E.source == A.source
 
     if b is None:
-        b = A.source.random(seed=seed)
+        with new_rng(0):
+            b = A.source.random()
     else:
         assert b in A.source
 

src/pymor/algorithms/hapod.py

@@ -10,6 +10,7 @@
 
 from pymor.algorithms.pod import pod
 from pymor.core.logger import getLogger
+from pymor.tools.random import spawn_rng
 
 
 class Node:
@@ -180,7 +181,7 @@ async def hapod_step(node):
 
         if node.children:
             modes, svals, snap_counts = zip(
-                *await asyncio.gather(*(hapod_step(c) for c in node.children))
+                *await asyncio.gather(*(spawn_rng(hapod_step(c)) for c in node.children))
             )
             for m, sv in zip(modes, svals):
                 m.scal(sv)
@@ -221,7 +222,7 @@ def main():
         nonlocal result
         result = asyncio.run(hapod_step(tree))
     result = None
-    hapod_thread = Thread(target=main)
+    hapod_thread = Thread(target=spawn_rng(main))
     hapod_thread.start()
     hapod_thread.join()
     return result
@@ -433,7 +434,7 @@ def __init__(self, executor, max_workers=None):
 
     def submit(self, f, *args):
         future = asyncio.get_event_loop().create_future()
-        self.queue.append((future, f, args))
+        self.queue.append((future, spawn_rng(f), args))
         asyncio.get_event_loop().create_task(self.run_task())
         return future
 

src/pymor/algorithms/lradi.py

@@ -12,18 +12,17 @@
 from pymor.core.defaults import defaults
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator, InverseOperator
-from pymor.tools.random import get_random_state
+from pymor.tools.random import new_rng
 from pymor.vectorarrays.constructions import cat_arrays
 
 
 @defaults('lradi_tol', 'lradi_maxiter', 'lradi_shifts', 'projection_shifts_init_maxiter',
-          'projection_shifts_init_seed', 'projection_shifts_subspace_columns',
+          'projection_shifts_subspace_columns',
           'wachspress_large_ritz_num', 'wachspress_small_ritz_num', 'wachspress_tol')
 def lyap_lrcf_solver_options(lradi_tol=1e-10,
                              lradi_maxiter=500,
                              lradi_shifts='projection_shifts',
                              projection_shifts_init_maxiter=20,
-                             projection_shifts_init_seed=None,
                              projection_shifts_subspace_columns=6,
                              wachspress_large_ritz_num=50,
                              wachspress_small_ritz_num=25,
@@ -40,8 +39,6 @@ def lyap_lrcf_solver_options(lradi_tol=1e-10,
         See :func:`solve_lyap_lrcf`.
     projection_shifts_init_maxiter
         See :func:`projection_shifts_init`.
-    projection_shifts_init_seed
-        See :func:`projection_shifts_init`.
     projection_shifts_subspace_columns
         See :func:`projection_shifts`.
     wachspress_large_ritz_num
@@ -62,7 +59,6 @@ def lyap_lrcf_solver_options(lradi_tol=1e-10,
                       'shift_options':
                       {'projection_shifts': {'type': 'projection_shifts',
                                              'init_maxiter': projection_shifts_init_maxiter,
-                                             'init_seed': projection_shifts_init_seed,
                                              'subspace_columns': projection_shifts_subspace_columns},
                        'wachspress_shifts': {'type': 'wachspress_shifts',
                                              'large_ritz_num': wachspress_large_ritz_num,
@@ -194,14 +190,15 @@ def projection_shifts_init(A, E, B, shift_options):
     shifts
         A |NumPy array| containing a set of stable shift parameters.
     """
-    random_state = get_random_state(seed=shift_options['init_seed'])
+    rng = new_rng(0)
     for i in range(shift_options['init_maxiter']):
         Q = gram_schmidt(B, atol=0, rtol=0)
         shifts = spla.eigvals(A.apply2(Q, Q), E.apply2(Q, Q))
         shifts = shifts[shifts.real < 0]
         if shifts.size == 0:
             # use random subspace instead of span{B} (with same dimensions)
-            B = B.random(len(B), distribution='normal', random_state=random_state)
+            with rng:
+                B = B.random(len(B), distribution='normal')
         else:
             return shifts
     raise RuntimeError('Could not generate initial shifts for low-rank ADI iteration.')

src/pymor/algorithms/lrradi.py

@@ -12,16 +12,15 @@
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator
 from pymor.algorithms.gram_schmidt import gram_schmidt
-from pymor.tools.random import get_random_state
+from pymor.tools.random import new_rng
 
 
 @defaults('lrradi_tol', 'lrradi_maxiter', 'lrradi_shifts', 'hamiltonian_shifts_init_maxiter',
-          'hamiltonian_shifts_init_seed', 'hamiltonian_shifts_subspace_columns')
+          'hamiltonian_shifts_subspace_columns')
 def ricc_lrcf_solver_options(lrradi_tol=1e-10,
                              lrradi_maxiter=500,
                              lrradi_shifts='hamiltonian_shifts',
                              hamiltonian_shifts_init_maxiter=20,
-                             hamiltonian_shifts_init_seed=None,
                              hamiltonian_shifts_subspace_columns=6):
     """Returns available Riccati equation solvers with default solver options.
 
@@ -35,8 +34,6 @@ def ricc_lrcf_solver_options(lrradi_tol=1e-10,
         See :func:`solve_ricc_lrcf`.
     hamiltonian_shifts_init_maxiter
         See :func:`hamiltonian_shifts_init`.
-    hamiltonian_shifts_init_seed
-        See :func:`hamiltonian_shifts_init`.
     hamiltonian_shifts_subspace_columns
         See :func:`hamiltonian_shifts`.
 
@@ -51,7 +48,6 @@ def ricc_lrcf_solver_options(lrradi_tol=1e-10,
                        'shift_options':
                        {'hamiltonian_shifts': {'type': 'hamiltonian_shifts',
                                                'init_maxiter': hamiltonian_shifts_init_maxiter,
-                                               'init_seed': hamiltonian_shifts_init_seed,
                                                'subspace_columns': hamiltonian_shifts_subspace_columns}}}}
 
 
@@ -229,7 +225,7 @@ def hamiltonian_shifts_init(A, E, B, C, shift_options):
     shifts
         A |NumPy array| containing a set of stable shift parameters.
     """
-    random_state = get_random_state(seed=shift_options['init_seed'])
+    rng = new_rng(0)
     for _ in range(shift_options['init_maxiter']):
         Q = gram_schmidt(C, atol=0, rtol=0)
         Ap = A.apply2(Q, Q)
@@ -249,7 +245,8 @@ def hamiltonian_shifts_init(A, E, B, C, shift_options):
         eigpairs = list(filter(lambda e: e[0].real < 0, eigpairs))
         if len(eigpairs) == 0:
             # use random subspace instead of span{C} (with same dimensions)
-            C = C.random(len(C), distribution='normal', random_state=random_state)
+            with rng:
+                C = C.random(len(C), distribution='normal')
             continue
         # find shift with most impact on convergence
         maxval = -1

src/pymor/algorithms/samdp.py

@@ -10,12 +10,13 @@
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator
 from pymor.operators.interface import Operator
+from pymor.tools.random import new_rng
 
 
 @defaults('which', 'tol', 'imagtol', 'conjtol', 'dorqitol', 'rqitol', 'maxrestart', 'krestart', 'init_shifts',
-          'rqi_maxiter', 'seed')
+          'rqi_maxiter')
 def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=1e-6, conjtol=1e-8,
-          dorqitol=1e-4, rqitol=1e-10, maxrestart=100, krestart=20, rqi_maxiter=10, seed=0):
+          dorqitol=1e-4, rqitol=1e-10, maxrestart=100, krestart=20, rqi_maxiter=10):
     """Compute the dominant pole triplets and residues of the transfer function of an LTI system.
 
     This function uses the subspace accelerated dominant pole (SAMDP) algorithm as described in
@@ -70,8 +71,6 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
         Maximum dimension of search space before performing a restart.
     rqi_maxiter
         Maximum number of iterations for the two-sided Rayleigh quotient iteration.
-    seed
-        Random seed which is used for computing the initial shift and random restarts.
 
     Returns
     -------
@@ -106,7 +105,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
     k = 0
     nrestart = 0
     nr_converged = 0
-    np.random.seed(seed)
+    rng = new_rng(0)
 
     X = A.source.empty()
     Q = A.source.empty()
@@ -121,7 +120,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
     poles = np.empty(0)
 
     if init_shifts is None:
-        st = np.random.uniform() * 10.j
+        st = rng.uniform() * 10.j
         shift_nr = 0
         nr_shifts = 0
     else:
@@ -164,7 +163,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
         SH, UR, URt, res = _select_max_eig(H, G, X, V, B_defl, C_defl, which)
 
         if np.all(res < np.finfo(float).eps):
-            st = np.random.uniform() * 10.j
+            st = rng.uniform() * 10.j
             found = False
         else:
             found = True
@@ -287,7 +286,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
                     if found:
                         st = SH[0, 0]
                     else:
-                        st = np.random.uniform() * 10.j
+                        st = rng.uniform() * 10.j
 
                     if shift_nr < nr_shifts:
                         st = shifts[shift_nr]

src/pymor/basic.py

@@ -76,7 +76,7 @@
                                   SOBTReductor)
 from pymor.reductors.sor_irka import SORIRKAReductor
 
-from pymor.tools.random import default_random_state
+from pymor.tools.random import get_rng, new_rng
 
 from pymor.vectorarrays.constructions import cat_arrays
 from pymor.vectorarrays.list import ListVectorSpace

src/pymor/bindings/dunegdt.py

@@ -125,8 +125,8 @@ def real_zero_vector(self):
     def real_full_vector(self, value):
         return DuneXTVector(self.dune_vector_type(self.dim, value))
 
-    def real_random_vector(self, distribution, random_state, **kwargs):
-        values = _create_random_values(self.dim, distribution, random_state, **kwargs)
+    def real_random_vector(self, distribution, **kwargs):
+        values = _create_random_values(self.dim, distribution, **kwargs)
         return self.real_vector_from_numpy(values)
 
     def real_vector_from_numpy(self, data, ensure_copy=False):

src/pymor/bindings/fenics.py

@@ -157,9 +157,9 @@ def real_full_vector(self, value):
         impl += value
         return FenicsVector(impl)
 
-    def real_random_vector(self, distribution, random_state, **kwargs):
+    def real_random_vector(self, distribution, **kwargs):
         impl = df.Function(self.V).vector()
-        values = _create_random_values(impl.local_size(), distribution, random_state, **kwargs)
+        values = _create_random_values(impl.local_size(), distribution, **kwargs)
         impl[:] = np.ascontiguousarray(values)
         return FenicsVector(impl)
 

src/pymor/parallel/ipython.py

@@ -147,7 +147,8 @@ def __init__(self, num_engines=None, **kwargs):
         else:
             self.view = self.client[:]
         self.logger.info(f'Connected to {len(self.view)} engines')
-        self.view.map_sync(_setup_worker, range(len(self.view)))
+        from pymor.tools.random import get_seed_seq
+        self.view.map_sync(_setup_worker, get_seed_seq().spawn(len(self.view)))
         self._remote_objects_created = Counter()
 
         if defaults.defaults_changes() > 0:
@@ -202,15 +203,13 @@ def _worker_call_function(function, loop, args, kwargs):
         return function(*args, **kwargs)
 
 
-def _setup_worker(worker_id):
+def _setup_worker(seed_seq):
     global _remote_objects
     _remote_objects = {}
-    # ensure that each worker starts with a different RandomState
-    from pymor.tools import random
-    import numpy as np
-    state = random.default_random_state()
-    new_state = np.random.RandomState(state.randint(0, 2**16) + worker_id)
-    random._default_random_state = new_state
+    # ensure that each worker starts with a different yet deterministically
+    # initialized rng
+    from pymor.tools.random import new_rng
+    new_rng(seed_seq).install()
 
 
 def _push_object(remote_id, obj):

src/pymor/parallel/mpi.py

@@ -8,6 +8,7 @@
 
 from pymor.parallel.basic import WorkerPoolBase
 from pymor.tools import mpi
+from pymor.tools.random import get_seed_seq
 
 
 class MPIPool(WorkerPoolBase):
@@ -18,6 +19,7 @@ def __init__(self):
         self.logger.info(f'Connected to {mpi.size} ranks')
         self._payload = mpi.call(mpi.function_call_manage, _setup_worker)
         self._apply(os.chdir, os.getcwd())
+        self._map(_setup_rng, [[[s] for s in get_seed_seq().spawn(mpi.size)]])
 
     def __del__(self):
         mpi.call(mpi.remove_object, self._payload)
@@ -90,15 +92,15 @@ def _worker_map_function(payload, function, **kwargs):
 
 
 def _setup_worker():
-    # ensure that each worker starts with a different RandomState
-    if not mpi.rank0:
-        from pymor.tools import random
-        import numpy as np
-        state = random.default_random_state()
-        new_state = np.random.RandomState(state.randint(0, 2**16) + mpi.rank)
-        random._default_random_state = new_state
     return [None]
 
 
+def _setup_rng(seed_seq):
+    # ensure that each worker starts with a different yet deterministically
+    # initialized rng
+    from pymor.tools.random import new_rng
+    new_rng(seed_seq).install()
+
+
 def _push_object(obj):
     return obj