From 5b6e219e78c0c006452ab4a0059d14951f83cba9 Mon Sep 17 00:00:00 2001
From: Margaret Duff <margaret.duff@stfc.ac.uk>
Date: Tue, 13 Aug 2024 15:53:06 +0000
Subject: [PATCH 1/6] Added documentation for CGLS including tolerance
 documentation

---
 .../cil/optimisation/algorithms/CGLS.py       | 64 +++++++++++--------
 1 file changed, 39 insertions(+), 25 deletions(-)

diff --git a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
index fe6fab863..64386c5fd 100644
--- a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
@@ -25,33 +25,41 @@
 
 class CGLS(Algorithm):
 
-    r'''Conjugate Gradient Least Squares algorithm
+    r'''Conjugate Gradient Least Squares (CGLS) algorithm
+    
+    The Conjugate Gradient Least Squares (CGLS) algorithm is commonly used for solving large systems of linear equations, due to its fast convergence.
 
     Problem:
 
     .. math::
 
-      \min || A x - b ||^2_2
-
-    |
-
-    Parameters :
-
-      :parameter operator : Linear operator for the inverse problem
-      :parameter initial : Initial guess ( Default initial = 0)
-      :parameter data : Acquired data to reconstruct
-      :parameter tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm
-
-    Reference:
-        https://web.stanford.edu/group/SOL/software/cgls/
+      \min_x || A x - b ||^2_2
+      
+      
+    Note
+    ----
+    By default, this algorithm will terminate if the value of :math:`||A^T(Ax-b)||_2 < tol*||A^T(Ax_0-b)||_2` where 'tol' is set to default as '1e-6', :math:`x` is the current iterate and :math:`x_0` is the initial value. 
+    It will also terminate if the algorithm begins to diverge i.e. if :math:`||x||_2>1/{tol}`. 
+    
+    By setting 'tolerance' to be '0' you can prevent the algorithm stopping on either of these criteria. 
+
+    Parameters
+    ------------
+    operator : Operator
+        Linear operator for the inverse problem
+    initial : (optional) DataContainer in the domain of the operator, default is a DataContainer filled with zeros. 
+        Initial guess 
+    data : Data Container in the range of the operator 
+        Acquired data to reconstruct
+    tolerance: float, default 1e-6 
+        Tolerance/ Stopping Criterion to end CGLS algorithm.  
+
+    Reference
+    ---------
+    https://web.stanford.edu/group/SOL/software/cgls/
     '''
     def __init__(self, initial=None, operator=None, data=None, tolerance=1e-6, **kwargs):
         '''initialisation of the algorithm
-
-        :param operator : Linear operator for the inverse problem
-        :param initial : Initial guess ( Default initial = 0)
-        :param data : Acquired data to reconstruct
-        :param tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm
         '''
         super(CGLS, self).__init__(**kwargs)
 
@@ -61,13 +69,19 @@ def __init__(self, initial=None, operator=None, data=None, tolerance=1e-6, **kwa
             self.set_up(initial=initial, operator=operator, data=data, tolerance=tolerance)
 
     def set_up(self, initial, operator, data, tolerance=1e-6):
-        '''initialisation of the algorithm
-
-        :param operator: Linear operator for the inverse problem
-        :param initial: Initial guess ( Default initial = 0)
-        :param data: Acquired data to reconstruct
-        :param tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm
+        r'''initialisation of the algorithm
+        Parameters
+        ------------
+        operator : Operator
+            Linear operator for the inverse problem
+        initial : (optional) DataContainer in the domain of the operator, default is a DataContainer filled with zeros. 
+            Initial guess 
+        data : Data Container in the range of the operator 
+            Acquired data to reconstruct
+        tolerance: float, default 1e-6 
+            Tolerance/ Stopping Criterion to end the CGLS algorithm
         '''
+        
         log.info("%s setting up", self.__class__.__name__)
         self.x = initial.copy()
         self.operator = operator

From f902afc9bb0428d43c1d21286fc14f0b6249d324 Mon Sep 17 00:00:00 2001
From: Margaret Duff <margaret.duff@stfc.ac.uk>
Date: Thu, 15 Aug 2024 11:16:18 +0000
Subject: [PATCH 2/6] New unit tests for CGLS, deprecated tolerance and two new
 callbacks

---
 .../cil/optimisation/algorithms/CGLS.py       |  37 +++--
 .../cil/optimisation/utilities/callbacks.py   |  43 +++++
 Wrappers/Python/test/test_algorithms.py       | 155 +++++++++++++++---
 3 files changed, 195 insertions(+), 40 deletions(-)

diff --git a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
index 64386c5fd..cec86b61b 100644
--- a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
@@ -19,6 +19,7 @@
 from cil.optimisation.algorithms import Algorithm
 import numpy
 import logging
+import warnings 
 
 log = logging.getLogger(__name__)
 
@@ -36,13 +37,6 @@ class CGLS(Algorithm):
       \min_x || A x - b ||^2_2
       
       
-    Note
-    ----
-    By default, this algorithm will terminate if the value of :math:`||A^T(Ax-b)||_2 < tol*||A^T(Ax_0-b)||_2` where 'tol' is set to default as '1e-6', :math:`x` is the current iterate and :math:`x_0` is the initial value. 
-    It will also terminate if the algorithm begins to diverge i.e. if :math:`||x||_2>1/{tol}`. 
-    
-    By setting 'tolerance' to be '0' you can prevent the algorithm stopping on either of these criteria. 
-
     Parameters
     ------------
     operator : Operator
@@ -51,24 +45,33 @@ class CGLS(Algorithm):
         Initial guess 
     data : Data Container in the range of the operator 
         Acquired data to reconstruct
-    tolerance: float, default 1e-6 
-        Tolerance/ Stopping Criterion to end CGLS algorithm.  
+        
+    Note
+    -----
+    Passing tolerance directly to CGLS is being deprecated. Instead we recommend using the callback functionality: https://tomographicimaging.github.io/CIL/nightly/optimisation/#callbacks and in particular the CGLSEarlyStopping callback replicated the old behaviour.
 
     Reference
     ---------
     https://web.stanford.edu/group/SOL/software/cgls/
     '''
-    def __init__(self, initial=None, operator=None, data=None, tolerance=1e-6, **kwargs):
+    def __init__(self, initial=None, operator=None, data=None,  **kwargs):
         '''initialisation of the algorithm
         '''
+        #We are deprecating tolerance 
+        self.tolerance=kwargs.pop("tolerance", None)
+        if self.tolerance is not None:
+            warnings.warn( stacklevel=2, category=DeprecationWarning, message="Passing tolerance directly to CGLS is being deprecated. Instead we recommend using the callback functionality: https://tomographicimaging.github.io/CIL/nightly/optimisation/#callbacks and in particular the CGLSEarlyStopping callback replicated the old behaviour")
+        else:
+            self.tolerance = 0
+        
         super(CGLS, self).__init__(**kwargs)
 
         if initial is None and operator is not None:
             initial = operator.domain_geometry().allocate(0)
         if initial is not None and operator is not None and data is not None:
-            self.set_up(initial=initial, operator=operator, data=data, tolerance=tolerance)
+            self.set_up(initial=initial, operator=operator, data=data) 
 
-    def set_up(self, initial, operator, data, tolerance=1e-6):
+    def set_up(self, initial, operator, data):
         r'''initialisation of the algorithm
         Parameters
         ------------
@@ -78,14 +81,12 @@ def set_up(self, initial, operator, data, tolerance=1e-6):
             Initial guess 
         data : Data Container in the range of the operator 
             Acquired data to reconstruct
-        tolerance: float, default 1e-6 
-            Tolerance/ Stopping Criterion to end the CGLS algorithm
+
         '''
         
         log.info("%s setting up", self.__class__.__name__)
         self.x = initial.copy()
         self.operator = operator
-        self.tolerance = tolerance
 
         self.r = data - self.operator.direct(self.x)
         self.s = self.operator.adjoint(self.r)
@@ -124,7 +125,7 @@ def update(self):
         #self.p = self.s + self.beta * self.p
         self.p.sapyb(self.beta, self.s, 1, out=self.p)
 
-        self.normx = self.x.norm()
+        self.normx = self.x.norm()#Can be deprecated with tolerance 
 
 
     def update_objective(self):
@@ -133,10 +134,10 @@ def update_objective(self):
             raise StopIteration()
         self.loss.append(a)
 
-    def should_stop(self):
+    def should_stop(self): #can be deprecated with tolerance 
         return self.flag() or super().should_stop()
 
-    def flag(self):
+    def flag(self): #can be deprecated with tolerance 
         '''returns whether the tolerance has been reached'''
         flag  = (self.norms <= self.norms0 * self.tolerance) or (self.normx * self.tolerance >= 1)
 
diff --git a/Wrappers/Python/cil/optimisation/utilities/callbacks.py b/Wrappers/Python/cil/optimisation/utilities/callbacks.py
index 2c357d9a8..b5a79a7f5 100644
--- a/Wrappers/Python/cil/optimisation/utilities/callbacks.py
+++ b/Wrappers/Python/cil/optimisation/utilities/callbacks.py
@@ -3,6 +3,7 @@
 
 from tqdm.auto import tqdm as tqdm_auto
 from tqdm.std import tqdm as tqdm_std
+import numpy as np
 
 
 class Callback(ABC):
@@ -134,3 +135,45 @@ class LogfileCallback(TextProgressCallback):
     def __init__(self, log_file, mode='a', **kwargs):
         self.fd = open(log_file, mode=mode)
         super().__init__(file=self.fd, **kwargs)
+        
+class EarlyStoppingObjectiveValue(Callback):
+    '''Callback that stops iterations if the change in the objective value is less than a provided threshold value.
+
+    Parameters
+    ----------
+    threshold: float, default 1e-6 
+
+    Note
+    -----
+    This callback only compares the last two calculated objective values. If `update_objective_interval` is greater than 1, the objective value is not calculated at each iteration. 
+    
+        '''
+    def __init__(self, threshold=1e-6):
+        self.threshold=threshold
+    
+    
+    def __call__(self, algorithm):
+        if len(algorithm.loss)>=2:
+            if np.abs(algorithm.loss[-1]-algorithm.loss[-2])<self.threshold:
+                raise StopIteration
+                
+class CGLSEarlyStopping(Callback):
+    '''Callback to work with CGLS. It causes the algorithm to terminate if the value of :math:`||A^T(Ax-b)||_2 < tol*||A^T(Ax_0-b)||_2` where 'tol' is set to default as '1e-6', :math:`x` is the current iterate and :math:`x_0` is the initial value. 
+    It will also terminate if the algorithm begins to diverge i.e. if :math:`||x||_2>1/{tol}`. 
+    Parameters
+    ----------
+    tolerance: float, default 1e-6 
+ 
+    '''
+    def __init__(self, tolerance=1e-6):
+        self.tolerance=tolerance
+    
+    
+    def __call__(self, algorithm):
+        self.normx = algorithm.x.norm()
+        
+        if (algorithm.norms <= algorithm.norms0 * self.tolerance) or (algorithm.normx * self.tolerance >= 1):
+            print('Tolerance is reached: {}'.format(self.tolerance))
+            raise StopIteration
+            
+        
diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py
index 7ae29e170..888a59d7e 100644
--- a/Wrappers/Python/test/test_algorithms.py
+++ b/Wrappers/Python/test/test_algorithms.py
@@ -65,6 +65,8 @@
 from testclass import CCPiTestClass
 from utils import has_astra, initialise_tests
 
+from unittest.mock import MagicMock
+
 log = logging.getLogger(__name__)
 initialise_tests()
 
@@ -579,30 +581,103 @@ def test_with_cvxpy(self):
 
 
 class TestCGLS(CCPiTestClass):
-    def test_CGLS(self):
-        ig = ImageGeometry(10, 2)
+    def setUp(self):
+        self.ig = ImageGeometry(10, 2)
         np.random.seed(2)
-        initial = ig.allocate(1.)
-        b = ig.allocate('random')
-        identity = IdentityOperator(ig)
-
-        alg = CGLS(initial=initial, operator=identity, data=b)
-
-        np.testing.assert_array_equal(
-            initial.as_array(), alg.solution.as_array())
-
-        alg.max_iteration = 200
-        alg.run(20, verbose=0)
-        self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array())
-
-        alg = CGLS(initial=initial, operator=identity, data=b,
-                   max_iteration=200, update_objective_interval=2)
-        self.assertTrue(alg.max_iteration == 200)
-        self.assertTrue(alg.update_objective_interval == 2)
-        alg.run(20, verbose=0)
-        self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array())
+        self.initial = self.ig.allocate(1.)
+        self.data = self.ig.allocate('random')
+        self.operator = IdentityOperator(self.ig)
+        self.alg = CGLS(initial=self.initial, operator=self.operator, data=self.data,
+                   update_objective_interval=2)
+
+        
+    def test_initialization_with_default_tolerance(self):
+
+        self.assertEqual(self.alg.tolerance, 0)
+    
+    def test_initialization_with_custom_tolerance(self):
+        with self.assertWarns(DeprecationWarning):
+            alg = CGLS(initial=self.initial, operator=self.operator, data=self.data, tolerance=0.01)
+        self.assertEqual(alg.tolerance, 0.01)
+        
+    def test_set_up(self):
+        # Test the set_up method
+        
+        self.alg.set_up(initial=self.initial, operator=self.operator, data=self.data)
+        
+        # Check if internal variables are set up correctly
+        self.assertNumpyArrayEqual(self.alg.x.as_array(), self.initial.as_array())
+        self.assertEqual(self.alg.operator, self.operator)
+        self.assertNumpyArrayEqual(self.alg.r.as_array(), (self.data-self.initial).as_array())
+        self.assertNumpyArrayEqual(self.alg.s.as_array(), (self.data-self.initial).as_array())
+        self.assertNumpyArrayEqual(self.alg.p.as_array(), (self.data-self.initial).as_array())
+        self.assertTrue(self.alg.configured)
 
+    def test_update(self):
 
+        self.alg.set_up(initial=self.mock_initial, operator=self.mock_operator, data=self.mock_data)
+        
+        self.alg.update()
+        
+        self.mock_operator.direct.assert_called_with(self.mock_data, out=self.alg.q)
+        self.mock_operator.adjoint.assert_called_with(self.alg.r, out=self.alg.s)
+        
+    def test_convergence(self):
+
+        
+        self.alg.run(20, verbose=0)
+        self.assertNumpyArrayAlmostEqual(self.alg.x.as_array(), self.data.as_array())
+        
+    def test_should_stop_flag_false(self):
+        # Mocking norms to ensure tolerance isn't reached
+        self.alg.run(2)
+        self.alg.norms = 10
+        self.alg.norms0 = 99
+        self.alg.tolerance = 0.1
+        self.alg.normx = 0.1
+
+        self.assertFalse(self.alg.flag())
+        
+    def test_should_stop_flag_true(self):
+        # Mocking norms to ensure tolerance is reached
+        self.alg.run(4)
+        self.alg.norms = 1
+        self.alg.norms0 = 10
+        self.alg.tolerance = 0.1
+        self.alg.normx = 10
+
+        self.assertTrue(self.alg.flag())
+        
+    def test_update_objective(self):
+        # Mocking squared_norm to return a finite value
+
+        self.alg.r=MagicMock()
+        self.alg.r.squared_norm.return_value = 1.0
+
+
+        self.alg.update_objective()
+
+        # Ensure the loss list is updated
+        self.assertEqual(self.alg.loss, [1.0])
+        
+    def test_update_objective_with_nan_raises_stop_iteration(self):
+        # Mocking squared_norm to return NaN
+        self.alg.r=MagicMock()
+        self.alg.r.squared_norm.return_value = np.nan
+        
+        with self.assertRaises(StopIteration):
+            self.alg.update_objective()
+            
+    def test_update(self):
+        self.alg.gamma=4
+        
+        self.alg.update()
+        norm=(self.data-self.initial).squared_norm()
+        alpha=4/norm
+        self.assertNumpyArrayEqual(self.alg.x.as_array(), (self.initial+alpha*(self.data-self.initial)).as_array())
+        self.assertNumpyArrayEqual(self.alg.r.as_array(), (self.data - self.initial-alpha*(self.data-self.initial)).as_array())
+        beta= ((self.data - self.initial-alpha*(self.data-self.initial)).norm()**2)/4
+        self.assertNumpyArrayEqual(self.alg.p.as_array(), ((self.data - self.initial-alpha*(self.data-self.initial))+beta*(self.data-self.initial)).as_array())
 class TestPDHG(CCPiTestClass):
 
     def test_PDHG_Denoising(self):
@@ -1298,7 +1373,43 @@ def __call__(self, algorithm: Algorithm):
         algo.run(20, callbacks=[EarlyStopping()])
         self.assertEqual(algo.iteration, 15)
 
-
+    def test_CGLSEarlyStopping(self):
+        ig = ImageGeometry(10, 2)
+        np.random.seed(2)
+        initial = ig.allocate(1.)
+        data = ig.allocate('random')
+        operator = IdentityOperator(ig)
+        alg = CGLS(initial=initial, operator=operator, data=data,
+                   update_objective_interval=2)
+        # Mocking norms to ensure tolerance isn't reached
+        
+        alg.norms = 10
+        alg.norms0 = 99
+        alg.normx = 0.1
+        callbacks.CGLSEarlyStopping(0.1)(alg)
+        
+        # Mocking norms to ensure tolerance is reached
+        alg.norms = 1
+        alg.norms0 = 100
+     
+        with self.assertRaises(StopIteration):
+            callbacks.CGLSEarlyStopping(0.1)(alg)
+
+    def test_EarlyStoppingObjectiveValue(self):
+        ig = ImageGeometry(10, 2)
+        np.random.seed(2)
+        alg = MagicMock()
+        alg.loss=[5]
+        callbacks.EarlyStoppingObjectiveValue(0.1)(alg)
+        alg.loss=[]
+        callbacks.EarlyStoppingObjectiveValue(0.1)(alg)
+        alg.loss=[5,10]
+        callbacks.EarlyStoppingObjectiveValue(0.1)(alg)
+        alg.loss=[5, 5.001]
+        with self.assertRaises(StopIteration):
+            callbacks.EarlyStoppingObjectiveValue(0.1)(alg)
+        
+        
 class TestADMM(unittest.TestCase):
     def setUp(self):
         ig = ImageGeometry(2, 3, 2)

From 533e2c56175fa22cee3e965892376f0e2bbed30a Mon Sep 17 00:00:00 2001
From: Margaret Duff <margaret.duff@stfc.ac.uk>
Date: Thu, 15 Aug 2024 15:31:42 +0000
Subject: [PATCH 3/6] Added comments on depracation to the unit tests

---
 Wrappers/Python/test/test_algorithms.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py
index 888a59d7e..afc1d89f7 100644
--- a/Wrappers/Python/test/test_algorithms.py
+++ b/Wrappers/Python/test/test_algorithms.py
@@ -591,11 +591,11 @@ def setUp(self):
                    update_objective_interval=2)
 
         
-    def test_initialization_with_default_tolerance(self):
+    def test_initialization_with_default_tolerance(self):#can be deprecated when tolerance is deprecated in CGLS
 
         self.assertEqual(self.alg.tolerance, 0)
     
-    def test_initialization_with_custom_tolerance(self):
+    def test_initialization_with_custom_tolerance(self):#can be deprecated when tolerance is deprecated in CGLS
         with self.assertWarns(DeprecationWarning):
             alg = CGLS(initial=self.initial, operator=self.operator, data=self.data, tolerance=0.01)
         self.assertEqual(alg.tolerance, 0.01)
@@ -628,7 +628,7 @@ def test_convergence(self):
         self.alg.run(20, verbose=0)
         self.assertNumpyArrayAlmostEqual(self.alg.x.as_array(), self.data.as_array())
         
-    def test_should_stop_flag_false(self):
+    def test_should_stop_flag_false(self): #can be deprecated when tolerance is deprecated in CGLS
         # Mocking norms to ensure tolerance isn't reached
         self.alg.run(2)
         self.alg.norms = 10
@@ -638,7 +638,7 @@ def test_should_stop_flag_false(self):
 
         self.assertFalse(self.alg.flag())
         
-    def test_should_stop_flag_true(self):
+    def test_should_stop_flag_true(self): #can be deprecated when tolerance is deprecated in CGLS
         # Mocking norms to ensure tolerance is reached
         self.alg.run(4)
         self.alg.norms = 1

From 760523f883526626b7a71934cebeb7df210c5cc4 Mon Sep 17 00:00:00 2001
From: Margaret Duff <margaret.duff@stfc.ac.uk>
Date: Fri, 16 Aug 2024 15:20:40 +0000
Subject: [PATCH 4/6] Removed 'verbose_output()' from CGLS as it is being
 deprecated and causing issue #1894

---
 Wrappers/Python/cil/optimisation/algorithms/CGLS.py | 2 --
 Wrappers/Python/test/test_algorithms.py             | 6 ++++++
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
index cec86b61b..fe5293825 100644
--- a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
@@ -143,8 +143,6 @@ def flag(self): #can be deprecated with tolerance
 
         if flag:
             self.update_objective()
-            if self.iteration > self._iteration[-1]:
-                print (self.verbose_output())
             print('Tolerance is reached: {}'.format(self.tolerance))
 
         return flag
diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py
index afc1d89f7..88d7baed7 100644
--- a/Wrappers/Python/test/test_algorithms.py
+++ b/Wrappers/Python/test/test_algorithms.py
@@ -648,6 +648,12 @@ def test_should_stop_flag_true(self): #can be deprecated when tolerance is depre
 
         self.assertTrue(self.alg.flag())
         
+    def test_tolerance_reached_immediately(self): #can be deprecated when tolerance is deprecated in CGLS
+        alg = CGLS(initial=self.operator.domain_geometry().allocate(0), operator=self.operator, data=self.operator.domain_geometry().allocate(0))
+        alg.run(2)
+
+        
+        
     def test_update_objective(self):
         # Mocking squared_norm to return a finite value
 

From 151c933e4b4e57f500f9994f02b265f4f1fd90f4 Mon Sep 17 00:00:00 2001
From: Margaret Duff <margaret.duff@stfc.ac.uk>
Date: Tue, 20 Aug 2024 12:28:05 +0000
Subject: [PATCH 5/6] Gemma and Jakob's comments

---
 .../cil/optimisation/algorithms/CGLS.py       | 14 ++++----
 .../cil/optimisation/utilities/callbacks.py   | 30 +++++++++++------
 Wrappers/Python/test/test_algorithms.py       | 33 +++++++++++++++----
 3 files changed, 54 insertions(+), 23 deletions(-)

diff --git a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
index fe5293825..f420f88cb 100644
--- a/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/cil/optimisation/algorithms/CGLS.py
@@ -43,7 +43,7 @@ class CGLS(Algorithm):
         Linear operator for the inverse problem
     initial : (optional) DataContainer in the domain of the operator, default is a DataContainer filled with zeros. 
         Initial guess 
-    data : Data Container in the range of the operator 
+    data : DataContainer in the range of the operator 
         Acquired data to reconstruct
         
     Note
@@ -54,7 +54,7 @@ class CGLS(Algorithm):
     ---------
     https://web.stanford.edu/group/SOL/software/cgls/
     '''
-    def __init__(self, initial=None, operator=None, data=None,  **kwargs):
+    def __init__(self, initial=None, operator=None, data=None, **kwargs):
         '''initialisation of the algorithm
         '''
         #We are deprecating tolerance 
@@ -72,14 +72,14 @@ def __init__(self, initial=None, operator=None, data=None,  **kwargs):
             self.set_up(initial=initial, operator=operator, data=data) 
 
     def set_up(self, initial, operator, data):
-        r'''initialisation of the algorithm
+        r'''Initialisation of the algorithm
         Parameters
         ------------
         operator : Operator
             Linear operator for the inverse problem
         initial : (optional) DataContainer in the domain of the operator, default is a DataContainer filled with zeros. 
             Initial guess 
-        data : Data Container in the range of the operator 
+        data : DataContainer in the range of the operator 
             Acquired data to reconstruct
 
         '''
@@ -125,7 +125,7 @@ def update(self):
         #self.p = self.s + self.beta * self.p
         self.p.sapyb(self.beta, self.s, 1, out=self.p)
 
-        self.normx = self.x.norm()#Can be deprecated with tolerance 
+        self.normx = self.x.norm()# TODO: Deprecated, remove when CGLS tolerance is removed
 
 
     def update_objective(self):
@@ -134,10 +134,10 @@ def update_objective(self):
             raise StopIteration()
         self.loss.append(a)
 
-    def should_stop(self): #can be deprecated with tolerance 
+    def should_stop(self): # TODO: Deprecated, remove when CGLS tolerance is removed
         return self.flag() or super().should_stop()
 
-    def flag(self): #can be deprecated with tolerance 
+    def flag(self): # TODO: Deprecated, remove when CGLS tolerance is removed
         '''returns whether the tolerance has been reached'''
         flag  = (self.norms <= self.norms0 * self.tolerance) or (self.normx * self.tolerance >= 1)
 
diff --git a/Wrappers/Python/cil/optimisation/utilities/callbacks.py b/Wrappers/Python/cil/optimisation/utilities/callbacks.py
index b5a79a7f5..8ba79e1aa 100644
--- a/Wrappers/Python/cil/optimisation/utilities/callbacks.py
+++ b/Wrappers/Python/cil/optimisation/utilities/callbacks.py
@@ -145,7 +145,7 @@ class EarlyStoppingObjectiveValue(Callback):
 
     Note
     -----
-    This callback only compares the last two calculated objective values. If `update_objective_interval` is greater than 1, the objective value is not calculated at each iteration. 
+    This callback only compares the last two calculated objective values. f update_objective_interval is greater than 1, the objective value is not calculated at each iteration (which is the default behaviour), only every update_objective_interval iterations.
     
         '''
     def __init__(self, threshold=1e-6):
@@ -158,22 +158,32 @@ def __call__(self, algorithm):
                 raise StopIteration
                 
 class CGLSEarlyStopping(Callback):
-    '''Callback to work with CGLS. It causes the algorithm to terminate if the value of :math:`||A^T(Ax-b)||_2 < tol*||A^T(Ax_0-b)||_2` where 'tol' is set to default as '1e-6', :math:`x` is the current iterate and :math:`x_0` is the initial value. 
-    It will also terminate if the algorithm begins to diverge i.e. if :math:`||x||_2>1/{tol}`. 
+    '''Callback to work with CGLS. It causes the algorithm to terminate if  :math:`||A^T(Ax-b)||_2 < \epsilon||A^T(Ax_0-b)||_2` where `epsilon` is set to default as '1e-6', :math:`x` is the current iterate and :math:`x_0` is the initial value. 
+    It will also terminate if the algorithm begins to diverge i.e. if :math:`||x||_2> \Omega`, where `omega` is set to default as 1e6. 
     Parameters
     ----------
-    tolerance: float, default 1e-6 
- 
+    epsilon: float, default 1e-6 
+        Usually a small number: the algorithm to terminate if :math:`||A^T(Ax-b)||_2 < \epsilon||A^T(Ax_0-b)||_2`
+    lambda: float, default 1e6 
+        Usually a large number: the algorithm will terminate if  :math:`||x||_2> \Omega`
+        
+    Note
+    ---------
+    This callback is implemented to replicate the automatic behaviour of CGLS in CIL versions <=24. It also replicates the behaviour of https://web.stanford.edu/group/SOL/software/cgls/. 
     '''
-    def __init__(self, tolerance=1e-6):
-        self.tolerance=tolerance
+    def __init__(self, epsilon=1e-6, omega=1e6):
+        self.epsilon=epsilon
+        self.omega=omega
     
     
     def __call__(self, algorithm):
-        self.normx = algorithm.x.norm()
         
-        if (algorithm.norms <= algorithm.norms0 * self.tolerance) or (algorithm.normx * self.tolerance >= 1):
-            print('Tolerance is reached: {}'.format(self.tolerance))
+        if (algorithm.norms <= algorithm.norms0 * self.epsilon):
+            print('The norm of the residual is less than {} times the norm of the initial residual and so the algorithm is terminated'.format(self.epsilon))
+            raise StopIteration
+        self.normx = algorithm.x.norm()
+        if algorithm.normx >= self.omega:
+            print('The norm of the solution is greater than {} and so the algorithm is terminated'.format(self.omega))
             raise StopIteration
             
         
diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py
index 88d7baed7..c1b0b0c0e 100644
--- a/Wrappers/Python/test/test_algorithms.py
+++ b/Wrappers/Python/test/test_algorithms.py
@@ -1387,19 +1387,40 @@ def test_CGLSEarlyStopping(self):
         operator = IdentityOperator(ig)
         alg = CGLS(initial=initial, operator=operator, data=data,
                    update_objective_interval=2)
-        # Mocking norms to ensure tolerance isn't reached
+       
+        #Test init
+        cb = callbacks.CGLSEarlyStopping(epsilon = 0.1, omega = 33)
+        self.assertEqual(cb.epsilon, 0.1)
+        self.assertEqual(cb.omega, 33)
+        
+        #Test default values
+        cb = callbacks.CGLSEarlyStopping()
+        self.assertEqual(cb.epsilon, 1e-6)
+        self.assertEqual(cb.omega, 1e6)
         
+        #Tests it doesn't stops iterations if the norm of the current residual is not less than epsilon times the norm of the original residual
         alg.norms = 10
         alg.norms0 = 99
-        alg.normx = 0.1
-        callbacks.CGLSEarlyStopping(0.1)(alg)
+        callbacks.CGLSEarlyStopping(epsilon = 0.1)(alg)
         
-        # Mocking norms to ensure tolerance is reached
+        #Test it stops iterations if the norm of the current residual is less than epsilon times the norm of the original residual
         alg.norms = 1
         alg.norms0 = 100
-     
         with self.assertRaises(StopIteration):
-            callbacks.CGLSEarlyStopping(0.1)(alg)
+            callbacks.CGLSEarlyStopping(epsilon = 0.1)(alg)
+           
+        #Test it doesn't stop iterations if the norm of x is smaller than omega 
+        alg.norms = 10
+        alg.norms0 = 99
+        alg.x = data
+        callbacks.CGLSEarlyStopping(epsilon = 0.1)(alg)
+        
+        #Test it stops iterations if the norm of x is larger than omega
+        alg.norms = 10
+        alg.norms0 = 99
+        alg.x = data
+        with self.assertRaises(StopIteration):
+            callbacks.CGLSEarlyStopping(epsilon = 0.1, omega = 0.33)(alg)
 
     def test_EarlyStoppingObjectiveValue(self):
         ig = ImageGeometry(10, 2)

From 75e1f25c7e638f84533e322d5aac2d4497875a95 Mon Sep 17 00:00:00 2001
From: Margaret Duff <margaret.duff@stfc.ac.uk>
Date: Thu, 22 Aug 2024 08:42:59 +0000
Subject: [PATCH 6/6] Jakob's final comments and change log

---
 CHANGELOG.md                                           |  4 ++++
 .../Python/cil/optimisation/utilities/callbacks.py     | 10 +++++-----
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index eee67f4cc..e8fd2af5c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,12 +3,16 @@
     - Added SVRG and LSVRG stochastic functions (#1625)
     - Added SAG and SAGA stochastic functions (#1624)
     - Allow `SumFunction` with 1 item (#1857)
+    - Added callback `optimisation.utilities.callbacks.EarlyStoppingObjectiveValue` which stops iterations if an algorithm objective changes less than a provided threshold (#1892)
+    - Added callback `optimisation.utilities.callbacks.CGLSEarlyStopping` which replicates the automatic behaviour of CGLS in CIL versions <=24. (#1892)
   - Enhancements:
     - Use ravel instead of flat in KullbackLeibler numba backend (#1874)
     - Upgrade Python wrapper (#1873, #1875)
   - Bug fixes:
     - `ImageData` removes dimensions of size 1 from the input array. This fixes an issue where single slice reconstructions from 3D data would fail due to shape mismatches (#1885)
     - Make Binner accept accelerated=False (#1887)
+  - Changes that break backwards compatibility:
+    - CGLS will no longer automatically stop iterations once a default tolerance is reached. The option to pass `tolerance` will be deprecated to be replaced by `optimisation.utilities.callbacks` (#1892)
 
 * 24.1.0
   - New Features:
diff --git a/Wrappers/Python/cil/optimisation/utilities/callbacks.py b/Wrappers/Python/cil/optimisation/utilities/callbacks.py
index 8ba79e1aa..a5c136695 100644
--- a/Wrappers/Python/cil/optimisation/utilities/callbacks.py
+++ b/Wrappers/Python/cil/optimisation/utilities/callbacks.py
@@ -145,7 +145,7 @@ class EarlyStoppingObjectiveValue(Callback):
 
     Note
     -----
-    This callback only compares the last two calculated objective values. f update_objective_interval is greater than 1, the objective value is not calculated at each iteration (which is the default behaviour), only every update_objective_interval iterations.
+    This callback only compares the last two calculated objective values. If `update_objective_interval` is greater than 1, the objective value is not calculated at each iteration (which is the default behaviour), only every `update_objective_interval` iterations.
     
         '''
     def __init__(self, threshold=1e-6):
@@ -159,16 +159,16 @@ def __call__(self, algorithm):
                 
 class CGLSEarlyStopping(Callback):
     '''Callback to work with CGLS. It causes the algorithm to terminate if  :math:`||A^T(Ax-b)||_2 < \epsilon||A^T(Ax_0-b)||_2` where `epsilon` is set to default as '1e-6', :math:`x` is the current iterate and :math:`x_0` is the initial value. 
-    It will also terminate if the algorithm begins to diverge i.e. if :math:`||x||_2> \Omega`, where `omega` is set to default as 1e6. 
+    It will also terminate if the algorithm begins to diverge i.e. if :math:`||x||_2> \omega`, where `omega` is set to default as 1e6. 
     Parameters
     ----------
     epsilon: float, default 1e-6 
         Usually a small number: the algorithm to terminate if :math:`||A^T(Ax-b)||_2 < \epsilon||A^T(Ax_0-b)||_2`
-    lambda: float, default 1e6 
-        Usually a large number: the algorithm will terminate if  :math:`||x||_2> \Omega`
+    omega: float, default 1e6 
+        Usually a large number: the algorithm will terminate if  :math:`||x||_2> \omega`
         
     Note
-    ---------
+    -----
     This callback is implemented to replicate the automatic behaviour of CGLS in CIL versions <=24. It also replicates the behaviour of https://web.stanford.edu/group/SOL/software/cgls/. 
     '''
     def __init__(self, epsilon=1e-6, omega=1e6):