refactor

DualTVDD/src/DualTVDD.jl

@@ -7,6 +7,7 @@ include("fd.jl")
 
 using LinearAlgebra: I
 using .FD
+using .FD: parts
 
 # export
 
@@ -60,7 +61,10 @@ Fixed point iteration (Chambolle).
 """
 # TODO: investigate matrix multiplication performance with GridFunctions
 function fpiter(grid, p, wg, A, α, β, τ)
-    x = FD.GridFunction(FD.parent(grid), (A'*A + β*I) \ (FD.divergence_z(FD.zero_continuation(p)) + wg).data[:])
+    # applying B is a global operation
+    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..., :))
     ## this misses out on the right boundary
     #q = FD.gradient(FD.GridFunction(grid, view(x.data, grid.indices..., :)))
@@ -90,38 +94,73 @@ end
 
 
 """
-dditer(grid, p, A, f; α, β, τ, ϵ, ρ)
+    dditer(metagrid, p, A, f; α, β, τ, ϵ, ρ)
 
 Domain decomposition algorithm, calling `fpalg` for each part.
 """
-function dditer(grid, p, A, f; α, β, τ, ϵ, ρ)
+function dditer(metagrid, p, A, f; α, β, τ, ϵ, ρ)
+    # contribution from previous step
     p̂ = (1-ρ) * p
 
-    i = 0
-    for part in FD.parts(grid)
-        gridv = view(grid, Tuple(part)...)
+    for part in parts(metagrid)
+        grid = view(metagrid, Tuple(part)...)
 
         fv = view(f, Tuple(part)...)
 
         p̂v = view(p̂, Tuple(part)...)
         pv = view(p, Tuple(part)...)
 
-        θ = partition_function(grid, Tuple(part)...)
+        θ = partition_function(metagrid, Tuple(part)...)
         θv = view(θ, Tuple(part)...)
 
-        w = FD.GridFunction(grid, A' * f.data[:]) + FD.divergence(FD.pwmul(FD.GridFunction(grid, 1 .- θ), p))
+        w = FD.GridFunction(metagrid, A' * f.data[:]) +
+            FD.divergence(FD.pwmul(FD.GridFunction(metagrid, 1 .- θ), p))
         wv = view(w, Tuple(part)...)
 
-        p̂v .+= ρ .* fpalg(gridv, w, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
-
-        i += 1
+        p̂v .+= ρ .* fpalg(grid, w, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
     end
     p .= p̂
 
-    u = p2u(grid, p, A, f, α=α, β=β)
+    # reconstruct primal solution and plot (FIXME: hack)
+    u = p2u(metagrid, p, A, f, α=α, β=β)
     FD.my_plot(u)
     #FD.my_plot(FD.pw_norm(p))
 end
 
+using LinearAlgebra: norm
+using SparseArrays: spdiagm
+"""
+    ddalg()
+
+temporary run function
+"""
+function ddalg()
+    grid = MetaGrid((0:0.01:1, 0:0.01:1), (1, 1), 5)
+
+    f = GridFunction(x -> float(norm(x .- [0.5, 0.5]) < 0.4), 1, grid)
+    g = GridFunction(x -> norm(x), 2, grid)
+
+    A = I
+    #A = spdiagm(0 => ones(length(grid)), 1 => 0.3 * ones(length(grid)-1))
+    #f = GridFunction(grid, A * f.data[:])
+
+    "TV parameter"
+    α=1
+    "L2 parameter"
+    β=0
+    "chambolle weighing"
+    τ=1/8
+    "residual error tolerance for local problem"
+    ϵ=1e-9
+    "???"
+    ρ=0.5
+
+    p = GridFunction(x -> 0, 2, grid)
+
+    for n = 1:10
+      println("iteration $n")
+      dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
+    end
+end
 
 end # module

DualTVDD/src/fd.jl

@@ -25,6 +25,7 @@ using .FDiff
 # export
 export Grid, MetaGrid, GridView, GridFunction
 export grid
+export ⊙
 
 
 # abstract types
@@ -71,8 +72,8 @@ Grid(domain::FDStepRange...) = Grid(domain)
 
 # TODO: somehow display(grid) is incorrect for d = 1 ?
 Base.show(io::IO, grid::Grid) =
-    print("Grid from $(first.(grid.domain)) to $(last.(grid.domain)) ",
-          "with spacing $(step.(grid.domain))")
+    print(io, "Grid from $(first.(grid.domain)) to $(last.(grid.domain)) ",
+              "with spacing $(step.(grid.domain))")
 
 Base.size(grid::Grid) = ntuple(i -> grid.domain[i].len, ndims(grid))
 
@@ -207,8 +208,8 @@ function GridFunction(grid::AbstractGrid, ::UndefInitializer, ::Val{n} = Val(1))
     GridFunction{typeof(grid), eltype(data), ndims(grid), n, typeof(data), ndims(grid) + 1}(grid, data)
 end
 
-#Base.show(io::IO, grid::GridFunction) =
-#    print("GridFunction on $(length(grid)) grid points")
+Base.show(io::IO, f::GridFunction) =
+    print(io, "$(rangedim(f))d-GridFunction on $(f.grid), $(length(f)) dofs")
 
 Base.size(f::GridFunction) = (length(f.data),)
 #Base.axes(f::GridFunction) = (Base.OneTo(length(f.data)),)
@@ -390,6 +391,7 @@ end
 Pointwise multiplication of two gridfunctions.
 """
 ⊙(f, g) = pwmul(f, g)
+pwmul(f::GridFunction{G,T,d,n}, g::GridFunction{G,T,d,n}) where {G,T,d,n} = f .* g
 function pwmul(f::GridFunction{G,T,d,1}, g::GridFunction{G,T,d}) where {G,T,d}
     # TODO: maybe comprehension? Or make broadcast work
     data = similar(g.data)