Skip to content

Commit acfba10

Browse files
author
Alexis Montoison
committed
Optimized version of the artificial compressed sensing problems
1 parent 24fa6a1 commit acfba10

File tree

6 files changed

+231
-26
lines changed

6 files changed

+231
-26
lines changed

‎fft_example_1D.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ function fft_example_1D(Nt::Int; gpu::Bool=false, rdft::Bool=false, check::Bool=
2929
centers = centering(DFTdim, DFTsize, missing_prob)
3030
radius = 1
3131
index_missing, z_zero = punching(DFTdim, DFTsize, centers, radius, y)
32+
# println("length(index_missing) = ", length(index_missing))
3233
end
3334

3435
M_perptz = M_perp_tz_wei(DFTdim, DFTsize, z_zero) # M_perptz

‎fft_example_2D.jl‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -20,13 +20,14 @@ function fft_example_2D(Nt::Int, Ns::Int; gpu::Bool=false, rdft::Bool=false, che
2020

2121
# randomly generate missing indices
2222
if check
23-
index_missing_Cartesian = Int[]
23+
index_missing = Int[]
2424
z_zero = y
2525
else
2626
missing_prob = 0.15
2727
centers = centering(DFTdim, DFTsize, missing_prob)
2828
radius = 1
29-
index_missing_Cartesian, z_zero = punching(DFTdim, DFTsize, centers, radius, y)
29+
index_missing, z_zero = punching(DFTdim, DFTsize, centers, radius, y)
30+
# println("length(index_missing) = ", length(index_missing))
3031
end
3132

3233
# unify parameters for barrier method
@@ -42,7 +43,7 @@ function fft_example_2D(Nt::Int, Ns::Int; gpu::Bool=false, rdft::Bool=false, che
4243
eps_barrier = 1e-6
4344
mu_barrier = 10
4445

45-
parameters = FFTParameters(DFTdim, DFTsize, M_perptz, lambda, index_missing_Cartesian, alpha_LS, gamma_LS, eps_NT, mu_barrier, eps_barrier)
46+
parameters = FFTParameters(DFTdim, DFTsize, M_perptz, lambda, index_missing, alpha_LS, gamma_LS, eps_NT, mu_barrier, eps_barrier)
4647

4748
t_init = 1
4849
beta_init = zeros(prod(DFTsize))

‎fft_example_3D.jl‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -21,13 +21,14 @@ function fft_example_3D(N1::Int, N2::Int, N3::Int; gpu::Bool=false, rdft::Bool=f
2121

2222
# randomly generate missing indices
2323
if check
24-
index_missing_Cartesian = Int[]
24+
index_missing = Int[]
2525
z_zero = y
2626
else
2727
missing_prob = 0.15
2828
centers = centering(DFTdim, DFTsize, missing_prob)
2929
radius = 1
30-
index_missing_Cartesian, z_zero = punching(DFTdim, DFTsize, centers, radius, y)
30+
index_missing, z_zero = punching(DFTdim, DFTsize, centers, radius, y)
31+
# println("length(index_missing) = ", length(index_missing))
3132
end
3233

3334
# unify parameters for barrier method
@@ -43,7 +44,7 @@ function fft_example_3D(N1::Int, N2::Int, N3::Int; gpu::Bool=false, rdft::Bool=f
4344
eps_barrier = 1e-6
4445
mu_barrier = 10
4546

46-
parameters = FFTParameters(DFTdim, DFTsize, M_perptz, lambda, index_missing_Cartesian, alpha_LS, gamma_LS, eps_NT, mu_barrier, eps_barrier)
47+
parameters = FFTParameters(DFTdim, DFTsize, M_perptz, lambda, index_missing, alpha_LS, gamma_LS, eps_NT, mu_barrier, eps_barrier)
4748

4849
t_init = 1
4950
beta_init = zeros(prod(DFTsize))
@@ -67,7 +68,7 @@ function fft_example_3D(N1::Int, N2::Int, N3::Int; gpu::Bool=false, rdft::Bool=f
6768
nlp_scaling=false,
6869
dual_initialized=true,
6970
richardson_max_iter=0,
70-
tol=1e-6,
71+
tol=1e-8,
7172
richardson_tol=Inf,
7273
)
7374
results = ipm_solve!(solver)

‎fft_model.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -90,7 +90,8 @@ end
9090
include("kkt.jl")
9191
include("fft_wei.jl")
9292
include("fft_utils.jl")
93-
include("punching_centering.jl")
93+
# include("punching_centering.jl")
94+
include("punching_centering_v2.jl")
9495

9596
function NLPModels.cons!(nlp::FFTNLPModel, x::AbstractVector, c::AbstractVector)
9697
increment!(nlp, :neval_cons)

‎punching_centering_v2.jl‎

Lines changed: 201 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,201 @@
1+
## punching
2+
function punching(DFTdim, DFTsize, centers, radius, data)
3+
if radius == 1
4+
return punching_optimized(DFTdim, DFTsize, centers, data)
5+
end
6+
if DFTdim == 1
7+
index_missing = punching1D(DFTsize, centers, radius)
8+
data[index_missing] .= 0
9+
elseif DFTdim == 2
10+
index_missing = punching2D(DFTsize, centers, radius)
11+
data[index_missing] .= 0
12+
else
13+
index_missing = punching3D(DFTsize, centers, radius)
14+
data[index_missing] .= 0
15+
end
16+
return index_missing, data
17+
end
18+
19+
function punching1D(DFTsize, centers, radius)
20+
N = prod(DFTsize)
21+
index_missing = Vector{CartesianIndex{1}}(undef, N)
22+
pos = 0
23+
for i = 1:DFTsize[1]
24+
for center in centers
25+
if abs(i - center[1]) <= radius
26+
pos = pos + 1
27+
index_missing[pos] = CartesianIndex{1}(i)
28+
end
29+
end
30+
end
31+
resize!(index_missing, pos)
32+
return index_missing
33+
end
34+
35+
function punching2D(DFTsize, centers, radius)
36+
N = prod(DFTsize)
37+
index_missing = Vector{CartesianIndex{2}}(undef, N)
38+
pos = 0
39+
for i = 1:DFTsize[1]
40+
for j = 1:DFTsize[2]
41+
for center in centers
42+
if (center[1] - i)^2 + (center[2] - j)^2 <= radius^2
43+
pos = pos + 1
44+
index_missing[pos] = CartesianIndex{2}(i, j)
45+
end
46+
end
47+
end
48+
end
49+
resize!(index_missing, pos)
50+
return index_missing
51+
end
52+
53+
function punching3D(DFTsize, centers, radius)
54+
N = prod(DFTsize)
55+
index_missing = Vector{CartesianIndex{3}}(undef, N)
56+
pos = 0
57+
for i = 1:DFTsize[1]
58+
for j = 1:DFTsize[2]
59+
for k = 1:DFTsize[3]
60+
for center in centers
61+
if (center[1] - i)^2 + (center[2] - j)^2 + (center[3] - k)^2 <= radius^2
62+
pos = pos + 1
63+
index_missing[pos] = CartesianIndex{3}(i, j, k)
64+
end
65+
end
66+
end
67+
end
68+
end
69+
resize!(index_missing, pos)
70+
return index_missing
71+
end
72+
73+
## centering
74+
function centering(DFTdim, DFTsize, missing_prob)
75+
if DFTdim == 1
76+
return center_1d(DFTsize, missing_prob)
77+
elseif DFTdim == 2
78+
return center_2d(DFTsize, missing_prob)
79+
else
80+
return center_3d(DFTsize, missing_prob)
81+
end
82+
end
83+
84+
function center_1d(DFTsize, missing_prob)
85+
N = prod(DFTsize)
86+
n = N*missing_prob/3
87+
stepsize = ceil(N/n) |> Int
88+
centers = CartesianIndices((1:stepsize:N))
89+
return centers
90+
end
91+
92+
function center_2d(DFTsize, missing_prob)
93+
N = prod(DFTsize)
94+
Nt = DFTsize[1]
95+
Ns = DFTsize[2]
96+
n = (N*missing_prob/5)^(1/2)
97+
stepsize1 = ceil(Nt/n) |> Int
98+
stepsize2 = ceil(Ns/n) |> Int
99+
centers = CartesianIndices((1:stepsize1:Nt, 1:stepsize2:Ns))
100+
return centers
101+
end
102+
103+
function center_3d(DFTsize, missing_prob)
104+
N = prod(DFTsize)
105+
N1 = DFTsize[1]
106+
N2 = DFTsize[2]
107+
N3 = DFTsize[3]
108+
n = (N*missing_prob/7)^(1/3)
109+
stepsize1 = ceil(N1/n) |> Int
110+
stepsize2 = ceil(N2/n) |> Int
111+
stepsize3 = ceil(N3/n) |> Int
112+
centers = CartesianIndices((1:stepsize1:N1, 1:stepsize2:N2, 1:stepsize3:N3))
113+
return centers
114+
end
115+
116+
## punching_optimized
117+
function punching_optimized(DFTdim, DFTsize, centers, data)
118+
if DFTdim == 1
119+
index_missing = punching_optimized_1D(DFTsize, centers)
120+
data[index_missing] .= 0
121+
elseif DFTdim == 2
122+
index_missing = punching_optimized_2D(DFTsize, centers)
123+
data[index_missing] .= 0
124+
else
125+
index_missing = punching_optimized_3D(DFTsize, centers)
126+
data[index_missing] .= 0
127+
end
128+
return index_missing, data
129+
end
130+
131+
function punching_optimized_1D(DFTsize, centers)
132+
ncenters = prod(centers |> size)
133+
index_missing = Vector{CartesianIndex{1}}(undef, 3*ncenters)
134+
Nx = DFTsize[1]
135+
pos = 0
136+
for center in centers
137+
for i in center[1]
138+
for i2 = i-1:i+1
139+
if 1 <= i2 <= Nx
140+
pos = pos + 1
141+
index_missing[pos] = CartesianIndex{1}(i2)
142+
end
143+
end
144+
end
145+
end
146+
resize!(index_missing, pos)
147+
return index_missing
148+
end
149+
150+
function punching_optimized_2D(DFTsize, centers)
151+
ncenters = prod(centers |> size)
152+
index_missing = Vector{CartesianIndex{2}}(undef, 5*ncenters)
153+
Nx = DFTsize[1]
154+
Ny = DFTsize[2]
155+
pos = 0
156+
for center in centers
157+
for i in center[1]
158+
for j in center[2]
159+
for i2 = i-1:i+1
160+
for j2 = j-1:j+1
161+
if (1 <= i2 <= Nx) && (1 <= j2 <= Ny) && (abs(i2 - i) + abs(j2 - j) <= 1)
162+
pos = pos + 1
163+
index_missing[pos] = CartesianIndex{2}(i2, j2)
164+
end
165+
end
166+
end
167+
end
168+
end
169+
end
170+
resize!(index_missing, pos)
171+
return index_missing
172+
end
173+
174+
function punching_optimized_3D(DFTsize, centers)
175+
ncenters = prod(centers |> size)
176+
index_missing = Vector{CartesianIndex{3}}(undef, 7*ncenters)
177+
Nx = DFTsize[1]
178+
Ny = DFTsize[2]
179+
Nz = DFTsize[3]
180+
pos = 0
181+
for center in centers
182+
for i in center[1]
183+
for j in center[2]
184+
for k in center[3]
185+
for i2 = i-1:i+1
186+
for j2 = j-1:j+1
187+
for k2 = k-1:k+1
188+
if (1 <= i2 <= Nx) && (1 <= j2 <= Ny) && (1 <= k2 <= Nz) && (abs(i2 - i) + abs(j2 - j) + abs(k2 - k) <= 1)
189+
pos = pos + 1
190+
index_missing[pos] = CartesianIndex{3}(i2, j2, k2)
191+
end
192+
end
193+
end
194+
end
195+
end
196+
end
197+
end
198+
end
199+
resize!(index_missing, pos)
200+
return index_missing
201+
end

‎unit_tests.jl‎

Lines changed: 18 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -26,14 +26,14 @@ if dim1
2626
@testset "1D -- CPU -- rdft=$rdft -- $N" begin
2727
nlp, solver, results = fft_example_1D(N; gpu=false, rdft, check=true)
2828

29-
z2 = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z, Int[]; rdft)
29+
z2 = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
3030
@test z2 ≈ z
3131

32-
res1 = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z; rdft)
32+
res1 = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z, nlp.fft_timer, nlp.mapping_timer; rdft)
3333
@test norm(res1) ≈ norm(z)
3434
@test res1_wei ≈ res1
3535

36-
res2 = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z, Int[]; rdft)
36+
res2 = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
3737
@test norm(res2) ≈ norm(z)
3838
@test res2_wei ≈ res2
3939
end
@@ -44,14 +44,14 @@ if dim1
4444

4545
z_gpu = CuArray(z)
4646

47-
z2_gpu = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z_gpu, Int[]; rdft)
47+
z2_gpu = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z_gpu, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
4848
@test z2_gpu ≈ z_gpu
4949

50-
res1_gpu = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z_gpu; rdft)
50+
res1_gpu = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z_gpu, nlp.fft_timer, nlp.mapping_timer; rdft)
5151
@test norm(res1_gpu) ≈ norm(z_gpu)
5252
@test res1_wei ≈ collect(res1_gpu)
5353

54-
res2_gpu = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z_gpu, Int[]; rdft)
54+
res2_gpu = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 1, (N,), z_gpu, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
5555
@test norm(res2_gpu) ≈ norm(z_gpu)
5656
@test res2_wei ≈ collect(res2_gpu)
5757
end
@@ -78,14 +78,14 @@ if dim2
7878
@testset "2D -- CPU -- rdft=$rdft -- $N1 × $N2" begin
7979
nlp, solver, results = fft_example_2D(N1, N2; gpu=false, rdft, check=true)
8080

81-
z2 = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z, Int[]; rdft)
81+
z2 = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
8282
@test z2 ≈ z
8383

84-
res1 = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), reshape(z, (N1, N2)); rdft)
84+
res1 = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), reshape(z, (N1, N2)), nlp.fft_timer, nlp.mapping_timer; rdft)
8585
@test norm(res1) ≈ norm(z)
8686
@test res1_wei ≈ res1
8787

88-
res2 = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z, Int[]; rdft)
88+
res2 = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
8989
@test norm(res2) ≈ norm(z)
9090
@test res2_wei ≈ res2
9191
end
@@ -96,14 +96,14 @@ if dim2
9696

9797
z_gpu = CuArray(z)
9898

99-
z2_gpu = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z_gpu, Int[]; rdft)
99+
z2_gpu = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z_gpu, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
100100
@test z2_gpu ≈ z_gpu
101101

102-
res1_gpu = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), reshape(z_gpu, (N1, N2)); rdft)
102+
res1_gpu = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), reshape(z_gpu, (N1, N2)), nlp.fft_timer, nlp.mapping_timer; rdft)
103103
@test norm(res1_gpu) ≈ norm(z_gpu)
104104
@test res1_wei ≈ collect(res1_gpu)
105105

106-
res2_gpu = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z_gpu, Int[]; rdft)
106+
res2_gpu = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 2, (N1, N2), z_gpu, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
107107
@test norm(res2_gpu) ≈ norm(z_gpu)
108108
@test res2_wei ≈ collect(res2_gpu)
109109
end
@@ -130,14 +130,14 @@ if dim3
130130
@testset "3D -- CPU -- rdft=$rdft -- $N1 × $N2 × $N3" begin
131131
nlp, solver, results = fft_example_3D(N1, N2, N3; gpu=false, rdft, check=true)
132132

133-
z2 = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z, Int[]; rdft)
133+
z2 = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
134134
@test z2 ≈ z
135135

136-
res1 = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), reshape(z, (N1, N2, N3)); rdft)
136+
res1 = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), reshape(z, (N1, N2, N3)), nlp.fft_timer, nlp.mapping_timer; rdft)
137137
@test norm(res1) ≈ norm(z)
138138
@test res1_wei ≈ res1
139139

140-
res2 = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z, Int[]; rdft)
140+
res2 = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
141141
@test norm(res2) ≈ norm(z)
142142
@test res2_wei ≈ res2
143143
end
@@ -148,14 +148,14 @@ if dim3
148148

149149
z_gpu = CuArray(z)
150150

151-
z2_gpu = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z_gpu, Int[]; rdft)
151+
z2_gpu = M_perpt_M_perp_vec(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z_gpu, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
152152
@test z2_gpu ≈ z_gpu
153153

154-
res1_gpu = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), reshape(z_gpu, (N1, N2, N3)); rdft)
154+
res1_gpu = M_perp_tz(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), reshape(z_gpu, (N1, N2, N3)), nlp.fft_timer, nlp.mapping_timer; rdft)
155155
@test norm(res1_gpu) ≈ norm(z_gpu)
156156
@test res1_wei ≈ collect(res1_gpu)
157157

158-
res2_gpu = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z_gpu, Int[]; rdft)
158+
res2_gpu = M_perp_beta(nlp.buffer_real, nlp.buffer_complex1, nlp.buffer_complex2, nlp.op, 3, (N1, N2, N3), z_gpu, Int[], nlp.fft_timer, nlp.mapping_timer; rdft)
159159
@test norm(res2_gpu) ≈ norm(z_gpu)
160160
@test res2_wei ≈ collect(res2_gpu)
161161
end

0 commit comments

Comments
 (0)