From c2cb31c0432dadf9f0257c6ad4a87a669a6fa6b2 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Thu, 16 Jan 2020 14:57:15 +0000
Subject: add force recalculation of norm (#487)
* add force recalculation of norm
closes #394
better handling and documentation of Operator norm method.
* added docstring and test for dot_test
---
.../ccpi/optimisation/operators/LinearOperator.py | 25 ++++++++++++---
.../Python/ccpi/optimisation/operators/Operator.py | 21 +++++++++++--
Wrappers/Python/test/test_Operator.py | 36 ++++++++++++++++++++--
3 files changed, 74 insertions(+), 8 deletions(-)
(limited to 'Wrappers/Python')
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index 5ff5b01..60a4c2a 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -72,14 +72,22 @@ class LinearOperator(Operator):
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)
+ '''Returns the norm of the LinearOperator as calculated by the PowerMethod
+
+ :param iterations: number of iterations to run
+ :type iterations: int, optional, default = 25
+ :param x_init: starting point for the iteration in the operator domain
+ :type x_init: same type as domain, a subclass of :code:`DataContainer`, optional, None
+ :parameter force: forces the recalculation of the norm
+ :type force: boolean, default :code:`False`
+ '''
+ x0 = kwargs.get('x_init', None)
iterations = kwargs.get('iterations', 25)
s1, sall, svec = LinearOperator.PowerMethod(self, iterations, x_init=x0)
return s1
@staticmethod
- def dot_test(operator, domain_init=None, range_init=None, verbose=False):
+ def dot_test(operator, domain_init=None, range_init=None, verbose=False, **kwargs):
r'''Does a dot linearity test on the operator
Evaluates if the following equivalence holds
@@ -88,10 +96,18 @@ class LinearOperator(Operator):
Ax\times y = y \times A^Tx
+ The equivalence is tested within a user specified precision
+
+ .. code::
+
+ abs(desired-actual) < 1.5 * 10**(-decimal)
+
:param operator: operator to test
:param range_init: optional initialisation container in the operator range
:param domain_init: optional initialisation container in the operator domain
:returns: boolean, True if the test is passed.
+ :param decimal: desired precision
+ :type decimal: int, optional, default 4
'''
if range_init is None:
y = operator.range_geometry().allocate('random_int')
@@ -108,8 +124,9 @@ class LinearOperator(Operator):
b = by.dot(x)
if verbose:
print ('Left hand side {}, \nRight hand side {}'.format(a, b))
+ print ("decimal ", kwargs.get('decimal', 4))
try:
- numpy.testing.assert_almost_equal(abs((a-b)/a), 0, decimal=4)
+ numpy.testing.assert_almost_equal(abs((a-b)/a), 0., decimal=kwargs.get('decimal',4))
return True
except AssertionError as ae:
print (ae)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
index 3cd0899..14d53e8 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
@@ -33,8 +33,25 @@ class Operator(object):
'''Returns the application of the Operator on x'''
raise NotImplementedError
def norm(self, **kwargs):
- '''Returns the norm of the Operator'''
- if self.__norm is None:
+ '''Returns the norm of the Operator
+
+ Calling norm triggers the calculation of the norm of the operator. Normally this
+ is a computationally expensive task, therefore we store the result of norm into
+ a member of the class. If the calculation has already run, following calls to
+ norm just return the saved member.
+ It is possible to force recalculation by setting the optional force parameter. Notice that
+ norm doesn't take notice of how many iterations or of the initialisation of the PowerMethod,
+ so in case you want to recalculate by setting a higher number of iterations or changing the
+ starting point or both you need to set :code:`force=True`
+
+ :param iterations: number of iterations to run
+ :type iterations: int, optional, default = 25
+ :param x_init: starting point for the iteration in the operator domain
+ :type x_init: same type as domain, a subclass of :code:`DataContainer`, optional, default None
+ :parameter force: forces the recalculation of the norm
+ :type force: boolean, default :code:`False`
+ '''
+ if self.__norm is None or kwargs.get('force', False):
self.__norm = self.calculate_norm(**kwargs)
return self.__norm
def calculate_norm(self, **kwargs):
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 5c659a4..57dd41f 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -191,6 +191,7 @@ class TestOperator(CCPiTestClass):
def test_Norm(self):
print ("test_BlockOperator")
##
+ numpy.random.seed(1)
N, M = 200, 300
ig = ImageGeometry(N, M)
@@ -200,8 +201,26 @@ class TestOperator(CCPiTestClass):
t1 = timer()
norm2 = G.norm()
t2 = timer()
- print ("Norm dT1 {} dT2 {}".format(t1-t0,t2-t1))
+ norm3 = G.norm(force=True)
+ t3 = timer()
+ print ("Norm dT1 {} dT2 {} dT3 {}".format(t1-t0,t2-t1, t3-t2))
self.assertLess(t2-t1, t1-t0)
+ self.assertLess(t2-t1, t3-t2)
+
+ numpy.random.seed(1)
+ t4 = timer()
+ norm4 = G.norm(iterations=50, force=True)
+ t5 = timer()
+ self.assertLess(t2-t1, t5-t4)
+
+ numpy.random.seed(1)
+ t4 = timer()
+ norm5 = G.norm(x_init=ig.allocate('random'), iterations=50, force=True)
+ t5 = timer()
+ self.assertLess(t2-t1, t5-t4)
+ for n in [norm, norm2, norm3, norm4, norm5]:
+ print ("norm {}", format(n))
+
class TestGradients(CCPiTestClass):
def setUp(self):
@@ -301,9 +320,22 @@ class TestGradients(CCPiTestClass):
lhs3 = E3.direct(u3).dot(w3)
rhs3 = u3.dot(E3.adjoint(w3))
numpy.testing.assert_almost_equal(lhs3, rhs3)
+ self.assertAlmostEqual(lhs3, rhs3)
+ print (lhs3, rhs3, abs((rhs3-lhs3)/rhs3) , 1.5 * 10**(-4), abs((rhs3-lhs3)/rhs3) < 1.5 * 10**(-4))
+ self.assertTrue( LinearOperator.dot_test(E3, range_init = w3, domain_init=u3) )
+ def test_dot_test(self):
+ Grad3 = Gradient(self.ig3, correlation = 'Space', backend='numpy')
+
# self.assertAlmostEqual(lhs3, rhs3)
- # self.assertTrue( LinearOperator.dot_test(E3) )
+ self.assertTrue( LinearOperator.dot_test(Grad3 , verbose=True))
+ self.assertTrue( LinearOperator.dot_test(Grad3 , decimal=6, verbose=True))
+ def test_dot_test2(self):
+ Grad3 = Gradient(self.ig3, correlation = 'SpaceChannel', backend='c')
+
+ # self.assertAlmostEqual(lhs3, rhs3)
+ self.assertTrue( LinearOperator.dot_test(Grad3 , verbose=True))
+ self.assertTrue( LinearOperator.dot_test(Grad3 , decimal=6, verbose=True))
--
cgit v1.2.3