try precond

3 files changed, 479 insertions, 0 deletions

diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_DiagPrecond.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_DiagPrecond.py
new file mode 100644
index 0000000..d65478c
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_DiagPrecond.py
@@ -0,0 +1,208 @@
+# -*- 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 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
+
+#   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.
+
+
+""" 
+Total Variation Denoising using PDHG algorithm
+
+Problem:     min_x \alpha * ||\nabla x||_{1} + || x - g ||_{2}^{2}
+
+             \nabla: Gradient operator 
+             g: Noisy Data with Gaussian Noise
+             \alpha: Regularization parameter
+
+"""
+
+from ccpi.framework import ImageData, ImageGeometry
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
+from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
+                      MixedL21Norm, BlockFunction
+
+
+# Create phantom for TV Gaussian denoising
+N = 400
+
+data = np.zeros((N,N))
+data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+data = ImageData(data)
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+ag = ig
+
+# Create noisy data. Add Gaussian noise
+np.random.seed(10)
+noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.05, size=ig.shape) )
+
+# Regularisation Parameter
+alpha = 2
+
+method = '1'
+
+if method == '0':
+
+    # Create operators
+    op1 = Gradient(ig)
+    op2 = Identity(ig, ag)
+
+    # Create BlockOperator
+    operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+    # Create functions
+      
+    f1 = alpha * MixedL21Norm()
+    f2 = 0.5 * L2NormSquared(b = noisy_data)    
+    f = BlockFunction(f1, f2)  
+                                      
+    g = ZeroFunction()
+    
+else:
+    
+    # Without the "Block Framework"
+    operator =  Gradient(ig)
+    f =  alpha * MixedL21Norm()
+    g =  0.5 * L2NormSquared(b = noisy_data)
+        
+    
+diag_precon =  False
+
+
+if diag_precon:
+    
+    def tau_sigma_precond(operator):
+        
+        tau = 1/operator.sum_abs_col()
+        sigma = 1/operator.sum_abs_row()
+        
+        sigma[0].as_array()[sigma[0].as_array()==np.inf]=0
+        sigma[1].as_array()[sigma[1].as_array()==np.inf]=0
+        
+        return tau, sigma
+
+    tau, sigma = tau_sigma_precond(operator)
+             
+else:
+    # Compute operator Norm
+    normK = operator.norm()
+    
+    # Primal & dual stepsizes
+    sigma = 1
+    tau = 1/(sigma*normK**2)
+
+
+# Setup and run the PDHG algorithm
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
+pdhg.max_iteration = 3000
+pdhg.update_objective_interval = 200
+pdhg.run(3000, verbose=False)
+
+#%%
+plt.figure(figsize=(15,15))
+plt.subplot(3,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(3,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(3,1,3)
+plt.imshow(pdhg.get_output().as_array())
+plt.title('TV Reconstruction')
+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), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
+
+
+#%% Check with CVX solution
+
+from ccpi.optimisation.operators import SparseFiniteDiff
+
+try:
+    from cvxpy import *
+    cvx_not_installable = True
+except ImportError:
+    cvx_not_installable = False
+
+
+if cvx_not_installable:
+
+    ##Construct problem    
+    u = Variable(ig.shape)
+    
+    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+    
+    # Define Total Variation as a regulariser
+    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+    fidelity = 0.5 * sum_squares(u - noisy_data.as_array())
+    
+    # choose solver
+    if 'MOSEK' in installed_solvers():
+        solver = MOSEK
+    else:
+        solver = SCS  
+        
+    obj =  Minimize( regulariser +  fidelity)
+    prob = Problem(obj)
+    result = prob.solve(verbose = True, solver = MOSEK)
+    
+    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
+        
+    plt.figure(figsize=(15,15))
+    plt.subplot(3,1,1)
+    plt.imshow(pdhg.get_output().as_array())
+    plt.title('PDHG solution')
+    plt.colorbar()
+    plt.subplot(3,1,2)
+    plt.imshow(u.value)
+    plt.title('CVX solution')
+    plt.colorbar()
+    plt.subplot(3,1,3)
+    plt.imshow(diff_cvx)
+    plt.title('Difference')
+    plt.colorbar()
+    plt.show()    
+    
+    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
+    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
+    plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'Truth')
+    
+    plt.legend()
+    plt.title('Middle Line Profiles')
+    plt.show()
+            
+    print('Primal Objective (CVX) {} '.format(obj.value))
+    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
+
+
+
+
+
diff --git a/Wrappers/Python/wip/Demos/check_blockOperator_sum_row_cols.py b/Wrappers/Python/wip/Demos/check_blockOperator_sum_row_cols.py
new file mode 100644
index 0000000..bdb2c38
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/check_blockOperator_sum_row_cols.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May  3 13:10:09 2019
+
+@author: evangelos
+"""
+
+from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff, BlockOperator, Gradient
+from ccpi.framework import ImageGeometry, AcquisitionGeometry, BlockDataContainer, ImageData
+from ccpi.astra.ops import AstraProjectorSimple
+
+from scipy import sparse
+import numpy as np
+
+N = 3
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+u = ig.allocate('random_int')
+
+# Compare FiniteDiff with SparseFiniteDiff
+
+DY = FiniteDiff(ig, direction = 0, bnd_cond = 'Neumann')
+DX = FiniteDiff(ig, direction = 1, bnd_cond = 'Neumann')
+
+DXu = DX.direct(u)
+DYu = DY.direct(u)
+
+DX_sparse = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+DY_sparse = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+
+DXu_sparse = DX_sparse.direct(u)
+DYu_sparse = DY_sparse.direct(u)
+
+#np.testing.assert_array_almost_equal(DYu.as_array(), DYu_sparse.as_array(), decimal=4)
+#np.testing.assert_array_almost_equal(DXu.as_array(), DXu_sparse.as_array(), decimal=4)
+
+#%%  Tau/ Sigma
+
+A1 = DY_sparse.matrix()
+A2 = DX_sparse.matrix()
+A3 = sparse.eye(np.prod(ig.shape))
+
+sum_rows1 = np.array(np.sum(abs(A1), axis=1))
+sum_rows2 = np.array(np.sum(abs(A2), axis=1))
+sum_rows3 = np.array(np.sum(abs(A3), axis=1))
+
+sum_cols1 = np.array(np.sum(abs(A1), axis=0))
+sum_cols2 = np.array(np.sum(abs(A2), axis=0))
+sum_cols3 = np.array(np.sum(abs(A2), axis=0))
+
+# Check if Grad sum row/cols is OK
+Grad = Gradient(ig)
+
+Sum_Block_row = Grad.sum_abs_row()
+Sum_Block_col = Grad.sum_abs_col()
+
+tmp1 = BlockDataContainer( ImageData(np.reshape(sum_rows1, ig.shape, order='F')),\
+                           ImageData(np.reshape(sum_rows2, ig.shape, order='F')))
+
+
+#np.testing.assert_array_almost_equal(tmp1[0].as_array(), Sum_Block_row[0].as_array(), decimal=4)
+#np.testing.assert_array_almost_equal(tmp1[1].as_array(), Sum_Block_row[1].as_array(), decimal=4)
+
+tmp2 = ImageData(np.reshape(sum_cols1 + sum_cols2, ig.shape, order='F'))
+
+#np.testing.assert_array_almost_equal(tmp2.as_array(), Sum_Block_col.as_array(), decimal=4)
+
+
+#%% BlockOperator with Gradient, Identity
+
+Id = Identity(ig)
+Block_GrId = BlockOperator(Grad, Id, shape=(2,1))
+
+
+Sum_Block_GrId_row = Block_GrId.sum_abs_row()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Wrappers/Python/wip/Demos/check_precond.py b/Wrappers/Python/wip/Demos/check_precond.py
new file mode 100644
index 0000000..8cf95fa
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/check_precond.py
@@ -0,0 +1,182 @@
+# -*- 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 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
+
+#   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 ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import BlockOperator, Gradient
+from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
+                      MixedL21Norm, BlockFunction
+
+from ccpi.astra.ops import AstraProjectorSimple
+
+# Create phantom for TV 2D tomography 
+N = 75
+x = np.zeros((N,N))
+x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+
+data = ImageData(x)
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+
+detectors = N
+angles = np.linspace(0, np.pi, N, dtype=np.float32)
+
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
+sin = Aop.direct(data)
+
+# Create noisy data
+np.random.seed(10)
+n1 = np.random.random(sin.shape)
+noisy_data = sin + ImageData(5*n1)
+
+#%%
+
+# Regularisation Parameter
+alpha = 50
+
+# Create operators
+op1 = Gradient(ig)
+op2 = Aop
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+
+
+# Create functions
+      
+f1 = alpha * MixedL21Norm()
+f2 = L2NormSquared(b=noisy_data)    
+f = BlockFunction(f1, f2)  
+                                      
+g = ZeroFunction()
+    
+diag_precon =  True
+
+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:
+    # Compute operator Norm
+    normK = operator.norm()
+    # Primal & dual stepsizes
+    sigma = 10
+    tau = 1/(sigma*normK**2)
+
+
+# Setup and run the PDHG algorithm
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
+pdhg.max_iteration = 1000
+pdhg.update_objective_interval = 200
+pdhg.run(1000)
+
+#%% Check with CVX solution
+
+from ccpi.optimisation.operators import SparseFiniteDiff
+import astra
+import numpy
+
+try:
+    from cvxpy import *
+    cvx_not_installable = True
+except ImportError:
+    cvx_not_installable = False
+
+
+if cvx_not_installable:
+    
+    ##Construct problem    
+    u = Variable(N*N)
+    
+    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+
+    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+    
+    # create matrix representation for Astra operator
+
+    vol_geom = astra.create_vol_geom(N, N)
+    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
+
+    proj_id = astra.create_projector('line', proj_geom, vol_geom)
+
+    matrix_id = astra.projector.matrix(proj_id)
+
+    ProjMat = astra.matrix.get(matrix_id)
+    
+    fidelity = sum_squares( ProjMat * u - noisy_data.as_array().ravel()) 
+    #constraints = [q>=fidelity]
+#    constraints = [u>=0]
+        
+    solver = MOSEK
+    obj =  Minimize( regulariser +  fidelity)
+    prob = Problem(obj)
+    result = prob.solve(verbose = True, solver = solver)    
+
+
+#%%
+    
+plt.figure(figsize=(15,15))
+plt.subplot(2,2,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+
+plt.subplot(2,2,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+
+plt.subplot(2,2,3)
+plt.imshow(pdhg.get_output().as_array())
+plt.title('PDHG Reconstruction')
+
+plt.subplot(2,2,4)
+plt.imshow(np.reshape(u.value, ig.shape))
+plt.title('CVX Reconstruction')
+
+plt.show()
+
+#%%
+plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
+plt.plot(np.linspace(0,N,N), np.reshape(u.value, ig.shape)[int(N/2),:], label = 'CVX')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
+
+
+
+
+
+
+
+