Kullback leibler example

7 files changed, 285 insertions, 16 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index df53e57..f25cdbf 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -139,6 +139,8 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
         x_old = x
         y_old = y   
                 
+
+
 #        if isinstance(f, FunctionOperatorComposition):
 #        p1 = f(x) + g(x)
 #        else:
diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
new file mode 100644
index 0000000..18af154
--- /dev/null
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -0,0 +1,82 @@
+# -*- 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.
+
+import numpy
+from ccpi.optimisation.functions import Function
+from ccpi.optimisation.functions.ScaledFunction import ScaledFunction
+from ccpi.framework import DataContainer, ImageData, ImageGeometry  
+
+class KullbackLeibler(Function):
+    
+    def __init__(self,data,**kwargs):
+        
+        super(KullbackLeibler, self).__init__()
+        
+        self.b = data        
+        self.bnoise = kwargs.get('bnoise', 0)
+
+        self.sum_value = self.b + self.bnoise        
+        if  (self.sum_value.as_array()<0).any():
+            self.sum_value = numpy.inf
+                                                
+    def __call__(self, x):
+        
+        if self.sum_value==numpy.inf:
+            return numpy.inf
+        else:
+            return numpy.sum( x.as_array() - self.b.as_array() * numpy.log(self.sum_value.as_array()))
+
+        
+    def gradient(self, x):
+        
+        #TODO Division check
+        return 1 - self.b/(x + self.bnoise)
+    
+    def convex_conjugate(self, x, out=None):
+        pass
+    
+    def proximal(self, x, tau, out=None):
+        
+        z = x + tau * self.bnoise
+        return (z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()
+
+    
+    def proximal_conjugate(self, x, tau, out=None):
+        pass
+        
+        
+    
+
+if __name__ == '__main__':
+    
+    N, M = 2,3
+    ig  = ImageGeometry(N, M)
+    data = ImageData(numpy.random.randint(-10, 100, size=(M, N)))
+    x = ImageData(numpy.random.randint(-10, 100, size=(M, N)))
+    
+    bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N)))
+    
+    f = KullbackLeibler(data, bnoise=bnoise)
+    print(f.sum_value)
+    
+    print(f(x))
+
+
+    
+        
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/functions/__init__.py b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
index 9dbb505..65e8848 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
@@ -9,4 +9,5 @@ from .BlockFunction import BlockFunction
 from .FunctionOperatorComposition import FunctionOperatorComposition
 from .MixedL21Norm import MixedL21Norm
 from .IndicatorBox import IndicatorBox
+from .KullbackLeibler import KullbackLeibler
 from .Norm2Sq import Norm2sq
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 484dc61..6c080bb 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -235,6 +235,7 @@ class BlockOperator(Operator):
             tmp = sum(res)
             return ImageData(tmp)
         else:
+            
             return BlockDataContainer(*res)
         
     def sum_abs_col(self):
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
index d54db9b..5e318ff 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -70,7 +70,7 @@ class SparseFiniteDiff():
     
     def sum_abs_col(self):
         
-        res = 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, 'F') )
         res[res==0]=1
         return ImageData(res)
         
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
index e9a85cc..dea8e5c 100644
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
+++ b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py
@@ -13,7 +13,7 @@ from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, Acquisi
 import numpy as np                           
 import matplotlib.pyplot as plt
 
-from ccpi.optimisation.algorithms import PDHG
+from ccpi.optimisation.algorithms import PDHG, PDHG_old
 
 from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
 from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \
@@ -107,26 +107,41 @@ normK = operator.norm()
 
 ## Primal & dual stepsizes
 
-sigma = 10
+sigma = 1
 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 = 20
+#sigma = 1/normK
+#tau = 1/normK
 
-pdhg.run(5000)
+opt = {'niter':2000}
 
-#%%
-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()
+res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+ 
+plt.figure(figsize=(5,5))
+plt.imshow(res.as_array())
+plt.colorbar()
 plt.show()
 
+#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 = 20
+#
+#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()
+
 
 #%%
 plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py
new file mode 100644
index 0000000..9fad6f8
--- /dev/null
+++ b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py
@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Feb 22 14:53:03 2019
+
+@author: evangelos
+"""
+
+from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer
+
+import numpy as np                           
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG, PDHG_old
+
+from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
+from ccpi.optimisation.functions import ZeroFun, KullbackLeibler, \
+                      MixedL21Norm, FunctionOperatorComposition, BlockFunction
+                 
+
+from skimage.util import random_noise
+
+
+
+# ############################################################################
+# Create phantom for TV denoising
+
+N = 100
+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
+
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+ag = ig
+
+# Create noisy data. Add Gaussian noise
+n1 = random_noise(data, mode = 'poisson')
+noisy_data = ImageData(n1)
+
+plt.imshow(noisy_data.as_array())
+plt.colorbar()
+plt.show()
+
+#%%
+
+# Regularisation Parameter
+alpha = 10
+
+#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
+method = '1'
+if method == '0':
+
+    # Create operators
+    op1 = Gradient(ig)
+    op2 = Identity(ig, ag)
+
+    # Form Composite Operator
+    operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+    #### Create functions
+#    f = FunctionComposition_new(operator, mixed_L12Norm(alpha), \
+#                                    L2NormSq(0.5, b = noisy_data) )    
+    
+    f1 = alpha * MixedL21Norm()
+    f2 = KullbackLeibler(b = noisy_data)
+    
+    f = BlockFunction(f1, f2 )                                        
+    g = ZeroFun()
+    
+else:
+    
+    ###########################################################################
+    #         No Composite #
+    ###########################################################################
+    operator = Gradient(ig)
+    f = alpha *  FunctionOperatorComposition(operator, MixedL21Norm())
+    g = KullbackLeibler(noisy_data)
+    ###########################################################################
+#%%
+    
+# Compute operator Norm
+normK = operator.norm()
+print ("normK", normK)
+# Primal & dual stepsizes
+#sigma = 1
+#tau = 1/(sigma*normK**2)
+
+sigma = 1/normK
+tau = 1/normK
+
+opt = {'niter':2000}
+
+res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+ 
+plt.figure(figsize=(5,5))
+plt.imshow(res.as_array())
+plt.colorbar()
+plt.show()
+
+#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+#pdhg.max_iteration = 2000
+#pdhg.update_objective_interval = 10
+#
+#pdhg.run(2000)
+
+    
+
+#sol = pdhg.get_output().as_array()
+##sol = result.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[int(N/2),:], label = 'GTruth')
+#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
+#plt.legend()
+#plt.show()
+
+
+#%% Compare with cvx
+
+#try_cvx = input("Do you want CVX comparison (0/1)")
+#
+#if try_cvx=='0':
+#
+#    from cvxpy import *
+#    import sys
+#    sys.path.insert(0,'/Users/evangelos/Desktop/Projects/CCPi/CCPi-Framework/Wrappers/Python/ccpi/optimisation/cvx_scripts')
+#    from cvx_functions import TV_cvx
+#
+#    u = Variable((N, N))
+#    fidelity = pnorm( u - noisy_data.as_array(),1)
+#    regulariser = alpha * TV_cvx(u)
+#    solver = MOSEK
+#    obj =  Minimize( regulariser +  fidelity)
+#    constraints = []
+#    prob = Problem(obj, constraints)
+#
+#    # Choose solver (SCS is fast but less accurate than MOSEK)
+#    result = prob.solve(verbose = True, solver = solver)
+#
+#    print('Objective value is {} '.format(obj.value))
+#
+#    diff_pdhg_cvx = np.abs(u.value - res.as_array())
+#    plt.imshow(diff_pdhg_cvx)
+#    plt.colorbar()
+#    plt.title('|CVX-PDHG|')        
+#    plt.show()
+#
+#    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
+#    plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'PDHG')
+#    plt.legend()
+#    plt.show()
+#
+#else:
+#    print('No CVX solution available')
+
+
+
+