tweak synthetic experiments

scripts/run.jl

@@ -5,5 +5,26 @@ const outpath = "/home/stev47/stuff/phd/data"
 
 ctx = Util.Context(datapath, outpath)
 
+# convergence
+
+# algorithm comparison
 ctx(experiment_convergence_rate, "fem/convergence/rate")
+# L1 convergence, difference in algorithm
+
+# L1 non-convergence / oscillations
+
+
+# adaptivity
+
+# residual estimator, discrete g, good
+# pd est, L2-TV, discrete g, good?
+# pd est, zero
+# pd est, L1-TV, bad
+
+
+# applications
+
+# denoising
+# adaptive
+# optical flow
 ctx(experiment_optflow_middlebury_all, "fem/optflow/middlebury")

scripts/run_experiments.jl

@@ -857,12 +857,12 @@ function denoise_pd(ctx)
 
     # semi-implicit primal dual parameters
     mu = st.alpha2 + st.beta # T = I, S = I
-    #mu /= 100 # kind of arbitrary?
+    mu /= 100 # kind of arbitrary?
 
     tau = 1e-1
-    L = 100
+    L = norm_gradient(mesh)
     sigma = inv(tau * L^2)
-    theta = 1.
+    theta = NaN # initialized later
 
     #project_l2_lagrange!(st.g, ctx.params.g_arr)
     interpolate!(st.g, x -> evaluate_bilinear(ctx.params.g_arr, x))
@@ -909,6 +909,7 @@ function denoise_pd(ctx)
 
             log!(ctx.params.df; k,
                 norm_step = norm_step_,
+                norm_residual = norm_residual(st),
                 #est = estimate_error(st),
                 primal_energy = primal_energy(st))
 
@@ -978,7 +979,7 @@ function experiment_convergence_difference(ctx)
     algparams = (
         init_zero = true,
         alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
-        gamma1 = 1e-10, gamma2 = 1e-10,
+        gamma1 = 1e-7, gamma2 = 1e-7,
         tol = 1e-10,
         max_iters = 10000,
     )
@@ -990,12 +991,16 @@ function experiment_convergence_difference(ctx)
     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)
+    p3 = plot(df1.k, df1.norm_residual; axis = :log)
+    p4 = plot(df2.k, df2.norm_residual; axis = :log)
     display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
 
     display(st1.u.data)
+    display(norm_residual(st1))
+    display(primal_energy(st1))
     display(st2.u.data)
+    display(norm_residual(st2))
+    display(primal_energy(st2))
 
     return
 end
@@ -1005,11 +1010,12 @@ function experiment_convergence_oscillations(ctx)
     mesh = init_grid(g_arr;)
 
     algparams = (
+        #init_zero = true,
         #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
+        alpha1 = 0.5, alpha2 = 0., lambda = 1., beta = 1., # oscillating
         # DP0 for p1 has less oscillations
-        gamma1 = 1e-5, gamma2 = 1e-5,
+        gamma1 = 1e-7, gamma2 = 1e-7,
         tol = 1e-10,
         max_iters = 10000,
     )
@@ -1026,7 +1032,9 @@ function experiment_convergence_oscillations(ctx)
     display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
 
     display(st1.u.data)
+    display(primal_energy(st1))
     display(st2.u.data)
+    display(primal_energy(st2))
 
     return
 end