From 97a45ae06f9bc95f9df640158e16f522d356392d Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 21 May 2019 15:44:59 +0100 Subject: Fix upper/lower case inconsistency to allow astra demo to run --- Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py | 4 ++-- Wrappers/Python/ccpi/optimisation/functions/__init__.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) (limited to 'Wrappers') 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 -- cgit v1.2.3 From 6794aaaa388f622cc8cd18c56c7345dca016163b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 23 May 2019 16:04:50 +0100 Subject: update pretty print --- .../ccpi/optimisation/algorithms/Algorithm.py | 44 +++++++++++++++++----- 1 file changed, 34 insertions(+), 10 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index e119a1c..0778838 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -149,10 +149,7 @@ class Algorithm(object): i = 0 if verbose: - print ("{:>9} {:>10} {:>11} {:>20}".format('Iter', - 'Max Iter', - 's/Iter', - 'Objective Value')) + print (self.verbose_header()) for _ in self: if (self.iteration -1) % self.update_objective_interval == 0: if verbose: @@ -173,12 +170,39 @@ class Algorithm(object): t = sum(timing)/len(timing) el = [ self.iteration-1, self.max_iteration, - "{:.3f} s/it".format(t), + "{:.3f}".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) - else: - out = "{:>9} {:>10} {:>11} {:>20.5e}".format(*el) + string = self.objective_to_string() + out = "{:>9} {:>10} {:>13} {}".format(*el[:-1] , string) return out + + def objective_to_string(self): + el = self.get_last_objective()[0] + 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 = 1 + 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} {:>13} {:>20}".format('Iter', + 'Max Iter', + 'Time/Iter', + 'Objective') + out += "{:>9} {:>10} {:>13} {:>20}".format('', + '', + '[s]', + '') + return out \ No newline at end of file -- cgit v1.2.3 From 31410bc217c7de4659e5df35f3423fdbbd04f59a Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 23 May 2019 16:20:08 +0100 Subject: updated print --- Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 0778838..78fe196 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -147,10 +147,9 @@ class Algorithm(object): if self.should_stop(): print ("Stop cryterion has been reached.") i = 0 - - if verbose: - print (self.verbose_header()) for _ in self: + if i == 0 and verbose: + print (self.verbose_header()) if (self.iteration -1) % self.update_objective_interval == 0: if verbose: #print ("Iteration {:>7} max: {:>7}, = {}".format(self.iteration-1, @@ -178,7 +177,7 @@ class Algorithm(object): return out def objective_to_string(self): - el = self.get_last_objective()[0] + 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]) @@ -186,7 +185,7 @@ class Algorithm(object): string = "{:>20.5e}".format(el) return string def verbose_header(self): - el = 1 + el = self.get_last_objective() if type(el) == list: out = "{:>9} {:>10} {:>13} {:>13} {:>13} {:>15}\n".format('Iter', 'Max Iter', @@ -197,7 +196,7 @@ class Algorithm(object): '[s]', 'Objective' , 'Objective', 'Gap') else: - out = "{:>9} {:>10} {:>13} {:>20}".format('Iter', + out = "{:>9} {:>10} {:>13} {:>20}\n".format('Iter', 'Max Iter', 'Time/Iter', 'Objective') -- cgit v1.2.3 From 8787adfef38ded4852b53a6d7812bb7d4f888d86 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 3 Jun 2019 11:18:23 -0400 Subject: fix import and calls in tests --- Wrappers/Python/test/test_functions.py | 6 +++--- Wrappers/Python/test/test_run_test.py | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py index af419c7..523008d 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(A, L2NormSquared(b=b)) 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..a6c13f4 100755 --- a/Wrappers/Python/test/test_run_test.py +++ b/Wrappers/Python/test/test_run_test.py @@ -8,7 +8,7 @@ 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 @@ -80,7 +80,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,7 +144,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 @@ -218,7 +218,7 @@ 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() @@ -286,7 +286,7 @@ 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] -- cgit v1.2.3 From d688da2c586e3165bc366772668ecb008d03eca3 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 4 Jun 2019 04:57:46 -0400 Subject: fix print for python27 --- Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index 78fe196..9f54c1e 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -167,13 +167,12 @@ class Algorithm(object): t = 0 else: t = sum(timing)/len(timing) - el = [ self.iteration-1, - self.max_iteration, - "{:.3f}".format(t), - self.get_last_objective() ] - - string = self.objective_to_string() - out = "{:>9} {:>10} {:>13} {}".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): @@ -204,4 +203,4 @@ class Algorithm(object): '', '[s]', '') - return out \ No newline at end of file + return out -- cgit v1.2.3 From 1580b6fbaebc77649c2eebf89ea7a9ebda8ba91c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 6 Jun 2019 12:16:53 +0000 Subject: Norm refactor (#299) * Fix upper/lower case inconsistency to allow astra demo to run * fix import and calls in tests * add storage in Operator Class of the operator norm scaled operator is not an Operator * add test for norm * reduced iteration in PowerMethod --- .../ccpi/optimisation/operators/BlockOperator.py | 2 +- .../optimisation/operators/BlockScaledOperator.py | 2 +- .../operators/FiniteDifferenceOperator.py | 2 +- .../optimisation/operators/GradientOperator.py | 2 +- .../optimisation/operators/IdentityOperator.py | 2 +- .../Python/ccpi/optimisation/operators/Operator.py | 8 +++++++ .../operators/SymmetrizedGradientOperator.py | 2 +- .../ccpi/optimisation/operators/ZeroOperator.py | 2 +- Wrappers/Python/test/test_Operator.py | 26 +++++++++++++++++----- 9 files changed, 36 insertions(+), 12 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index 1d77510..10fe34c 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -104,7 +104,7 @@ class BlockOperator(Operator): index = row*self.shape[1]+col return self.operators[index] - def norm(self): + def calculate_norm(self): norm = [op.norm()**2 for op in self.operators] return numpy.sqrt(sum(norm)) diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py index aeb6c53..dc8cbf7 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py @@ -49,7 +49,7 @@ class BlockScaledOperator(ScaledOperator): return self.scalar * self.operator.adjoint(x, out=out) else: raise TypeError('Operator is not linear') - def norm(self): + def calculate_norm(self): return numpy.abs(self.scalar) * self.operator.norm() def range_geometry(self): return self.operator.range_geometry() diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py index 2c42532..89f3549 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -315,7 +315,7 @@ class FiniteDiff(LinearOperator): '''Returns the domain geometry''' return self.gm_domain - def norm(self): + def calculate_norm(self): x0 = self.gm_domain.allocate('random_int') self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0) return self.s1 diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index e0b8a32..475022e 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -86,7 +86,7 @@ class Gradient(LinearOperator): def range_geometry(self): return self.gm_range - def norm(self): + def calculate_norm(self): x0 = self.gm_domain.allocate('random') self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0) diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py index a58a296..61428ae 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): return 1.0 def domain_geometry(self): diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py index 2d2089b..0cac586 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -8,6 +8,9 @@ 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 @@ -16,6 +19,11 @@ class Operator(object): raise NotImplementedError def norm(self): '''Returns the norm of the Operator''' + if self.__norm is None: + self.__norm = self.calculate_norm() + return self.__norm + def calculate_norm(self): + '''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/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index 3fc43b8..da73de6 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -102,7 +102,7 @@ class SymmetrizedGradient(Gradient): def range_dim(self): return self.gm_range - def norm(self): + 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) diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py index a7c5f09..2331856 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -29,7 +29,7 @@ class ZeroOp(Operator): else: return ImageData(np.zeros(self.gm_domain)) - def norm(self): + def calculate_norm(self): return 0 def domain_dim(self): diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 5c97545..2fc45a5 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) @@ -134,7 +137,20 @@ class TestOperator(CCPiTestClass): 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 +191,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 +374,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 +443,7 @@ class TestBlockOperator(unittest.TestCase): self.assertGreater(t1,t2) def test_BlockOperatorLinearValidity(self): - + print ("test_BlockOperatorLinearValidity") M, N = 3, 4 ig = ImageGeometry(M, N) -- cgit v1.2.3 From 68d716ff14264bfaac8c9311d435530689382e7e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 7 Jun 2019 13:25:15 +0000 Subject: Merge norm ref (#301) * Fix upper/lower case inconsistency to allow astra demo to run * fix import and calls in tests * add storage in Operator Class of the operator norm scaled operator is not an Operator * add test for norm * simplified calls to norm for LinearOperators, added default in base class * add pillow as dependency * resolve conflict --- .../ccpi/optimisation/operators/BlockOperator.py | 4 +-- .../optimisation/operators/BlockScaledOperator.py | 4 +-- .../operators/FiniteDifferenceOperator.py | 9 ++---- .../optimisation/operators/GradientOperator.py | 5 ---- .../optimisation/operators/IdentityOperator.py | 2 +- .../ccpi/optimisation/operators/LinearOperator.py | 26 +++++------------- .../optimisation/operators/LinearOperatorMatrix.py | 21 ++------------ .../Python/ccpi/optimisation/operators/Operator.py | 6 ++-- .../ccpi/optimisation/operators/ScaledOperator.py | 4 +-- .../optimisation/operators/SparseFiniteDiff.py | 2 +- .../operators/SymmetrizedGradientOperator.py | 32 +--------------------- .../ccpi/optimisation/operators/ZeroOperator.py | 4 +-- Wrappers/Python/conda-recipe/meta.yaml | 1 + Wrappers/Python/test/test_Operator.py | 11 ++++---- 14 files changed, 33 insertions(+), 98 deletions(-) (limited to 'Wrappers') 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]) -- cgit v1.2.3