From 7dd757baddcfd3a71b100d55b371ebdf7111d698 Mon Sep 17 00:00:00 2001 From: scivision Date: Tue, 26 Mar 2024 00:58:30 -0400 Subject: [PATCH] matlab: modernize --- +airtools/kaczmarz.m | 6 ++++-- +airtools/logmart.m | 6 +++--- +airtools/maxent.m | 4 ++-- +airtools/unitTest.m | 29 ++++++++++++++++---------- src/airtools/maxent.py | 1 - src/airtools/picard.py | 1 - src/airtools/tests/conftest.py | 4 +++- src/airtools/tests/matlab_engine.py | 13 ------------ src/airtools/tests/test_matlab.py | 32 ++++++++++++++++++++--------- 9 files changed, 52 insertions(+), 44 deletions(-) delete mode 100644 src/airtools/tests/matlab_engine.py diff --git a/+airtools/kaczmarz.m b/+airtools/kaczmarz.m index b0678d0..a250cbc 100644 --- a/+airtools/kaczmarz.m +++ b/+airtools/kaczmarz.m @@ -64,8 +64,10 @@ end % The sizes of A, b and x must match. -if size(b,1) ~= m || size(b,2) ~= 1 - error('The sizes of A and b do not match') +if size(b,1) ~= m + error("rows of A " + int2str(m) + " and columns of b " + int2str(size(b,1)) + " do not match") +elseif size(b,2) ~= 1 + error('b must be a column vector') elseif size(x0,1) ~= n || size(x0,2) ~= 1 error('The size of x0 does not match the problem') end diff --git a/+airtools/logmart.m b/+airtools/logmart.m index 3646d63..f3e3767 100644 --- a/+airtools/logmart.m +++ b/+airtools/logmart.m @@ -30,9 +30,9 @@ y(y<=1e-8) = 1e-8; if isempty(x0) - x=(A'*y)./sum(A(:)); - xA=A*x; - x=x.*max(y(:))/max(xA(:)); + x = (A'*y) ./ sum(A(:)); + xA = A*x; + x = x.*max(y(:)) / max(xA(:)); else x=x0; end diff --git a/+airtools/maxent.m b/+airtools/maxent.m index 679e1cc..2b0afcf 100644 --- a/+airtools/maxent.m +++ b/+airtools/maxent.m @@ -32,7 +32,7 @@ tau0 = 1e-3; % Initial threshold used in secant root finder. % Initialization. -[m,n] = size(A); x_lambda = zeros(n,length(lambda)); F = zeros(maxit,1); +n = size(A,2); x_lambda = zeros(n,length(lambda)); F = zeros(maxit,1); if (min(lambda) <= 0) error('Regularization parameter lambda must be positive') end @@ -51,7 +51,7 @@ % Start the nonlinear CG iteration here. delta_x = x; dF = 1; it = 0; phi0 = p'*g; - while (norm(delta_x) > minstep*norm(x) && dF > flat && it < maxit && phi0 < 0) + while (all(norm(delta_x) > minstep*norm(x), 'all') && dF > flat && it < maxit && all(phi0 < 0, 'all')) it = it + 1; % Compute some CG quantities. diff --git a/+airtools/unitTest.m b/+airtools/unitTest.m index 96c919f..fea69e4 100644 --- a/+airtools/unitTest.m +++ b/+airtools/unitTest.m @@ -50,8 +50,8 @@ function test_condition(tc, name) %% well-posed, well-conditioned problem? A = tc.TestData.(name).A; b = tc.TestData.(name).b; -[U,s] = csvd(A); -picard(U, s, b); +[U,s] = airtools.csvd(A); +airtools.picard(U, s, b); disp("Condition #: " + num2str(cond(A))) end @@ -70,7 +70,7 @@ function test_pseudoinv(tc, name) x_true = tc.TestData.(name).x_true; x_pinv = pinv(A)*b; -tc.verifyEqual(x_pinv, x_true, 'RelTol', 0.005) +tc.verifyEqual(x_pinv, x_true, RelTol=0.005) end function test_logmart(tc, name) @@ -80,19 +80,26 @@ function test_logmart(tc, name) tc.assumeGreaterThanOrEqual(b, 0) -x_logmart = logmart(b,A); -tc.verifyEqual(x_logmart, x_true, 'RelTol', 0.1) +x_logmart = airtools.logmart(b,A); +tc.verifyEqual(x_logmart, x_true, RelTol=0.1) end function test_maxent(tc, name) A = tc.TestData.(name).A; b = tc.TestData.(name).b; x_true = tc.TestData.(name).x_true; -% x_python = py.airtools.maxent.maxent(A,b,0.00002) -% py.numpy.testing.assert_array_almost_equal(x_python,x_true) -x_maxent = maxent(A,b,0.001); -tc.verifyEqual(x_maxent, x_true, 'RelTol', 0.05, 'maxent') +x_maxent = airtools.maxent(A,b,0.001); +tc.verifyEqual(x_maxent, x_true, RelTol=0.05) +end + +function test_maxent_python(tc, name) +A = tc.TestData.(name).A; +b = tc.TestData.(name).b; +x_true = tc.TestData.(name).x_true; + +x_python = py.airtools.maxent.maxent(A,b,0.00002) +tc.verifyEqual(x_python, x_true) end function test_kart(tc,name) @@ -102,8 +109,8 @@ function test_kart(tc,name) % x_python = py.airtools.kaczmarz.kaczmarz(A,b,200)[0] % py.numpy.testing.assert_array_almost_equal(x_python,x_true) % -x_kaczmarz = kaczmarz(A,b,250); -tc.verifyEqual(x_kaczmarz, x_true, 'RelTol', 0.05, 'kaczmarz ART') +x_kaczmarz = airtools.kaczmarz(A,b,250); +tc.verifyEqual(x_kaczmarz, x_true, RelTol=0.05) end end diff --git a/src/airtools/maxent.py b/src/airtools/maxent.py index 69a36df..bb941e5 100755 --- a/src/airtools/maxent.py +++ b/src/airtools/maxent.py @@ -76,7 +76,6 @@ def maxent(A, b, lamb, w=None, x0=None) -> tuple: # Treat each lambda separately. for j in range(Nlambda): - # Prepare for nonlinear CG iteration. l2 = lamb[j] ** 2.0 x = x0 diff --git a/src/airtools/picard.py b/src/airtools/picard.py index 7d01e29..7071ca5 100644 --- a/src/airtools/picard.py +++ b/src/airtools/picard.py @@ -29,7 +29,6 @@ def picard(U, s, b, d=0) -> tuple: - n, ps = np.atleast_2d(s).T.shape beta = np.abs(np.asfortranarray(U[:, :n]).T.dot(b)) diff --git a/src/airtools/tests/conftest.py b/src/airtools/tests/conftest.py index 75a849f..d9e286a 100644 --- a/src/airtools/tests/conftest.py +++ b/src/airtools/tests/conftest.py @@ -15,7 +15,9 @@ def matrices(name: str) -> np.ndarray: [0.126491, 0.357771, 1.41421, 4.0], ] ), - "fiedler": np.array([[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]]), + "fiedler": np.array( + [[0, 1, 2, 3], [1, 0, 1, 2], [2, 1, 0, 1], [3, 2, 1, 0]], dtype=np.float64 + ), "hilbert": np.array( [ [1.0, 1 / 2, 1 / 3, 1 / 4], diff --git a/src/airtools/tests/matlab_engine.py b/src/airtools/tests/matlab_engine.py deleted file mode 100644 index d7b02f9..0000000 --- a/src/airtools/tests/matlab_engine.py +++ /dev/null @@ -1,13 +0,0 @@ -from pathlib import Path -import functools - -import matlab.engine - - -@functools.cache -def matlab_engine(): - """ - only cached because used by Pytest in multiple tests - """ - - return matlab.engine.start_matlab("-nojvm") diff --git a/src/airtools/tests/test_matlab.py b/src/airtools/tests/test_matlab.py index b806693..600a504 100755 --- a/src/airtools/tests/test_matlab.py +++ b/src/airtools/tests/test_matlab.py @@ -1,4 +1,5 @@ from pathlib import Path +import functools import numpy as np @@ -8,7 +9,7 @@ import airtools try: - from .matlab_engine import matlab_engine + import matlab.engine matlab_skip = False except (ImportError, RuntimeError): @@ -26,18 +27,27 @@ used = ("identity", "fiedler") +@functools.cache +def matlab_engine(): + """ + cached to use in multiple tests without restarting + """ + eng = matlab.engine.start_matlab("-nojvm") + eng.addpath(eng.genpath(str(Rmatlab)), nargout=0) + return eng + + @pytest.mark.skipif(matlab_skip, reason="Matlab Engine not available") @pytest.mark.parametrize("name", used) def test_maxent(matrices, name): - eng = matlab_engine() - eng.addpath(eng.genpath(str(Rmatlab)), nargout=0) A = matrices b = A @ x lamb = 2.5e-5 - x_matlab = eng.airtools.maxent(A, b, lamb).squeeze() + x_matlab = eng.airtools.maxent(A, b, lamb) + x_matlab = np.asarray(x_matlab).squeeze() assert x_matlab == approx(x, rel=0.01) @@ -48,9 +58,7 @@ def test_maxent(matrices, name): @pytest.mark.skipif(matlab_skip, reason="Matlab Engine not available") @pytest.mark.parametrize("name", used) def test_kaczmarz(matrices, name): - eng = matlab_engine() - eng.addpath(eng.genpath(str(Rmatlab)), nargout=0) A = matrices b = A @ x @@ -58,7 +66,12 @@ def test_kaczmarz(matrices, name): lamb = 1.0 x0 = np.zeros_like(x) - x_matlab = eng.airtools.kaczmarz(A, b, max_iter, x0, {"lambda": lamb}).squeeze() + print("kaczmarz: A.shape", A.shape) + print("kaczmarz: b.shape", b.shape) + print("kaczmarz: x0.shape", x0.shape) + + x_matlab = eng.airtools.kaczmarz(A, b[:, None], max_iter, x0[:, None], {"lambda": lamb}) + x_matlab = np.asarray(x_matlab).squeeze() assert x_matlab == approx(x, rel=0.01) x_est = airtools.kaczmarz(A, b, x0=x0, max_iter=max_iter, lamb=lamb)[0] @@ -68,9 +81,7 @@ def test_kaczmarz(matrices, name): @pytest.mark.skipif(matlab_skip, reason="Matlab Engine not available") @pytest.mark.parametrize("name", used) def test_logmart(matrices, name): - eng = matlab_engine() - eng.addpath(eng.genpath(str(Rmatlab)), nargout=0) A = matrices @@ -79,7 +90,8 @@ def test_logmart(matrices, name): max_iter = 2000 sigma = 1.0 - x_matlab = eng.airtools.logmart(b, A, relax, [], sigma, max_iter).squeeze() + x_matlab = eng.airtools.logmart(b, A, relax, eng.double.empty(), sigma, max_iter) + x_matlab = np.asarray(x_matlab).squeeze() assert x_matlab == approx(x, rel=0.01) x_est = airtools.logmart(A, b, relax=relax, sigma=sigma, max_iter=max_iter)[0]