diff options
14 files changed, 33 insertions, 98 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index 10fe34c..1577f06 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -104,8 +104,8 @@ class BlockOperator(Operator): index = row*self.shape[1]+col return self.operators[index] - def calculate_norm(self): - norm = [op.norm()**2 for op in self.operators] + def calculate_norm(self, **kwargs): + norm = [op.norm(**kwargs)**2 for op in self.operators] return numpy.sqrt(sum(norm)) def direct(self, x, out=None): diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py index dc8cbf7..74ba9e4 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py @@ -49,8 +49,8 @@ class BlockScaledOperator(ScaledOperator): return self.scalar * self.operator.adjoint(x, out=out) else: raise TypeError('Operator is not linear') - def calculate_norm(self): - return numpy.abs(self.scalar) * self.operator.norm() + def norm(self, **kwargs): + return numpy.abs(self.scalar) * self.operator.norm(**kwargs) def range_geometry(self): return self.operator.range_geometry() def domain_geometry(self): diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py index 89f3549..05278ef 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -314,13 +314,8 @@ class FiniteDiff(LinearOperator): def domain_geometry(self): '''Returns the domain geometry''' return self.gm_domain - - def calculate_norm(self): - x0 = self.gm_domain.allocate('random_int') - self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0) - return self.s1 - - + + if __name__ == '__main__': from ccpi.framework import ImageGeometry diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index 475022e..8aa7f5b 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -85,12 +85,7 @@ class Gradient(LinearOperator): def range_geometry(self): return self.gm_range - - def calculate_norm(self): - x0 = self.gm_domain.allocate('random') - self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0) - return self.s1 def __rmul__(self, scalar): return ScaledOperator(self, scalar) diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py index 61428ae..f084504 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py @@ -35,7 +35,7 @@ class Identity(LinearOperator): else: out.fill(x) - def calculate_norm(self): + def calculate_norm(self, **kwargs): return 1.0 def domain_geometry(self): diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py index f304cc6..55eb692 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py @@ -29,7 +29,7 @@ class LinearOperator(Operator): # Initialise random
if x_init is None:
- x0 = operator.domain_geometry().allocate(ImageGeometry.RANDOM_INT)
+ x0 = operator.domain_geometry().allocate(type(operator.domain_geometry()).RANDOM_INT)
else:
x0 = x_init.copy()
@@ -45,23 +45,11 @@ class LinearOperator(Operator): x1.multiply((1.0/x1norm), out=x0)
return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
- @staticmethod
- def PowerMethodNonsquare(op,numiters , x_init=None):
- '''Power method to calculate iteratively the Lipschitz constant'''
-
- 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 = x1.norm()
- s[it] = (x1*x0).sum() / (x0.squared_norm())
- x0 = (1.0/x1norm)*x1
- return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
+ def calculate_norm(self, **kwargs):
+ '''Returns the norm of the LinearOperator as calculated by the PowerMethod'''
+ x0 = kwargs.get('x0', None)
+ iterations = kwargs.get('iterations', 25)
+ s1, sall, svec = LinearOperator.PowerMethod(self, iterations, x_init=x0)
+ return s1
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py index 62e22e0..8cc4c71 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py @@ -30,22 +30,7 @@ class LinearOperatorMatrix(LinearOperator): def size(self): return self.A.shape - def get_max_sing_val(self): + def calculate_norm(self, **kwargs): # If unknown, compute and store. If known, simply return it. - if self.s1 is None: - self.s1 = svds(self.A,1,return_singular_vectors=False)[0] - return self.s1 - else: - return self.s1 - def allocate_direct(self): - '''allocates the memory to hold the result of adjoint''' - #numpy.dot(self.A.transpose(),x.as_array()) - M_A, N_A = self.A.shape - out = numpy.zeros((N_A,1)) - return DataContainer(out) - def allocate_adjoint(self): - '''allocate the memory to hold the result of direct''' - #numpy.dot(self.A.transpose(),x.as_array()) - M_A, N_A = self.A.shape - out = numpy.zeros((M_A,1)) - return DataContainer(out) + return svds(self.A,1,return_singular_vectors=False)[0] + diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py index 0cac586..c1adf62 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -17,12 +17,12 @@ class Operator(object): def direct(self,x, out=None):
'''Returns the application of the Operator on x'''
raise NotImplementedError
- def norm(self):
+ def norm(self, **kwargs):
'''Returns the norm of the Operator'''
if self.__norm is None:
- self.__norm = self.calculate_norm()
+ self.__norm = self.calculate_norm(**kwargs)
return self.__norm
- def calculate_norm(self):
+ def calculate_norm(self, **kwargs):
'''Calculates the norm of the Operator'''
raise NotImplementedError
def range_geometry(self):
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py index ba0049e..f7be5de 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py @@ -40,8 +40,8 @@ class ScaledOperator(object): out *= self.scalar else: raise TypeError('Operator is not linear') - def norm(self): - return numpy.abs(self.scalar) * self.operator.norm() + def norm(self, **kwargs): + return numpy.abs(self.scalar) * self.operator.norm(**kwargs) def range_geometry(self): return self.operator.range_geometry() def domain_geometry(self): diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py index 5e318ff..eb51cc3 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py @@ -10,7 +10,7 @@ import scipy.sparse as sp import numpy as np from ccpi.framework import ImageData -class SparseFiniteDiff(): +class SparseFiniteDiff(object): def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'): diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index da73de6..1897040 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -75,42 +75,12 @@ class SymmetrizedGradient(Gradient): tmp += FiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).direct(x.get_item(j)) res.append(tmp) return res - - - -# for - -# tmp = numpy.zeros(self.gm_domain) -# -# tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) + \ -# FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) -# -# tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) + \ -# FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1]) -# -# return type(x)(tmp) - - def alloc_domain_dim(self): - return ImageData(numpy.zeros(self.gm_domain)) - - def alloc_range_dim(self): - return ImageData(numpy.zeros(self.range_dim)) - + def domain_dim(self): return self.gm_domain def range_dim(self): return self.gm_range - - def calculate_norm(self): -# return np.sqrt(4*len(self.domainDim())) - #TODO this takes time for big ImageData - # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12) - x0 = ImageData(np.random.random_sample(self.domain_dim())) - self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0) - return self.s1 - - if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py index 2331856..a2a94b7 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -29,8 +29,8 @@ class ZeroOp(Operator): else: return ImageData(np.zeros(self.gm_domain)) - def calculate_norm(self): - return 0 + def calculate_norm(self, **kwargs): + return 0. def domain_dim(self): return self.gm_domain diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index dd3238e..163d326 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -36,6 +36,7 @@ requirements: - scipy #- matplotlib - h5py + - pillow about: home: http://www.ccpi.ac.uk diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 2fc45a5..d890e46 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -121,19 +121,20 @@ class TestOperator(CCPiTestClass): uid = Id.domain_geometry().allocate(ImageGeometry.RANDOM_INT, seed=1) a = LinearOperator.PowerMethod(Id, niter, uid) - b = LinearOperator.PowerMethodNonsquare(Id, niter, uid) - + #b = LinearOperator.PowerMethodNonsquare(Id, niter, uid) + b = LinearOperator.PowerMethod(Id, niter) print ("Edo impl", a[0]) - print ("old impl", b[0]) + print ("None impl", b[0]) #self.assertAlmostEqual(a[0], b[0]) self.assertNumpyArrayAlmostEqual(a[0],b[0],decimal=6) a = LinearOperator.PowerMethod(G, niter, uid) - b = LinearOperator.PowerMethodNonsquare(G, niter, uid) + b = LinearOperator.PowerMethod(G, niter) + #b = LinearOperator.PowerMethodNonsquare(G, niter, uid) print ("Edo impl", a[0]) - print ("old impl", b[0]) + #print ("old impl", b[0]) self.assertNumpyArrayAlmostEqual(a[0],b[0],decimal=2) #self.assertAlmostEqual(a[0], b[0]) |