diff options
69 files changed, 8520 insertions, 566 deletions
diff --git a/Wrappers/Python/build/lib/ccpi/contrib/optimisation/__init__.py b/Wrappers/Python/build/lib/ccpi/contrib/optimisation/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/contrib/optimisation/__init__.py diff --git a/Wrappers/Python/build/lib/ccpi/contrib/optimisation/algorithms/__init__.py b/Wrappers/Python/build/lib/ccpi/contrib/optimisation/algorithms/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/contrib/optimisation/algorithms/__init__.py diff --git a/Wrappers/Python/build/lib/ccpi/contrib/optimisation/algorithms/spdhg.py b/Wrappers/Python/build/lib/ccpi/contrib/optimisation/algorithms/spdhg.py new file mode 100644 index 0000000..263a7cd --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/contrib/optimisation/algorithms/spdhg.py @@ -0,0 +1,338 @@ +# Copyright 2018 Matthias Ehrhardt, 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.
+
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy
+
+from ccpi.optimisation.funcs import Function
+from ccpi.framework import ImageData
+from ccpi.framework import AcquisitionData
+
+
+class spdhg():
+ """Computes a saddle point with a stochastic PDHG.
+
+ This means, a solution (x*, y*), y* = (y*_1, ..., y*_n) such that
+
+ (x*, y*) in arg min_x max_y sum_i=1^n <y_i, A_i> - f*[i](y_i) + g(x)
+
+ where g : X -> IR_infty and f[i] : Y[i] -> IR_infty are convex, l.s.c. and
+ proper functionals. For this algorithm, they all may be non-smooth and no
+ strong convexity is assumed.
+
+ Parameters
+ ----------
+ f : list of functions
+ Functionals Y[i] -> IR_infty that all have a convex conjugate with a
+ proximal operator, i.e.
+ f[i].convex_conj.prox(sigma[i]) : Y[i] -> Y[i].
+ g : function
+ Functional X -> IR_infty that has a proximal operator, i.e.
+ g.prox(tau) : X -> X.
+ A : list of functions
+ Operators A[i] : X -> Y[i] that possess adjoints: A[i].adjoint
+ x : primal variable, optional
+ By default equals 0.
+ y : dual variable, optional
+ Part of a product space. By default equals 0.
+ z : variable, optional
+ Adjoint of dual variable, z = A^* y. By default equals 0 if y = 0.
+ tau : scalar / vector / matrix, optional
+ Step size for primal variable. Note that the proximal operator of g
+ has to be well-defined for this input.
+ sigma : scalar, optional
+ Scalar / vector / matrix used as step size for dual variable. Note that
+ the proximal operator related to f (see above) has to be well-defined
+ for this input.
+ prob : list of scalars, optional
+ Probabilities prob[i] that a subset i is selected in each iteration.
+ If fun_select is not given, then the sum of all probabilities must
+ equal 1.
+ A_norms : list of scalars, optional
+ Norms of the operators in A. Can be used to determine the step sizes
+ tau and sigma and the probabilities prob.
+ fun_select : function, optional
+ Function that selects blocks at every iteration IN -> {1,...,n}. By
+ default this is serial sampling, fun_select(k) selects an index
+ i \in {1,...,n} with probability prob[i].
+
+ References
+ ----------
+ [CERS2018] A. Chambolle, M. J. Ehrhardt, P. Richtarik and C.-B. Schoenlieb,
+ *Stochastic Primal-Dual Hybrid Gradient Algorithm with Arbitrary Sampling
+ and Imaging Applications*. SIAM Journal on Optimization, 28(4), 2783-2808
+ (2018) http://doi.org/10.1007/s10851-010-0251-1
+
+ [E+2017] M. J. Ehrhardt, P. J. Markiewicz, P. Richtarik, J. Schott,
+ A. Chambolle and C.-B. Schoenlieb, *Faster PET reconstruction with a
+ stochastic primal-dual hybrid gradient method*. Wavelets and Sparsity XVII,
+ 58 (2017) http://doi.org/10.1117/12.2272946.
+
+ [EMS2018] M. J. Ehrhardt, P. J. Markiewicz and C.-B. Schoenlieb, *Faster
+ PET Reconstruction with Non-Smooth Priors by Randomization and
+ Preconditioning*. (2018) ArXiv: http://arxiv.org/abs/1808.07150
+ """
+
+ def __init__(self, f, g, A, x=None, y=None, z=None, tau=None, sigma=None,
+ prob=None, A_norms=None, fun_select=None):
+ # fun_select is optional and by default performs serial sampling
+
+ if x is None:
+ x = A[0].allocate_direct(0)
+
+ if y is None:
+ if z is not None:
+ raise ValueError('y and z have to be defaulted together')
+
+ y = [Ai.allocate_adjoint(0) for Ai in A]
+ z = 0 * x.copy()
+
+ else:
+ if z is None:
+ raise ValueError('y and z have to be defaulted together')
+
+ if A_norms is not None:
+ if tau is not None or sigma is not None or prob is not None:
+ raise ValueError('Either A_norms or (tau, sigma, prob) must '
+ 'be given')
+
+ tau = 1 / sum(A_norms)
+ sigma = [1 / nA for nA in A_norms]
+ prob = [nA / sum(A_norms) for nA in A_norms]
+
+ #uniform prob, needs different sigma and tau
+ #n = len(A)
+ #prob = [1./n] * n
+
+ if fun_select is None:
+ if prob is None:
+ raise ValueError('prob was not determined')
+
+ def fun_select(k):
+ return [int(numpy.random.choice(len(A), 1, p=prob))]
+
+ self.iter = 0
+ self.x = x
+
+ self.y = y
+ self.z = z
+
+ self.f = f
+ self.g = g
+ self.A = A
+ self.tau = tau
+ self.sigma = sigma
+ self.prob = prob
+ self.fun_select = fun_select
+
+ # Initialize variables
+ self.z_relax = z.copy()
+ self.tmp = self.x.copy()
+
+ def update(self):
+ # select block
+ selected = self.fun_select(self.iter)
+
+ # update primal variable
+ #tmp = (self.x - self.tau * self.z_relax).as_array()
+ #self.x.fill(self.g.prox(tmp, self.tau))
+ self.tmp = - self.tau * self.z_relax
+ self.tmp += self.x
+ self.x = self.g.prox(self.tmp, self.tau)
+
+ # update dual variable and z, z_relax
+ self.z_relax = self.z.copy()
+ for i in selected:
+ # save old yi
+ y_old = self.y[i].copy()
+
+ # y[i]= prox(tmp)
+ tmp = y_old + self.sigma[i] * self.A[i].direct(self.x)
+ self.y[i] = self.f[i].convex_conj.prox(tmp, self.sigma[i])
+
+ # update adjoint of dual variable
+ dz = self.A[i].adjoint(self.y[i] - y_old)
+ self.z += dz
+
+ # compute extrapolation
+ self.z_relax += (1 + 1 / self.prob[i]) * dz
+
+ self.iter += 1
+
+
+## Functions
+
+class KullbackLeibler(Function):
+ def __init__(self, data, background):
+ self.data = data
+ self.background = background
+ self.__offset = None
+
+ def __call__(self, x):
+ """Return the KL-diveregnce in the point ``x``.
+
+ If any components of ``x`` is non-positive, the value is positive
+ infinity.
+
+ Needs one extra array of memory of the size of `prior`.
+ """
+
+ # define short variable names
+ y = self.data
+ r = self.background
+
+ # Compute
+ # sum(x + r - y + y * log(y / (x + r)))
+ # = sum(x - y * log(x + r)) + self.offset
+ # Assume that
+ # x + r > 0
+
+ # sum the result up
+ obj = numpy.sum(x - y * numpy.log(x + r)) + self.offset()
+
+ if numpy.isnan(obj):
+ # In this case, some element was less than or equal to zero
+ return numpy.inf
+ else:
+ return obj
+
+ @property
+ def convex_conj(self):
+ """The convex conjugate functional of the KL-functional."""
+ return KullbackLeiblerConvexConjugate(self.data, self.background)
+
+ def offset(self):
+ """The offset which is independent of the unknown."""
+
+ if self.__offset is None:
+ tmp = self.domain.element()
+
+ # define short variable names
+ y = self.data
+ r = self.background
+
+ tmp = self.domain.element(numpy.maximum(y, 1))
+ tmp = r - y + y * numpy.log(tmp)
+
+ # sum the result up
+ self.__offset = numpy.sum(tmp)
+
+ return self.__offset
+
+# def __repr__(self):
+# """to be added???"""
+# """Return ``repr(self)``."""
+ # return '{}({!r}, {!r}, {!r})'.format(self.__class__.__name__,
+ ## self.domain, self.data,
+ # self.background)
+
+
+class KullbackLeiblerConvexConjugate(Function):
+ """The convex conjugate of Kullback-Leibler divergence functional.
+
+ Notes
+ -----
+ The functional :math:`F^*` with prior :math:`g>0` is given by:
+
+ .. math::
+ F^*(x)
+ =
+ \\begin{cases}
+ \\sum_{i} \left( -g_i \ln(1 - x_i) \\right)
+ & \\text{if } x_i < 1 \\forall i
+ \\\\
+ +\\infty & \\text{else}
+ \\end{cases}
+
+ See Also
+ --------
+ KullbackLeibler : convex conjugate functional
+ """
+
+ def __init__(self, data, background):
+ self.data = data
+ self.background = background
+
+ def __call__(self, x):
+ y = self.data
+ r = self.background
+
+ tmp = numpy.sum(- x * r - y * numpy.log(1 - x))
+
+ if numpy.isnan(tmp):
+ # In this case, some element was larger than or equal to one
+ return numpy.inf
+ else:
+ return tmp
+
+
+ def prox(self, x, tau, out=None):
+ # Let y = data, r = background, z = x + tau * r
+ # Compute 0.5 * (z + 1 - sqrt((z - 1)**2 + 4 * tau * y))
+ # Currently it needs 3 extra copies of memory.
+
+ if out is None:
+ out = x.copy()
+
+ # define short variable names
+ try: # this should be standard SIRF/CIL mode
+ y = self.data.as_array()
+ r = self.background.as_array()
+ x = x.as_array()
+
+ try:
+ taua = tau.as_array()
+ except:
+ taua = tau
+
+ z = x + tau * r
+
+ out.fill(0.5 * (z + 1 - numpy.sqrt((z - 1) ** 2 + 4 * taua * y)))
+
+ return out
+
+ except: # e.g. for NumPy
+ y = self.data
+ r = self.background
+
+ try:
+ taua = tau.as_array()
+ except:
+ taua = tau
+
+ z = x + tau * r
+
+ out[:] = 0.5 * (z + 1 - numpy.sqrt((z - 1) ** 2 + 4 * taua * y))
+
+ return out
+
+ @property
+ def convex_conj(self):
+ return KullbackLeibler(self.data, self.background)
+
+
+def mult(x, y):
+ try:
+ xa = x.as_array()
+ except:
+ xa = x
+
+ out = y.clone()
+ out.fill(xa * y.as_array())
+
+ return out
diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py new file mode 100644 index 0000000..166014b --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py @@ -0,0 +1,484 @@ + # -*- coding: utf-8 -*- +""" +Created on Tue Mar 5 16:04:45 2019 + +@author: ofn77899 +""" +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + +import numpy +from numbers import Number +import functools +from ccpi.framework import DataContainer +#from ccpi.framework import AcquisitionData, ImageData +#from ccpi.optimisation.operators import Operator, LinearOperator + +class BlockDataContainer(object): + '''Class to hold DataContainers as column vector + + Provides basic algebra between BlockDataContainer's, DataContainer's and + subclasses and Numbers + + 1) algebra between `BlockDataContainer`s will be element-wise, only if + the shape of the 2 `BlockDataContainer`s is the same, otherwise it + will fail + 2) algebra between `BlockDataContainer`s and `list` or `numpy array` will + work as long as the number of `rows` and element of the arrays match, + indipendently on the fact that the `BlockDataContainer` could be nested + 3) algebra between `BlockDataContainer` and one `DataContainer` is possible. + It will require that all the `DataContainers` in the block to be + compatible with the `DataContainer` we want to algebra with. Should we + require that the `DataContainer` is the same type? Like `ImageData` or `AcquisitionData`? + 4) algebra between `BlockDataContainer` and a `Number` is possible and it + will be done with each element of the `BlockDataContainer` even if nested + + A = [ [B,C] , D] + A * 3 = [ 3 * [B,C] , 3* D] = [ [ 3*B, 3*C] , 3*D ] + + ''' + ADD = 'add' + SUBTRACT = 'subtract' + MULTIPLY = 'multiply' + DIVIDE = 'divide' + POWER = 'power' + __array_priority__ = 1 + __container_priority__ = 2 + def __init__(self, *args, **kwargs): + '''''' + self.containers = args + self.index = 0 + shape = kwargs.get('shape', None) + if shape is None: + shape = (len(args),1) +# shape = (len(args),1) + self.shape = shape + + n_elements = functools.reduce(lambda x,y: x*y, shape, 1) + if len(args) != n_elements: + raise ValueError( + 'Dimension and size do not match: expected {} got {}' + .format(n_elements, len(args))) + + + def __iter__(self): + '''BlockDataContainer is Iterable''' + return self + def next(self): + '''python2 backwards compatibility''' + return self.__next__() + def __next__(self): + try: + out = self[self.index] + except IndexError as ie: + raise StopIteration() + self.index+=1 + return out + + def is_compatible(self, other): + '''basic check if the size of the 2 objects fit''' + + if isinstance(other, Number): + return True + elif isinstance(other, (list, numpy.ndarray)) : + for ot in other: + if not isinstance(ot, (Number,\ + numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\ + numpy.float, numpy.float16, numpy.float32, numpy.float64, \ + numpy.complex)): + raise ValueError('List/ numpy array can only contain numbers {}'\ + .format(type(ot))) + return len(self.containers) == len(other) + elif issubclass(other.__class__, DataContainer): + ret = True + for i, el in enumerate(self.containers): + if isinstance(el, BlockDataContainer): + a = el.is_compatible(other) + else: + a = el.shape == other.shape + ret = ret and a + return ret + #return self.get_item(0).shape == other.shape + return len(self.containers) == len(other.containers) + + def get_item(self, row): + if row > self.shape[0]: + raise ValueError('Requested row {} > max {}'.format(row, self.shape[0])) + return self.containers[row] + + def __getitem__(self, row): + return self.get_item(row) + + def add(self, other, *args, **kwargs): + '''Algebra: add method of BlockDataContainer with number/DataContainer or BlockDataContainer + + :param: other (number, DataContainer or subclasses or BlockDataContainer + :param: out (optional): provides a placehold for the resul. + ''' + out = kwargs.get('out', None) + if out is not None: + self.binary_operations(BlockDataContainer.ADD, other, *args, **kwargs) + else: + return self.binary_operations(BlockDataContainer.ADD, other, *args, **kwargs) + def subtract(self, other, *args, **kwargs): + '''Algebra: subtract method of BlockDataContainer with number/DataContainer or BlockDataContainer + + :param: other (number, DataContainer or subclasses or BlockDataContainer + :param: out (optional): provides a placehold for the resul. + ''' + out = kwargs.get('out', None) + if out is not None: + self.binary_operations(BlockDataContainer.SUBTRACT, other, *args, **kwargs) + else: + return self.binary_operations(BlockDataContainer.SUBTRACT, other, *args, **kwargs) + def multiply(self, other, *args, **kwargs): + '''Algebra: multiply method of BlockDataContainer with number/DataContainer or BlockDataContainer + + :param: other (number, DataContainer or subclasses or BlockDataContainer + :param: out (optional): provides a placehold for the resul. + ''' + out = kwargs.get('out', None) + if out is not None: + self.binary_operations(BlockDataContainer.MULTIPLY, other, *args, **kwargs) + else: + return self.binary_operations(BlockDataContainer.MULTIPLY, other, *args, **kwargs) + def divide(self, other, *args, **kwargs): + '''Algebra: divide method of BlockDataContainer with number/DataContainer or BlockDataContainer + + :param: other (number, DataContainer or subclasses or BlockDataContainer + :param: out (optional): provides a placehold for the resul. + ''' + out = kwargs.get('out', None) + if out is not None: + self.binary_operations(BlockDataContainer.DIVIDE, other, *args, **kwargs) + else: + return self.binary_operations(BlockDataContainer.DIVIDE, other, *args, **kwargs) + + + def binary_operations(self, operation, other, *args, **kwargs): + '''Algebra: generic method of algebric operation with BlockDataContainer with number/DataContainer or BlockDataContainer + + Provides commutativity with DataContainer and subclasses, i.e. this + class's reverse algebric methods take precedence w.r.t. direct algebric + methods of DataContainer and subclasses. + + This method is not to be used directly + ''' + if not self.is_compatible(other): + raise ValueError('Incompatible for divide') + out = kwargs.get('out', None) + if isinstance(other, Number) or issubclass(other.__class__, DataContainer): + # try to do algebra with one DataContainer. Will raise error if not compatible + kw = kwargs.copy() + res = [] + for i,el in enumerate(self.containers): + if operation == BlockDataContainer.ADD: + op = el.add + elif operation == BlockDataContainer.SUBTRACT: + op = el.subtract + elif operation == BlockDataContainer.MULTIPLY: + op = el.multiply + elif operation == BlockDataContainer.DIVIDE: + op = el.divide + elif operation == BlockDataContainer.POWER: + op = el.power + else: + raise ValueError('Unsupported operation', operation) + if out is not None: + kw['out'] = out.get_item(i) + op(other, *args, **kw) + else: + res.append(op(other, *args, **kw)) + if out is not None: + return + else: + return type(self)(*res, shape=self.shape) + elif isinstance(other, (list, numpy.ndarray, BlockDataContainer)): + # try to do algebra with one DataContainer. Will raise error if not compatible + kw = kwargs.copy() + res = [] + if isinstance(other, BlockDataContainer): + the_other = other.containers + else: + the_other = other + for i,zel in enumerate(zip ( self.containers, the_other) ): + el = zel[0] + ot = zel[1] + if operation == BlockDataContainer.ADD: + op = el.add + elif operation == BlockDataContainer.SUBTRACT: + op = el.subtract + elif operation == BlockDataContainer.MULTIPLY: + op = el.multiply + elif operation == BlockDataContainer.DIVIDE: + op = el.divide + elif operation == BlockDataContainer.POWER: + op = el.power + else: + raise ValueError('Unsupported operation', operation) + if out is not None: + kw['out'] = out.get_item(i) + op(ot, *args, **kw) + else: + res.append(op(ot, *args, **kw)) + if out is not None: + return + else: + return type(self)(*res, shape=self.shape) + return type(self)(*[ operation(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape) + else: + raise ValueError('Incompatible type {}'.format(type(other))) + + + def power(self, other, *args, **kwargs): + if not self.is_compatible(other): + raise ValueError('Incompatible for power') + out = kwargs.get('out', None) + if isinstance(other, Number): + return type(self)(*[ el.power(other, *args, **kwargs) for el in self.containers], shape=self.shape) + elif isinstance(other, list) or isinstance(other, numpy.ndarray): + return type(self)(*[ el.power(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape) + return type(self)(*[ el.power(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape) + + def maximum(self,other, *args, **kwargs): + if not self.is_compatible(other): + raise ValueError('Incompatible for maximum') + out = kwargs.get('out', None) + if isinstance(other, Number): + return type(self)(*[ el.maximum(other, *args, **kwargs) for el in self.containers], shape=self.shape) + elif isinstance(other, list) or isinstance(other, numpy.ndarray): + return type(self)(*[ el.maximum(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape) + return type(self)(*[ el.maximum(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape) + + ## unary operations + def abs(self, *args, **kwargs): + return type(self)(*[ el.abs(*args, **kwargs) for el in self.containers], shape=self.shape) + def sign(self, *args, **kwargs): + return type(self)(*[ el.sign(*args, **kwargs) for el in self.containers], shape=self.shape) + def sqrt(self, *args, **kwargs): + return type(self)(*[ el.sqrt(*args, **kwargs) for el in self.containers], shape=self.shape) + def conjugate(self, out=None): + return type(self)(*[el.conjugate() for el in self.containers], shape=self.shape) + + ## reductions + + def sum(self, *args, **kwargs): + return numpy.sum([ el.sum(*args, **kwargs) for el in self.containers]) + + def squared_norm(self): + y = numpy.asarray([el.squared_norm() for el in self.containers]) + return y.sum() + + + def norm(self): + return numpy.sqrt(self.squared_norm()) + + def pnorm(self, p=2): + + if p==1: + return sum(self.abs()) + elif p==2: + return sum([el*el for el in self.containers]).sqrt() + else: + return ValueError('Not implemented') + + def copy(self): + '''alias of clone''' + return self.clone() + def clone(self): + return type(self)(*[el.copy() for el in self.containers], shape=self.shape) + def fill(self, other): + if isinstance (other, BlockDataContainer): + if not self.is_compatible(other): + raise ValueError('Incompatible containers') + for el,ot in zip(self.containers, other.containers): + el.fill(ot) + else: + return ValueError('Cannot fill with object provided {}'.format(type(other))) + + def __add__(self, other): + return self.add( other ) + # __radd__ + + def __sub__(self, other): + return self.subtract( other ) + # __rsub__ + + def __mul__(self, other): + return self.multiply(other) + # __rmul__ + + def __div__(self, other): + return self.divide(other) + # __rdiv__ + def __truediv__(self, other): + return self.divide(other) + + def __pow__(self, other): + return self.power(other) + # reverse operand + def __radd__(self, other): + '''Reverse addition + + to make sure that this method is called rather than the __mul__ of a numpy array + the class constant __array_priority__ must be set > 0 + https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__ + ''' + return self + other + # __radd__ + + def __rsub__(self, other): + '''Reverse subtraction + + to make sure that this method is called rather than the __mul__ of a numpy array + the class constant __array_priority__ must be set > 0 + https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__ + ''' + return (-1 * self) + other + # __rsub__ + + def __rmul__(self, other): + '''Reverse multiplication + + to make sure that this method is called rather than the __mul__ of a numpy array + the class constant __array_priority__ must be set > 0 + https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__ + ''' + return self * other + # __rmul__ + + def __rdiv__(self, other): + '''Reverse division + + to make sure that this method is called rather than the __mul__ of a numpy array + the class constant __array_priority__ must be set > 0 + https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__ + ''' + return pow(self / other, -1) + # __rdiv__ + def __rtruediv__(self, other): + '''Reverse truedivision + + to make sure that this method is called rather than the __mul__ of a numpy array + the class constant __array_priority__ must be set > 0 + https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__ + ''' + return self.__rdiv__(other) + + def __rpow__(self, other): + '''Reverse power + + to make sure that this method is called rather than the __mul__ of a numpy array + the class constant __array_priority__ must be set > 0 + https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__ + ''' + return other.power(self) + + def __iadd__(self, other): + '''Inline addition''' + if isinstance (other, BlockDataContainer): + for el,ot in zip(self.containers, other.containers): + el += ot + elif isinstance(other, Number): + for el in self.containers: + el += other + elif isinstance(other, list) or isinstance(other, numpy.ndarray): + if not self.is_compatible(other): + raise ValueError('Incompatible for __iadd__') + for el,ot in zip(self.containers, other): + el += ot + return self + # __iadd__ + + def __isub__(self, other): + '''Inline subtraction''' + if isinstance (other, BlockDataContainer): + for el,ot in zip(self.containers, other.containers): + el -= ot + elif isinstance(other, Number): + for el in self.containers: + el -= other + elif isinstance(other, list) or isinstance(other, numpy.ndarray): + if not self.is_compatible(other): + raise ValueError('Incompatible for __isub__') + for el,ot in zip(self.containers, other): + el -= ot + return self + # __isub__ + + def __imul__(self, other): + '''Inline multiplication''' + if isinstance (other, BlockDataContainer): + for el,ot in zip(self.containers, other.containers): + el *= ot + elif isinstance(other, Number): + for el in self.containers: + el *= other + elif isinstance(other, list) or isinstance(other, numpy.ndarray): + if not self.is_compatible(other): + raise ValueError('Incompatible for __imul__') + for el,ot in zip(self.containers, other): + el *= ot + return self + # __imul__ + + def __idiv__(self, other): + '''Inline division''' + if isinstance (other, BlockDataContainer): + for el,ot in zip(self.containers, other.containers): + el /= ot + elif isinstance(other, Number): + for el in self.containers: + el /= other + elif isinstance(other, list) or isinstance(other, numpy.ndarray): + if not self.is_compatible(other): + raise ValueError('Incompatible for __idiv__') + for el,ot in zip(self.containers, other): + el /= ot + return self + # __rdiv__ + def __itruediv__(self, other): + '''Inline truedivision''' + return self.__idiv__(other) + + def dot(self, other): +# + tmp = [ self.containers[i].dot(other.containers[i]) for i in range(self.shape[0])] + return sum(tmp) + + + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry, BlockGeometry + import numpy + + N, M = 2, 3 + ig = ImageGeometry(N, M) + BG = BlockGeometry(ig, ig) + + U = BG.allocate('random_int') + V = BG.allocate('random_int') + + + print ("test sum BDC " ) + w = U[0].as_array() + U[1].as_array() + w1 = sum(U).as_array() + numpy.testing.assert_array_equal(w, w1) + + print ("test sum BDC " ) + z = numpy.sqrt(U[0].as_array()**2 + U[1].as_array()**2) + z1 = sum(U**2).sqrt().as_array() + numpy.testing.assert_array_equal(z, z1) + + z2 = U.pnorm(2) + + zzz = U.dot(V) + + + + + + diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py new file mode 100644 index 0000000..ed44d99 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py @@ -0,0 +1,80 @@ +from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy
+from numbers import Number
+import functools
+from ccpi.framework import BlockDataContainer
+#from ccpi.optimisation.operators import Operator, LinearOperator
+
+class BlockGeometry(object):
+ '''Class to hold Geometry as column vector'''
+ #__array_priority__ = 1
+ def __init__(self, *args, **kwargs):
+ ''''''
+ self.geometries = args
+ self.index = 0
+
+ shape = (len(args),1)
+ self.shape = shape
+
+ n_elements = functools.reduce(lambda x,y: x*y, shape, 1)
+ if len(args) != n_elements:
+ raise ValueError(
+ 'Dimension and size do not match: expected {} got {}'
+ .format(n_elements, len(args)))
+
+
+ def get_item(self, index):
+ '''returns the Geometry in the BlockGeometry located at position index'''
+ return self.geometries[index]
+
+ def allocate(self, value=0, dimension_labels=None, **kwargs):
+
+ symmetry = kwargs.get('symmetry',False)
+ containers = [geom.allocate(value) for geom in self.geometries]
+
+ if symmetry == True:
+
+ # for 2x2
+ # [ ig11, ig12\
+ # ig21, ig22]
+
+ # Row-wise Order
+
+ if len(containers)==4:
+ containers[1]=containers[2]
+
+ # for 3x3
+ # [ ig11, ig12, ig13\
+ # ig21, ig22, ig23\
+ # ig31, ig32, ig33]
+
+ elif len(containers)==9:
+ containers[1]=containers[3]
+ containers[2]=containers[6]
+ containers[5]=containers[7]
+
+ # for 4x4
+ # [ ig11, ig12, ig13, ig14\
+ # ig21, ig22, ig23, ig24\
+ # ig31, ig32, ig33, ig34
+ # ig41, ig42, ig43, ig44]
+
+ elif len(containers) == 16:
+ containers[1]=containers[4]
+ containers[2]=containers[8]
+ containers[3]=containers[12]
+ containers[6]=containers[9]
+ containers[7]=containers[10]
+ containers[11]=containers[15]
+
+
+
+
+ return BlockDataContainer(*containers)
+
+
+
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/Algorithm.py new file mode 100644 index 0000000..12cbabc --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/Algorithm.py @@ -0,0 +1,185 @@ +# -*- 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 2019 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 time +from numbers import Integral + +class Algorithm(object): + '''Base class for iterative algorithms + + provides the minimal infrastructure. + Algorithms are iterables so can be easily run in a for loop. They will + stop as soon as the stop cryterion is met. + The user is required to implement the set_up, __init__, update and + and update_objective methods + + A courtesy method run is available to run n iterations. The method accepts + a callback function that receives the current iteration number and the actual objective + value and can be used to trigger print to screens and other user interactions. The run + method will stop when the stopping cryterion is met. + ''' + + def __init__(self): + '''Constructor + + Set the minimal number of parameters: + iteration: current iteration number + max_iteration: maximum number of iterations + memopt: whether to use memory optimisation () + timing: list to hold the times it took to run each iteration + update_objectice_interval: the interval every which we would save the current + objective. 1 means every iteration, 2 every 2 iteration + and so forth. This is by default 1 and should be increased + when evaluating the objective is computationally expensive. + ''' + self.iteration = 0 + self.__max_iteration = 0 + self.__loss = [] + self.memopt = False + self.timing = [] + self.update_objective_interval = 1 + def set_up(self, *args, **kwargs): + '''Set up the algorithm''' + raise NotImplementedError() + def update(self): + '''A single iteration of the algorithm''' + raise NotImplementedError() + + def should_stop(self): + '''default stopping cryterion: number of iterations + + The user can change this in concrete implementatition of iterative algorithms.''' + return self.max_iteration_stop_cryterion() + + def max_iteration_stop_cryterion(self): + '''default stop cryterion for iterative algorithm: max_iteration reached''' + return self.iteration >= self.max_iteration + def __iter__(self): + '''Algorithm is an iterable''' + return self + def next(self): + '''Algorithm is an iterable + + python2 backwards compatibility''' + return self.__next__() + def __next__(self): + '''Algorithm is an iterable + + calling this method triggers update and update_objective + ''' + if self.should_stop(): + raise StopIteration() + else: + time0 = time.time() + self.update() + self.timing.append( time.time() - time0 ) + if self.iteration % self.update_objective_interval == 0: + self.update_objective() + self.iteration += 1 + def get_output(self): + '''Returns the solution found''' + return self.x + def get_last_loss(self): + '''Returns the last stored value of the loss function + + if update_objective_interval is 1 it is the value of the objective at the current + iteration. If update_objective_interval > 1 it is the last stored value. + ''' + return self.__loss[-1] + def get_last_objective(self): + '''alias to get_last_loss''' + return self.get_last_loss() + def update_objective(self): + '''calculates the objective with the current solution''' + raise NotImplementedError() + @property + def loss(self): + '''returns the list of the values of the objective during the iteration + + The length of this list may be shorter than the number of iterations run when + the update_objective_interval > 1 + ''' + return self.__loss + @property + def objective(self): + '''alias of loss''' + return self.loss + @property + def max_iteration(self): + '''gets the maximum number of iterations''' + return self.__max_iteration + @max_iteration.setter + def max_iteration(self, value): + '''sets the maximum number of iterations''' + assert isinstance(value, int) + self.__max_iteration = value + @property + def update_objective_interval(self): + return self.__update_objective_interval + @update_objective_interval.setter + def update_objective_interval(self, value): + if isinstance(value, Integral): + if value >= 1: + self.__update_objective_interval = value + else: + raise ValueError('Update objective interval must be an integer >= 1') + else: + raise ValueError('Update objective interval must be an integer >= 1') + def run(self, iterations, verbose=True, callback=None): + '''run n iterations and update the user with the callback if specified''' + if self.should_stop(): + print ("Stop cryterion has been reached.") + i = 0 + +# print("Iteration {:<5} Primal {:<5} Dual {:<5} PDgap".format('','','')) + for _ in self: + + + if self.iteration % self.update_objective_interval == 0: + + if callback is not None: + callback(self.iteration, self.get_last_objective(), self.x) + + else: + + if verbose: + +# if verbose and self.iteration % self.update_objective_interval == 0: + #pass + # \t for tab +# print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\ +# format(self.iteration, self.max_iteration,'', \ +# self.get_last_objective()[0],'',\ +# self.get_last_objective()[1],'',\ +# self.get_last_objective()[2])) + + + print ("Iteration {}/{}, {}".format(self.iteration, + self.max_iteration, self.get_last_objective()) ) + + #print ("Iteration {}/{}, Primal, Dual, PDgap = {}".format(self.iteration, + # self.max_iteration, self.get_last_objective()) ) + + +# else: +# if callback is not None: +# callback(self.iteration, self.get_last_objective(), self.x) + i += 1 + if i == iterations: + break + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/FISTA.py new file mode 100644 index 0000000..ee51049 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/FISTA.py @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 21 11:07:30 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.functions import ZeroFunction +import numpy + +class FISTA(Algorithm): + '''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 + ''' + + def __init__(self, **kwargs): + '''initialisation can be done at creation time if all + proper variables are passed or later with set_up''' + super(FISTA, self).__init__() + self.f = None + self.g = None + self.invL = None + self.t_old = 1 + args = ['x_init', 'f', 'g', 'opt'] + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(kwargs['x_init'], + f=kwargs['f'], + g=kwargs['g'], + opt=kwargs['opt']) + + def set_up(self, x_init, f=None, g=None, opt=None): + + # default inputs + if f is None: + self.f = ZeroFunction() + else: + self.f = f + if g is None: + g = ZeroFunction() + self.g = g + else: + self.g = g + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-4, 'memopt':False} + + self.tol = opt['tol'] if 'tol' in opt.keys() else 1e-4 + memopt = opt['memopt'] if 'memopt' in opt.keys() else False + self.memopt = memopt + + # initialization + if memopt: + self.y = x_init.copy() + self.x_old = x_init.copy() + self.x = x_init.copy() + self.u = x_init.copy() + else: + self.x_old = x_init.copy() + self.y = x_init.copy() + + #timing = numpy.zeros(max_iter) + #criter = numpy.zeros(max_iter) + + + self.invL = 1/f.L + + self.t_old = 1 + + def update(self): + # algorithm loop + #for it in range(0, max_iter): + + if self.memopt: + # u = y - invL*f.grad(y) + # store the result in x_old + self.f.gradient(self.y, out=self.u) + self.u.__imul__( -self.invL ) + self.u.__iadd__( self.y ) + # x = g.prox(u,invL) + self.g.proximal(self.u, self.invL, out=self.x) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + # y = x + (t_old-1)/t*(x-x_old) + self.x.subtract(self.x_old, out=self.y) + self.y.__imul__ ((self.t_old-1)/self.t) + self.y.__iadd__( self.x ) + + self.x_old.fill(self.x) + self.t_old = self.t + + + else: + u = self.y - self.invL*self.f.gradient(self.y) + + self.x = self.g.proximal(u,self.invL) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + self.y = self.x + (self.t_old-1)/self.t*(self.x-self.x_old) + + self.x_old = self.x.copy() + self.t_old = self.t + + def update_objective(self): + self.loss.append( self.f(self.x) + self.g(self.x) )
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py new file mode 100644 index 0000000..0f5e8ef --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py @@ -0,0 +1,236 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Feb 4 16:18:06 2019 + +@author: evangelos +""" +from ccpi.optimisation.algorithms import Algorithm +from ccpi.framework import ImageData, DataContainer +import numpy as np +import numpy +import time +from ccpi.optimisation.operators import BlockOperator +from ccpi.framework import BlockDataContainer +from ccpi.optimisation.functions import FunctionOperatorComposition +import matplotlib.pyplot as plt + +class PDHG(Algorithm): + '''Primal Dual Hybrid Gradient''' + + def __init__(self, **kwargs): + super(PDHG, self).__init__() + self.f = kwargs.get('f', None) + 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.memopt = kwargs.get('memopt', False) + + if self.f is not None and self.operator is not None and \ + self.g is not None: + print ("Calling from creator") + self.set_up(self.f, + self.operator, + self.g, + self.tau, + self.sigma) + + def set_up(self, f, g, operator, tau = None, sigma = None, opt = None, **kwargs): + # algorithmic parameters + + if sigma is None and tau is None: + raise ValueError('Need sigma*tau||K||^2<1') + + + self.x_old = self.operator.domain_geometry().allocate() + self.y_old = self.operator.range_geometry().allocate() + + self.xbar = self.x_old.copy() + self.x_tmp = self.x_old.copy() + self.x = self.x_old.copy() + + self.y_tmp = self.y_old.copy() + self.y = self.y_tmp.copy() + + + + #self.y = self.y_old.copy() + + + #if self.memopt: + # self.y_tmp = self.y_old.copy() + # self.x_tmp = self.x_old.copy() + + + # relaxation parameter + self.theta = 1 + + def update(self): + + if self.memopt: + # Gradient descent, Dual problem solution + # self.y_old += self.sigma * self.operator.direct(self.xbar) + self.operator.direct(self.xbar, out=self.y_tmp) + self.y_tmp *= self.sigma + self.y_tmp += self.y_old + + #self.y = self.f.proximal_conjugate(self.y_old, self.sigma) + self.f.proximal_conjugate(self.y_tmp, self.sigma, out=self.y) + + # Gradient ascent, Primal problem solution + self.operator.adjoint(self.y, out=self.x_tmp) + self.x_tmp *= -1*self.tau + self.x_tmp += self.x_old + + self.g.proximal(self.x_tmp, self.tau, out=self.x) + + #Update + self.x.subtract(self.x_old, out=self.xbar) + self.xbar *= self.theta + self.xbar += self.x + + self.x_old.fill(self.x) + self.y_old.fill(self.y) + + + else: + # Gradient descent, Dual problem solution + self.y_old += self.sigma * self.operator.direct(self.xbar) + self.y = self.f.proximal_conjugate(self.y_old, self.sigma) + + # Gradient ascent, Primal problem solution + self.x_old -= self.tau * self.operator.adjoint(self.y) + self.x = self.g.proximal(self.x_old, self.tau) + + #Update + self.x.subtract(self.x_old, out=self.xbar) + self.xbar *= self.theta + self.xbar += self.x + + self.x_old.fill(self.x) + self.y_old.fill(self.y) + + #xbar = x + theta * (x - x_old) +# self.xbar.fill(self.x) +# self.xbar -= self.x_old +# self.xbar *= self.theta +# self.xbar += self.x + +# self.x_old.fill(self.x) +# self.y_old.fill(self.y) + + + + def update_objective(self): + + p1 = self.f(self.operator.direct(self.x)) + self.g(self.x) + d1 = -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(-1*self.operator.adjoint(self.y))) + + self.loss.append([p1,d1,p1-d1]) + + + +def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-6, 'niter': 500, 'show_iter': 100, \ + 'memopt': False} + + if sigma is None and tau is None: + raise ValueError('Need sigma*tau||K||^2<1') + + niter = opt['niter'] if 'niter' 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 + show_iter = opt['show_iter'] if 'show_iter' in opt.keys() else False + stop_crit = opt['stop_crit'] if 'stop_crit' in opt.keys() else False + + x_old = operator.domain_geometry().allocate() + y_old = operator.range_geometry().allocate() + + xbar = x_old.copy() + x_tmp = x_old.copy() + x = x_old.copy() + + y_tmp = y_old.copy() + y = y_tmp.copy() + + + # relaxation parameter + theta = 1 + + t = time.time() + + primal = [] + dual = [] + pdgap = [] + + + for i in range(niter): + + + if not memopt: + + y_tmp = y_old + sigma * operator.direct(xbar) + y = f.proximal_conjugate(y_tmp, sigma) + + x_tmp = x_old - tau*operator.adjoint(y) + x = g.proximal(x_tmp, tau) + + x.subtract(x_old, out=xbar) + xbar *= theta + xbar += x + + if i%50==0: + + p1 = f(operator.direct(x)) + g(x) + d1 = - ( f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) ) + primal.append(p1) + dual.append(d1) + pdgap.append(p1-d1) + print(p1, d1, p1-d1) + + x_old.fill(x) + y_old.fill(y) + + + else: + + operator.direct(xbar, out = y_tmp) + y_tmp *= sigma + y_tmp += y_old + f.proximal_conjugate(y_tmp, sigma, out=y) + + operator.adjoint(y, out = x_tmp) + x_tmp *= -1*tau + x_tmp += x_old + + g.proximal(x_tmp, tau, out = x) + + x.subtract(x_old, out=xbar) + xbar *= theta + xbar += x + + if i%50==0: + + p1 = f(operator.direct(x)) + g(x) + d1 = - ( f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) ) + primal.append(p1) + dual.append(d1) + pdgap.append(p1-d1) + print(p1, d1, p1-d1) + + x_old.fill(x) + y_old.fill(y) + + + + + t_end = time.time() + + return x, t_end - t, primal, dual, pdgap + + + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/funcs.py b/Wrappers/Python/build/lib/ccpi/optimisation/funcs.py new file mode 100644 index 0000000..b2b9791 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/funcs.py @@ -0,0 +1,265 @@ +# -*- 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) ) diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py new file mode 100644 index 0000000..0836324 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py @@ -0,0 +1,221 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 10:01:31 2019 + +@author: evangelos +""" + +from ccpi.optimisation.functions import Function +from ccpi.framework import BlockDataContainer +from numbers import Number + +class BlockFunction(Function): + + '''BlockFunction acts as a separable sum function, i.e., + + f = [f_1,...,f_n] + + f([x_1,...,x_n]) = f_1(x_1) + .... + f_n(x_n) + + ''' + def __init__(self, *functions): + + self.functions = functions + self.length = len(self.functions) + + super(BlockFunction, self).__init__() + + def __call__(self, x): + + '''Evaluates the BlockFunction at a BlockDataContainer x + + :param: x (BlockDataContainer): must have as many rows as self.length + + returns sum(f_i(x_i)) + ''' + + if self.length != x.shape[0]: + raise ValueError('BlockFunction and BlockDataContainer have incompatible size') + t = 0 + for i in range(x.shape[0]): + t += self.functions[i](x.get_item(i)) + return t + + def convex_conjugate(self, x): + + ''' Evaluate convex conjugate of BlockFunction at x + + returns sum(f_i^{*}(x_i)) + + ''' + t = 0 + for i in range(x.shape[0]): + t += self.functions[i].convex_conjugate(x.get_item(i)) + return t + + + def proximal_conjugate(self, x, tau, out = None): + + ''' Evaluate Proximal Operator of tau * f(\cdot) at x + + prox_{tau*f}(x) = sum_{i} prox_{tau*f_{i}}(x_{i}) + + + ''' + + if out is not None: + if isinstance(tau, Number): + for i in range(self.length): + self.functions[i].proximal_conjugate(x.get_item(i), tau, out=out.get_item(i)) + else: + for i in range(self.length): + self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(i),out=out.get_item(i)) + + else: + + out = [None]*self.length + if isinstance(tau, Number): + for i in range(self.length): + out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau) + else: + for i in range(self.length): + out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(i)) + + return BlockDataContainer(*out) + + + def proximal(self, x, tau, out = None): + + ''' Evaluate Proximal Operator of tau * f^{*}(\cdot) at x + + prox_{tau*f^{*}}(x) = sum_{i} prox_{tau*f^{*}_{i}}(x_{i}) + + + ''' + + if out is None: + + out = [None]*self.length + if isinstance(tau, Number): + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau) + else: + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i)) + + return BlockDataContainer(*out) + + else: + if isinstance(tau, Number): + for i in range(self.length): + self.functions[i].proximal(x.get_item(i), tau, out[i]) + else: + for i in range(self.length): + self.functions[i].proximal(x.get_item(i), tau.get_item(i), out[i]) + + + + def gradient(self,x, out=None): + + ''' Evaluate gradient of f at x: f'(x) + + returns: BlockDataContainer [f_{1}'(x_{1}), ... , f_{n}'(x_{n})] + + ''' + + out = [None]*self.length + for i in range(self.length): + out[i] = self.functions[i].gradient(x.get_item(i)) + + return BlockDataContainer(*out) + + + +if __name__ == '__main__': + + M, N, K = 2,3,5 + + from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm + from ccpi.framework import ImageGeometry, BlockGeometry + from ccpi.optimisation.operators import Gradient, Identity, BlockOperator + import numpy + import numpy as np + + + ig = ImageGeometry(M, N) + BG = BlockGeometry(ig, ig) + + u = ig.allocate('random_int') + B = BlockOperator( Gradient(ig), Identity(ig) ) + + U = B.direct(u) + b = ig.allocate('random_int') + + f1 = 10 * MixedL21Norm() + f2 = 0.5 * L2NormSquared(b=b) + + f = BlockFunction(f1, f2) + tau = 0.3 + + print( " without out " ) + res_no_out = f.proximal_conjugate( U, tau) + res_out = B.range_geometry().allocate() + f.proximal_conjugate( U, tau, out = res_out) + + numpy.testing.assert_array_almost_equal(res_no_out[0][0].as_array(), \ + res_out[0][0].as_array(), decimal=4) + + numpy.testing.assert_array_almost_equal(res_no_out[0][1].as_array(), \ + res_out[0][1].as_array(), decimal=4) + + numpy.testing.assert_array_almost_equal(res_no_out[1].as_array(), \ + res_out[1].as_array(), decimal=4) + + + + ########################################################################## + + + + + + + +# zzz = B.range_geometry().allocate('random_int') +# www = B.range_geometry().allocate() +# www.fill(zzz) + +# res[0].fill(z) + + + + +# f.proximal_conjugate(z, sigma, out = res) + +# print(z1[0][0].as_array()) +# print(res[0][0].as_array()) + + + + +# U = BG.allocate('random_int') +# RES = BG.allocate() +# f = BlockFunction(f1, f2) +# +# z = f.proximal_conjugate(U, 0.2) +# f.proximal_conjugate(U, 0.2, out = RES) +# +# print(z[0].as_array()) +# print(RES[0].as_array()) +# +# print(z[1].as_array()) +# print(RES[1].as_array()) + + + + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition.py new file mode 100644 index 0000000..8895f3d --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 09:55:36 2019 + +@author: evangelos +""" + +from ccpi.optimisation.functions import Function +from ccpi.optimisation.functions import ScaledFunction + + +class FunctionOperatorComposition(Function): + + ''' Function composition with Operator, i.e., f(Ax) + + A: operator + f: function + + ''' + + def __init__(self, function, operator): + + super(FunctionOperatorComposition, self).__init__() + + self.function = function + self.operator = operator + self.L = function.L * operator.norm()**2 + + + def __call__(self, x): + + ''' Evaluate FunctionOperatorComposition at x + + returns f(Ax) + + ''' + + return self.function(self.operator.direct(x)) + + def gradient(self, x, out=None): +# + ''' Gradient takes into account the Operator''' + if out is None: + return self.operator.adjoint(self.function.gradient(self.operator.direct(x))) + else: + tmp = self.operator.range_geometry().allocate() + self.operator.direct(x, out=tmp) + self.function.gradient(tmp, out=tmp) + self.operator.adjoint(tmp, out=out) + + + + + #TODO do not know if we need it + #def call_adjoint(self, x): + # + # return self.function(self.operator.adjoint(x)) + + + #def convex_conjugate(self, x): + # + # ''' convex_conjugate does not take into account the Operator''' + # return self.function.convex_conjugate(x) + + + + + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + from ccpi.optimisation.operators import Gradient + from ccpi.optimisation.functions import L2NormSquared + + M, N, K = 2,3 + ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N) + + G = Gradient(ig) + alpha = 0.5 + + f = L2NormSquared() + f_comp = FunctionOperatorComposition(G, alpha * f) + x = ig.allocate('random_int') + print(f_comp.gradient(x).shape + + ) + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition_old.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition_old.py new file mode 100644 index 0000000..70511bb --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition_old.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 09:55:36 2019 + +@author: evangelos +""" + +from ccpi.optimisation.functions import Function +from ccpi.optimisation.functions import ScaledFunction + + +class FunctionOperatorComposition(Function): + + ''' Function composition with Operator, i.e., f(Ax) + + A: operator + f: function + + ''' + + def __init__(self, operator, function): + + super(FunctionOperatorComposition, self).__init__() + self.function = function + self.operator = operator + alpha = 1 + + if isinstance (function, ScaledFunction): + alpha = function.scalar + self.L = 2 * alpha * operator.norm()**2 + + + def __call__(self, x): + + ''' Evaluate FunctionOperatorComposition at x + + returns f(Ax) + + ''' + + return self.function(self.operator.direct(x)) + + #TODO do not know if we need it + def call_adjoint(self, x): + + return self.function(self.operator.adjoint(x)) + + + def convex_conjugate(self, x): + + ''' convex_conjugate does not take into account the Operator''' + return self.function.convex_conjugate(x) + + def proximal(self, x, tau, out=None): + + '''proximal does not take into account the Operator''' + if out is None: + return self.function.proximal(x, tau) + else: + self.function.proximal(x, tau, out=out) + + + def proximal_conjugate(self, x, tau, out=None): + + ''' proximal conjugate does not take into account the Operator''' + if out is None: + return self.function.proximal_conjugate(x, tau) + else: + self.function.proximal_conjugate(x, tau, out=out) + + def gradient(self, x, out=None): + + ''' Gradient takes into account the Operator''' + if out is None: + return self.operator.adjoint( + self.function.gradient(self.operator.direct(x)) + ) + else: + self.operator.adjoint( + self.function.gradient(self.operator.direct(x), + out=out) + ) + +
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/KullbackLeibler.py new file mode 100644 index 0000000..b53f669 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/KullbackLeibler.py @@ -0,0 +1,175 @@ +# -*- 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-2019 Evangelos Papoutsellis 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.optimisation.functions import Function +from ccpi.optimisation.functions.ScaledFunction import ScaledFunction +from ccpi.framework import ImageData, ImageGeometry +import functools + +class KullbackLeibler(Function): + + ''' Assume that data > 0 + + ''' + + def __init__(self,data, **kwargs): + + super(KullbackLeibler, self).__init__() + + self.b = data + self.bnoise = kwargs.get('bnoise', 0) + + + def __call__(self, x): + + res_tmp = numpy.zeros(x.shape) + + tmp = x + self.bnoise + ind = x.as_array()>0 + + res_tmp[ind] = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind]) + + return res_tmp.sum() + + def log(self, datacontainer): + '''calculates the in-place log of the datacontainer''' + if not functools.reduce(lambda x,y: x and y>0, + datacontainer.as_array().ravel(), True): + raise ValueError('KullbackLeibler. Cannot calculate log of negative number') + datacontainer.fill( numpy.log(datacontainer.as_array()) ) + + + def gradient(self, x, out=None): + + #TODO Division check + if out is None: + return 1 - self.b/(x + self.bnoise) + else: + x.add(self.bnoise, out=out) + self.b.divide(out, out=out) + out.subtract(1, out=out) + out *= -1 + + def convex_conjugate(self, x): + + tmp = self.b/(1-x) + ind = tmp.as_array()>0 + + return (self.b.as_array()[ind] * (numpy.log(tmp.as_array()[ind])-1)).sum() + + + def proximal(self, x, tau, out=None): + + if out is None: + return 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b ) .sqrt() ) + else: + + tmp = 0.5 *( (x - self.bnoise - tau) + + ( (x + self.bnoise - tau)**2 + 4*tau*self.b ) .sqrt() + ) + x.add(self.bnoise, out=out) + out -= tau + out *= out + tmp = self.b * (4 * tau) + out.add(tmp, out=out) + out.sqrt(out=out) + + x.subtract(self.bnoise, out=tmp) + tmp -= tau + + out += tmp + + out *= 0.5 + + def proximal_conjugate(self, x, tau, out=None): + + + if out is None: + z = x + tau * self.bnoise + return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) + else: +# z = x + tau * self.bnoise +# out.fill( 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) ) + + tmp1 = x + tau * self.bnoise - 1 + tmp2 = tmp1 + 2 + + self.b.multiply(4*tau, out=out) + tmp1.multiply(tmp1, out=tmp1) + out += tmp1 + out.sqrt(out=out) + + out *= -1 + out += tmp2 + out *= 0.5 + + + + def __rmul__(self, scalar): + + ''' Multiplication of L2NormSquared with a scalar + + Returns: ScaledFunction + + ''' + + return ScaledFunction(self, scalar) + + + + +if __name__ == '__main__': + + + from ccpi.framework import ImageGeometry + import numpy + + N, M = 2,3 + ig = ImageGeometry(N, M) + data = ImageData(numpy.random.randint(-10, 10, size=(M, N))) + x = ImageData(numpy.random.randint(-10, 10, size=(M, N))) + + bnoise = ImageData(numpy.random.randint(-10, 10, size=(M, N))) + + f = KullbackLeibler(data) + + print(f(x)) + +# numpy.random.seed(10) +# +# +# x = numpy.random.randint(-10, 10, size = (2,3)) +# b = numpy.random.randint(1, 10, size = (2,3)) +# +# ind1 = x>0 +# +# res = x[ind1] - b * numpy.log(x[ind1]) +# +## ind = x>0 +# +## y = x[ind] +# +# +# +# +# + + + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py new file mode 100644 index 0000000..e73da93 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py @@ -0,0 +1,303 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.optimisation.functions import Function +from ccpi.optimisation.functions.ScaledFunction import ScaledFunction +from ccpi.optimisation.functions import FunctionOperatorComposition + +class L2NormSquared(Function): + + ''' + + Cases: a) f(x) = \|x\|^{2}_{2} + + b) f(x) = ||x - b||^{2}_{2} + + ''' + + def __init__(self, **kwargs): + + super(L2NormSquared, self).__init__() + self.b = kwargs.get('b',None) + self.L = 2 + + def __call__(self, x): + + ''' Evaluate L2NormSquared at x: f(x) ''' + + y = x + if self.b is not None: + y = x - self.b + try: + return y.squared_norm() + except AttributeError as ae: + # added for compatibility with SIRF + return (y.norm()**2) + + def gradient(self, x, out=None): + + ''' Evaluate gradient of L2NormSquared at x: f'(x) ''' + + if out is not None: + + out.fill(x) + if self.b is not None: + out -= self.b + out *= 2 + + else: + + y = x + if self.b is not None: + y = x - self.b + return 2*y + + + def convex_conjugate(self, x): + + ''' Evaluate convex conjugate of L2NormSquared at x: f^{*}(x)''' + + tmp = 0 + + if self.b is not None: + tmp = (x * self.b).sum() + + return (1./4.) * x.squared_norm() + tmp + + + def proximal(self, x, tau, out = None): + + ''' Evaluate Proximal Operator of tau * f(\cdot) at x: + + prox_{tau*f(\cdot)}(x) = \argmin_{z} \frac{1}{2}|| z - x ||^{2}_{2} + tau * f(z) + + ''' + + if out is None: + + if self.b is None: + return x/(1+2*tau) + else: + + tmp = x.subtract(self.b) + tmp /= (1+2*tau) + tmp += self.b + return tmp + + else: + out.fill(x) + if self.b is not None: + out -= self.b + out /= (1+2*tau) + if self.b is not None: + out += self.b + + + + + def proximal_conjugate(self, x, tau, out=None): + + ''' Evaluate Proximal Operator of tau * f^{*}(\cdot) at x (i.e., the convex conjugate of f) : + + prox_{tau*f(\cdot)}(x) = \argmin_{z} \frac{1}{2}|| z - x ||^{2}_{2} + tau * f^{*}(z) + + ''' + + if out is None: + if self.b is not None: + return (x - tau*self.b)/(1 + tau/2) + else: + return x/(1 + tau/2) + else: + if self.b is not None: + x.subtract(tau*self.b, out=out) + out.divide(1+tau/2, out=out) + else: + x.divide(1 + tau/2, out=out) + + def __rmul__(self, scalar): + + ''' Multiplication of L2NormSquared with a scalar + + Returns: ScaledFunction + + ''' + + return ScaledFunction(self, scalar) + + + def composition(self, operator): + + return FunctionOperatorComposition(operator) + + + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + import numpy + # TESTS for L2 and scalar * L2 + + M, N, K = 2,3,5 + ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K) + u = ig.allocate('random_int') + b = ig.allocate('random_int') + + # check grad/call no data + f = L2NormSquared() + a1 = f.gradient(u) + a2 = 2 * u + numpy.testing.assert_array_almost_equal(a1.as_array(), a2.as_array(), decimal=4) + numpy.testing.assert_equal(f(u), u.squared_norm()) + + # check grad/call with data + f1 = L2NormSquared(b=b) + b1 = f1.gradient(u) + b2 = 2 * (u-b) + + numpy.testing.assert_array_almost_equal(b1.as_array(), b2.as_array(), decimal=4) + numpy.testing.assert_equal(f1(u), (u-b).squared_norm()) + + #check convex conjuagate no data + c1 = f.convex_conjugate(u) + c2 = 1/4 * u.squared_norm() + numpy.testing.assert_equal(c1, c2) + + #check convex conjuagate with data + d1 = f1.convex_conjugate(u) + d2 = (1/4) * u.squared_norm() + (u*b).sum() + numpy.testing.assert_equal(d1, d2) + + # check proximal no data + tau = 5 + e1 = f.proximal(u, tau) + e2 = u/(1+2*tau) + numpy.testing.assert_array_almost_equal(e1.as_array(), e2.as_array(), decimal=4) + + # check proximal with data + tau = 5 + h1 = f1.proximal(u, tau) + h2 = (u-b)/(1+2*tau) + b + numpy.testing.assert_array_almost_equal(h1.as_array(), h2.as_array(), decimal=4) + + # check proximal conjugate no data + tau = 0.2 + k1 = f.proximal_conjugate(u, tau) + k2 = u/(1 + tau/2 ) + numpy.testing.assert_array_almost_equal(k1.as_array(), k2.as_array(), decimal=4) + + # check proximal conjugate with data + l1 = f1.proximal_conjugate(u, tau) + l2 = (u - tau * b)/(1 + tau/2 ) + numpy.testing.assert_array_almost_equal(l1.as_array(), l2.as_array(), decimal=4) + + + # check scaled function properties + + # scalar + scalar = 100 + f_scaled_no_data = scalar * L2NormSquared() + f_scaled_data = scalar * L2NormSquared(b=b) + + # call + numpy.testing.assert_equal(f_scaled_no_data(u), scalar*f(u)) + numpy.testing.assert_equal(f_scaled_data(u), scalar*f1(u)) + + # grad + numpy.testing.assert_array_almost_equal(f_scaled_no_data.gradient(u).as_array(), scalar*f.gradient(u).as_array(), decimal=4) + numpy.testing.assert_array_almost_equal(f_scaled_data.gradient(u).as_array(), scalar*f1.gradient(u).as_array(), decimal=4) + + # conj + numpy.testing.assert_almost_equal(f_scaled_no_data.convex_conjugate(u), \ + f.convex_conjugate(u/scalar) * scalar, decimal=4) + + numpy.testing.assert_almost_equal(f_scaled_data.convex_conjugate(u), \ + scalar * f1.convex_conjugate(u/scalar), decimal=4) + + # proximal + numpy.testing.assert_array_almost_equal(f_scaled_no_data.proximal(u, tau).as_array(), \ + f.proximal(u, tau*scalar).as_array()) + + + numpy.testing.assert_array_almost_equal(f_scaled_data.proximal(u, tau).as_array(), \ + f1.proximal(u, tau*scalar).as_array()) + + + # proximal conjugate + numpy.testing.assert_array_almost_equal(f_scaled_no_data.proximal_conjugate(u, tau).as_array(), \ + (u/(1 + tau/(2*scalar) )).as_array(), decimal=4) + + numpy.testing.assert_array_almost_equal(f_scaled_data.proximal_conjugate(u, tau).as_array(), \ + ((u - tau * b)/(1 + tau/(2*scalar) )).as_array(), decimal=4) + + + + print( " ####### check without out ######### " ) + + + u_out_no_out = ig.allocate('random_int') + res_no_out = f_scaled_data.proximal_conjugate(u_out_no_out, 0.5) + print(res_no_out.as_array()) + + print( " ####### check with out ######### " ) + + res_out = ig.allocate() + f_scaled_data.proximal_conjugate(u_out_no_out, 0.5, out = res_out) + + print(res_out.as_array()) + + numpy.testing.assert_array_almost_equal(res_no_out.as_array(), \ + res_out.as_array(), decimal=4) + + + + ig1 = ImageGeometry(2,3) + + tau = 0.1 + + u = ig1.allocate('random_int') + b = ig1.allocate('random_int') + + scalar = 0.5 + f_scaled = scalar * L2NormSquared(b=b) + f_noscaled = L2NormSquared(b=b) + + + res1 = f_scaled.proximal(u, tau) + res2 = f_noscaled.proximal(u, tau*scalar) + +# res2 = (u + tau*b)/(1+tau) + + numpy.testing.assert_array_almost_equal(res1.as_array(), \ + res2.as_array(), decimal=4) + + + + + + + + + + + + + + + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/MixedL21Norm.py new file mode 100644 index 0000000..e8f6da4 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/MixedL21Norm.py @@ -0,0 +1,159 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.optimisation.functions import Function, ScaledFunction +from ccpi.framework import BlockDataContainer + +import functools + +class MixedL21Norm(Function): + + + ''' + f(x) = ||x||_{2,1} = \sum |x|_{2} + ''' + + def __init__(self, **kwargs): + + super(MixedL21Norm, self).__init__() + self.SymTensor = kwargs.get('SymTensor',False) + + def __call__(self, x): + + ''' Evaluates L2,1Norm at point x + + :param: x is a BlockDataContainer + + ''' + if not isinstance(x, BlockDataContainer): + raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) + + tmp = [ el**2 for el in x.containers ] + res = sum(tmp).sqrt().sum() + + return res + + def gradient(self, x, out=None): + return ValueError('Not Differentiable') + + def convex_conjugate(self,x): + + ''' This is the Indicator function of ||\cdot||_{2, \infty} + which is either 0 if ||x||_{2, \infty} or \infty + ''' + + return 0.0 + + #tmp = [ el**2 for el in x.containers ] + #print(sum(tmp).sqrt().as_array().max()) + #return sum(tmp).sqrt().as_array().max() + + def proximal(self, x, tau, out=None): + + ''' + For this we need to define a MixedL2,2 norm acting on BDC, + different form L2NormSquared which acts on DC + + ''' + pass + + def proximal_conjugate(self, x, tau, out=None): + + + if out is None: + tmp = [ el*el for el in x.containers] + res = sum(tmp).sqrt().maximum(1.0) + frac = [el/res for el in x.containers] + return BlockDataContainer(*frac) + + + #TODO this is slow, why??? +# return x.divide(x.pnorm().maximum(1.0)) + else: + + res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ) + res = res1.sqrt().maximum(1.0) + x.divide(res, out=out) + +# x.divide(sum([el*el for el in x.containers]).sqrt().maximum(1.0), out=out) + #TODO this is slow, why ??? +# x.divide(x.pnorm().maximum(1.0), out=out) + + + def __rmul__(self, scalar): + + ''' Multiplication of L2NormSquared with a scalar + + Returns: ScaledFunction + + ''' + return ScaledFunction(self, scalar) + + +# +if __name__ == '__main__': + + M, N, K = 2,3,5 + from ccpi.framework import BlockGeometry + import numpy + + ig = ImageGeometry(M, N) + + BG = BlockGeometry(ig, ig) + + U = BG.allocate('random_int') + + # Define no scale and scaled + f_no_scaled = MixedL21Norm() + f_scaled = 0.5 * MixedL21Norm() + + # call + + a1 = f_no_scaled(U) + a2 = f_scaled(U) + print(a1, 2*a2) + + + print( " ####### check without out ######### " ) + + + u_out_no_out = BG.allocate('random_int') + res_no_out = f_scaled.proximal_conjugate(u_out_no_out, 0.5) + print(res_no_out[0].as_array()) + + print( " ####### check with out ######### " ) +# + res_out = BG.allocate() + f_scaled.proximal_conjugate(u_out_no_out, 0.5, out = res_out) +# + print(res_out[0].as_array()) +# + numpy.testing.assert_array_almost_equal(res_no_out[0].as_array(), \ + res_out[0].as_array(), decimal=4) + + numpy.testing.assert_array_almost_equal(res_no_out[1].as_array(), \ + res_out[1].as_array(), decimal=4) +# + + + + + + + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ScaledFunction.py new file mode 100644 index 0000000..1db223b --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ScaledFunction.py @@ -0,0 +1,150 @@ +# -*- 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-2019 Evangelos Papoutsellis 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.
+from numbers import Number
+import numpy
+
+class ScaledFunction(object):
+
+ '''ScaledFunction
+
+ A class to represent the scalar multiplication of an Function with a scalar.
+ It holds a function and a scalar. Basically it returns the multiplication
+ of the product of the function __call__, convex_conjugate and gradient with the scalar.
+ For the rest it behaves like the function it holds.
+
+ Args:
+ function (Function): a Function or BlockOperator
+ scalar (Number): a scalar multiplier
+ Example:
+ The scaled operator behaves like the following:
+
+ '''
+ def __init__(self, function, scalar):
+ super(ScaledFunction, self).__init__()
+
+ if not isinstance (scalar, Number):
+ raise TypeError('expected scalar: got {}'.format(type(scalar)))
+ self.scalar = scalar
+ self.function = function
+
+ if self.function.L is not None:
+ self.L = self.scalar * self.function.L
+
+ def __call__(self,x, out=None):
+ '''Evaluates the function at x '''
+ return self.scalar * self.function(x)
+
+ def convex_conjugate(self, x):
+ '''returns the convex_conjugate of the scaled function '''
+ return self.scalar * self.function.convex_conjugate(x/self.scalar)
+
+ def gradient(self, x, out=None):
+ '''Returns the gradient of the function at x, if the function is differentiable'''
+ if out is None:
+ return self.scalar * self.function.gradient(x)
+ else:
+ self.function.gradient(x, out=out)
+ out *= self.scalar
+
+ def proximal(self, x, tau, out=None):
+ '''This returns the proximal operator for the function at x, tau
+ '''
+ if out is None:
+ return self.function.proximal(x, tau*self.scalar)
+ else:
+ self.function.proximal(x, tau*self.scalar, out = out)
+
+ def proximal_conjugate(self, x, tau, out = None):
+ '''This returns the proximal operator for the function at x, tau
+ '''
+ if out is None:
+ return self.scalar * self.function.proximal_conjugate(x/self.scalar, tau/self.scalar)
+ else:
+ self.function.proximal_conjugate(x/self.scalar, tau/self.scalar, out=out)
+ out *= self.scalar
+
+ def grad(self, x):
+ '''Alias of gradient(x,None)'''
+ warnings.warn('''This method will disappear in following
+ versions of the CIL. Use gradient instead''', DeprecationWarning)
+ return self.gradient(x, out=None)
+
+ def prox(self, x, tau):
+ '''Alias of proximal(x, tau, None)'''
+ warnings.warn('''This method will disappear in following
+ versions of the CIL. Use proximal instead''', DeprecationWarning)
+ return self.proximal(x, out=None)
+
+
+
+if __name__ == '__main__':
+
+ from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm
+ from ccpi.framework import ImageGeometry, BlockGeometry
+
+ M, N, K = 2,3,5
+ ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K)
+
+ u = ig.allocate('random_int')
+ b = ig.allocate('random_int')
+
+ BG = BlockGeometry(ig, ig)
+ U = BG.allocate('random_int')
+
+ f2 = 0.5 * L2NormSquared(b=b)
+ f1 = 30 * MixedL21Norm()
+ tau = 0.355
+
+ res_no_out1 = f1.proximal_conjugate(U, tau)
+ res_no_out2 = f2.proximal_conjugate(u, tau)
+
+
+# print( " ######## with out ######## ")
+ res_out1 = BG.allocate()
+ res_out2 = ig.allocate()
+
+ f1.proximal_conjugate(U, tau, out = res_out1)
+ f2.proximal_conjugate(u, tau, out = res_out2)
+
+
+ numpy.testing.assert_array_almost_equal(res_no_out1[0].as_array(), \
+ res_out1[0].as_array(), decimal=4)
+
+ numpy.testing.assert_array_almost_equal(res_no_out2.as_array(), \
+ res_out2.as_array(), decimal=4)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/ZeroFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ZeroFunction.py new file mode 100644 index 0000000..a019815 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ZeroFunction.py @@ -0,0 +1,55 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.optimisation.functions import Function +from ccpi.framework import BlockDataContainer + +class ZeroFunction(Function): + + ''' ZeroFunction: f(x) = 0 + + + ''' + + def __init__(self): + super(ZeroFunction, self).__init__() + + def __call__(self,x): + return 0 + + def convex_conjugate(self, x): + + ''' This is the support function sup <x, x^{*}> which in fact is the + indicator function for the set = {0} + So 0 if x=0, or inf if x neq 0 + ''' + return x.maximum(0).sum() + + + def proximal(self, x, tau, out=None): + if out is None: + return x.copy() + else: + out.fill(x) + + def proximal_conjugate(self, x, tau, out = None): + if out is None: + return 0 + else: + return 0 diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/BlockOperator.py new file mode 100644 index 0000000..c8bd546 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/BlockOperator.py @@ -0,0 +1,417 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 14 12:36:40 2019 + +@author: ofn77899 +""" +#from ccpi.optimisation.ops import Operator +import numpy +from numbers import Number +import functools +from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer, DataContainer +from ccpi.optimisation.operators import Operator, LinearOperator +from ccpi.optimisation.operators.BlockScaledOperator import BlockScaledOperator +from ccpi.framework import BlockGeometry + +class BlockOperator(Operator): + '''A Block matrix containing Operators + + The Block Framework is a generic strategy to treat variational problems in the + following form: + + .. math:: + + min Regulariser + Fidelity + + + BlockOperators have a generic shape M x N, and when applied on an + Nx1 BlockDataContainer, will yield and Mx1 BlockDataContainer. + Notice: BlockDatacontainer are only allowed to have the shape of N x 1, with + N rows and 1 column. + + User may specify the shape of the block, by default is a row vector + + Operators in a Block are required to have the same domain column-wise and the + same range row-wise. + ''' + __array_priority__ = 1 + def __init__(self, *args, **kwargs): + ''' + Class creator + + Note: + Do not include the `self` parameter in the ``Args`` section. + + Args: + :param: vararg (Operator): Operators in the block. + :param: shape (:obj:`tuple`, optional): If shape is passed the Operators in + vararg are considered input in a row-by-row fashion. + Shape and number of Operators must match. + + Example: + BlockOperator(op0,op1) results in a row block + BlockOperator(op0,op1,shape=(1,2)) results in a column block + ''' + self.operators = args + shape = kwargs.get('shape', None) + if shape is None: + shape = (len(args),1) + self.shape = shape + n_elements = functools.reduce(lambda x,y: x*y, shape, 1) + if len(args) != n_elements: + raise ValueError( + 'Dimension and size do not match: expected {} got {}' + .format(n_elements,len(args))) + # test if operators are compatible + if not self.column_wise_compatible(): + raise ValueError('Operators in each column must have the same domain') + if not self.row_wise_compatible(): + raise ValueError('Operators in each row must have the same range') + + def column_wise_compatible(self): + '''Operators in a Block should have the same domain per column''' + rows, cols = self.shape + compatible = True + for col in range(cols): + column_compatible = True + for row in range(1,rows): + dg0 = self.get_item(row-1,col).domain_geometry() + dg1 = self.get_item(row,col).domain_geometry() + column_compatible = dg0.__dict__ == dg1.__dict__ and column_compatible + compatible = compatible and column_compatible + return compatible + + def row_wise_compatible(self): + '''Operators in a Block should have the same range per row''' + rows, cols = self.shape + compatible = True + for row in range(rows): + row_compatible = True + for col in range(1,cols): + dg0 = self.get_item(row,col-1).range_geometry() + dg1 = self.get_item(row,col).range_geometry() + row_compatible = dg0.__dict__ == dg1.__dict__ and row_compatible + compatible = compatible and row_compatible + return compatible + + def get_item(self, row, col): + '''returns the Operator at specified row and col''' + if row > self.shape[0]: + raise ValueError('Requested row {} > max {}'.format(row, self.shape[0])) + if col > self.shape[1]: + raise ValueError('Requested col {} > max {}'.format(col, self.shape[1])) + + index = row*self.shape[1]+col + return self.operators[index] + + def norm(self): + norm = [op.norm()**2 for op in self.operators] + return numpy.sqrt(sum(norm)) + + def direct(self, x, out=None): + '''Direct operation for the BlockOperator + + BlockOperator work on BlockDataContainer, but they will work on DataContainers + and inherited classes by simple wrapping the input in a BlockDataContainer of shape (1,1) + ''' + + if not isinstance (x, BlockDataContainer): + x_b = BlockDataContainer(x) + else: + x_b = x + shape = self.get_output_shape(x_b.shape) + res = [] + + if out is None: + + for row in range(self.shape[0]): + for col in range(self.shape[1]): + if col == 0: + prod = self.get_item(row,col).direct(x_b.get_item(col)) + else: + prod += self.get_item(row,col).direct(x_b.get_item(col)) + res.append(prod) + return BlockDataContainer(*res, shape=shape) + + else: + + tmp = self.range_geometry().allocate() + for row in range(self.shape[0]): + for col in range(self.shape[1]): + if col == 0: + self.get_item(row,col).direct( + x_b.get_item(col), + out=out.get_item(row)) + else: + a = out.get_item(row) + self.get_item(row,col).direct( + x_b.get_item(col), + out=tmp.get_item(row)) + a += tmp.get_item(row) + + def adjoint(self, x, out=None): + '''Adjoint operation for the BlockOperator + + BlockOperator may contain both LinearOperator and Operator + This method exists in BlockOperator as it is not known what type of + Operator it will contain. + + BlockOperator work on BlockDataContainer, but they will work on DataContainers + and inherited classes by simple wrapping the input in a BlockDataContainer of shape (1,1) + + Raises: ValueError if the contained Operators are not linear + ''' + if not self.is_linear(): + raise ValueError('Not all operators in Block are linear.') + if not isinstance (x, BlockDataContainer): + x_b = BlockDataContainer(x) + else: + x_b = x + shape = self.get_output_shape(x_b.shape, adjoint=True) + if out is None: + res = [] + for col in range(self.shape[1]): + for row in range(self.shape[0]): + if row == 0: + prod = self.get_item(row, col).adjoint(x_b.get_item(row)) + else: + prod += self.get_item(row, col).adjoint(x_b.get_item(row)) + res.append(prod) + if self.shape[1]==1: + return ImageData(*res) + else: + return BlockDataContainer(*res, shape=shape) + else: + #tmp = self.domain_geometry().allocate() + + for col in range(self.shape[1]): + for row in range(self.shape[0]): + if row == 0: + if issubclass(out.__class__, DataContainer): + self.get_item(row, col).adjoint( + x_b.get_item(row), + out=out) + else: + op = self.get_item(row,col) + self.get_item(row, col).adjoint( + x_b.get_item(row), + out=out.get_item(col)) + else: + if issubclass(out.__class__, DataContainer): + out += self.get_item(row,col).adjoint( + x_b.get_item(row)) + else: + a = out.get_item(col) + a += self.get_item(row,col).adjoint( + x_b.get_item(row), + ) + def is_linear(self): + '''returns whether all the elements of the BlockOperator are linear''' + return functools.reduce(lambda x, y: x and y.is_linear(), self.operators, True) + + def get_output_shape(self, xshape, adjoint=False): + '''returns the shape of the output BlockDataContainer + + A(N,M) direct u(M,1) -> N,1 + A(N,M)^T adjoint u(N,1) -> M,1 + ''' + rows , cols = self.shape + xrows, xcols = xshape + if xcols != 1: + raise ValueError('BlockDataContainer cannot have more than 1 column') + if adjoint: + if rows != xrows: + raise ValueError('Incompatible shapes {} {}'.format(self.shape, xshape)) + return (cols,xcols) + if cols != xrows: + raise ValueError('Incompatible shapes {} {}'.format((rows,cols), xshape)) + return (rows,xcols) + + def __rmul__(self, scalar): + '''Defines the left multiplication with a scalar + + Args: scalar (number or iterable containing numbers): + + Returns: a block operator with Scaled Operators inside''' + if isinstance (scalar, list) or isinstance(scalar, tuple) or \ + isinstance(scalar, numpy.ndarray): + if len(scalar) != len(self.operators): + raise ValueError('dimensions of scalars and operators do not match') + scalars = scalar + else: + scalars = [scalar for _ in self.operators] + # create a list of ScaledOperator-s + ops = [ v * op for v,op in zip(scalars, self.operators)] + #return BlockScaledOperator(self, scalars ,shape=self.shape) + return type(self)(*ops, shape=self.shape) + @property + def T(self): + '''Return the transposed of self + + input in a row-by-row''' + newshape = (self.shape[1], self.shape[0]) + oplist = [] + for col in range(newshape[1]): + for row in range(newshape[0]): + oplist.append(self.get_item(col,row)) + return type(self)(*oplist, shape=newshape) + + def domain_geometry(self): + '''returns the domain of the BlockOperator + + If the shape of the BlockOperator is (N,1) the domain is a ImageGeometry or AcquisitionGeometry. + Otherwise it is a BlockGeometry. + ''' + if self.shape[1] == 1: + # column BlockOperator + return self.get_item(0,0).domain_geometry() + else: + # get the geometries column wise + # we need only the geometries from the first row + # since it is compatible from __init__ + tmp = [] + for i in range(self.shape[1]): + tmp.append(self.get_item(0,i).domain_geometry()) + return BlockGeometry(*tmp) + + #shape = (self.shape[0], 1) + #return BlockGeometry(*[el.domain_geometry() for el in self.operators], + # shape=self.shape) + + def range_geometry(self): + '''returns the range of the BlockOperator''' + + tmp = [] + for i in range(self.shape[0]): + tmp.append(self.get_item(i,0).range_geometry()) + return BlockGeometry(*tmp) + + + #shape = (self.shape[1], 1) + #return BlockGeometry(*[el.range_geometry() for el in self.operators], + # shape=shape) + + def sum_abs_row(self): + + res = [] + for row in range(self.shape[0]): + for col in range(self.shape[1]): + if col == 0: + prod = self.get_item(row,col).sum_abs_row() + else: + prod += self.get_item(row,col).sum_abs_row() + res.append(prod) + + if self.shape[1]==1: + tmp = sum(res) + return ImageData(tmp) + else: + + return BlockDataContainer(*res) + + def sum_abs_col(self): + + res = [] + for row in range(self.shape[0]): + for col in range(self.shape[1]): + if col == 0: + prod = self.get_item(row, col).sum_abs_col() + else: + prod += self.get_item(row, col).sum_abs_col() + res.append(prod) + + return BlockDataContainer(*res) + + + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + from ccpi.optimisation.operators import Gradient, Identity, \ + SparseFiniteDiff, SymmetrizedGradient, ZeroOperator + + + M, N = 4, 3 + ig = ImageGeometry(M, N) + arr = ig.allocate('random_int') + + G = Gradient(ig) + Id = Identity(ig) + + B = BlockOperator(G, Id) + + print(B.sum_abs_row()) +# + Gx = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + Gy = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + + d1 = abs(Gx.matrix()).toarray().sum(axis=0) + d2 = abs(Gy.matrix()).toarray().sum(axis=0) + d3 = abs(Id.matrix()).toarray().sum(axis=0) + + + d_res = numpy.reshape(d1 + d2 + d3, ig.shape, 'F') + + print(d_res) +# + z1 = abs(Gx.matrix()).toarray().sum(axis=1) + z2 = abs(Gy.matrix()).toarray().sum(axis=1) + z3 = abs(Id.matrix()).toarray().sum(axis=1) +# + z_res = BlockDataContainer(BlockDataContainer(ImageData(numpy.reshape(z2, ig.shape, 'F')),\ + ImageData(numpy.reshape(z1, ig.shape, 'F'))),\ + ImageData(numpy.reshape(z3, ig.shape, 'F'))) +# + ttt = B.sum_abs_col() +# + #TODO this is not working +# numpy.testing.assert_array_almost_equal(z_res[0][0].as_array(), ttt[0][0].as_array(), decimal=4) +# numpy.testing.assert_array_almost_equal(z_res[0][1].as_array(), ttt[0][1].as_array(), decimal=4) +# numpy.testing.assert_array_almost_equal(z_res[1].as_array(), ttt[1].as_array(), decimal=4) + + + u = ig.allocate('random_int') + + z1 = B.direct(u) + res = B.range_geometry().allocate() + + B.direct(u, out = res) + + + + ########################################################################### + # Block Operator for TGV reconstruction + + M, N = 2,3 + ig = ImageGeometry(M, N) + ag = ig + + op11 = Gradient(ig) + op12 = Identity(op11.range_geometry()) + + op22 = SymmetrizedGradient(op11.domain_geometry()) + + op21 = ZeroOperator(ig, op22.range_geometry()) + + + op31 = Identity(ig, ag) + op32 = ZeroOperator(op22.domain_geometry(), ag) + + operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) + + z1 = operator.domain_geometry() + z2 = operator.range_geometry() + + print(z1.shape) + print(z2.shape) + + + + + + + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/FiniteDifferenceOperator.py new file mode 100644 index 0000000..f459334 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -0,0 +1,372 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:51:17 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import LinearOperator +from ccpi.framework import ImageData, BlockDataContainer +import numpy as np + +class FiniteDiff(LinearOperator): + + # Works for Neum/Symmetric & periodic boundary conditions + # TODO add central differences??? + # TODO not very well optimised, too many conditions + # TODO add discretisation step, should get that from imageGeometry + + # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x'] + # Grad_order = ['channels', 'direction_y', 'direction_x'] + # Grad_order = ['direction_z', 'direction_y', 'direction_x'] + # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x'] + + def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'): + '''''' + super(FiniteDiff, self).__init__() + '''FIXME: domain and range should be geometries''' + self.gm_domain = gm_domain + self.gm_range = gm_range + + self.direction = direction + self.bnd_cond = bnd_cond + + # Domain Geometry = Range Geometry if not stated + if self.gm_range is None: + self.gm_range = self.gm_domain + # check direction and "length" of geometry + if self.direction + 1 > len(self.gm_domain.shape): + raise ValueError('Gradient directions more than geometry domain') + + #self.voxel_size = kwargs.get('voxel_size',1) + # this wrongly assumes a homogeneous voxel size +# self.voxel_size = self.gm_domain.voxel_size_x + + + def direct(self, x, out=None): + + x_asarr = x.as_array() + x_sz = len(x.shape) + + if out is None: + out = np.zeros_like(x_asarr) + else: + out = out.as_array() + out[:]=0 + + + + ######################## Direct for 2D ############################### + if x_sz == 2: + + if self.direction == 1: + + np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = out[:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0], x_asarr[:,-1], out = out[:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 0: + + np.subtract( x_asarr[1:], x_asarr[0:-1], out = out[0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:], x_asarr[-1,:], out = out[-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Direct for 3D ############################### + elif x_sz == 3: + + if self.direction == 0: + + np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = out[0:-1,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = out[-1,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + + np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = out[:,0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = out[:,-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + + if self.direction == 2: + + np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = out[:,:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = out[:,:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Direct for 4D ############################### + elif x_sz == 4: + + if self.direction == 0: + np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = out[0:-1,:,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = out[-1,:,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = out[:,0:-1,:,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = out[:,-1,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = out[:,:,0:-1,:] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = out[:,:,-1,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 3: + np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = out[:,:,:,0:-1] ) + + if self.bnd_cond == 'Neumann': + pass + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = out[:,:,:,-1] ) + else: + raise ValueError('No valid boundary conditions') + + else: + raise NotImplementedError + +# res = out #/self.voxel_size + return type(x)(out) + + + def adjoint(self, x, out=None): + + x_asarr = x.as_array() + x_sz = len(x.shape) + + if out is None: + out = np.zeros_like(x_asarr) + else: + out = out.as_array() + out[:]=0 + + + ######################## Adjoint for 2D ############################### + if x_sz == 2: + + if self.direction == 1: + + np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = out[:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0], 0, out = out[:,0] ) + np.subtract( -x_asarr[:,-2], 0, out = out[:,-1] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0], x_asarr[:,-1], out = out[:,0] ) + + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 0: + + np.subtract( x_asarr[1:,:], x_asarr[0:-1,:], out = out[1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:], 0, out = out[0,:] ) + np.subtract( -x_asarr[-2,:], 0, out = out[-1,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:], x_asarr[-1,:], out = out[0,:] ) + + else: + raise ValueError('No valid boundary conditions') + + ######################## Adjoint for 3D ############################### + elif x_sz == 3: + + if self.direction == 0: + + np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = out[1:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:,:], 0, out = out[0,:,:] ) + np.subtract( -x_asarr[-2,:,:], 0, out = out[-1,:,:] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = out[0,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = out[:,1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0,:], 0, out = out[:,0,:] ) + np.subtract( -x_asarr[:,-2,:], 0, out = out[:,-1,:] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = out[:,0,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = out[:,:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,0], 0, out = out[:,:,0] ) + np.subtract( -x_asarr[:,:,-2], 0, out = out[:,:,-1] ) + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = out[:,:,0] ) + else: + raise ValueError('No valid boundary conditions') + + ######################## Adjoint for 4D ############################### + elif x_sz == 4: + + if self.direction == 0: + np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = out[1:,:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[0,:,:,:], 0, out = out[0,:,:,:] ) + np.subtract( -x_asarr[-2,:,:,:], 0, out = out[-1,:,:,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = out[0,:,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 1: + np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = out[:,1:,:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,0,:,:], 0, out = out[:,0,:,:] ) + np.subtract( -x_asarr[:,-2,:,:], 0, out = out[:,-1,:,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = out[:,0,:,:] ) + else: + raise ValueError('No valid boundary conditions') + + + if self.direction == 2: + np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = out[:,:,1:,:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,0,:], 0, out = out[:,:,0,:] ) + np.subtract( -x_asarr[:,:,-2,:], 0, out = out[:,:,-1,:] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = out[:,:,0,:] ) + else: + raise ValueError('No valid boundary conditions') + + if self.direction == 3: + np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = out[:,:,:,1:] ) + + if self.bnd_cond == 'Neumann': + np.subtract( x_asarr[:,:,:,0], 0, out = out[:,:,:,0] ) + np.subtract( -x_asarr[:,:,:,-2], 0, out = out[:,:,:,-1] ) + + elif self.bnd_cond == 'Periodic': + np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = out[:,:,:,0] ) + else: + raise ValueError('No valid boundary conditions') + + else: + raise NotImplementedError + + out *= -1 #/self.voxel_size + return type(x)(out) + + def range_geometry(self): + '''Returns the range geometry''' + return self.gm_range + + def domain_geometry(self): + '''Returns the domain geometry''' + return self.gm_domain + + def norm(self): + x0 = self.gm_domain.allocate('random_int') + self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0) + return self.s1 + + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + import numpy + + N, M = 2, 3 + + ig = ImageGeometry(N, M) + + + FD = FiniteDiff(ig, direction = 1, bnd_cond = 'Neumann') + u = FD.domain_geometry().allocate('random_int') + + res = FD.domain_geometry().allocate() + res1 = FD.range_geometry().allocate() + FD.direct(u, out=res) + + z = FD.direct(u) +# print(z.as_array(), res.as_array()) + + for i in range(10): +# + z1 = FD.direct(u) + FD.direct(u, out=res) + + u = ig.allocate('random_int') + res = u + z1 = u + numpy.testing.assert_array_almost_equal(z1.as_array(), \ + res.as_array(), decimal=4) + +# print(z1.as_array(), res.as_array()) + z2 = FD.adjoint(z1) + FD.adjoint(z1, out=res1) + numpy.testing.assert_array_almost_equal(z2.as_array(), \ + res1.as_array(), decimal=4) + + + + + + + + +# w = G.range_geometry().allocate('random_int') + + + +
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py new file mode 100644 index 0000000..e0b8a32 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py @@ -0,0 +1,242 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:50:04 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator +from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer +import numpy +from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff + +#%% + +class Gradient(LinearOperator): + + def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): + + super(Gradient, self).__init__() + + self.gm_domain = gm_domain # Domain of Grad Operator + + self.correlation = kwargs.get('correlation','Space') + + if self.correlation=='Space': + if self.gm_domain.channels>1: + self.gm_range = BlockGeometry(*[self.gm_domain for _ in range(self.gm_domain.length-1)] ) + self.ind = numpy.arange(1,self.gm_domain.length) + else: + self.gm_range = BlockGeometry(*[self.gm_domain for _ in range(self.gm_domain.length) ] ) + self.ind = numpy.arange(self.gm_domain.length) + elif self.correlation=='SpaceChannels': + if self.gm_domain.channels>1: + self.gm_range = BlockGeometry(*[self.gm_domain for _ in range(self.gm_domain.length)]) + self.ind = range(self.gm_domain.length) + else: + raise ValueError('No channels to correlate') + + self.bnd_cond = bnd_cond + + # Call FiniteDiff operator + + self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond) + + + def direct(self, x, out=None): + + + if out is not None: + + for i in range(self.gm_range.shape[0]): + self.FD.direction = self.ind[i] + self.FD.direct(x, out = out[i]) + else: + tmp = self.gm_range.allocate() + for i in range(tmp.shape[0]): + self.FD.direction=self.ind[i] + tmp.get_item(i).fill(self.FD.direct(x)) + return tmp + + def adjoint(self, x, out=None): + + if out is not None: + + tmp = self.gm_domain.allocate() + for i in range(x.shape[0]): + self.FD.direction=self.ind[i] + self.FD.adjoint(x.get_item(i), out = tmp) + if i == 0: + out.fill(tmp) + else: + out += tmp + else: + tmp = self.gm_domain.allocate() + for i in range(x.shape[0]): + self.FD.direction=self.ind[i] + + tmp += self.FD.adjoint(x.get_item(i)) + return tmp + + + def domain_geometry(self): + return self.gm_domain + + def range_geometry(self): + return self.gm_range + + def norm(self): + + x0 = self.gm_domain.allocate('random') + self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0) + return self.s1 + + def __rmul__(self, scalar): + return ScaledOperator(self, scalar) + + ########################################################################### + ############### For preconditioning ###################################### + ########################################################################### + def matrix(self): + + tmp = self.gm_range.allocate() + + mat = [] + for i in range(tmp.shape[0]): + + spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond) + mat.append(spMat.matrix()) + + return BlockDataContainer(*mat) + + + def sum_abs_row(self): + + tmp = self.gm_range.allocate() + res = self.gm_domain.allocate() + for i in range(tmp.shape[0]): + spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond) + res += spMat.sum_abs_row() + return res + + def sum_abs_col(self): + + tmp = self.gm_range.allocate() + res = [] + for i in range(tmp.shape[0]): + spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond) + res.append(spMat.sum_abs_col()) + return BlockDataContainer(*res) + + +if __name__ == '__main__': + + + from ccpi.optimisation.operators import Identity, BlockOperator + + + M, N = 2, 3 + ig = ImageGeometry(M, N) + arr = ig.allocate('random_int' ) + + # check direct of Gradient and sparse matrix + G = Gradient(ig) + G_sp = G.matrix() + + res1 = G.direct(arr) + res1y = numpy.reshape(G_sp[0].toarray().dot(arr.as_array().flatten('F')), ig.shape, 'F') + + print(res1[0].as_array()) + print(res1y) + + res1x = numpy.reshape(G_sp[1].toarray().dot(arr.as_array().flatten('F')), ig.shape, 'F') + + print(res1[1].as_array()) + print(res1x) + + #check sum abs row + conc_spmat = numpy.abs(numpy.concatenate( (G_sp[0].toarray(), G_sp[1].toarray() ))) + print(numpy.reshape(conc_spmat.sum(axis=0), ig.shape, 'F')) + print(G.sum_abs_row().as_array()) + + print(numpy.reshape(conc_spmat.sum(axis=1), ((2,) + ig.shape), 'F')) + + print(G.sum_abs_col()[0].as_array()) + print(G.sum_abs_col()[1].as_array()) + + # Check Blockoperator sum abs col and row + + op1 = Gradient(ig) + op2 = Identity(ig) + + B = BlockOperator( op1, op2) + + Brow = B.sum_abs_row() + Bcol = B.sum_abs_col() + + concB = numpy.concatenate( (numpy.abs(numpy.concatenate( (G_sp[0].toarray(), G_sp[1].toarray() ))), op2.matrix().toarray())) + + print(numpy.reshape(concB.sum(axis=0), ig.shape, 'F')) + print(Brow.as_array()) + + print(numpy.reshape(concB.sum(axis=1)[0:12], ((2,) + ig.shape), 'F')) + print(Bcol[1].as_array()) + + +# print(numpy.concatene(G_sp[0].toarray()+ )) +# print(G_sp[1].toarray()) +# +# d1 = G.sum_abs_row() +# print(d1.as_array()) +# +# d2 = G_neum.sum_abs_col() +## print(d2) +# +# +# ########################################################### + a = BlockDataContainer( BlockDataContainer(arr, arr), arr) + b = BlockDataContainer( BlockDataContainer(arr+5, arr+3), arr+2) + c = a/b + + print(c[0][0].as_array(), (arr/(arr+5)).as_array()) + print(c[0][1].as_array(), (arr/(arr+3)).as_array()) + print(c[1].as_array(), (arr/(arr+2)).as_array()) + + + a1 = BlockDataContainer( arr, BlockDataContainer(arr, arr)) +# +# c1 = arr + a +# c2 = arr + a +# c2 = a1 + arr + + from ccpi.framework import ImageGeometry +# from ccpi.optimisation.operators import Gradient +# + N, M = 2, 3 +# + ig = ImageGeometry(N, M) +# + G = Gradient(ig) +# + u = G.domain_geometry().allocate('random_int') + w = G.range_geometry().allocate('random_int') + + + print( "################ without out #############") + + print( (G.direct(u)*w).sum(), (u*G.adjoint(w)).sum() ) + + + print( "################ with out #############") + + res = G.range_geometry().allocate() + res1 = G.domain_geometry().allocate() + G.direct(u, out = res) + G.adjoint(w, out = res1) + + print( (res*w).sum(), (u*res1).sum() ) + + + + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperator.py new file mode 100644 index 0000000..f304cc6 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperator.py @@ -0,0 +1,67 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 15:57:52 2019
+
+@author: ofn77899
+"""
+
+from ccpi.optimisation.operators import Operator
+from ccpi.framework import ImageGeometry
+import numpy
+
+
+class LinearOperator(Operator):
+ '''A Linear Operator that maps from a space X <-> Y'''
+ def __init__(self):
+ super(LinearOperator, self).__init__()
+ def is_linear(self):
+ '''Returns if the operator is linear'''
+ return True
+ def adjoint(self,x, out=None):
+ '''returns the adjoint/inverse operation
+
+ only available to linear operators'''
+ raise NotImplementedError
+
+ @staticmethod
+ def PowerMethod(operator, iterations, x_init=None):
+ '''Power method to calculate iteratively the Lipschitz constant'''
+
+ # Initialise random
+ if x_init is None:
+ x0 = operator.domain_geometry().allocate(ImageGeometry.RANDOM_INT)
+ else:
+ x0 = x_init.copy()
+
+ x1 = operator.domain_geometry().allocate()
+ y_tmp = operator.range_geometry().allocate()
+ s = numpy.zeros(iterations)
+ # Loop
+ for it in numpy.arange(iterations):
+ operator.direct(x0,out=y_tmp)
+ operator.adjoint(y_tmp,out=x1)
+ x1norm = x1.norm()
+ s[it] = x1.dot(x0) / x0.squared_norm()
+ x1.multiply((1.0/x1norm), out=x0)
+ return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
+
+ @staticmethod
+ def PowerMethodNonsquare(op,numiters , x_init=None):
+ '''Power method to calculate iteratively the Lipschitz constant'''
+
+ if x_init is None:
+ x0 = op.allocate_direct()
+ x0.fill(numpy.random.randn(*x0.shape))
+ else:
+ x0 = x_init.copy()
+
+ s = numpy.zeros(numiters)
+ # Loop
+ for it in numpy.arange(numiters):
+ x1 = op.adjoint(op.direct(x0))
+ x1norm = x1.norm()
+ s[it] = (x1*x0).sum() / (x0.squared_norm())
+ x0 = (1.0/x1norm)*x1
+ return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
+
+
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperatorMatrix.py new file mode 100644 index 0000000..62e22e0 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperatorMatrix.py @@ -0,0 +1,51 @@ +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 +from numbers import Number +from ccpi.optimisation.operators import LinearOperator +class LinearOperatorMatrix(LinearOperator): + def __init__(self,A): + self.A = A + self.s1 = None # Largest singular value, initially unknown + super(LinearOperatorMatrix, self).__init__() + + 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, 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 + + def get_max_sing_val(self): + # If unknown, compute and store. If known, simply return it. + if self.s1 is None: + self.s1 = svds(self.A,1,return_singular_vectors=False)[0] + 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/build/lib/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SymmetrizedGradientOperator.py new file mode 100644 index 0000000..243809b --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -0,0 +1,244 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 1 22:53:55 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import Gradient, Operator, LinearOperator, ScaledOperator +from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer +import numpy +from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff + + +class SymmetrizedGradient(Gradient): + + ''' Symmetrized Gradient, denoted by E: V --> W + where V is the Range of the Gradient Operator + and W is the Range of the Symmetrized Gradient. + ''' + + + def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): + + super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) + + ''' + Domain of SymGrad is the Range of Gradient + ''' + + self.gm_domain = self.gm_range + self.bnd_cond = bnd_cond + + self.channels = self.gm_range.get_item(0).channels + + tmp_gm = len(self.gm_domain.geometries)*self.gm_domain.geometries + + self.gm_range = BlockGeometry(*tmp_gm) + + self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond) + + if self.gm_domain.shape[0]==2: + self.order_ind = [0,2,1,3] + else: + self.order_ind = [0,3,6,1,4,7,2,5,8] + + + def direct(self, x, out=None): + + if out is None: + + tmp = [] + for i in range(self.gm_domain.shape[0]): + for j in range(x.shape[0]): + self.FD.direction = i + tmp.append(self.FD.adjoint(x.get_item(j))) + + tmp1 = [tmp[i] for i in self.order_ind] + + res = [0.5 * sum(x) for x in zip(tmp, tmp1)] + + return BlockDataContainer(*res) + + else: + + ind = 0 + for i in range(self.gm_domain.shape[0]): + for j in range(x.shape[0]): + self.FD.direction = i + self.FD.adjoint(x.get_item(j), out=out[ind]) + ind+=1 + out1 = BlockDataContainer(*[out[i] for i in self.order_ind]) + out.fill( 0.5 * (out + out1) ) + + + def adjoint(self, x, out=None): + + if out is None: + + tmp = [None]*self.gm_domain.shape[0] + i = 0 + + for k in range(self.gm_domain.shape[0]): + tmp1 = 0 + for j in range(self.gm_domain.shape[0]): + self.FD.direction = j + tmp1 += self.FD.direct(x[i]) + i+=1 + tmp[k] = tmp1 + return BlockDataContainer(*tmp) + + + else: + + tmp = self.gm_domain.allocate() + i = 0 + for k in range(self.gm_domain.shape[0]): + tmp1 = 0 + for j in range(self.gm_domain.shape[0]): + self.FD.direction = j + self.FD.direct(x[i], out=tmp[j]) + i+=1 + tmp1+=tmp[j] + out[k].fill(tmp1) +# tmp = self.adjoint(x) +# out.fill(tmp) + + + def domain_geometry(self): + return self.gm_domain + + def range_geometry(self): + return self.gm_range + + def norm(self): + +# TODO need dot method for BlockDataContainer +# return numpy.sqrt(4*self.gm_domain.shape[0]) + +# x0 = self.gm_domain.allocate('random') + self.s1, sall, svec = LinearOperator.PowerMethod(self, 50) + return self.s1 + + + +if __name__ == '__main__': + + ########################################################################### + ## Symmetrized Gradient Tests + from ccpi.framework import DataContainer + from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff + import numpy as np + + N, M = 2, 3 + K = 2 + C = 2 + + ig1 = ImageGeometry(N, M) + ig2 = ImageGeometry(N, M, channels=C) + + E1 = SymmetrizedGradient(ig1, correlation = 'Space', bnd_cond='Neumann') + + try: + E1 = SymmetrizedGradient(ig1, correlation = 'SpaceChannels', bnd_cond='Neumann') + except: + print("No Channels to correlate") + + E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Neumann') + + print(E1.domain_geometry().shape, E1.range_geometry().shape) + print(E2.domain_geometry().shape, E2.range_geometry().shape) + + #check Linear operator property + + u1 = E1.domain_geometry().allocate('random_int') + u2 = E2.domain_geometry().allocate('random_int') + + # Need to allocate random_int at the Range of SymGradient + + #a1 = ig1.allocate('random_int') + #a2 = ig1.allocate('random_int') + #a3 = ig1.allocate('random_int') + + #a4 = ig1.allocate('random_int') + #a5 = ig1.allocate('random_int') + #a6 = ig1.allocate('random_int') + + # TODO allocate has to create this symmetry by default!!!!! + #w1 = BlockDataContainer(*[a1, a2, \ + # a2, a3]) + w1 = E1.range_geometry().allocate('random_int',symmetry=True) + + LHS = (E1.direct(u1) * w1).sum() + RHS = (u1 * E1.adjoint(w1)).sum() + + numpy.testing.assert_equal(LHS, RHS) + + u2 = E2.gm_domain.allocate('random_int') + + #aa1 = ig2.allocate('random_int') + #aa2 = ig2.allocate('random_int') + #aa3 = ig2.allocate('random_int') + #aa4 = ig2.allocate('random_int') + #aa5 = ig2.allocate('random_int') + #aa6 = ig2.allocate('random_int') + + #w2 = BlockDataContainer(*[aa1, aa2, aa3, \ + # aa2, aa4, aa5, \ + # aa3, aa5, aa6]) + w2 = E2.range_geometry().allocate('random_int',symmetry=True) + + + LHS1 = (E2.direct(u2) * w2).sum() + RHS1 = (u2 * E2.adjoint(w2)).sum() + + numpy.testing.assert_equal(LHS1, RHS1) + + out = E1.range_geometry().allocate() + E1.direct(u1, out=out) + a1 = E1.direct(u1) + numpy.testing.assert_array_equal(a1[0].as_array(), out[0].as_array()) + numpy.testing.assert_array_equal(a1[1].as_array(), out[1].as_array()) + numpy.testing.assert_array_equal(a1[2].as_array(), out[2].as_array()) + numpy.testing.assert_array_equal(a1[3].as_array(), out[3].as_array()) + + + out1 = E1.domain_geometry().allocate() + E1.adjoint(w1, out=out1) + b1 = E1.adjoint(w1) + + LHS_out = (out * w1).sum() + RHS_out = (u1 * out1).sum() + print(LHS_out, RHS_out) + + + out2 = E2.range_geometry().allocate() + E2.direct(u2, out=out2) + a2 = E2.direct(u2) + + out21 = E2.domain_geometry().allocate() + E2.adjoint(w2, out=out21) + b2 = E2.adjoint(w2) + + LHS_out = (out2 * w2).sum() + RHS_out = (u2 * out21).sum() + print(LHS_out, RHS_out) + + + out = E1.range_geometry().allocate() + E1.direct(u1, out=out) + E1.adjoint(out, out=out1) + + print(E1.norm()) + print(E2.norm()) + + + + + + +# +# +# +#
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/ZeroOperator.py new file mode 100644 index 0000000..8168f0b --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/ZeroOperator.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 6 19:25:53 2019 + +@author: evangelos +""" + +import numpy as np +from ccpi.framework import ImageData +from ccpi.optimisation.operators import LinearOperator + +class ZeroOperator(LinearOperator): + + def __init__(self, gm_domain, gm_range=None): + + super(ZeroOperator, self).__init__() + + self.gm_domain = gm_domain + self.gm_range = gm_range + if self.gm_range is None: + self.gm_range = self.gm_domain + + + def direct(self,x,out=None): + if out is None: + return self.gm_range.allocate() + else: + out.fill(self.gm_range.allocate()) + + def adjoint(self,x, out=None): + if out is None: + return self.gm_domain.allocate() + else: + out.fill(self.gm_domain.allocate()) + + def norm(self): + return 0 + + def domain_geometry(self): + return self.gm_domain + + def range_geometry(self): + return self.gm_range
\ No newline at end of file diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/__init__.py new file mode 100644 index 0000000..23222d4 --- /dev/null +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/__init__.py @@ -0,0 +1,23 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 15:56:27 2019
+
+@author: ofn77899
+"""
+
+from .Operator import Operator
+from .LinearOperator import LinearOperator
+from .ScaledOperator import ScaledOperator
+from .BlockOperator import BlockOperator
+from .BlockScaledOperator import BlockScaledOperator
+
+from .SparseFiniteDiff import SparseFiniteDiff
+from .ShrinkageOperator import ShrinkageOperator
+
+from .FiniteDifferenceOperator import FiniteDiff
+from .GradientOperator import Gradient
+from .SymmetrizedGradientOperator import SymmetrizedGradient
+from .IdentityOperator import Identity
+from .ZeroOperator import ZeroOperator
+from .LinearOperatorMatrix import LinearOperatorMatrix
+
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index 7236e0e..dbe7d0a 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -762,12 +762,16 @@ class DataContainer(object): def norm(self): '''return the euclidean norm of the DataContainer viewed as a vector''' return numpy.sqrt(self.squared_norm()) + + def dot(self, other, *args, **kwargs): '''return the inner product of 2 DataContainers viewed as vectors''' method = kwargs.get('method', 'reduce') + if method not in ['numpy','reduce']: raise ValueError('dot: specified method not valid. Expecting numpy or reduce got {} '.format( method)) + if self.shape == other.shape: # return (self*other).sum() if method == 'numpy': @@ -784,6 +788,7 @@ class DataContainer(object): else: raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape)) + diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py index e65bc89..4d4843c 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 + class CGLS(Algorithm): '''Conjugate Gradient Least Squares algorithm @@ -83,4 +84,4 @@ class CGLS(Algorithm): self.d = s + beta*self.d def update_objective(self): - self.loss.append(self.r.squared_norm()) + self.loss.append(self.r.squared_norm())
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py index 8ea2b6c..ee51049 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -65,10 +65,10 @@ class FISTA(Algorithm): # initialization if memopt: - self.y = x_init.clone() - self.x_old = x_init.clone() - self.x = x_init.clone() - self.u = x_init.clone() + self.y = x_init.copy() + self.x_old = x_init.copy() + self.x = x_init.copy() + self.u = x_init.copy() else: self.x_old = x_init.copy() self.y = x_init.copy() diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py index 70511bb..8895f3d 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py @@ -19,16 +19,13 @@ class FunctionOperatorComposition(Function): ''' - def __init__(self, operator, function): + def __init__(self, function, operator): super(FunctionOperatorComposition, self).__init__() + self.function = function self.operator = operator - alpha = 1 - - if isinstance (function, ScaledFunction): - alpha = function.scalar - self.L = 2 * alpha * operator.norm()**2 + self.L = function.L * operator.norm()**2 def __call__(self, x): @@ -39,47 +36,57 @@ class FunctionOperatorComposition(Function): ''' - return self.function(self.operator.direct(x)) + return self.function(self.operator.direct(x)) + + def gradient(self, x, out=None): +# + ''' Gradient takes into account the Operator''' + if out is None: + return self.operator.adjoint(self.function.gradient(self.operator.direct(x))) + else: + tmp = self.operator.range_geometry().allocate() + self.operator.direct(x, out=tmp) + self.function.gradient(tmp, out=tmp) + self.operator.adjoint(tmp, out=out) - #TODO do not know if we need it - def call_adjoint(self, x): - return self.function(self.operator.adjoint(x)) + + #TODO do not know if we need it + #def call_adjoint(self, x): + # + # return self.function(self.operator.adjoint(x)) - def convex_conjugate(self, x): - - ''' convex_conjugate does not take into account the Operator''' - return self.function.convex_conjugate(x) - def proximal(self, x, tau, out=None): - - '''proximal does not take into account the Operator''' - if out is None: - return self.function.proximal(x, tau) - else: - self.function.proximal(x, tau, out=out) - + #def convex_conjugate(self, x): + # + # ''' convex_conjugate does not take into account the Operator''' + # return self.function.convex_conjugate(x) - def proximal_conjugate(self, x, tau, out=None): + - ''' proximal conjugate does not take into account the Operator''' - if out is None: - return self.function.proximal_conjugate(x, tau) - else: - self.function.proximal_conjugate(x, tau, out=out) - def gradient(self, x, out=None): - - ''' Gradient takes into account the Operator''' - if out is None: - return self.operator.adjoint( - self.function.gradient(self.operator.direct(x)) - ) - else: - self.operator.adjoint( - self.function.gradient(self.operator.direct(x), - out=out) - ) + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + from ccpi.optimisation.operators import Gradient + from ccpi.optimisation.functions import L2NormSquared + + M, N, K = 2,3 + ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N) + + G = Gradient(ig) + alpha = 0.5 + + f = L2NormSquared() + f_comp = FunctionOperatorComposition(G, alpha * f) + x = ig.allocate('random_int') + print(f_comp.gradient(x).shape + + ) + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition_old.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition_old.py new file mode 100644 index 0000000..70511bb --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition_old.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 8 09:55:36 2019 + +@author: evangelos +""" + +from ccpi.optimisation.functions import Function +from ccpi.optimisation.functions import ScaledFunction + + +class FunctionOperatorComposition(Function): + + ''' Function composition with Operator, i.e., f(Ax) + + A: operator + f: function + + ''' + + def __init__(self, operator, function): + + super(FunctionOperatorComposition, self).__init__() + self.function = function + self.operator = operator + alpha = 1 + + if isinstance (function, ScaledFunction): + alpha = function.scalar + self.L = 2 * alpha * operator.norm()**2 + + + def __call__(self, x): + + ''' Evaluate FunctionOperatorComposition at x + + returns f(Ax) + + ''' + + return self.function(self.operator.direct(x)) + + #TODO do not know if we need it + def call_adjoint(self, x): + + return self.function(self.operator.adjoint(x)) + + + def convex_conjugate(self, x): + + ''' convex_conjugate does not take into account the Operator''' + return self.function.convex_conjugate(x) + + def proximal(self, x, tau, out=None): + + '''proximal does not take into account the Operator''' + if out is None: + return self.function.proximal(x, tau) + else: + self.function.proximal(x, tau, out=out) + + + def proximal_conjugate(self, x, tau, out=None): + + ''' proximal conjugate does not take into account the Operator''' + if out is None: + return self.function.proximal_conjugate(x, tau) + else: + self.function.proximal_conjugate(x, tau, out=out) + + def gradient(self, x, out=None): + + ''' Gradient takes into account the Operator''' + if out is None: + return self.operator.adjoint( + self.function.gradient(self.operator.direct(x)) + ) + else: + self.operator.adjoint( + self.function.gradient(self.operator.direct(x), + out=out) + ) + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 36ba146..cf1a244 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -111,11 +111,11 @@ class KullbackLeibler(Function): self.b.multiply(4*tau, out=out) z_m.multiply(z_m, out=z_m) out += z_m + out.sqrt(out=out) - z_m.sqrt(out=z_m) - z_m += 2 + out *= -1 - out += z_m + out += tmp2 out *= 0.5 @@ -131,23 +131,6 @@ class KullbackLeibler(Function): return ScaledFunction(self, scalar) - - -if __name__ == '__main__': - - - from ccpi.framework import ImageGeometry - import numpy - - N, M = 2,3 - ig = ImageGeometry(N, M) - data = ImageData(numpy.random.randint(-10, 10, size=(M, N))) - x = ImageData(numpy.random.randint(-10, 10, size=(M, N))) - - bnoise = ImageData(numpy.random.randint(-10, 10, size=(M, N))) - - f = KullbackLeibler(data) - print(f(x)) diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 8740376..b77d472 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -19,6 +19,7 @@ from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction +from ccpi.optimisation.functions import FunctionOperatorComposition class L2NormSquared(Function): @@ -97,8 +98,6 @@ class L2NormSquared(Function): tmp /= (1+2*tau) tmp += self.b return tmp -# return (x-self.b)/(1+2*tau) + self.b - else: if self.b is not None: @@ -107,6 +106,7 @@ class L2NormSquared(Function): out += self.b else: x.divide((1+2*tau), out=out) + def proximal_conjugate(self, x, tau, out=None): @@ -139,9 +139,9 @@ class L2NormSquared(Function): return ScaledFunction(self, scalar) - def operator_composition(self, operator): + def composition(self, operator): - return FunctionOperatorComposition + return FunctionOperatorComposition(operator) diff --git a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py index 7caeab2..1db223b 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py @@ -59,7 +59,8 @@ class ScaledFunction(object): if out is None:
return self.scalar * self.function.gradient(x)
else:
- out.fill( self.scalar * self.function.gradient(x) )
+ self.function.gradient(x, out=out)
+ out *= self.scalar
def proximal(self, x, tau, out=None):
'''This returns the proximal operator for the function at x, tau
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index e0b8a32..6ffaf70 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -111,7 +111,7 @@ class Gradient(LinearOperator): return BlockDataContainer(*mat) - def sum_abs_row(self): + def sum_abs_col(self): tmp = self.gm_range.allocate() res = self.gm_domain.allocate() @@ -120,7 +120,7 @@ class Gradient(LinearOperator): res += spMat.sum_abs_row() return res - def sum_abs_col(self): + def sum_abs_row(self): tmp = self.gm_range.allocate() res = [] diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py index a58a296..a853b8d 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py @@ -50,7 +50,7 @@ class Identity(LinearOperator): def sum_abs_row(self): - return self.gm_domain.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))) + return self.gm_range.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))) def sum_abs_col(self): diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py index 5e318ff..c5c2ec9 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py @@ -65,13 +65,13 @@ class SparseFiniteDiff(): def sum_abs_row(self): res = np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')) - res[res==0]=1 + #res[res==0]=0 return ImageData(res) def sum_abs_col(self): res = np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F') ) - res[res==0]=1 + #res[res==0]=0 return ImageData(res) if __name__ == '__main__': diff --git a/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py b/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py new file mode 100644 index 0000000..eb62761 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py @@ -0,0 +1,117 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + + +""" +Compare FISTA & PDHG classes + + +Problem: min_x alpha * ||\grad x ||^{2}_{2} + || x - g ||_{1} + + A: Projection operator + g: Sinogram + +""" + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import FISTA, PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity +from ccpi.optimisation.functions import L2NormSquared, L1Norm, \ + FunctionOperatorComposition, BlockFunction, ZeroFunction + +from skimage.util import random_noise + + +# Create Ground truth and noisy data + +N = 100 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 5 + + +# Setup and run the FISTA algorithm +operator = Gradient(ig) +fidelity = L1Norm(b=noisy_data) +regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator) + +x_init = ig.allocate() +opt = {'memopt':True} +fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt) +fista.max_iteration = 2000 +fista.update_objective_interval = 50 +fista.run(2000, verbose=False) + +# Setup and run the PDHG algorithm +op1 = Gradient(ig) +op2 = Identity(ig, ag) + +operator = BlockOperator(op1, op2, shape=(2,1) ) +f = BlockFunction(alpha * L2NormSquared(), fidelity) +g = ZeroFunction() + +normK = operator.norm() + +sigma = 1 +tau = 1/(sigma*normK**2) + +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000, verbose=False) + +#%% +# Show results + +plt.figure(figsize=(15,15)) + +plt.subplot(1,2,1) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') + +plt.subplot(1,2,2) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA reconstruction') + +plt.show() + +diff1 = pdhg.get_output() - fista.get_output() + +plt.imshow(diff1.abs().as_array()) +plt.title('Diff PDHG vs CGLS') +plt.colorbar() +plt.show() + + diff --git a/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py b/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py new file mode 100644 index 0000000..c877018 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py @@ -0,0 +1,198 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + + +""" +Compare Least Squares minimization problem using FISTA, PDHG, CGLS classes +and Astra Built-in CGLS + +Problem: min_x || A x - g ||_{2}^{2} + + A: Projection operator + g: Sinogram + +""" + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA + +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition +from ccpi.astra.ops import AstraProjectorSimple +import astra + +# Create Ground truth phantom and Sinogram + +N = 128 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D', angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = sin + +# Setup and run Astra CGLS algorithm + +vol_geom = astra.create_vol_geom(N, N) +proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) + +proj_id = astra.create_projector('strip', proj_geom, vol_geom) + +# Create a sinogram from a phantom +sinogram_id, sinogram = astra.create_sino(x, proj_id) + +# Create a data object for the reconstruction +rec_id = astra.data2d.create('-vol', vol_geom) + +cgls_astra = astra.astra_dict('CGLS') +cgls_astra['ReconstructionDataId'] = rec_id +cgls_astra['ProjectionDataId'] = sinogram_id +cgls_astra['ProjectorId'] = proj_id + +# Create the algorithm object from the configuration structure +alg_id = astra.algorithm.create(cgls_astra) + +astra.algorithm.run(alg_id, 1000) + +recon_cgls_astra = astra.data2d.get(rec_id) + +# Setup and run the CGLS algorithm + +x_init = ig.allocate() + +cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) +cgls.max_iteration = 1000 +cgls.update_objective_interval = 200 +cgls.run(1000) + +# Setup and run the PDHG algorithm + +operator = Aop +f = L2NormSquared(b = noisy_data) +g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 1000 +pdhg.update_objective_interval = 200 +pdhg.run(1000) + +# Setup and run the FISTA algorithm + +fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop) +regularizer = ZeroFunction() + +opt = {'memopt':True} +fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt) +fista.max_iteration = 1000 +fista.update_objective_interval = 200 +fista.run(1000, verbose=True) + +#%% Show results + +plt.figure(figsize=(15,15)) +plt.suptitle('Reconstructions ') + +plt.subplot(2,2,1) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') + +plt.subplot(2,2,2) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA reconstruction') + +plt.subplot(2,2,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') + +plt.subplot(2,2,4) +plt.imshow(recon_cgls_astra) +plt.title('CGLS astra') + +#%% +diff1 = pdhg.get_output() - cgls.get_output() +diff2 = fista.get_output() - cgls.get_output() +diff3 = ImageData(recon_cgls_astra) - cgls.get_output() + +print( diff1.squared_norm()) +print( diff2.squared_norm()) +print( diff3.squared_norm()) + +plt.figure(figsize=(15,15)) + +plt.subplot(3,1,1) +plt.imshow(diff1.abs().as_array()) +plt.title('Diff PDHG vs CGLS') +plt.colorbar() + +plt.subplot(3,1,2) +plt.imshow(diff2.abs().as_array()) +plt.title('Diff FISTA vs CGLS') +plt.colorbar() + +plt.subplot(3,1,3) +plt.imshow(diff3.abs().as_array()) +plt.title('Diff CLGS astra vs CGLS') +plt.colorbar() + + + + + + + + + + + + + + + + + + + +# +# +# +# +# +# +# +# diff --git a/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py b/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py new file mode 100644 index 0000000..3155654 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py @@ -0,0 +1,127 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS + +from ccpi.optimisation.operators import Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition +from skimage.util import random_noise +from ccpi.astra.ops import AstraProjectorSimple + +#%% + +N = 128 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = sin + +x_init = ig.allocate() + +## Setup and run the CGLS algorithm +cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) +cgls.max_iteration = 500 +cgls.update_objective_interval = 50 +cgls.run(500, verbose=True) + +# Create BlockOperator +operator = Aop +f = 0.5 * L2NormSquared(b = noisy_data) +g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +sigma = 0.1 +tau = 1/(sigma*normK**2) +# +# +## Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +#%% + +diff = pdhg.get_output() - cgls.get_output() +print( diff.norm()) +# +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(diff.abs().as_array()) +plt.title('Difference reconstruction') +plt.colorbar() +plt.show() + + + + + + + + + + + + + + + + + + + + + + +# +# +# +# +# +# +# +# diff --git a/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py b/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py new file mode 100644 index 0000000..984fca4 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py @@ -0,0 +1,152 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +""" +Compare Tikhonov with PDHG, CGLS classes + + +Problem: min_x alpha * ||\grad x ||^{2}_{2} + || A x - g ||_{2}^{2} + + A: Projection operator + g: Sinogram + +""" + + +from ccpi.framework import ImageData, ImageGeometry, \ + AcquisitionGeometry, BlockDataContainer, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS +from ccpi.optimisation.operators import BlockOperator, Gradient + +from ccpi.optimisation.functions import ZeroFunction, BlockFunction, L2NormSquared +from ccpi.astra.ops import AstraProjectorSimple + +# Create Ground truth phantom and Sinogram + +N = 128 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D', angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,3,ig.shape)) + +# Setup and run the CGLS algorithm + +alpha = 50 +Grad = Gradient(ig) + +# Form Tikhonov as a Block CGLS structure + +#%% +op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1)) +block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate()) + +x_init = ig.allocate() + +cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data) +cgls.max_iteration = 1000 +cgls.update_objective_interval = 200 +cgls.run(1000) + +#%% +#Setup and run the PDHG algorithm + +# Create BlockOperator +op_PDHG = BlockOperator(Grad, Aop, shape=(2,1) ) + +# Create functions + +f1 = 0.5 * alpha**2 * L2NormSquared() +f2 = 0.5 * L2NormSquared(b = noisy_data) +f = BlockFunction(f1, f2) +g = ZeroFunction() + +## Compute operator Norm +normK = op_PDHG.norm() + +## Primal & dual stepsizes +sigma = 10 +tau = 1/(sigma*normK**2) + +pdhg = PDHG(f=f,g=g,operator=op_PDHG, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 1000 +pdhg.update_objective_interval = 200 +pdhg.run(1000) + +#%% +# Show results + +plt.figure(figsize=(15,15)) + +plt.subplot(1,2,1) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') + +plt.subplot(1,2,2) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') + +plt.show() + +diff1 = pdhg.get_output() - cgls.get_output() + +plt.imshow(diff1.abs().as_array()) +plt.title('Diff PDHG vs CGLS') +plt.colorbar() +plt.show() + +#%% + +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') +plt.plot(np.linspace(0,N,N), cgls.get_output().as_array()[int(N/2),:], label = 'CGLS') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + + + + + + + + +# +# +# +# +# +# +# +# diff --git a/Wrappers/Python/wip/Demos/FISTA_vs_CGLS.py b/Wrappers/Python/wip/Demos/FISTA_vs_CGLS.py new file mode 100644 index 0000000..2dcaa89 --- /dev/null +++ b/Wrappers/Python/wip/Demos/FISTA_vs_CGLS.py @@ -0,0 +1,119 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import FISTA, CGLS + +from ccpi.optimisation.operators import Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition +from skimage.util import random_noise +from ccpi.astra.ops import AstraProjectorSimple + +#%% + +N = 75 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = sin + +fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop) +regularizer = ZeroFunction() + +x_init = ig.allocate() + +## Setup and run the FISTA algorithm +opt = {'tol': 1e-4, 'memopt':True} +fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt) +fista.max_iteration = 500 +fista.update_objective_interval = 50 +fista.run(500, verbose=True) + +## Setup and run the CGLS algorithm +cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) +cgls.max_iteration = 500 +cgls.update_objective_interval = 50 +cgls.run(500, verbose=True) + +diff = fista.get_output() - cgls.get_output() + + +#%% +print( diff.norm()) + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA reconstruction') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(diff.abs().as_array()) +plt.title('Difference reconstruction') +plt.colorbar() +plt.show() + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Wrappers/Python/wip/Demos/FISTA_vs_PDHG.py b/Wrappers/Python/wip/Demos/FISTA_vs_PDHG.py new file mode 100644 index 0000000..b7777ef --- /dev/null +++ b/Wrappers/Python/wip/Demos/FISTA_vs_PDHG.py @@ -0,0 +1,120 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import FISTA, PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity +from ccpi.optimisation.functions import L2NormSquared, L1Norm, \ + MixedL21Norm, FunctionOperatorComposition, BlockFunction, ZeroFunction + +from skimage.util import random_noise + +# Create phantom for TV Gaussian denoising +N = 100 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Add Gaussian noise +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 5 + +operator = Gradient(ig) + +fidelity = L1Norm(b=noisy_data) +regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator) + +x_init = ig.allocate() + +## Setup and run the PDHG algorithm +opt = {'tol': 1e-4, 'memopt':True} +fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt) +fista.max_iteration = 2000 +fista.update_objective_interval = 50 +fista.run(2000, verbose=True) + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(fista.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() + +# Compare with PDHG +# Create operators +op1 = Gradient(ig) +op2 = Identity(ig, ag) + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) +f = BlockFunction(alpha * L2NormSquared(), fidelity) +g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() +# +## Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) +# +# +## Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) +# +#%% +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(np.abs(pdhg.get_output().as_array()-fista.get_output().as_array())) +plt.title('Diff FISTA-PDHG') +plt.colorbar() +plt.show() + + diff --git a/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py new file mode 100644 index 0000000..e67bdb1 --- /dev/null +++ b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py @@ -0,0 +1,164 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData +from astropy.io import fits +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, L2NormSquared,\ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple + + +# load IMAT file +#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_141.fits' +filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_564.fits' + +sino_handler = fits.open(filename_sino) +sino_tmp = numpy.array(sino_handler[0].data, dtype=float) +# reorder sino coordinate: channels, angles, detectors +sinogram = numpy.rollaxis(sino_tmp, 2) +sino_handler.close() +#%% +# white beam data +sinogram_wb = sinogram.sum(axis=0) + +pixh = sinogram_wb.shape[1] # detectors +pixv = sinogram_wb.shape[1] # detectors + +# WhiteBeam Geometry +igWB = ImageGeometry(voxel_num_x = pixh, voxel_num_y = pixv) + +# Load Golden angles +with open("golden_angles.txt") as f: + angles_string = [line.rstrip() for line in f] + angles = numpy.array(angles_string).astype(float) +agWB = AcquisitionGeometry('parallel', '2D', angles * numpy.pi / 180, pixh) +op_WB = AstraProjectorSimple(igWB, agWB, 'gpu') +sinogram_aqdata = AcquisitionData(sinogram_wb, agWB) + +# BackProjection +result_bp = op_WB.adjoint(sinogram_aqdata) + +plt.imshow(result_bp.subset(channel=50).array) +plt.title('BackProjection') +plt.show() + + + +#%% + +# Regularisation Parameter +alpha = 2000 + +# Create operators +op1 = Gradient(igWB) +op2 = op_WB + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +# Create functions + +f1 = alpha * MixedL21Norm() +f2 = KullbackLeibler(sinogram_aqdata) +#f2 = L2NormSquared(b = sinogram_aqdata) +f = BlockFunction(f1, f2) + +g = ZeroFunction() + +diag_precon = False + +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_row() + sigma = 1/ operator.sum_abs_col() + + return tau, sigma + + tau, sigma = tau_sigma_precond(operator) + +else: + # Compute operator Norm + normK = operator.norm() + print ("normK", normK) + # Primal & dual stepsizes + sigma = 0.1 + tau = 1/(sigma*normK**2) + +#%% + + +## Primal & dual stepsizes +#sigma = 0.1 +#tau = 1/(sigma*normK**2) +# +# +## Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 10000 +pdhg.update_objective_interval = 500 + +def circ_mask(h, w, center=None, radius=None): + + if center is None: # use the middle of the image + center = [int(w/2), int(h/2)] + if radius is None: # use the smallest distance between the center and image walls + radius = min(center[0], center[1], w-center[0], h-center[1]) + + Y, X = numpy.ogrid[:h, :w] + dist_from_center = numpy.sqrt((X - center[0])**2 + (Y-center[1])**2) + + mask = dist_from_center <= radius + return mask + +def show_result(niter, objective, solution): + + mask = circ_mask(pixh, pixv, center=None, radius = 220) # 55 with 141, + plt.imshow(solution.as_array() * mask) + plt.colorbar() + plt.title("Iter: {}".format(niter)) + plt.show() + + + print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\ + format(niter, pdhg.max_iteration,'', \ + objective[0],'',\ + objective[1],'',\ + objective[2])) + +pdhg.run(10000, callback = show_result) + +#%% + +mask = circ_mask(pixh, pixv, center=None, radius = 210) # 55 with 141, +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(pdhg.get_output().as_array() * mask) +plt.title('Ground Truth') +plt.colorbar() +plt.show() diff --git a/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt new file mode 100644 index 0000000..95ce73a --- /dev/null +++ b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt @@ -0,0 +1,186 @@ +0 +0.9045 +1.809 +2.368 +3.2725 +4.736 +5.6405 +6.1995 +7.104 +8.5675 +9.472 +10.9356 +11.8401 +12.3991 +13.3036 +14.7671 +15.6716 +16.2306 +17.1351 +18.0396 +18.5986 +19.5031 +20.9666 +21.8711 +22.4301 +23.3346 +24.7981 +25.7026 +27.1661 +28.0706 +28.6297 +29.5342 +30.9977 +31.9022 +32.4612 +33.3657 +34.8292 +35.7337 +37.1972 +38.1017 +38.6607 +39.5652 +41.0287 +41.9332 +42.4922 +43.3967 +44.3012 +44.8602 +45.7647 +47.2283 +48.1328 +48.6918 +49.5963 +51.0598 +51.9643 +53.4278 +54.3323 +54.8913 +55.7958 +57.2593 +58.1638 +58.7228 +59.6273 +60.5318 +61.0908 +61.9953 +63.4588 +64.3633 +64.9224 +65.8269 +67.2904 +68.1949 +69.6584 +70.5629 +71.1219 +72.0264 +73.4899 +74.3944 +74.9534 +75.8579 +77.3214 +78.2259 +79.6894 +80.5939 +81.1529 +82.0574 +83.521 +84.4255 +84.9845 +85.889 +86.7935 +87.3525 +88.257 +89.7205 +90.625 +91.184 +92.0885 +93.552 +94.4565 +95.92 +96.8245 +97.3835 +98.288 +99.7516 +100.656 +101.215 +102.12 +103.583 +104.488 +105.951 +106.856 +107.415 +108.319 +109.783 +110.687 +111.246 +112.151 +113.055 +113.614 +114.519 +115.982 +116.887 +117.446 +118.35 +119.814 +120.718 +122.182 +123.086 +123.645 +124.55 +126.013 +126.918 +127.477 +128.381 +129.286 +129.845 +130.749 +132.213 +133.117 +133.676 +134.581 +136.044 +136.949 +138.412 +139.317 +139.876 +140.78 +142.244 +143.148 +143.707 +144.612 +146.075 +146.98 +148.443 +149.348 +149.907 +150.811 +152.275 +153.179 +153.738 +154.643 +155.547 +156.106 +157.011 +158.474 +159.379 +159.938 +160.842 +162.306 +163.21 +164.674 +165.578 +166.137 +167.042 +168.505 +169.41 +169.969 +170.873 +172.337 +173.242 +174.705 +175.609 +176.168 +177.073 +178.536 +179.441 diff --git a/Wrappers/Python/wip/Demos/LeastSq_CGLS_FISTA_PDHG.py b/Wrappers/Python/wip/Demos/LeastSq_CGLS_FISTA_PDHG.py new file mode 100644 index 0000000..97c71ba --- /dev/null +++ b/Wrappers/Python/wip/Demos/LeastSq_CGLS_FISTA_PDHG.py @@ -0,0 +1,154 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA + +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition +from ccpi.astra.ops import AstraProjectorSimple + +#%% + +N = 68 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = sin + + +#%% +############################################################################### +## Setup and run the CGLS algorithm + +x_init = ig.allocate() +cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) +cgls.max_iteration = 500 +cgls.update_objective_interval = 50 +cgls.run(500, verbose=True) + +#%% +plt.imshow(cgls.get_output().as_array()) +#%% +############################################################################### +## Setup and run the PDHG algorithm + +operator = Aop +f = L2NormSquared(b = noisy_data) +g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +sigma = 0.1 +tau = 1/(sigma*normK**2) + +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + + +#%% +############################################################################### +## Setup and run the FISTA algorithm + +fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop) +regularizer = ZeroFunction() + +## Setup and run the FISTA algorithm +opt = {'memopt':True} +fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt) +fista.max_iteration = 2000 +fista.update_objective_interval = 200 +fista.run(2000, verbose=True) + +#%% Show results + +diff1 = pdhg.get_output() - cgls.get_output() +diff2 = fista.get_output() - cgls.get_output() + +print( diff1.norm()) +print( diff2.norm()) + +plt.figure(figsize=(10,10)) +plt.subplot(2,3,1) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') +plt.subplot(2,3,2) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') +plt.subplot(2,3,3) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA reconstruction') +plt.subplot(2,3,4) +plt.imshow(diff1.abs().as_array()) +plt.title('Diff PDHG vs CGLS') +plt.colorbar() +plt.subplot(2,3,5) +plt.imshow(diff2.abs().as_array()) +plt.title('Diff FISTA vs CGLS') +plt.colorbar() +plt.show() + + + + + + + + + + + + + + + + + + + + + + +# +# +# +# +# +# +# +# diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py index 26578bb..49d4db6 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py +++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py @@ -51,7 +51,7 @@ detectors = N angles = np.linspace(0, np.pi, N, dtype=np.float32) ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') +Aop = AstraProjectorSimple(ig, ag, 'gpu') sin = Aop.direct(data) # Create noisy data. Apply Poisson noise diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py index 1a3e0df..860e76e 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py @@ -17,6 +17,28 @@ # See the License for the specific language governing permissions and # limitations under the License. + +""" + +Total Variation Denoising using PDHG algorithm: + + min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) + + +Problem: min_x \alpha * ||\nabla x||_{1} + || x - g ||_{2}^{2} + + \nabla: Gradient operator + g: Noisy Data with Gaussian Noise + \alpha: Regularization parameter + + Method = 0: K = [ \nabla, + Identity] + + Method = 1: K = \nabla + + +""" + from ccpi.framework import ImageData, ImageGeometry import numpy as np @@ -29,10 +51,9 @@ from ccpi.optimisation.operators import BlockOperator, Identity, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction -from skimage.util import random_noise # Create phantom for TV Gaussian denoising -N = 300 +N = 200 data = np.zeros((N,N)) data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 @@ -42,11 +63,11 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) ag = ig # Create noisy data. Add Gaussian noise -n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10) -noisy_data = ImageData(n1) +np.random.seed(10) +noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.05, size=ig.shape) ) # Regularisation Parameter -alpha = 0.5 +alpha = 2 method = '0' @@ -70,11 +91,10 @@ if method == '0': else: # Without the "Block Framework" - operator = Gradient(ig) + operator = Gradient(ig) f = alpha * MixedL21Norm() g = 0.5 * L2NormSquared(b = noisy_data) - # Compute operator Norm normK = operator.norm() @@ -82,13 +102,11 @@ normK = operator.norm() sigma = 1 tau = 1/(sigma*normK**2) - # Setup and run the PDHG algorithm pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(2000) - +pdhg.max_iteration = 3000 +pdhg.update_objective_interval = 200 +pdhg.run(3000, verbose=False) plt.figure(figsize=(15,15)) plt.subplot(3,1,1) @@ -104,7 +122,7 @@ plt.imshow(pdhg.get_output().as_array()) plt.title('TV Reconstruction') plt.colorbar() plt.show() -## + plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') plt.legend() @@ -112,7 +130,7 @@ plt.title('Middle Line Profiles') plt.show() -##%% Check with CVX solution +#%% Check with CVX solution from ccpi.optimisation.operators import SparseFiniteDiff @@ -143,7 +161,7 @@ if cvx_not_installable: obj = Minimize( regulariser + fidelity) prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) + result = prob.solve(verbose = True, solver = MOSEK) diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) @@ -164,6 +182,8 @@ if cvx_not_installable: plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'Truth') + plt.legend() plt.title('Middle Line Profiles') plt.show() @@ -171,7 +191,3 @@ if cvx_not_installable: print('Primal Objective (CVX) {} '.format(obj.value)) print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) - - - - diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_DiagPrecond.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_DiagPrecond.py new file mode 100644 index 0000000..d65478c --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_DiagPrecond.py @@ -0,0 +1,208 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + + +""" +Total Variation Denoising using PDHG algorithm + +Problem: min_x \alpha * ||\nabla x||_{1} + || x - g ||_{2}^{2} + + \nabla: Gradient operator + g: Noisy Data with Gaussian Noise + \alpha: Regularization parameter + +""" + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + + +# Create phantom for TV Gaussian denoising +N = 400 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Add Gaussian noise +np.random.seed(10) +noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.05, size=ig.shape) ) + +# Regularisation Parameter +alpha = 2 + +method = '1' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * MixedL21Norm() + f2 = 0.5 * L2NormSquared(b = noisy_data) + f = BlockFunction(f1, f2) + + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = 0.5 * L2NormSquared(b = noisy_data) + + +diag_precon = False + + +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_col() + sigma = 1/operator.sum_abs_row() + + sigma[0].as_array()[sigma[0].as_array()==np.inf]=0 + sigma[1].as_array()[sigma[1].as_array()==np.inf]=0 + + return tau, sigma + + tau, sigma = tau_sigma_precond(operator) + +else: + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 3000 +pdhg.update_objective_interval = 200 +pdhg.run(3000, verbose=False) + +#%% +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + +#%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(ig.shape) + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = MOSEK) + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output().as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(u.value) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'Truth') + + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py index ccdabb2..3c295f5 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py @@ -17,6 +17,27 @@ # See the License for the specific language governing permissions and # limitations under the License. +""" + +Total Variation Denoising using PDHG algorithm: + + min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + \int x - g * log(x) + + \nabla: Gradient operator + g: Noisy Data with Poisson Noise + \alpha: Regularization parameter + + Method = 0: K = [ \nabla, + Identity] + + Method = 1: K = \nabla + + +""" + from ccpi.framework import ImageData, ImageGeometry import numpy as np @@ -124,63 +145,63 @@ plt.show() #%% Check with CVX solution -#from ccpi.optimisation.operators import SparseFiniteDiff -# -#try: -# from cvxpy import * -# cvx_not_installable = True -#except ImportError: -# cvx_not_installable = False -# -# -#if cvx_not_installable: -# -# ##Construct problem -# u1 = Variable(ig.shape) -# q = Variable() -# -# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') -# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') -# -# # Define Total Variation as a regulariser -# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0)) -# -# fidelity = sum( u1 - multiply(noisy_data.as_array(), log(u1)) ) -# constraints = [q>= fidelity, u1>=0] -# -# solver = ECOS -# obj = Minimize( regulariser + q) -# prob = Problem(obj, constraints) -# result = prob.solve(verbose = True, solver = solver) -# -# -# diff_cvx = numpy.abs( pdhg.get_output().as_array() - u1.value ) -# -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(pdhg.get_output().as_array()) -# plt.title('PDHG solution') -# plt.colorbar() -# plt.subplot(3,1,2) -# plt.imshow(u1.value) -# plt.title('CVX solution') -# plt.colorbar() -# plt.subplot(3,1,3) -# plt.imshow(diff_cvx) -# plt.title('Difference') -# plt.colorbar() -# plt.show() -# -# plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX') -# plt.legend() -# plt.title('Middle Line Profiles') -# plt.show() -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) -# -# -# -# -# +from ccpi.optimisation.operators import SparseFiniteDiff + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u1 = Variable(ig.shape) + q = Variable() + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0)) + + fidelity = sum( u1 - multiply(noisy_data.as_array(), log(u1)) ) + constraints = [q>= fidelity, u1>=0] + + solver = ECOS + obj = Minimize( regulariser + q) + prob = Problem(obj, constraints) + result = prob.solve(verbose = True, solver = solver) + + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - u1.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output().as_array()) + plt.title('PDHG solution') + plt.colorbar() + plt.subplot(3,1,2) + plt.imshow(u1.value) + plt.title('CVX solution') + plt.colorbar() + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') + plt.colorbar() + plt.show() + + plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') + plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX') + plt.legend() + plt.title('Middle Line Profiles') + plt.show() + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py index 484fe42..f5d4ce4 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py @@ -17,6 +17,27 @@ # See the License for the specific language governing permissions and # limitations under the License. +""" + +Total Variation Denoising using PDHG algorithm: + + min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + ||x-g||_{1} + + \nabla: Gradient operator + g: Noisy Data with Salt & Pepper Noise + \alpha: Regularization parameter + + Method = 0: K = [ \nabla, + Identity] + + Method = 1: K = \nabla + + +""" + from ccpi.framework import ImageData, ImageGeometry import numpy as np @@ -48,7 +69,7 @@ noisy_data = ImageData(n1) # Regularisation Parameter alpha = 2 -method = '1' +method = '0' if method == '0': diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py index 0711e91..87d5328 100644 --- a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py @@ -31,6 +31,25 @@ from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ from ccpi.astra.ops import AstraProjectorSimple +""" + +Total Variation Denoising using PDHG algorithm: + + min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + int A x -g log(Ax + \eta) + + \nabla: Gradient operator + + A: Projection Matrix + g: Noisy sinogram corrupted with Poisson Noise + + \eta: Background Noise + \alpha: Regularization parameter + +""" + # Create phantom for TV 2D tomography N = 75 x = np.zeros((N,N)) @@ -70,12 +89,25 @@ f = BlockFunction(f1, f2) g = ZeroFunction() -# Compute operator Norm -normK = operator.norm() +diag_precon = True -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_row() + sigma = 1/ operator.sum_abs_col() + + return tau, sigma + + tau, sigma = tau_sigma_precond(operator) + +else: + # Compute operator Norm + normK = operator.norm() + # Primal & dual stepsizes + sigma = 10 + tau = 1/(sigma*normK**2) # Setup and run the PDHG algorithm @@ -84,6 +116,8 @@ pdhg.max_iteration = 2000 pdhg.update_objective_interval = 50 pdhg.run(2000) + +#%% plt.figure(figsize=(15,15)) plt.subplot(3,1,1) plt.imshow(data.as_array()) diff --git a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py index 3f275e2..041d4ee 100644 --- a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py +++ b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py @@ -47,7 +47,7 @@ noisy_data = ImageData(n1) # Regularisation Parameter alpha = 4 -method = '1' +method = '0' if method == '0': diff --git a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py index 5c03362..f17c4fe 100644 --- a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py +++ b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py @@ -43,7 +43,7 @@ detectors = N angles = np.linspace(0, np.pi, N, dtype=np.float32) ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') +Aop = AstraProjectorSimple(ig, ag, 'gpu') sin = Aop.direct(data) # Create noisy data. Apply Gaussian noise diff --git a/Wrappers/Python/wip/Demos/PDHG_vs_CGLS.py b/Wrappers/Python/wip/Demos/PDHG_vs_CGLS.py new file mode 100644 index 0000000..3155654 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_vs_CGLS.py @@ -0,0 +1,127 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS + +from ccpi.optimisation.operators import Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition +from skimage.util import random_noise +from ccpi.astra.ops import AstraProjectorSimple + +#%% + +N = 128 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = sin + +x_init = ig.allocate() + +## Setup and run the CGLS algorithm +cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) +cgls.max_iteration = 500 +cgls.update_objective_interval = 50 +cgls.run(500, verbose=True) + +# Create BlockOperator +operator = Aop +f = 0.5 * L2NormSquared(b = noisy_data) +g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +sigma = 0.1 +tau = 1/(sigma*normK**2) +# +# +## Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +#%% + +diff = pdhg.get_output() - cgls.get_output() +print( diff.norm()) +# +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(diff.abs().as_array()) +plt.title('Difference reconstruction') +plt.colorbar() +plt.show() + + + + + + + + + + + + + + + + + + + + + + +# +# +# +# +# +# +# +# diff --git a/Wrappers/Python/wip/Demos/check_blockOperator_sum_row_cols.py b/Wrappers/Python/wip/Demos/check_blockOperator_sum_row_cols.py new file mode 100644 index 0000000..bdb2c38 --- /dev/null +++ b/Wrappers/Python/wip/Demos/check_blockOperator_sum_row_cols.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri May 3 13:10:09 2019 + +@author: evangelos +""" + +from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff, BlockOperator, Gradient +from ccpi.framework import ImageGeometry, AcquisitionGeometry, BlockDataContainer, ImageData +from ccpi.astra.ops import AstraProjectorSimple + +from scipy import sparse +import numpy as np + +N = 3 +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +u = ig.allocate('random_int') + +# Compare FiniteDiff with SparseFiniteDiff + +DY = FiniteDiff(ig, direction = 0, bnd_cond = 'Neumann') +DX = FiniteDiff(ig, direction = 1, bnd_cond = 'Neumann') + +DXu = DX.direct(u) +DYu = DY.direct(u) + +DX_sparse = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +DY_sparse = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + +DXu_sparse = DX_sparse.direct(u) +DYu_sparse = DY_sparse.direct(u) + +#np.testing.assert_array_almost_equal(DYu.as_array(), DYu_sparse.as_array(), decimal=4) +#np.testing.assert_array_almost_equal(DXu.as_array(), DXu_sparse.as_array(), decimal=4) + +#%% Tau/ Sigma + +A1 = DY_sparse.matrix() +A2 = DX_sparse.matrix() +A3 = sparse.eye(np.prod(ig.shape)) + +sum_rows1 = np.array(np.sum(abs(A1), axis=1)) +sum_rows2 = np.array(np.sum(abs(A2), axis=1)) +sum_rows3 = np.array(np.sum(abs(A3), axis=1)) + +sum_cols1 = np.array(np.sum(abs(A1), axis=0)) +sum_cols2 = np.array(np.sum(abs(A2), axis=0)) +sum_cols3 = np.array(np.sum(abs(A2), axis=0)) + +# Check if Grad sum row/cols is OK +Grad = Gradient(ig) + +Sum_Block_row = Grad.sum_abs_row() +Sum_Block_col = Grad.sum_abs_col() + +tmp1 = BlockDataContainer( ImageData(np.reshape(sum_rows1, ig.shape, order='F')),\ + ImageData(np.reshape(sum_rows2, ig.shape, order='F'))) + + +#np.testing.assert_array_almost_equal(tmp1[0].as_array(), Sum_Block_row[0].as_array(), decimal=4) +#np.testing.assert_array_almost_equal(tmp1[1].as_array(), Sum_Block_row[1].as_array(), decimal=4) + +tmp2 = ImageData(np.reshape(sum_cols1 + sum_cols2, ig.shape, order='F')) + +#np.testing.assert_array_almost_equal(tmp2.as_array(), Sum_Block_col.as_array(), decimal=4) + + +#%% BlockOperator with Gradient, Identity + +Id = Identity(ig) +Block_GrId = BlockOperator(Grad, Id, shape=(2,1)) + + +Sum_Block_GrId_row = Block_GrId.sum_abs_row() + + + + + + + + + + + + + + diff --git a/Wrappers/Python/wip/Demos/check_precond.py b/Wrappers/Python/wip/Demos/check_precond.py new file mode 100644 index 0000000..8cf95fa --- /dev/null +++ b/Wrappers/Python/wip/Demos/check_precond.py @@ -0,0 +1,182 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple + +# Create phantom for TV 2D tomography +N = 75 +x = np.zeros((N,N)) +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 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +# Create noisy data +np.random.seed(10) +n1 = np.random.random(sin.shape) +noisy_data = sin + ImageData(5*n1) + +#%% + +# Regularisation Parameter +alpha = 50 + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + + + +# Create functions + +f1 = alpha * MixedL21Norm() +f2 = L2NormSquared(b=noisy_data) +f = BlockFunction(f1, f2) + +g = ZeroFunction() + +diag_precon = True + +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_row() + sigma = 1/ operator.sum_abs_col() + + return tau, sigma + + tau, sigma = tau_sigma_precond(operator) + +else: + # Compute operator Norm + normK = operator.norm() + # Primal & dual stepsizes + sigma = 10 + tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 1000 +pdhg.update_objective_interval = 200 +pdhg.run(1000) + +#%% Check with CVX solution + +from ccpi.optimisation.operators import SparseFiniteDiff +import astra +import numpy + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(N*N) + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + + # create matrix representation for Astra operator + + vol_geom = astra.create_vol_geom(N, N) + proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) + + proj_id = astra.create_projector('line', proj_geom, vol_geom) + + matrix_id = astra.projector.matrix(proj_id) + + ProjMat = astra.matrix.get(matrix_id) + + fidelity = sum_squares( ProjMat * u - noisy_data.as_array().ravel()) + #constraints = [q>=fidelity] +# constraints = [u>=0] + + solver = MOSEK + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + +#%% + +plt.figure(figsize=(15,15)) +plt.subplot(2,2,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') + +plt.subplot(2,2,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') + +plt.subplot(2,2,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG Reconstruction') + +plt.subplot(2,2,4) +plt.imshow(np.reshape(u.value, ig.shape)) +plt.title('CVX Reconstruction') + +plt.show() + +#%% +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') +plt.plot(np.linspace(0,N,N), np.reshape(u.value, ig.shape)[int(N/2),:], label = 'CVX') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + + + + + + + diff --git a/Wrappers/Python/wip/Demos/fista_test.py b/Wrappers/Python/wip/Demos/fista_test.py new file mode 100644 index 0000000..dd1f6fa --- /dev/null +++ b/Wrappers/Python/wip/Demos/fista_test.py @@ -0,0 +1,127 @@ +# -*- 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-2019 Evangelos Papoutsellis 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. + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import FISTA, PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity +from ccpi.optimisation.functions import L2NormSquared, L1Norm, \ + MixedL21Norm, FunctionOperatorComposition, BlockFunction, ZeroFunction + +from skimage.util import random_noise + +# Create phantom for TV Gaussian denoising +N = 100 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Add Gaussian noise +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 5 + +operator = Gradient(ig) + +#fidelity = L1Norm(b=noisy_data) +#regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator) + +fidelity = FunctionOperatorComposition(alpha * MixedL21Norm(), operator) +regulariser = 0.5 * L2NormSquared(b = noisy_data) + +x_init = ig.allocate() + +## Setup and run the PDHG algorithm +opt = {'tol': 1e-4, 'memopt':True} +fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt) +fista.max_iteration = 2000 +fista.update_objective_interval = 50 +fista.run(2000, verbose=True) + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(fista.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() + +# Compare with PDHG +method = '0' +# +if method == '0': +# +# # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) +# +# # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + f = BlockFunction(alpha * L2NormSquared(), fidelity) + g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() +# +## Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) +# +# +## Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) +# +#%% +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(np.abs(pdhg.get_output().as_array()-fista.get_output().as_array())) +plt.title('Diff FISTA-PDHG') +plt.colorbar() +plt.show() + + diff --git a/Wrappers/Python/wip/old_demos/demo_colourbay.py b/Wrappers/Python/wip/old_demos/demo_colourbay.py new file mode 100644 index 0000000..5dbf2e1 --- /dev/null +++ b/Wrappers/Python/wip/old_demos/demo_colourbay.py @@ -0,0 +1,137 @@ +# This script demonstrates how to load a mat-file with UoM colour-bay data +# into the CIL optimisation framework and run (simple) multichannel +# reconstruction methods. + +# All third-party imports. +import numpy +from scipy.io import loadmat +import matplotlib.pyplot as plt + +# All own imports. +from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.astra.ops import AstraProjectorMC +from ccpi.optimisation.algs import CGLS, FISTA +from ccpi.optimisation.funcs import Norm2sq, Norm1 + +# Load full data and permute to expected ordering. Change path as necessary. +# The loaded X has dims 80x60x80x150, which is pix x angle x pix x channel. +# Permute (numpy.transpose) puts into our default ordering which is +# (channel, angle, vertical, horizontal). + +pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/ColourBay/spectral_data_sets/CarbonPd/' +filename = 'carbonPd_full_sinogram_stripes_removed.mat' + +X = loadmat(pathname + filename) +X = numpy.transpose(X['SS'],(3,1,2,0)) + +# Store geometric variables for reuse +num_channels = X.shape[0] +num_pixels_h = X.shape[3] +num_pixels_v = X.shape[2] +num_angles = X.shape[1] + +# Display a single projection in a single channel +plt.imshow(X[100,5,:,:]) +plt.title('Example of a projection image in one channel' ) +plt.show() + +# Set angles to use +angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False) + +# Define full 3D acquisition geometry and data container. +# Geometric info is taken from the txt-file in the same dir as the mat-file +ag = AcquisitionGeometry('cone', + '3D', + angles, + pixel_num_h=num_pixels_h, + pixel_size_h=0.25, + pixel_num_v=num_pixels_v, + pixel_size_v=0.25, + dist_source_center=233.0, + dist_center_detector=245.0, + channels=num_channels) +data = AcquisitionData(X, geometry=ag) + +# Reduce to central slice by extracting relevant parameters from data and its +# geometry. Perhaps create function to extract central slice automatically? +data2d = data.subset(vertical=40) +ag2d = AcquisitionGeometry('cone', + '2D', + ag.angles, + pixel_num_h=ag.pixel_num_h, + pixel_size_h=ag.pixel_size_h, + pixel_num_v=1, + pixel_size_v=ag.pixel_size_h, + dist_source_center=ag.dist_source_center, + dist_center_detector=ag.dist_center_detector, + channels=ag.channels) +data2d.geometry = ag2d + +# Set up 2D Image Geometry. +# First need the geometric magnification to scale the voxel size relative +# to the detector pixel size. +mag = (ag.dist_source_center + ag.dist_center_detector)/ag.dist_source_center +ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, + voxel_num_y=ag2d.pixel_num_h, + voxel_size_x=ag2d.pixel_size_h/mag, + voxel_size_y=ag2d.pixel_size_h/mag, + channels=X.shape[0]) + +# Create GPU multichannel projector/backprojector operator with ASTRA. +Aall = AstraProjectorMC(ig2d,ag2d,'gpu') + +# Compute and simple backprojction and display one channel as image. +Xbp = Aall.adjoint(data2d) +plt.imshow(Xbp.subset(channel=100).array) +plt.show() + +# Set initial guess ImageData with zeros for algorithms, and algorithm options. +x_init = ImageData(numpy.zeros((num_channels,num_pixels_v,num_pixels_h)), + geometry=ig2d, + dimension_labels=['channel','horizontal_y','horizontal_x']) +opt_CGLS = {'tol': 1e-4, 'iter': 5} + +# Run CGLS algorithm and display one channel. +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aall, data2d, opt_CGLS) + +plt.imshow(x_CGLS.subset(channel=100).array) +plt.title('CGLS') +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS Criterion vs iterations') +plt.show() + +# Create least squares object instance with projector, test data and a constant +# coefficient of 0.5. Note it is least squares over all channels. +f = Norm2sq(Aall,data2d,c=0.5) + +# Options for FISTA algorithm. +opt = {'tol': 1e-4, 'iter': 100} + +# Run FISTA for least squares without regularization and display one channel +# reconstruction as image. +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) + +plt.imshow(x_fista0.subset(channel=100).array) +plt.title('FISTA LS') +plt.show() + +plt.semilogy(criter0) +plt.title('FISTA LS Criterion vs iterations') +plt.show() + +# Set up 1-norm regularisation (over all channels), solve with FISTA, and +# display one channel of reconstruction. +lam = 0.1 +g0 = Norm1(lam) + +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) + +plt.imshow(x_fista1.subset(channel=100).array) +plt.title('FISTA LS+1') +plt.show() + +plt.semilogy(criter1) +plt.title('FISTA LS+1 Criterion vs iterations') +plt.show()
\ No newline at end of file diff --git a/Wrappers/Python/wip/old_demos/demo_compare_cvx.py b/Wrappers/Python/wip/old_demos/demo_compare_cvx.py new file mode 100644 index 0000000..27b1c97 --- /dev/null +++ b/Wrappers/Python/wip/old_demos/demo_compare_cvx.py @@ -0,0 +1,306 @@ + +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, Norm2 + +from ccpi.optimisation.ops import LinearOperatorMatrix, TomoIdentity +from ccpi.optimisation.ops import Identity +from ccpi.optimisation.ops import FiniteDiff2D + +# Requires CVXPY, see http://www.cvxpy.org/ +# CVXPY can be installed in anaconda using +# conda install -c cvxgrp cvxpy libgcc + +# Whether to use or omit CVXPY +use_cvxpy = True +if use_cvxpy: + from cvxpy import * + +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 +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]) + +if use_cvxpy: + # 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) + +# Plot criterion curve to see FISTA converge to same value as CVX. +iternum = np.arange(1,1001) +plt.figure() +plt.loglog(iternum[[0,-1]],[objective0.value, objective0.value], label='CVX LS') +plt.loglog(iternum,criter0,label='FISTA LS') +plt.legend() +plt.show() + +# Create 1-norm object instance +g1 = Norm1(lam) + +g1(x_init) +x_rand = DataContainer(np.reshape(np.random.rand(n),(n,1))) +x_rand2 = DataContainer(np.reshape(np.random.rand(n-1),(n-1,1))) +v = g1.prox(x_rand,0.02) +#vv = g1.prox(x_rand2,0.02) +vv = v.copy() +vv *= 0 +print (">>>>>>>>>>vv" , vv.as_array()) +vv.fill(v) +print (">>>>>>>>>>fill" , vv.as_array()) +g1.proximal(x_rand, 0.02, out=vv) +print (">>>>>>>>>>v" , v.as_array()) +print (">>>>>>>>>>gradient" , vv.as_array()) + +print (">>>>>>>>>>" , (v-vv).as_array()) +import sys +#sys.exit(0) +# 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]) + +if use_cvxpy: + # 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) + +# Now try another algorithm FBPD for same problem: +x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init,Identity(), None, f, g1) +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[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1') +plt.loglog(iternum,criter1,label='FISTA LS+1') +plt.legend() +plt.show() + +plt.figure() +plt.loglog(iternum[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1') +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 = 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 + +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 = y.array + 0.1*np.random.randn(N, N) + +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))) + +# 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]) + +#plt.imshow(x_fista1_denoise.as_array()) +#plt.title('FISTA LS+1') +#plt.show() + +# Now denoise LS + 1-norm with FBPD +x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise, \ + criterfbpd1_denoise = FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) +print(x_fbpd1_denoise) +print(criterfbpd1_denoise[-1]) + +#plt.imshow(x_fbpd1_denoise.as_array()) +#plt.title('FBPD LS+1') +#plt.show() + +if use_cvxpy: + # 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) + +x1_cvx = x1_denoise.value +x1_cvx.shape = (N,N) + + + +#plt.imshow(x1_cvx) +#plt.title('CVX LS+1') +#plt.show() + +fig = plt.figure() +plt.subplot(1,4,1) +plt.imshow(y.array) +plt.title("LS+1") +plt.subplot(1,4,2) +plt.imshow(x_fista1_denoise.as_array()) +plt.title("fista") +plt.subplot(1,4,3) +plt.imshow(x_fbpd1_denoise.as_array()) +plt.title("fbpd") +plt.subplot(1,4,4) +plt.imshow(x1_cvx) +plt.title("cvx") +plt.show() + +############################################################## +# Now TV with FBPD and Norm2 +lam_tv = 0.1 +gtv = TV2D(lam_tv) +norm2 = Norm2(lam_tv) +op = FiniteDiff2D() +#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, op, None, \ + f_denoise, norm2 ,opt=opt_tv) +print(x_fbpdtv_denoise) +print(criterfbpdtv_denoise[-1]) + +plt.imshow(x_fbpdtv_denoise.as_array()) +plt.title('FBPD TV') +#plt.show() + +if use_cvxpy: + # Compare to CVXPY + + # Construct the problem. + xtv_denoise = Variable((N,N)) + #print (xtv_denoise.value.shape) + 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) + +plt.imshow(xtv_denoise.value) +plt.title('CVX TV') +#plt.show() + +fig = plt.figure() +plt.subplot(1,3,1) +plt.imshow(y.array) +plt.title("TV2D") +plt.subplot(1,3,2) +plt.imshow(x_fbpdtv_denoise.as_array()) +plt.title("fbpd tv denoise") +plt.subplot(1,3,3) +plt.imshow(xtv_denoise.value) +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') diff --git a/Wrappers/Python/wip/old_demos/demo_gradient_descent.py b/Wrappers/Python/wip/old_demos/demo_gradient_descent.py new file mode 100755 index 0000000..4d6647e --- /dev/null +++ b/Wrappers/Python/wip/old_demos/demo_gradient_descent.py @@ -0,0 +1,295 @@ + +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, Norm2 + +from ccpi.optimisation.ops import LinearOperatorMatrix, TomoIdentity +from ccpi.optimisation.ops import Identity +from ccpi.optimisation.ops import FiniteDiff2D + +# Requires CVXPY, see http://www.cvxpy.org/ +# CVXPY can be installed in anaconda using +# conda install -c cvxgrp cvxpy libgcc + +# Whether to use or omit CVXPY + +import numpy as np +import matplotlib.pyplot as plt + +class Algorithm(object): + def __init__(self, *args, **kwargs): + pass + def set_up(self, *args, **kwargs): + raise NotImplementedError() + def update(self): + raise NotImplementedError() + + def should_stop(self): + raise NotImplementedError() + + def __iter__(self): + return self + + def __next__(self): + if self.should_stop(): + raise StopIteration() + else: + self.update() + +class GradientDescent(Algorithm): + x = None + rate = 0 + objective_function = None + regulariser = None + iteration = 0 + stop_cryterion = 'max_iter' + __max_iteration = 0 + __loss = [] + def __init__(self, **kwargs): + args = ['x_init', 'objective_function', 'rate'] + present = True + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(x_init=kwargs['x_init'], + objective_function=kwargs['objective_function'], + rate=kwargs['rate']) + + def should_stop(self): + return self.iteration >= self.max_iteration + + def set_up(self, x_init, objective_function, rate): + self.x = x_init.copy() + self.x_update = x_init.copy() + self.objective_function = objective_function + self.rate = rate + self.__loss.append(objective_function(x_init)) + + def update(self): + + self.objective_function.gradient(self.x, out=self.x_update) + self.x_update *= -self.rate + self.x += self.x_update + self.__loss.append(self.objective_function(self.x)) + self.iteration += 1 + + def get_output(self): + return self.x + def get_current_loss(self): + return self.__loss[-1] + @property + def loss(self): + return self.__loss + @property + def max_iteration(self): + return self.__max_iteration + @max_iteration.setter + def max_iteration(self, value): + assert isinstance(value, int) + self.__max_iteration = value + + + + + +# 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]) + +gd = GradientDescent(x_init=x_init, objective_function=f, rate=0.001) +gd.max_iteration = 5000 + +for i,el in enumerate(gd): + if i%100 == 0: + print ("\rIteration {} Loss: {}".format(gd.iteration, + gd.get_current_loss())) + + +#%% + + +# +#if use_cvxpy: +# # 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) +# +## Plot criterion curve to see FISTA converge to same value as CVX. +#iternum = np.arange(1,1001) +#plt.figure() +#plt.loglog(iternum[[0,-1]],[objective0.value, objective0.value], label='CVX LS') +#plt.loglog(iternum,criter0,label='FISTA LS') +#plt.legend() +#plt.show() +# +## Create 1-norm object instance +#g1 = Norm1(lam) +# +#g1(x_init) +#x_rand = DataContainer(np.reshape(np.random.rand(n),(n,1))) +#x_rand2 = DataContainer(np.reshape(np.random.rand(n-1),(n-1,1))) +#v = g1.prox(x_rand,0.02) +##vv = g1.prox(x_rand2,0.02) +#vv = v.copy() +#vv *= 0 +#print (">>>>>>>>>>vv" , vv.as_array()) +#vv.fill(v) +#print (">>>>>>>>>>fill" , vv.as_array()) +#g1.proximal(x_rand, 0.02, out=vv) +#print (">>>>>>>>>>v" , v.as_array()) +#print (">>>>>>>>>>gradient" , vv.as_array()) +# +#print (">>>>>>>>>>" , (v-vv).as_array()) +#import sys +##sys.exit(0) +## 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]) +# +#if use_cvxpy: +# # 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) +# +## Now try another algorithm FBPD for same problem: +#x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init,Identity(), None, f, g1) +#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[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1') +#plt.loglog(iternum,criter1,label='FISTA LS+1') +#plt.legend() +#plt.show() +# +#plt.figure() +#plt.loglog(iternum[[0,-1]],[objective1.value, objective1.value], label='CVX LS+1') +#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 = 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 + +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 = y.array + 0.1*np.random.randn(N, N) + +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))) + +# 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]) + +f_2 = +gd = GradientDescent(x_init=x_init_denoise, + objective_function=f, rate=0.001) +gd.max_iteration = 5000 + +for i,el in enumerate(gd): + if i%100 == 0: + print ("\rIteration {} Loss: {}".format(gd.iteration, + gd.get_current_loss())) + +plt.imshow(gd.get_output().as_array()) +plt.title('GD image') +plt.show() + diff --git a/Wrappers/Python/wip/old_demos/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/old_demos/demo_imat_multichan_RGLTK.py new file mode 100644 index 0000000..8370c78 --- /dev/null +++ b/Wrappers/Python/wip/old_demos/demo_imat_multichan_RGLTK.py @@ -0,0 +1,151 @@ +# This script demonstrates how to load IMAT fits data +# into the CIL optimisation framework and run reconstruction methods. +# +# Demo to reconstruct energy-discretized channels of IMAT data + +# needs dxchange: conda install -c conda-forge dxchange +# needs astropy: conda install astropy + + +# All third-party imports. +import numpy as np +import matplotlib.pyplot as plt +from dxchange.reader import read_fits +from astropy.io import fits + +# All own imports. +from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData, DataContainer +from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple +from ccpi.optimisation.algs import CGLS, FISTA +from ccpi.optimisation.funcs import Norm2sq, Norm1 +from ccpi.plugins.regularisers import FGP_TV + +# set main parameters here +n = 512 +totalAngles = 250 # total number of projection angles +# spectral discretization parameter +num_average = 145 # channel discretization frequency (total number of averaged channels) +numChannels = 2970 # 2970 +totChannels = round(numChannels/num_average) # the resulting number of channels +Projections_stack = np.zeros((num_average,n,n),dtype='uint16') +ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32') + +######################################################################### +print ("Loading the data...") +MainPath = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/' # path to data +pathname0 = '{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/Sample/') +counterFileName = 4675 +# A main loop over all available angles +for ll in range(0,totalAngles,1): + pathnameData = '{!s}{!s}{!s}/'.format(pathname0,'angle',str(ll)) + filenameCurr = '{!s}{!s}{!s}'.format('IMAT0000',str(counterFileName),'_Tomo_test_000_') + counterT = 0 + # loop over reduced channels (discretized) + for i in range(0,totChannels,1): + sumCount = 0 + # loop over actual channels to obtain averaged one + for j in range(0,num_average,1): + if counterT < 10: + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'0000',str(counterT)) + if ((counterT >= 10) & (counterT < 100)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'000',str(counterT)) + if ((counterT >= 100) & (counterT < 1000)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'00',str(counterT)) + if ((counterT >= 1000) & (counterT < 10000)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'0',str(counterT)) + try: + Projections_stack[j,:,:] = read_fits(outfile) + except: + print("Fits is corrupted, skipping no.", counterT) + sumCount -= 1 + counterT += 1 + sumCount += 1 + AverageProj=np.sum(Projections_stack,axis=0)/sumCount # averaged projection over "num_average" channels + ProjAngleChannels[ll,i,:,:] = AverageProj + print("Angle is processed", ll) + counterFileName += 1 +######################################################################### + +flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')) +nonzero = flat1 > 0 +# Apply flat field and take negative log +for ll in range(0,totalAngles,1): + for i in range(0,totChannels,1): + ProjAngleChannels[ll,i,nonzero] = ProjAngleChannels[ll,i,nonzero]/flat1[nonzero] + +eqzero = ProjAngleChannels == 0 +ProjAngleChannels[eqzero] = 1 +ProjAngleChannels_NormLog = -np.log(ProjAngleChannels) # normalised and neg-log data + +# extact sinogram over energy channels +selectedVertical_slice = 256 +sino_all_channels = ProjAngleChannels_NormLog[:,:,:,selectedVertical_slice] +# Set angles to use +angles = np.linspace(-np.pi,np.pi,totalAngles,endpoint=False) + +# set the geometry +ig = ImageGeometry(n,n) +ag = AcquisitionGeometry('parallel', + '2D', + angles, + n,1) +Aop = AstraProjectorSimple(ig, ag, 'gpu') + + +# loop to reconstruct energy channels +REC_chan = np.zeros((totChannels, n, n), 'float32') +for i in range(0,totChannels,1): + sino_channel = sino_all_channels[:,i,:] # extract a sinogram for i-th channel + + print ("Initial guess") + x_init = ImageData(geometry=ig) + + # Create least squares object instance with projector and data. + print ("Create least squares object instance with projector and data.") + f = Norm2sq(Aop,DataContainer(sino_channel),c=0.5) + + print ("Run FISTA-TV for least squares") + lamtv = 5 + opt = {'tol': 1e-4, 'iter': 200} + g_fgp = FGP_TV(lambdaReg = lamtv, + iterationsTV=50, + tolerance=1e-6, + methodTV=0, + nonnegativity=0, + printing=0, + device='gpu') + + x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp, opt) + REC_chan[i,:,:] = x_fista_fgp.array + """ + plt.figure() + plt.subplot(121) + plt.imshow(x_fista_fgp.array, vmin=0, vmax=0.05) + plt.title('FISTA FGP TV') + plt.subplot(122) + plt.semilogy(criter_fgp) + plt.show() + """ + """ + print ("Run CGLS for least squares") + opt = {'tol': 1e-4, 'iter': 20} + x_init = ImageData(geometry=ig) + x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, DataContainer(sino_channel), opt=opt) + + plt.figure() + plt.imshow(x_CGLS.array,vmin=0, vmax=0.05) + plt.title('CGLS') + plt.show() + """ +# Saving images into fits using astrapy if required +counter = 0 +filename = 'FISTA_TV_imat_slice' +for i in range(totChannels): + im = REC_chan[i,:,:] + add_val = np.min(im[:]) + im += abs(add_val) + im = im/np.max(im[:])*65535 + outfile = '{!s}_{!s}_{!s}.fits'.format(filename,str(selectedVertical_slice),str(counter)) + hdu = fits.PrimaryHDU(np.uint16(im)) + hdu.writeto(outfile, overwrite=True) + counter += 1
\ No newline at end of file diff --git a/Wrappers/Python/wip/old_demos/demo_imat_whitebeam.py b/Wrappers/Python/wip/old_demos/demo_imat_whitebeam.py new file mode 100644 index 0000000..e0d213e --- /dev/null +++ b/Wrappers/Python/wip/old_demos/demo_imat_whitebeam.py @@ -0,0 +1,138 @@ +# This script demonstrates how to load IMAT fits data +# into the CIL optimisation framework and run reconstruction methods. +# +# This demo loads the summedImg files which are the non-spectral images +# resulting from summing projections over all spectral channels. + +# needs dxchange: conda install -c conda-forge dxchange +# needs astropy: conda install astropy + + +# All third-party imports. +import numpy +from scipy.io import loadmat +import matplotlib.pyplot as plt +from dxchange.reader import read_fits + +# All own imports. +from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple +from ccpi.optimisation.algs import CGLS, FISTA +from ccpi.optimisation.funcs import Norm2sq, Norm1 + +# Load and display a couple of summed projection as examples +pathname0 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/' +filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' + +data0 = read_fits(pathname0 + filename0) + +pathname10 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/' +filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' + +data10 = read_fits(pathname10 + filename10) + +# Load a flat field (more are available, should we average over them?) +flat1 = read_fits('/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') + +# Apply flat field and display after flat-field correction and negative log +data0_rel = numpy.zeros(numpy.shape(flat1), dtype = float) +nonzero = flat1 > 0 +data0_rel[nonzero] = data0[nonzero] / flat1[nonzero] +data10_rel = numpy.zeros(numpy.shape(flat1), dtype = float) +data10_rel[nonzero] = data10[nonzero] / flat1[nonzero] + +plt.imshow(data0_rel) +plt.colorbar() +plt.show() + +plt.imshow(-numpy.log(data0_rel)) +plt.colorbar() +plt.show() + +plt.imshow(data10_rel) +plt.colorbar() +plt.show() + +plt.imshow(-numpy.log(data10_rel)) +plt.colorbar() +plt.show() + +# Set up for loading all summed images at 250 angles. +pathname = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/' +filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' + +# Dimensions +num_angles = 250 +imsize = 512 + +# Initialise array +data = numpy.zeros((num_angles,imsize,imsize)) + +# Load only 0-249, as 250 is at repetition of zero degrees just like image 0 +for i in range(0,250): + curimfile = (pathname + filename).format(i, i+4675) + data[i,:,:] = read_fits(curimfile) + +# Apply flat field and take negative log +nonzero = flat1 > 0 +for i in range(0,250): + data[i,nonzero] = data[i,nonzero]/flat1[nonzero] + +eqzero = data == 0 +data[eqzero] = 1 + +data_rel = -numpy.log(data) + +# Permute order to get: angles, vertical, horizontal, as default in framework. +data_rel = numpy.transpose(data_rel,(0,2,1)) + +# Set angles to use +angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False) + +# Create 3D acquisition geometry and acquisition data +ag = AcquisitionGeometry('parallel', + '3D', + angles, + pixel_num_h=imsize, + pixel_num_v=imsize) +b = AcquisitionData(data_rel, geometry=ag) + +# Reduce to single (noncentral) slice by extracting relevant parameters from data and its +# geometry. Perhaps create function to extract central slice automatically? +b2d = b.subset(vertical=128) +ag2d = AcquisitionGeometry('parallel', + '2D', + ag.angles, + pixel_num_h=ag.pixel_num_h) +b2d.geometry = ag2d + +# Create 2D image geometry +ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, + voxel_num_y=ag2d.pixel_num_h) + +# Create GPU projector/backprojector operator with ASTRA. +Aop = AstraProjectorSimple(ig2d,ag2d,'gpu') + +# Demonstrate operator is working by applying simple backprojection. +z = Aop.adjoint(b2d) +plt.imshow(z.array) +plt.title('Simple backprojection') +plt.colorbar() +plt.show() + +# Set initial guess ImageData with zeros for algorithms, and algorithm options. +x_init = ImageData(numpy.zeros((imsize,imsize)), + geometry=ig2d) +opt_CGLS = {'tol': 1e-4, 'iter': 20} + +# Run CGLS algorithm and display reconstruction. +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b2d, opt_CGLS) + +plt.imshow(x_CGLS.array) +plt.title('CGLS') +plt.colorbar() +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS Criterion vs iterations') +plt.show()
\ No newline at end of file diff --git a/Wrappers/Python/wip/old_demos/demo_memhandle.py b/Wrappers/Python/wip/old_demos/demo_memhandle.py new file mode 100755 index 0000000..db48d73 --- /dev/null +++ b/Wrappers/Python/wip/old_demos/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() diff --git a/Wrappers/Python/wip/old_demos/demo_test_sirt.py b/Wrappers/Python/wip/old_demos/demo_test_sirt.py new file mode 100644 index 0000000..6f5a44d --- /dev/null +++ b/Wrappers/Python/wip/old_demos/demo_test_sirt.py @@ -0,0 +1,176 @@ +# This demo illustrates how to use the SIRT algorithm without and with +# nonnegativity and box constraints. The ASTRA 2D projectors are used. + +# First make all imports +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \ + AcquisitionData +from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT +from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D, IndicatorBox +from ccpi.astra.ops import AstraProjectorSimple + +import numpy as np +import matplotlib.pyplot as plt + +# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case +test_case = 1 + +# 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 = 128 +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() + +# Set up AcquisitionGeometry object to hold the parameters of the measurement +# setup geometry: # Number of angles, the actual angles from 0 to +# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector +# pixel relative to an object pixel, the number of detector pixels, and the +# source-origin and origin-detector distance (here the origin-detector distance +# set to 0 to simulate a "virtual detector" with same detector pixel size as +# object pixel size). +angles_num = 20 +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) + ag = AcquisitionGeometry('parallel', + '2D', + angles, + det_num,det_w) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) + ag = AcquisitionGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) +else: + NotImplemented + +# Set up Operator object combining the ImageGeometry and AcquisitionGeometry +# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU. +Aop = AstraProjectorSimple(ig, ag, 'gpu') + +# Forward and backprojection are available as methods direct and adjoint. Here +# generate test data b and do simple backprojection to obtain z. +b = Aop.direct(Phantom) +z = Aop.adjoint(b) + +plt.imshow(b.array) +plt.title('Simulated data') +plt.show() + +plt.imshow(z.array) +plt.title('Backprojected data') +plt.show() + +# Using the test data b, different reconstruction methods can now be set up as +# demonstrated in the rest of this file. In general all methods need an initial +# guess and some algorithm options to be set: +x_init = ImageData(np.zeros(x.shape),geometry=ig) +opt = {'tol': 1e-4, 'iter': 1000} + +# First a CGLS reconstruction can be done: +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) + +plt.imshow(x_CGLS.array) +plt.title('CGLS') +plt.colorbar() +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS criterion') +plt.show() + +# A SIRT unconstrained reconstruction can be done: similarly: +x_SIRT, it_SIRT, timing_SIRT, criter_SIRT = SIRT(x_init, Aop, b, opt) + +plt.imshow(x_SIRT.array) +plt.title('SIRT unconstrained') +plt.colorbar() +plt.show() + +plt.semilogy(criter_SIRT) +plt.title('SIRT unconstrained criterion') +plt.show() + +# A SIRT nonnegativity constrained reconstruction can be done using the +# additional input "constraint" set to a box indicator function with 0 as the +# lower bound and the default upper bound of infinity: +x_SIRT0, it_SIRT0, timing_SIRT0, criter_SIRT0 = SIRT(x_init, Aop, b, opt, + constraint=IndicatorBox(lower=0)) + +plt.imshow(x_SIRT0.array) +plt.title('SIRT nonneg') +plt.colorbar() +plt.show() + +plt.semilogy(criter_SIRT0) +plt.title('SIRT nonneg criterion') +plt.show() + +# A SIRT reconstruction with box constraints on [0,1] can also be done: +x_SIRT01, it_SIRT01, timing_SIRT01, criter_SIRT01 = SIRT(x_init, Aop, b, opt, + constraint=IndicatorBox(lower=0,upper=1)) + +plt.imshow(x_SIRT01.array) +plt.title('SIRT box(0,1)') +plt.colorbar() +plt.show() + +plt.semilogy(criter_SIRT01) +plt.title('SIRT box(0,1) criterion') +plt.show() + +# The indicator function can also be used with the FISTA algorithm to do +# least squares with nonnegativity constraint. + +# Create least squares object instance with projector, test data and a constant +# coefficient of 0.5: +f = Norm2sq(Aop,b,c=0.5) + +# Run FISTA for least squares without constraints +x_fista, it, timing, criter = FISTA(x_init, f, None,opt) + +plt.imshow(x_fista.array) +plt.title('FISTA Least squares') +plt.show() + +plt.semilogy(criter) +plt.title('FISTA Least squares criterion') +plt.show() + +# Run FISTA for least squares with nonnegativity constraint +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, IndicatorBox(lower=0),opt) + +plt.imshow(x_fista0.array) +plt.title('FISTA Least squares nonneg') +plt.show() + +plt.semilogy(criter0) +plt.title('FISTA Least squares nonneg criterion') +plt.show() + +# Run FISTA for least squares with box constraint [0,1] +x_fista01, it01, timing01, criter01 = FISTA(x_init, f, IndicatorBox(lower=0,upper=1),opt) + +plt.imshow(x_fista01.array) +plt.title('FISTA Least squares box(0,1)') +plt.show() + +plt.semilogy(criter01) +plt.title('FISTA Least squares box(0,1) criterion') +plt.show()
\ No newline at end of file diff --git a/Wrappers/Python/wip/old_demos/multifile_nexus.py b/Wrappers/Python/wip/old_demos/multifile_nexus.py new file mode 100755 index 0000000..d1ad438 --- /dev/null +++ b/Wrappers/Python/wip/old_demos/multifile_nexus.py @@ -0,0 +1,307 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug 15 16:00:53 2018
+
+@author: ofn77899
+"""
+
+import os
+from ccpi.io.reader import NexusReader
+
+from sys import getsizeof
+
+import matplotlib.pyplot as plt
+
+from ccpi.framework import DataProcessor, DataContainer
+from ccpi.processors import Normalizer
+from ccpi.processors import CenterOfRotationFinder
+import numpy
+
+
+class averager(object):
+ def __init__(self):
+ self.reset()
+
+ def reset(self):
+ self.N = 0
+ self.avg = 0
+ self.min = 0
+ self.max = 0
+ self.var = 0
+ self.ske = 0
+ self.kur = 0
+
+ def add_reading(self, val):
+
+ if (self.N == 0):
+ self.avg = val;
+ self.min = val;
+ self.max = val;
+ elif (self.N == 1):
+ #//set min/max
+ self.max = val if val > self.max else self.max
+ self.min = val if val < self.min else self.min
+
+
+ thisavg = (self.avg + val)/2
+ #// initial value is in avg
+ self.var = (self.avg - thisavg)*(self.avg-thisavg) + (val - thisavg) * (val-thisavg)
+ self.ske = self.var * (self.avg - thisavg)
+ self.kur = self.var * self.var
+ self.avg = thisavg
+ else:
+ self.max = val if val > self.max else self.max
+ self.min = val if val < self.min else self.min
+
+ M = self.N
+
+ #// b-factor =(<v>_N + v_(N+1)) / (N+1)
+ #float b = (val -avg) / (M+1);
+ b = (val -self.avg) / (M+1)
+
+ self.kur = self.kur + (M *(b*b*b*b) * (1+M*M*M))- (4* b * self.ske) + (6 * b *b * self.var * (M-1))
+
+ self.ske = self.ske + (M * b*b*b *(M-1)*(M+1)) - (3 * b * self.var * (M-1))
+
+ #//var = var * ((M-1)/M) + ((val - avg)*(val - avg)/(M+1)) ;
+ self.var = self.var * ((M-1)/M) + (b * b * (M+1))
+
+ self.avg = self.avg * (M/(M+1)) + val / (M+1)
+
+ self.N += 1
+
+ def stats(self, vector):
+ i = 0
+ while i < vector.size:
+ self.add_reading(vector[i])
+ i+=1
+
+avg = averager()
+a = numpy.linspace(0,39,40)
+avg.stats(a)
+print ("average" , avg.avg, a.mean())
+print ("variance" , avg.var, a.var())
+b = a - a.mean()
+b *= b
+b = numpy.sqrt(sum(b)/(a.size-1))
+print ("std" , numpy.sqrt(avg.var), b)
+#%%
+
+class DataStatMoments(DataProcessor):
+ '''Normalization based on flat and dark
+
+ This processor read in a AcquisitionData and normalises it based on
+ the instrument reading with and without incident photons or neutrons.
+
+ Input: AcquisitionData
+ Parameter: 2D projection with flat field (or stack)
+ 2D projection with dark field (or stack)
+ Output: AcquisitionDataSetn
+ '''
+
+ def __init__(self, axis, skewness=False, kurtosis=False, offset=0):
+ kwargs = {
+ 'axis' : axis,
+ 'skewness' : skewness,
+ 'kurtosis' : kurtosis,
+ 'offset' : offset,
+ }
+ #DataProcessor.__init__(self, **kwargs)
+ super(DataStatMoments, self).__init__(**kwargs)
+
+
+ def check_input(self, dataset):
+ #self.N = dataset.get_dimension_size(self.axis)
+ return True
+
+ @staticmethod
+ def add_sample(dataset, N, axis, stats=None, offset=0):
+ #dataset = self.get_input()
+ if (N == 0):
+ # get the axis number along to calculate the stats
+
+
+ #axs = dataset.get_dimension_size(self.axis)
+ # create a placeholder for the output
+ if stats is None:
+ ax = dataset.get_dimension_axis(axis)
+ shape = [dataset.shape[i] for i in range(len(dataset.shape)) if i != ax]
+ # output has 4 components avg, std, skewness and kurtosis + last avg+ (val-thisavg)
+ shape.insert(0, 4+2)
+ stats = numpy.zeros(shape)
+
+ stats[0] = dataset.subset(**{axis:N-offset}).array[:]
+
+ #avg = val
+ elif (N == 1):
+ # val
+ stats[5] = dataset.subset(**{axis:N-offset}).array
+ stats[4] = stats[0] + stats[5]
+ stats[4] /= 2 # thisavg
+ stats[5] -= stats[4] # (val - thisavg)
+
+ #float thisavg = (avg + val)/2;
+
+ #// initial value is in avg
+ #var = (avg - thisavg)*(avg-thisavg) + (val - thisavg) * (val-thisavg);
+ stats[1] = stats[5] * stats[5] + stats[5] * stats[5]
+ #skewness = var * (avg - thisavg);
+ stats[2] = stats[1] * stats[5]
+ #kurtosis = var * var;
+ stats[3] = stats[1] * stats[1]
+ #avg = thisavg;
+ stats[0] = stats[4]
+ else:
+
+ #float M = (float)N;
+ M = N
+ #val
+ stats[4] = dataset.subset(**{axis:N-offset}).array
+ #// b-factor =(<v>_N + v_(N+1)) / (N+1)
+ stats[5] = stats[4] - stats[0]
+ stats[5] /= (M+1) # b factor
+ #float b = (val -avg) / (M+1);
+
+ #kurtosis = kurtosis + (M *(b*b*b*b) * (1+M*M*M))- (4* b * skewness) + (6 * b *b * var * (M-1));
+ #if self.kurtosis:
+ # stats[3] += (M * stats[5] * stats[5] * stats[5] * stats[5]) - \
+ # (4 * stats[5] * stats[2]) + \
+ # (6 * stats[5] * stats[5] * stats[1] * (M-1))
+
+ #skewness = skewness + (M * b*b*b *(M-1)*(M+1)) - (3 * b * var * (M-1));
+ #if self.skewness:
+ # stats[2] = stats[2] + (M * stats[5]* stats[5] * stats[5] * (M-1)*(M-1) ) -\
+ # 3 * stats[5] * stats[1] * (M-1)
+ #//var = var * ((M-1)/M) + ((val - avg)*(val - avg)/(M+1)) ;
+ #var = var * ((M-1)/M) + (b * b * (M+1));
+ stats[1] = ((M-1)/M) * stats[1] + (stats[5] * stats[5] * (M+1))
+
+ #avg = avg * (M/(M+1)) + val / (M+1)
+ stats[0] = stats[0] * (M/(M+1)) + stats[4] / (M+1)
+
+ N += 1
+ return stats , N
+
+
+ def process(self):
+
+ data = self.get_input()
+
+ #stats, i = add_sample(0)
+ N = data.get_dimension_size(self.axis)
+ ax = data.get_dimension_axis(self.axis)
+ stats = None
+ i = 0
+ while i < N:
+ stats , i = DataStatMoments.add_sample(data, i, self.axis, stats, offset=self.offset)
+ self.offset += N
+ labels = ['StatisticalMoments']
+
+ labels += [data.dimension_labels[i] \
+ for i in range(len(data.dimension_labels)) if i != ax]
+ y = DataContainer( stats[:4] , False,
+ dimension_labels=labels)
+ return y
+
+directory = r'E:\Documents\Dataset\CCPi\Nexus_test'
+data_path="entry1/instrument/pco1_hw_hdf_nochunking/data"
+
+reader = NexusReader(os.path.join( os.path.abspath(directory) , '74331.nxs'))
+
+print ("read flat")
+read_flat = NexusReader(os.path.join( os.path.abspath(directory) , '74240.nxs'))
+read_flat.data_path = data_path
+flatsslice = read_flat.get_acquisition_data_whole()
+avg = DataStatMoments('angle')
+
+avg.set_input(flatsslice)
+flats = avg.get_output()
+
+ave = averager()
+ave.stats(flatsslice.array[:,0,0])
+
+print ("avg" , ave.avg, flats.array[0][0][0])
+print ("var" , ave.var, flats.array[1][0][0])
+
+print ("read dark")
+read_dark = NexusReader(os.path.join( os.path.abspath(directory) , '74243.nxs'))
+read_dark.data_path = data_path
+
+## darks are very many so we proceed in batches
+total_size = read_dark.get_projection_dimensions()[0]
+print ("total_size", total_size)
+
+batchsize = 40
+if batchsize > total_size:
+ batchlimits = [batchsize * (i+1) for i in range(int(total_size/batchsize))] + [total_size-1]
+else:
+ batchlimits = [total_size]
+#avg.N = 0
+avg.offset = 0
+N = 0
+for batch in range(len(batchlimits)):
+ print ("running batch " , batch)
+ bmax = batchlimits[batch]
+ bmin = bmax-batchsize
+
+ darksslice = read_dark.get_acquisition_data_batch(bmin,bmax)
+ if batch == 0:
+ #create stats
+ ax = darksslice.get_dimension_axis('angle')
+ shape = [darksslice.shape[i] for i in range(len(darksslice.shape)) if i != ax]
+ # output has 4 components avg, std, skewness and kurtosis + last avg+ (val-thisavg)
+ shape.insert(0, 4+2)
+ print ("create stats shape ", shape)
+ stats = numpy.zeros(shape)
+ print ("N" , N)
+ #avg.set_input(darksslice)
+ i = bmin
+ while i < bmax:
+ stats , i = DataStatMoments.add_sample(darksslice, i, 'angle', stats, bmin)
+ print ("{0}-{1}-{2}".format(bmin, i , bmax ) )
+
+darks = stats
+#%%
+
+fig = plt.subplot(2,2,1)
+fig.imshow(flats.subset(StatisticalMoments=0).array)
+fig = plt.subplot(2,2,2)
+fig.imshow(numpy.sqrt(flats.subset(StatisticalMoments=1).array))
+fig = plt.subplot(2,2,3)
+fig.imshow(darks[0])
+fig = plt.subplot(2,2,4)
+fig.imshow(numpy.sqrt(darks[1]))
+plt.show()
+
+
+#%%
+norm = Normalizer(flat_field=flats.array[0,200,:], dark_field=darks[0,200,:])
+#norm.set_flat_field(flats.array[0,200,:])
+#norm.set_dark_field(darks.array[0,200,:])
+norm.set_input(reader.get_acquisition_data_slice(200))
+
+n = Normalizer.normalize_projection(norm.get_input().as_array(), flats.array[0,200,:], darks[0,200,:], 1e-5)
+#dn_n= Normalizer.estimate_normalised_error(norm.get_input().as_array(), flats.array[0,200,:], darks[0,200,:],
+# numpy.sqrt(flats.array[1,200,:]), numpy.sqrt(darks[1,200,:]))
+#%%
+
+
+#%%
+fig = plt.subplot(2,1,1)
+
+
+fig.imshow(norm.get_input().as_array())
+fig = plt.subplot(2,1,2)
+fig.imshow(n)
+
+#fig = plt.subplot(3,1,3)
+#fig.imshow(dn_n)
+
+
+plt.show()
+
+
+
+
+
+
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index 1be3dfa..e123739 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -90,7 +90,7 @@ g = ZeroFunction() normK = operator.norm() ## Primal & dual stepsizes -diag_precon = False +diag_precon = True if diag_precon: diff --git a/Wrappers/Python/wip/test_pdhg_gap.py b/Wrappers/Python/wip/test_pdhg_gap.py deleted file mode 100644 index 6c7ccc9..0000000 --- a/Wrappers/Python/wip/test_pdhg_gap.py +++ /dev/null @@ -1,140 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Apr 2 12:26:24 2019 - -@author: vaggelis -""" - - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \ - MixedL21Norm, BlockFunction, ScaledFunction - -from ccpi.astra.ops import AstraProjectorSimple -from skimage.util import random_noise - - -#%%############################################################################### -# Create phantom for TV tomography - -#import os -#import tomophantom -#from tomophantom import TomoP2D -#from tomophantom.supp.qualitymetrics import QualityTools - -#model = 1 # select a model number from the library -#N = 150 # set dimension of the phantom -## one can specify an exact path to the parameters file -## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' -#path = os.path.dirname(tomophantom.__file__) -#path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -##This will generate a N_size x N_size phantom (2D) -#phantom_2D = TomoP2D.Model(model, N, path_library2D) -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -#data = ImageData(phantom_2D, geometry=ig) - -N = 150 -x = np.zeros((N,N)) -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 - -data = ImageData(x) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - - -detectors = 150 -angles = np.linspace(0,np.pi,100) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -plt.imshow(sin.as_array()) -plt.title('Sinogram') -plt.colorbar() -plt.show() - -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - - -#%% Works only with Composite Operator Structure of PDHG - -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Form Composite Operator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -alpha = 50 -f = BlockFunction( alpha * MixedL21Norm(), \ - 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFun() - -# Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes - -sigma = 10 -tau = 1/(sigma*normK**2) - -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 100 - -pdhg.run(5000) -#%% - -opt = {'niter':2000} - -res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) - -#%% -sol = pdhg.get_output().as_array() -sol_old = res[0].as_array() -fig = plt.figure(figsize=(20,10)) -plt.subplot(1,3,1) -plt.imshow(noisy_data.as_array()) -#plt.colorbar() -plt.subplot(1,3,2) -plt.imshow(sol) -#plt.colorbar() - -plt.subplot(1,3,3) -plt.imshow(sol_old) -plt.show() - -plt.imshow(np.abs(sol-sol_old)) -plt.colorbar() -plt.show() - - -# -# -##%% -#plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon') -#plt.legend() -#plt.show() - - diff --git a/Wrappers/Python/wip/test_profile.py b/Wrappers/Python/wip/test_profile.py deleted file mode 100644 index f14c0c3..0000000 --- a/Wrappers/Python/wip/test_profile.py +++ /dev/null @@ -1,84 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Mon Apr 8 13:57:46 2019 - -@author: evangelos -""" - -# profile direct, adjoint, gradient - -from ccpi.framework import ImageGeometry -from ccpi.optimisation.operators import Gradient, BlockOperator, Identity -from ccpi.optimisation.functions import MixedL21Norm, L2NormSquared, BlockFunction -import numpy - -N, M, K = 2, 3, 2 - -ig = ImageGeometry(N, M) -b = ig.allocate('random_int') - -G = Gradient(ig) -Id = Identity(ig) - -#operator = BlockOperator(G, Id) -operator = G - -f1 = MixedL21Norm() -f2 = L2NormSquared(b = b) - -f = BlockFunction( f1, f2) - - -x_old = operator.domain_geometry().allocate() -y_old = operator.range_geometry().allocate('random_int') - - -xbar = operator.domain_geometry().allocate('random_int') - -x_tmp = x_old.copy() -x = x_old.copy() - -y_tmp = operator.range_geometry().allocate() -y = y_old.copy() - -y1 = y.copy() - -sigma = 20 - -for i in range(100): - - operator.direct(xbar, out = y_tmp) - y_tmp *= sigma - y_tmp += y_old - - - y_tmp1 = sigma * operator.direct(xbar) + y_old - - print(i) - print(" y_old :", y_old[0].as_array(), "\n") - print(" y_tmp[0] :", y_tmp[0].as_array(),"\n") - print(" y_tmp1[0] :", y_tmp1[0].as_array()) - - - numpy.testing.assert_array_equal(y_tmp[0].as_array(), \ - y_tmp1[0].as_array()) - - numpy.testing.assert_array_equal(y_tmp[1].as_array(), \ - y_tmp1[1].as_array()) - - - y1 = f.proximal_conjugate(y_tmp1, sigma) - f.proximal_conjugate(y_tmp, sigma, y) - - - - - - - - - - - -
\ No newline at end of file diff --git a/Wrappers/Python/wip/test_symGrad_method1.py b/Wrappers/Python/wip/test_symGrad_method1.py deleted file mode 100644 index 36adee1..0000000 --- a/Wrappers/Python/wip/test_symGrad_method1.py +++ /dev/null @@ -1,178 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, \ - Gradient, SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -#def dt(steps): -# return steps[-1] - steps[-2] - -# Create phantom for TGV Gaussian denoising - -N = 3 - -data = np.zeros((N,N)) - -x1 = np.linspace(0, int(N/2), N) -x2 = np.linspace(int(N/2), 0., N) -xv, yv = np.meshgrid(x1, x2) - -xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T - -data = xv -data = data/data.max() - -plt.imshow(data) -plt.show() - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.005, seed=10) -noisy_data = ImageData(n1) - - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy data') -plt.show() - -alpha, beta = 0.2, 1 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '1' - - -# Create operators -op11 = Gradient(ig) -op12 = Identity(op11.range_geometry()) - -op22 = SymmetrizedGradient(op11.domain_geometry()) - -op21 = ZeroOperator(ig, op22.range_geometry()) - - -op31 = Identity(ig, ag) -op32 = ZeroOperator(op22.domain_geometry(), ag) - -operator1 = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) - - -f1 = alpha * MixedL21Norm() -f2 = beta * MixedL21Norm() -f3 = ZeroFunction() -f_B3 = BlockFunction(f1, f2, f3) -g_B3 = ZeroFunction() - - - -# Create operators -op11 = Gradient(ig) -op12 = Identity(op11.range_geometry()) - -op22 = SymmetrizedGradient(op11.domain_geometry()) - -op21 = ZeroOperator(ig, op22.range_geometry()) - -operator2 = BlockOperator(op11, -1*op12, \ - op21, op22, \ - shape=(2,2) ) - -#f1 = alpha * MixedL21Norm() -#f2 = beta * MixedL21Norm() -f_B2 = BlockFunction(f1, f2) -g_B2 = 0.5 * L2NormSquared(b = noisy_data) - - -#%% - -x_old1 = operator1.domain_geometry().allocate('random_int') -y_old1 = operator1.range_geometry().allocate() - -xbar1 = x_old1.copy() -x_tmp1 = x_old1.copy() -x1 = x_old1.copy() - -y_tmp1 = y_old1.copy() -y1 = y_tmp1.copy() - -x_old2 = x_old1.copy() -y_old2 = operator2.range_geometry().allocate() - -xbar2 = x_old2.copy() -x_tmp2 = x_old2.copy() -x2 = x_old2.copy() - -y_tmp2 = y_old2.copy() -y2 = y_tmp2.copy() - -sigma = 0.4 -tau = 0.4 - -y_tmp1 = y_old1 + sigma * operator1.direct(xbar1) -y_tmp2 = y_old2 + sigma * operator2.direct(xbar2) - -numpy.testing.assert_array_equal(y_tmp1[0][0].as_array(), y_tmp2[0][0].as_array()) -numpy.testing.assert_array_equal(y_tmp1[0][1].as_array(), y_tmp2[0][1].as_array()) -numpy.testing.assert_array_equal(y_tmp1[1][0].as_array(), y_tmp2[1][0].as_array()) -numpy.testing.assert_array_equal(y_tmp1[1][1].as_array(), y_tmp2[1][1].as_array()) - - -y1 = f_B3.proximal_conjugate(y_tmp1, sigma) -y2 = f_B2.proximal_conjugate(y_tmp2, sigma) - -numpy.testing.assert_array_equal(y1[0][0].as_array(), y2[0][0].as_array()) -numpy.testing.assert_array_equal(y1[0][1].as_array(), y2[0][1].as_array()) -numpy.testing.assert_array_equal(y1[1][0].as_array(), y2[1][0].as_array()) -numpy.testing.assert_array_equal(y1[1][1].as_array(), y2[1][1].as_array()) - - -x_tmp1 = x_old1 - tau * operator1.adjoint(y1) -x_tmp2 = x_old2 - tau * operator2.adjoint(y2) - -numpy.testing.assert_array_equal(x_tmp1[0].as_array(), x_tmp2[0].as_array()) - - - - - - - - - - - -############################################################################## -#x_1 = operator1.domain_geometry().allocate('random_int') -# -#x_2 = BlockDataContainer(x_1[0], x_1[1]) -# -#res1 = operator1.direct(x_1) -#res2 = operator2.direct(x_2) -# -#print(res1[0][0].as_array(), res2[0][0].as_array()) -#print(res1[0][1].as_array(), res2[0][1].as_array()) -# -#print(res1[1][0].as_array(), res2[1][0].as_array()) -#print(res1[1][1].as_array(), res2[1][1].as_array()) -# -##res1 = op11.direct(x1[0]) - op12.direct(x1[1]) -##res2 = op21.direct(x1[0]) - op22.direct(x1[1]) |