precond with Tomo

14 files changed, 302 insertions, 57 deletions

diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py
index 8e55b67..21ef3f0 100644
--- a/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py
@@ -52,6 +52,11 @@ 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]

+        

         if isinstance(other, Number):

             return True   

         elif isinstance(other, list):

diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py
index d336305..0f43155 100644
--- a/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py
+++ b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py
@@ -27,6 +27,10 @@ class BlockGeometry(object):
             raise ValueError(

                     'Dimension and size do not match: expected {} got {}'

                     .format(n_elements, len(args)))

+            

+    def get_item(self, index):

+        '''returns the Geometry in the BlockGeometry located at position index'''

+        return self.geometries[index]            

 

     def allocate(self, value=0, dimension_labels=None):

         containers = [geom.allocate(value) for geom in self.geometries]

diff --git a/Wrappers/Python/build/lib/ccpi/framework/__init__.py b/Wrappers/Python/build/lib/ccpi/framework/__init__.py
index 66e2f56..229edb5 100644
--- a/Wrappers/Python/build/lib/ccpi/framework/__init__.py
+++ b/Wrappers/Python/build/lib/ccpi/framework/__init__.py
@@ -15,6 +15,7 @@ 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

diff --git a/Wrappers/Python/build/lib/ccpi/framework/framework.py b/Wrappers/Python/build/lib/ccpi/framework/framework.py
index ae9faf7..07c2ead 100644
--- a/Wrappers/Python/build/lib/ccpi/framework/framework.py
+++ b/Wrappers/Python/build/lib/ccpi/framework/framework.py
@@ -29,6 +29,7 @@ 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]
@@ -496,6 +497,7 @@ class DataContainer(object):
     ## algebra 
     def __add__(self, other, *args, **kwargs):
         out = kwargs.get('out', None)
+        
         if issubclass(type(other), DataContainer):    
             if self.check_dimensions(other):
                 out = self.as_array() + other.as_array()
@@ -601,6 +603,7 @@ class DataContainer(object):
                                deep_copy=True, 
                                dimension_labels=self.dimension_labels,
                                geometry=self.geometry)
+            
         else:
             raise TypeError('Cannot {0} DataContainer with {1}'.format("multiply" ,
                             type(other)))
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py
index 043fe38..d0e27ae 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py
@@ -8,11 +8,14 @@ Created on Mon Feb  4 16:18:06 2019
 from ccpi.optimisation.algorithms import Algorithm
 from ccpi.framework import ImageData
 import numpy as np
-#import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
 import time
 from ccpi.optimisation.operators import BlockOperator
 from ccpi.framework import BlockDataContainer
 
+
+import matplotlib.pyplot as plt
+
 class PDHG(Algorithm):
     '''Primal Dual Hybrid Gradient'''
 
@@ -69,9 +72,10 @@ class PDHG(Algorithm):
         self.xbar *= self.theta
         self.xbar += self.x
                         
-        self.x_old.fill(self.x)
-        self.y_old.fill(self.y)
-        #self.y_old = y.copy()
+#        self.x_old.fill(self.x)
+#        self.y_old.fill(self.y)
+        self.y_old = self.y.copy()
+        self.x_old = self.x.copy()
         #self.y = self.y_old
 
     def update_objective(self):
@@ -80,3 +84,72 @@ class PDHG(Algorithm):
         ])
 
 
+
+def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
+        
+    # algorithmic parameters
+    if opt is None: 
+        opt = {'tol': 1e-6, 'niter': 500, 'show_iter': 100, \
+               'memopt': False} 
+        
+    if sigma is None and tau is None:
+        raise ValueError('Need sigma*tau||K||^2<1') 
+                
+    niter = opt['niter'] if 'niter' in opt.keys() else 1000
+    tol = opt['tol'] if 'tol' in opt.keys() else 1e-4
+    memopt = opt['memopt'] if 'memopt' in opt.keys() else False  
+    show_iter = opt['show_iter'] if 'show_iter' in opt.keys() else False 
+    stop_crit = opt['stop_crit'] if 'stop_crit' in opt.keys() else False 
+
+
+    x_old = operator.domain_geometry().allocate()
+    y_old = operator.range_geometry().allocate()       
+        
+    
+    xbar = x_old
+    x_tmp = x_old
+    x = x_old
+    
+    y_tmp = y_old
+    y = y_tmp
+        
+    # relaxation parameter
+    theta = 1
+    
+    t = time.time()
+    
+    objective = []
+    
+    
+    for i in range(niter):
+        
+        # Gradient descent, Dual problem solution
+        y_tmp = y_old + sigma * operator.direct(xbar)
+        y = f.proximal_conjugate(y_tmp, sigma)
+
+        # Gradient ascent, Primal problem solution
+        x_tmp = x_old - tau * operator.adjoint(y)
+        x = g.proximal(x_tmp, tau)
+        
+        #Update
+        xbar = x + theta * (x - x_old)
+                                
+        x_old = x
+        y_old = y   
+                
+        if i%100==0:
+            
+            primal = f(operator.direct(x)) + g(x)
+            dual = -(f.convex_conjugate(y) + g(-1*operator.adjoint(y)))
+            print( i, primal, dual, primal-dual)
+            
+#            plt.imshow(x.as_array())
+#            plt.show()
+#            print(f(operator.direct(x)) + g(x), i)
+                         
+    t_end = time.time()        
+        
+    return x, t_end - t, objective
+
+
+
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py
index 443bc78..f562973 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py
@@ -28,3 +28,5 @@ from .GradientDescent import GradientDescent
 from .FISTA import FISTA
 from .FBPD import FBPD
 from .PDHG import PDHG
+from .PDHG import PDHG_old
+
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py
index 70216a9..81c16cd 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py
@@ -10,6 +10,7 @@ import numpy as np
 #from ccpi.optimisation.funcs import Function
 from ccpi.optimisation.functions import Function
 from ccpi.framework import BlockDataContainer
+from numbers import Number
 
 class BlockFunction(Function):
     '''A Block vector of Functions
@@ -52,16 +53,24 @@ class BlockFunction(Function):
     def proximal_conjugate(self, x, tau, out = None):
         '''proximal_conjugate does not take into account the BlockOperator'''
         out = [None]*self.length
-        for i in range(self.length):
-            out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau)
+        if isinstance(tau, Number):
+            for i in range(self.length):
+                out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau)
+        else:
+            for i in range(self.length):
+                out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(i))
 
         return BlockDataContainer(*out) 
     
     def proximal(self, x, tau, out = None):
         '''proximal does not take into account the BlockOperator'''
         out = [None]*self.length
-        for i in range(self.length):
-            out[i] = self.functions[i].proximal(x.get_item(i), tau)
+        if isinstance(tau, Number):
+            for i in range(self.length):
+                out[i] = self.functions[i].proximal(x.get_item(i), tau)
+        else:
+            for i in range(self.length):
+                out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i))                        
 
         return BlockDataContainer(*out)     
     
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py
index 5489d92..889d703 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py
@@ -1,12 +1,21 @@
 # -*- 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
 
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb  7 13:10:56 2019
+#   Copyright 2018-2019 Evangelos Papoutsellis and 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
 
-@author: evangelos
-"""
+#   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.
 
 import numpy
 from ccpi.optimisation.functions import Function
@@ -71,14 +80,15 @@ class L2NormSquared(Function):
             
         tmp = 0
         if self.b is not None:
-            tmp = (self.b * x).sum()
+#            tmp = (self.b * x).sum()
+            tmp = (x * self.b).sum()
             
         if out is None:
             # FIXME: this is a number
-            return (1/4) * x.squared_norm() + tmp
+            return (1./4.) * x.squared_norm() + tmp
         else:
             # FIXME: this is a DataContainer
-            out.fill((1/4) * x.squared_norm() + tmp)
+            out.fill((1./4.) * x.squared_norm() + tmp)
                     
 
     def proximal(self, x, tau, out = None):
@@ -108,7 +118,8 @@ class L2NormSquared(Function):
         
         if out is None:
             if self.b is not None:
-                return (x - tau*self.b)/(1 + tau/2) 
+                # change the order cannot add ImageData + NestedBlock
+                return (-1* tau*self.b + x)/(1 + tau/2) 
             else:
                 return x/(1 + tau/2 )
         else:
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py
index ec14b8f..e00de0c 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py
@@ -8,9 +8,9 @@ Created on Fri Mar  1 22:50:04 2019
 
 from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator
 from ccpi.optimisation.ops import PowerMethodNonsquare
-from ccpi.framework import ImageData, ImageGeometry, BlockGeometry
+from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
 import numpy 
-from ccpi.optimisation.operators import FiniteDiff
+from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
 
 #%%
 
@@ -45,7 +45,6 @@ class Gradient(LinearOperator):
         
         tmp = self.gm_range.allocate()
         
-        
         for i in range(tmp.shape[0]):
             tmp.get_item(i).fill(FiniteDiff(self.gm_domain, direction = self.ind[i], bnd_cond = self.bnd_cond).direct(x))
         return tmp    
@@ -73,6 +72,115 @@ class Gradient(LinearOperator):
     def __rmul__(self, scalar):
         return ScaledOperator(self, scalar) 
     
+    
+    def matrix(self):
+        
+        tmp = self.gm_range.allocate()
+        
+        mat = []
+        for i in range(tmp.shape[0]):
+            
+            spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond)
+            mat.append(spMat.matrix())
+    
+        return BlockDataContainer(*mat)    
+
+
+    def sum_abs_row(self):
+        
+        tmp = self.gm_range.allocate()
+        res = self.gm_domain.allocate()
+        for i in range(tmp.shape[0]):
+            spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond)
+            res += spMat.sum_abs_row()
+        return res
+    
+    def sum_abs_col(self):
+        
+        tmp = self.gm_range.allocate()
+        res = []
+        for i in range(tmp.shape[0]):
+            spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond)
+            res.append(spMat.sum_abs_col())
+        return BlockDataContainer(*res)
+   
+       
 if __name__ == '__main__':
     
-    pass
+    
+    from ccpi.optimisation.operators import Identity, BlockOperator
+    
+    M, N = 2, 3
+    ig = ImageGeometry(M, N)
+    arr = ig.allocate('random_int' )
+    
+    # check direct of Gradient and sparse matrix
+    G = Gradient(ig)
+    G_sp = G.matrix()
+    
+    res1 = G.direct(arr)
+    res1y = numpy.reshape(G_sp[0].toarray().dot(arr.as_array().flatten('F')), ig.shape, 'F')
+    
+    print(res1[0].as_array())
+    print(res1y)
+    
+    res1x = numpy.reshape(G_sp[1].toarray().dot(arr.as_array().flatten('F')), ig.shape, 'F')
+    
+    print(res1[1].as_array())
+    print(res1x)    
+    
+    #check sum abs row
+    conc_spmat = numpy.abs(numpy.concatenate( (G_sp[0].toarray(), G_sp[1].toarray() )))
+    print(numpy.reshape(conc_spmat.sum(axis=0), ig.shape, 'F'))    
+    print(G.sum_abs_row().as_array())
+    
+    print(numpy.reshape(conc_spmat.sum(axis=1), ((2,) + ig.shape), 'F'))
+    
+    print(G.sum_abs_col()[0].as_array())
+    print(G.sum_abs_col()[1].as_array())   
+    
+    # Check Blockoperator sum abs col and row
+    
+    op1 = Gradient(ig)
+    op2 = Identity(ig)
+    
+    B = BlockOperator( op1, op2)
+    
+    Brow = B.sum_abs_row()
+    Bcol = B.sum_abs_col()
+    
+    concB = numpy.concatenate( (numpy.abs(numpy.concatenate( (G_sp[0].toarray(), G_sp[1].toarray() ))), op2.matrix().toarray()))
+    
+    print(numpy.reshape(concB.sum(axis=0), ig.shape, 'F'))
+    print(Brow.as_array())
+    
+    print(numpy.reshape(concB.sum(axis=1)[0:12], ((2,) + ig.shape), 'F'))
+    print(Bcol[1].as_array())    
+    
+        
+#    print(numpy.concatene(G_sp[0].toarray()+ ))
+#    print(G_sp[1].toarray())
+#    
+#    d1 = G.sum_abs_row()
+#    print(d1.as_array())
+#    
+#    d2 = G_neum.sum_abs_col()
+##    print(d2)    
+#    
+#    
+#    ###########################################################
+    a = BlockDataContainer( BlockDataContainer(arr, arr), arr)
+    b = BlockDataContainer( BlockDataContainer(arr+5, arr+3), arr+2)
+    c = a/b
+    
+    print(c[0][0].as_array(), (arr/(arr+5)).as_array())
+    print(c[0][1].as_array(), (arr/(arr+3)).as_array())
+    print(c[1].as_array(), (arr/(arr+2)).as_array())
+    
+    
+    a1 = BlockDataContainer( arr, BlockDataContainer(arr, arr))
+#    
+#    c1 = arr + a
+#    c2 = arr + a
+#    c2 = a1 + arr
+#    
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py
index 52c7c3b..a58a296 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py
@@ -50,11 +50,11 @@ class Identity(LinearOperator):
     
     def sum_abs_row(self):
         
-        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
+        return self.gm_domain.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
  
     def sum_abs_col(self):
         
-        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F')))
+        return self.gm_domain.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F')))
             
     
 if __name__ == '__main__':
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py
index 0fb5efb..1b88cba 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -64,11 +64,15 @@ class SparseFiniteDiff():
     
     def sum_abs_row(self):
         
-        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
+        res = np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))
+        res[res==0]=1
+        return ImageData(res)
     
     def sum_abs_col(self):
         
-        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F')))
+        res = np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'C'))
+        res[res==0]=1
+        return ImageData(res)
         
 if __name__ == '__main__':
     
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
index a28c0bf..f562973 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
@@ -29,3 +29,4 @@ from .FISTA import FISTA
 from .FBPD import FBPD
 from .PDHG import PDHG
 from .PDHG import PDHG_old
+
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
index 2e0b9f4..d5c021d 100644
--- a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
+++ b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
@@ -74,11 +74,12 @@ else:
     ###########################################################################
 #%%
 
-diag_precon = False 
+diag_precon = True 
 
 if diag_precon:
     
     def tau_sigma_precond(operator):
+        
         tau = 1/operator.sum_abs_row()
         sigma = 1/ operator.sum_abs_col()
     
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
index 8bd63c2..9feb05b 100644
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py
+++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
@@ -26,10 +26,10 @@ from skimage.util import random_noise
 #%%###############################################################################
 # Create phantom for TV tomography
 
-import os
-import tomophantom
-from tomophantom import TomoP2D
-from tomophantom.supp.qualitymetrics import QualityTools
+#import os
+#import tomophantom
+#from tomophantom import TomoP2D
+#from tomophantom.supp.qualitymetrics import QualityTools
 
 #model = 1 # select a model number from the library
 #N = 150 # set dimension of the phantom
@@ -55,7 +55,7 @@ detectors = 150
 angles = np.linspace(0,np.pi,100)
 
 ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-Aop = AstraProjectorSimple(ig, ag, 'gpu')
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
 sin = Aop.direct(data)
 
 plt.imshow(sin.as_array())
@@ -95,32 +95,55 @@ g = ZeroFun()
 normK = operator.norm()
 
 ## Primal & dual stepsizes
-
-sigma = 10
-tau = 1/(sigma*normK**2)
-
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 5000
-pdhg.update_objective_interval = 250
-
-pdhg.run(5000)
-
-#%%
-sol = pdhg.get_output().as_array()
-fig = plt.figure()
-plt.subplot(1,2,1)
-plt.imshow(noisy_data.as_array())
-#plt.colorbar()
-plt.subplot(1,2,2)
-plt.imshow(sol)
-#plt.colorbar()
-plt.show()
-
+diag_precon = False 
+
+if diag_precon:
+    
+    def tau_sigma_precond(operator):
+        
+        tau = 1/operator.sum_abs_row()
+        sigma = 1/ operator.sum_abs_col()
+    
+        return tau, sigma
+
+    tau, sigma = tau_sigma_precond(operator)
+    
+else:    
+    sigma = 10
+    tau = 1/(sigma*normK**2)
+
+#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+#pdhg.max_iteration = 5000
+#pdhg.update_objective_interval = 250
+#
+#pdhg.run(5000)
+    
+opt = {'niter':2000}
+#
+res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+
+aaa = res[0].as_array()
+#    
+plt.imshow(aaa)
+plt.colorbar()
+plt.show()    
 
 #%%
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
-plt.legend()
-plt.show()
+#sol = pdhg.get_output().as_array()
+#fig = plt.figure()
+#plt.subplot(1,2,1)
+#plt.imshow(noisy_data.as_array())
+##plt.colorbar()
+#plt.subplot(1,2,2)
+#plt.imshow(sol)
+##plt.colorbar()
+#plt.show()
+#
+#
+##%%
+#plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
+#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
+#plt.legend()
+#plt.show()