diff options
26 files changed, 156 insertions, 224 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 0ba2c55..c62d0ea 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -51,6 +51,7 @@ class Algorithm(object): self.__max_iteration = kwargs.get('max_iteration', 0) self.__loss = [] self.memopt = False + self.configured = False self.timing = [] self.update_objective_interval = kwargs.get('update_objective_interval', 1) def set_up(self, *args, **kwargs): @@ -86,6 +87,8 @@ class Algorithm(object): raise StopIteration() else: time0 = time.time() + if not self.configured: + raise ValueError('Algorithm not configured correctly. Please run set_up.') self.update() self.timing.append( time.time() - time0 ) if self.iteration % self.update_objective_interval == 0: @@ -147,13 +150,9 @@ class Algorithm(object): if self.should_stop(): print ("Stop cryterion has been reached.") i = 0 - - if verbose: - print ("{:>9} {:>10} {:>11} {:>20}".format('Iter', - 'Max Iter', - 's/Iter', - 'Objective Value')) for _ in self: + if i == 0 and verbose: + print (self.verbose_header()) if (self.iteration -1) % self.update_objective_interval == 0: if verbose: print (self.verbose_output()) @@ -170,14 +169,40 @@ class Algorithm(object): t = 0 else: t = sum(timing)/len(timing) - el = [ self.iteration-1, - self.max_iteration, - "{:.3f} s/it".format(t), - self.get_last_objective() ] - - if type(el[-1] ) == list: - string = functools.reduce(lambda x,y: x+' {:>15.5e}'.format(y), el[-1],'') - out = "{:>9} {:>10} {:>11} {}".format(*el[:-1] , string) + out = "{:>9} {:>10} {:>13} {}".format( + self.iteration-1, + self.max_iteration, + "{:.3f}".format(t), + self.objective_to_string() + ) + return out + + def objective_to_string(self): + el = self.get_last_objective() + if type(el) == list: + string = functools.reduce(lambda x,y: x+' {:>13.5e}'.format(y), el[:-1],'') + string += '{:>15.5e}'.format(el[-1]) + else: + string = "{:>20.5e}".format(el) + return string + def verbose_header(self): + el = self.get_last_objective() + if type(el) == list: + out = "{:>9} {:>10} {:>13} {:>13} {:>13} {:>15}\n".format('Iter', + 'Max Iter', + 'Time/Iter', + 'Primal' , 'Dual', 'Primal-Dual') + out += "{:>9} {:>10} {:>13} {:>13} {:>13} {:>15}".format('', + '', + '[s]', + 'Objective' , 'Objective', 'Gap') else: - out = "{:>9} {:>10} {:>11} {:>20.5e}".format(*el) + out = "{:>9} {:>10} {:>13} {:>20}\n".format('Iter', + 'Max Iter', + 'Time/Iter', + 'Objective') + out += "{:>9} {:>10} {:>13} {:>20}".format('', + '', + '[s]', + '') return out diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py index 4d4843c..6b610a0 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -23,6 +23,7 @@ Created on Thu Feb 21 11:11:23 2019 """ from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.functions import Norm2Sq class CGLS(Algorithm): @@ -59,6 +60,9 @@ class CGLS(Algorithm): # self.normr2 = sum(self.normr2) #self.normr2 = numpy.sqrt(self.normr2) #print ("set_up" , self.normr2) + n = Norm2Sq(operator, self.data) + self.loss.append(n(x_init)) + self.configured = True def update(self): @@ -84,4 +88,4 @@ class CGLS(Algorithm): self.d = s + beta*self.d def update_objective(self): - self.loss.append(self.r.squared_norm())
\ No newline at end of file + self.loss.append(self.r.squared_norm()) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py deleted file mode 100644 index aa07359..0000000 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py +++ /dev/null @@ -1,86 +0,0 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2019 Edoardo Pasca - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -""" -Created on Thu Feb 21 11:09:03 2019 - -@author: ofn77899 -""" - -from ccpi.optimisation.algorithms import Algorithm -from ccpi.optimisation.functions import ZeroFunction - -class FBPD(Algorithm): - '''FBPD Algorithm - - Parameters: - x_init: initial guess - f: constraint - g: data fidelity - h: regularizer - opt: additional algorithm - ''' - constraint = None - data_fidelity = None - regulariser = None - def __init__(self, **kwargs): - pass - def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\ - regulariser=None, opt=None): - - # default inputs - if constraint is None: - self.constraint = ZeroFun() - else: - self.constraint = constraint - if data_fidelity is None: - data_fidelity = ZeroFun() - else: - self.data_fidelity = data_fidelity - if regulariser is None: - self.regulariser = ZeroFun() - else: - self.regulariser = regulariser - - # algorithmic parameters - - - # step-sizes - self.tau = 2 / (self.data_fidelity.L + 2) - self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L - - self.inv_sigma = 1/self.sigma - - # initialization - self.x = x_init - self.y = operator.direct(self.x) - - - def update(self): - - # primal forward-backward step - x_old = self.x - self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) ) - self.x = self.constraint.prox(self.x, self.tau); - - # dual forward-backward step - self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old); - self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma); - - # time and criterion - self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x)) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py index 419958a..c8fd0d4 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -29,12 +29,13 @@ class FISTA(Algorithm): proper variables are passed or later with set_up''' super(FISTA, self).__init__() self.f = kwargs.get('f', None) - self.g = kwargs.get('g', None) + self.g = kwargs.get('g', ZeroFunction()) self.x_init = kwargs.get('x_init',None) self.invL = None self.t_old = 1 - if self.f is not None and self.g is not None: - print ("Calling from creator") + if self.x_init is not None and \ + self.f is not None and self.g is not None: + print ("FISTA set_up called from creator") self.set_up(self.x_init, self.f, self.g) @@ -56,6 +57,8 @@ class FISTA(Algorithm): self.invL = 1/f.L self.t_old = 1 + self.update_objective() + self.configured = True def update(self): diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py index 14763c5..34bf954 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py @@ -40,7 +40,7 @@ class GradientDescent(Algorithm): if k in args: args.pop(args.index(k)) if len(args) == 0: - return self.set_up(x_init=kwargs['x_init'], + self.set_up(x_init=kwargs['x_init'], objective_function=kwargs['objective_function'], rate=kwargs['rate']) @@ -61,6 +61,7 @@ class GradientDescent(Algorithm): self.memopt = False if self.memopt: self.x_update = x_init.copy() + self.configured = True def update(self): '''Single iteration''' diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 39b092b..3afd8b0 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -23,10 +23,16 @@ class PDHG(Algorithm): self.operator = kwargs.get('operator', None) self.g = kwargs.get('g', None) self.tau = kwargs.get('tau', None) - self.sigma = kwargs.get('sigma', None) + self.sigma = kwargs.get('sigma', 1.) + if self.f is not None and self.operator is not None and \ self.g is not None: + if self.tau is None: + # Compute operator Norm + normK = self.operator.norm() + # Primal & dual stepsizes + self.tau = 1/(self.sigma*normK**2) print ("Calling from creator") self.set_up(self.f, self.g, @@ -57,6 +63,8 @@ class PDHG(Algorithm): # relaxation parameter self.theta = 1 + self.update_objective() + self.configured = True def update(self): diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py index 30584d4..c73d323 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py @@ -59,6 +59,7 @@ class SIRT(Algorithm): # Set up scaling matrices D and M. self.M = 1/self.operator.direct(self.operator.domain_geometry().allocate(value=1.0)) self.D = 1/self.operator.adjoint(self.operator.range_geometry().allocate(value=1.0)) + self.configured = True def update(self): @@ -67,8 +68,9 @@ class SIRT(Algorithm): self.x += self.relax_par * (self.D*self.operator.adjoint(self.M*self.r)) - if self.constraint != None: - self.x = self.constraint.prox(self.x,None) + if self.constraint is not None: + self.x = self.constraint.proximal(self.x,None) + # self.constraint.proximal(self.x,None, out=self.x) def update_objective(self): - self.loss.append(self.r.squared_norm())
\ No newline at end of file + self.loss.append(self.r.squared_norm()) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py index 2dbacfc..8f255f3 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -27,7 +27,5 @@ from .CGLS import CGLS from .SIRT import SIRT from .GradientDescent import GradientDescent from .FISTA import FISTA -from .FBPD import FBPD from .PDHG import PDHG -from .PDHG import PDHG_old diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py index b553d7c..0b6bb0b 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py +++ b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py @@ -21,7 +21,7 @@ import numpy import warnings # Define a class for squared 2-norm -class Norm2sq(Function): +class Norm2Sq(Function): ''' f(x) = c*||A*x-b||_2^2 @@ -38,7 +38,7 @@ class Norm2sq(Function): ''' def __init__(self,A,b,c=1.0,memopt=False): - super(Norm2sq, self).__init__() + super(Norm2Sq, self).__init__() self.A = A # Should be an operator, default identity self.b = b # Default zero DataSet? diff --git a/Wrappers/Python/ccpi/optimisation/functions/__init__.py b/Wrappers/Python/ccpi/optimisation/functions/__init__.py index a82ee3e..c0eab31 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py @@ -10,4 +10,4 @@ from .FunctionOperatorComposition import FunctionOperatorComposition from .MixedL21Norm import MixedL21Norm from .IndicatorBox import IndicatorBox from .KullbackLeibler import KullbackLeibler -from .Norm2Sq import Norm2sq +from .Norm2Sq import Norm2Sq diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index c8bd546..5f04363 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 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 aeb6c53..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 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 f459334..876f45f 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 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 d98961b..cd58b7d 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -109,13 +109,6 @@ class Gradient(LinearOperator): def range_geometry(self): return self.gm_range - - def norm(self, **kwargs): - - x0 = self.gm_domain.allocate('random') - iterations = kwargs.get('iterations', 10) - self.s1, sall, svec = LinearOperator.PowerMethod(self, iterations, 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 a853b8d..8f35373 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 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 2d2089b..c1adf62 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -8,14 +8,22 @@ from ccpi.optimisation.operators.ScaledOperator import ScaledOperator class Operator(object):
'''Operator that maps from a space X -> Y'''
+ def __init__(self, **kwargs):
+ self.__norm = None
+
def is_linear(self):
'''Returns if the operator is linear'''
return False
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(**kwargs)
+ return self.__norm
+ def calculate_norm(self, **kwargs):
+ '''Calculates the norm of the Operator'''
raise NotImplementedError
def range_geometry(self):
'''Returns the range of the Operator: Y space'''
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 c5c2ec9..cb76dce 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 243809b..205f7e1 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -111,17 +111,6 @@ class SymmetrizedGradient(Gradient): def range_geometry(self): return self.gm_range - - def norm(self): - -# TODO need dot method for BlockDataContainer -# return numpy.sqrt(4*self.gm_domain.shape[0]) - -# x0 = self.gm_domain.allocate('random') - self.s1, sall, svec = LinearOperator.PowerMethod(self, 50) - return self.s1 - - if __name__ == '__main__': @@ -241,4 +230,4 @@ if __name__ == '__main__': # # # -#
\ No newline at end of file +# diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py index 8168f0b..67808de 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -34,8 +34,8 @@ class ZeroOperator(LinearOperator): else: out.fill(self.gm_domain.allocate()) - def norm(self): - return 0 + def calculate_norm(self, **kwargs): + return 0. def domain_geometry(self): return self.gm_domain diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 5c97545..d890e46 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -51,6 +51,7 @@ class CCPiTestClass(unittest.TestCase): class TestOperator(CCPiTestClass): def test_ScaledOperator(self): + print ("test_ScaledOperator") ig = ImageGeometry(10,20,30) img = ig.allocate() scalar = 0.5 @@ -58,7 +59,8 @@ class TestOperator(CCPiTestClass): numpy.testing.assert_array_equal(scalar * img.as_array(), sid.direct(img).as_array()) - def test_TomoIdentity(self): + def test_Identity(self): + print ("test_Identity") ig = ImageGeometry(10,20,30) img = ig.allocate() self.assertTrue(img.shape == (30,20,10)) @@ -107,9 +109,10 @@ class TestOperator(CCPiTestClass): self.assertBlockDataContainerEqual(res, w) def test_PowerMethod(self): + print ("test_BlockOperator") N, M = 200, 300 - niter = 1000 + niter = 10 ig = ImageGeometry(N, M) Id = Identity(ig) @@ -118,23 +121,37 @@ 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]) + def test_Norm(self): + print ("test_BlockOperator") + ## + N, M = 200, 300 + ig = ImageGeometry(N, M) + G = Gradient(ig) + t0 = timer() + norm = G.norm() + t1 = timer() + norm2 = G.norm() + t2 = timer() + print ("Norm dT1 {} dT2 {}".format(t1-t0,t2-t1)) + self.assertLess(t2-t1, t1-t0) @@ -175,7 +192,7 @@ class TestBlockOperator(unittest.TestCase): self.assertTrue(res) def test_BlockOperator(self): - + print ("test_BlockOperator") M, N = 3, 4 ig = ImageGeometry(M, N) @@ -358,7 +375,7 @@ class TestBlockOperator(unittest.TestCase): print("Z1", Z1[0][1].as_array()) print("RES1", RES1[0][1].as_array()) def test_timedifference(self): - + print ("test_timedifference") M, N ,W = 100, 512, 512 ig = ImageGeometry(M, N, W) arr = ig.allocate('random_int') @@ -427,7 +444,7 @@ class TestBlockOperator(unittest.TestCase): self.assertGreater(t1,t2) def test_BlockOperatorLinearValidity(self): - + print ("test_BlockOperatorLinearValidity") M, N = 3, 4 ig = ImageGeometry(M, N) diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 669804e..8c398f4 100755 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -13,11 +13,11 @@ from ccpi.framework import AcquisitionData from ccpi.framework import ImageGeometry from ccpi.framework import AcquisitionGeometry from ccpi.optimisation.operators import Identity -from ccpi.optimisation.funcs import Norm2sq +from ccpi.optimisation.functions import Norm2Sq, ZeroFunction, \ + L2NormSquared, FunctionOperatorComposition from ccpi.optimisation.algorithms import GradientDescent from ccpi.optimisation.algorithms import CGLS from ccpi.optimisation.algorithms import FISTA -from ccpi.optimisation.algorithms import FBPD @@ -50,11 +50,13 @@ class TestAlgorithms(unittest.TestCase): identity = Identity(ig) - norm2sq = Norm2sq(identity, b) + norm2sq = Norm2Sq(identity, b) + rate = 0.3 + rate = norm2sq.L / 3. alg = GradientDescent(x_init=x_init, objective_function=norm2sq, - rate=0.3) + rate=rate) alg.max_iteration = 20 alg.run(20, verbose=True) self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) @@ -80,20 +82,19 @@ class TestAlgorithms(unittest.TestCase): b = x_init.copy() # fill with random numbers b.fill(numpy.random.random(x_init.shape)) - x_init = ImageData(geometry=ig) - x_init.fill(numpy.random.random(x_init.shape)) - + x_init = ig.allocate(ImageGeometry.RANDOM) identity = Identity(ig) - norm2sq = Norm2sq(identity, b) - norm2sq.L = 2 * norm2sq.c * identity.norm()**2 + #### it seems FISTA does not work with Nowm2Sq + # norm2sq = Norm2Sq(identity, b) + # norm2sq.L = 2 * norm2sq.c * identity.norm()**2 + norm2sq = FunctionOperatorComposition(L2NormSquared(b=b), identity) opt = {'tol': 1e-4, 'memopt':False} - alg = FISTA(x_init=x_init, f=norm2sq, g=None, opt=opt) + print ("initial objective", norm2sq(x_init)) + alg = FISTA(x_init=x_init, f=norm2sq, g=ZeroFunction()) alg.max_iteration = 2 alg.run(20, verbose=True) self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) - alg.run(20, verbose=True) - self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index 082548b..021ad99 100644 --- a/Wrappers/Python/test/test_functions.py +++ b/Wrappers/Python/test/test_functions.py @@ -20,7 +20,7 @@ from ccpi.optimisation.operators import Gradient from ccpi.optimisation.functions import L2NormSquared from ccpi.optimisation.functions import L1Norm, MixedL21Norm -from ccpi.optimisation.functions import Norm2sq +from ccpi.optimisation.functions import Norm2Sq from ccpi.optimisation.functions import ZeroFunction from ccpi.optimisation.functions import FunctionOperatorComposition @@ -298,7 +298,7 @@ class TestFunction(unittest.TestCase): b = ig.allocate(ImageGeometry.RANDOM_INT) A = 0.5 * Identity(ig) - old_chisq = Norm2sq(A, b, 1.0) + old_chisq = Norm2Sq(A, b, 1.0) new_chisq = FunctionOperatorComposition(L2NormSquared(b=b),A) yold = old_chisq(u) @@ -364,4 +364,4 @@ class TestFunction(unittest.TestCase): proxc = f.proximal_conjugate(x,1.2) f.proximal_conjugate(x, 1.2, out=out) - numpy.testing.assert_array_equal(proxc.as_array(), out.as_array())
\ No newline at end of file + numpy.testing.assert_array_equal(proxc.as_array(), out.as_array()) diff --git a/Wrappers/Python/test/test_run_test.py b/Wrappers/Python/test/test_run_test.py index 9b4d53b..22d4a78 100755 --- a/Wrappers/Python/test/test_run_test.py +++ b/Wrappers/Python/test/test_run_test.py @@ -8,10 +8,11 @@ from ccpi.framework import ImageGeometry from ccpi.framework import AcquisitionGeometry from ccpi.optimisation.algorithms import FISTA #from ccpi.optimisation.algs import FBPD -from ccpi.optimisation.functions import Norm2sq +from ccpi.optimisation.functions import Norm2Sq from ccpi.optimisation.functions import ZeroFunction -from ccpi.optimisation.funcs import Norm1 -from ccpi.optimisation.funcs import Norm2 +from ccpi.optimisation.functions import L1Norm +# This was removed +#from ccpi.optimisation.funcs import Norm2 from ccpi.optimisation.operators import LinearOperatorMatrix from ccpi.optimisation.operators import Identity @@ -80,7 +81,7 @@ class TestAlgorithms(unittest.TestCase): lam = 10 opt = {'memopt': True} # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) + f = Norm2Sq(A, b, c=0.5, memopt=True) g0 = ZeroFunction() # Initial guess @@ -144,14 +145,14 @@ class TestAlgorithms(unittest.TestCase): lam = 10 opt = {'memopt': True} # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) + f = Norm2Sq(A, b, c=0.5, memopt=True) g0 = ZeroFunction() # Initial guess x_init = DataContainer(np.zeros((n, 1))) # Create 1-norm object instance - g1 = Norm1(lam) + g1 = lam * L1Norm() g1(x_init) g1.prox(x_init, 0.02) @@ -218,14 +219,14 @@ class TestAlgorithms(unittest.TestCase): x_init = DataContainer(np.random.randn(n, 1)) # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) + f = Norm2Sq(A, b, c=0.5, memopt=True) f.L = LinearOperator.PowerMethod(A, 25, x_init)[0] print ("Lipschitz", f.L) g0 = ZeroFun() # Create 1-norm object instance - g1 = Norm1(lam) + g1 = lam * L1Norm() # Compare to CVXPY @@ -286,13 +287,13 @@ class TestAlgorithms(unittest.TestCase): y.array = y.array + 0.1*np.random.randn(N, N) # Data fidelity term - f_denoise = Norm2sq(I, y, c=0.5, memopt=True) + f_denoise = Norm2Sq(I, y, c=0.5, memopt=True) x_init = ImageData(geometry=ig) f_denoise.L = LinearOperator.PowerMethod(I, 25, x_init)[0] # 1-norm regulariser lam1_denoise = 1.0 - g1_denoise = Norm1(lam1_denoise) + g1_denoise = lam1_denoise * L1Norm() # Initial guess x_init_denoise = ImageData(np.zeros((N, N))) |