From d9d7d55bdfc68553edcf95da3b993afc35fa4570 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Mon, 28 Feb 2022 19:44:41 +0100
Subject: [PATCH] update denoise_approximation experiment

---
 scripts/run_experiments.jl | 48 +++++++++++++++++++++++---------------
 1 file changed, 29 insertions(+), 19 deletions(-)

diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 8a316a2..c573932 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -693,7 +693,7 @@ function estimate_res!(st::L1L2TVState)
             fj = mesh.cells[j, cell]
 
             gdofs = space.dofmap[:, 1, cell]
-            st.est.data[gdofs] .+= (isSnabla ? hfacet : inv(hfacet)) *
+            st.est.data[gdofs] .+=  (isSnabla ? hfacet : inv(hfacet)) *
                 norm2(facetv[(fi, fj)]) * intel
         end
     end
@@ -921,8 +921,8 @@ function denoise_pd(ctx)
                 #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
+            #norm_residual_ < params.eps_newton && norm_step_ < params.eps_newton && break
+            haskey(params, :eps_newton) && norm_step_ < params.eps_newton && break
             haskey(params, :max_iters) && k >= params.max_iters && break
         end
 
@@ -941,7 +941,7 @@ function experiment_convergence_rate(ctx)
     algparams = (
         alpha1=0., alpha2=30., lambda=1., beta=0.,
         gamma1=1e-2, gamma2=1e-3,
-        tol = 1e-10,
+        eps_newton = 1e-10,
         max_iters = 10000,
     )
     algctx = Util.Context(ctx; g_arr, mesh, algparams...)
@@ -988,7 +988,7 @@ function experiment_convergence_difference(ctx)
         init_zero = true,
         alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
         gamma1 = 1e-7, gamma2 = 1e-7,
-        tol = 1e-10,
+        eps_newton = 1e-10,
         max_iters = 10000,
     )
     algctx = Util.Context(ctx; g_arr, mesh, algparams...)
@@ -1024,7 +1024,7 @@ function experiment_convergence_oscillations(ctx)
         alpha1 = 0.5, alpha2 = 0., lambda = 1., beta = 1., # oscillating
         # DP0 for p1 has less oscillations
         gamma1 = 1e-7, gamma2 = 1e-7,
-        tol = 1e-10,
+        eps_newton = 1e-10,
         max_iters = 10000,
     )
     algctx = Util.Context(ctx; g_arr, mesh, algparams...)
@@ -1053,8 +1053,8 @@ function denoise_approximation(ctx)
     m = 1
     T(tdata, u) = u
     T(::typeof(adjoint), tdata, v) = v
-    S(u, nablau) = u
-    #S(u, nablau) = nablau
+    #S(u, nablau) = u
+    S(u, nablau) = nablau'
 
     st = L1L2TVState{m}(ctx.params.mesh;
         T, tdata = nothing, S,
@@ -1066,8 +1066,15 @@ function denoise_approximation(ctx)
     w = FeFunction(st.u.space, name = "w")
     w.data .= 0
 
-    #project_l2_lagrange!(st.g, img)
-    interpolate!(st.g, x -> evaluate_bilinear(ctx.params.g_arr, x))
+    #project_image! = project_l2_lagrange!
+    project_image!(f, img) = interpolate!(f, x -> evaluate_bilinear(img, x))
+
+    function interpolate_data!()
+        project_image!(st.g, ctx.params.g_arr)
+        interpolate!(st.g, x -> float(norm(x .- 1., Inf) <= 3/4))
+    end
+
+    interpolate_data!()
     st.u.data .= st.g.data
     #st.u.data .= rand(size(st.u.data))
 
@@ -1107,13 +1114,14 @@ function denoise_approximation(ctx)
 
         k = 0
         k_refine = 0
-        while true && ndofs(st.u.space) < 1000 && k < 300
+        while true && ndofs(st.u.space) < 2000 && k < 300
+            k += 1
+
             step!(st)
+
             solve_primal!(w, st)
 
             norm_step_ = norm_step(st) * domain_factor
-
-            k += 1
             norm_residual_ = norm_residual(st) * domain_factor
             estimate_res!(st)
 
@@ -1126,7 +1134,7 @@ function denoise_approximation(ctx)
                 dual_energy = dual_energy(st),
                 est = estimate_error(st),
             )
-            if haskey(ctx.params, :tol) && norm_step_ < ctx.params.tol
+            if haskey(ctx.params, :eps_newton) && norm_step_ < ctx.params.eps_newton
                 k += 1
                 k_refine += 1
                 norm_residual_ = norm_residual(st) * domain_factor
@@ -1151,6 +1159,8 @@ function denoise_approximation(ctx)
                     fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
                 w = fs.w
 
+                interpolate_data!()
+
                 norm_residual_ = norm_residual(st) * domain_factor
                 estimate_res!(st)
 
@@ -1173,18 +1183,18 @@ end
 
 function experiment_approximation(ctx)
     phi = (sqrt(5) + 1) / 2
-    #g_arr = from_img([1. 0.; 0. 0.])
     g_arr = from_img([0. 0.; 1. 0.])
+    #g_arr = rand(30, 30)
     mesh = init_grid(g_arr;)
 
     df = DataFrame()
     denoise_approximation(Util.Context(ctx; name = "test", df,
         g_arr, mesh,
-        alpha1 = 0., alpha2 = 20., lambda = 1., beta = 0.,
+        alpha1 = 0., alpha2 = 100., lambda = 1., beta = 1e-5,
         #alpha1 = 1. - inv(phi), alpha2 = 0., lambda = 0., beta = 1.,
-        #alpha1 = 1., alpha2 = 0., lambda = 1., beta = 1e-5,
-        gamma1 = 1e-5, gamma2 = 1e-5,
-        tol = 1e-10, adaptive = true,
+        #alpha1 = 2., alpha2 = 0., lambda = 1., beta = 1e-5,
+        gamma1 = 1e-4, gamma2 = 1e-4,
+        eps_newton = 1e-5, adaptive = true,
     ))
 end
 
-- 
GitLab