summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py35
-rw-r--r--Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py59
-rwxr-xr-xWrappers/Python/wip/pdhg_TV_denoising.py69
-rw-r--r--Wrappers/Python/wip/pdhg_tv_denoising_poisson.py66
4 files changed, 178 insertions, 51 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index a165e55..bc080f8 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -166,15 +166,15 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
y_old.fill(y)
- if i%100==0:
-
- p1 = f(operator.direct(x)) + g(x)
- d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) )
- primal.append(p1)
- dual.append(d1)
- pdgap.append(p1-d1)
- print(p1, d1, p1-d1)
-
+# if i%10==0:
+##
+# p1 = f(operator.direct(x)) + g(x)
+## print(p1)
+# d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) )
+# primal.append(p1)
+# dual.append(d1)
+# pdgap.append(p1-d1)
+# print(p1, d1, p1-d1)
else:
@@ -196,14 +196,15 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
x_old.fill(x)
y_old.fill(y)
- if i%100==0:
-
- p1 = f(operator.direct(x)) + g(x)
- d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) )
- primal.append(p1)
- dual.append(d1)
- pdgap.append(p1-d1)
- print(p1, d1, p1-d1)
+# if i%10==0:
+##
+# p1 = f(operator.direct(x)) + g(x)
+## print(p1)
+# d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) )
+# primal.append(p1)
+# dual.append(d1)
+# pdgap.append(p1-d1)
+# print(p1, d1, p1-d1)
t_end = time.time()
diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
index 40dddd7..7889cad 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -39,19 +39,14 @@ class KullbackLeibler(Function):
def __call__(self, x):
# TODO check
-
- self.sum_value = x + self.bnoise
- if (self.sum_value.as_array()<0).any():
- self.sum_value = numpy.inf
-
- if self.sum_value==numpy.inf:
- return numpy.inf
- else:
- tmp = self.sum_value.as_array()
- return (x - self.b * ImageData( numpy.log(tmp))).sum()
-
-# return numpy.sum( x.as_array() - self.b.as_array() * numpy.log(self.sum_value.as_array()))
+
+ tmp = x + self.bnoise
+ ind = tmp.as_array()>0
+ res = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind])
+
+ return sum(res)
+
def gradient(self, x, out=None):
@@ -64,10 +59,12 @@ class KullbackLeibler(Function):
def convex_conjugate(self, x):
- tmp = self.b.as_array()/( 1 - x.as_array() )
+ tmp = self.b/( 1 - x )
+ ind = tmp.as_array()>0
+
+ sel
return (self.b * ( ImageData( numpy.log(tmp) ) - 1 ) - self.bnoise * (x - 1)).sum()
-# return self.b * ( ImageData(numpy.log(self.b/(1-x)) - 1 )) - self.bnoise * (x - 1)
def proximal(self, x, tau, out=None):
@@ -105,19 +102,39 @@ class KullbackLeibler(Function):
if __name__ == '__main__':
+
+ from ccpi.framework import ImageGeometry
+ import numpy
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)))
+ data = ImageData(numpy.random.randint(-10, 10, size=(M, N)))
+ x = ImageData(numpy.random.randint(-10, 10, size=(M, N)))
- f = KullbackLeibler(data, bnoise=bnoise)
- print(f.sum_value)
+ bnoise = ImageData(numpy.random.randint(-10, 10, size=(M, N)))
- print(f(x))
+ f = KullbackLeibler(data)
+ print(f(x))
+
+# numpy.random.seed(10)
+#
+#
+# x = numpy.random.randint(-10, 10, size = (2,3))
+# b = numpy.random.randint(1, 10, size = (2,3))
+#
+# ind1 = x>0
+#
+# res = x[ind1] - b * numpy.log(x[ind1])
+#
+## ind = x>0
+#
+## y = x[ind]
+#
+#
+#
+#
+#
\ No newline at end of file
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py
index e6d38e9..0e7f628 100755
--- a/Wrappers/Python/wip/pdhg_TV_denoising.py
+++ b/Wrappers/Python/wip/pdhg_TV_denoising.py
@@ -8,7 +8,8 @@ Created on Fri Feb 22 14:53:03 2019
from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer
-import numpy as np
+import numpy as np
+import numpy
import matplotlib.pyplot as plt
from ccpi.optimisation.algorithms import PDHG, PDHG_old
@@ -89,13 +90,12 @@ t1 = timer()
res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt)
t2 = timer()
-print(" Run memopt")
t3 = timer()
res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1)
t4 = timer()
-
-#%%
+#
+##%%
plt.figure(figsize=(15,15))
plt.subplot(3,1,1)
plt.imshow(res.as_array())
@@ -115,10 +115,67 @@ plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt')
plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt')
plt.legend()
plt.show()
-
+#
print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3))
diff = (res1 - res).abs().as_array().max()
-
+#
print(" Max of abs difference is {}".format(diff))
+#%% 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 = solver)
+
+ diff_cvx = numpy.abs( res.as_array() - u.value )
+
+ # Show result
+ plt.figure(figsize=(15,15))
+ plt.subplot(3,1,1)
+ plt.imshow(res.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()
+
+ print('Primal Objective (CVX) {} '.format(obj.value))
+ print('Primal Objective (PDHG) {} '.format(primal[-1]))
+
+
+
+
diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py
index fbe0d9b..16e6b43 100644
--- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py
+++ b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py
@@ -6,7 +6,8 @@ Created on Fri Feb 22 14:53:03 2019
@author: evangelos
"""
-import numpy as np
+import numpy as np
+import numpy
import matplotlib.pyplot as plt
from ccpi.framework import ImageData, ImageGeometry
@@ -46,6 +47,7 @@ plt.show()
alpha = 2
#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
+
method = '0'
if method == '0':
@@ -57,7 +59,7 @@ if method == '0':
operator = BlockOperator(op1, op2, shape=(2,1) )
f1 = alpha * MixedL21Norm()
- f2 = KullbackLeibler(noisy_data)
+ f2 = 0.5 * KullbackLeibler(noisy_data)
f = BlockFunction(f1, f2 )
g = ZeroFunction()
@@ -69,9 +71,8 @@ else:
###########################################################################
operator = Gradient(ig)
f = alpha * MixedL21Norm()
- g = KullbackLeibler(noisy_data)
+ g = 0.5 * KullbackLeibler(noisy_data)
###########################################################################
-#%%
# Compute operator Norm
normK = operator.norm()
@@ -87,13 +88,11 @@ t1 = timer()
res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt)
t2 = timer()
-print(" Run memopt")
-
t3 = timer()
res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1)
t4 = timer()
-#%%
+
plt.figure(figsize=(15,15))
plt.subplot(3,1,1)
plt.imshow(res.as_array())
@@ -120,6 +119,59 @@ diff = (res1 - res).abs().as_array().max()
print(" Max of abs difference is {}".format(diff))
+#%% 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())
+
+
+ solver = SCS
+ obj = Minimize( regulariser + fidelity)
+ prob = Problem(obj)
+ result = prob.solve(verbose = True, solver = solver)
+
+ diff_cvx = numpy.abs( res.as_array() - u.value )
+
+ # Show result
+ plt.figure(figsize=(15,15))
+ plt.subplot(3,1,1)
+ plt.imshow(res.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()
+
+ print('Primal Objective (CVX) {} '.format(obj.value))
+ print('Primal Objective (PDHG) {} '.format(primal[-1]))
+
+
+
#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
#pdhg.max_iteration = 2000
#pdhg.update_objective_interval = 10
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 09:40:24 +0000