Merge remote-tracking branch 'origin/fix_test' into update_fix_test

Conflicts:
	Wrappers/Python/ccpi/framework/__init__.py
	Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
	Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
	Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
	Wrappers/Python/test/test_run_test.py

11 files changed, 403 insertions, 50 deletions

diff --git a/Wrappers/Python/ccpi/framework/VectorData.py b/Wrappers/Python/ccpi/framework/VectorData.py
new file mode 100755
index 0000000..d0fed10
--- /dev/null
+++ b/Wrappers/Python/ccpi/framework/VectorData.py
@@ -0,0 +1,58 @@
+# -*- 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, 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..072ac84 100755
--- a/Wrappers/Python/ccpi/framework/__init__.py
+++ b/Wrappers/Python/ccpi/framework/__init__.py
@@ -1,28 +1,29 @@
-# -*- 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
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())