added tests

src/pymortests/algorithms/rand_la.py

-# This file is part of the pyMOR project (https://www.pymor.org).
-# Copyright pyMOR developers and contributors. All rights reserved.
-# License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
-
-import numpy as np
-import pytest
-from hypothesis import assume, settings, HealthCheck
-from hypothesis.strategies import sampled_from
-
-from pymor.algorithms.basic import almost_equal
-from pymor.algorithms.rand_la import rrf 
-from pymor.algorithms.basic import contains_zero_vector
-from pymor.core.logger import log_levels
-from pymor.vectorarrays.numpy import NumpyVectorSpace
-from pymortests.base import runmodule
-from pymortests.strategies import given_vector_arrays
-
-methods = [rrf]
-
-def test_random_generalized_svd_without_products_vectorArray(vector_array, method):
-    pass 
-
-def test_random_generalized_svd_without_products_Operator(operator, method):
-    pass 
-
-def test_random_generalized_svd_with_products_Operator(operator, method):
-    pass 
-
-def test_random_generalized_svd_with_products_vectorArray(vector_array, method): 
-    pass 
-
-def test_random_generalized_svd_source_and_range_quadratic(vectorarray, method): 
-    pass 
-
-def test_p_not_to_large(vectorarray, method):
-    pass 
-
-def test_not_many_modes(method):
-    pass 
-
-def test_random_ghep_single_pass(operator, method):
-    pass 
-
-def test_random_ghep_NO_single_pass(operator, method): 
-    pass
-
-
-
-
-if __name__ == "__main__":
-    runmodule(filename=__file__)

src/pymortests/algorithms/svd_va.py

@@ -63,6 +63,7 @@ def test_not_too_many_modes(method):
     U, s, V = method(vec_array, atol=0, rtol=0)
     assert len(U) == len(s) == len(V) == 1
 
-
 if __name__ == "__main__":
     runmodule(filename=__file__)
+
+

src/pymortests/rand_la.py

@@ -3,45 +3,104 @@
 # License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
 
 import numpy as np
+import scipy as sp
 from numpy.random import uniform
 
-from pymor.algorithms.rand_la import rrf, adaptive_rrf
+from pymor.algorithms.rand_la import rrf, adaptive_rrf, random_ghep, random_generalized_svd
 from pymor.operators.numpy import NumpyMatrixOperator
 from pymor.operators.constructions import VectorArrayOperator
 
 
-np.random.seed(0)
-A = uniform(low=-1.0, high=1.0, size=(100, 100))
-A = A.dot(A.T)
-range_product = NumpyMatrixOperator(A)
+def test_adaptive_rrf():
+    np.random.seed(0)
+    A = uniform(low=-1.0, high=1.0, size=(100, 100))
+    A = A @ A.T
+    range_product = NumpyMatrixOperator(A)
+
+    np.random.seed(1)
+    B = uniform(low=-1.0, high=1.0, size=(10, 10))
+    B = B.dot(B.T)
+    source_product = NumpyMatrixOperator(B)
 
-np.random.seed(1)
-A = uniform(low=-1.0, high=1.0, size=(10, 10))
-A = A.dot(A.T)
-source_product = NumpyMatrixOperator(A)
+    C = range_product.range.random(10, seed=10)
+    op = VectorArrayOperator(C)
 
-B = range_product.range.random(10, seed=10)
-op = VectorArrayOperator(B)
+    D = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
+    op_complex = VectorArrayOperator(D)
 
-C = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
-op_complex = VectorArrayOperator(C)
+    Q1 = adaptive_rrf(op, source_product, range_product)
+    assert Q1 in op.range
+
+    Q2 = adaptive_rrf(op_complex, iscomplex=True)
+    assert np.iscomplexobj(Q2.to_numpy())
+    assert Q2 in op.range
 
 
 def test_rrf():
-    Q = rrf(op, source_product, range_product)
-    assert Q in op.range
-    assert len(Q) == 8
+    np.random.seed(2)
+    A = uniform(low=-1.0, high=1.0, size=(100, 100))
+    A = A @ A.T
+    range_product = NumpyMatrixOperator(A)
 
-    Q = rrf(op_complex, iscomplex=True)
-    assert np.iscomplexobj(Q.to_numpy())
-    assert Q in op.range
-    assert len(Q) == 8
+    np.random.seed(3)
+    B = uniform(low=-1.0, high=1.0, size=(10, 10))
+    B = B @B.T
+    source_product = NumpyMatrixOperator(B)
 
+    C = range_product.range.random(10, seed=10)
+    op = VectorArrayOperator(C)
 
-def test_adaptive_rrf():
-    B = adaptive_rrf(op, source_product, range_product)
-    assert B in op.range
+    D = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
+    op_complex = VectorArrayOperator(D)
+
+    Q1 = rrf(op, source_product, range_product)
+    assert Q1 in op.range
+    assert len(Q1) == 8
+
+    Q2 = rrf(op_complex, iscomplex=True)
+    assert np.iscomplexobj(Q2.to_numpy())
+    assert Q2 in op.range
+    assert len(Q2) == 8
+
+
+def test_random_generalized_svd():
+    np.random.seed(4)
+    E = uniform(low=-1.0, high=1.0, size=(5, 5)) 
+    E_op = NumpyMatrixOperator(E)
+
+    modes = 3
+    U, s, Vh = random_generalized_svd(E_op, modes=modes, p=1)
+    U_real, s_real, Vh_real = sp.linalg.svd(E)
+    
+    assert abs(np.linalg.norm(s-s_real[:modes])) <= 1e-2
+    assert len(U) ==  modes
+    assert len(Vh) == modes
+    assert len(s) == modes
+    assert U in E_op.range
+    assert Vh in E_op.source
+
+
+def test_random_ghep():
+    np.random.seed(5)
+    D = uniform(low=-1.0, high=1.0, size=(5, 5)) 
+    D = D @ D.T
+    D_op = NumpyMatrixOperator(D)
+
+    modes = 3
+    w1, V1 = random_ghep(D_op, modes=modes, p=1, single_pass=False)
+    w2, V2 = random_ghep(D_op, modes=modes, p=1, single_pass=True)
+    w_real, V_real = sp.linalg.eigh(D)
+    w_real = w_real[::-1]
+    V_real = V_real[:,::-1]
+
+    assert abs(np.linalg.norm(w1-w_real[:modes])) <= 1e-2
+    assert abs(np.linalg.norm(w2-w_real[:modes])) <= 1
+
+    for i in range(0,modes):
+        assert np.linalg.norm(abs(V1.to_numpy()[i,:]) - abs(V_real[:,i])) <= 1
+    for i in range(0,modes):
+        assert np.linalg.norm(abs(V2.to_numpy()[i,:]) - abs(V_real[:,i])) <= 1
 
-    B = adaptive_rrf(op_complex, iscomplex=True)
-    assert np.iscomplexobj(B.to_numpy())
-    assert B in op.range
+    assert len(w1) == modes
+    assert len(V1) == modes 
+    assert V1.dim == D_op.source.dim