From 6e6603854c5ceaf9a12c887991f9c1f3b335e5e8 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan.hilb@mathematik.uni-stuttgart.de>
Date: Mon, 9 Sep 2019 13:39:56 +0200
Subject: [PATCH] refactor
---
DualTVDD/src/DualTVDD.jl | 63 ++++++++++++++++++++++++++++++++--------
DualTVDD/src/fd.jl | 10 ++++---
2 files changed, 57 insertions(+), 16 deletions(-)
diff --git a/DualTVDD/src/DualTVDD.jl b/DualTVDD/src/DualTVDD.jl
index 1150f71..b28348c 100644
--- a/DualTVDD/src/DualTVDD.jl
+++ b/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
diff --git a/DualTVDD/src/fd.jl b/DualTVDD/src/fd.jl
index 6005035..6829c9e 100644
--- a/DualTVDD/src/fd.jl
+++ b/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)
--
GitLab