added bool configured to Algorithm

Algorithms will not run if not configured  (meaning set_up has been run).
closes #304

7 files changed, 47 insertions, 32 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
index 4fbf83b..c62d0ea 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
@@ -51,6 +51,7 @@ class Algorithm(object):
         self.__max_iteration = kwargs.get('max_iteration', 0)
         self.__loss = []
         self.memopt = False
+        self.configured = False
         self.timing = []
         self.update_objective_interval = kwargs.get('update_objective_interval', 1)
     def set_up(self, *args, **kwargs):
@@ -86,6 +87,8 @@ class Algorithm(object):
             raise StopIteration()
         else:
             time0 = time.time()
+            if not self.configured:
+                raise ValueError('Algorithm not configured correctly. Please run set_up.')
             self.update()
             self.timing.append( time.time() - time0 )
             if self.iteration % self.update_objective_interval == 0:
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
index 4d4843c..6b610a0 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
@@ -23,6 +23,7 @@ Created on Thu Feb 21 11:11:23 2019
 """
 
 from ccpi.optimisation.algorithms import Algorithm
+from ccpi.optimisation.functions import Norm2Sq
 
 class CGLS(Algorithm):
 
@@ -59,6 +60,9 @@ class CGLS(Algorithm):
         #    self.normr2 = sum(self.normr2)
         #self.normr2 = numpy.sqrt(self.normr2)
         #print ("set_up" , self.normr2)
+        n = Norm2Sq(operator, self.data)
+        self.loss.append(n(x_init))
+        self.configured = True
 
     def update(self):
 
@@ -84,4 +88,4 @@ class CGLS(Algorithm):
         self.d = s + beta*self.d
 
     def update_objective(self):
-        self.loss.append(self.r.squared_norm())
\ No newline at end of file
+        self.loss.append(self.r.squared_norm())
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py
index aa07359..269438c 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py
@@ -35,52 +35,48 @@ class FBPD(Algorithm):
       h: regularizer
       opt: additional algorithm 
     '''
-    constraint = None
-    data_fidelity = None
-    regulariser = None
     def __init__(self, **kwargs):
-        pass
-    def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\
-         regulariser=None, opt=None):
+        super(FBPD, self).__init__()
+        self.f = kwargs.get('f', None)
+        self.g = kwargs.get('g', ZeroFunction())
+        self.g = kwargs.get('h', ZeroFunction())
+        self.operator = kwargs.get('operator', None)
+        self.x_init = kwargs.get('x_init',None)
+        if self.x_init is not None and self.operator is not None:
+            self.set_up(self.x_init, self.operator, self.f, self.g, self.h)
+
+    def set_up(self, x_init, operator, constraint, data_fidelity, 
+               regulariser, opt=None):
 
-        # default inputs
-        if constraint    is None: 
-            self.constraint    = ZeroFun()
-        else:
-            self.constraint = constraint
-        if data_fidelity is None:
-            data_fidelity = ZeroFun()
-        else:
-            self.data_fidelity = data_fidelity
-        if regulariser   is None:
-            self.regulariser   = ZeroFun()
-        else:
-            self.regulariser = regulariser
         
         # algorithmic parameters
         
         
         # step-sizes
-        self.tau   = 2 / (self.data_fidelity.L + 2)
-        self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L
+        self.tau   = 2 / (data_fidelity.L + 2)
+        self.sigma = (1/self.tau - data_fidelity.L/2) / regulariser.L
         
         self.inv_sigma = 1/self.sigma
     
         # initialization
         self.x = x_init
         self.y = operator.direct(self.x)
+        self.update_objective()
+        self.configured = True
         
     
     def update(self):
     
         # primal forward-backward step
         x_old = self.x
-        self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) )
-        self.x = self.constraint.prox(self.x, self.tau);
+        self.x = self.x - self.tau * ( self.g.gradient(self.x) + self.operator.adjoint(self.y) )
+        self.x = self.f.proximal(self.x, self.tau)
     
         # dual forward-backward step
-        self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old);
-        self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma);   
+        self.y += self.sigma * self.operator.direct(2*self.x - x_old);
+        self.y -= self.sigma * self.h.proximal(self.inv_sigma*self.y, self.inv_sigma)   
 
         # time and criterion
-        self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x))
+
+    def update_objective(self):
+        self.loss.append(self.f(self.x) + self.g(self.x) + self.h(self.operator.direct(self.x)))
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
index 3c7a8d1..647ae98 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
@@ -58,6 +58,7 @@ class FISTA(Algorithm):
         
         self.t_old = 1
         self.update_objective()
+        self.configured = True
             
     def update(self):
 
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
index 14763c5..34bf954 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
@@ -40,7 +40,7 @@ class GradientDescent(Algorithm):
             if k in args:
                 args.pop(args.index(k))
         if len(args) == 0:
-            return self.set_up(x_init=kwargs['x_init'],
+            self.set_up(x_init=kwargs['x_init'],
                                objective_function=kwargs['objective_function'],
                                rate=kwargs['rate'])
     
@@ -61,6 +61,7 @@ class GradientDescent(Algorithm):
             self.memopt = False
         if self.memopt:
             self.x_update = x_init.copy()
+        self.configured = True
 
     def update(self):
         '''Single iteration'''
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index 39b092b..3afd8b0 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -23,10 +23,16 @@ class PDHG(Algorithm):
         self.operator = kwargs.get('operator', None)
         self.g        = kwargs.get('g', None)
         self.tau      = kwargs.get('tau', None)
-        self.sigma    = kwargs.get('sigma', None)
+        self.sigma    = kwargs.get('sigma', 1.)
+        
         
         if self.f is not None and self.operator is not None and \
            self.g is not None:
+            if self.tau is None:
+                # Compute operator Norm
+                normK = self.operator.norm()
+                # Primal & dual stepsizes
+                self.tau = 1/(self.sigma*normK**2)
             print ("Calling from creator")
             self.set_up(self.f,
                         self.g,
@@ -57,6 +63,8 @@ class PDHG(Algorithm):
 
         # relaxation parameter
         self.theta = 1
+        self.update_objective()
+        self.configured = True
 
     def update(self):
         
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
index 30584d4..c73d323 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
@@ -59,6 +59,7 @@ class SIRT(Algorithm):
         # Set up scaling matrices D and M.
         self.M = 1/self.operator.direct(self.operator.domain_geometry().allocate(value=1.0))        
         self.D = 1/self.operator.adjoint(self.operator.range_geometry().allocate(value=1.0))
+        self.configured = True
 
 
     def update(self):
@@ -67,8 +68,9 @@ class SIRT(Algorithm):
         
         self.x += self.relax_par * (self.D*self.operator.adjoint(self.M*self.r))
         
-        if self.constraint != None:
-            self.x = self.constraint.prox(self.x,None)
+        if self.constraint is not None:
+            self.x = self.constraint.proximal(self.x,None)
+            # self.constraint.proximal(self.x,None, out=self.x)
 
     def update_objective(self):
-        self.loss.append(self.r.squared_norm())
\ No newline at end of file
+        self.loss.append(self.r.squared_norm())