From 6100586a53b7e6dccbeb74d9dd657f47991d7e64 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Wed, 9 Feb 2022 17:02:32 +0100
Subject: [PATCH] add convergence experiments
---
scripts/run_experiments.jl | 130 +++++++++++++++++++++++++++----------
1 file changed, 97 insertions(+), 33 deletions(-)
diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 170bed8..954c193 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -794,10 +794,14 @@ function denoise_pd(ctx)
#project_l2_lagrange!(st.g, ctx.params.g_arr)
interpolate!(st.g, x -> evaluate_bilinear(ctx.params.g_arr, x))
- st.u.data .= st.g.data
+ if haskey(ctx.params, :init_zero) && ctx.params.init_zero
+ st.u.data .= 0
+ else
+ st.u.data .= st.g.data
+ end
- save_denoise(st, i) =
- output(st, "output/$(st.name)_$(lpad(i, 5, '0')).vtu",
+ save_step(st, i) =
+ output(st, joinpath(ctx.outdir, "output_$(ctx.params.algorithm)_$(lpad(i, 5, '0')).vtu"),
st.g, st.u, st.p1, st.p2)
log!(x::Nothing; kwargs...) = x
@@ -808,42 +812,41 @@ function denoise_pd(ctx)
println(row)
end
- #pvd = paraview_collection("output/$(st.name).pvd")
- #pvd[0] = save_denoise(st, 0)
+ pvd = paraview_collection(joinpath(ctx.outdir, "output_$(ctx.params.algorithm).pvd")) do pvd
+ k = 0
+ println("primal energy: $(primal_energy(st))")
- k = 0
- println("primal energy: $(primal_energy(st))")
+ while true
+ k += 1
+ if ctx.params.algorithm == :pd1
+ # no step size control
+ step_pd2!(st; sigma, tau, theta)
+ elseif ctx.params.algorithm == :pd2
+ theta = 1 / sqrt(1 + 2 * mu * tau)
+ tau *= theta
+ sigma /= theta
+ step_pd2!(st; sigma, tau, theta)
+ elseif ctx.params.algorithm == :newton
+ step!(st)
+ end
+ #pvd[k] = save_step(st, k)
- while true
- k += 1
- if ctx.params.algorithm == :pd1
- # no step size control
- step_pd2!(st; sigma, tau, theta)
- elseif ctx.params.algorithm == :pd2
- theta = 1 / sqrt(1 + 2 * mu * tau)
- tau *= theta
- sigma /= theta
- step_pd2!(st; sigma, tau, theta)
- elseif ctx.params.algorithm == :newton
- step!(st)
- end
- #pvd[k] = save_denoise(st, k)
+ domain_factor = 1 / sqrt(area(mesh))
+ norm_step_ = norm_step(st) * domain_factor
+ #estimate_pd!(st)
- domain_factor = 1 / sqrt(area(mesh))
- norm_step_ = norm_step(st) * domain_factor
- #estimate_pd!(st)
+ log!(ctx.params.df; k,
+ norm_step = norm_step_,
+ #est = estimate_error(st),
+ primal_energy = primal_energy(st))
- log!(ctx.params.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
- #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
+ pvd[1] = save_step(st, 0)
end
- #pvd[1] = save_denoise(st, 1)
- #vtk_save(pvd)
return st
end
@@ -896,6 +899,66 @@ function experiment_convergence_rate(ctx)
energy_min, algparams...)
end
+function experiment_convergence_difference(ctx)
+ g_arr = from_img([0. 0.; 1. 0.])
+ mesh = init_grid(g_arr;)
+
+ algparams = (
+ init_zero = true,
+ alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
+ gamma1 = 1e-10, gamma2 = 1e-10,
+ tol = 1e-10,
+ max_iters = 10000,
+ )
+ algctx = Util.Context(ctx; g_arr, mesh, algparams...)
+
+ df1 = DataFrame()
+ df2 = DataFrame()
+ st1 = denoise_pd(Util.Context(algctx; algorithm = :pd2, df = df1));
+ st2 = denoise_pd(Util.Context(algctx; algorithm = :newton, df = df2));
+ p1 = plot(df1.k, df1.norm_step; axis = :log)
+ p2 = plot(df2.k, df2.norm_step; axis = :log)
+ p3 = plot(df1.k, df1.primal_energy; xaxis = :log)
+ p4 = plot(df2.k, df2.primal_energy; xaxis = :log)
+ display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
+
+ display(st1.u.data)
+ display(st2.u.data)
+
+ return
+end
+
+function experiment_convergence_oscillations(ctx)
+ g_arr = from_img([0. 0.; 1. 0.])
+ mesh = init_grid(g_arr;)
+
+ algparams = (
+ #alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
+ #alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
+ alpha1 = 1., alpha2 = 0., lambda = 0., beta = 1., # oscillating
+ # DP0 for p1 has less oscillations
+ gamma1 = 1e-5, gamma2 = 1e-5,
+ tol = 1e-10,
+ max_iters = 10000,
+ )
+ algctx = Util.Context(ctx; g_arr, mesh, algparams...)
+
+ df1 = DataFrame()
+ df2 = DataFrame()
+ st1 = denoise_pd(Util.Context(algctx; algorithm = :pd2, df = df1));
+ st2 = denoise_pd(Util.Context(algctx; algorithm = :newton, df = df2));
+ p1 = plot(df1.k, df1.norm_step; axis = :log)
+ p2 = plot(df2.k, df2.norm_step; axis = :log)
+ p3 = plot(df1.k, df1.primal_energy; xaxis = :log)
+ p4 = plot(df2.k, df2.primal_energy; xaxis = :log)
+ display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
+
+ display(st1.u.data)
+ display(st2.u.data)
+
+ return
+end
+
function denoise_approximation(ctx)
domain_factor = 1 / sqrt(area(ctx.params.mesh))
@@ -1016,6 +1079,7 @@ function denoise_approximation(ctx)
return st
end
+
function experiment_approximation(ctx)
#g_arr = from_img([1. 0.; 0. 0.])
g_arr = from_img([0. 0.; 1. 0.])
--
GitLab