From 9a8c03aff07d63f2c635c75fade0da835eb1fb98 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 16 Apr 2019 15:28:55 +0100 Subject: fix new implementation of PowerMethod fixes different values found in #224 --- .../ccpi/optimisation/operators/LinearOperator.py | 39 ++++++---------------- 1 file changed, 11 insertions(+), 28 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py index bc18f9b..f304cc6 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py @@ -24,59 +24,42 @@ class LinearOperator(Operator): raise NotImplementedError @staticmethod - def PowerMethod(operator, iterations, x0=None): + def PowerMethod(operator, iterations, x_init=None): '''Power method to calculate iteratively the Lipschitz constant''' - # Initialise random - if x0 is None: - #x0 = op.create_image_data() + # Initialise random + if x_init is None: x0 = operator.domain_geometry().allocate(ImageGeometry.RANDOM_INT) - + else: + x0 = x_init.copy() + x1 = operator.domain_geometry().allocate() y_tmp = operator.range_geometry().allocate() s = numpy.zeros(iterations) # Loop for it in numpy.arange(iterations): - #x1 = operator.adjoint(operator.direct(x0)) operator.direct(x0,out=y_tmp) operator.adjoint(y_tmp,out=x1) x1norm = x1.norm() - #s[it] = (x1*x0).sum() / (x0.squared_norm()) s[it] = x1.dot(x0) / x0.squared_norm() - #x0 = (1.0/x1norm)*x1 - #x1 *= (1.0 / x1norm) - #x0.fill(x1) x1.multiply((1.0/x1norm), out=x0) return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 @staticmethod - def PowerMethodNonsquare(op,numiters , x0=None): - # Initialise random - # Jakob's - # inputsize , outputsize = op.size() - #x0 = ImageContainer(numpy.random.randn(*inputsize) - # Edo's - #vg = ImageGeometry(voxel_num_x=inputsize[0], - # voxel_num_y=inputsize[1], - # voxel_num_z=inputsize[2]) - # - #x0 = ImageData(geometry = vg, dimension_labels=['vertical','horizontal_y','horizontal_x']) - #print (x0) - #x0.fill(numpy.random.randn(*x0.shape)) + def PowerMethodNonsquare(op,numiters , x_init=None): + '''Power method to calculate iteratively the Lipschitz constant''' - if x0 is None: - #x0 = op.create_image_data() + if x_init is None: x0 = op.allocate_direct() x0.fill(numpy.random.randn(*x0.shape)) + else: + x0 = x_init.copy() s = numpy.zeros(numiters) # Loop for it in numpy.arange(numiters): x1 = op.adjoint(op.direct(x0)) - #x1norm = numpy.sqrt((x1*x1).sum()) x1norm = x1.norm() - #print ("x0 **********" ,x0) - #print ("x1 **********" ,x1) s[it] = (x1*x0).sum() / (x0.squared_norm()) x0 = (1.0/x1norm)*x1 return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 -- cgit v1.2.3