summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Wrappers/Python/ccpi/framework.py14
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algs.py5
-rwxr-xr-xWrappers/Python/ccpi/optimisation/funcs.py2
-rwxr-xr-xWrappers/Python/ccpi/optimisation/ops.py5
-rwxr-xr-xWrappers/Python/conda-recipe/run_test.py28
5 files changed, 38 insertions, 16 deletions
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index 318eb92..e1a2dff 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -3,7 +3,7 @@
# Visual Analytics and Imaging System Group of the Science Technology
# Facilities Council, STFC
-# Copyright 2018 Edoardo Pasca
+# Copyright 2018-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.
@@ -26,6 +26,7 @@ import numpy
import sys
from datetime import timedelta, datetime
import warnings
+from functools import reduce
def find_key(dic, val):
"""return the key of dictionary dic given the value"""
@@ -720,10 +721,15 @@ class DataContainer(object):
## reductions
def sum(self, out=None, *args, **kwargs):
return self.as_array().sum(*args, **kwargs)
- def norm(self):
- '''return the norm of the DataContainer'''
- y = self.as_array()
+ def squared_norm(self):
+ '''return the squared euclidean norm of the DataContainer viewed as a vector'''
+ shape = self.shape
+ size = reduce(lambda x,y:x*y, shape, 1)
+ y = numpy.reshape(self.as_array(), (size, ))
return numpy.dot(y, y.conjugate())
+ def norm(self):
+ '''return the euclidean norm of the DataContainer viewed as a vector'''
+ return numpy.sqrt(self.squared_norm())
class ImageData(DataContainer):
'''DataContainer for holding 2D or 3D DataContainer'''
diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py
index a736277..15638a9 100755
--- a/Wrappers/Python/ccpi/optimisation/algs.py
+++ b/Wrappers/Python/ccpi/optimisation/algs.py
@@ -158,9 +158,8 @@ def FBPD(x_init, operator=None, constraint=None, data_fidelity=None,\
memopt = opt['memopts'] if 'memopts' in opt.keys() else False
# step-sizes
- tau = 2 / (data_fidelity.L + 2);
- sigma = (1/tau - data_fidelity.L/2) / regulariser.L;
-
+ tau = 2 / (data_fidelity.L + 2)
+ sigma = (1/tau - data_fidelity.L/2) / regulariser.L
inv_sigma = 1/sigma
# initialization
diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py
index c7a6143..9b9fc36 100755
--- a/Wrappers/Python/ccpi/optimisation/funcs.py
+++ b/Wrappers/Python/ccpi/optimisation/funcs.py
@@ -217,10 +217,10 @@ class ZeroFun(Function):
class Norm1(Function):
def __init__(self,gamma):
+ super(Norm1, self).__init__()
self.gamma = gamma
self.L = 1
self.sign_x = None
- super(Norm1, self).__init__()
def __call__(self,x,out=None):
if out is None:
diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py
index 450b084..b2f996d 100755
--- a/Wrappers/Python/ccpi/optimisation/ops.py
+++ b/Wrappers/Python/ccpi/optimisation/ops.py
@@ -209,10 +209,11 @@ def PowerMethodNonsquare(op,numiters , x0=None):
# Loop
for it in numpy.arange(numiters):
x1 = op.adjoint(op.direct(x0))
- x1norm = numpy.sqrt((x1*x1).sum())
+ #x1norm = numpy.sqrt((x1*x1).sum())
+ x1norm = x1.norm()
#print ("x0 **********" ,x0)
#print ("x1 **********" ,x1)
- s[it] = (x1*x0).sum() / (x0*x0).sum()
+ 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/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py
index 5bf6538..b52af55 100755
--- a/Wrappers/Python/conda-recipe/run_test.py
+++ b/Wrappers/Python/conda-recipe/run_test.py
@@ -19,6 +19,7 @@ from ccpi.optimisation.funcs import Norm2
from ccpi.optimisation.ops import LinearOperatorMatrix
from ccpi.optimisation.ops import TomoIdentity
from ccpi.optimisation.ops import Identity
+from ccpi.optimisation.ops import PowerMethodNonsquare
import numpy.testing
@@ -494,6 +495,9 @@ class TestDataContainer(unittest.TestCase):
except AssertionError as err:
res = False
print(err)
+ print("expected " , second)
+ print("actual " , first)
+
self.assertTrue(res)
def test_DataContainerChaining(self):
dc = self.create_DataContainer(256,256,256,1)
@@ -501,7 +505,13 @@ class TestDataContainer(unittest.TestCase):
dc.add(9,out=dc)\
.subtract(1,out=dc)
self.assertEqual(1+9-1,dc.as_array().flatten()[0])
-
+ def test_reduction(self):
+ print ("test reductions")
+ dc = self.create_DataContainer(2,2,2,value=1)
+ sqnorm = dc.squared_norm()
+ norm = dc.norm()
+ self.assertEqual(sqnorm, 8.0)
+ numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7)
@@ -643,7 +653,8 @@ class TestAlgorithms(unittest.TestCase):
else:
self.assertTrue(cvx_not_installable)
- def test_FBPD_Norm1_cvx(self):
+ def skip_test_FBPD_Norm1_cvx(self):
+ print ("test_FBPD_Norm1_cvx")
if not cvx_not_installable:
opt = {'memopt': True}
# Problem data.
@@ -663,12 +674,15 @@ class TestAlgorithms(unittest.TestCase):
# Regularization parameter
lam = 10
opt = {'memopt': True}
+ # Initial guess
+ x_init = DataContainer(np.random.randn(n, 1))
+
# Create object instances with the test data A and b.
f = Norm2sq(A, b, c=0.5, memopt=True)
+ f.L = PowerMethodNonsquare(A, 25, x_init)[0]
+ print ("Lipschitz", f.L)
g0 = ZeroFun()
- # Initial guess
- x_init = DataContainer(np.zeros((n, 1)))
# Create 1-norm object instance
g1 = Norm1(lam)
@@ -710,7 +724,7 @@ class TestAlgorithms(unittest.TestCase):
# 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.
- def test_FISTA_denoise_cvx(self):
+ def skip_test_FISTA_denoise_cvx(self):
if not cvx_not_installable:
opt = {'memopt': True}
N = 64
@@ -733,6 +747,8 @@ class TestAlgorithms(unittest.TestCase):
# Data fidelity term
f_denoise = Norm2sq(I, y, c=0.5, memopt=True)
+ x_init = ImageData(geometry=ig)
+ f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0]
# 1-norm regulariser
lam1_denoise = 1.0
@@ -759,7 +775,7 @@ class TestAlgorithms(unittest.TestCase):
# Compare to CVXPY
# Construct the problem.
- x1_denoise = Variable(N**2, 1)
+ x1_denoise = Variable(N**2)
objective1_denoise = Minimize(
0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1))
prob1_denoise = Problem(objective1_denoise)
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-04 15:06:33 +0000