summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xWrappers/Python/ccpi/framework/Vector.py96
-rwxr-xr-xWrappers/Python/ccpi/framework/VectorData.py59
-rwxr-xr-xWrappers/Python/ccpi/framework/VectorGeometry.py69
-rwxr-xr-xWrappers/Python/ccpi/framework/__init__.py59
-rwxr-xr-xWrappers/Python/ccpi/framework/framework.py2
-rwxr-xr-xWrappers/Python/ccpi/optimisation/algorithms/CGLS.py48
-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_run_test.py3
-rw-r--r--Wrappers/Python/wip/compare_CGLS_algos.py28
-rwxr-xr-xWrappers/Python/wip/fix_test.py167
12 files changed, 502 insertions, 50 deletions
diff --git a/Wrappers/Python/ccpi/framework/Vector.py b/Wrappers/Python/ccpi/framework/Vector.py
new file mode 100755
index 0000000..b9ce486
--- /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.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/VectorData.py b/Wrappers/Python/ccpi/framework/VectorData.py
new file mode 100755
index 0000000..47ac0dd
--- /dev/null
+++ b/Wrappers/Python/ccpi/framework/VectorData.py
@@ -0,0 +1,59 @@
+# -*- 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
+from ccpi.framework import VectorGeometry
+
+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.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)
diff --git a/Wrappers/Python/ccpi/framework/VectorGeometry.py b/Wrappers/Python/ccpi/framework/VectorGeometry.py
new file mode 100755
index 0000000..255d2a0
--- /dev/null
+++ b/Wrappers/Python/ccpi/framework/VectorGeometry.py
@@ -0,0 +1,69 @@
+# -*- 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.VectorData import VectorData
+
+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 \ No newline at end of file
diff --git a/Wrappers/Python/ccpi/framework/__init__.py b/Wrappers/Python/ccpi/framework/__init__.py
index 8926897..66eb94d 100755
--- a/Wrappers/Python/ccpi/framework/__init__.py
+++ b/Wrappers/Python/ccpi/framework/__init__.py
@@ -1,28 +1,31 @@
-# -*- 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 .VectorGeometry import VectorGeometry
+#from .VectorData import VectorData
+#from .pippo import pippo, pupi
+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..8474d89 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):
@@ -56,6 +57,8 @@ class CGLS(Algorithm):
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,45 @@ 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 ('cannot update solution')
+ raise StopIteration()
+ alpha = self.normr2/norm
+ if alpha == 0.:
+ print ('cannot update solution')
+ raise StopIteration()
+ self.d *= alpha
+ Ad *= alpha
+ self.r -= Ad
+ if numpy.isnan(self.r.as_array()).any():
+ print ("some nan")
+ raise StopIteration()
+ 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_run_test.py b/Wrappers/Python/test/test_run_test.py
index 22d4a78..4f66da1 100755
--- a/Wrappers/Python/test/test_run_test.py
+++ b/Wrappers/Python/test/test_run_test.py
@@ -13,6 +13,7 @@ from ccpi.optimisation.functions import ZeroFunction
from ccpi.optimisation.functions import L1Norm
# This was removed
#from ccpi.optimisation.funcs import Norm2
+#from ccpi.optimisation.funcs import Norm1
from ccpi.optimisation.operators import LinearOperatorMatrix
from ccpi.optimisation.operators import Identity
@@ -82,6 +83,7 @@ class TestAlgorithms(unittest.TestCase):
opt = {'memopt': True}
# Create object instances with the test data A and b.
f = Norm2Sq(A, b, c=0.5, memopt=True)
+ f.L = LinearOperator.PowerMethod(A, 10)
g0 = ZeroFunction()
# Initial guess
@@ -124,6 +126,7 @@ class TestAlgorithms(unittest.TestCase):
self.assertTrue(cvx_not_installable)
def test_FISTA_Norm1_cvx(self):
+ print ("test_FISTA_Norm1_cvx")
if not cvx_not_installable:
try:
opt = {'memopt': True}
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..316606e
--- /dev/null
+++ b/Wrappers/Python/wip/fix_test.py
@@ -0,0 +1,167 @@
+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 = 4
+n = 10
+np.random.seed(1)
+Amat = np.asarray( np.random.randn(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(2, 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())
+
+
+# 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 (objective, f(solution))
+
+fa = FISTA(x_init=x_init, f=f, g=g1)
+fa.max_iteration = 100
+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 = 100
+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 = 2
+
+#cgls.should_stop = stop_criterion(cgls)
+cgls.run(10, 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())
generated by cgit v1.2.3 (git 2.25.1) at 2025-11-04 02:32:17 +0000