-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfigure5.py
118 lines (97 loc) · 3.8 KB
/
figure5.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
"""
Executing this file reproduces figure 5 from the paper.
For fixed regularization parameter alpha, the error between
Tikhonov regularization and Standard-, Nyström-, and SVD-EKI is plotted for varying sample size J
"""
import matplotlib
matplotlib.rcParams['text.usetex'] = True
from matplotlib import pyplot as plt
import numpy as np
from inversion import *
# YOU CAN ADAPT THESE PARAMETERS
use_ray = True # toggles parallelization; set False for final experiment.
alpha = 0.03
snr = 10. # desired signal-to-noise ratio
scaling_factor = 0.25 # determines the dimension; set to 0.4 for final experiment
h = 0.01
def compute_figure5():
"""
Performs computations for figure 5 and stores the results in "out" as "figure5_sizes.csv", "figure5_eki.csv",
"figure5_nys.csv" and "figure5_svd.csv".
: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.")
# determine a good value for the regularization parameter alpha using the discrepancy principle
# and Tikhonov regularization
options = {"parallel": use_ray}
x_tik = tikhonov(fwd=fwd, y=y_hat, x0=x0, c0_root=c0_root, alpha=alpha, options=options)
# apply EKI
traj_std = []
traj_nys = []
traj_svd = []
sizes = [100, 500, 1000, 1500, 2000, 2500, 3000, 5000, 8000]
options["c0_root"] = c0_root
options["c0_eigenvalues"] = evals
options["c0_eigenvectors"] = evecs
for j in sizes:
options["j"] = j
print("Sample size: ", j)
options["sampling"] = "standard"
x_std = direct_eki(fwd=fwd, y=y_hat, x0=x0, c0=c0, alpha=alpha, options=options)
traj_std.append(x_std)
options["sampling"] = "nystroem"
x_nys = direct_eki(fwd=fwd, y=y_hat, x0=x0, c0=c0, alpha=alpha, options=options)
traj_nys.append(x_nys)
options["sampling"] = "svd"
x_svd = direct_eki(fwd=fwd, y=y_hat, x0=x0, c0=c0, alpha=alpha, options=options)
traj_svd.append(x_svd)
# compute approximation errors
def approximation_error(x_hat):
return np.linalg.norm(x_hat - x_tik[:,np.newaxis], axis=0) / np.linalg.norm(x)
traj_std = np.array(traj_std).T
traj_nys = np.array(traj_nys).T
traj_svd = np.array(traj_svd).T
e_eki = approximation_error(traj_std)
e_nys = approximation_error(traj_nys)
e_svd = approximation_error(traj_svd)
# Store results
np.savetxt("out/figure5_sizes.csv", sizes, delimiter=",")
np.savetxt("out/figure5_eki.csv", e_eki, delimiter=",")
np.savetxt("out/figure5_nys.csv", e_nys, delimiter=",")
np.savetxt("out/figure5_svd.csv", e_svd, delimiter=",")
def plot_figure5():
"""
Plots figure 5 using csv-files and stores it as 'out/figure5.png'.
:return:
"""
sizes = np.loadtxt("out/figure5_sizes.csv", delimiter=",")
e_eki = np.loadtxt("out/figure5_eki.csv", delimiter=",")
e_nys = np.loadtxt("out/figure5_nys.csv", delimiter=",")
e_svd = np.loadtxt("out/figure5_svd.csv", delimiter=",")
# plotting
plt.rc("font", size=15) # Larger font size
plt.plot(sizes, e_eki, 'ro--', label="Standard-EKI")
plt.plot(sizes, e_nys, 'bx--', label="Nyström-EKI")
plt.plot(sizes, e_svd, 'gv--', label="SVD-EKI")
plt.xlabel(r"$J$")
plt.ylabel(r"$e_\mathrm{app}$")
plt.legend(loc="upper right")
plt.savefig("out/figure5.png", bbox_inches='tight')
compute_figure5()
plot_figure5()