add spdhg + Kullback Leibler + example without subsets

examples/Python/PET/interactive/spdhg_reconstruction.py

@@ -113,105 +113,47 @@ def make_cylindrical_FOV(image):
 # Since we are not using a penalty, or prior in this example, it
 # defaults to using MLEM, but we will modify it to OSEM
 
-import pCIL
+import pCIL  # the code from this module needs to be imported somehow differently
 
-data = acquired_data.as_array()
+data = acquired_data
 background = 0 * data.copy()
 
-f = pCIL.KullbackLeibler(data, background)
+f = [pCIL.KullbackLeibler(data, background)]
 g = pCIL.ZeroFun()
 
-init_image = image.asarray().clone()
+init_image = 0 * image.copy()
 z = init_image.copy()
-y = data.copy()
-A = am
+y = [0 * data.copy()]
+A = [am]
 
 niter = 10
 
-L = pCIL.PowerMethodNonsquare(A, 10)
+L = pCIL.PowerMethodNonsquare(A[0], 10, x0=image.copy())
 L *= 1.05
 tau = 1 / L
 sigma = [1 / L]
 
 # Selection function
-def fun_select(x):
-    return [int(numpy.random.choice(len(A), 1, p=1 / len(A)))]
-
+def fun_select(k):
+    n = len(A)
+    return [int(numpy.random.choice(n, 1, p=[1./n] * n))]
+
+init_image = 0 * image.copy()   
 recon = pCIL.spdhg(init_image, y, z, f, g, A, tau, sigma, fun_select)
 
+# %%
 for i in range(niter):
     print(i)
     recon.update()
 
+# a = am.forward(x, subset_num=0, num_subsets=10)
+# b = am.forward(x)
+# b.as_array().shape
 #%%  create initial image
 # we could just use a uniform image but here we will create a disk with a different
 # initial value (this will help the display later on)
-init_image=image.clone()
-init_image.fill(cmax/4)
-make_cylindrical_FOV(init_image)
-# display
-idata = init_image.as_array()
+idata = recon.x.as_array()
 slice=idata.shape[0]/2;
 plt.figure()
 imshow(idata[slice,:,:],[0,cmax], 'initial image');
-
-#%% reconstruct the image 
-reconstructed_image=init_image.clone()
-# set up the reconstructor
-recon.set_up(reconstructed_image)
-# do actual recon
-recon.reconstruct(reconstructed_image)
-
-#%% bitmap display of images
-reconstructed_array=reconstructed_image.as_array()
-
-plt.figure();
-plt.subplot(1,2,1);
-imshow(image_array[slice,:,:,], [0,cmax*1.2],'emission image');
-plt.subplot(1,2,2);
-imshow(reconstructed_array[slice,:,:,], [0,cmax*1.2], 'reconstructed image');
-
-
-#%% Generate a noisy realisation of the data
-noisy_array=numpy.random.poisson(acquisition_array).astype('float64')
-print(' Maximum counts in the data: %d' % noisy_array.max())
-# stuff into a new AcquisitionData object
-noisy_data = acquired_data.clone()
-noisy_data.fill(noisy_array);
-
-#%% Display bitmaps of the middle sinogram
-plt.figure()
-plt.subplot(1,2,1);
-imshow(acquisition_array[slice,:,:,], [0,acquisition_array.max()], 'original');
-plt.subplot(1,2,2);
-imshow(noisy_array[slice,:,:,], [0,acquisition_array.max()], 'noisy');
-
-#%% reconstruct the noisy data
-obj_fun.set_acquisition_data(noisy_data)
-# We could save the data to file if we wanted to, but currently we don't.
-# recon.set_output_filename_prefix('reconstructedImage_noisydata')
-
-noisy_reconstructed_image=init_image.clone()
-recon.reconstruct(noisy_reconstructed_image)
-#%% bitmap display of images
-noisy_reconstructed_array=noisy_reconstructed_image.as_array()
-
-plt.figure();
-plt.subplot(1,2,1);
-imshow(reconstructed_array[slice,:,:,], [0,cmax*1.2], 'no noise');
-plt.subplot(1,2,2);
-imshow(noisy_reconstructed_array[slice,:,:,], [0,cmax*1.2], 'with noise');
-
-#%% run same reconstruction but saving images and objective function values every sub-iteration
-num_subiters = 64;
-all_osem_images = numpy.ndarray(shape=(num_subiters+1,) + idata.shape );
-current_image = init_image.clone()
-osem_objective_function_values = [ obj_fun.value(current_image) ]
-all_osem_images[0,:,:,:] =  current_image.as_array();
-for iter in range(1, num_subiters+1):
-    recon.update(current_image);
-
-    obj_fun_value = obj_fun.value(current_image);
-    osem_objective_function_values.append(obj_fun_value);
-    all_osem_images[iter,:,:,:] =  current_image.as_array();
 

src/pCIL/__init__.py

@@ -0,0 +1,2 @@
+from functions import *
+from spdhg import *

src/pCIL/functions.py

@@ -0,0 +1,211 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jul 27 11:08:29 2018
+
+@author: sirfuser
+"""
+
+import numpy
+
+class Function(object):
+    def __init__(self):
+        pass
+
+    def __call__(self, x):
+        raise NotImplementedError
+
+    def grad(self, x):
+        raise NotImplementedError
+
+    def prox(self, x, tau, out=None):
+        raise NotImplementedError
+
+    @property
+    def convex_conj():
+        raise NotImplementedError
+
+
+class ZeroFun(Function):
+
+    def __init__(self):
+        super(ZeroFun, self).__init__()
+
+    def __call__(self, x):
+        return 0
+
+    def prox(self, x, tau):
+        return x
+
+
+class KullbackLeibler(Function):
+    def __init__(self, data, background):
+        self.data = data
+        self.background = background
+        self.__offset = None
+
+    def __call__(self, x):
+        """Return the KL-diveregnce in the point ``x``.
+
+        If any components of ``x`` is non-positive, the value is positive
+        infinity.
+
+        Needs one extra array of memory of the size of `prior`.
+        """
+
+        # define short variable names
+        y = self.data
+        r = self.background
+
+        # Compute
+        #   sum(x + r - y + y * log(y / (x + r)))
+        # = sum(x - y * log(x + r)) + self.offset
+        # Assume that
+        #   x + r > 0
+
+        # sum the result up
+        obj = numpy.sum(x - y * numpy.log(x + r)) + self.offset()
+
+        if numpy.isnan(obj):
+            # In this case, some element was less than or equal to zero
+            return numpy.inf
+        else:
+            return obj
+
+    @property
+    def convex_conj(self):
+        """The convex conjugate functional of the KL-functional."""
+        return KullbackLeiblerConvexConjugate(self.data, self.background)
+
+    def offset(self):
+        """The offset which is independent of the unknown."""
+
+        if self.__offset is None:
+            tmp = self.domain.element()
+
+            # define short variable names
+            y = self.data
+            r = self.background
+
+            tmp = self.domain.element(numpy.maximum(y, 1))
+            tmp = r - y + y * numpy.log(tmp)
+
+            # sum the result up
+            self.__offset = numpy.sum(tmp)
+
+        return self.__offset
+
+#     def __repr__(self):
+#         """to be added???"""
+#         """Return ``repr(self)``."""
+        # return '{}({!r}, {!r}, {!r})'.format(self.__class__.__name__,
+          ##                                    self.domain, self.data,
+           #                                   self.background)
+
+
+class KullbackLeiblerConvexConjugate(Function):
+    """The convex conjugate of Kullback-Leibler divergence functional.
+
+    Notes
+    -----
+    The functional :math:`F^*` with prior :math:`g>0` is given by:
+
+    .. math::
+        F^*(x)
+        =
+        \\begin{cases}
+            \\sum_{i} \left( -g_i \ln(1 - x_i) \\right)
+            & \\text{if } x_i < 1 \\forall i
+            \\\\
+            +\\infty & \\text{else}
+        \\end{cases}
+
+    See Also
+    --------
+    KullbackLeibler : convex conjugate functional
+    """
+
+    def __init__(self, data, background):
+        self.data = data
+        self.background = background
+
+    def __call__(self, x):
+        y = self.data
+        r = self.background
+
+        tmp = numpy.sum(- x * r - y * numpy.log(1 - x))
+
+        if numpy.isnan(tmp):
+            # In this case, some element was larger than or equal to one
+            return numpy.inf
+        else:
+            return tmp
+
+
+    def prox(self, x, tau, out=None):
+        # Let y = data, r = background, z = x + tau * r
+        # Compute 0.5 * (z + 1 - sqrt((z - 1)**2 + 4 * tau * y))
+        # Currently it needs 3 extra copies of memory.
+
+        if out is None:
+            out = x.clone()
+
+        # define short variable names
+        y = self.data.as_array()
+        r = self.background.as_array()
+        x = x.as_array()
+
+        try:
+            taua = tau.as_array()
+        except:
+            taua = tau
+
+        z = x + tau * r
+
+        out.fill(0.5 * (z + 1 - numpy.sqrt((z - 1) ** 2 + 4 * taua * y)))
+
+        return out 
+
+    @property
+    def convex_conj(self):
+        return KullbackLeibler(self.data, self.background)
+
+
+def mult(x, y):
+    try:
+        xa = x.as_array()
+    except:
+        xa = x
+
+    out = y.clone()
+    out.fill(xa * y.as_array())
+
+    return out
+
+def PowerMethodNonsquare(op, numiters, x0=None):
+    # Initialise random
+    # Jakob's
+    #inputsize = op.size()[1]
+    #x0 = ImageContainer(numpy.random.randn(*inputsize)
+    # Edo's
+    #vg = ImageGeometry(voxel_num_x=inputsize[0],
+    #                   voxel_num_y=inputsize[1], 
+    #                   voxel_num_z=inputsize[2])
+    #
+    #x0 = ImageData(geometry = vg, dimension_labels=['vertical','horizontal_y','horizontal_x'])
+    #print (x0)
+    #x0.fill(numpy.random.randn(*x0.shape))
+
+    if x0 is None:
+        x0 = op.create_image_data()
+
+    #s = numpy.zeros(numiters)
+    # Loop
+    for it in numpy.arange(numiters):
+        x1 = op.adjoint(op.direct(x0))
+        x1norm = numpy.sqrt(x1.dot(x1))
+        #print ("x0 **********" ,x0)
+        #print ("x1 **********" ,x1)
+        #s[it] = x1.dot(x0) / x0.dot(x0)
+        x0 = (1.0/x1norm)*x1
+    return numpy.sqrt(x1norm)