summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers')
-rwxr-xr-xWrappers/Python/ccpi/framework/Vector.py96
-rwxr-xr-xWrappers/Python/ccpi/framework/__init__.py56
-rwxr-xr-xWrappers/Python/ccpi/framework/framework.py2
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/CGLS.py49
-rwxr-xr-xWrappers/Python/ccpi/optimisation/functions/Norm2Sq.py1
-rwxr-xr-xWrappers/Python/ccpi/optimisation/functions/ScaledFunction.py3
-rw-r--r--Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py17
-rwxr-xr-xWrappers/Python/test/test_DataContainer.py53
-rwxr-xr-xWrappers/Python/test/test_run_test.py19
-rw-r--r--Wrappers/Python/wip/compare_CGLS_algos.py28
-rwxr-xr-xWrappers/Python/wip/fix_test.py208
11 files changed, 450 insertions, 82 deletions
diff --git a/Wrappers/Python/ccpi/framework/Vector.py b/Wrappers/Python/ccpi/framework/Vector.py
new file mode 100755
index 0000000..542cb7a
--- /dev/null
+++ b/Wrappers/Python/ccpi/framework/Vector.py
@@ -0,0 +1,96 @@
+# -*- 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 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
+import sys
+from datetime import timedelta, datetime
+import warnings
+from functools import reduce
+from numbers import Number
+from ccpi.framework import DataContainer
+
+class VectorData(DataContainer):
+ def __init__(self, array=None, **kwargs):
+ self.geometry = kwargs.get('geometry', None)
+ self.dtype = kwargs.get('dtype', numpy.float32)
+
+ if self.geometry is None:
+ if array is None:
+ raise ValueError('Please specify either a geometry or an array')
+ else:
+ if len(array.shape) > 1:
+ raise ValueError('Incompatible size: expected 1D got {}'.format(array.shape))
+ out = array
+ self.geometry = VectorGeometry(array.shape[0])
+ self.length = self.geometry.length
+ else:
+ self.length = self.geometry.length
+
+ if array is None:
+ out = numpy.zeros((self.length,), dtype=self.dtype)
+ else:
+ if self.length == array.shape[0]:
+ out = array
+ else:
+ raise ValueError('Incompatible size: expecting {} got {}'.format((self.length,), array.shape))
+ deep_copy = True
+ super(VectorData, self).__init__(out, deep_copy, None)
+
+class VectorGeometry(object):
+ RANDOM = 'random'
+ RANDOM_INT = 'random_int'
+
+ def __init__(self,
+ length):
+
+ self.length = length
+ self.shape = (length, )
+
+
+ def clone(self):
+ '''returns a copy of VectorGeometry'''
+ return VectorGeometry(self.length)
+
+ def allocate(self, value=0, **kwargs):
+ '''allocates an VectorData according to the size expressed in the instance'''
+ self.dtype = kwargs.get('dtype', numpy.float32)
+ out = VectorData(geometry=self, dtype=self.dtype)
+ if isinstance(value, Number):
+ if value != 0:
+ out += value
+ else:
+ if value == VectorGeometry.RANDOM:
+ seed = kwargs.get('seed', None)
+ if seed is not None:
+ numpy.random.seed(seed)
+ out.fill(numpy.random.random_sample(self.shape))
+ elif value == VectorGeometry.RANDOM_INT:
+ seed = kwargs.get('seed', None)
+ if seed is not None:
+ numpy.random.seed(seed)
+ max_value = kwargs.get('max_value', 100)
+ out.fill(numpy.random.randint(max_value,size=self.shape))
+ else:
+ raise ValueError('Value {} unknown'.format(value))
+ return out
diff --git a/Wrappers/Python/ccpi/framework/__init__.py b/Wrappers/Python/ccpi/framework/__init__.py
index 8926897..3de27ed 100755
--- a/Wrappers/Python/ccpi/framework/__init__.py
+++ b/Wrappers/Python/ccpi/framework/__init__.py
@@ -1,28 +1,28 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Mar 5 16:00:18 2019
-
-@author: ofn77899
-"""
-from __future__ import absolute_import
-from __future__ import division
-from __future__ import print_function
-from __future__ import unicode_literals
-
-import numpy
-import sys
-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
-from .framework import find_key, message
-from .framework import DataProcessor
-from .framework import AX, PixelByPixelDataProcessor, CastDataContainer
-from .BlockDataContainer import BlockDataContainer
-from .BlockGeometry import BlockGeometry
-
-from .TestData import TestData
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 5 16:00:18 2019
+
+@author: ofn77899
+"""
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy
+import sys
+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
+from .framework import find_key, message
+from .framework import DataProcessor
+from .framework import AX, PixelByPixelDataProcessor, CastDataContainer
+from .BlockDataContainer import BlockDataContainer
+from .BlockGeometry import BlockGeometry
+from .TestData import TestData
+from .Vector import VectorGeometry, VectorData
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index e906ca6..fd0e0d8 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -29,7 +29,6 @@ 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]
@@ -679,6 +678,7 @@ class DataContainer(object):
def pixel_wise_binary(self, pwop, x2, *args, **kwargs):
out = kwargs.get('out', None)
+
if out is None:
if isinstance(x2, (int, float, complex)):
out = pwop(self.as_array() , x2 , *args, **kwargs )
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
index 6b610a0..15acc31 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
@@ -24,6 +24,7 @@ Created on Thu Feb 21 11:11:23 2019
from ccpi.optimisation.algorithms import Algorithm
from ccpi.optimisation.functions import Norm2Sq
+import numpy
class CGLS(Algorithm):
@@ -49,13 +50,15 @@ class CGLS(Algorithm):
def set_up(self, x_init, operator , data ):
self.r = data.copy()
- self.x = x_init.copy()
+ self.x = x_init * 0
self.operator = operator
self.d = operator.adjoint(self.r)
self.normr2 = self.d.squared_norm()
+
+ self.s = self.operator.domain_geometry().allocate()
#if isinstance(self.normr2, Iterable):
# self.normr2 = sum(self.normr2)
#self.normr2 = numpy.sqrt(self.normr2)
@@ -65,7 +68,8 @@ class CGLS(Algorithm):
self.configured = True
def update(self):
-
+ self.update_new()
+ def update_old(self):
Ad = self.operator.direct(self.d)
#norm = (Ad*Ad).sum()
#if isinstance(norm, Iterable):
@@ -87,5 +91,44 @@ class CGLS(Algorithm):
self.normr2 = normr2_new
self.d = s + beta*self.d
+ def update_new(self):
+
+ Ad = self.operator.direct(self.d)
+ norm = Ad.squared_norm()
+ if norm == 0.:
+ print ('norm = 0, cannot update solution')
+ print ("self.d norm", self.d.squared_norm(), self.d.as_array())
+ raise StopIteration()
+ alpha = self.normr2/norm
+ if alpha == 0.:
+ print ('alpha = 0, cannot update solution')
+ raise StopIteration()
+ self.d *= alpha
+ Ad *= alpha
+ self.r -= Ad
+
+ self.x += self.d
+
+ self.operator.adjoint(self.r, out=self.s)
+ s = self.s
+
+ normr2_new = s.squared_norm()
+
+ beta = normr2_new/self.normr2
+ self.normr2 = normr2_new
+ self.d *= (beta/alpha)
+ self.d += s
+
def update_objective(self):
- self.loss.append(self.r.squared_norm())
+ a = self.r.squared_norm()
+ if a is numpy.nan:
+ raise StopIteration()
+ self.loss.append(a)
+
+# def should_stop(self):
+# if self.iteration > 0:
+# x = self.get_last_objective()
+# a = x > 0
+# return self.max_iteration_stop_cryterion() or (not a)
+# else:
+# return False
diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
index 0b6bb0b..8e77f56 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
@@ -88,7 +88,6 @@ class Norm2Sq(Function):
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)
diff --git a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
index 1db223b..8bf502a 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
@@ -18,6 +18,7 @@
# limitations under the License.
from numbers import Number
import numpy
+import warnings
class ScaledFunction(object):
@@ -89,7 +90,7 @@ class ScaledFunction(object):
'''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)
+ return self.proximal(x, tau, out=None)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
index 8cc4c71..90ef938 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
@@ -2,14 +2,17 @@ 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 VectorData
+from ccpi.framework import VectorGeometry
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
+ M_A, N_A = self.A.shape
+ self.gm_domain = VectorGeometry(N_A)
+ self.gm_range = VectorGeometry(M_A)
self.s1 = None # Largest singular value, initially unknown
super(LinearOperatorMatrix, self).__init__()
@@ -18,15 +21,13 @@ class LinearOperatorMatrix(LinearOperator):
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
@@ -34,3 +35,7 @@ class LinearOperatorMatrix(LinearOperator):
# If unknown, compute and store. If known, simply return it.
return svds(self.A,1,return_singular_vectors=False)[0]
+ def domain_geometry(self):
+ return self.gm_domain
+ def range_geometry(self):
+ return self.gm_range
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 4c53df8..16f7b86 100755
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -201,21 +201,25 @@ class TestDataContainer(unittest.TestCase):
self.assertNumpyArrayEqual(out.as_array(), ds2.as_array())
ds0 = ds
- steps.append(timer())
- ds0.add(2, out=out)
- steps.append(timer())
- print("ds0.add(2,out=out)", dt(steps), 3, ds0.as_array()[0][0][0])
- self.assertEqual(3., out.as_array()[0][0][0])
-
- dt1 = dt(steps)
- steps.append(timer())
- ds3 = ds0.add(2)
- steps.append(timer())
- print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0])
- dt2 = dt(steps)
- self.assertLess(dt1, dt2)
+ dt1 = 0
+ dt2 = 0
+ for i in range(10):
+ steps.append(timer())
+ ds0.add(2, out=out)
+ steps.append(timer())
+ print("ds0.add(2,out=out)", dt(steps), 3, ds0.as_array()[0][0][0])
+ self.assertEqual(3., out.as_array()[0][0][0])
+
+ dt1 += dt(steps)/10
+ steps.append(timer())
+ ds3 = ds0.add(2)
+ steps.append(timer())
+ print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0])
+ dt2 += dt(steps)/10
self.assertNumpyArrayEqual(out.as_array(), ds3.as_array())
+ self.assertLess(dt1, dt2)
+
def binary_subtract(self):
print("Test binary subtract")
@@ -324,16 +328,19 @@ class TestDataContainer(unittest.TestCase):
ds = DataContainer(a, False, ['X', 'Y', 'Z'])
ds1 = ds.copy()
- steps.append(timer())
- ds.divide(ds1, out=ds)
- steps.append(timer())
- t1 = dt(steps)
- print("ds.divide(ds1, out=ds)", dt(steps))
- steps.append(timer())
- ds2 = ds.divide(ds1)
- steps.append(timer())
- t2 = dt(steps)
- print("ds2 = ds.divide(ds1)", dt(steps))
+ t1 = 0
+ t2 = 0
+ for i in range(10):
+ steps.append(timer())
+ ds.divide(ds1, out=ds)
+ steps.append(timer())
+ t1 += dt(steps)/10.
+ print("ds.divide(ds1, out=ds)", dt(steps))
+ steps.append(timer())
+ ds2 = ds.divide(ds1)
+ steps.append(timer())
+ t2 += dt(steps)/10.
+ print("ds2 = ds.divide(ds1)", dt(steps))
self.assertLess(t1, t2)
self.assertEqual(ds.as_array()[0][0][0], 1.)
diff --git a/Wrappers/Python/test/test_run_test.py b/Wrappers/Python/test/test_run_test.py
index 22d4a78..81ee738 100755
--- a/Wrappers/Python/test/test_run_test.py
+++ b/Wrappers/Python/test/test_run_test.py
@@ -10,9 +10,8 @@ from ccpi.optimisation.algorithms import FISTA
#from ccpi.optimisation.algs import FBPD
from ccpi.optimisation.functions import Norm2Sq
from ccpi.optimisation.functions import ZeroFunction
+# from ccpi.optimisation.funcs import Norm1
from ccpi.optimisation.functions import L1Norm
-# This was removed
-#from ccpi.optimisation.funcs import Norm2
from ccpi.optimisation.operators import LinearOperatorMatrix
from ccpi.optimisation.operators import Identity
@@ -138,8 +137,11 @@ class TestAlgorithms(unittest.TestCase):
# A = Identity()
# Change n to equal to m.
-
- b = DataContainer(bmat)
+ vgb = VectorGeometry(m)
+ vgx = VectorGeometry(n)
+ b = vgb.allocate()
+ b.fill(bmat)
+ #b = DataContainer(bmat)
# Regularization parameter
lam = 10
@@ -149,10 +151,11 @@ class TestAlgorithms(unittest.TestCase):
g0 = ZeroFunction()
# Initial guess
- x_init = DataContainer(np.zeros((n, 1)))
+ #x_init = DataContainer(np.zeros((n, 1)))
+ x_init = vgx.allocate()
# Create 1-norm object instance
- g1 = lam * L1Norm()
+ g1 = Norm1(lam)
g1(x_init)
g1.prox(x_init, 0.02)
@@ -226,7 +229,7 @@ class TestAlgorithms(unittest.TestCase):
# Create 1-norm object instance
- g1 = lam * L1Norm()
+ g1 = Norm1(lam)
# Compare to CVXPY
@@ -293,7 +296,7 @@ class TestAlgorithms(unittest.TestCase):
# 1-norm regulariser
lam1_denoise = 1.0
- g1_denoise = lam1_denoise * L1Norm()
+ g1_denoise = Norm1(lam1_denoise)
# Initial guess
x_init_denoise = ImageData(np.zeros((N, N)))
diff --git a/Wrappers/Python/wip/compare_CGLS_algos.py b/Wrappers/Python/wip/compare_CGLS_algos.py
index 119752c..52f3f31 100644
--- a/Wrappers/Python/wip/compare_CGLS_algos.py
+++ b/Wrappers/Python/wip/compare_CGLS_algos.py
@@ -12,6 +12,8 @@ from ccpi.optimisation.algorithms import CGLS as CGLSalg
import numpy as np
import matplotlib.pyplot as plt
+from ccpi.optimisation.functions import Norm2Sq
+
# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
test_case = 1
@@ -101,25 +103,29 @@ x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt)
#plt.title('CGLS criterion')
#plt.show()
+f = Norm2Sq(Aop, b, c=1.)
+def callback(it, objective, solution):
+ print (objective, f(solution))
# 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'])
+CGLS_alg.max_iteration = 500
+CGLS_alg.update_objective_interval = 10
+CGLS_alg.run(300, callback=callback)
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.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()
+plt.figure()
+plt.semilogy(CGLS_alg.objective)
+plt.title('CGLS criterion')
+plt.show()
print(criter_CGLS)
print(CGLS_alg.objective)
diff --git a/Wrappers/Python/wip/fix_test.py b/Wrappers/Python/wip/fix_test.py
new file mode 100755
index 0000000..b1006c0
--- /dev/null
+++ b/Wrappers/Python/wip/fix_test.py
@@ -0,0 +1,208 @@
+import numpy as np
+import numpy
+from ccpi.optimisation.operators import *
+from ccpi.optimisation.algorithms import *
+from ccpi.optimisation.functions import *
+from ccpi.framework import *
+
+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 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) )
+
+opt = {'memopt': True}
+# Problem data.
+m = 500
+n = 200
+
+# if m < n then the problem is under-determined and algorithms will struggle to find a solution.
+# One approach is to add regularisation
+
+#np.random.seed(1)
+Amat = np.asarray( np.random.randn(m, n), dtype=numpy.float32)
+Amat = np.asarray( np.random.random_integers(1,10, (m, n)), dtype=numpy.float32)
+#Amat = np.asarray(np.eye(m), dtype=np.float32) * 2
+A = LinearOperatorMatrix(Amat)
+bmat = np.asarray( np.random.randn(m), dtype=numpy.float32)
+#bmat *= 0
+#bmat += 2
+print ("bmat", bmat.shape)
+print ("A", A.A)
+#bmat.shape = (bmat.shape[0], 1)
+
+# A = Identity()
+# Change n to equal to m.
+vgb = VectorGeometry(m)
+vgx = VectorGeometry(n)
+b = vgb.allocate(VectorGeometry.RANDOM_INT, dtype=numpy.float32)
+# b.fill(bmat)
+#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=1., memopt=True)
+#f = FunctionOperatorComposition(A, L2NormSquared(b=bmat))
+g0 = ZeroFunction()
+
+#f.L = 30.003
+x_init = vgx.allocate(VectorGeometry.RANDOM, dtype=numpy.float32)
+x_initcgls = x_init.copy()
+
+a = VectorData(x_init.as_array(), deep_copy=True)
+
+assert id(x_init.as_array()) != id(a.as_array())
+
+
+#f.L = LinearOperator.PowerMethod(A, 25, x_init)[0]
+#print ('f.L', f.L)
+rate = (1 / f.L) / 6
+#f.L *= 12
+
+# Initial guess
+#x_init = DataContainer(np.zeros((n, 1)))
+print ('x_init', x_init.as_array())
+print ('b', b.as_array())
+# Create 1-norm object instance
+g1_new = lam * L1Norm()
+g1 = Norm1(lam)
+
+g1 = ZeroFunction()
+#g1(x_init)
+x = g1.prox(x_init, 1/f.L )
+print ("g1.proximal ", x.as_array())
+
+x = g1.prox(x_init, 0.03 )
+print ("g1.proximal ", x.as_array())
+x = g1_new.proximal(x_init, 0.03 )
+print ("g1.proximal ", x.as_array())
+
+x1 = vgx.allocate(VectorGeometry.RANDOM, dtype=numpy.float32)
+pippo = vgx.allocate()
+
+print ("x_init", x_init.as_array())
+print ("x1", x1.as_array())
+a = x_init.subtract(x1, out=pippo)
+
+print ("pippo", pippo.as_array())
+print ("x_init", x_init.as_array())
+print ("x1", x1.as_array())
+
+y = A.direct(x_init)
+y *= 0
+A.direct(x_init, out=y)
+
+# Combine with least squares and solve using generic FISTA implementation
+#x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt)
+def callback(it, objective, solution):
+ print ("callback " , it , objective, f(solution))
+
+fa = FISTA(x_init=x_init, f=f, g=g1)
+fa.max_iteration = 1000
+fa.update_objective_interval = int( fa.max_iteration / 10 )
+fa.run(fa.max_iteration, callback = None, verbose=True)
+
+gd = GradientDescent(x_init=x_init, objective_function=f, rate = rate )
+gd.max_iteration = 5000
+gd.update_objective_interval = int( gd.max_iteration / 10 )
+gd.run(gd.max_iteration, callback = None, verbose=True)
+
+
+
+cgls = CGLS(x_init= x_initcgls, operator=A, data=b)
+cgls.max_iteration = 1000
+cgls.update_objective_interval = int( cgls.max_iteration / 10 )
+
+#cgls.should_stop = stop_criterion(cgls)
+cgls.run(cgls.max_iteration, callback = callback, verbose=True)
+
+
+
+# Print for comparison
+print("FISTA least squares plus 1-norm solution and objective value:")
+print(fa.get_output().as_array())
+print(fa.get_last_objective())
+
+print ("data ", b.as_array())
+print ('FISTA ', A.direct(fa.get_output()).as_array())
+print ('GradientDescent', A.direct(gd.get_output()).as_array())
+print ('CGLS ', A.direct(cgls.get_output()).as_array())
+
+cond = numpy.linalg.cond(A.A)
+
+print ("cond" , cond)
+
+#%%
+try:
+ import cvxpy as cp
+ # Construct the problem.
+ x = cp.Variable(n)
+ objective = cp.Minimize(cp.sum_squares(A.A*x - bmat))
+ prob = cp.Problem(objective)
+ # The optimal objective is returned by prob.solve().
+ result = prob.solve(solver = cp.MOSEK)
+
+ print ('CGLS ', cgls.get_output().as_array())
+ print ('CVX ', x.value)
+
+ print ('FISTA ', fa.get_output().as_array())
+ print ('GD ', gd.get_output().as_array())
+except ImportError as ir:
+ pass
+
+ #%%
+
+
+
+
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-29 04:55:09 +0000