added test for autocorrelation function for gibbs sampler

ver-python/fraction_missing_information.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Mar 26 19:25:12 2020
+
+@author: [email protected]
+"""
+
+import numpy as np
+
+def fraction_missing_pixel(distribution, syst_matrix, avg_scattered, projector):
+    """
+    Asymptotic value for the Bayesian fraction of missing information.
+    
+    Possibly this function is not correct. Some theoretical analysis 
+    is necessary.
+    
+    Parameters
+    ----------
+    distribution : TYPE
+        DESCRIPTION.
+    syst_matrix : TYPE
+        DESCRIPTION.
+    projector : TYPE
+        DESCRIPTION.
+
+    Returns
+    -------
+    None.
+
+    """
+    npixels = np.sqrt(syst_matrix.shape[1]).astype(int)
+
+    projector_vec = np.reshape(projector, (npixels**2, 1))
+    distribution_vec = np.reshape(distribution, (npixels**2, 1))
+    scattered_vec = np.reshape(avg_scattered, (syst_matrix.shape[0],1))
+
+    vec_a = syst_matrix.dot(projector_vec)
+    pixel_intensity = np.sum(np.multiply(distribution_vec, projector_vec))
+
+    normalization = syst_matrix.dot(distribution_vec) + scattered_vec
+
+    multiplier_to_sensitivity = np.ones((syst_matrix.shape[0], 1)) - np.divide(pixel_intensity*vec_a, 
+                                                      normalization)
+
+    main_nominator = np.sum(np.multiply(vec_a, multiplier_to_sensitivity))
+    main_denominator = np.sum(vec_a)
+
+    return (1-1./(1. + main_nominator/main_denominator))
+
+def fraction_missing_linear(distribution, syst_matrix, projector):
+    """
+    
+
+    Parameters
+    ----------
+    distribution : TYPE
+        DESCRIPTION.
+    syst_matrix : TYPE
+        DESCRIPTION.
+    projector : TYPE
+        DESCRIPTION.
+
+    Returns
+    -------
+    None.
+
+    """
+    return None

ver-python/gibbs_sampler_emission.py

@@ -40,9 +40,9 @@ def gibbs_sampler_emission_gamma(sinogram_counts, syst_matrix, avg_scattered,
     syst_matrix : matrix of type (nshift*nphi, dim_im*dim_im) : system matrix of observation 
                                                                 data
     avg_scattered : matrix of type (nshift, nphi) : average number of scattered photons per LOR
-    gamma_prior_params : tuple of type (alpha, beta) : parameters for Gamma-distribution prior
+    gamma_prior_params : tuple of type (float, float) : parameters (shape, rate) for Gamma-distribution prior
+    init_point : matrix of type (dim_im, dim_im) : initial guess for the distribution
     burn_in_size : integer : number of iterations for burn-in 
-    init_point : matrix of type (dim_im, dim_im) : initial guess for the isotope distribution
     sample_size : integer : number of samples
     
 
@@ -84,7 +84,7 @@ def gibbs_sampler_emission_gamma(sinogram_counts, syst_matrix, avg_scattered,
                                                                                                   keepdims=True)
             multinomial_prob_matrix = np.append(multinomial_prob_matrix, multinomial_prob_matrix_last_column, axis=1)
 
-        # generate backprojected data n_ij, scattered
+        # generate backprojected data n_ij, scattered photons
             for i in range(nshift*nphi):
                 backprojection_data[i, :] = np.random.multinomial(sinogram_counts_vec[i], multinomial_prob_matrix[i, :])[:-1] 
         # end of Step 1
@@ -99,10 +99,9 @@ def gibbs_sampler_emission_gamma(sinogram_counts, syst_matrix, avg_scattered,
         # save the sample
             if (iteration > burn_in_size-1):
                 output_array[:, :, iteration-burn_in_size] = np.reshape(current_density_vec, (npixels, npixels))
-
-        # end of loop 
+    # end of iteration loop 
 
-        # clean memory in a loop (not efficient)
+        # clean memory in a loop (not very efficient, better set variables before loop)
             del(backprojection_data)
             del(multinomial_prob_matrix)
             del(multinomial_prob_matrix_nominator)

ver-python/test_gibbs_sampler_emission.py

@@ -34,14 +34,14 @@
 system_matrix = radon_transform2d_xray_matrix(64, 64, 64, 1.0)
 
 sinogram_vector = system_matrix.dot(np.reshape(image, (64*64, 1)))
-noise_sinogram_vector = generate_noise_emission(sinogram_vector, 200, 10)
+noise_sinogram_vector = generate_noise_emission(sinogram_vector, 20000000, 10)
 noise_sinogram = np.reshape(noise_sinogram_vector, (64, 64))
 
 avg_scattered = 10*np.ones((64,64))
 gamma_prior_params = (1.,1.)
 init_point = np.ones((64,64))
 burn_in_size = 0
-sample_size = 100
+sample_size = 1000
 
 posterior_sample = gibbs_sampler_emission_gamma(noise_sinogram, system_matrix, avg_scattered, gamma_prior_params, init_point, burn_in_size, sample_size)
 
@@ -55,3 +55,51 @@
 plt.title("Posterior standard deviation")
 plt.colorbar()
 
+# autocorrelation test
+
+# create image with one "hot pixel"
+hot_image = image 
+hot_image[32,22] = 10. # create hot pixel
+
+# for comparison we choose also a "cold" pixel
+cold_projector = np.zeros((64, 64)) 
+cold_projector[32,32] = 1.  # projector for "cold" pixel
+
+hot_projector = np.zeros((64, 64))
+hot_projector[32,22] = 1. # projector for "hot" pixel
+
+hot_sinogram_vector = system_matrix.dot(np.reshape(hot_image, (64*64, 1)))
+hot_noise_sinogram_vector = generate_noise_emission(hot_sinogram_vector, 2e6, 10)
+hot_noise_sinogram = np.reshape(hot_noise_sinogram_vector, (64, 64))
+
+avg_scattered = 10*np.ones((64,64))
+gamma_prior_params = (1.,1.)
+init_point = np.ones((64,64))
+burn_in_size = 0
+sample_size = 1e3
+
+hot_posterior_sample = gibbs_sampler_emission_gamma(hot_noise_sinogram, system_matrix, avg_scattered, gamma_prior_params, init_point, burn_in_size, sample_size)
+
+autocorrelation_array_hot = np.array([])
+autocorrelation_array_cold = np.array([])
+
+for i in range(sample_size):
+    autocorrelation_array_hot = np.append(autocorrelation_array_hot, np.sum(np.multiply(hot_projector, hot_posterior_sample[:, :, i])[:]))
+    autocorrelation_array_cold = np.append(autocorrelation_array_cold, np.sum(np.multiply(cold_projector, hot_posterior_sample[:, :, i])[:]))
+
+# autocorrelation function
+def acf(x, length=100):
+    return np.array([1]+[np.corrcoef(x[:-i], x[i:])[0,1]  \
+        for i in range(1, length)])
+
+autocorrelation_lag_array_hot = acf(autocorrelation_array_hot)
+autocorrelation_lag_array_cold = acf(autocorrelation_array_cold)
+
+fig3 = plt.plot(autocorrelation_lag_array_hot)
+plt.plot(autocorrelation_lag_array_cold)
+plt.legend(["Hot pixel", "Cold pixel"])
+plt.xlabel("lag")
+plt.ylabel("correlation")
+plt.title("autocorrelation function")
+
+