From 973748987533e7491efff3ea6aa6067b5e5fc632 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Mon, 2 Aug 2021 19:02:19 +0200
Subject: [PATCH] add first experiment

---
 Manifest.toml |   6 +++
 Project.toml  |   1 +
 src/image.jl  |  17 ++++++-
 src/run.jl    | 127 +++++++++++++++++++++++++++++++++++---------------
 4 files changed, 113 insertions(+), 38 deletions(-)

diff --git a/Manifest.toml b/Manifest.toml
index cd18d24..af12a4e 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -33,6 +33,12 @@ git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 version = "0.11.0"
 
+[[Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
+git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.12.8"
+
 [[CommonSubexpressions]]
 deps = ["MacroTools", "Test"]
 git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
diff --git a/Project.toml b/Project.toml
index d7d50c6..7aef229 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/src/image.jl b/src/image.jl
index ab14b37..a8bebd6 100644
--- a/src/image.jl
+++ b/src/image.jl
@@ -1,4 +1,4 @@
-export interpolate_bilinear
+export interpolate_bilinear, halve
 
 eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...]
 
@@ -28,3 +28,18 @@ function warp_backwards(img, u)
     end
     return res
 end
+
+function halve(img)
+    all(iseven.(size(img))) || throw(ArgumentError("uneven size"))
+
+    res = similar(img, size(img) .÷ 2)
+    fill!(res, zero(eltype(res)))
+
+    for i in CartesianIndices(img)
+        j = CartesianIndex((Tuple(i) .+ 1) .÷ 2)
+        res[j] += img[i]
+    end
+    res ./= 2 ^ ndims(img)
+
+    return res
+end
diff --git a/src/run.jl b/src/run.jl
index ad33d2b..2478e3e 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -1,6 +1,7 @@
 export myrun, denoise, denoise_pd, inpaint, optflow, solve_primal!, estimate!, loadimg, saveimg
 
 using LinearAlgebra: norm
+using Colors: Gray
 
 # avoid world-age-issues by preloading ColorTypes
 import ColorTypes
@@ -43,7 +44,7 @@ function L1L2TVContext(name, mesh, m; T, tdata, S,
     Vest = FeSpace(mesh, DP0(), (1,))
     Vg = FeSpace(mesh, P1(), (1,))
     Vu = FeSpace(mesh, P1(), (m,))
-    Vp1 = FeSpace(mesh, P1(), (1,))
+    Vp1 = FeSpace(mesh, DP0(), (1,))
     Vp2 = FeSpace(mesh, DP1(), (m, d))
 
     est = FeFunction(Vest, name="est")
@@ -145,13 +146,15 @@ function step!(ctx::L1L2TVContext)
 	ctx.u, ctx.nablau, ctx.p2, ctx.du, ctx.nabladu)
 
     # newton update
-    ctx.u.data .+= ctx.du.data
-    ctx.p1.data .+= ctx.dp1.data
-    ctx.p2.data .+= ctx.dp2.data
+    theta = 1.
+    ctx.u.data .+= theta * ctx.du.data
+    ctx.p1.data .+= theta * ctx.dp1.data
+    ctx.p2.data .+= theta * ctx.dp2.data
 
     # reproject p1
     function p1_project!(p1, alpha1)
 	p1.space.element == DP0() || p1.space.element == P1() ||
+	    p1.space.element == DP1() ||
 	    throw(ArgumentError("element unsupported"))
 	p1.data .= clamp.(p1.data, -alpha1, alpha1)
     end
@@ -286,6 +289,7 @@ function step_pd2!(ctx::L1L2TVContext; sigma, tau, theta = 1.)
     # reproject p1
     function p1_project!(p1, alpha1)
 	p1.space.element == DP0() || p1.space.element == P1() ||
+	    p1.space.element == DP1() ||
 	    throw(ArgumentError("element unsupported"))
 	p1.data .= clamp.(p1.data, -alpha1, alpha1)
     end
@@ -365,11 +369,11 @@ function estimate!(ctx::L1L2TVContext)
     # FIXME: sign?
     function estf(x_; g, u, p1, p2, nablau, w, nablaw, tdata)
 	alpha1part = iszero(ctx.alpha1) ? 0. : ctx.alpha1 * (
-	    huber(norm(ctx.T(tdata, u) - g), ctx.gamma1) +
+	    huber(norm(ctx.T(tdata, u) - g), ctx.gamma1) -
 	    dot(ctx.T(tdata, u) - g, p1 / ctx.alpha1) +
 	    ctx.gamma1 / 2 * norm(p1 / ctx.alpha1)^2)
 	lambdapart = iszero(ctx.lambda) ? 0. : ctx.lambda * (
-	    huber(norm(nablau), ctx.gamma2) +
+	    huber(norm(nablau), ctx.gamma2) -
 	    dot(nablau, p2 / ctx.lambda) +
 	    ctx.gamma2 / 2 * norm(p2 / ctx.lambda)^2)
 	bpart = 1 / 2 * (
@@ -448,7 +452,7 @@ end
 
 toimg(arr) = Gray.(clamp.(arr, 0., 1.))
 loadimg(x) = reverse(transpose(Float64.(FileIO.load(x))); dims = 2)
-saveimg(io, x) = FileIO.save(io, toimg(x))
+saveimg(io, x) = FileIO.save(io, transpose(reverse(toimg(x); dims = 2)))
 
 function primal_energy(ctx::L1L2TVContext)
     function integrand(x; g, u, nablau, tdata)
@@ -469,7 +473,18 @@ function norm_residual(ctx::L1L2TVContext)
     w = FeFunction(ctx.u.space)
     solve_primal!(w, ctx)
     w.data .-= ctx.u.data
-    return norm_l2(w)
+    upart2 = norm_l2(w)^2
+
+    function integrand(x; g, u, nablau, p1, p2, tdata)
+	p1part = p1 * max(ctx.gamma1, norm(ctx.T(tdata, u) - g)) -
+	    ctx.alpha1 * (ctx.T(tdata, u) - g)
+	p2part = p2 * max(ctx.gamma2, norm(nablau)) -
+	    ctx.lambda * nablau
+	return norm(p1part)^2 + norm(p2part)^2
+    end
+    ppart2 = integrate(ctx.mesh, integrand; ctx.g, ctx.u, ctx.nablau, ctx.p1, ctx.p2, ctx.tdata)
+
+    return sqrt(upart2 + ppart2)
 end
 
 function denoise(img; name, params...)
@@ -527,7 +542,8 @@ function denoise(img; name, params...)
     return ctx
 end
 
-function denoise_pd(img; name, params...)
+
+function denoise_pd(img; df=nothing, name, algorithm, params...)
     m = 1
     mesh = init_grid(img; type=:vertex)
     #mesh = init_grid(img, 5, 5)
@@ -539,19 +555,30 @@ function denoise_pd(img; name, params...)
 
     # semi-implicit primal dual parameters
     gamma = ctx.alpha2 + ctx.beta # T = I, S = I
-    sigma = 1e-2
-    tau = 1e-2
+    gamma /= 100 # kind of arbitrary?
+
+    tau = 1e-1
+    L = 100
+    sigma = inv(tau * L^2)
     theta = 1.
 
-    project_img!(ctx.g, img)
-    #interpolate!(ctx.g, x -> interpolate_bilinear(img, x))
-    #m = (size(img) .- 1) ./ 2 .+ 1
-    #interpolate!(ctx.g, x -> norm(x .- m) < norm(m .- 1) / 3)
+    #project_img!(ctx.g, img)
+    interpolate!(ctx.g, x -> interpolate_bilinear(img, x))
     ctx.u.data .= ctx.g.data
 
     save_denoise(ctx, i) =
 	output(ctx, "output/$(ctx.name)_$(lpad(i, 5, '0')).vtu",
-	    ctx.g, ctx.u, ctx.p1, ctx.p2, ctx.est, ctx.du)
+	    ctx.g, ctx.u, ctx.p1, ctx.p2)
+
+    log!(x::Nothing; kwargs...) = x
+    function log!(df::DataFrame; k, norm_step, norm_residual)
+	push!(df, (;
+	    k,
+	    primal_energy = primal_energy(ctx),
+	    norm_step,
+	    norm_residual))
+	println(NamedTuple(last(df)))
+    end
 
     pvd = paraview_collection("output/$(ctx.name).pvd")
     pvd[0] = save_denoise(ctx, 0)
@@ -559,44 +586,70 @@ function denoise_pd(img; name, params...)
     k = 0
     println("primal energy: $(primal_energy(ctx))")
 
-    algorithm = :pd2
     while true
 	k += 1
-	#println("")
-	if algorithm == :pd2 #|| primal_energy(ctx) < 2409
-	    println("step_pd: theta = $theta, tau = $tau, sigma = $sigma")
-	    #theta = 1 / sqrt(1 + 2 * gamma * tau)
-	    #tau *= theta
-	    #sigma /= theta
-	    step_pd2!(ctx; sigma, tau, theta)
-	elseif algorithm == :pd1
+	if algorithm == :pd1
 	    # no step size control
-	    step_pd!(ctx; sigma, tau, theta)
+	    step_pd2!(ctx; sigma, tau, theta)
+	elseif algorithm == :pd2
+	    theta = 1 / sqrt(1 + 2 * gamma * tau)
+	    tau *= theta
+	    sigma /= theta
+	    step_pd2!(ctx; sigma, tau, theta)
 	elseif algorithm == :newton
 	    step!(ctx)
 	end
-	#estimate!(ctx)
 	#pvd[k] = save_denoise(ctx, k)
-	println()
 
 	domain_factor = 1 / sqrt(area(mesh))
 	norm_step_ = norm_step(ctx) * domain_factor
-	residual = -1.
-	residual = norm_residual(ctx) * domain_factor
+	norm_residual_ = norm_residual(ctx) * domain_factor
 
-	#println("ndofs: $(ndofs(ctx.u.space)), est: $(norm_l2(ctx.est)))")
-	println("primal energy: $(primal_energy(ctx))")
-	println("norm_step: $(norm_step_)")
-	println("norm_residual: $(residual)")
+	log!(df; k, norm_step = norm_step_, norm_residual = norm_residual_)
 
-	#return ctx
-	#residual <= 1e-3 && break
-	#k >= 5 && break
+	norm_residual_ < 1e-6 && norm_step_ < 1e-6 && break
+	norm_step_ < 1e-13 && break
     end
+    pvd[1] = save_denoise(ctx, 1)
     vtk_save(pvd)
     return ctx
 end
 
+"""
+    logfilter(dt; a=20)
+
+Reduce dataset such that the resulting dataset would have approximately
+equidistant datapoints on a log-scale. `a` controls density.
+"""
+logfilter(dt; a=20) = filter(:k => k -> round(a*log(k+1)) - round(a*log(k)) > 0, dt)
+
+function experiment1(img)
+    path = "data/pd-comparison"
+
+    img = loadimg("data/denoising/input.png")
+
+    algparams = (alpha1=0., alpha2=30., lambda=1., beta=0., gamma1=1e-3, gamma2=1e-3)
+
+    df1 = DataFrame()
+    df2 = DataFrame()
+    df3 = DataFrame()
+
+    denoise_pd(img; name="test", algorithm=:pd1, df = df1, algparams...);
+    denoise_pd(img; name="test", algorithm=:pd2, df = df2, algparams...);
+    denoise_pd(img; name="test", algorithm=:newton, df = df3, algparams...);
+
+    energy_min = min(minimum(df1.primal_energy), minimum(df2.primal_energy),
+		     minimum(df3.primal_energy))
+
+    df1.primal_energy .-= energy_min
+    df2.primal_energy .-= energy_min
+    df3.primal_energy .-= energy_min
+
+    CSV.write(joinpath(path, "semiimplicit.csv"), logfilter(df1))
+    CSV.write(joinpath(path, "semiimplicit-accelerated.csv"), logfilter(df2))
+    CSV.write(joinpath(path, "newton.csv"), logfilter(df3))
+end
+
 function inpaint(img, imgmask; name, params...)
     size(img) == size(imgmask) ||
 	throw(ArgumentError("non-matching dimensions"))
-- 
GitLab