Merge pull request #55 from CUQI-DTU/sprint27_all_paper_2_notebooks

Sprint27 all paper 2 notebooks

Author: amal-ghamdi

EIT/EIT.ipynb

@@ -2,13 +2,16 @@
 import dolfin as dl
 import matplotlib.pyplot as plt
 import os
+import sys
 from cuqi.geometry import Continuous1D
 from cuqi.model import PDEModel
 from cuqi.distribution import Gaussian, JointDistribution
 from cuqi.sampler import MH
 from cuqipy_fenics.geometry import FEniCSContinuous, MaternKLExpansion,\
                                    FEniCSMappedGeometry 
 from cuqipy_fenics.pde import SteadyStateLinearFEniCSPDE
+import cuqi
+import cuqipy_fenics
 
 # Function for extracting the indices of the boundary nodes
 def extract_boundary_dofs_indices(solution_space):
@@ -55,8 +58,33 @@ def Heaviside(func):
     return G_Heavi
 
 if __name__ == "__main__":
+    # Parse command line arguments: noise level, number of samples,
+    # number of burn-in samples, thinning factor, and random seed.
+    # Parse noise level which is passed as a command line argument. Only
+    # 5, 10, and 20 percent noise levels are supported.
+
+    if len(sys.argv) != 6:
+        print("Usage: python EIT.py <noise_percent> <num_samples> <num_burnin> <thinning_factor> <random_seed>")
+        sys.exit(1)
+
+    noise_percent = int(sys.argv[1])
+    Ns = int(sys.argv[2])
+    Nb = int(sys.argv[3])
+    Nt = int(sys.argv[4])
+    seed = int(sys.argv[5])
+
+    # Check if noise level is supported
+    if noise_percent not in [5, 10, 20]:
+        print("Only 5, 10, and 20 percent noise levels are supported")
+        sys.exit(1)
+    print("Running EIT with noise level: ", noise_percent, "%")
+
+    # Print cuqi and cuqipy_fenics version
+    print("cuqi version: ", cuqi.__version__)
+    print("cuqipy_fenics version: ", cuqipy_fenics.__version__)
+
     # Fix the random seed for reproducibility 
-    np.random.seed(2)
+    np.random.seed(seed)
 
     #%% 1 setting up FEniCS
     # loading computational mesh
@@ -165,7 +193,6 @@ def observation4(sigma, v4):
 
     #%% 4 Creating the posterior distribution
     # loading signal from file
-    noise_percent = 5
     obs_data = np.load('./data/obs_circular_inclusion_2_'+str(noise_percent)+'per_noise.npz')
     b_exact = obs_data['b_exact']
     s_noise_list = np.sqrt(obs_data['sigma2']) # read the noise variance and
@@ -228,7 +255,7 @@ def observation4(sigma, v4):
     sampler = MH(posterior)
 
     # Sampling using the Metropolis-Hastings sampler
-    posterior_samples = sampler.sample_adapt(1000000)
+    posterior_samples = sampler.sample_adapt(Ns)
 
     #%% 6 visualization
     # plotting prior samples
@@ -246,7 +273,7 @@ def observation4(sigma, v4):
     plt.savefig("plot_prior_samples"+str(noise_percent)+".png")
 
     # plotting posterior samples
-    idx = np.random.permutation(1000000) # create randomized index
+    idx = np.random.permutation(Ns) # create randomized index
     f, axes = plt.subplots(1,3)
     plt.sca(axes[0])
     posterior_samples.plot(idx[0],subplots=False)
@@ -258,7 +285,7 @@ def observation4(sigma, v4):
     plt.savefig("plot_posterior_samples"+str(noise_percent)+".png")
 
     # burn-thin the samples
-    posterior_samples = posterior_samples.burnthin(200000, 4000)
+    posterior_samples = posterior_samples.burnthin(Nb, Nt)
 
     # plotting the mean
     f, axes = plt.subplots(1,2)

EIT/EIT.py

@@ -1,8 +1,14 @@
 # Instructions for running the EIT Bayesian Inverse Problem
 
-To generate the EIT results shown in the paper, run the script `EIT.py` for each of the different noise levels: $5\\%$, $10\\%$, $20\\%$. Then run the script `plot_paper_figures.py` to generate the figures 6, 7, 8, and 9 in the paper. To set the noise level, change the variable `noise_percent` in the script `EIT.py` to be either `5`, `10`, or `20`. The script `EIT.py` will then generate the sampling results for the corresponding noise level and save it in the folder `/stat/`. The script `plot_paper_figures.py` will then use these results to generate the paper figures.
+Full instructions for running the EIT Bayesian Inverse Problem are provided in 
+the notebook `EIT.ipynb`. 
 
-Observed data for each noise-levels are available in `/data/`:
-  - `/data/obs_circular_inclusion_2_5per_noise.npz` with 5 percent noise-level.
-  - `/data/obs_circular_inclusion_2_10per_noise.npz` with 10 percent noise-level.
-  - `/data/obs_circular_inclusion_2_20per_noise.npz` with 20 percent noise-level.
+Observed data are available in `./data/`:
+  - `./data/obs_circular_inclusion_2_5per_noise.npz` with 5 percent noise-level.
+  - `./data/obs_circular_inclusion_2_10per_noise.npz` with 10 percent noise-level.
+  - `./data/obs_circular_inclusion_2_20per_noise.npz` with 20 percent noise-level.
+
+Pre-computed results are available in `./stat_paper/`:
+  - `./results/stat_circular_inclusion_2_5per_noise_thinned.npz` for the 5 percent noise-level case.
+  - `./results/stat_circular_inclusion_2_10per_noise_thinned.npz` for the 10 percent noise-level case.
+  - `./results/stat_circular_inclusion_2_20per_noise_thinned.npz` for the 20 percent noise-level case.

EIT/README.md

@@ -116,14 +116,16 @@ def eval(self,values,x):
     subfigs[1].subplots_adjust(wspace=0, bottom=0.2, top=0.8)
 
 
-def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
+def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3,
+                 prior_samples_idx_to_plot=[0, 2], 
+                 posterior_samples_idx_to_plot=[0, 1]):
     cm_to_in = 1/2.54
     fig = plt.figure( figsize=(17.8*cm_to_in, 10.5*cm_to_in),layout='constrained')
     subfigs = fig.subfigures(2, 1)
 
     axes = subfigs[0].subplots(1,4,sharey=True)
     plt.sca(axes[0])
-    prior_samples.plot([ 0 ], subplots=False)
+    prior_samples.plot([prior_samples_idx_to_plot[0]], subplots=False)
     axes[0].set_xlim([-1.1,1.1])
     axes[0].set_ylim([-1.1,1.1])
     axes[0].set_xticks([-1,0,1])
@@ -134,7 +136,7 @@ def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
     axes[0].set_ylabel(r'$\xi_2$')
 
     plt.sca(axes[1])
-    cuqi_samples1.plot([ 0 ], subplots=False)
+    cuqi_samples1.plot([posterior_samples_idx_to_plot[0]], subplots=False)
     axes[1].set_ylabel('')
     axes[1].set_xlim([-1.1,1.1])
     axes[1].set_ylim([-1.1,1.1])
@@ -143,7 +145,7 @@ def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
     axes[1].set_xlabel(r'$\xi_1$')
 
     plt.sca(axes[2])
-    cuqi_samples2.plot([ 0 ], subplots=False)
+    cuqi_samples2.plot([posterior_samples_idx_to_plot[0]], subplots=False)
     axes[2].set_ylabel('')
     axes[2].set_xlim([-1.1,1.1])
     axes[2].set_ylim([-1.1,1.1])
@@ -152,9 +154,7 @@ def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
     axes[2].set_xlabel(r'$\xi_1$')
 
     plt.sca(axes[3])
-    cuqi_samples3.plot([ 0 ], subplots=False)
-    #prior_samples.plot([ 3 ], subplots=False)
-    #axes[0,3].set_yticks([])
+    cuqi_samples3.plot([posterior_samples_idx_to_plot[0]], subplots=False)
     axes[3].set_ylabel('')
     axes[3].set_xlim([-1.1,1.1])
     axes[3].set_ylim([-1.1,1.1])
@@ -164,7 +164,7 @@ def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
 
     axes = subfigs[1].subplots(1,4,sharey=True)
     plt.sca(axes[0])
-    prior_samples.plot([ 2 ], subplots=False)
+    prior_samples.plot([prior_samples_idx_to_plot[1]], subplots=False)
     axes[0].set_xlim([-1.1,1.1])
     axes[0].set_ylim([-1.1,1.1])
     axes[0].set_xticks([-1,0,1])
@@ -175,7 +175,7 @@ def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
     axes[0].set_ylabel(r'$\xi_2$')
 
     plt.sca(axes[1])
-    cuqi_samples1.plot([ -10 ], subplots=False)
+    cuqi_samples1.plot([posterior_samples_idx_to_plot[1]], subplots=False)
     axes[1].set_ylabel('')
     axes[1].set_xlim([-1.1,1.1])
     axes[1].set_ylim([-1.1,1.1])
@@ -184,7 +184,7 @@ def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
     axes[1].set_xlabel(r'$\xi_1$')
 
     plt.sca(axes[2])
-    cuqi_samples2.plot([ -10 ], subplots=False)
+    cuqi_samples2.plot([posterior_samples_idx_to_plot[1]], subplots=False)
     axes[2].set_ylabel('')
     axes[2].set_xlim([-1.1,1.1])
     axes[2].set_ylim([-1.1,1.1])
@@ -193,7 +193,7 @@ def plot_figure7(prior_samples, cuqi_samples1, cuqi_samples2, cuqi_samples3):
     axes[2].set_xlabel(r'$\xi_1$')
 
     plt.sca(axes[3])
-    cuqi_samples3.plot([ -10 ], subplots=False)
+    cuqi_samples3.plot([posterior_samples_idx_to_plot[1]], subplots=False)
     axes[3].set_ylabel('')
     axes[3].set_xlim([-1.1,1.1])
     axes[3].set_ylim([-1.1,1.1])