summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers')
-rw-r--r--Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py5
-rw-r--r--Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py4
-rw-r--r--Wrappers/Python/build/lib/ccpi/framework/__init__.py1
-rw-r--r--Wrappers/Python/build/lib/ccpi/framework/framework.py3
-rw-r--r--Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py81
-rw-r--r--Wrappers/Python/build/lib/ccpi/optimisation/algorithms/__init__.py2
-rw-r--r--Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py17
-rw-r--r--Wrappers/Python/build/lib/ccpi/optimisation/functions/L2NormSquared.py31
-rw-r--r--Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py116
-rw-r--r--Wrappers/Python/build/lib/ccpi/optimisation/operators/IdentityOperator.py4
-rw-r--r--Wrappers/Python/build/lib/ccpi/optimisation/operators/SparseFiniteDiff.py8
-rw-r--r--Wrappers/Python/ccpi/optimisation/algorithms/__init__.py1
-rw-r--r--Wrappers/Python/wip/pdhg_TV_denoising_precond.py3
-rw-r--r--Wrappers/Python/wip/pdhg_TV_tomography2D.py83
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()
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 23:51:06 +0000