diff options
Diffstat (limited to 'Wrappers/Python')
6 files changed, 160 insertions, 112 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index d71f597..a6c9545 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -240,4 +240,4 @@ if __name__ == '__main__': - + diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py index 09e550e..1c2c43f 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -177,4 +177,4 @@ if __name__ == '__main__': - +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 92e0116..f5108ba 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -280,4 +280,4 @@ if __name__ == '__main__': numpy.testing.assert_array_almost_equal(res1.as_array(), \ res2.as_array(), decimal=4) - +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py index 6e2f49a..0dd7d4b 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -194,7 +194,7 @@ class FiniteDiff(LinearOperator): def adjoint(self, x, out=None): - + x_asarr = x.as_array() x_sz = len(x.shape) outnone = False diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index 8d14cf8..c85abfa 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -21,15 +21,12 @@ from __future__ import print_function from __future__ import unicode_literals -from ccpi.optimisation.operators import Gradient, Operator, LinearOperator,\ - ScaledOperator -from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, \ - BlockDataContainer -import numpy +from ccpi.optimisation.operators import LinearOperator +from ccpi.framework import BlockGeometry, BlockDataContainer from ccpi.optimisation.operators import FiniteDiff -class SymmetrizedGradient(Gradient): +class SymmetrizedGradient(LinearOperator): r'''Symmetrized Gradient Operator: E: V -> W @@ -49,21 +46,23 @@ class SymmetrizedGradient(Gradient): | ''' + CORRELATION_SPACE = "Space" + CORRELATION_SPACECHANNEL = "SpaceChannels" def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): - super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) + super(SymmetrizedGradient, self).__init__() - self.gm_domain = self.gm_range + self.gm_domain = gm_domain self.bnd_cond = bnd_cond - - self.channels = self.gm_range.get_item(0).channels - + self.correlation = kwargs.get('correlation',SymmetrizedGradient.CORRELATION_SPACE) + tmp_gm = len(self.gm_domain.geometries)*self.gm_domain.geometries self.gm_range = BlockGeometry(*tmp_gm) - self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond) + # Define FD operator. We need one geometry from the BlockGeometry of the domain + self.FD = FiniteDiff(self.gm_domain.get_item(0), direction = 0, bnd_cond = self.bnd_cond) if self.gm_domain.shape[0]==2: self.order_ind = [0,2,1,3] @@ -130,8 +129,6 @@ class SymmetrizedGradient(Gradient): i+=1 tmp1+=tmp[j] out[k].fill(tmp1) -# tmp = self.adjoint(x) -# out.fill(tmp) def domain_geometry(self): @@ -144,110 +141,59 @@ if __name__ == '__main__': ########################################################################### ## Symmetrized Gradient Tests - from ccpi.framework import DataContainer - from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff + from ccpi.framework import ImageGeometry + from ccpi.optimisation.operators import Gradient import numpy as np N, M = 2, 3 K = 2 C = 2 - ig1 = ImageGeometry(N, M) - ig2 = ImageGeometry(N, M, channels=C) + ########################################################################### + # 2D geometry no channels + ig = ImageGeometry(N, M) + Grad = Gradient(ig) - E1 = SymmetrizedGradient(ig1, correlation = 'Space', bnd_cond='Neumann') + E1 = SymmetrizedGradient(Grad.range_geometry()) + np.testing.assert_almost_equal(E1.norm(), np.sqrt(8), 1e-5) - try: - E1 = SymmetrizedGradient(ig1, correlation = 'SpaceChannels', bnd_cond='Neumann') - except: - print("No Channels to correlate") - - E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Neumann') - print(E1.domain_geometry().shape, E1.range_geometry().shape) - print(E2.domain_geometry().shape, E2.range_geometry().shape) - - #check Linear operator property - u1 = E1.domain_geometry().allocate('random_int') - u2 = E2.domain_geometry().allocate('random_int') - - # Need to allocate random_int at the Range of SymGradient - - #a1 = ig1.allocate('random_int') - #a2 = ig1.allocate('random_int') - #a3 = ig1.allocate('random_int') - - #a4 = ig1.allocate('random_int') - #a5 = ig1.allocate('random_int') - #a6 = ig1.allocate('random_int') - - # TODO allocate has to create this symmetry by default!!!!! - #w1 = BlockDataContainer(*[a1, a2, \ - # a2, a3]) - w1 = E1.range_geometry().allocate('random_int',symmetry=True) - - LHS = (E1.direct(u1) * w1).sum() - RHS = (u1 * E1.adjoint(w1)).sum() - - numpy.testing.assert_equal(LHS, RHS) + w1 = E1.range_geometry().allocate('random_int', symmetry = True) + - u2 = E2.gm_domain.allocate('random_int') + lhs = E1.direct(u1).dot(w1) + rhs = u1.dot(E1.adjoint(w1)) + np.testing.assert_almost_equal(lhs, rhs) - #aa1 = ig2.allocate('random_int') - #aa2 = ig2.allocate('random_int') - #aa3 = ig2.allocate('random_int') - #aa4 = ig2.allocate('random_int') - #aa5 = ig2.allocate('random_int') - #aa6 = ig2.allocate('random_int') + ########################################################################### + # 2D geometry with channels + ig2 = ImageGeometry(N, M, channels = C) + Grad2 = Gradient(ig2, correlation = 'Space') - #w2 = BlockDataContainer(*[aa1, aa2, aa3, \ - # aa2, aa4, aa5, \ - # aa3, aa5, aa6]) - w2 = E2.range_geometry().allocate('random_int',symmetry=True) - + E2 = SymmetrizedGradient(Grad2.range_geometry()) + np.testing.assert_almost_equal(E2.norm(), np.sqrt(12), 1e-6) - LHS1 = (E2.direct(u2) * w2).sum() - RHS1 = (u2 * E2.adjoint(w2)).sum() - - numpy.testing.assert_equal(LHS1, RHS1) - - out = E1.range_geometry().allocate() - E1.direct(u1, out=out) - a1 = E1.direct(u1) - numpy.testing.assert_array_equal(a1[0].as_array(), out[0].as_array()) - numpy.testing.assert_array_equal(a1[1].as_array(), out[1].as_array()) - numpy.testing.assert_array_equal(a1[2].as_array(), out[2].as_array()) - numpy.testing.assert_array_equal(a1[3].as_array(), out[3].as_array()) - - - out1 = E1.domain_geometry().allocate() - E1.adjoint(w1, out=out1) - b1 = E1.adjoint(w1) - - LHS_out = (out * w1).sum() - RHS_out = (u1 * out1).sum() - print(LHS_out, RHS_out) - - - out2 = E2.range_geometry().allocate() - E2.direct(u2, out=out2) - a2 = E2.direct(u2) - - out21 = E2.domain_geometry().allocate() - E2.adjoint(w2, out=out21) - b2 = E2.adjoint(w2) - - LHS_out = (out2 * w2).sum() - RHS_out = (u2 * out21).sum() - print(LHS_out, RHS_out) - - - out = E1.range_geometry().allocate() - E1.direct(u1, out=out) - E1.adjoint(out, out=out1) - - print(E1.norm()) - print(E2.norm()) - - + print(E2.domain_geometry().shape, E2.range_geometry().shape) + u2 = E2.domain_geometry().allocate('random_int') + w2 = E2.range_geometry().allocate('random_int', symmetry = True) +# + lhs2 = E2.direct(u2).dot(w2) + rhs2 = u2.dot(E2.adjoint(w2)) + np.testing.assert_almost_equal(lhs2, rhs2) + + ########################################################################### + # 3D geometry no channels + ig3 = ImageGeometry(N, M, K) + Grad3 = Gradient(ig3, correlation = 'Space') + + E3 = SymmetrizedGradient(Grad3.range_geometry()) + np.testing.assert_almost_equal(E3.norm(), np.sqrt(12), 1e-6) + + print(E3.domain_geometry().shape, E3.range_geometry().shape) + u3 = E3.domain_geometry().allocate('random_int') + w3 = E3.range_geometry().allocate('random_int', symmetry = True) +# + lhs3 = E3.direct(u3).dot(w3) + rhs3 = u3.dot(E3.adjoint(w3)) + np.testing.assert_almost_equal(lhs3, rhs3)
\ No newline at end of file diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index b26bb5d..3372b9b 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -18,13 +18,14 @@ import unittest from ccpi.framework import ImageGeometry, ImageData, BlockDataContainer, DataContainer from ccpi.optimisation.operators import BlockOperator, BlockScaledOperator,\ - FiniteDiff + FiniteDiff, SymmetrizedGradient import numpy from timeit import default_timer as timer from ccpi.framework import ImageGeometry from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff from ccpi.optimisation.operators import LinearOperator + def dt(steps): return steps[-1] - steps[-2] @@ -170,6 +171,107 @@ class TestOperator(CCPiTestClass): print ("Norm dT1 {} dT2 {}".format(t1-t0,t2-t1)) self.assertLess(t2-t1, t1-t0) +class TestGradients(CCPiTestClass): + def setUp(self): + N, M = 20, 30 + K = 20 + C = 20 + self.decimal = 1 + self.iterations = 50 + ########################################################################### + # 2D geometry no channels + self.ig = ImageGeometry(N, M) + self.ig2 = ImageGeometry(N, M, channels = C) + self.ig3 = ImageGeometry(N, M, K) + + def test_SymmetrizedGradient1a(self): + ########################################################################### + ## Symmetrized Gradient Tests + print ("Test SymmetrizedGradient") + ########################################################################### + # 2D geometry no channels + # ig = ImageGeometry(N, M) + Grad = Gradient(self.ig) + + E1 = SymmetrizedGradient(Grad.range_geometry()) + numpy.testing.assert_almost_equal(E1.norm(iterations=self.iterations), numpy.sqrt(8), decimal = self.decimal) + + def test_SymmetrizedGradient1b(self): + ########################################################################### + ## Symmetrized Gradient Tests + print ("Test SymmetrizedGradient") + ########################################################################### + # 2D geometry no channels + # ig = ImageGeometry(N, M) + Grad = Gradient(self.ig) + + E1 = SymmetrizedGradient(Grad.range_geometry()) + u1 = E1.domain_geometry().allocate('random_int') + w1 = E1.range_geometry().allocate('random_int', symmetry = True) + + + lhs = E1.direct(u1).dot(w1) + rhs = u1.dot(E1.adjoint(w1)) + # self.assertAlmostEqual(lhs, rhs) + numpy.testing.assert_almost_equal(lhs, rhs) + + def test_SymmetrizedGradient2(self): + ########################################################################### + # 2D geometry with channels + # ig2 = ImageGeometry(N, M, channels = C) + Grad2 = Gradient(self.ig2, correlation = 'Space') + + E2 = SymmetrizedGradient(Grad2.range_geometry()) + + u2 = E2.domain_geometry().allocate('random_int') + w2 = E2.range_geometry().allocate('random_int', symmetry = True) + # + lhs2 = E2.direct(u2).dot(w2) + rhs2 = u2.dot(E2.adjoint(w2)) + + numpy.testing.assert_almost_equal(lhs2, rhs2) + + def test_SymmetrizedGradient2a(self): + ########################################################################### + # 2D geometry with channels + # ig2 = ImageGeometry(N, M, channels = C) + Grad2 = Gradient(self.ig2, correlation = 'Space') + + E2 = SymmetrizedGradient(Grad2.range_geometry()) + numpy.testing.assert_almost_equal(E2.norm(iterations=self.iterations), + numpy.sqrt(8), decimal = self.decimal) + + + def test_SymmetrizedGradient3a(self): + ########################################################################### + # 3D geometry no channels + #ig3 = ImageGeometry(N, M, K) + Grad3 = Gradient(self.ig3, correlation = 'Space') + + E3 = SymmetrizedGradient(Grad3.range_geometry()) + + norm1 = E3.norm() + norm2 = E3.calculate_norm(iterations=100) + print (norm1,norm2) + numpy.testing.assert_almost_equal(norm2, numpy.sqrt(12), decimal = self.decimal) + + def test_SymmetrizedGradient3b(self): + ########################################################################### + # 3D geometry no channels + #ig3 = ImageGeometry(N, M, K) + Grad3 = Gradient(self.ig3, correlation = 'Space') + + E3 = SymmetrizedGradient(Grad3.range_geometry()) + + u3 = E3.domain_geometry().allocate('random_int') + w3 = E3.range_geometry().allocate('random_int', symmetry = True) + # + lhs3 = E3.direct(u3).dot(w3) + rhs3 = u3.dot(E3.adjoint(w3)) + numpy.testing.assert_almost_equal(lhs3, rhs3) + # self.assertAlmostEqual(lhs3, rhs3) + # self.assertTrue( LinearOperator.dot_test(E3) ) + |