From fb12a13702d19f5809933e85a13f63333c6b5a99 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Wed, 8 Sep 2021 12:20:55 +0200
Subject: [PATCH] finalize primal-dual comparison experiment
---
scripts/run_experiments.jl | 121 ++++++++++++++++++++++---------------
scripts/util.jl | 11 ++--
2 files changed, 77 insertions(+), 55 deletions(-)
diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index d1064c2..de1f869 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -410,6 +410,9 @@ function estimate!(ctx::L1L2TVContext)
nablau = nabla(ctx.u), w, nablaw = nabla(w), ctx.tdata)
end
+estimate_error(st::L1L2TVContext) =
+ sqrt(sum(st.est.data) / area(st.mesh))
+
# minimal Dörfler marking
function mark(ctx::L1L2TVContext; theta=0.5)
n = ncells(ctx.mesh)
@@ -529,113 +532,132 @@ function denoise(img; name, params...)
end
-function denoise_pd(ctx, img; df=nothing, name, algorithm, params_...)
+function denoise_pd(st, img; df=nothing, name, algorithm, params_...)
params = NamedTuple(params_)
m = 1
img = from_img(img) # coord flip
- mesh = init_grid(img; type=:vertex)
+ mesh = init_grid(img)
#mesh = init_grid(img, 5, 5)
T(tdata, u) = u
S(u, nablau) = u
- ctx = L1L2TVContext(name, mesh, m;
+ st = L1L2TVContext(name, mesh, m;
T, tdata = nothing, S,
params.alpha1, params.alpha2, params.lambda, params.beta,
params.gamma1, params.gamma2)
# semi-implicit primal dual parameters
- gamma = ctx.alpha2 + ctx.beta # T = I, S = I
- gamma /= 100 # kind of arbitrary?
+ mu = st.alpha2 + st.beta # T = I, S = I
+ #mu /= 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))
- ctx.u.data .= ctx.g.data
+ #project_img!(st.g, img)
+ interpolate!(st.g, x -> interpolate_bilinear(img, x))
+ st.u.data .= st.g.data
- save_denoise(ctx, i) =
- output(ctx, "output/$(ctx.name)_$(lpad(i, 5, '0')).vtu",
- ctx.g, ctx.u, ctx.p1, ctx.p2)
+ save_denoise(st, i) =
+ output(st, "output/$(st.name)_$(lpad(i, 5, '0')).vtu",
+ st.g, st.u, st.p1, st.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)))
+ function log!(df::DataFrame; kwargs...)
+ row = NamedTuple(kwargs)
+ push!(df, row)
+ #println(df)
+ println(row)
end
- #pvd = paraview_collection("output/$(ctx.name).pvd")
- #pvd[0] = save_denoise(ctx, 0)
+ #pvd = paraview_collection("output/$(st.name).pvd")
+ #pvd[0] = save_denoise(st, 0)
k = 0
- println("primal energy: $(primal_energy(ctx))")
+ println("primal energy: $(primal_energy(st))")
while true
k += 1
if algorithm == :pd1
# no step size control
- step_pd2!(ctx; sigma, tau, theta)
+ step_pd2!(st; sigma, tau, theta)
elseif algorithm == :pd2
- theta = 1 / sqrt(1 + 2 * gamma * tau)
+ theta = 1 / sqrt(1 + 2 * mu * tau)
tau *= theta
sigma /= theta
- step_pd2!(ctx; sigma, tau, theta)
+ step_pd2!(st; sigma, tau, theta)
elseif algorithm == :newton
- step!(ctx)
+ step!(st)
end
- #pvd[k] = save_denoise(ctx, k)
+ #pvd[k] = save_denoise(st, k)
domain_factor = 1 / sqrt(area(mesh))
- norm_step_ = norm_step(ctx) * domain_factor
- norm_residual_ = norm_residual(ctx) * domain_factor
+ norm_step_ = norm_step(st) * domain_factor
+ #estimate!(st)
- log!(df; k, norm_step = norm_step_, norm_residual = norm_residual_)
+ log!(df; k,
+ norm_step = norm_step_,
+ #est = estimate_error(st),
+ primal_energy = primal_energy(st))
#norm_residual_ < params.tol && norm_step_ < params.tol && break
haskey(params, :tol) && norm_step_ < params.tol && break
haskey(params, :max_iters) && k >= params.max_iters && break
end
- #pvd[1] = save_denoise(ctx, 1)
+ #pvd[1] = save_denoise(st, 1)
#vtk_save(pvd)
- return ctx
+ return st
end
function experiment_pd_comparison(ctx)
img = loadimg(joinpath(ctx.indir, "input.png"))
- img = from_img(img) # coord flip
+ #img = [0. 0.; 1. 0.]
algparams = (
alpha1=0., alpha2=30., lambda=1., beta=0.,
- gamma1=1e-3, gamma2=1e-3,
- tol = 1e-6, max_iters = 50,
+ gamma1=1e-2, gamma2=1e-3,
+ tol = 1e-10,
+ max_iters = 10000,
)
df1 = DataFrame()
df2 = DataFrame()
df3 = DataFrame()
- denoise_pd(ctx, img; name="test", algorithm=:pd1, df = df1, algparams...);
- denoise_pd(ctx, img; name="test", algorithm=:pd2, df = df2, algparams...);
- denoise_pd(ctx, img; name="test", algorithm=:newton, df = df3, algparams...);
+ st1 = denoise_pd(ctx, img; name="test",
+ algorithm=:pd1, df = df1, algparams...);
+ st2 = denoise_pd(ctx, img; name="test",
+ algorithm=:pd2, df = df2, algparams...);
+ st3 = denoise_pd(ctx, 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(ctx.outdir, "semiimplicit.csv"), logfilter(df1))
- CSV.write(joinpath(ctx.outdir, "semiimplicit-accelerated.csv"), logfilter(df2))
- CSV.write(joinpath(ctx.outdir, "newton.csv"), logfilter(df3))
+ df1[!, :energy_min] .= energy_min
+ df2[!, :energy_min] .= energy_min
+ df3[!, :energy_min] .= energy_min
+
+ CSV.write(joinpath(ctx.outdir, "semi-implicit.csv"),
+ logfilter(df1))
+ CSV.write(joinpath(ctx.outdir, "semi-implicit-accelerated.csv"),
+ logfilter(df2))
+ CSV.write(joinpath(ctx.outdir, "newton.csv"),
+ logfilter(df3))
+
+ saveimg(joinpath(ctx.outdir, "input.png"),
+ img)
+ saveimg(joinpath(ctx.outdir, "semi-implicit.png"),
+ to_img(sample(st1.u)))
+ saveimg(joinpath(ctx.outdir, "semi-implicit-accelerated.png"),
+ to_img(sample(st2.u)))
+ saveimg(joinpath(ctx.outdir, "newton.png"),
+ to_img(sample(st3.u)))
+
+ savedata(joinpath(ctx.outdir, "data.tex");
+ energy_min, algparams...)
end
function denoise_approximation(ctx)
@@ -662,7 +684,8 @@ function denoise_approximation(ctx)
function log!(df::DataFrame; kwargs...)
row = NamedTuple(kwargs)
push!(df, row)
- println(row)
+ println(df)
+ #println(row)
end
pvd_path = joinpath(ctx.outdir, "$(ctx.params.name).pvd")
@@ -716,9 +739,9 @@ function experiment_approximation(ctx)
denoise_approximation(Util.Context(ctx; name = "test", df,
img, mesh,
- #alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
- alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
- gamma1 = 1e-5, gamma2 = 1e-3,
+ alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
+ #alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
+ gamma1 = 1e-5, gamma2 = 1e-5,
tol = 1e-10, adaptive = true,
))
end
diff --git a/scripts/util.jl b/scripts/util.jl
index f57f0d7..5eff824 100644
--- a/scripts/util.jl
+++ b/scripts/util.jl
@@ -44,13 +44,12 @@ end
# LaTeX Data Output
-# TODO: don't depend on ctx.path and use filename only
-savedata(ctx::Context, filename; data...) =
- open(joinpath(ctx.outdir, filename); write=true) do io
- _savedata(io, ctx.path; data...)
+savedata(filename; data...) =
+ open(filename; write=true) do io
+ _savedata(io; data...)
end
-function _savedata(io, path; data...)
+function _savedata(io; data...)
isvalidname(x) = contains(string(x), r"^[^/\\]+$")
# define tex command to access data
@@ -67,7 +66,7 @@ function _savedata(io, path; data...)
# actual data encoded as tex-commands
write(io, "\\expandafter\\def")
- write(io, "\\csname $(texcmd)/$(path)/$(key)\\endcsname{")
+ write(io, "\\csname $(texcmd)/\\DataPrefix/$(key)\\endcsname{")
show(io, MIME("text/latex"), value)
write(io, "}\n")
end
--
GitLab