diff options
Diffstat (limited to 'Wrappers/Python')
| -rw-r--r-- | Wrappers/Python/ccpi/framework.py | 244 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algs.py | 78 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/funcs.py | 139 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/ops.py | 114 | ||||
| -rw-r--r-- | Wrappers/Python/conda-recipe/meta.yaml | 4 | ||||
| -rwxr-xr-x | Wrappers/Python/conda-recipe/run_test.py | 536 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_compare_cvx.py | 25 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/demo_memhandle.py | 193 |
8 files changed, 1200 insertions, 133 deletions
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 0c2432f..a34ca37 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -33,6 +33,15 @@ def find_key(dic, val): """return the key of dictionary dic given the value""" return [k for k, v in dic.items() if v == val][0] +def message(cls, msg, *args): + msg = "{0}: " + msg + for i in range(len(args)): + msg += " {%d}" %(i+1) + args = list(args) + args.insert(0, cls.__name__ ) + + return msg.format(*args ) + class ImageGeometry(object): @@ -211,7 +220,7 @@ class DataContainer(object): if type(array) == numpy.ndarray: if deep_copy: - self.array = array[:] + self.array = array.copy() else: self.array = array else: @@ -335,12 +344,17 @@ class DataContainer(object): def fill(self, array, **dimension): '''fills the internal numpy array with the one provided''' if dimension == {}: - if numpy.shape(array) != numpy.shape(self.array): - raise ValueError('Cannot fill with the provided array.' + \ - 'Expecting {0} got {1}'.format( - numpy.shape(self.array), - numpy.shape(array))) - self.array = array[:] + if issubclass(type(array), DataContainer) or\ + issubclass(type(array), numpy.ndarray): + if array.shape != self.shape: + raise ValueError('Cannot fill with the provided array.' + \ + 'Expecting {0} got {1}'.format( + self.shape,array.shape)) + if issubclass(type(array), DataContainer): + numpy.copyto(self.array, array.array) + else: + #self.array[:] = array + numpy.copyto(self.array, array) else: command = 'self.array[' @@ -360,8 +374,10 @@ class DataContainer(object): def check_dimensions(self, other): return self.shape == other.shape - - def __add__(self, other): + + ## algebra + + def __add__(self, other , *args, out=None, **kwargs): if issubclass(type(other), DataContainer): if self.check_dimensions(other): out = self.as_array() + other.as_array() @@ -407,7 +423,6 @@ class DataContainer(object): return self.__div__(other) def __div__(self, other): - print ("calling __div__") if issubclass(type(other), DataContainer): if self.check_dimensions(other): out = self.as_array() / other.as_array() @@ -470,40 +485,6 @@ class DataContainer(object): type(other))) # __mul__ - - #def __abs__(self): - # operation = FM.OPERATION.ABS - # return self.callFieldMath(operation, None, self.mask, self.maskOnValue) - # __abs__ - - def abs(self): - out = numpy.abs(self.as_array() ) - return type(self)(out, - deep_copy=True, - dimension_labels=self.dimension_labels, - geometry=self.geometry) - - def maximum(self,otherscalar): - out = numpy.maximum(self.as_array(),otherscalar) - return type(self)(out, - deep_copy=True, - dimension_labels=self.dimension_labels, - geometry=self.geometry) - - def minimum(self,otherscalar): - out = numpy.minimum(self.as_array(),otherscalar) - return type(self)(out, - deep_copy=True, - dimension_labels=self.dimension_labels, - geometry=self.geometry) - - def sign(self): - out = numpy.sign(self.as_array() ) - return type(self)(out, - deep_copy=True, - dimension_labels=self.dimension_labels, - geometry=self.geometry) - # reverse operand def __radd__(self, other): return self + other @@ -530,7 +511,7 @@ class DataContainer(object): return type(self)(fother ** self.array , dimension_labels=self.dimension_labels, geometry=self.geometry) - elif issubclass(other, DataContainer): + elif issubclass(type(other), DataContainer): if self.check_dimensions(other): return type(self)(other.as_array() ** self.array , dimension_labels=self.dimension_labels, @@ -539,27 +520,59 @@ class DataContainer(object): raise ValueError('Dimensions do not match') # __rpow__ - def sum(self): - return self.as_array().sum() - # in-place arithmetic operators: # (+=, -=, *=, /= , //=, + # must return self + def __iadd__(self, other): - return self + other + if isinstance(other, (int, float)) : + #print ("__iadd__", self.array.shape) + numpy.add(self.array, other, out=self.array) + elif issubclass(type(other), DataContainer): + if self.check_dimensions(other): + numpy.add(self.array, other.array, out=self.array) + else: + raise ValueError('Dimensions do not match') + return self # __iadd__ def __imul__(self, other): - return self * other + if isinstance(other, (int, float)) : + #print ("__imul__", self.array.shape) + #print ("type(self)", type(self)) + #print ("self.array", self.array, other) + arr = self.as_array() + #print ("arr", arr) + numpy.multiply(arr, other, out=arr) + elif issubclass(type(other), DataContainer): + if self.check_dimensions(other): + numpy.multiply(self.array, other.array, out=self.array) + else: + raise ValueError('Dimensions do not match') + return self # __imul__ def __isub__(self, other): - return self - other + if isinstance(other, (int, float)) : + numpy.subtract(self.array, other, out=self.array) + elif issubclass(type(other), DataContainer): + if self.check_dimensions(other): + numpy.subtract(self.array, other.array, out=self.array) + else: + raise ValueError('Dimensions do not match') + return self # __isub__ def __idiv__(self, other): - print ("call __idiv__") - return self / other + if isinstance(other, (int, float)) : + numpy.divide(self.array, other, out=self.array) + elif issubclass(type(other), DataContainer): + if self.check_dimensions(other): + numpy.divide(self.array, other.array, out=self.array) + else: + raise ValueError('Dimensions do not match') + return self # __idiv__ def __str__ (self, representation=False): @@ -578,9 +591,6 @@ class DataContainer(object): dimension_labels=self.dimension_labels, deep_copy=True, geometry=self.geometry ) - def copy(self): - '''alias of clone''' - return self.clone() def get_data_axes_order(self,new_order=None): '''returns the axes label of self as a list @@ -617,13 +627,127 @@ class DataContainer(object): .format(len(self.shape),len(new_order))) - - + def copy(self): + '''alias of clone''' + return self.clone() + + ## binary operations + + def pixel_wise_binary(self,pwop, x2 , *args, out=None, **kwargs): + if out is None: + if isinstance(x2, (int, float, complex)): + out = pwop(self.as_array() , x2 , *args, **kwargs ) + elif issubclass(type(x2) , DataContainer): + out = pwop(self.as_array() , x2.as_array() , *args, **kwargs ) + return type(self)(out, + deep_copy=False, + dimension_labels=self.dimension_labels, + geometry=self.geometry) + + + elif issubclass(type(out), DataContainer) and issubclass(type(x2), DataContainer): + if self.check_dimensions(out) and self.check_dimensions(x2): + pwop(self.as_array(), x2.as_array(), *args, out=out.as_array(), **kwargs ) + return type(self)(out.as_array(), + deep_copy=False, + dimension_labels=self.dimension_labels, + geometry=self.geometry) + else: + raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) + elif issubclass(type(out), DataContainer) and isinstance(x2, (int,float,complex)): + if self.check_dimensions(out): + pwop(self.as_array(), x2, *args, out=out.as_array(), **kwargs ) + return type(self)(out.as_array(), + deep_copy=False, + dimension_labels=self.dimension_labels, + geometry=self.geometry) + else: + raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) + elif issubclass(type(out), numpy.ndarray): + if self.array.shape == out.shape and self.array.dtype == out.dtype: + pwop(self.as_array(), x2 , *args, out=out, **kwargs) + return type(self)(out, + deep_copy=False, + dimension_labels=self.dimension_labels, + geometry=self.geometry) + else: + raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) + + def add(self, other , *args, out=None, **kwargs): + return self.pixel_wise_binary(numpy.add, other, *args, out=out, **kwargs) + + def subtract(self, other , *args, out=None, **kwargs): + return self.pixel_wise_binary(numpy.subtract, other, *args, out=out, **kwargs) + + def multiply(self, other , *args, out=None, **kwargs): + return self.pixel_wise_binary(numpy.multiply, other, *args, out=out, **kwargs) + + def divide(self, other , *args, out=None, **kwargs): + return self.pixel_wise_binary(numpy.divide, other, *args, out=out, **kwargs) + + def power(self, other , *args, out=None, **kwargs): + return self.pixel_wise_binary(numpy.power, other, *args, out=out, **kwargs) + + def maximum(self,x2, out=None): + return self.pixel_wise_binary(numpy.maximum, x2=x2, out=out) + + ## unary operations + + def pixel_wise_unary(self,pwop, *args, out=None, **kwargs): + if out is None: + out = pwop(self.as_array() , *args, **kwargs ) + return type(self)(out, + deep_copy=False, + dimension_labels=self.dimension_labels, + geometry=self.geometry) + elif issubclass(type(out), DataContainer): + if self.check_dimensions(out): + pwop(self.as_array(), *args, out=out.as_array(), **kwargs ) + return type(self)(out.as_array(), + deep_copy=False, + dimension_labels=self.dimension_labels, + geometry=self.geometry) + else: + raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) + elif issubclass(type(out), numpy.ndarray): + if self.array.shape == out.shape and self.array.dtype == out.dtype: + pwop(self.as_array(), *args, out=out, **kwargs) + return type(self)(out, + deep_copy=False, + dimension_labels=self.dimension_labels, + geometry=self.geometry) + else: + raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) + + def abs(self, out=None): + return self.pixel_wise_unary(numpy.abs, out=out) + + def sign(self, out=None): + #out = numpy.sign(self.as_array() ) + #return type(self)(out, + # deep_copy=False, + # dimension_labels=self.dimension_labels, + # geometry=self.geometry) + return self.pixel_wise_unary(numpy.sign , out=out) + + def sqrt(self, out=None): + return self.pixel_wise_unary(numpy.sqrt, out=out) + + #def __abs__(self): + # operation = FM.OPERATION.ABS + # return self.callFieldMath(operation, None, self.mask, self.maskOnValue) + # __abs__ + + ## reductions + def sum(self): + return self.as_array().sum() + + class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' def __init__(self, array = None, - deep_copy=True, + deep_copy=False, dimension_labels=None, **kwargs): @@ -1131,4 +1255,4 @@ if __name__ == '__main__': pixel_num_h=5 , channels=2) sino = AcquisitionData(geometry=sgeometry) sino2 = sino.clone() -
\ No newline at end of file + diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index 570df4b..4a6a383 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -43,22 +43,22 @@ def FISTA(x_init, f=None, g=None, opt=None): # 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'] + 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 - x_old = x_init - y = x_init; + 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) @@ -66,22 +66,44 @@ def FISTA(x_init, f=None, g=None, opt=None): 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() - - u = y - invL*f.grad(y) - - x = g.prox(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 - t_old = t + 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.grad(y) + + x = g.prox(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 @@ -93,7 +115,7 @@ def FISTA(x_init, f=None, g=None, opt=None): #print(it, 'out of', 10, 'iterations', end='\r'); - criter = criter[0:it+1]; + #criter = criter[0:it+1]; timing = numpy.cumsum(timing[0:it+1]); return x, it, timing, criter @@ -127,6 +149,7 @@ def FBPD(x_init, f=None, g=None, h=None, opt=None): 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 / (g.L + 2); @@ -141,6 +164,9 @@ def FBPD(x_init, f=None, g=None, h=None, opt=None): timing = numpy.zeros(max_iter) criter = numpy.zeros(max_iter) + + + # algorithm loop for it in range(0, max_iter): diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py index 4a0a139..490d742 100755 --- a/Wrappers/Python/ccpi/optimisation/funcs.py +++ b/Wrappers/Python/ccpi/optimisation/funcs.py @@ -19,16 +19,31 @@ from ccpi.optimisation.ops import Identity, FiniteDiff2D import numpy +from ccpi.framework import DataContainer +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 Function(object): def __init__(self): self.op = Identity() - def __call__(self,x): return 0 - def grad(self, x): return 0 - def prox(self, x, tau): return x - def gradient(self, x): return self.grad(x) - def proximal(self, x, tau): return self.prox(x, tau) + def __call__(self,x, out=None): return 0 + def grad(self, x): return 0 + def prox(self, x, tau): return x + def gradient(self, x, out=None): return self.grad(x) + def proximal(self, x, tau, out=None): return self.prox(x, tau) class Norm2(Function): @@ -39,10 +54,25 @@ class Norm2(Function): self.gamma = gamma; self.direction = direction; - def __call__(self, x): + def __call__(self, x, out=None): - xx = numpy.sqrt(numpy.sum(numpy.square(x.as_array()), self.direction, + 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 @@ -55,6 +85,30 @@ class Norm2(Function): 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) ) + + class TV2D(Norm2): @@ -81,10 +135,16 @@ class Norm2sq(Function): ''' - def __init__(self,A,b,c=1.0): + def __init__(self,A,b,c=1.0,memopt=False): self.A = A # Should be an operator, default identity self.b = b # Default zero DataSet? self.c = c # Default 1. + self.memopt = memopt + if memopt: + #self.direct_placehold = A.adjoint(b) + self.direct_placehold = A.allocate_direct() + self.adjoint_placehold = A.allocate_adjoint() + # Compute the Lipschitz parameter from the operator. # Initialise to None instead and only call when needed. @@ -93,11 +153,31 @@ class Norm2sq(Function): def grad(self,x): #return 2*self.c*self.A.adjoint( self.A.direct(x) - self.b ) - return 2.0*self.c*self.A.adjoint( self.A.direct(x) - self.b ) + return (2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b ) def __call__(self,x): #return self.c* np.sum(np.square((self.A.direct(x) - self.b).ravel())) - return self.c*( ( (self.A.direct(x)-self.b)**2).sum() ) + #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 + + 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.adjoint_placehold) + self.adjoint_placehold.__isub__( self.b ) + self.A.adjoint(self.adjoint_placehold, out=self.direct_placehold) + self.direct_placehold.__imul__(2.0 * self.c) + # can this be avoided? + out.fill(self.direct_placehold) + else: + return self.grad(x) + class ZeroFun(Function): @@ -111,20 +191,33 @@ class ZeroFun(Function): return 0 def prox(self,x,tau): - return x + return x.copy() + + def proximal(self, x, tau, out=None): + if out is None: + return self.prox(x, tau) + else: + out.fill(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): - # Do nothing self.gamma = gamma self.L = 1 + self.sign_x = None super(Norm1, self).__init__() - def __call__(self,x): - return self.gamma*(x.abs().sum()) + 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() @@ -153,3 +246,21 @@ class IndicatorBox(Function): 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 ) + diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index 668b07e..c6775a8 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -19,55 +19,106 @@ import numpy from scipy.sparse.linalg import svds +from ccpi.framework import DataContainer +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry # Maybe operators need to know what types they take as inputs/outputs # to not just use generic DataContainer class Operator(object): - def direct(self,x): + def direct(self,x, out=None): return x - def adjoint(self,x): + def adjoint(self,x, out=None): return x def size(self): # To be defined for specific class raise NotImplementedError def get_max_sing_val(self): raise NotImplementedError + def allocate_direct(self): + raise NotImplementedError + def allocate_adjoint(self): + raise NotImplementedError class Identity(Operator): def __init__(self): self.s1 = 1.0 super(Identity, self).__init__() - def direct(self,x): - return x + def direct(self,x,out=None): + if out is None: + return x.copy() + else: + out.fill(x) - def adjoint(self,x): - return x + def adjoint(self,x, out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def size(self): + return NotImplemented + + def get_max_sing_val(self): + return self.s1 + +class TomoIdentity(Operator): + def __init__(self, geometry, **kwargs): + self.s1 = 1.0 + self.geometry = geometry + super(TomoIdentity, self).__init__() + + def direct(self,x,out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def adjoint(self,x, out=None): + if out is None: + return x.copy() + else: + out.fill(x) def size(self): return NotImplemented def get_max_sing_val(self): return self.s1 + def allocate_direct(self): + if issubclass(type(self.geometry), ImageGeometry): + return ImageData(geometry=self.geometry) + elif issubclass(type(self.geometry), AcquisitionGeometry): + return AcquisitionData(geometry=self.geometry) + else: + raise ValueError("Wrong geometry type: expected ImageGeometry of AcquisitionGeometry, got ", type(self.geometry)) + def allocate_adjoint(self): + return self.allocate_direct() + + class FiniteDiff2D(Operator): def __init__(self): self.s1 = 8.0 super(FiniteDiff2D, self).__init__() - def direct(self,x): + def direct(self,x, out=None): '''Forward differences with Neumann BC.''' + # FIXME this seems to be working only with numpy arrays d1 = numpy.zeros_like(x.as_array()) d1[:,:-1] = x.as_array()[:,1:] - x.as_array()[:,:-1] d2 = numpy.zeros_like(x.as_array()) d2[:-1,:] = x.as_array()[1:,:] - x.as_array()[:-1,:] d = numpy.stack((d1,d2),0) - return type(x)(d,geometry=x.geometry) + return type(x)(d,False,geometry=x.geometry) - def adjoint(self,x): + def adjoint(self,x, out=None): '''Backward differences, Neumann BC.''' Nrows = x.get_dimension_size('horizontal_x') Ncols = x.get_dimension_size('horizontal_x') @@ -76,12 +127,16 @@ class FiniteDiff2D(Operator): Nchannels = x.get_dimension_size('channel') zer = numpy.zeros((Nrows,1)) xxx = x.as_array()[0,:,:-1] - h = numpy.concatenate((zer,xxx), 1) - numpy.concatenate((xxx,zer), 1) + # + h = numpy.concatenate((zer,xxx), 1) + h -= numpy.concatenate((xxx,zer), 1) zer = numpy.zeros((1,Ncols)) xxx = x.as_array()[1,:-1,:] - v = numpy.concatenate((zer,xxx), 0) - numpy.concatenate((xxx,zer), 0) - return type(x)(h + v,geometry=x.geometry) + # + v = numpy.concatenate((zer,xxx), 0) + v -= numpy.concatenate((xxx,zer), 0) + return type(x)(h + v, False, geometry=x.geometry) def size(self): return NotImplemented @@ -132,7 +187,7 @@ def PowerMethodNonsquareOld(op,numiters): def PowerMethodNonsquare(op,numiters): # Initialise random # Jakob's - #inputsize = op.size()[1] + inputsize , outputsize = op.size() #x0 = ImageContainer(numpy.random.randn(*inputsize) # Edo's #vg = ImageGeometry(voxel_num_x=inputsize[0], @@ -143,7 +198,10 @@ def PowerMethodNonsquare(op,numiters): #print (x0) #x0.fill(numpy.random.randn(*x0.shape)) - x0 = op.create_image_data() + + #x0 = op.create_image_data() + x0 = op.allocate_direct() + x0.fill(numpy.random.randn(*x0.shape)) s = numpy.zeros(numiters) # Loop @@ -162,11 +220,19 @@ class LinearOperatorMatrix(Operator): self.s1 = None # Largest singular value, initially unknown super(LinearOperatorMatrix, self).__init__() - def direct(self,x): - return type(x)(numpy.dot(self.A,x.as_array())) + def direct(self,x, out=None): + if out is None: + return type(x)(numpy.dot(self.A,x.as_array())) + else: + numpy.dot(self.A, x.as_array(), out=out.as_array()) + - def adjoint(self,x): - return type(x)(numpy.dot(self.A.transpose(),x.as_array())) + def adjoint(self,x, out=None): + if out is None: + return type(x)(numpy.dot(self.A.transpose(),x.as_array())) + else: + numpy.dot(self.A.transpose(),x.as_array(), out=out.as_array()) + def size(self): return self.A.shape @@ -178,3 +244,15 @@ class LinearOperatorMatrix(Operator): return self.s1 else: return self.s1 + def allocate_direct(self): + '''allocates the memory to hold the result of adjoint''' + #numpy.dot(self.A.transpose(),x.as_array()) + M_A, N_A = self.A.shape + out = numpy.zeros((N_A,1)) + return DataContainer(out) + def allocate_adjoint(self): + '''allocate the memory to hold the result of direct''' + #numpy.dot(self.A.transpose(),x.as_array()) + M_A, N_A = self.A.shape + out = numpy.zeros((M_A,1)) + return DataContainer(out) diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index 8575b8a..1b3d910 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -20,6 +20,10 @@ requirements: - numpy - scipy - matplotlib + run: + - python + - numpy + - cvxpy about: home: http://www.ccpi.ac.uk diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py index 6a20d05..7df5c07 100755 --- a/Wrappers/Python/conda-recipe/run_test.py +++ b/Wrappers/Python/conda-recipe/run_test.py @@ -1,8 +1,33 @@ import unittest
import numpy
-from ccpi.framework import DataContainer, ImageData, AcquisitionData, \
- ImageGeometry, AcquisitionGeometry
+import numpy as np
+from ccpi.framework import DataContainer
+from ccpi.framework import ImageData
+from ccpi.framework import AcquisitionData
+from ccpi.framework import ImageGeometry
+from ccpi.framework import AcquisitionGeometry
import sys
+from timeit import default_timer as timer
+from ccpi.optimisation.algs import FISTA
+from ccpi.optimisation.algs import FBPD
+from ccpi.optimisation.funcs import Norm2sq
+from ccpi.optimisation.funcs import ZeroFun
+from ccpi.optimisation.funcs import Norm1
+from ccpi.optimisation.funcs import TV2D
+
+from ccpi.optimisation.ops import LinearOperatorMatrix
+from ccpi.optimisation.ops import TomoIdentity
+
+from cvxpy import *
+
+
+def aid(x):
+ # This function returns the memory
+ # block address of an array.
+ return x.__array_interface__['data'][0]
+
+def dt (steps):
+ return steps[-1] - steps[-2]
class TestDataContainer(unittest.TestCase):
@@ -21,6 +46,282 @@ class TestDataContainer(unittest.TestCase): self.assertEqual(sys.getrefcount(a),3)
self.assertEqual(ds.dimension_labels , {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'})
+ def testGb_creation_nocopy(self):
+ X,Y,Z = 512,512,512
+ steps = [timer()]
+ a = numpy.ones((X,Y,Z), dtype='float32')
+ steps.append(timer())
+ t0 = dt(steps)
+ print("test clone")
+ #print("a refcount " , sys.getrefcount(a))
+ ds = DataContainer(a, False, ['X', 'Y','Z'])
+ #print("a refcount " , sys.getrefcount(a))
+ self.assertEqual(sys.getrefcount(a),3)
+ ds1 = ds.copy()
+ self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array()))
+ ds1 = ds.clone()
+ self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array()))
+
+
+ def testInlineAlgebra(self):
+ print ("Test Inline Algebra")
+ X,Y,Z = 1024,512,512
+ steps = [timer()]
+ a = numpy.ones((X,Y,Z), dtype='float32')
+ steps.append(timer())
+ t0 = dt(steps)
+ print(t0)
+ #print("a refcount " , sys.getrefcount(a))
+ ds = DataContainer(a, False, ['X', 'Y','Z'])
+ #ds.__iadd__( 2 )
+ ds += 2
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],3.)
+ #ds.__isub__( 2 )
+ ds -= 2
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],1.)
+ #ds.__imul__( 2 )
+ ds *= 2
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],2.)
+ #ds.__idiv__( 2 )
+ ds /= 2
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],1.)
+
+ ds1 = ds.copy()
+ #ds1.__iadd__( 1 )
+ ds1 += 1
+ #ds.__iadd__( ds1 )
+ ds += ds1
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],3.)
+ #ds.__isub__( ds1 )
+ ds -= ds1
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],1.)
+ #ds.__imul__( ds1 )
+ ds *= ds1
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],2.)
+ #ds.__idiv__( ds1 )
+ ds /= ds1
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],1.)
+
+
+ def test_unary_operations(self):
+ print ("Test unary operations")
+ X,Y,Z = 1024,512,512
+ steps = [timer()]
+ a = -numpy.ones((X,Y,Z), dtype='float32')
+ steps.append(timer())
+ t0 = dt(steps)
+ print(t0)
+ #print("a refcount " , sys.getrefcount(a))
+ ds = DataContainer(a, False, ['X', 'Y','Z'])
+
+ ds.sign(out=ds)
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],-1.)
+
+ ds.abs(out=ds)
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],1.)
+
+ ds.__imul__( 2 )
+ ds.sqrt(out=ds)
+ steps.append(timer())
+ print(dt(steps))
+ self.assertEqual(ds.as_array()[0][0][0],numpy.sqrt(2., dtype='float32'))
+
+
+
+ def test_binary_operations(self):
+ self.binary_add()
+ self.binary_subtract()
+ self.binary_multiply()
+ self.binary_divide()
+
+ def binary_add(self):
+ print ("Test binary add")
+ X,Y,Z = 512,512,512
+ steps = [timer()]
+ a = numpy.ones((X,Y,Z), dtype='float32')
+ steps.append(timer())
+ t0 = dt(steps)
+
+ #print("a refcount " , sys.getrefcount(a))
+ ds = DataContainer(a, False, ['X', 'Y','Z'])
+ ds1 = ds.copy()
+
+ steps.append(timer())
+ ds.add(ds1, out=ds)
+ steps.append(timer())
+ t1 = dt(steps)
+ print("ds.add(ds1, out=ds)",dt(steps))
+ steps.append(timer())
+ ds2 = ds.add(ds1)
+ steps.append(timer())
+ t2 = dt(steps)
+ print("ds2 = ds.add(ds1)",dt(steps))
+
+ self.assertLess(t1,t2)
+ self.assertEqual(ds.as_array()[0][0][0] , 2.)
+
+ ds0 = ds
+ ds0.add(2,out=ds0)
+ steps.append(timer())
+ print ("ds0.add(2,out=ds0)", dt(steps), 3 , ds0.as_array()[0][0][0])
+ self.assertEqual(4., ds0.as_array()[0][0][0])
+
+ dt1 = dt(steps)
+ ds3 = ds0.add(2)
+ steps.append(timer())
+ print ("ds3 = ds0.add(2)", dt(steps), 5 , ds3.as_array()[0][0][0])
+ dt2 = dt(steps)
+ self.assertLess(dt1,dt2)
+
+ def binary_subtract(self):
+ print ("Test binary subtract")
+ X,Y,Z = 512,512,512
+ steps = [timer()]
+ a = numpy.ones((X,Y,Z), dtype='float32')
+ steps.append(timer())
+ t0 = dt(steps)
+
+ #print("a refcount " , sys.getrefcount(a))
+ ds = DataContainer(a, False, ['X', 'Y','Z'])
+ ds1 = ds.copy()
+
+ steps.append(timer())
+ ds.subtract(ds1, out=ds)
+ steps.append(timer())
+ t1 = dt(steps)
+ print("ds.subtract(ds1, out=ds)",dt(steps))
+ self.assertEqual(0., ds.as_array()[0][0][0])
+
+ steps.append(timer())
+ ds2 = ds.subtract(ds1)
+ self.assertEqual(-1., ds2.as_array()[0][0][0])
+
+ steps.append(timer())
+ t2 = dt(steps)
+ print("ds2 = ds.subtract(ds1)",dt(steps))
+
+ self.assertLess(t1,t2)
+
+ del ds1
+ ds0 = ds.copy()
+ steps.append(timer())
+ ds0.subtract(2,out=ds0)
+ #ds0.__isub__( 2 )
+ steps.append(timer())
+ print ("ds0.subtract(2,out=ds0)", dt(steps), -2. , ds0.as_array()[0][0][0])
+ self.assertEqual(-2., ds0.as_array()[0][0][0])
+
+ dt1 = dt(steps)
+ ds3 = ds0.subtract(2)
+ steps.append(timer())
+ print ("ds3 = ds0.subtract(2)", dt(steps), 0. , ds3.as_array()[0][0][0])
+ dt2 = dt(steps)
+ self.assertLess(dt1,dt2)
+ self.assertEqual(-2., ds0.as_array()[0][0][0])
+ self.assertEqual(-4., ds3.as_array()[0][0][0])
+
+ def binary_multiply(self):
+ print ("Test binary multiply")
+ X,Y,Z = 1024,512,512
+ steps = [timer()]
+ a = numpy.ones((X,Y,Z), dtype='float32')
+ steps.append(timer())
+ t0 = dt(steps)
+
+ #print("a refcount " , sys.getrefcount(a))
+ ds = DataContainer(a, False, ['X', 'Y','Z'])
+ ds1 = ds.copy()
+
+ steps.append(timer())
+ ds.multiply(ds1, out=ds)
+ steps.append(timer())
+ t1 = dt(steps)
+ print("ds.multiply(ds1, out=ds)",dt(steps))
+ steps.append(timer())
+ ds2 = ds.multiply(ds1)
+ steps.append(timer())
+ t2 = dt(steps)
+ print("ds2 = ds.multiply(ds1)",dt(steps))
+
+ self.assertLess(t1,t2)
+
+ ds0 = ds
+ ds0.multiply(2,out=ds0)
+ steps.append(timer())
+ print ("ds0.multiply(2,out=ds0)", dt(steps), 2. , ds0.as_array()[0][0][0])
+ self.assertEqual(2., ds0.as_array()[0][0][0])
+
+ dt1 = dt(steps)
+ ds3 = ds0.multiply(2)
+ steps.append(timer())
+ print ("ds3 = ds0.multiply(2)", dt(steps), 4. , ds3.as_array()[0][0][0])
+ dt2 = dt(steps)
+ self.assertLess(dt1,dt2)
+ self.assertEqual(4., ds3.as_array()[0][0][0])
+ self.assertEqual(2., ds.as_array()[0][0][0])
+
+ def binary_divide(self):
+ print ("Test binary divide")
+ X,Y,Z = 1024,512,512
+ steps = [timer()]
+ a = numpy.ones((X,Y,Z), dtype='float32')
+ steps.append(timer())
+ t0 = dt(steps)
+
+ #print("a refcount " , sys.getrefcount(a))
+ ds = DataContainer(a, False, ['X', 'Y','Z'])
+ ds1 = ds.copy()
+
+ steps.append(timer())
+ ds.divide(ds1, out=ds)
+ steps.append(timer())
+ t1 = dt(steps)
+ print("ds.divide(ds1, out=ds)",dt(steps))
+ steps.append(timer())
+ ds2 = ds.divide(ds1)
+ steps.append(timer())
+ t2 = dt(steps)
+ print("ds2 = ds.divide(ds1)",dt(steps))
+
+ self.assertLess(t1,t2)
+ self.assertEqual(ds.as_array()[0][0][0] , 1.)
+
+ ds0 = ds
+ ds0.divide(2,out=ds0)
+ steps.append(timer())
+ print ("ds0.divide(2,out=ds0)", dt(steps), 0.5 , ds0.as_array()[0][0][0])
+ self.assertEqual(0.5, ds0.as_array()[0][0][0])
+
+ dt1 = dt(steps)
+ ds3 = ds0.divide(2)
+ steps.append(timer())
+ print ("ds3 = ds0.divide(2)", dt(steps), 0.25 , ds3.as_array()[0][0][0])
+ dt2 = dt(steps)
+ self.assertLess(dt1,dt2)
+ self.assertEqual(.25, ds3.as_array()[0][0][0])
+ self.assertEqual(.5, ds.as_array()[0][0][0])
+
+
def test_creation_copy(self):
shape = (2,3,4,5)
size = shape[0]
@@ -147,6 +448,235 @@ class TestDataContainer(unittest.TestCase): res = False
print (err)
self.assertTrue(res)
+
+
+ def assertNumpyArrayAlmostEqual(self, first, second, decimal=6):
+ res = True
+ try:
+ numpy.testing.assert_array_almost_equal(first, second, decimal)
+ except AssertionError as err:
+ res = False
+ print (err)
+ self.assertTrue(res)
+
+class TestAlgorithms(unittest.TestCase):
+ def assertNumpyArrayEqual(self, first, second):
+ res = True
+ try:
+ numpy.testing.assert_array_equal(first, second)
+ except AssertionError as err:
+ res = False
+ print (err)
+ self.assertTrue(res)
+
+
+ def assertNumpyArrayAlmostEqual(self, first, second, decimal=6):
+ res = True
+ try:
+ numpy.testing.assert_array_almost_equal(first, second, decimal)
+ except AssertionError as err:
+ res = False
+ print (err)
+ self.assertTrue(res)
+ def test_FISTA(self):
+ # Problem data.
+ m = 30
+ n = 20
+ np.random.seed(1)
+ Amat = np.random.randn(m, n)
+ A = LinearOperatorMatrix(Amat)
+ bmat = np.random.randn(m)
+ bmat.shape = (bmat.shape[0],1)
+
+ # A = Identity()
+ # Change n to equal to m.
+
+ b = DataContainer(bmat)
+
+ # Regularization parameter
+ lam = 10
+ opt = {'memopt':True}
+ # Create object instances with the test data A and b.
+ f = Norm2sq(A,b,c=0.5, memopt=True)
+ g0 = ZeroFun()
+
+ # Initial guess
+ x_init = DataContainer(np.zeros((n,1)))
+
+ f.grad(x_init)
+
+ # Run FISTA for least squares plus zero function.
+ x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0 , opt=opt)
+
+ # Print solution and final objective/criterion value for comparison
+ print("FISTA least squares plus zero function solution and objective value:")
+ print(x_fista0.array)
+ print(criter0[-1])
+
+ # Compare to CVXPY
+
+ # Construct the problem.
+ x0 = Variable(n)
+ objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]) )
+ prob0 = Problem(objective0)
+
+ # The optimal objective is returned by prob.solve().
+ result0 = prob0.solve(verbose=False,solver=SCS,eps=1e-9)
+
+ # The optimal solution for x is stored in x.value and optimal objective value
+ # is in result as well as in objective.value
+ print("CVXPY least squares plus zero function solution and objective value:")
+ print(x0.value)
+ print(objective0.value)
+ self.assertNumpyArrayAlmostEqual(
+ numpy.squeeze(x_fista0.array),x0.value,6)
+ # Create 1-norm object instance
+ g1 = Norm1(lam)
+
+ g1(x_init)
+ g1.prox(x_init,0.02)
+
+ # Combine with least squares and solve using generic FISTA implementation
+ x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1,opt=opt)
+
+ # Print for comparison
+ print("FISTA least squares plus 1-norm solution and objective value:")
+ print(x_fista1)
+ print(criter1[-1])
+
+ # Compare to CVXPY
+
+ # Construct the problem.
+ x1 = Variable(n)
+ objective1 = Minimize(0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1,1) )
+ prob1 = Problem(objective1)
+
+ # The optimal objective is returned by prob.solve().
+ result1 = prob1.solve(verbose=False,solver=SCS,eps=1e-9)
+
+ # The optimal solution for x is stored in x.value and optimal objective value
+ # is in result as well as in objective.value
+ print("CVXPY least squares plus 1-norm solution and objective value:")
+ print(x1.value)
+ print(objective1.value)
+
+ self.assertNumpyArrayAlmostEqual(
+ numpy.squeeze(x_fista1.array),x1.value,6)
+
+ # Now try another algorithm FBPD for same problem:
+ x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, None, f, g1)
+ print(x_fbpd1)
+ print(criterfbpd1[-1])
+
+ self.assertNumpyArrayAlmostEqual(
+ numpy.squeeze(x_fbpd1.array),x1.value,6)
+ # Plot criterion curve to see both FISTA and FBPD converge to same value.
+ # Note that FISTA is very efficient for 1-norm minimization so it beats
+ # FBPD in this test by a lot. But FBPD can handle a larger class of problems
+ # than FISTA can.
+
+
+ # Now try 1-norm and TV denoising with FBPD, first 1-norm.
+
+ # Set up phantom size NxN by creating ImageGeometry, initialising the
+ # ImageData object with this geometry and empty array and finally put some
+ # data into its array, and display as image.
+ N = 64
+ ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+ Phantom = ImageData(geometry=ig)
+
+ x = Phantom.as_array()
+ x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+ x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+
+
+ # Identity operator for denoising
+ I = TomoIdentity(ig)
+
+ # Data and add noise
+ y = I.direct(Phantom)
+ y.array = y.array + 0.1*np.random.randn(N, N)
+
+
+ # Data fidelity term
+ f_denoise = Norm2sq(I,y,c=0.5,memopt=True)
+
+ # 1-norm regulariser
+ lam1_denoise = 1.0
+ g1_denoise = Norm1(lam1_denoise)
+
+ # Initial guess
+ x_init_denoise = ImageData(np.zeros((N,N)))
+
+ # Combine with least squares and solve using generic FISTA implementation
+ x_fista1_denoise, it1_denoise, timing1_denoise, criter1_denoise = FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt)
+
+ print(x_fista1_denoise)
+ print(criter1_denoise[-1])
+
+
+ # Now denoise LS + 1-norm with FBPD
+ x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise, criterfbpd1_denoise = FBPD(x_init_denoise, None, f_denoise, g1_denoise)
+ print(x_fbpd1_denoise)
+ print(criterfbpd1_denoise[-1])
+
+
+ # Compare to CVXPY
+
+ # Construct the problem.
+ x1_denoise = Variable(N**2,1)
+ objective1_denoise = Minimize(0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise,1) )
+ prob1_denoise = Problem(objective1_denoise)
+
+ # The optimal objective is returned by prob.solve().
+ result1_denoise = prob1_denoise.solve(verbose=False,solver=SCS,eps=1e-12)
+
+ # The optimal solution for x is stored in x.value and optimal objective value
+ # is in result as well as in objective.value
+ print("CVXPY least squares plus 1-norm solution and objective value:")
+ print(x1_denoise.value)
+ print(objective1_denoise.value)
+ self.assertNumpyArrayAlmostEqual(
+ x_fista1_denoise.array.flatten(),x1_denoise.value,5)
+
+ self.assertNumpyArrayAlmostEqual(
+ x_fbpd1_denoise.array.flatten(),x1_denoise.value,5)
+ x1_cvx = x1_denoise.value
+ x1_cvx.shape = (N,N)
+
+
+ # Now TV with FBPD
+ lam_tv = 0.1
+ gtv = TV2D(lam_tv)
+ gtv(gtv.op.direct(x_init_denoise))
+
+ opt_tv = {'tol': 1e-4, 'iter': 10000}
+
+ x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, criterfbpdtv_denoise = FBPD(x_init_denoise, None, f_denoise, gtv,opt=opt_tv)
+ print(x_fbpdtv_denoise)
+ print(criterfbpdtv_denoise[-1])
+
+
+
+ # Compare to CVXPY
+
+ # Construct the problem.
+ xtv_denoise = Variable((N,N))
+ objectivetv_denoise = Minimize(0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise) )
+ probtv_denoise = Problem(objectivetv_denoise)
+
+ # The optimal objective is returned by prob.solve().
+ resulttv_denoise = probtv_denoise.solve(verbose=False,solver=SCS,eps=1e-12)
+
+ # The optimal solution for x is stored in x.value and optimal objective value
+ # is in result as well as in objective.value
+ print("CVXPY least squares plus 1-norm solution and objective value:")
+ print(xtv_denoise.value)
+ print(objectivetv_denoise.value)
+
+ self.assertNumpyArrayAlmostEqual(
+ x_fbpdtv_denoise.as_array(),xtv_denoise.value,1)
+
# =============================================================================
# def test_upper(self):
# self.assertEqual('foo'.upper(), 'FOO')
@@ -164,4 +694,4 @@ class TestDataContainer(unittest.TestCase): # =============================================================================
if __name__ == '__main__':
- unittest.main()
\ No newline at end of file + unittest.main()
diff --git a/Wrappers/Python/wip/demo_compare_cvx.py b/Wrappers/Python/wip/demo_compare_cvx.py index cbfe50e..2df11c1 100644 --- a/Wrappers/Python/wip/demo_compare_cvx.py +++ b/Wrappers/Python/wip/demo_compare_cvx.py @@ -3,7 +3,7 @@ from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataCo from ccpi.optimisation.algs import FISTA, FBPD, CGLS from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D -from ccpi.optimisation.ops import LinearOperatorMatrix, Identity +from ccpi.optimisation.ops import LinearOperatorMatrix, TomoIdentity # Requires CVXPY, see http://www.cvxpy.org/ # CVXPY can be installed in anaconda using @@ -33,9 +33,9 @@ b = DataContainer(bmat) # Regularization parameter lam = 10 - +opt = {'memopt':True} # Create object instances with the test data A and b. -f = Norm2sq(A,b,c=0.5) +f = Norm2sq(A,b,c=0.5, memopt=True) g0 = ZeroFun() # Initial guess @@ -44,7 +44,7 @@ x_init = DataContainer(np.zeros((n,1))) f.grad(x_init) # Run FISTA for least squares plus zero function. -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0) +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0 , opt=opt) # Print solution and final objective/criterion value for comparison print("FISTA least squares plus zero function solution and objective value:") @@ -56,7 +56,7 @@ if use_cvxpy: # Construct the problem. x0 = Variable(n) - objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat) ) + objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]) ) prob0 = Problem(objective0) # The optimal objective is returned by prob.solve(). @@ -83,7 +83,7 @@ g1(x_init) g1.prox(x_init,0.02) # Combine with least squares and solve using generic FISTA implementation -x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1) +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1,opt=opt) # Print for comparison print("FISTA least squares plus 1-norm solution and objective value:") @@ -95,7 +95,7 @@ if use_cvxpy: # Construct the problem. x1 = Variable(n) - objective1 = Minimize(0.5*sum_squares(Amat*x1 - bmat) + lam*norm(x1,1) ) + objective1 = Minimize(0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1,1) ) prob1 = Problem(objective1) # The optimal objective is returned by prob.solve(). @@ -147,7 +147,7 @@ plt.title('Phantom image') plt.show() # Identity operator for denoising -I = Identity() +I = TomoIdentity(ig) # Data and add noise y = I.direct(Phantom) @@ -158,7 +158,7 @@ plt.title('Noisy image') plt.show() # Data fidelity term -f_denoise = Norm2sq(I,y,c=0.5) +f_denoise = Norm2sq(I,y,c=0.5,memopt=True) # 1-norm regulariser lam1_denoise = 1.0 @@ -168,7 +168,7 @@ g1_denoise = Norm1(lam1_denoise) x_init_denoise = ImageData(np.zeros((N,N))) # Combine with least squares and solve using generic FISTA implementation -x_fista1_denoise, it1_denoise, timing1_denoise, criter1_denoise = FISTA(x_init_denoise, f_denoise, g1_denoise) +x_fista1_denoise, it1_denoise, timing1_denoise, criter1_denoise = FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) print(x_fista1_denoise) print(criter1_denoise[-1]) @@ -229,7 +229,8 @@ if use_cvxpy: # Compare to CVXPY # Construct the problem. - xtv_denoise = Variable(N,N) + xtv_denoise = Variable((N,N)) + print (xtv_denoise.shape) objectivetv_denoise = Minimize(0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise) ) probtv_denoise = Problem(objectivetv_denoise) @@ -247,4 +248,4 @@ plt.title('CVX TV') plt.show() plt.loglog([0,opt_tv['iter']], [objectivetv_denoise.value,objectivetv_denoise.value], label='CVX TV') -plt.loglog(criterfbpdtv_denoise, label='FBPD TV')
\ No newline at end of file +plt.loglog(criterfbpdtv_denoise, label='FBPD TV') diff --git a/Wrappers/Python/wip/demo_memhandle.py b/Wrappers/Python/wip/demo_memhandle.py new file mode 100755 index 0000000..db48d73 --- /dev/null +++ b/Wrappers/Python/wip/demo_memhandle.py @@ -0,0 +1,193 @@ + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer +from ccpi.optimisation.algs import FISTA, FBPD, CGLS +from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D + +from ccpi.optimisation.ops import LinearOperatorMatrix, Identity +from ccpi.optimisation.ops import TomoIdentity + +# Requires CVXPY, see http://www.cvxpy.org/ +# CVXPY can be installed in anaconda using +# conda install -c cvxgrp cvxpy libgcc + + +import numpy as np +import matplotlib.pyplot as plt + +# Problem data. +m = 30 +n = 20 +np.random.seed(1) +Amat = np.random.randn(m, n) +A = LinearOperatorMatrix(Amat) +bmat = np.random.randn(m) +bmat.shape = (bmat.shape[0],1) + + + +# A = Identity() +# Change n to equal to m. + +b = DataContainer(bmat) + +# Regularization parameter +lam = 10 + +# Create object instances with the test data A and b. +f = Norm2sq(A,b,c=0.5, memopt=True) +g0 = ZeroFun() + +# Initial guess +x_init = DataContainer(np.zeros((n,1))) + +f.grad(x_init) +opt = {'memopt': True} +# Run FISTA for least squares plus zero function. +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0) +x_fista0_m, it0_m, timing0_m, criter0_m = FISTA(x_init, f, g0, opt=opt) + +iternum = [i for i in range(len(criter0))] +# Print solution and final objective/criterion value for comparison +print("FISTA least squares plus zero function solution and objective value:") +print(x_fista0.array) +print(criter0[-1]) + + +# Plot criterion curve to see FISTA converge to same value as CVX. +#iternum = np.arange(1,1001) +plt.figure() +plt.loglog(iternum,criter0,label='FISTA LS') +plt.loglog(iternum,criter0_m,label='FISTA LS memopt') +plt.legend() +plt.show() +#%% +# Create 1-norm object instance +g1 = Norm1(lam) + +g1(x_init) +g1.prox(x_init,0.02) + +# Combine with least squares and solve using generic FISTA implementation +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1) +x_fista1_m, it1_m, timing1_m, criter1_m = FISTA(x_init, f, g1, opt=opt) +iternum = [i for i in range(len(criter1))] +# Print for comparison +print("FISTA least squares plus 1-norm solution and objective value:") +print(x_fista1) +print(criter1[-1]) + + +# Now try another algorithm FBPD for same problem: +x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, None, f, g1) +iternum = [i for i in range(len(criterfbpd1))] +print(x_fbpd1) +print(criterfbpd1[-1]) + +# Plot criterion curve to see both FISTA and FBPD converge to same value. +# Note that FISTA is very efficient for 1-norm minimization so it beats +# FBPD in this test by a lot. But FBPD can handle a larger class of problems +# than FISTA can. +plt.figure() +plt.loglog(iternum,criter1,label='FISTA LS+1') +plt.loglog(iternum,criter1_m,label='FISTA LS+1 memopt') +plt.legend() +plt.show() + +plt.figure() +plt.loglog(iternum,criter1,label='FISTA LS+1') +plt.loglog(iternum,criterfbpd1,label='FBPD LS+1') +plt.legend() +plt.show() +#%% +# Now try 1-norm and TV denoising with FBPD, first 1-norm. + +# Set up phantom size NxN by creating ImageGeometry, initialising the +# ImageData object with this geometry and empty array and finally put some +# data into its array, and display as image. +N = 1000 +ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = ImageData(geometry=ig) + +x = Phantom.as_array() +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 + +plt.imshow(x) +plt.title('Phantom image') +plt.show() + +# Identity operator for denoising +I = TomoIdentity(ig) + +# Data and add noise +y = I.direct(Phantom) +y.array += 0.1*np.random.randn(N, N) + +plt.figure() +plt.imshow(y.array) +plt.title('Noisy image') +plt.show() + +# Data fidelity term +f_denoise = Norm2sq(I,y,c=0.5, memopt=True) + +# 1-norm regulariser +lam1_denoise = 1.0 +g1_denoise = Norm1(lam1_denoise) + +# Initial guess +x_init_denoise = ImageData(np.zeros((N,N))) +opt = {'memopt': False, 'iter' : 50} +# Combine with least squares and solve using generic FISTA implementation +print ("no memopt") +x_fista1_denoise, it1_denoise, timing1_denoise, \ + criter1_denoise = FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) +opt = {'memopt': True, 'iter' : 50} +print ("yes memopt") +x_fista1_denoise_m, it1_denoise_m, timing1_denoise_m, \ + criter1_denoise_m = FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) + +print(x_fista1_denoise) +print(criter1_denoise[-1]) + +plt.figure() +plt.imshow(x_fista1_denoise.as_array()) +plt.title('FISTA LS+1') +plt.show() + +plt.figure() +plt.imshow(x_fista1_denoise_m.as_array()) +plt.title('FISTA LS+1 memopt') +plt.show() + +plt.figure() +plt.loglog(iternum,criter1_denoise,label='FISTA LS+1') +plt.loglog(iternum,criter1_denoise_m,label='FISTA LS+1 memopt') +plt.legend() +plt.show() +#%% +# Now denoise LS + 1-norm with FBPD +x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise, criterfbpd1_denoise = FBPD(x_init_denoise, None, f_denoise, g1_denoise) +print(x_fbpd1_denoise) +print(criterfbpd1_denoise[-1]) + +plt.figure() +plt.imshow(x_fbpd1_denoise.as_array()) +plt.title('FBPD LS+1') +plt.show() + + +# Now TV with FBPD +lam_tv = 0.1 +gtv = TV2D(lam_tv) +gtv(gtv.op.direct(x_init_denoise)) + +opt_tv = {'tol': 1e-4, 'iter': 10000} + +x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, criterfbpdtv_denoise = FBPD(x_init_denoise, None, f_denoise, gtv,opt=opt_tv) +print(x_fbpdtv_denoise) +print(criterfbpdtv_denoise[-1]) + +plt.imshow(x_fbpdtv_denoise.as_array()) +plt.title('FBPD TV') +plt.show() |