diff options
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py | 7 | ||||
| -rwxr-xr-x | Wrappers/Python/test/test_run_test.py | 23 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/fix_test.py | 66 |
3 files changed, 57 insertions, 39 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py index 8474d89..4faffad 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -50,7 +50,7 @@ class CGLS(Algorithm): def set_up(self, x_init, operator , data ): self.r = data.copy() - self.x = x_init.copy() + self.x = x_init * 0 self.operator = operator self.d = operator.adjoint(self.r) @@ -96,11 +96,12 @@ class CGLS(Algorithm): Ad = self.operator.direct(self.d) norm = Ad.squared_norm() if norm == 0.: - print ('cannot update solution') + print ('norm = 0, cannot update solution') + print ("self.d norm", self.d.squared_norm(), self.d.as_array()) raise StopIteration() alpha = self.normr2/norm if alpha == 0.: - print ('cannot update solution') + print ('alpha = 0, cannot update solution') raise StopIteration() self.d *= alpha Ad *= alpha diff --git a/Wrappers/Python/test/test_run_test.py b/Wrappers/Python/test/test_run_test.py index 4f66da1..a0db9cb 100755 --- a/Wrappers/Python/test/test_run_test.py +++ b/Wrappers/Python/test/test_run_test.py @@ -10,10 +10,9 @@ from ccpi.optimisation.algorithms import FISTA #from ccpi.optimisation.algs import FBPD from ccpi.optimisation.functions import Norm2Sq from ccpi.optimisation.functions import ZeroFunction +# from ccpi.optimisation.funcs import Norm1 from ccpi.optimisation.functions import L1Norm -# This was removed -#from ccpi.optimisation.funcs import Norm2 -#from ccpi.optimisation.funcs import Norm1 +from ccpi.optimisation.funcs import Norm2 from ccpi.optimisation.operators import LinearOperatorMatrix from ccpi.optimisation.operators import Identity @@ -83,7 +82,6 @@ class TestAlgorithms(unittest.TestCase): opt = {'memopt': True} # Create object instances with the test data A and b. f = Norm2Sq(A, b, c=0.5, memopt=True) - f.L = LinearOperator.PowerMethod(A, 10) g0 = ZeroFunction() # Initial guess @@ -126,7 +124,6 @@ class TestAlgorithms(unittest.TestCase): self.assertTrue(cvx_not_installable) def test_FISTA_Norm1_cvx(self): - print ("test_FISTA_Norm1_cvx") if not cvx_not_installable: try: opt = {'memopt': True} @@ -141,8 +138,11 @@ class TestAlgorithms(unittest.TestCase): # A = Identity() # Change n to equal to m. - - b = DataContainer(bmat) + vgb = VectorGeometry(m) + vgx = VectorGeometry(n) + b = vgb.allocate() + b.fill(bmat) + #b = DataContainer(bmat) # Regularization parameter lam = 10 @@ -152,10 +152,11 @@ class TestAlgorithms(unittest.TestCase): g0 = ZeroFunction() # Initial guess - x_init = DataContainer(np.zeros((n, 1))) + #x_init = DataContainer(np.zeros((n, 1))) + x_init = vgx.allocate() # Create 1-norm object instance - g1 = lam * L1Norm() + g1 = Norm1(lam) g1(x_init) g1.prox(x_init, 0.02) @@ -229,7 +230,7 @@ class TestAlgorithms(unittest.TestCase): # Create 1-norm object instance - g1 = lam * L1Norm() + g1 = Norm1(lam) # Compare to CVXPY @@ -296,7 +297,7 @@ class TestAlgorithms(unittest.TestCase): # 1-norm regulariser lam1_denoise = 1.0 - g1_denoise = lam1_denoise * L1Norm() + g1_denoise = Norm1(lam1_denoise) # Initial guess x_init_denoise = ImageData(np.zeros((N, N))) diff --git a/Wrappers/Python/wip/fix_test.py b/Wrappers/Python/wip/fix_test.py index 7f0124f..b1006c0 100755 --- a/Wrappers/Python/wip/fix_test.py +++ b/Wrappers/Python/wip/fix_test.py @@ -61,10 +61,15 @@ class Norm1(Function): opt = {'memopt': True} # Problem data. -m = 10 -n = 10 -np.random.seed(1) +m = 500 +n = 200 + +# if m < n then the problem is under-determined and algorithms will struggle to find a solution. +# One approach is to add regularisation + +#np.random.seed(1) Amat = np.asarray( np.random.randn(m, n), dtype=numpy.float32) +Amat = np.asarray( np.random.random_integers(1,10, (m, n)), dtype=numpy.float32) #Amat = np.asarray(np.eye(m), dtype=np.float32) * 2 A = LinearOperatorMatrix(Amat) bmat = np.asarray( np.random.randn(m), dtype=numpy.float32) @@ -78,8 +83,8 @@ print ("A", A.A) # Change n to equal to m. vgb = VectorGeometry(m) vgx = VectorGeometry(n) -b = vgb.allocate(0, dtype=numpy.float32) -b.fill(bmat) +b = vgb.allocate(VectorGeometry.RANDOM_INT, dtype=numpy.float32) +# b.fill(bmat) #b = DataContainer(bmat) # Regularization parameter @@ -133,14 +138,17 @@ print ("pippo", pippo.as_array()) print ("x_init", x_init.as_array()) print ("x1", x1.as_array()) +y = A.direct(x_init) +y *= 0 +A.direct(x_init, out=y) # Combine with least squares and solve using generic FISTA implementation #x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt) def callback(it, objective, solution): - print (objective, f(solution)) + print ("callback " , it , objective, f(solution)) fa = FISTA(x_init=x_init, f=f, g=g1) -fa.max_iteration = 100 +fa.max_iteration = 1000 fa.update_objective_interval = int( fa.max_iteration / 10 ) fa.run(fa.max_iteration, callback = None, verbose=True) @@ -149,12 +157,16 @@ gd.max_iteration = 5000 gd.update_objective_interval = int( gd.max_iteration / 10 ) gd.run(gd.max_iteration, callback = None, verbose=True) + + cgls = CGLS(x_init= x_initcgls, operator=A, data=b) cgls.max_iteration = 1000 -cgls.update_objective_interval = 100 +cgls.update_objective_interval = int( cgls.max_iteration / 10 ) #cgls.should_stop = stop_criterion(cgls) -cgls.run(1000, callback = callback, verbose=True) +cgls.run(cgls.max_iteration, callback = callback, verbose=True) + + # Print for comparison print("FISTA least squares plus 1-norm solution and objective value:") @@ -166,25 +178,29 @@ print ('FISTA ', A.direct(fa.get_output()).as_array()) print ('GradientDescent', A.direct(gd.get_output()).as_array()) print ('CGLS ', A.direct(cgls.get_output()).as_array()) +cond = numpy.linalg.cond(A.A) -#%% - -import cvxpy as cp -# Construct the problem. -x = cp.Variable(n) -objective = cp.Minimize(cp.sum_squares(A.A*x - bmat)) -prob = cp.Problem(objective) -# The optimal objective is returned by prob.solve(). -result = prob.solve(solver = cp.MOSEK) - -print ('CGLS ', cgls.get_output().as_array()) -print ('CVX ', x.value) - -print ('FISTA ', fa.get_output().as_array()) -print ('GD ', gd.get_output().as_array()) - +print ("cond" , cond) #%% +try: + import cvxpy as cp + # Construct the problem. + x = cp.Variable(n) + objective = cp.Minimize(cp.sum_squares(A.A*x - bmat)) + prob = cp.Problem(objective) + # The optimal objective is returned by prob.solve(). + result = prob.solve(solver = cp.MOSEK) + + print ('CGLS ', cgls.get_output().as_array()) + print ('CVX ', x.value) + + print ('FISTA ', fa.get_output().as_array()) + print ('GD ', gd.get_output().as_array()) +except ImportError as ir: + pass + + #%% |