delete algs/funcs

2 files changed, 0 insertions, 572 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py
deleted file mode 100755
index f5ba85e..0000000
--- a/Wrappers/Python/ccpi/optimisation/algs.py
+++ /dev/null
@@ -1,307 +0,0 @@
-# -*- coding: utf-8 -*-
-#   This work is part of the Core Imaging Library developed by
-#   Visual Analytics and Imaging System Group of the Science Technology
-#   Facilities Council, STFC
-
-#   Copyright 2018 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca
-
-#   Licensed under the Apache License, Version 2.0 (the "License");
-#   you may not use this file except in compliance with the License.
-#   You may obtain a copy of the License at
-
-#       http://www.apache.org/licenses/LICENSE-2.0
-
-#   Unless required by applicable law or agreed to in writing, software
-#   distributed under the License is distributed on an "AS IS" BASIS,
-#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-#   See the License for the specific language governing permissions and
-#   limitations under the License.
-
-import numpy
-import time
-
-
-
-
-def FISTA(x_init, f=None, g=None, opt=None):
-    '''Fast Iterative Shrinkage-Thresholding Algorithm
-    
-    Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding 
-    algorithm for linear inverse problems. 
-    SIAM journal on imaging sciences,2(1), pp.183-202.
-    
-    Parameters:
-      x_init: initial guess
-      f: data fidelity
-      g: regularizer
-      h:
-      opt: additional algorithm 
-    '''
-    # default inputs
-    if f   is None: f = ZeroFun()
-    if g   is None: g = ZeroFun()
-    
-    # algorithmic parameters
-    if opt is None: 
-        opt = {'tol': 1e-4, 'iter': 1000, 'memopt':False}
-    
-    max_iter = opt['iter'] if 'iter' in opt.keys() else 1000
-    tol = opt['tol'] if 'tol' in opt.keys() else 1e-4
-    memopt = opt['memopt'] if 'memopt' in opt.keys() else False
-        
-        
-    # initialization
-    if memopt:
-        y = x_init.clone()
-        x_old = x_init.clone()
-        x = x_init.clone()
-        u = x_init.clone()
-    else:
-        x_old = x_init
-        y = x_init;
-    
-    timing = numpy.zeros(max_iter)
-    criter = numpy.zeros(max_iter)
-    
-    invL = 1/f.L
-    
-    t_old = 1
-    
-#    c = f(x_init) + g(x_init)
-
-    # algorithm loop
-    for it in range(0, max_iter):
-    
-        time0 = time.time()
-        if memopt:
-            # u = y - invL*f.grad(y)
-            # store the result in x_old
-            f.gradient(y, out=u)
-            u.__imul__( -invL )
-            u.__iadd__( y )
-            # x = g.prox(u,invL)
-            g.proximal(u, invL, out=x)
-            
-            t = 0.5*(1 + numpy.sqrt(1 + 4*(t_old**2)))
-            
-            # y = x + (t_old-1)/t*(x-x_old)
-            x.subtract(x_old, out=y)
-            y.__imul__ ((t_old-1)/t)
-            y.__iadd__( x )
-            
-            x_old.fill(x)
-            t_old = t
-            
-            
-        else:
-            u = y - invL*f.gradient(y)
-            
-            x = g.proximal(u,invL)
-            
-            t = 0.5*(1 + numpy.sqrt(1 + 4*(t_old**2)))
-            
-            y = x + (t_old-1)/t*(x-x_old)
-            
-            x_old = x.copy()
-            t_old = t
-        
-        # time and criterion
-#        timing[it] = time.time() - time0
-#        criter[it] = f(x) + g(x);
-        
-        # stopping rule
-        #if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10:
-        #   break
-    
-        #print(it, 'out of', 10, 'iterations', end='\r');
-
-    #criter = criter[0:it+1];
-#    timing = numpy.cumsum(timing[0:it+1]);
-    
-    return x #, it, timing, criter
-
-def FBPD(x_init, operator=None, constraint=None, data_fidelity=None,\
-         regulariser=None, opt=None):
-    '''FBPD Algorithm
-    
-    Parameters:
-      x_init: initial guess
-      f: constraint
-      g: data fidelity
-      h: regularizer
-      opt: additional algorithm 
-    '''
-    # default inputs
-    if constraint    is None: constraint    = ZeroFun()
-    if data_fidelity is None: data_fidelity = ZeroFun()
-    if regulariser   is None: regulariser   = ZeroFun()
-    
-    # algorithmic parameters
-    if opt is None: 
-        opt = {'tol': 1e-4, 'iter': 1000}
-    else:
-        try:
-            max_iter = opt['iter']
-        except KeyError as ke:
-            opt[ke] = 1000
-        try:
-            opt['tol'] = 1000
-        except KeyError as ke:
-            opt[ke] = 1e-4
-    tol = opt['tol']
-    max_iter = opt['iter']
-    memopt = opt['memopts'] if 'memopts' in opt.keys() else False
-    
-    # step-sizes
-    tau   = 2 / (data_fidelity.L + 2)
-    sigma = (1/tau - data_fidelity.L/2) / regulariser.L
-    inv_sigma = 1/sigma
-
-    # initialization
-    x = x_init
-    y = operator.direct(x);
-    
-    timing = numpy.zeros(max_iter)
-    criter = numpy.zeros(max_iter)
-
-    
-    
-        
-    # algorithm loop
-    for it in range(0, max_iter):
-    
-        t = time.time()
-    
-        # primal forward-backward step
-        x_old = x;
-        x = x - tau * ( data_fidelity.grad(x) + operator.adjoint(y) );
-        x = constraint.prox(x, tau);
-    
-        # dual forward-backward step
-        y = y + sigma * operator.direct(2*x - x_old);
-        y = y - sigma * regulariser.prox(inv_sigma*y, inv_sigma);   
-
-        # time and criterion
-        timing[it] = time.time() - t
-        criter[it] = constraint(x) + data_fidelity(x) + regulariser(operator.direct(x))
-           
-        # stopping rule
-        #if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10:
-        #   break
-
-    criter = criter[0:it+1]
-    timing = numpy.cumsum(timing[0:it+1])
-    
-    return x, it, timing, criter
-
-def CGLS(x_init, operator , data , opt=None):
-    '''Conjugate Gradient Least Squares algorithm
-    
-    Parameters:
-      x_init: initial guess
-      operator: operator for forward/backward projections
-      data: data to operate on
-      opt: additional algorithm 
-    '''
-    
-    if opt is None: 
-        opt = {'tol': 1e-4, 'iter': 1000}
-    else:
-        try:
-            max_iter = opt['iter']
-        except KeyError as ke:
-            opt[ke] = 1000
-        try:
-            opt['tol'] = 1000
-        except KeyError as ke:
-            opt[ke] = 1e-4
-    tol = opt['tol']
-    max_iter = opt['iter']
-    
-    r = data.copy()
-    x = x_init.copy()
-    
-    d = operator.adjoint(r)
-    
-    normr2 = (d**2).sum()
-    
-    timing = numpy.zeros(max_iter)
-    criter = numpy.zeros(max_iter)
-
-    # algorithm loop
-    for it in range(0, max_iter):
-    
-        t = time.time()
-        
-        Ad = operator.direct(d)
-        alpha = normr2/( (Ad**2).sum() )
-        x  = x + alpha*d
-        r  = r - alpha*Ad
-        s  = operator.adjoint(r)
-        
-        normr2_new = (s**2).sum()
-        beta = normr2_new/normr2
-        normr2 = normr2_new
-        d = s + beta*d
-        
-        # time and criterion
-        timing[it] = time.time() - t
-        criter[it] = (r**2).sum()
-    
-    return x, it, timing, criter
-
-def SIRT(x_init, operator , data , opt=None, constraint=None):
-    '''Simultaneous Iterative Reconstruction Technique
-    
-    Parameters:
-      x_init: initial guess
-      operator: operator for forward/backward projections
-      data: data to operate on
-      opt: additional algorithm 
-      constraint: func of Indicator type specifying convex constraint.
-    '''
-    
-    if opt is None: 
-        opt = {'tol': 1e-4, 'iter': 1000}
-    else:
-        try:
-            max_iter = opt['iter']
-        except KeyError as ke:
-            opt[ke] = 1000
-        try:
-            opt['tol'] = 1000
-        except KeyError as ke:
-            opt[ke] = 1e-4
-    tol = opt['tol']
-    max_iter = opt['iter']
-    
-    x = x_init.clone()
-    
-    timing = numpy.zeros(max_iter)
-    criter = numpy.zeros(max_iter)
-    
-    # Relaxation parameter must be strictly between 0 and 2. For now fix at 1.0
-    relax_par = 1.0
-    
-    # Set up scaling matrices D and M.
-    M = 1/operator.direct(operator.domain_geometry().allocate(value=1.0))
-    D = 1/operator.adjoint(operator.range_geometry().allocate(value=1.0))
-    
-    # algorithm loop
-    for it in range(0, max_iter):
-        t = time.time()
-        r = data - operator.direct(x)
-        
-        x = x + relax_par * (D*operator.adjoint(M*r))
-
-        if constraint != None:
-            x = constraint.prox(x,None)
-        
-        timing[it] = time.time() - t
-        if it > 0:
-            criter[it-1] = (r**2).sum()
-    
-    r = data - operator.direct(x)
-    criter[it] = (r**2).sum()
-    return x, it, timing,  criter
-    
diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py
deleted file mode 100755
index b2b9791..0000000
--- a/Wrappers/Python/ccpi/optimisation/funcs.py
+++ /dev/null
@@ -1,265 +0,0 @@
-# -*- coding: utf-8 -*-
-#   This work is part of the Core Imaging Library developed by
-#   Visual Analytics and Imaging System Group of the Science Technology
-#   Facilities Council, STFC
-
-#   Copyright 2018 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca
-
-#   Licensed under the Apache License, Version 2.0 (the "License");
-#   you may not use this file except in compliance with the License.
-#   You may obtain a copy of the License at
-
-#       http://www.apache.org/licenses/LICENSE-2.0
-
-#   Unless required by applicable law or agreed to in writing, software
-#   distributed under the License is distributed on an "AS IS" BASIS,
-#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-#   See the License for the specific language governing permissions and
-#   limitations under the License.
-
-import numpy
-from ccpi.framework import DataContainer
-import warnings
-from ccpi.optimisation.functions import Function
-def isSizeCorrect(data1 ,data2):
-    if issubclass(type(data1), DataContainer) and \
-       issubclass(type(data2), DataContainer):
-        # check dimensionality
-        if data1.check_dimensions(data2):
-            return True
-    elif issubclass(type(data1) , numpy.ndarray) and \
-         issubclass(type(data2) , numpy.ndarray):
-        return data1.shape == data2.shape
-    else:
-        raise ValueError("{0}: getting two incompatible types: {1} {2}"\
-                         .format('Function', type(data1), type(data2)))
-    return False
-class Norm2(Function):
-    
-    def __init__(self, 
-                 gamma=1.0, 
-                 direction=None):
-        super(Norm2, self).__init__()
-        self.gamma     = gamma;
-        self.direction = direction; 
-    
-    def __call__(self, x, out=None):
-        
-        if out is None:
-            xx = numpy.sqrt(numpy.sum(numpy.square(x.as_array()), self.direction,
-                                  keepdims=True))
-        else:
-            if isSizeCorrect(out, x):
-                # check dimensionality
-                if issubclass(type(out), DataContainer):
-                    arr = out.as_array()
-                    numpy.square(x.as_array(), out=arr)
-                    xx = numpy.sqrt(numpy.sum(arr, self.direction, keepdims=True))
-                        
-                elif issubclass(type(out) , numpy.ndarray):
-                    numpy.square(x.as_array(), out=out)
-                    xx = numpy.sqrt(numpy.sum(out, self.direction, keepdims=True))
-            else:
-                raise ValueError ('Wrong size: x{0} out{1}'.format(x.shape,out.shape) )
-        
-        p  = numpy.sum(self.gamma*xx)        
-        
-        return p
-    
-    def prox(self, x, tau):
-
-        xx = numpy.sqrt(numpy.sum( numpy.square(x.as_array()), self.direction, 
-                                  keepdims=True ))
-        xx = numpy.maximum(0, 1 - tau*self.gamma / xx)
-        p  = x.as_array() * xx
-        
-        return type(x)(p,geometry=x.geometry)
-    def proximal(self, x, tau, out=None):
-        if out is None:
-            return self.prox(x,tau)
-        else:
-            if isSizeCorrect(out, x):
-                # check dimensionality
-                if issubclass(type(out), DataContainer):
-                    numpy.square(x.as_array(), out = out.as_array())
-                    xx = numpy.sqrt(numpy.sum( out.as_array() , self.direction, 
-                                  keepdims=True ))
-                    xx = numpy.maximum(0, 1 - tau*self.gamma / xx)
-                    x.multiply(xx, out= out.as_array())
-                    
-                        
-                elif issubclass(type(out) , numpy.ndarray):
-                    numpy.square(x.as_array(), out=out)
-                    xx = numpy.sqrt(numpy.sum(out, self.direction, keepdims=True))
-                    
-                    xx = numpy.maximum(0, 1 - tau*self.gamma / xx)
-                    x.multiply(xx, out= out)
-            else:
-                raise ValueError ('Wrong size: x{0} out{1}'.format(x.shape,out.shape) )
-        
-
-        
-
-# Define a class for squared 2-norm
-class Norm2sq(Function):
-    '''
-    f(x) = c*||A*x-b||_2^2
-    
-    which has 
-    
-    grad[f](x) = 2*c*A^T*(A*x-b)
-    
-    and Lipschitz constant
-    
-    L = 2*c*||A||_2^2 = 2*s1(A)^2
-    
-    where s1(A) is the largest singular value of A.
-    
-    '''
-    
-    def __init__(self,A,b,c=1.0,memopt=False):
-        super(Norm2sq, self).__init__()
-    
-        self.A = A  # Should be an operator, default identity
-        self.b = b  # Default zero DataSet?
-        self.c = c  # Default 1.
-        if memopt:
-            try:
-                self.range_tmp = A.range_geometry().allocate()
-                self.domain_tmp = A.domain_geometry().allocate()
-                self.memopt = True
-            except NameError as ne:
-                warnings.warn(str(ne))
-                self.memopt = False
-            except NotImplementedError as nie:
-                print (nie)
-                warnings.warn(str(nie))
-                self.memopt = False
-        else:
-            self.memopt = False
-        
-        # Compute the Lipschitz parameter from the operator if possible
-        # Leave it initialised to None otherwise
-        try:
-            self.L = 2.0*self.c*(self.A.norm()**2)
-        except AttributeError as ae:
-            pass
-        except NotImplementedError as noe:
-            pass
-        
-    #def grad(self,x):
-    #    return self.gradient(x, out=None)
-
-    def __call__(self,x):
-        #return self.c* np.sum(np.square((self.A.direct(x) - self.b).ravel()))
-        #if out is None:
-        #    return self.c*( ( (self.A.direct(x)-self.b)**2).sum() )
-        #else:
-        y = self.A.direct(x)
-        y.__isub__(self.b)
-        #y.__imul__(y)
-        #return y.sum() * self.c
-        try:
-            return y.squared_norm() * self.c
-        except AttributeError as ae:
-            # added for compatibility with SIRF 
-            return (y.norm()**2) * self.c
-    
-    def gradient(self, x, out = None):
-        if self.memopt:
-            #return 2.0*self.c*self.A.adjoint( self.A.direct(x) - self.b )
-            
-            self.A.direct(x, out=self.range_tmp)
-            self.range_tmp -= self.b 
-            self.A.adjoint(self.range_tmp, out=out)
-            #self.direct_placehold.multiply(2.0*self.c, out=out)
-            out *= (self.c * 2.0)
-        else:
-            return (2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b )
-
-
-
-# Box constraints indicator function. Calling returns 0 if argument is within 
-# the box. The prox operator is projection onto the box. Only implements one 
-# scalar lower and one upper as constraint on all elements. Should generalise 
-# to vectors to allow different constraints one elements.
-class IndicatorBox(Function):
-    
-    def __init__(self,lower=-numpy.inf,upper=numpy.inf):
-        # Do nothing
-        super(IndicatorBox, self).__init__()
-        self.lower = lower
-        self.upper = upper
-        
-    
-    def __call__(self,x):
-        
-        if (numpy.all(x.array>=self.lower) and 
-            numpy.all(x.array <= self.upper) ):
-            val = 0
-        else:
-            val = numpy.inf
-        return val
-    
-    def prox(self,x,tau=None):
-        return  (x.maximum(self.lower)).minimum(self.upper)
-    
-    def proximal(self, x, tau, out=None):
-        if out is None:
-            return self.prox(x, tau)
-        else:
-            if not x.shape == out.shape:
-                raise ValueError('Norm1 Incompatible size:',
-                                 x.shape, out.shape)
-            #(x.abs() - tau*self.gamma).maximum(0) * x.sign()
-            x.abs(out = out)
-            out.__isub__(tau*self.gamma)
-            out.maximum(0, out=out)
-            if self.sign_x is None or not x.shape == self.sign_x.shape:
-                self.sign_x = x.sign()
-            else:
-                x.sign(out=self.sign_x)
-                
-            out.__imul__( self.sign_x )
-
-# A more interesting example, least squares plus 1-norm minimization.
-# Define class to represent 1-norm including prox function
-class Norm1(Function):
-    
-    def __init__(self,gamma):
-        super(Norm1, self).__init__()
-        self.gamma = gamma
-        self.L = 1
-        self.sign_x = None
-    
-    def __call__(self,x,out=None):
-        if out is None:
-            return self.gamma*(x.abs().sum())
-        else:
-            if not x.shape == out.shape:
-                raise ValueError('Norm1 Incompatible size:',
-                                 x.shape, out.shape)
-            x.abs(out=out)
-            return out.sum() * self.gamma
-    
-    def prox(self,x,tau):
-        return (x.abs() - tau*self.gamma).maximum(0) * x.sign()
-    
-    def proximal(self, x, tau, out=None):
-        if out is None:
-            return self.prox(x, tau)
-        else:
-            if isSizeCorrect(x,out):
-                # check dimensionality
-                if issubclass(type(out), DataContainer):
-                    v = (x.abs() - tau*self.gamma).maximum(0)
-                    x.sign(out=out)
-                    out *= v
-                    #out.fill(self.prox(x,tau))    
-                elif issubclass(type(out) , numpy.ndarray):
-                    v = (x.abs() - tau*self.gamma).maximum(0)
-                    out[:] = x.sign()
-                    out *= v
-                    #out[:] = self.prox(x,tau)
-            else:
-                raise ValueError ('Wrong size: x{0} out{1}'.format(x.shape,out.shape) )