diff options
Diffstat (limited to 'Wrappers/Python')
14 files changed, 302 insertions, 57 deletions
diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py index 8e55b67..21ef3f0 100644 --- a/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py @@ -52,6 +52,11 @@ class BlockDataContainer(object): def is_compatible(self, other):
'''basic check if the size of the 2 objects fit'''
+
+ for i in range(len(self.containers)):
+ if type(self.containers[i])==type(self):
+ self = self.containers[i]
+
if isinstance(other, Number):
return True
elif isinstance(other, list):
diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py index d336305..0f43155 100644 --- a/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py +++ b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py @@ -27,6 +27,10 @@ class BlockGeometry(object): 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):
containers = [geom.allocate(value) for geom in self.geometries]
diff --git a/Wrappers/Python/build/lib/ccpi/framework/__init__.py b/Wrappers/Python/build/lib/ccpi/framework/__init__.py index 66e2f56..229edb5 100644 --- a/Wrappers/Python/build/lib/ccpi/framework/__init__.py +++ b/Wrappers/Python/build/lib/ccpi/framework/__init__.py @@ -15,6 +15,7 @@ from datetime import timedelta, datetime import warnings
from functools import reduce
+
from .framework import DataContainer
from .framework import ImageData, AcquisitionData
from .framework import ImageGeometry, AcquisitionGeometry
diff --git a/Wrappers/Python/build/lib/ccpi/framework/framework.py b/Wrappers/Python/build/lib/ccpi/framework/framework.py index ae9faf7..07c2ead 100644 --- a/Wrappers/Python/build/lib/ccpi/framework/framework.py +++ b/Wrappers/Python/build/lib/ccpi/framework/framework.py @@ -29,6 +29,7 @@ import warnings from functools import reduce from numbers import Number + def find_key(dic, val): """return the key of dictionary dic given the value""" return [k for k, v in dic.items() if v == val][0] @@ -496,6 +497,7 @@ class DataContainer(object): ## algebra def __add__(self, other, *args, **kwargs): out = kwargs.get('out', None) + if issubclass(type(other), DataContainer): if self.check_dimensions(other): out = self.as_array() + other.as_array() @@ -601,6 +603,7 @@ class DataContainer(object): deep_copy=True, dimension_labels=self.dimension_labels, geometry=self.geometry) + else: raise TypeError('Cannot {0} DataContainer with {1}'.format("multiply" , type(other))) diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py index 043fe38..d0e27ae 100644 --- a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py @@ -8,11 +8,14 @@ Created on Mon Feb 4 16:18:06 2019 from ccpi.optimisation.algorithms import Algorithm from ccpi.framework import ImageData import numpy as np -#import matplotlib.pyplot as plt +import matplotlib.pyplot as plt import time from ccpi.optimisation.operators import BlockOperator from ccpi.framework import BlockDataContainer + +import matplotlib.pyplot as plt + class PDHG(Algorithm): '''Primal Dual Hybrid Gradient''' @@ -69,9 +72,10 @@ class PDHG(Algorithm): self.xbar *= self.theta self.xbar += self.x - self.x_old.fill(self.x) - self.y_old.fill(self.y) - #self.y_old = y.copy() +# self.x_old.fill(self.x) +# self.y_old.fill(self.y) + self.y_old = self.y.copy() + self.x_old = self.x.copy() #self.y = self.y_old def update_objective(self): @@ -80,3 +84,72 @@ class PDHG(Algorithm): ]) + +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 + x_tmp = x_old + x = x_old + + y_tmp = y_old + y = y_tmp + + # relaxation parameter + theta = 1 + + t = time.time() + + objective = [] + + + for i in range(niter): + + # Gradient descent, Dual problem solution + y_tmp = y_old + sigma * operator.direct(xbar) + y = f.proximal_conjugate(y_tmp, sigma) + + # Gradient ascent, Primal problem solution + x_tmp = x_old - tau * operator.adjoint(y) + x = g.proximal(x_tmp, tau) + + #Update + xbar = x + theta * (x - x_old) + + x_old = x + y_old = y + + if i%100==0: + + primal = f(operator.direct(x)) + g(x) + dual = -(f.convex_conjugate(y) + g(-1*operator.adjoint(y))) + print( i, primal, dual, primal-dual) + +# plt.imshow(x.as_array()) +# plt.show() +# print(f(operator.direct(x)) + g(x), i) + + t_end = time.time() + + return x, t_end - t, objective + + + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py index 443bc78..f562973 100644 --- a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py @@ -28,3 +28,5 @@ from .GradientDescent import GradientDescent from .FISTA import FISTA from .FBPD import FBPD from .PDHG import PDHG +from .PDHG import PDHG_old + diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py index 70216a9..81c16cd 100644 --- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py @@ -10,6 +10,7 @@ import numpy as np #from ccpi.optimisation.funcs import Function from ccpi.optimisation.functions import Function from ccpi.framework import BlockDataContainer +from numbers import Number class BlockFunction(Function): '''A Block vector of Functions @@ -52,16 +53,24 @@ class BlockFunction(Function): def proximal_conjugate(self, x, tau, out = None): '''proximal_conjugate does not take into account the BlockOperator''' out = [None]*self.length - for i in range(self.length): - out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau) + 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): '''proximal does not take into account the BlockOperator''' out = [None]*self.length - for i in range(self.length): - out[i] = self.functions[i].proximal(x.get_item(i), tau) + 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) diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py index 5489d92..889d703 100644 --- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py @@ -1,12 +1,21 @@ # -*- 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 -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Feb 7 13:10:56 2019 +# 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 -@author: evangelos -""" +# 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 @@ -71,14 +80,15 @@ class L2NormSquared(Function): tmp = 0 if self.b is not None: - tmp = (self.b * x).sum() +# tmp = (self.b * x).sum() + tmp = (x * self.b).sum() if out is None: # FIXME: this is a number - return (1/4) * x.squared_norm() + tmp + return (1./4.) * x.squared_norm() + tmp else: # FIXME: this is a DataContainer - out.fill((1/4) * x.squared_norm() + tmp) + out.fill((1./4.) * x.squared_norm() + tmp) def proximal(self, x, tau, out = None): @@ -108,7 +118,8 @@ class L2NormSquared(Function): if out is None: if self.b is not None: - return (x - tau*self.b)/(1 + tau/2) + # change the order cannot add ImageData + NestedBlock + return (-1* tau*self.b + x)/(1 + tau/2) else: return x/(1 + tau/2 ) else: diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py index ec14b8f..e00de0c 100644 --- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/build/lib/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,115 @@ 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 + + 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 +# diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py index 52c7c3b..a58a296 100644 --- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py @@ -50,11 +50,11 @@ class Identity(LinearOperator): def sum_abs_row(self): - return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))) + return self.gm_domain.allocate(1)#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'))) + return self.gm_domain.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F'))) if __name__ == '__main__': diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py index 0fb5efb..1b88cba 100644 --- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py +++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py @@ -64,11 +64,15 @@ class SparseFiniteDiff(): def sum_abs_row(self): - return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))) + res = np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')) + res[res==0]=1 + return ImageData(res) def sum_abs_col(self): - return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F'))) + res = np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'C')) + res[res==0]=1 + return ImageData(res) if __name__ == '__main__': diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py index a28c0bf..f562973 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -29,3 +29,4 @@ from .FISTA import FISTA from .FBPD import FBPD from .PDHG import PDHG from .PDHG import PDHG_old + diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py index 2e0b9f4..d5c021d 100644 --- a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py @@ -74,11 +74,12 @@ else: ########################################################################### #%% -diag_precon = False +diag_precon = True if diag_precon: def tau_sigma_precond(operator): + tau = 1/operator.sum_abs_row() sigma = 1/ operator.sum_abs_col() diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index 8bd63c2..9feb05b 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -26,10 +26,10 @@ 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 +#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 @@ -55,7 +55,7 @@ detectors = 150 angles = np.linspace(0,np.pi,100) ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'gpu') +Aop = AstraProjectorSimple(ig, ag, 'cpu') sin = Aop.direct(data) plt.imshow(sin.as_array()) @@ -95,32 +95,55 @@ g = ZeroFun() 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 = 5000 -pdhg.update_objective_interval = 250 - -pdhg.run(5000) - -#%% -sol = pdhg.get_output().as_array() -fig = plt.figure() -plt.subplot(1,2,1) -plt.imshow(noisy_data.as_array()) -#plt.colorbar() -plt.subplot(1,2,2) -plt.imshow(sol) -#plt.colorbar() -plt.show() - +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: + sigma = 10 + tau = 1/(sigma*normK**2) + +#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) +#pdhg.max_iteration = 5000 +#pdhg.update_objective_interval = 250 +# +#pdhg.run(5000) + +opt = {'niter':2000} +# +res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) + +aaa = res[0].as_array() +# +plt.imshow(aaa) +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() +#sol = pdhg.get_output().as_array() +#fig = plt.figure() +#plt.subplot(1,2,1) +#plt.imshow(noisy_data.as_array()) +##plt.colorbar() +#plt.subplot(1,2,2) +#plt.imshow(sol) +##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() |