add convergence experiments

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))
+    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,9 +812,7 @@ 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))")
 
@@ -827,7 +829,7 @@ function denoise_pd(ctx)
             elseif ctx.params.algorithm == :newton
                 step!(st)
             end
-	#pvd[k] = save_denoise(st, k)
+            #pvd[k] = save_step(st, k)
 
             domain_factor = 1 / sqrt(area(mesh))
             norm_step_ = norm_step(st) * domain_factor
@@ -842,8 +844,9 @@ function denoise_pd(ctx)
             haskey(params, :tol) && norm_step_ < params.tol && break
             haskey(params, :max_iters) && k >= params.max_iters && break
         end
-    #pvd[1] = save_denoise(st, 1)
-    #vtk_save(pvd)
+
+        pvd[1] = save_step(st, 0)
+    end
     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.])