diff --git a/pytests/test_spgl1.py b/pytests/test_spgl1.py index ca2c823..64b8bfb 100755 --- a/pytests/test_spgl1.py +++ b/pytests/test_spgl1.py @@ -3,7 +3,7 @@ from numpy.testing import assert_array_almost_equal from scipy.sparse import spdiags, csr_matrix -from spgl1.spgl1 import norm_l12nn_primal, norm_l12nn_dual, \ +from spgl1.spgl1 import oneprojector, norm_l12nn_primal, norm_l12nn_dual, \ norm_l12nn_project from spgl1 import spg_lasso, spg_bp, spg_bpdn, spg_mmv @@ -11,33 +11,57 @@ from spgl1.lsqr import lsqr -# dense matrix +# dense matrix (dtype=np.float64) par1 = {'n': 50, 'm': 50, 'k': 14, - 'sparse': False} # square + 'dtype': np.float64, 'sparse': False} # square par2 = {'n': 128, 'm': 50, 'k': 14, - 'sparse': False} # underdetermined + 'dtype': np.float64, 'sparse': False} # underdetermined par3 = {'n': 50, 'm': 100, 'k': 14, - 'sparse': False} # overdetermined + 'dtype': np.float64, 'sparse': False} # overdetermined -# sparse matrix +# sparse matrix (dtype=np.float64) par1_sp = {'n': 50, 'm': 50, 'k': 14, - 'sparse': True} # square + 'dtype': np.float64, 'sparse': True} # square par2_sp = {'n': 128, 'm': 50, 'k': 14, - 'sparse': True} # underdetermined + 'dtype': np.float64, 'sparse': True} # underdetermined par3_sp = {'n': 50, 'm': 100, 'k': 14, - 'sparse': True} # overdetermined + 'dtype': np.float64, 'sparse': True} # overdetermined -np.random.seed(1) +# dense matrix (dtype=np.float32) +par1_f32 = {'n': 50, 'm': 50, 'k': 14, + 'dtype': np.float32, 'sparse': False} # square +par2_f32 = {'n': 128, 'm': 50, 'k': 14, + 'dtype': np.float32, 'sparse': False} # underdetermined +par3_f32 = {'n': 50, 'm': 100, 'k': 14, + 'dtype': np.float32, 'sparse': False} # overdetermined + + +@pytest.mark.parametrize("par", [(par1), (par1_f32)]) +def test_oneprojector(par): + """One projector with tau=1 + + minimize_x ||b-x||_2 subject to ||x||_1 <= 1. + """ + np.random.seed(0) + + b = np.random.normal(0., 1., par['n']).astype(dtype=par['dtype']) + d = 1. + bp = oneprojector(b, d, 1.) + assert bp.dtype == par['dtype'] + assert np.linalg.norm(bp, 1) < 1. @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_lasso(par): """LASSO problem for ||x||_1 <= pi: minimize ||Ax-b||_2 subject to ||x||_1 <= 3.14159... """ + np.random.seed(0) + n = par['n'] m = par['m'] k = par['k'] @@ -47,11 +71,12 @@ def test_lasso(par): A = A1.copy() if par['sparse']: A = csr_matrix(A) + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up vector b, and run solver @@ -59,21 +84,26 @@ def test_lasso(par): tau = np.pi xinv, resid, _, _ = spg_lasso(A, b, tau, verbosity=0) - assert np.linalg.norm(xinv, 1) - np.pi < 1e-10 + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(xinv, 1) - np.pi < 1e-5 # Run solver with subspace minimization xinv, resid, _, _ = spg_lasso(A, b, tau, subspace_min=True, verbosity=0) - assert np.linalg.norm(xinv, 1) - np.pi < 1e-10 + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(xinv, 1) - np.pi < 1e-5 @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_bp(par): """Basis pursuit (BP) problem: minimize ||x||_1 subject to Ax = b """ + np.random.seed(0) + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -82,31 +112,37 @@ def test_bp(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up vector b, and run solver b = A.dot(x) xinv, _, _, _ = spg_bp(A, b, verbosity=0) + assert xinv.dtype == par['dtype'] assert_array_almost_equal(x, xinv, decimal=3) # Run solver with subspace minimization xinv, _, _, _ = spg_bp(A, b, subspace_min=True, verbosity=0) + assert xinv.dtype == par['dtype'] assert_array_almost_equal(x, xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par3), - (par1_sp), (par3_sp)]) + (par1_sp), (par3_sp), + (par1_f32), (par3_f32)]) def test_bpdn(par): """Basis pursuit denoise (BPDN) problem: minimize ||x||_1 subject to ||Ax - b||_2 <= 0.1 """ + np.random.seed(0) + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -115,32 +151,42 @@ def test_bpdn(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up vector b, and run solver - b = A.dot(x) + np.random.randn(n) * 0.075 + b = A.dot(x) + np.random.randn(n).astype(par['dtype']) * 0.075 + sigma = 0.10 xinv, resid, _, _ = spg_bpdn(A, b, sigma, iter_lim=1000, verbosity=0) - assert np.linalg.norm(resid) < sigma*1.1 # need to check why resid is slighly bigger than sigma + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(resid) < sigma * 1.1 # need to check why resid is slighly bigger than sigma # Run solver with subspace minimization xinv, _, _, _ = spg_bpdn(A, b, sigma, subspace_min=True, verbosity=0) - assert np.linalg.norm(resid) < sigma*1.1 # need to check why resid is slighly bigger than sigma + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(resid) < sigma * 1.1 # need to check why resid is slighly bigger than sigma @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_bp_complex(par): """Basis pursuit (BP) problem for complex variables: minimize ||x||_1 subject to Ax = b """ + np.random.seed(0) + + dtypecomplex = (np.ones(1, dtype=par['dtype']) + + 1j * np.ones(1, dtype=par['dtype'])).dtype + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -149,25 +195,29 @@ def test_bp_complex(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(dtypecomplex) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m, dtype=np.complex) + x = np.zeros(m, dtype=dtypecomplex) x[p] = np.random.randn(k) + 1j * np.random.randn(k) # Set up vector b, and run solver b = A.dot(x) xinv, _, _, _ = spg_bp(A, b, verbosity=0) + assert xinv.dtype == dtypecomplex assert_array_almost_equal(x, xinv, decimal=3) # Run solver with subspace minimization xinv, _, _, _ = spg_bp(A, b, subspace_min=True, verbosity=0) + assert xinv.dtype == dtypecomplex assert_array_almost_equal(x, xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_weighted_bp(par): """Weighted Basis pursuit (WBP) problem: @@ -178,6 +228,8 @@ def test_weighted_bp(par): minimize ||Wx||_1 subject to Ax = b """ + np.random.seed(0) + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -186,31 +238,36 @@ def test_weighted_bp(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up weights w and vector b - w = 0.1*np.random.rand(m) + 0.1 # Weights + w = 0.1 * np.random.rand(m).astype(dtype=par['dtype']) + 0.1 # Weights b = A.dot(x / w) # Signal xinv, _, _, _ = spg_bp(A, b, iter_lim=1000, weights=w, verbosity=0) # Reconstructed solution, with weighting xinv *= w + assert xinv.dtype == par['dtype'] assert_array_almost_equal(x, xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_multiple_measurements(par): """Multiple measurement vector (MMV) problem: minimize ||Y||_1,2 subject to AW^{-1}Y = B """ + np.random.seed(0) + # Create random m-by-n encoding matrix m = par['m'] n = par['n'] @@ -218,12 +275,14 @@ def test_multiple_measurements(par): l = 6 A = np.random.randn(m, n) + A = A.astype(par['dtype']) + p = np.random.permutation(n)[:k] - X = np.zeros((n, l)) + X = np.zeros((n, l), dtype=par['dtype']) X[p, :] = np.random.randn(k, l) - weights = 0.1 * np.random.rand(n) + 0.1 - W = 1 / weights * np.eye(n) + weights = 0.1 * np.random.rand(n).astype(par['dtype']) + 0.1 + W = 1 / weights * np.eye(n, dtype=par['dtype']) B = A.dot(W).dot(X) @@ -234,15 +293,20 @@ def test_multiple_measurements(par): X_w, _, _, _ = spg_mmv(A, B, 0, weights=weights, verbosity=0) X_w = spdiags(weights, 0, n, n).dot(X_w) - assert_array_almost_equal(X, X_uw, decimal=2) - assert_array_almost_equal(X, X_w, decimal=2) + assert X_uw.dtype == par['dtype'] + assert X_w.dtype == par['dtype'] + assert_array_almost_equal(X, X_uw, decimal=1) + assert_array_almost_equal(X, X_w, decimal=1) @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_multiple_measurements_nonnegative(par): """Multiple measurement vector (MMV) problem with non-negative norm: """ + np.random.seed(0) + # Create random m-by-n encoding matrix m = par['m'] n = par['n'] @@ -250,8 +314,10 @@ def test_multiple_measurements_nonnegative(par): l = 6 A = np.random.randn(m, n) + A = A.astype(par['dtype']) + p = np.random.permutation(n)[:k] - X = np.zeros((n, l)) + X = np.zeros((n, l), dtype=par['dtype']) X[p, :] = np.abs(np.random.randn(k, l)) B = A.dot(X) @@ -260,6 +326,7 @@ def test_multiple_measurements_nonnegative(par): primal_norm=norm_l12nn_primal, dual_norm=norm_l12nn_dual, iter_lim=20, verbosity=0) + assert Xnn.dtype == par['dtype'] assert np.any(Xnn < 0) == False @@ -276,6 +343,8 @@ def Aprodfun(A, x, mode): return np.dot(np.conj(A.T), x) return y + np.random.seed(0) + # Create random m-by-n encoding matrix m = par['m'] n = par['n'] diff --git a/spgl1/spgl1.py b/spgl1/spgl1.py index 2294b6c..ecec274 100755 --- a/spgl1/spgl1.py +++ b/spgl1/spgl1.py @@ -45,7 +45,7 @@ def __init__(self, A, nnz_idx, ebar, n): self.shape = (A.shape[0], self.nbar) self.dtype = A.dtype def _matvec(self, x): - y = np.zeros(self.n) + y = np.zeros(self.n, dtype=x.dtype) y[self.nnz_idd] = \ x - (1. / self.nbar) * np.dot(np.dot(np.conj(self.ebar), x), self.ebar) @@ -106,7 +106,7 @@ def _oneprojector_i(b, tau): b = b[idx] csb = np.cumsum(b) - tau - alpha = np.zeros(n+1) + alpha = np.zeros(n + 1) alpha[1:] = csb / (np.arange(n) + 1.0) alphaindex = np.where(alpha[1:] >= b)[0] if alphaindex.any(): @@ -160,11 +160,11 @@ def oneprojector(b, d, tau): Projects b onto the (weighted) one-norm ball of radius tau. If d=1 solves the problem:: - minimize_x ||b-x||_2 st ||x||_1 <= tau. + minimize_x ||b-x||_2 subject to ||x||_1 <= tau. else:: - minimize_x ||b-x||_2 st || Dx ||_1 <= tau. + minimize_x ||b-x||_2 subject to ||Dx||_1 <= tau. Parameters ---------- @@ -881,7 +881,7 @@ def spgl1(A, b, tau=0, sigma=0, x0=None, fid=None, verbosity=0, # Determine initial x and see if problem is complex realx = np.lib.isreal(A).all() and np.lib.isreal(b).all() if x0 is None: - x = np.zeros(n) + x = np.zeros(n, dtype=b.dtype) else: x = np.asarray(x0) @@ -1174,7 +1174,7 @@ def spgl1(A, b, tau=0, sigma=0, x0=None, fid=None, verbosity=0, # LSQR iterations successful. Take the subspace step. if istop != 4: # Push dx back into full space: dx = Z dx. - dx = np.zeros(n) + dx = np.zeros(n, dtype=x.dtype) dx[nnz_idx] = \ dxbar - (1/nebar)*np.dot(np.dot(np.conj(ebar.T), dxbar), dxbar) @@ -1288,7 +1288,7 @@ def spgl1(A, b, tau=0, sigma=0, x0=None, fid=None, verbosity=0, elif stat == EXIT_LEAST_SQUARES: _printf(fid, 'EXIT -- Found a least-squares solution') elif stat == EXIT_LINE_ERROR: - _printf(fid, 'ERROR EXIT -- Linesearch error (%i)' % lnerr) + _printf(fid, 'ERROR EXIT -- Linesearch error (%d)' % lnerr) elif stat == EXIT_SUBOPTIMAL_BP: _printf(fid, 'EXIT -- Found a suboptimal BP solution') elif stat == EXIT_MATVEC_LIMIT: