-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain_fourier.py
96 lines (74 loc) · 2.98 KB
/
main_fourier.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
from time import time
import matplotlib.pyplot as plt
import numpy as np
from src.operators.fourier_operator import FourierOperator
from src.solvers.fw import FW
def add_snr(y0, snr, N):
signal_power = np.mean(np.square(y0))
mse_db = 10 * np.log10(signal_power) - snr
mse = 10 ** (mse_db / 10)
w = np.random.normal(0, np.sqrt(mse / 2), (N, 2))
y = y0 + w
return y
def add_psnr(y0, psnr, N):
y0_max = np.max(np.abs(y0))
mse_db = 20 * np.log10(y0_max) - psnr
mse = 10 ** (mse_db / 10)
w = np.random.normal(0, np.sqrt(mse / 2), (N, 2))
y = y0 + w
return y
if __name__ == '__main__':
np.random.seed(1)
x0 = np.array([0.1, 0.25, 0.5, 0.7, 0.9])
# x0 = np.array([-.5,-.1,.1,.5])
# a0 = np.array([0.8,0.8,0.8,0.8])
# a0 = np.array([1, 1.5, 0.5, 2, 5])
# a0 = np.array([1, 15, 0.5, -3, 5])
a0 = np.array([1, 1, 1, 1, 1])
# x0 = np.array([0.2, 0.5, 0.8])
# a0 = np.array([1, 2, 1.5])
# x0 = np.random.uniform(-0.95, 0.95, 30)
# a0 = np.random.uniform(0.5, 3, 30)
# x0 = np.array([0.1, 0.25, 0.5, 0.51, 0.7, 0.75, 0.9, 0.92])
# a0 = np.array([1, 1, 1, 1, 1, 1, 1, 1])
# a0 = np.array([-1, 0.5, 1, 1, 1, 3, 1, 1])
# x0 = np.array([-0.89, -0.7, -0.68, -0.55, -0.46, - 0.24, -0.2, -0.05, 0.1, 0.25, 0.5, 0.51, 0.7, 0.75, 0.9, 0.92])
# a0 = np.array([3, 4.5, -1.5, -3, 4, 3, 1, 2.5, -1, -1.5, 1, 1, 1, 3, 1, 1])
# a0 = np.abs(a0)
bounds = np.array([0, 1])
# bounds = np.array([-1, 1])
N = 10 * len(x0)
freq_bounds = np.array([-10, 10])
forward_op = FourierOperator.get_RandomFourierOperator(x0, N, freq_bounds)
# Get measurements
y0 = forward_op(a0)
# add noise
psnr = 20
y = add_psnr(y0, psnr, N)
# y = add_snr(y0, psnr, N)
# Get lambda
lambda_max = max(abs((forward_op.adjoint(y))))
print("lambda_max = ", lambda_max)
lambda_ = 0.1 * lambda_max
x_dim = 1
lambdas = [0.001, 0.01, 0.02, 0.1]
options = {"initialization": "smoothing", "polyatomic": False, "swarm": False, "sliding": True, "positivity_constraint": False,
"max_iter": 20, "dual_certificate_tol": 1e-2, "smooth_sigma": 4, "n_particles": 10}
solver = FW(y, forward_op, lambda_, x_dim, bounds=bounds, verbose=False, show_progress=False, options=options)
t1 = time()
solver.fit()
print("Time: ", time() - t1)
solver.time_results()
solver.plot(x0, a0)
solver.flat_norm_results(x0, a0, lambdas)
solver.plot_solution(x0, a0)
options = {"initialization": "smoothing", "polyatomic": True, "swarm": False, "sliding": False, "positivity_constraint": False,
"max_iter": 20, "dual_certificate_tol": 1e-2, "smooth_sigma": 2.5}
solver = FW(y, forward_op, lambda_, x_dim, bounds=bounds, verbose=False, show_progress=False, options=options)
t1 = time()
solver.fit()
print("Time: ", time() - t1)
solver.time_results()
solver.plot(x0, a0)
solver.flat_norm_results(x0, a0, lambdas)
solver.plot_solution(x0, a0)