DRAFT: adaptive PDHG and SPDHG step sizes

Wrappers/Python/cil/framework/block.py

@@ -475,7 +475,6 @@ def binary_operations(self, operation, other, *args, **kwargs):
 
                     # As axpyb cannot return anything we `continue` to skip the rest of the code block
                     continue
-
                 else:
                     raise ValueError('Unsupported operation', operation)
                 if out is not None:
@@ -600,7 +599,7 @@ def __div__(self, other):
     # __rdiv__
     def __truediv__(self, other):
         return self.divide(other)
-
+    
     def __pow__(self, other):
         return self.power(other)
     # reverse operand

Wrappers/Python/cil/framework/data_container.py

@@ -567,6 +567,8 @@ def divide(self, other, *args, **kwargs):
             return _out
         return self.pixel_wise_binary(numpy.divide, other, *args, **kwargs)
 
+
+
     def power(self, other, *args, **kwargs):
         return self.pixel_wise_binary(numpy.power, other, *args, **kwargs)
 
@@ -575,6 +577,8 @@ def maximum(self, x2, *args, **kwargs):
 
     def minimum(self,x2, out=None, *args, **kwargs):
         return self.pixel_wise_binary(numpy.minimum, x2=x2, out=out, *args, **kwargs)
+
+
 
 
     def sapyb(self, a, y, b, out=None, num_threads=NUM_THREADS):

Wrappers/Python/cil/optimisation/algorithms/Algorithm.py

@@ -51,6 +51,18 @@ def __init__(self, update_objective_interval=1):
         self.update_objective_interval = update_objective_interval
         self.iter_string = 'Iter'
 
+    def _reset_iteration_state(self):
+        '''Resets the iteration counter and the objective/iteration history.
+
+        This is intended for internal use by step-size rules that re-run the
+        algorithm during set-up (e.g. Bayesian optimisation of the step sizes),
+        so that they do not need to reach into the algorithm's private state.
+        '''
+        self.iteration = -1
+        self.__loss = []
+        self._iteration = []
+        self._total_iterations = 1
+
     def set_up(self, *args, **kwargs):
         '''Set up the algorithm'''
         raise NotImplementedError
@@ -185,14 +197,14 @@ def loss(self):
     @property
     def update_objective_interval(self):
         '''gets the update_objective_interval'''
-        return self.__update_objective_interval
+        return self._update_objective_interval
 
     @update_objective_interval.setter
     def update_objective_interval(self, value):
         '''sets the update_objective_interval'''
         if not isinstance(value, Integral) or value < 0:
             raise ValueError('interval must be an integer >= 0')
-        self.__update_objective_interval = value
+        self._update_objective_interval = value
 
     def run(self, iterations=None, callbacks: Optional[List[Callback]] = None, verbose=1):
         r"""run upto :code:`iterations` with callbacks/logging.

Wrappers/Python/cil/optimisation/algorithms/PD3O.py

@@ -125,7 +125,6 @@ class PD3O(Algorithm):
         Yan, M. A New Primal–Dual Algorithm for Minimizing the Sum of Three Functions with a Linear Operator. J Sci Comput 76, 1698–1717 (2018). https://doi.org/10.1007/s10915-018-0680-3
      """
 
-
     def __init__(self, f, g, h, operator, delta=None, gamma=None, initial=None, **kwargs):
 
         super(PD3O, self).__init__(**kwargs)
@@ -139,7 +138,8 @@ def set_up(self, f, g, h, operator, delta=None, gamma=None, initial=None,**kwarg
         logging.info("{} setting up".format(self.__class__.__name__, ))
 
         if isinstance(f, ZeroFunction):
-            warnings.warn(" If f is the ZeroFunction, then PD3O = PDHG. Please use PDHG instead. Otherwise, select a relatively small parameter gamma ", UserWarning)
+            warnings.warn(
+                " If f is the ZeroFunction, then PD3O = PDHG. Please use PDHG instead. Otherwise, select a relatively small parameter gamma ", UserWarning)
             if gamma is None:
                 gamma = 1.0/operator.norm()
 
@@ -177,15 +177,13 @@ def set_up(self, f, g, h, operator, delta=None, gamma=None, initial=None,**kwarg
         self.s_old.sapyb(1, self.s, self.delta, out=self.s_old)
         self.h.proximal_conjugate(self.s_old, self.delta, out=self.s)
 
-
     def update(self):
         r""" Performs a single iteration of the PD3O algorithm
         """
 
         # Following equations 4 in https://link.springer.com/article/10.1007/s10915-018-0680-3
         # in this case order of proximal steps we recover the (primal) PDHG, when f=0
 
-
         tmp = self.x_old
         self.x_old = self.x
         self.x = tmp