-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfigure6.py
127 lines (105 loc) · 3.74 KB
/
figure6.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
"""
Executing this file reproduces figure 6 from the paper.
For fixed sample size J, the error between Tikhonov regularization and Standard-, Nyström-, and SVD-EKI is plotted
for varying regularization parameters alpha.
"""
import matplotlib
matplotlib.rcParams['text.usetex'] = True
import matplotlib.pyplot as plt
import numpy as np
from inversion import *
# YOU CAN ADAPT THESE PARAMETERS
use_ray = True
n_alpha = 30
j = 2000
snr = 10. # desired signal-to-noise ratio
alpha1 = 1. # minimal alpha
c = 0.8 # sampling error is computed for alpha, alpha/c0, alpha/c0^2, ...
h = 0.01
scaling_factor = 0.25 # determines the dimension of the problem
ALPHAS_FILE = "out/figure6_alphas.csv"
STD_FILE = "out/figure6_std.csv"
NYS_FILE = "out/figure6_nys.csv"
SVD_FILE = "out/figure6_svd.csv"
FIGNAME = "out/figure6.png"
def compute_figure6():
"""
Performs computations for figure6
:return:
"""
# obtain data
y_im, y_hat_im, x_im, fwd, delta = simulate_measurement(snr, scaling_factor)
n1, n2 = x_im.shape
x = x_im.flatten()
y_hat = y_hat_im.flatten()
n = x.size
m = y_hat.size
# Display information
print(f"Image size: {n1}x{n2}")
print(f"Parameter dimension: n={n}")
print(f"Measurement dimension: m={m}")
# Set up x0 and c0
x0 = np.zeros(n)
c0 = ornstein_uhlenbeck(n1, n2, h)
print("Computing SVD of c0...")
c0_root, evals, evecs = matrix_sqrt(c0)
print("...done.")
# create list of alphas
alphas = [alpha1]
alpha = alpha1
for i in range(n_alpha-1):
alpha = alpha * c
alphas.append(alpha)
# compute Tikhonov-regularized solutions
options = {"parallel": use_ray}
traj_tik = tikhonov_list(fwd=fwd, y=y_hat, x0=x0, c0_root=c0_root, alphas=alphas, options=options)
# computed direct EKI solutions
options["j"] = j
options["sampling"] = "standard"
options["c0_root"] = c0_root
print("Standard-EKI")
traj_std = eki_list(fwd=fwd, y=y_hat, x0=x0, c0=c0, alphas=alphas, options=options)
options["sampling"] = "nystroem"
print("Nystroem-EKI")
traj_nys = eki_list(fwd=fwd, y=y_hat, x0=x0, c0=c0, alphas=alphas, options=options)
options["sampling"] = "svd"
options["c0_eigenvalues"] = evals
options["c0_eigenvectors"] = evecs
print("SVD-EKI")
traj_svd = eki_list(fwd=fwd, y=y_hat, x0=x0, c0=c0, alphas=alphas, options=options)
# compute approximation errors
traj_tik = np.array(traj_tik).T
traj_std = np.array(traj_std).T
traj_nys = np.array(traj_nys).T
traj_svd = np.array(traj_svd).T
def approximation_error(traj_hat):
return np.linalg.norm(traj_hat - traj_tik, axis=0) / np.linalg.norm(x)
e_std = approximation_error(traj_std)
e_nys = approximation_error(traj_nys)
e_svd = approximation_error(traj_svd)
# store results in csv files
np.savetxt(ALPHAS_FILE, alphas, delimiter=",")
np.savetxt(STD_FILE, e_std, delimiter=",")
np.savetxt(NYS_FILE, e_nys, delimiter=",")
np.savetxt(SVD_FILE, e_svd, delimiter=",")
def plot_figure6():
"""
Creates the plot for figure 6 and stores it under `FIGNAME`.
"""
# Load arrays from csv-files.
alphas = np.loadtxt(ALPHAS_FILE, delimiter=",")
e_std = np.loadtxt(STD_FILE, delimiter=",")
e_nys = np.loadtxt(NYS_FILE, delimiter=",")
e_svd = np.loadtxt(SVD_FILE, delimiter=",")
# plotting
# plotting
plt.rc("font", size=15) # Larger font size
plt.plot(alphas, e_std, 'ro--', label="Standard-EKI")
plt.plot(alphas, e_nys, 'bx--', label="Nyström-EKI")
plt.plot(alphas, e_svd, 'gv--', label="SVD-EKI")
plt.xlabel(r"$\alpha$")
plt.ylabel(r"$e_\mathrm{app}$")
plt.legend(loc="upper right")
plt.savefig("out/figure6.png", bbox_inches='tight')
compute_figure6()
plot_figure6()