Reponse to Gemma's review

Wrappers/Python/cil/optimisation/operators/Operator.py

@@ -153,25 +153,9 @@ def adjoint(self, x, out=None):
     def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5,  return_all=False, method='auto'):
         r"""Power method or Power iteration algorithm 
 
-        The power method contains two different algorithms chosen by the `method` flag. 
-
-        In the case `method="direct_only"`, for operator, :math:`A`, the power method computes the iterations 
-        .. math::
-          x_{k+1} = A (\frac{x_k}{\|x_{k}\|}) \; (*)
-
-        initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. 
-
-        In the case `method="composed_with_adjoint"`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` 
-          returning the square root of this value, i.e. the iterations:
-          .. math::
-          x_{k+1} = A^*A (\frac{x_k}{\|x_{k}\|})
-
-        and returning  :math:`\sqrt{\|x_k\|}`.
-
-        The default flag is `method="auto"`, the algorithm checks to see if the `operator.domain_geometry() == operator.range_geometry()` and if so
-        uses the method "direct_only" and if not the method "composed_with_adjoint".
-
-
+        The Power method computes the largest (dominant) eigenvalue of a matrix in magnitude, e.g.,
+        absolute value in the real case and modulus in the complex case.        
+        
         Parameters
         ----------
 
@@ -184,10 +168,9 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5,  retur
             Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance.                    
         return_all: `boolean`, default = False
             Toggles the verbosity of the return
-        compose_with_adjoint: `boolean`, default None
-             Set this to `False` to apply the power method directly on the operator, :math: A :math:, and `True` to apply the power method to :math:A^TA:math: before taking the square root of the result.
-             Leave as default `None` to determine this from domain and range geometry.  
-
+        method: `string` one of `auto`, `composed_with_adjoint` and `direct_only`, default = `auto` 
+            The default `auto` lets the code choose the method, this can be specified with `direct_only` or `compose_with_adjoint`
+            
 
         Returns
         -------
@@ -200,6 +183,27 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5,  retur
             List of eigenvalues. Only returned if return_all is True.
         convergence: `boolean'
             Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True.
+        
+            
+        Note
+        -----
+        The power method contains two different algorithms chosen by the `method` flag. 
+
+        In the case `method="direct_only"`, for operator, :math:`A`, the power method computes the iterations 
+        .. math::
+          x_{k+1} = A (\frac{x_k}{\|x_{k}\|}) \; (*)
+
+        initialised with a random vector :math:`x_0` and returning the largest (dominant) eigenvalue in magnitude given by :math:`\|x_k\|`. 
+
+        In the case `method="composed_with_adjoint"`, the algorithm computes the largest (dominant) eigenvalue of :math:`A^{T}A` 
+          returning the square root of this value, i.e. the iterations:
+          .. math::
+          x_{k+1} = A^*A (\frac{x_k}{\|x_{k}\|})
+
+        and returning  :math:`\sqrt{\|x_k\|}`.
+
+        The default flag is `method="auto"`, the algorithm checks to see if the `operator.domain_geometry() == operator.range_geometry()` and if so
+        uses the method "direct_only" and if not the method "composed_with_adjoint".
 
         Examples
         --------    
@@ -216,20 +220,26 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5,  retur
         2.0005647295658866
 
         """
-        convergence_check = True
 
-        if method not in ["auto","direct_only","composed_with_adjoint"]:
-            ValueError(" The flag 'method' should take one of 'auto','direct_only','composed_with_adjoint'")
+
+        allowed_methods = ["auto","direct_only","composed_with_adjoint"]
+
+        if method not in allowed_methods:
+            raise ValueError("The argument 'method' can be set to one of {0} got {1}".format(allowed_methods, method))
 
+        apply_adjoint=True
+        if method == "direct_only":
+            apply_adjoint=False
         if method=="auto":
             try:
-                if operator.domain_geometry() == operator.range_geometry():
-                    method="direct_only"
-                else:
-                    method="composed_with_adjoint"
+                geometries_match = operator.domain_geometry() == operator.range_geometry()
+
             except AssertionError:
                 # catch AssertionError for SIRF objects https://github.com/SyneRBI/SIRF-SuperBuild/runs/5110228626?check_suite_focus=true#step:8:972
-                method="composed_with_adjoint"
+                    pass
+            else:
+                if geometries_match:
+                    apply_adjoint=False
 
 
         if initial is None:
@@ -245,19 +255,20 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5,  retur
 
         # initial guess for dominant eigenvalue
         eig_old = 1.
-
-        eig_list = []
+        if return_all:
+            eig_list = []
+            convergence_check = True
         diff = numpy.finfo('d').max
         i = 0
         while (i < max_iteration and diff > tolerance):
             operator.direct(x0, out=y_tmp)
 
-            if method=="direct_only":
+            if not apply_adjoint:
                 # swap datacontainer references
                 tmp = x0
                 x0 = y_tmp
                 y_tmp = tmp
-            elif method=="composed_with_adjoint":
+            else:
                 operator.adjoint(y_tmp, out=x0)
 
             # Get eigenvalue using Rayleigh quotient: denominator=1, due to normalization
@@ -270,14 +281,15 @@ def PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-5,  retur
             x0 /= x0_norm
 
             eig_new = numpy.abs(x0_norm)
-            if method=="composed_with_adjoint":
+            if apply_adjoint:
                 eig_new = numpy.sqrt(eig_new)
             diff = numpy.abs(eig_new - eig_old)
-            eig_list.append(eig_new)
+            if return_all:
+                eig_list.append(eig_new)
             eig_old = eig_new
             i += 1
 
-        if i == max_iteration:
+        if return_all and i == max_iteration:
             convergence_check = False
 
         if return_all: