diff options
Diffstat (limited to 'Wrappers/Python')
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py | 212 |
1 files changed, 105 insertions, 107 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index 25af156..3c15fa1 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -14,6 +14,12 @@ from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff class SymmetrizedGradient(Gradient): + ''' Symmetrized Gradient, denoted by E: V --> W + where V is the Range of the Gradient Operator + and W is the Range of the Symmetrized Gradient. + ''' + + def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) @@ -37,25 +43,7 @@ class SymmetrizedGradient(Gradient): self.order_ind = [0,2,1,3] else: self.order_ind = [0,3,6,1,4,7,2,5,8] - - -# if self.correlation=='Space': -# if self.channels>1: -# pass -# else: -## # 2D image ---> Dx v1, Dyv2, Dx -# tmp = self.gm_domain.geometries + (self.gm_domain.get_item(0),) -# self.gm_range = BlockGeometry(*tmp ) -# self.ind1 = range(self.gm_domain.get_item(0).length) -# self.ind2 = range(self.gm_domain.get_item(0).length-1, -1, -1) -## self.order = myorder = [0,1,2 3] -# -# elif self.correlation=='SpaceChannels': -# if self.channels>1: -# pass -# else: -# raise ValueError('No channels to correlate') - + def direct(self, x, out=None): @@ -71,10 +59,21 @@ class SymmetrizedGradient(Gradient): res = [0.5 * sum(x) for x in zip(tmp, tmp1)] - return BlockDataContainer(*res) - - + return BlockDataContainer(*res) + else: + + ind = 0 + for i in range(self.gm_domain.shape[0]): + for j in range(x.shape[0]): + self.FD.direction = i + self.FD.adjoint(x.get_item(j), out=out[ind]) + ind+=1 + out1 = [out[i] for i in self.order_ind] + res = [0.5 * sum(x) for x in zip(out, out1)] + out.fill(BlockDataContainer(*res)) + + def adjoint(self, x, out=None): if out is None: @@ -91,20 +90,20 @@ class SymmetrizedGradient(Gradient): tmp[k] = tmp1 return BlockDataContainer(*tmp) -# i = 0 -# -# tmp = self.gm_domain.allocate() -# for k in range(tmp.shape[0]): -# tmp1 = 0 -# for j in range(tmp.shape[0]): -# self.FD.direction = j -# tmp1 += self.FD.direct(x[i]) -# i+=1 -# tmp.get_item(k).fill(tmp1) -# return tmp - + else: - pass + + tmp = self.gm_domain.allocate() + i = 0 + for k in range(self.gm_domain.shape[0]): + tmp1 = 0 + for j in range(self.gm_domain.shape[0]): + self.FD.direction = j + self.FD.direct(x[i], out=tmp[j]) + i+=1 + tmp1+=tmp[j] + out[k].fill(tmp1) + def domain_geometry(self): @@ -117,6 +116,7 @@ class SymmetrizedGradient(Gradient): #TODO need dot method for BlockDataContainer return numpy.sqrt(4*self.gm_domain.shape[0]) + # x0 = self.gm_domain.allocate('random_int') # self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0) # return self.s1 @@ -126,114 +126,112 @@ class SymmetrizedGradient(Gradient): if __name__ == '__main__': ########################################################################### - ## Symmetrized Gradient + ## Symmetrized Gradient Tests from ccpi.framework import DataContainer from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff import numpy as np N, M = 2, 3 K = 2 + C = 2 ig1 = ImageGeometry(N, M) - ig2 = ImageGeometry(N, M, K) + ig2 = ImageGeometry(N, M, channels=C) E1 = SymmetrizedGradient(ig1, correlation = 'Space', bnd_cond='Neumann') - E2 = SymmetrizedGradient(ig2, correlation = 'Space', bnd_cond='Periodic') - print(E1.domain_geometry().shape, E1.range_geometry().shape) - print(E2.domain_geometry().shape, E2.range_geometry().shape) - - u1 = E1.gm_domain.allocate('random_int') - a1 = ig1.allocate('random_int') - a2 = ig1.allocate('random_int') - a3 = ig1.allocate('random_int') - w1 = BlockDataContainer(*[a1, a2, a3]) - w11 = BlockDataContainer(*[a1, a2, a2, a3]) + try: + E1 = SymmetrizedGradient(ig1, correlation = 'SpaceChannels', bnd_cond='Neumann') + except: + print("No Channels to correlate") - sym_direct = [None]*3 - - sym_direct[0] = FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').adjoint(u1[0]) - sym_direct[1] = FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').adjoint(u1[1]) - sym_direct[2] = 0.5 * (FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').adjoint(u1[0]) + \ - FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').adjoint(u1[1])) - - sym_direct1 = [None]*4 - - sym_direct1[0] = FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').adjoint(u1[0]) + 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) - sym_direct1[1] = 0.5 * (FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').adjoint(u1[1]) + \ - FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').adjoint(u1[0])) + #check Linear operator property - sym_direct1[2] = 0.5 * (FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').adjoint(u1[1]) + \ - FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').adjoint(u1[0])) - - sym_direct1[3] = FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').adjoint(u1[1]) + u1 = E1.domain_geometry().allocate('random_int') + u2 = E2.domain_geometry().allocate('random_int') - - - sym_adjoint = [None]*2 + # Need to allocate random_int at the Range of SymGradient - sym_adjoint[0] = FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').direct(w1[0]) + \ - FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').direct(w1[2]) - sym_adjoint[1] = FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').direct(w1[2]) + \ - FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').direct(w1[1]) - - sym_adjoint1 = [None]*2 + #a1 = ig1.allocate('random_int') + #a2 = ig1.allocate('random_int') + #a3 = ig1.allocate('random_int') - sym_adjoint1[0] = FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').direct(w11[0]) + \ - FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').direct(w11[1]) - sym_adjoint1[1] = FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann').direct(w11[2]) + \ - FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann').direct(w11[3]) - + #a4 = ig1.allocate('random_int') + #a5 = ig1.allocate('random_int') + #a6 = ig1.allocate('random_int') - LHS = (sym_direct[0] * w1[0] + \ - sym_direct[1] * w1[1] + \ - 2*sym_direct[2] * w1[2]).sum() + # TODO allocate has to create this symmetry by default!!!!! + #w1 = BlockDataContainer(*[a1, a2, \ + # a2, a3]) + w1 = E1.range_geometry().allocate('random_int',symmetry=True) - RHS = (u1[0]*sym_adjoint[0] + u1[1]*sym_adjoint[1]).sum() + LHS = (E1.direct(u1) * w1).sum() + RHS = (u1 * E1.adjoint(w1)).sum() - print(LHS, RHS) + numpy.testing.assert_equal(LHS, RHS) - LHS = (sym_direct1[0] * w11[0] + \ - sym_direct1[1] * w11[1] + \ - sym_direct1[2] * w11[2] + \ - sym_direct1[3] * w11[3] ).sum() + u2 = E2.gm_domain.allocate('random_int') - RHS = (u1[0]*sym_adjoint1[0] + u1[1]*sym_adjoint1[1]).sum() + #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') - print(LHS, RHS) + #w2 = BlockDataContainer(*[aa1, aa2, aa3, \ + # aa2, aa4, aa5, \ + # aa3, aa5, aa6]) + w2 = E2.range_geometry().allocate('random_int',symmetry=True) + + LHS1 = (E2.direct(u2) * w2).sum() + RHS1 = (u2 * E2.adjoint(w2)).sum() + numpy.testing.assert_equal(LHS1, RHS1) - a1 = (E1.direct(u1) * w11).sum() - b1 = (u1 * E1.adjoint(w11)).sum() - print(a1, b1) + 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()) - u2 = E2.gm_domain.allocate('random_int') + out1 = E1.domain_geometry().allocate() + E1.adjoint(w1, out=out1) + b1 = E1.adjoint(w1) - 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') + LHS_out = (out * w1).sum() + RHS_out = (u1 * out1).sum() + print(LHS_out, RHS_out) - w2 = BlockDataContainer(*[aa1, aa2, aa3, \ - aa2, aa4, aa5, \ - aa3, aa5, aa6]) - tmp1 = E2.direct(u2) - tmp2 = E2.adjoint(w2) - c1 = (tmp1 * w2).sum() - d1 = (u2 * tmp2).sum() + out2 = E2.range_geometry().allocate() + E2.direct(u2, out=out2) + a2 = E2.direct(u2) - print(c1, d1) + 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) -
\ No newline at end of file + +# +# +# +#
\ No newline at end of file |