From d2b119495c42182530c4f0329613c24b32d395fa Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Tue, 2 Apr 2019 17:47:18 +0100 Subject: add methods for precond pdhg --- .../optimisation/operators/GradientOperator.py | 57 ++++++++- .../optimisation/operators/SparseFiniteDiff.py | 140 +++++++++++++++++++++ .../Python/ccpi/optimisation/operators/__init__.py | 2 + 3 files changed, 195 insertions(+), 4 deletions(-) create mode 100644 Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index ec14b8f..87723f0 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -8,9 +8,9 @@ Created on Fri Mar 1 22:50:04 2019 from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator from ccpi.optimisation.ops import PowerMethodNonsquare -from ccpi.framework import ImageData, ImageGeometry, BlockGeometry +from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer import numpy -from ccpi.optimisation.operators import FiniteDiff +from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff #%% @@ -45,7 +45,6 @@ class Gradient(LinearOperator): tmp = self.gm_range.allocate() - for i in range(tmp.shape[0]): tmp.get_item(i).fill(FiniteDiff(self.gm_domain, direction = self.ind[i], bnd_cond = self.bnd_cond).direct(x)) return tmp @@ -73,6 +72,56 @@ class Gradient(LinearOperator): def __rmul__(self, scalar): return ScaledOperator(self, scalar) + + 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__': - pass + M, N = 2, 3 + ig = ImageGeometry(M, N) + arr = ig.allocate('random_int' ) + + G_neum = Gradient(ig) + + d = G_neum.matrix() + print(d[1]) + + d1 = G_neum.sum_abs_row() + print(d1.as_array()) + + d2 = G_neum.sum_abs_col() + print(d2) + + + + diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py new file mode 100644 index 0000000..0fb5efb --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Apr 2 14:06:15 2019 + +@author: vaggelis +""" + +import scipy.sparse as sp +import numpy as np +from ccpi.framework import ImageData + +class SparseFiniteDiff(): + + def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'): + + super(SparseFiniteDiff, self).__init__() + self.gm_domain = gm_domain + self.gm_range = gm_range + self.direction = direction + self.bnd_cond = bnd_cond + + if self.gm_range is None: + self.gm_range = self.gm_domain + + self.get_dims = [i for i in gm_domain.shape] + + if self.direction + 1 > len(self.gm_domain.shape): + raise ValueError('Gradient directions more than geometry domain') + + def matrix(self): + + i = self.direction + + mat = sp.spdiags(np.vstack([-np.ones((1,self.get_dims[i])),np.ones((1,self.get_dims[i]))]), [0,1], self.get_dims[i], self.get_dims[i], format = 'lil') + + if self.bnd_cond == 'Neumann': + mat[-1,:] = 0 + elif self.bnd_cond == 'Periodic': + mat[-1,0] = 1 + + tmpGrad = mat if i == 0 else sp.eye(self.get_dims[0]) + + for j in range(1, self.gm_domain.length): + + tmpGrad = sp.kron(mat, tmpGrad ) if j == i else sp.kron(sp.eye(self.get_dims[j]), tmpGrad ) + + return tmpGrad + + def T(self): + return self.matrix().T + + def direct(self, x): + + x_asarr = x.as_array() + res = np.reshape( self.matrix() * x_asarr.flatten('F'), self.gm_domain.shape, 'F') + return type(x)(res) + + def adjoint(self, x): + + x_asarr = x.as_array() + res = np.reshape( self.matrix().T * x_asarr.flatten('F'), self.gm_domain.shape, 'F') + return type(x)(res) + + def sum_abs_row(self): + + return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))) + + def sum_abs_col(self): + + return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F'))) + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + from ccpi.optimisation.operators import FiniteDiff + + # 2D + M, N= 2, 3 + ig = ImageGeometry(M, N) + arr = ig.allocate('random_int') + + for i in [0,1]: + + # Neumann + sFD_neum = SparseFiniteDiff(ig, direction=i, bnd_cond='Neumann') + G_neum = FiniteDiff(ig, direction=i, bnd_cond='Neumann') + + # Periodic + sFD_per = SparseFiniteDiff(ig, direction=i, bnd_cond='Periodic') + G_per = FiniteDiff(ig, direction=i, bnd_cond='Periodic') + + u_neum_direct = G_neum.direct(arr) + u_neum_sp_direct = sFD_neum.direct(arr) + np.testing.assert_array_almost_equal(u_neum_direct.as_array(), u_neum_sp_direct.as_array(), decimal=4) + + u_neum_adjoint = G_neum.adjoint(arr) + u_neum_sp_adjoint = sFD_neum.adjoint(arr) + np.testing.assert_array_almost_equal(u_neum_adjoint.as_array(), u_neum_sp_adjoint.as_array(), decimal=4) + + u_per_direct = G_neum.direct(arr) + u_per_sp_direct = sFD_neum.direct(arr) + np.testing.assert_array_almost_equal(u_per_direct.as_array(), u_per_sp_direct.as_array(), decimal=4) + + u_per_adjoint = G_per.adjoint(arr) + u_per_sp_adjoint = sFD_per.adjoint(arr) + np.testing.assert_array_almost_equal(u_per_adjoint.as_array(), u_per_sp_adjoint.as_array(), decimal=4) + + # 3D + M, N, K = 2, 3, 4 + ig3D = ImageGeometry(M, N, K) + arr3D = ig3D.allocate('random_int') + + for i in [0,1,2]: + + # Neumann + sFD_neum3D = SparseFiniteDiff(ig3D, direction=i, bnd_cond='Neumann') + G_neum3D = FiniteDiff(ig3D, direction=i, bnd_cond='Neumann') + + # Periodic + sFD_per3D = SparseFiniteDiff(ig3D, direction=i, bnd_cond='Periodic') + G_per3D = FiniteDiff(ig3D, direction=i, bnd_cond='Periodic') + + u_neum_direct3D = G_neum3D.direct(arr3D) + u_neum_sp_direct3D = sFD_neum3D.direct(arr3D) + np.testing.assert_array_almost_equal(u_neum_direct3D.as_array(), u_neum_sp_direct3D.as_array(), decimal=4) + + u_neum_adjoint3D = G_neum3D.adjoint(arr3D) + u_neum_sp_adjoint3D = sFD_neum3D.adjoint(arr3D) + np.testing.assert_array_almost_equal(u_neum_adjoint3D.as_array(), u_neum_sp_adjoint3D.as_array(), decimal=4) + + u_per_direct3D = G_neum3D.direct(arr3D) + u_per_sp_direct3D = sFD_neum3D.direct(arr3D) + np.testing.assert_array_almost_equal(u_per_direct3D.as_array(), u_per_sp_direct3D.as_array(), decimal=4) + + u_per_adjoint3D = G_per3D.adjoint(arr3D) + u_per_sp_adjoint3D = sFD_per3D.adjoint(arr3D) + np.testing.assert_array_almost_equal(u_per_adjoint3D.as_array(), u_per_sp_adjoint3D.as_array(), decimal=4) + + \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py index 1c09faf..57d89ad 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py @@ -17,3 +17,5 @@ from .GradientOperator import Gradient from .SymmetrizedGradientOperator import SymmetrizedGradient from .IdentityOperator import Identity from .ZeroOperator import ZeroOp + +from .SparseFiniteDiff import SparseFiniteDiff -- cgit v1.2.3