diff options
Diffstat (limited to 'Wrappers/Python')
6 files changed, 167 insertions, 33 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 54f8859..5489d92 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -31,8 +31,7 @@ class L2NormSquared(Function): super(L2NormSquared, self).__init__() self.b = kwargs.get('b',None) - def __call__(self, x, out=None): - + def __call__(self, x): ''' Evaluates L2NormSq at point x''' y = x @@ -44,33 +43,41 @@ class L2NormSquared(Function): # if out is None: # return x.squared_norm() # else: - return y.squared_norm() - + try: + return y.squared_norm() + except AttributeError as ae: + # added for compatibility with SIRF + return (y.norm()**2) + def gradient(self, x, out=None): - ''' Evaluates gradient of L2NormSq at point x''' - + if out is not None: + out.fill(x) + if self.b is not None: + out -= self.b + out *= 2 + else: + y = x if self.b is not None: # x.subtract(self.b, out=x) y = x - self.b - if out is None: - return 2*y - else: - return out.fill(2*y) + return 2*y + def convex_conjugate(self, x, out=None): - - ''' Evaluate convex conjugate of L2NormSq ''' + ''' Evaluate convex conjugate of L2NormSq''' tmp = 0 if self.b is not None: tmp = (self.b * x).sum() if out is None: + # FIXME: this is a number return (1/4) * x.squared_norm() + tmp else: + # FIXME: this is a DataContainer out.fill((1/4) * x.squared_norm() + tmp) @@ -86,10 +93,15 @@ class L2NormSquared(Function): else: return x/(1+2*tau) else: + out.fill(x) if self.b is not None: - out.fill((x - self.b)/(1+2*tau) + self.b) - else: - out.fill(x/(1+2*tau)) + out -= self.b + out /= (1+2*tau) + if self.b is not None: + out += self.b + #out.fill((x - self.b)/(1+2*tau) + self.b) + #else: + # out.fill(x/(1+2*tau)) def proximal_conjugate(self, x, tau, out=None): diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py index 9def741..8c28874 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py @@ -29,7 +29,7 @@ class ZeroFun(Function): if x.shape[0]==1: return x.maximum(0).sum() else: - if isinstance(x, CompositeDataContainer): + if isinstance(x, BlockDataContainer): return x.get_item(0).maximum(0).sum() + x.get_item(1).maximum(0).sum() else: return x.maximum(0).sum() + x.maximum(0).sum() diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index 1240b31..ee8f609 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -190,9 +190,15 @@ class BlockOperator(Operator): return type(self)(*ops, shape=self.shape) @property def T(self): - '''Return the transposed of self''' - shape = (self.shape[1], self.shape[0]) - return type(self)(*self.operators, shape=shape) + '''Return the transposed of self + + input in a row-by-row''' + newshape = (self.shape[1], self.shape[0]) + oplist = [] + for col in range(newshape[1]): + for row in range(newshape[0]): + oplist.append(self.get_item(col,row)) + return type(self)(*oplist, shape=newshape) def domain_geometry(self): '''returns the domain of the BlockOperator diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py index cdf15a7..2d2089b 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -4,7 +4,7 @@ Created on Tue Mar 5 15:55:56 2019 @author: ofn77899
"""
-from ccpi.optimisation.operators import ScaledOperator
+from ccpi.optimisation.operators.ScaledOperator import ScaledOperator
class Operator(object):
'''Operator that maps from a space X -> Y'''
diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py index 4eb84bb..e1c05fb 100644 --- a/Wrappers/Python/test/test_BlockOperator.py +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -39,6 +39,11 @@ class TestBlockOperator(unittest.TestCase): zero = numpy.zeros(X.get_item(0).shape) numpy.testing.assert_array_equal(Y.get_item(0).as_array(),len(x)+zero) + K2 = BlockOperator(*(ops+ops), shape=(3,2)) + Y = K2.T.direct(X) + # K.T (2,3) X (3,1) => output shape (2,1) + self.assertTrue(Y.shape == (2,1)) + try: # this should fail as the domain is not compatible ig = [ ImageGeometry(10,20,31) , \ diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index 384e351..3e5f26f 100644 --- a/Wrappers/Python/test/test_functions.py +++ b/Wrappers/Python/test/test_functions.py @@ -17,17 +17,20 @@ from ccpi.framework import BlockDataContainer from numbers import Number from ccpi.optimisation.operators import Gradient -from ccpi.optimisation.functions import SimpleL2NormSq -from ccpi.optimisation.functions import L2NormSq +#from ccpi.optimisation.functions import SimpleL2NormSq +from ccpi.optimisation.functions import L2NormSquared from ccpi.optimisation.functions import SimpleL1Norm from ccpi.optimisation.functions import L1Norm + +from ccpi.optimisation.funcs import Norm2sq # from ccpi.optimisation.functions.L2NormSquared import SimpleL2NormSq, L2NormSq # from ccpi.optimisation.functions.L1Norm import SimpleL1Norm, L1Norm -from ccpi.optimisation.functions import mixed_L12Norm +#from ccpi.optimisation.functions import mixed_L12Norm from ccpi.optimisation.functions import ZeroFun from ccpi.optimisation.functions import FunctionOperatorComposition -import unittest +import unittest +import numpy # @@ -46,12 +49,12 @@ class TestFunction(unittest.TestCase): operator = BlockOperator(op1, op2 , shape=(2,1) ) # Create functions - noisy_data = ImageData(np.random.randint(10, size=ag)) + noisy_data = ag.allocate(ImageGeometry.RANDOM_INT) - d = ImageData(np.random.randint(10, size=ag)) + d = ag.allocate(ImageGeometry.RANDOM_INT) alpha = 0.5 # scaled function - g = alpha * L2NormSq(b=noisy_data) + g = alpha * L2NormSquared(b=noisy_data) # Compare call of g a2 = alpha*(d - noisy_data).power(2).sum() @@ -63,13 +66,121 @@ class TestFunction(unittest.TestCase): self.assertEqual(a3, g.convex_conjugate(d)) #print( a3, g.convex_conjugate(d)) + def test_L2NormSquared(self): + # TESTS for L2 and scalar * L2 - def stest_ScaledFunctin(self): - ig = (N,N) - ag = ig - op1 = Gradient(ig) - op2 = Identity(ig, ag) - + M, N, K = 2,3,5 + ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K) + u = ig.allocate(ImageGeometry.RANDOM_INT) + b = ig.allocate(ImageGeometry.RANDOM_INT) + + # check grad/call no data + f = L2NormSquared() + a1 = f.gradient(u) + a2 = 2 * u + numpy.testing.assert_array_almost_equal(a1.as_array(), a2.as_array(), decimal=4) + numpy.testing.assert_equal(f(u), u.squared_norm()) + + # check grad/call with data + f1 = L2NormSquared(b=b) + b1 = f1.gradient(u) + b2 = 2 * (u-b) + + numpy.testing.assert_array_almost_equal(b1.as_array(), b2.as_array(), decimal=4) + numpy.testing.assert_equal(f1(u), (u-b).squared_norm()) + + #check convex conjuagate no data + c1 = f.convex_conjugate(u) + c2 = 1/4 * u.squared_norm() + numpy.testing.assert_equal(c1, c2) + + #check convex conjuagate with data + d1 = f1.convex_conjugate(u) + d2 = (1/4) * u.squared_norm() + (u*b).sum() + numpy.testing.assert_equal(d1, d2) + + # check proximal no data + tau = 5 + e1 = f.proximal(u, tau) + e2 = u/(1+2*tau) + numpy.testing.assert_array_almost_equal(e1.as_array(), e2.as_array(), decimal=4) + + # check proximal with data + tau = 5 + h1 = f1.proximal(u, tau) + h2 = (u-b)/(1+2*tau) + b + numpy.testing.assert_array_almost_equal(h1.as_array(), h2.as_array(), decimal=4) + + # check proximal conjugate no data + tau = 0.2 + k1 = f.proximal_conjugate(u, tau) + k2 = u/(1 + tau/2 ) + numpy.testing.assert_array_almost_equal(k1.as_array(), k2.as_array(), decimal=4) + + # check proximal conjugate with data + l1 = f1.proximal_conjugate(u, tau) + l2 = (u - tau * b)/(1 + tau/2 ) + numpy.testing.assert_array_almost_equal(l1.as_array(), l2.as_array(), decimal=4) + + + # check scaled function properties + + # scalar + scalar = 100 + f_scaled_no_data = scalar * L2NormSquared() + f_scaled_data = scalar * L2NormSquared(b=b) + + # call + numpy.testing.assert_equal(f_scaled_no_data(u), scalar*f(u)) + numpy.testing.assert_equal(f_scaled_data(u), scalar*f1(u)) + + # grad + numpy.testing.assert_array_almost_equal(f_scaled_no_data.gradient(u).as_array(), scalar*f.gradient(u).as_array(), decimal=4) + numpy.testing.assert_array_almost_equal(f_scaled_data.gradient(u).as_array(), scalar*f1.gradient(u).as_array(), decimal=4) + + # conj + numpy.testing.assert_almost_equal(f_scaled_no_data.convex_conjugate(u), \ + f.convex_conjugate(u/scalar) * scalar, decimal=4) + + numpy.testing.assert_almost_equal(f_scaled_data.convex_conjugate(u), \ + scalar * f1.convex_conjugate(u/scalar), decimal=4) + + # proximal + numpy.testing.assert_array_almost_equal(f_scaled_no_data.proximal(u, tau).as_array(), \ + f.proximal(u, tau*scalar).as_array()) + + + numpy.testing.assert_array_almost_equal(f_scaled_data.proximal(u, tau).as_array(), \ + f1.proximal(u, tau*scalar).as_array()) + + + # proximal conjugate + numpy.testing.assert_array_almost_equal(f_scaled_no_data.proximal_conjugate(u, tau).as_array(), \ + (u/(1 + tau/(2*scalar) )).as_array(), decimal=4) + + numpy.testing.assert_array_almost_equal(f_scaled_data.proximal_conjugate(u, tau).as_array(), \ + ((u - tau * b)/(1 + tau/(2*scalar) )).as_array(), decimal=4) + + + def test_Norm2sq_as_FunctionOperatorComposition(self): + M, N, K = 2,3,5 + ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K) + u = ig.allocate(ImageGeometry.RANDOM_INT) + b = ig.allocate(ImageGeometry.RANDOM_INT) + + A = 0.5 * Identity(ig) + old_chisq = Norm2sq(A, b, 1.0) + new_chisq = FunctionOperatorComposition(A, L2NormSquared(b=b)) + + yold = old_chisq(u) + ynew = new_chisq(u) + self.assertEqual(yold, ynew) + + yold = old_chisq.gradient(u) + ynew = new_chisq.gradient(u) + numpy.testing.assert_array_equal(yold.as_array(), ynew.as_array()) + + # # f1 = L2NormSq(alpha=1, b=noisy_data) |