add dual energy

scripts/run_experiments.jl

@@ -405,7 +405,7 @@ function estimate!(ctx::L1L2TVContext)
 end
 
 estimate_error(st::L1L2TVContext) =
-    sqrt(sum(st.est.data) / area(st.mesh))
+    sum(st.est.data) / area(st.mesh)
 
 # minimal Dörfler marking
 function mark(ctx::L1L2TVContext; theta=0.5)
@@ -445,6 +445,26 @@ function primal_energy(ctx::L1L2TVContext)
         nablau = nabla(ctx.u), ctx.tdata)
 end
 
+function dual_energy(st::L1L2TVContext)
+    # primal reconstruction
+    w = FeFunction(st.u.space)
+    solve_primal!(w, st)
+
+    # 1 / 2 * \|w\|_B^2 - alpha2 / 2 * \|g\| + <g, p_1> +
+    #   gamma1 / 2 / alpha1 * \|p_1\|^2 + gamma2 / 2 / lambda * \|p_2\|^2
+    function integrand(x; g, tdata, w, nablaw, p1, p2)
+        return 1 / 2 * (
+                st.alpha2 * dot(st.T(tdata, w), st.T(tdata, w)) +
+                st.beta * dot(st.S(w, nablaw), st.S(w, nablaw))) +
+            - st.alpha2 / 2 * dot(g, g) +
+            dot(g, p1) +
+            (iszero(st.alpha1) ? 0 : st.gamma1 / 2 / st.alpha1 * dot(p1, p1)) +
+            (iszero(st.lambda) ? 0 : st.gamma2 / 2 / st.lambda * dot(p2, p2))
+    end
+    return integrate(st.mesh, integrand; st.g, st.tdata,
+        w, nablaw = nabla(w), st.p1, st.p2)
+end
+
 norm_l2(f) = sqrt(integrate(f.space.mesh, (x; f) -> dot(f, f); f))
 
 norm_step(ctx::L1L2TVContext) =
@@ -693,6 +713,20 @@ function denoise_approximation(ctx)
             step!(st)
 
             norm_step_ = norm_step(st) * domain_factor
+
+            k += 1
+            norm_residual_ = norm_residual(st) * domain_factor
+            estimate!(st)
+
+            pvd[k] = save_vtk(st, k)
+            log!(ctx.params.df; k,
+                ndofs = ndofs(st.u.space),
+                hmax = mesh_size(st.mesh),
+                norm_step = norm_step_, norm_residual = norm_residual_,
+                primal_energy = primal_energy(st),
+                dual_energy = dual_energy(st),
+                est = estimate_error(st),
+            )
             if haskey(ctx.params, :tol) && norm_step_ < ctx.params.tol
                 k += 1
                 norm_residual_ = norm_residual(st) * domain_factor
@@ -704,19 +738,32 @@ function denoise_approximation(ctx)
                     hmax = mesh_size(st.mesh),
                     norm_step = norm_step_, norm_residual = norm_residual_,
                     primal_energy = primal_energy(st),
-                    est = sqrt(sum(st.est.data)) * domain_factor,
+                    dual_energy = dual_energy(st),
+                    est = estimate_error(st),
                 )
 
                 marked_cells = ctx.params.adaptive ?
                     mark(st; theta = 0.5) : # estimator + dörfler
                     Set(axes(st.mesh.cells, 2))
-
                 mesh, fs = refine(st.mesh, marked_cells;
                     st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
                 st = L1L2TVContext("run", mesh, st.d, st.m, T, nothing, S,
                     st.alpha1, st.alpha2, st.beta, st.lambda,
                     st.gamma1, st.gamma2,
                     fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
+
+                norm_residual_ = norm_residual(st) * domain_factor
+                estimate!(st)
+
+                pvd[k] = save_vtk(st, k)
+                log!(ctx.params.df; k,
+                    ndofs = ndofs(st.u.space),
+                    hmax = mesh_size(st.mesh),
+                    norm_step = norm_step_, norm_residual = norm_residual_,
+                    primal_energy = primal_energy(st),
+                    dual_energy = dual_energy(st),
+                    est = estimate_error(st),
+                )
             end
         end
     end
@@ -735,8 +782,8 @@ function experiment_approximation(ctx)
         img, mesh,
         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,
+        gamma1 = 1e-4, gamma2 = 1e-4,
+        tol = 1e-10, adaptive = false,
     ))
 end