diff options
Diffstat (limited to 'Wrappers/Python')
-rw-r--r-- | Wrappers/Python/ccpi/framework.py | 14 | ||||
-rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algs.py | 5 | ||||
-rwxr-xr-x | Wrappers/Python/ccpi/optimisation/funcs.py | 2 | ||||
-rwxr-xr-x | Wrappers/Python/ccpi/optimisation/ops.py | 5 | ||||
-rwxr-xr-x | Wrappers/Python/conda-recipe/run_test.py | 28 |
5 files changed, 38 insertions, 16 deletions
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 318eb92..e1a2dff 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -3,7 +3,7 @@ # Visual Analytics and Imaging System Group of the Science Technology # Facilities Council, STFC -# Copyright 2018 Edoardo Pasca +# Copyright 2018-2019 Edoardo Pasca # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -26,6 +26,7 @@ import numpy import sys from datetime import timedelta, datetime import warnings +from functools import reduce def find_key(dic, val): """return the key of dictionary dic given the value""" @@ -720,10 +721,15 @@ class DataContainer(object): ## reductions def sum(self, out=None, *args, **kwargs): return self.as_array().sum(*args, **kwargs) - def norm(self): - '''return the norm of the DataContainer''' - y = self.as_array() + def squared_norm(self): + '''return the squared euclidean norm of the DataContainer viewed as a vector''' + shape = self.shape + size = reduce(lambda x,y:x*y, shape, 1) + y = numpy.reshape(self.as_array(), (size, )) return numpy.dot(y, y.conjugate()) + def norm(self): + '''return the euclidean norm of the DataContainer viewed as a vector''' + return numpy.sqrt(self.squared_norm()) class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index a736277..15638a9 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -158,9 +158,8 @@ def FBPD(x_init, operator=None, constraint=None, data_fidelity=None,\ memopt = opt['memopts'] if 'memopts' in opt.keys() else False # step-sizes - tau = 2 / (data_fidelity.L + 2); - sigma = (1/tau - data_fidelity.L/2) / regulariser.L; - + tau = 2 / (data_fidelity.L + 2) + sigma = (1/tau - data_fidelity.L/2) / regulariser.L inv_sigma = 1/sigma # initialization diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py index c7a6143..9b9fc36 100755 --- a/Wrappers/Python/ccpi/optimisation/funcs.py +++ b/Wrappers/Python/ccpi/optimisation/funcs.py @@ -217,10 +217,10 @@ class ZeroFun(Function): class Norm1(Function): def __init__(self,gamma): + super(Norm1, self).__init__() self.gamma = gamma self.L = 1 self.sign_x = None - super(Norm1, self).__init__() def __call__(self,x,out=None): if out is None: diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index 450b084..b2f996d 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -209,10 +209,11 @@ def PowerMethodNonsquare(op,numiters , x0=None): # Loop for it in numpy.arange(numiters): x1 = op.adjoint(op.direct(x0)) - x1norm = numpy.sqrt((x1*x1).sum()) + #x1norm = numpy.sqrt((x1*x1).sum()) + x1norm = x1.norm() #print ("x0 **********" ,x0) #print ("x1 **********" ,x1) - s[it] = (x1*x0).sum() / (x0*x0).sum() + s[it] = (x1*x0).sum() / (x0.squared_norm()) x0 = (1.0/x1norm)*x1 return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py index 5bf6538..b52af55 100755 --- a/Wrappers/Python/conda-recipe/run_test.py +++ b/Wrappers/Python/conda-recipe/run_test.py @@ -19,6 +19,7 @@ from ccpi.optimisation.funcs import Norm2 from ccpi.optimisation.ops import LinearOperatorMatrix
from ccpi.optimisation.ops import TomoIdentity
from ccpi.optimisation.ops import Identity
+from ccpi.optimisation.ops import PowerMethodNonsquare
import numpy.testing
@@ -494,6 +495,9 @@ class TestDataContainer(unittest.TestCase): except AssertionError as err:
res = False
print(err)
+ print("expected " , second)
+ print("actual " , first)
+
self.assertTrue(res)
def test_DataContainerChaining(self):
dc = self.create_DataContainer(256,256,256,1)
@@ -501,7 +505,13 @@ class TestDataContainer(unittest.TestCase): dc.add(9,out=dc)\
.subtract(1,out=dc)
self.assertEqual(1+9-1,dc.as_array().flatten()[0])
-
+ def test_reduction(self):
+ print ("test reductions")
+ dc = self.create_DataContainer(2,2,2,value=1)
+ sqnorm = dc.squared_norm()
+ norm = dc.norm()
+ self.assertEqual(sqnorm, 8.0)
+ numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7)
@@ -643,7 +653,8 @@ class TestAlgorithms(unittest.TestCase): else:
self.assertTrue(cvx_not_installable)
- def test_FBPD_Norm1_cvx(self):
+ def skip_test_FBPD_Norm1_cvx(self):
+ print ("test_FBPD_Norm1_cvx")
if not cvx_not_installable:
opt = {'memopt': True}
# Problem data.
@@ -663,12 +674,15 @@ class TestAlgorithms(unittest.TestCase): # Regularization parameter
lam = 10
opt = {'memopt': True}
+ # Initial guess
+ x_init = DataContainer(np.random.randn(n, 1))
+
# Create object instances with the test data A and b.
f = Norm2sq(A, b, c=0.5, memopt=True)
+ f.L = PowerMethodNonsquare(A, 25, x_init)[0]
+ print ("Lipschitz", f.L)
g0 = ZeroFun()
- # Initial guess
- x_init = DataContainer(np.zeros((n, 1)))
# Create 1-norm object instance
g1 = Norm1(lam)
@@ -710,7 +724,7 @@ class TestAlgorithms(unittest.TestCase): # Set up phantom size NxN by creating ImageGeometry, initialising the
# ImageData object with this geometry and empty array and finally put some
# data into its array, and display as image.
- def test_FISTA_denoise_cvx(self):
+ def skip_test_FISTA_denoise_cvx(self):
if not cvx_not_installable:
opt = {'memopt': True}
N = 64
@@ -733,6 +747,8 @@ class TestAlgorithms(unittest.TestCase): # Data fidelity term
f_denoise = Norm2sq(I, y, c=0.5, memopt=True)
+ x_init = ImageData(geometry=ig)
+ f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0]
# 1-norm regulariser
lam1_denoise = 1.0
@@ -759,7 +775,7 @@ class TestAlgorithms(unittest.TestCase): # Compare to CVXPY
# Construct the problem.
- x1_denoise = Variable(N**2, 1)
+ x1_denoise = Variable(N**2)
objective1_denoise = Minimize(
0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1))
prob1_denoise = Problem(objective1_denoise)
|