changes for memopt option

6 files changed, 93 insertions, 50 deletions

diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
index 69b6931..8934f49 100755
--- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
@@ -53,9 +53,9 @@ class BlockDataContainer(object):
     def is_compatible(self, other):

         '''basic check if the size of the 2 objects fit'''

         

-        for i in range(len(self.containers)):

-            if type(self.containers[i])==type(self):

-                self = self.containers[i]

+#        for i in range(len(self.containers)):

+#            if type(self.containers[i])==type(self):

+#                self = self.containers[i]

         

         if isinstance(other, Number):

             return True   

@@ -341,3 +341,32 @@ class BlockDataContainer(object):
         '''Inline truedivision'''

         return self.__idiv__(other)

 

+if __name__ == '__main__':

+    

+    M, N, K = 2,3,5

+    from ccpi.framework import ImageGeometry, BlockGeometry, BlockDataContainer

+    

+    ig = ImageGeometry(N, M)

+    u = ig.allocate('random_int')

+    

+    BG = BlockGeometry(ig, ig)

+    U = BG.allocate('random_int')

+    

+    U_nested = BlockDataContainer(BlockDataContainer(u, u), u)

+    

+    

+    res1 = U + u

+    res2 = U_nested + u

+#    res2 = u + U

+

+    

+    

+    

+    

+    

+    

+    

+    

+    

+    

+    
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index 07c2ead..e03a29c 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -1109,6 +1109,7 @@ class AcquisitionData(DataContainer):
                 
             
 class DataProcessor(object):
+    
     '''Defines a generic DataContainer processor
     
     accepts DataContainer as inputs and 
@@ -1154,6 +1155,7 @@ class DataProcessor(object):
         raise NotImplementedError('Implement basic checks for input DataContainer')
         
     def get_output(self, out=None):
+        
         for k,v in self.__dict__.items():
             if v is None and k != 'output':
                 raise ValueError('Key {0} is None'.format(k))
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index 2a69857..5323e76 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -161,18 +161,20 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
         
         if not memopt:
             
-            y_tmp = y_old + sigma * operator.direct(xbar)
-            y = f.proximal_conjugate(y_tmp, sigma)
+            y_old += sigma * operator.direct(xbar)            
+            y = f.proximal_conjugate(y_old, sigma)
             
-            x_tmp = x_old - tau * operator.adjoint(y)
-            x = g.proximal(x_tmp, tau)
-            
-            xbar = x + theta * (x - x_old)
+            x_old -= tau*operator.adjoint(y)
+            x = g.proximal(x_old, tau)
+             
+            xbar.fill(x)
+            xbar -= x_old
+            xbar *= theta
+            xbar += x
             
-            x_old = x
-            y_old = y            
-        
-        
+            x_old.fill(x)
+            y_old.fill(y)            
+                               
         else:
             
             operator.direct(xbar, out = y_tmp)             
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
index 0ce7d8a..0c658a4 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -77,9 +77,7 @@ class MixedL21Norm(Function):
         pass
     
     def proximal_conjugate(self, x, tau, out=None): 
-        
-        
-        
+
         if self.SymTensor:
             
             param = [1]*x.shape[0]
@@ -91,38 +89,22 @@ class MixedL21Norm(Function):
             return res
         
         else:
-#             pass
 
-    
-# #                tmp2 = np.sqrt(x.as_array()[0]**2 +  x.as_array()[1]**2 +  2*x.as_array()[2]**2)/self.alpha
-# #                res = x.divide(ImageData(tmp2).maximum(1.0))                                
+                              
          if out is None:
             
-            tmp = [ el*el for el in x]
-            res = (sum(tmp).sqrt()).maximum(1.0) 
-            frac = [x[i]/res for i in range(x.shape[0])]
-            res = BlockDataContainer(*frac)
-                                                   
-            return res
+            tmp = [ el*el for el in x.containers]
+            res = (sum(tmp).sqrt()).maximum(1.0)                                                                
+            return x/res
         
          else:
-             
-            tmp = [ el*el for el in x]
-            res = (sum(tmp).sqrt()).maximum(1.0)             
-            out.fill(x/res)             
-             
-             
-             
-        #     tmp = [ el*el for el in x]
-        #     res = (sum(tmp).sqrt()).maximum(1.0) 
-        #     #frac = [x[i]/res for i in range(x.shape[0])]
-        #     for i in range(x.shape[0]):
-        #         a = out.get_item(i)
-        #         b = x.get_item(i)
-        #         b /= res
-        #         a.fill( b )
+                         
+            res = (sum(x**2).sqrt()).maximum(1.0)             
+            out.fill(x/res)
+
+            
+
 
-    
     def __rmul__(self, scalar):
         return ScaledFunction(self, scalar) 
 
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py
index 598acb0..22fee90 100755
--- a/Wrappers/Python/wip/pdhg_TV_denoising.py
+++ b/Wrappers/Python/wip/pdhg_TV_denoising.py
@@ -23,11 +23,13 @@ from timeit import default_timer as timer
 def dt(steps):
     return steps[-1] - steps[-2]
 
+#%%
+
 
 # ############################################################################
 # Create phantom for TV denoising
 
-N = 100
+N = 200
 
 data = np.zeros((N,N))
 data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
@@ -62,8 +64,8 @@ if method == '0':
 
     #### Create functions
       
-    f1 = MixedL21Norm()
-    f2 = L2NormSquared(b = noisy_data)    
+    f1 = alpha * MixedL21Norm()
+    f2 = 0.5 * L2NormSquared(b = noisy_data)    
     f = BlockFunction(f1, f2)  
                                       
     g = ZeroFun()
@@ -88,15 +90,18 @@ print ("normK", normK)
 sigma = 1
 tau = 1/(sigma*normK**2)
 
-<<<<<<< HEAD
 opt = {'niter':100}
 opt1 = {'niter':100, 'memopt': True}
 
+t1 = timer()
 res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+print(timer()-t1)
 
 print("with memopt \n")
 
+t2 = timer()
 res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
+print(timer()-t2)
 
 plt.figure(figsize=(5,5))
 plt.imshow(res.as_array())
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
index f06f166..159f2ea 100644
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py
+++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
@@ -21,6 +21,7 @@ from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \
 
 from ccpi.astra.ops import AstraProjectorSimple
 from skimage.util import random_noise
+from timeit import default_timer as timer
 
 
 #%%###############################################################################
@@ -118,15 +119,37 @@ else:
 #
 #pdhg.run(5000)
     
-opt = {'niter':2000}
-#
+opt = {'niter':300}
+opt1 = {'niter':300, 'memopt': True}
+
+
+t1 = timer()
 res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
 
+print(timer()-t1)
 plt.figure(figsize=(5,5))
 plt.imshow(res.as_array())
 plt.colorbar()
-plt.show() 
-
+plt.show()
+ 
+#%%
+print("with memopt \n")
+#
+t2 = timer()
+res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
+#print(timer()-t2)
+#
+#
+plt.figure(figsize=(5,5))
+plt.imshow(res1.as_array())
+plt.colorbar()
+plt.show()
+#
+#%%
+plt.figure(figsize=(5,5))
+plt.imshow(np.abs(res1.as_array()-res.as_array()))
+plt.colorbar()
+plt.show()
 #%%
 #sol = pdhg.get_output().as_array()
 #fig = plt.figure()