code: working chambolle dd

FD/src/FD.jl

@@ -225,11 +225,19 @@ grid(f::GridFunction) = f.grid
 Base.view(f::AbstractGridFunction{G}, idx...) where G<:MetaGrid =
     GridFunction(GridView(f.grid, indices(f.grid, idx...)), view(f.data, indices(f.grid, idx...)..., :))
 
-Base.parent(f::AbstractGridFunction{<:GridView}) =
+function zero_continuation(f::GridFunction{<:GridView})
+  dst = GridFunction(parent(grid(f)), zeros(size(parent(grid(f)))..., rangedim(f)), Val(rangedim(f)))
+  dst.data[grid(f).indices..., :] .= f.data
+  dst
+end
+
+Base.parent(f::GridFunction{<:GridView,<:Any,<:Any,<:Any,<:SubArray}) =
     GridFunction(parent(grid(f)), parent(f.data))
 
 
 
+
+
 function _check_shape(a::GridFunction, b::GridFunction)
   #TODO: why does fv.grid == gv.grid fail?
   a.grid === b.grid || throw(ArgumentError("grids need to match"))

chambollefp_dd.jl

@@ -36,35 +36,14 @@ end
 
 # TODO: investigate matrix multiplication performance with GridFunctions
 function fpiter(grid, p, wg, A, α, β, τ)
-  ## local variant
-  #w = FD.GridFunction(grid, view(wg.data, grid.indices..., :))
-  #x = FD.GridFunction(grid, (A'*A + β*I) \ (FD.divergence(p) + w).data[:])
-  #q = FD.gradient(x)
-
-  ## global variant
-  x = FD.GridFunction(FD.parent(grid), (A'*A + β*I) \ (FD.divergence(p) + wg).data[:])
-  #FD.my_plot(FD.divergence(p) + wg)
-
+  x = FD.GridFunction(FD.parent(grid), (A'*A + β*I) \ (FD.divergence_z(FD.zero_continuation(p)) + wg).data[:])
   q = FD.GridFunction(grid, view(FD.gradient(x).data, grid.indices..., :))
-  q2 = FD.gradient(FD.GridFunction(grid, view(x.data, grid.indices..., :)))
-
-  #display(x.data[10:13,51:54,1])
-  #display(q.data[10:13,52:53,2])
-  #display(q2.data[10:13,52:53,2])
-  #sleep(0.1)
-
+  ## this misses out on the right boundary
   #q = FD.gradient(FD.GridFunction(grid, view(x.data, grid.indices..., :)))
 
-
-  #FD.my_plot(FD.pw_norm(q))
-  #FD.my_plot(x)
-  #FD.my_plot(α)
-
   p_new = FD.GridFunction(grid, FD.pwmul(FD.GridFunction(grid, α ./ (α .+ τ .* FD.pw_norm(q))),
                                          FD.GridFunction(grid, p .+ τ .* q)))
 
-
-
   reserr = maximum(sum((p_new .- p).^2, dims=ndims(p)))
   (p_new, reserr)
 end
@@ -72,113 +51,47 @@ end
 
 function fpalg(grid, w, A, p = FD.GridFunction(x -> 0, ndims(grid), grid);
                α, β, τ = 1/8, ϵ = 1e-4)
-  local uimg
   reserr = 1
   k = 0
   while reserr > ϵ && k < 10e1
     k += 1
     (p, reserr) = fpiter(grid, p, w, A, α, β, τ)
-
-    #u = p2u(grid, p, A, f, α=α, β=β)
-    #@show k
-    #FD.my_plot(u)
   end
   return p
 end
 
 
 function dditer(grid, p, A, f; α, β, τ, ϵ, ρ)
-  #pzero = FD.GridFunction(x -> 0, ndims(grid), grid)
-
   p̂ = (1-ρ) * p
 
   i = 0
   for part in FD.parts(grid)
     gridv = view(grid, Tuple(part)...)
+
     fv = view(f, Tuple(part)...)
-    #p.data .= 0
+
     p̂v = view(p̂, Tuple(part)...)
     pv = view(p, Tuple(part)...)
-    #pzerov = view(pzero, Tuple(part)...)
 
     θ = theta(grid, Tuple(part)...)
     θv = view(θ, Tuple(part)...)
 
-    #y = FD.divergence(FD.pwmul(FD.GridFunction(grid, 1 .- θ), p))
-    #FD.my_plot(y)
-    #sleep(1)
-
     w = FD.GridFunction(grid, A' * f.data[:]) + FD.divergence(FD.pwmul(FD.GridFunction(grid, 1 .- θ), p))
     wv = view(w, Tuple(part)...)
 
-    #FD.my_plot(w)
-    #sleep(1)
-
-    #wv = A' * fv
-
-    #u = p2u(grid, p̂, A, f, α=α, β=β)
-    #uv = view(u, Tuple(part)...)
-
-    #FD.my_plot(uv)
-
-    #sleep(1)
-
     p̂v .+= ρ .* fpalg(gridv, w, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
-    #p̂v .= fpalg(gridv, wv, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
-
-  #FD.my_plot(FD.pw_norm(p))
-  #sleep(1)
-  #FD.my_plot(FD.pw_norm(f))
-  #sleep(1)
-  #FD.my_plot(FD.pw_norm(p̂))
-  #display(p̂.data[10:13,52:54,:])
-    #break
-    #FD.my_plot(FD.pw_norm(p̂))
-    #sleep(1)
-    #FD.my_plot(θ)
-    #display(p̂v.data)
-    #break
-    #sleep(1)
-
-    #u = p2u(grid, p, A, f, α=α*θ, β=β)
-    u = p2u(grid, p̂, A, f, α=α, β=β)
-    uv = view(u, Tuple(part)...)
-
-    #FD.my_plot(θv)
-
-    #sleep(1)
-    #FD.my_plot(FD.pw_norm(p̂v))
-    #@show maximum(u), minimum(u)
-    #@show maximum(p̂v), minimum(p̂v)
-    #sleep(1)
-    #FD.my_plot(FD.pw_norm(p̂))
-    #FD.my_plot(u)
-    #FD.my_plot(p̂v)
-    #y = FD.GridFunction(gridv, FD.gradient(θv).data[:,:,1])
-    #display(y.data)
-    #display(θv.data)
-    #FD.my_plot(FD.GridFunction(gridv, FD.gradient(θv).data[:,:,2]))
-
-    #display(FD.gradient(θv).data)
-    #break
-
-    #@show size(θv.data), size(p̂v.data)
-
-    #FD.my_plot(FD.pw_norm(p̂))
 
     i += 1
-    #break
-    #i > 3 && break
   end
   p .= p̂
 
   u = p2u(grid, p, A, f, α=α, β=β)
-  #FD.my_plot(u)
-  FD.my_plot(FD.pw_norm(p))
+  FD.my_plot(u)
+  #FD.my_plot(FD.pw_norm(p))
 end
 
 
-grid = FD.MetaGrid((0:0.01:1, 0:0.01:1), (1, 2), 5)
+grid = FD.MetaGrid((0:0.01:1, 0:0.01:1), (2, 2), 5)
 
 A = I
 #A = spdiagm(0 => ones(length(grid)), 1 => ones(length(grid)-1) )
@@ -195,7 +108,7 @@ g = FD.GridFunction(x -> norm(x), 2, grid)
 
 p = FD.GridFunction(x -> 0, 2, grid)
 
-for n = 1:1
+for n = 1:10
   println("iteration $n")
   dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
 end