add random_state and seed arguments to rrf and adaptive_rrf

src/pymor/algorithms/randrangefinder.py

@@ -9,11 +9,12 @@
 from pymor.algorithms.gram_schmidt import gram_schmidt
 from pymor.core.defaults import defaults
 from pymor.operators.interface import Operator
+from pymor.tools.random import get_random_state
 
 
 @defaults('tol', 'failure_tolerance', 'num_testvecs')
-def adaptive_rrf(A, source_product=None, range_product=None, tol=1e-4,
-                 failure_tolerance=1e-15, num_testvecs=20, lambda_min=None, iscomplex=False):
+def adaptive_rrf(A, source_product=None, range_product=None, tol=1e-4, failure_tolerance=1e-15, num_testvecs=20,
+                 lambda_min=None, iscomplex=False, random_state=None, seed=None):
     r"""Adaptive randomized range approximation of `A`.
 
     This is an implementation of Algorithm 1 in :cite:`BS18`.
@@ -47,6 +48,13 @@ def adaptive_rrf(A, source_product=None, range_product=None, tol=1e-4,
         If `None`, the smallest eigenvalue is computed using scipy.
     iscomplex
         If `True`, the random vectors are chosen complex.
+    random_state
+        :class:`~numpy.random.RandomState` to use for sampling.
+        If `None`, a new random state is generated using `seed`
+        as random seed, or the :func:`default <pymor.tools.random.default_random_state>`
+        random state is used.
+    seed
+        If not `None`, a new random state with this seed is used.
 
     Returns
     -------
@@ -57,11 +65,16 @@ def adaptive_rrf(A, source_product=None, range_product=None, tol=1e-4,
     assert range_product is None or isinstance(range_product, Operator)
     assert isinstance(A, Operator)
 
+    assert random_state is None or seed is None
+    random_state = get_random_state(random_state, seed)
+
     B = A.range.empty()
 
-    R = A.source.random(num_testvecs, distribution='normal')
     if iscomplex:
-        R += 1j*A.source.random(num_testvecs, distribution='normal')
+        R = A.source.random(2*num_testvecs, distribution='normal', random_state=random_state)
+        R = R[:num_testvecs] + 1j * R[num_testvecs:]
+    else:
+        R = A.source.random(num_testvecs, distribution='normal', random_state=random_state)
 
     if source_product is None:
         lambda_min = 1
@@ -82,9 +95,11 @@ def mvinv(v):
 
     while(maxnorm > testlimit):
         basis_length = len(B)
-        v = A.source.random(distribution='normal')
         if iscomplex:
-            v += 1j*A.source.random(distribution='normal')
+            v = A.source.random(2, distribution='normal', random_state=random_state)
+            v = v[0] + 1j * v[1]
+        else:
+            v = A.source.random(distribution='normal', random_state=random_state)
         B.append(A.apply(v))
         gram_schmidt(B, range_product, atol=0, rtol=0, offset=basis_length, copy=False)
         M -= B.lincomb(B.inner(M, range_product).T)
@@ -94,7 +109,7 @@ def mvinv(v):
 
 
 @defaults('q', 'l')
-def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False):
+def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False, random_state=None, seed=None):
     """Randomized range approximation of `A`.
 
     This is an implementation of Algorithm 4.4 in :cite:`HMT11`.
@@ -116,6 +131,13 @@ def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False):
         The block size of the normalized power iterations.
     iscomplex
         If `True`, the random vectors are chosen complex.
+    random_state
+        :class:`~numpy.random.RandomState` to use for sampling.
+        If `None`, a new random state is generated using `seed`
+        as random seed, or the :func:`default <pymor.tools.random.default_random_state>`
+        random state is used.
+    seed
+        If not `None`, a new random state with this seed is used.
 
     Returns
     -------
@@ -126,9 +148,14 @@ def rrf(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False):
     assert range_product is None or isinstance(range_product, Operator)
     assert isinstance(A, Operator)
 
-    R = A.source.random(l, distribution='normal')
+    assert random_state is None or seed is None
+    random_state = get_random_state(random_state, seed)
+
     if iscomplex:
-        R += 1j*A.source.random(l, distribution='normal')
+        R = A.source.random(2*l, distribution='normal', random_state=random_state)
+        R = R[:l] + 1j * R[l:]
+    else:
+        R = A.source.random(l, distribution='normal', random_state=random_state)
     Q = A.apply(R)
     gram_schmidt(Q, range_product, atol=0, rtol=0, copy=False)
 

src/pymor/vectorarrays/interface.py

@@ -153,7 +153,7 @@ def random(self, count=1, distribution='uniform', random_state=None, seed=None,
             as random seed, or the :func:`default <pymor.tools.random.default_random_state>`
             random state is used.
         seed
-            If not `None`, a new radom state with this seed is used.
+            If not `None`, a new random state with this seed is used.
         reserve
             Hint for the backend to which length the array will grow.
         """