correct space and parameter scaling in step_pd1

src/run.jl

@@ -43,7 +43,7 @@ function L1L2TVContext(name, mesh, m; T, tdata, S,
     Vest = FeSpace(mesh, DP0(), (1,))
     Vg = FeSpace(mesh, P1(), (1,))
     Vu = FeSpace(mesh, P1(), (m,))
-    Vp1 = FeSpace(mesh, DP0(), (1,))
+    Vp1 = FeSpace(mesh, P1(), (1,))
     Vp2 = FeSpace(mesh, DP1(), (m, d))
 
     est = FeFunction(Vest, name="est")
@@ -151,7 +151,8 @@ function step!(ctx::L1L2TVContext)
 
     # reproject p1
     function p1_project!(p1, alpha1)
-	p1.space.element::DP0
+	p1.space.element == DP0() || p1.space.element == P1() ||
+	    throw(ArgumentError("element unsupported"))
 	p1.data .= clamp.(p1.data, -alpha1, alpha1)
     end
     p1_project!(ctx.p1, ctx.alpha1)
@@ -175,11 +176,12 @@ end
 "
 function step_pd!(ctx::L1L2TVContext; sigma, tau, theta = 1.)
     # note: ignores gamma1, gamma2, beta and uses T = I, lambda = 1, m = 1!
+    # changed alpha2 -> alpha2 / 2
     # chambolle-pock require: sigma * tau * L^2 <= 1, L = |grad|
     if ctx.m != 1 || ctx.lambda != 1. || ctx.beta != 0.
 	error("unsupported parameters")
     end
-    beta = tau * ctx.alpha1 / (1 + 2 * tau * ctx.alpha2)
+    beta = tau * ctx.alpha1 / (1 + tau * ctx.alpha2)
 
     # u is P1
     # p2 is essentially DP0 (technically may be DP1)
@@ -215,11 +217,11 @@ function step_pd!(ctx::L1L2TVContext; sigma, tau, theta = 1.)
 	dot(z, phi)
 
     u_l(x, phi, nablaphi; u, g, p2) =
-	(dot(u + 2 * tau * ctx.alpha2 * g, phi) + tau * dot(p2, -nablaphi)) /
-	    (1 + 2 * tau * ctx.alpha2)
+	(dot(u + tau * ctx.alpha2 * g, phi) + tau * dot(p2, -nablaphi)) /
+	    (1 + tau * ctx.alpha2)
 
-    # z = 1 / (1 + 2 * tau * alpha2) *
-    #   (u + 2 * tau * alpha2 * g + tau * div(p))
+    # z = 1 / (1 + tau * alpha2) *
+    #   (u + tau * alpha2 * g + tau * div(p))
     z = FeFunction(ctx.u.space)
     A, b = assemble(z.space, u_a, u_l; ctx.g, ctx.u, ctx.p2)
     z.data .= A \ b
@@ -283,7 +285,8 @@ function step_pd2!(ctx::L1L2TVContext; sigma, tau, theta = 1.)
 
     # reproject p1
     function p1_project!(p1, alpha1)
-	p1.space.element::DP0
+	p1.space.element == DP0() || p1.space.element == P1() ||
+	    throw(ArgumentError("element unsupported"))
 	p1.data .= clamp.(p1.data, -alpha1, alpha1)
     end
     p1_project!(p1_next, ctx.alpha1)
@@ -538,7 +541,7 @@ function denoise_pd(img; name, params...)
     gamma = ctx.alpha2 + ctx.beta # T = I, S = I
     sigma = 1e-2
     tau = 1e-2
-    theta = 0.
+    theta = 1.
 
     project_img!(ctx.g, img)
     #interpolate!(ctx.g, x -> interpolate_bilinear(img, x))
@@ -556,16 +559,15 @@ function denoise_pd(img; name, params...)
     k = 0
     println("primal energy: $(primal_energy(ctx))")
 
-    algorithm = :newton
+    algorithm = :pd2
     while true
 	k += 1
 	#println("")
 	if algorithm == :pd2 #|| primal_energy(ctx) < 2409
 	    println("step_pd: theta = $theta, tau = $tau, sigma = $sigma")
-	    pd = true
-	    theta = 1 / sqrt(1 + 2 * gamma * tau)
-	    tau *= theta
-	    sigma /= theta
+	    #theta = 1 / sqrt(1 + 2 * gamma * tau)
+	    #tau *= theta
+	    #sigma /= theta
 	    step_pd2!(ctx; sigma, tau, theta)
 	elseif algorithm == :pd1
 	    # no step size control