From a9ec78c1669f460ebaf5227600f8b0c082cfcf56 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 30 Apr 2019 22:38:20 +0100 Subject: Fix dot product bug --- Wrappers/Python/ccpi/framework/framework.py | 3 +- Wrappers/Python/wip/compare_CGLS_algos.py | 127 ++++++++++++++++++++++++++++ 2 files changed, 129 insertions(+), 1 deletion(-) create mode 100644 Wrappers/Python/wip/compare_CGLS_algos.py (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index ffc91ae..e278795 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -765,7 +765,8 @@ class DataContainer(object): def dot(self, other, *args, **kwargs): '''return the inner product of 2 DataContainers viewed as vectors''' if self.shape == other.shape: - return numpy.dot(self.as_array().ravel(), other.as_array().ravel()) + return (self*other).sum() + #return numpy.dot(self.as_array().ravel(), other.as_array().ravel()) else: raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape)) diff --git a/Wrappers/Python/wip/compare_CGLS_algos.py b/Wrappers/Python/wip/compare_CGLS_algos.py new file mode 100644 index 0000000..333805d --- /dev/null +++ b/Wrappers/Python/wip/compare_CGLS_algos.py @@ -0,0 +1,127 @@ +# 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.astra.ops import AstraProjectorSimple + +from ccpi.optimisation.algorithms import CGLS as CGLSalg + +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.figure() +#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, 'cpu') + +# 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.figure() +#plt.imshow(b.array) +#plt.title('Simulated data') +#plt.show() + +#plt.figure() +#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': 7} + +# First a CGLS reconstruction using the function version of CGLS can be done: +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) + +#plt.figure() +#plt.imshow(x_CGLS.array) +#plt.title('CGLS') +#plt.colorbar() +#plt.show() + +#plt.figure() +#plt.semilogy(criter_CGLS) +#plt.title('CGLS criterion') +#plt.show() + + + +# Now CLGS using the algorithm class +CGLS_alg = CGLSalg() +CGLS_alg.set_up(x_init, Aop, b ) +CGLS_alg.max_iteration = 2000 +CGLS_alg.run(opt['iter']) +x_CGLS_alg = CGLS_alg.get_output() + +#plt.figure() +#plt.imshow(x_CGLS_alg.as_array()) +#plt.title('CGLS ALG') +#plt.colorbar() +#plt.show() + +#plt.figure() +#plt.semilogy(CGLS_alg.objective) +#plt.title('CGLS criterion') +#plt.show() + +print(criter_CGLS) +print(CGLS_alg.objective) + +print((x_CGLS - x_CGLS_alg).norm()) \ No newline at end of file -- cgit v1.2.3 From 6d1c24c8fc389365f2dd83ee480c63969b08ce9f Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:13:40 +0100 Subject: removed useless imports --- Wrappers/Python/ccpi/optimisation/algs.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index 89b5519..f5ba85e 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -20,13 +20,8 @@ import numpy import time -from ccpi.optimisation.functions import Function -from ccpi.optimisation.functions import ZeroFunction -from ccpi.framework import ImageData -from ccpi.framework import AcquisitionData -from ccpi.optimisation.spdhg import spdhg -from ccpi.optimisation.spdhg import KullbackLeibler -from ccpi.optimisation.spdhg import KullbackLeiblerConvexConjugate + + def FISTA(x_init, f=None, g=None, opt=None): '''Fast Iterative Shrinkage-Thresholding Algorithm -- cgit v1.2.3 From da581e6061ebe206e007fe4378cc9d449b67d791 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:14:02 +0100 Subject: reimplements dot product following discussion in #273 an implementation of the dot product is made where we rely on Python to choose an appropriate type for the result of the reduction (e.g. float64 if the data is float32) --- Wrappers/Python/ccpi/framework/framework.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index e278795..66420b9 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -764,13 +764,23 @@ class DataContainer(object): 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 self.shape == other.shape: - return (self*other).sum() - #return numpy.dot(self.as_array().ravel(), other.as_array().ravel()) + # return (self*other).sum() + if method == 'numpy': + return numpy.dot(self.as_array().ravel(), other.as_array()) + elif method == 'reduce': + # see https://github.com/vais-ral/CCPi-Framework/pull/273 + # notice that Python seems to be smart enough to use + # the appropriate type to hold the result of the reduction + sf = reduce(lambda x,y: x + y[0]*y[1], + zip(self.as_array().ravel(), + other.as_array().ravel()), + 0) + return sf else: raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape)) - - + -- cgit v1.2.3 From 67fa9b5a985912232d8b5aeeec12451ed55a1a7a Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:16:17 +0100 Subject: fix import --- Wrappers/Python/test/test_DataContainer.py | 5 +++++ Wrappers/Python/wip/compare_CGLS_algos.py | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py index 8e8ab87..e92d4c6 100755 --- a/Wrappers/Python/test/test_DataContainer.py +++ b/Wrappers/Python/test/test_DataContainer.py @@ -455,6 +455,11 @@ class TestDataContainer(unittest.TestCase): self.assertTrue(False) except ValueError as ve: self.assertTrue(True) + + print ("test dot numpy") + n0 = (ds0 * ds1).sum() + n1 = ds0.as_array().ravel().dot(ds1.as_array().ravel()) + self.assertEqual(n0, n1) diff --git a/Wrappers/Python/wip/compare_CGLS_algos.py b/Wrappers/Python/wip/compare_CGLS_algos.py index 333805d..119752c 100644 --- a/Wrappers/Python/wip/compare_CGLS_algos.py +++ b/Wrappers/Python/wip/compare_CGLS_algos.py @@ -5,7 +5,7 @@ from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \ AcquisitionData from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT -from ccpi.astra.ops import AstraProjectorSimple +from ccpi.astra.operators import AstraProjectorSimple from ccpi.optimisation.algorithms import CGLS as CGLSalg -- cgit v1.2.3 From 82d94d608ea639c0aa8aefb80cc97c5d8b1ba2cb Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 1 May 2019 16:29:46 +0100 Subject: raise error if method is not recognised --- Wrappers/Python/ccpi/framework/framework.py | 3 +++ 1 file changed, 3 insertions(+) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index 66420b9..7236e0e 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -765,6 +765,9 @@ class DataContainer(object): 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': -- cgit v1.2.3