finalize primal-dual comparison experiment

scripts/run_experiments.jl

@@ -410,6 +410,9 @@ function estimate!(ctx::L1L2TVContext)
         nablau = nabla(ctx.u), w, nablaw = nabla(w), ctx.tdata)
 end
 
+estimate_error(st::L1L2TVContext) =
+    sqrt(sum(st.est.data) / area(st.mesh))
+
 # minimal Dörfler marking
 function mark(ctx::L1L2TVContext; theta=0.5)
     n = ncells(ctx.mesh)
@@ -529,113 +532,132 @@ function denoise(img; name, params...)
 end
 
 
-function denoise_pd(ctx, img; df=nothing, name, algorithm, params_...)
+function denoise_pd(st, img; df=nothing, name, algorithm, params_...)
     params = NamedTuple(params_)
     m = 1
     img = from_img(img) # coord flip
-    mesh = init_grid(img; type=:vertex)
+    mesh = init_grid(img)
     #mesh = init_grid(img, 5, 5)
 
     T(tdata, u) = u
     S(u, nablau) = u
 
-    ctx = L1L2TVContext(name, mesh, m;
+    st = L1L2TVContext(name, mesh, m;
         T, tdata = nothing, S,
         params.alpha1, params.alpha2, params.lambda, params.beta,
         params.gamma1, params.gamma2)
 
     # semi-implicit primal dual parameters
-    gamma = ctx.alpha2 + ctx.beta # T = I, S = I
-    gamma /= 100 # kind of arbitrary?
+    mu = st.alpha2 + st.beta # T = I, S = I
+    #mu /= 100 # kind of arbitrary?
 
     tau = 1e-1
     L = 100
     sigma = inv(tau * L^2)
     theta = 1.
 
-    #project_img!(ctx.g, img)
-    interpolate!(ctx.g, x -> interpolate_bilinear(img, x))
-    ctx.u.data .= ctx.g.data
+    #project_img!(st.g, img)
+    interpolate!(st.g, x -> interpolate_bilinear(img, x))
+    st.u.data .= st.g.data
 
-    save_denoise(ctx, i) =
-	output(ctx, "output/$(ctx.name)_$(lpad(i, 5, '0')).vtu",
-	    ctx.g, ctx.u, ctx.p1, ctx.p2)
+    save_denoise(st, i) =
+	output(st, "output/$(st.name)_$(lpad(i, 5, '0')).vtu",
+	    st.g, st.u, st.p1, st.p2)
 
     log!(x::Nothing; kwargs...) = x
-    function log!(df::DataFrame; k, norm_step, norm_residual)
-	push!(df, (;
-	    k,
-	    primal_energy = primal_energy(ctx),
-	    norm_step,
-	    norm_residual))
-	println(NamedTuple(last(df)))
+    function log!(df::DataFrame; kwargs...)
+        row = NamedTuple(kwargs)
+        push!(df, row)
+        #println(df)
+	println(row)
     end
 
-    #pvd = paraview_collection("output/$(ctx.name).pvd")
-    #pvd[0] = save_denoise(ctx, 0)
+    #pvd = paraview_collection("output/$(st.name).pvd")
+    #pvd[0] = save_denoise(st, 0)
 
     k = 0
-    println("primal energy: $(primal_energy(ctx))")
+    println("primal energy: $(primal_energy(st))")
 
     while true
 	k += 1
 	if algorithm == :pd1
 	    # no step size control
-	    step_pd2!(ctx; sigma, tau, theta)
+	    step_pd2!(st; sigma, tau, theta)
 	elseif algorithm == :pd2
-	    theta = 1 / sqrt(1 + 2 * gamma * tau)
+	    theta = 1 / sqrt(1 + 2 * mu * tau)
 	    tau *= theta
 	    sigma /= theta
-	    step_pd2!(ctx; sigma, tau, theta)
+	    step_pd2!(st; sigma, tau, theta)
 	elseif algorithm == :newton
-	    step!(ctx)
+	    step!(st)
 	end
-	#pvd[k] = save_denoise(ctx, k)
+	#pvd[k] = save_denoise(st, k)
 
 	domain_factor = 1 / sqrt(area(mesh))
-	norm_step_ = norm_step(ctx) * domain_factor
-	norm_residual_ = norm_residual(ctx) * domain_factor
+	norm_step_ = norm_step(st) * domain_factor
+        #estimate!(st)
 
-	log!(df; k, norm_step = norm_step_, norm_residual = norm_residual_)
+        log!(df; k,
+            norm_step = norm_step_,
+            #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
         haskey(params, :max_iters) && k >= params.max_iters && break
     end
-    #pvd[1] = save_denoise(ctx, 1)
+    #pvd[1] = save_denoise(st, 1)
     #vtk_save(pvd)
-    return ctx
+    return st
 end
 
 
 function experiment_pd_comparison(ctx)
     img = loadimg(joinpath(ctx.indir, "input.png"))
-    img = from_img(img) # coord flip
+    #img = [0. 0.; 1. 0.]
 
     algparams = (
         alpha1=0., alpha2=30., lambda=1., beta=0.,
-        gamma1=1e-3, gamma2=1e-3,
-        tol = 1e-6, max_iters = 50,
+        gamma1=1e-2, gamma2=1e-3,
+        tol = 1e-10,
+        max_iters = 10000,
         )
 
     df1 = DataFrame()
     df2 = DataFrame()
     df3 = DataFrame()
 
-    denoise_pd(ctx, img; name="test", algorithm=:pd1, df = df1, algparams...);
-    denoise_pd(ctx, img; name="test", algorithm=:pd2, df = df2, algparams...);
-    denoise_pd(ctx, img; name="test", algorithm=:newton, df = df3, algparams...);
+    st1 = denoise_pd(ctx, img; name="test",
+        algorithm=:pd1, df = df1, algparams...);
+    st2 = denoise_pd(ctx, img; name="test",
+        algorithm=:pd2, df = df2, algparams...);
+    st3 = denoise_pd(ctx, img; name="test",
+        algorithm=:newton, df = df3, algparams...);
 
     energy_min = min(minimum(df1.primal_energy), minimum(df2.primal_energy),
 		     minimum(df3.primal_energy))
-
-    #df1.primal_energy .-= energy_min
-    #df2.primal_energy .-= energy_min
-    #df3.primal_energy .-= energy_min
-
-    CSV.write(joinpath(ctx.outdir, "semiimplicit.csv"), logfilter(df1))
-    CSV.write(joinpath(ctx.outdir, "semiimplicit-accelerated.csv"), logfilter(df2))
-    CSV.write(joinpath(ctx.outdir, "newton.csv"), logfilter(df3))
+    df1[!, :energy_min] .= energy_min
+    df2[!, :energy_min] .= energy_min
+    df3[!, :energy_min] .= energy_min
+
+    CSV.write(joinpath(ctx.outdir, "semi-implicit.csv"),
+        logfilter(df1))
+    CSV.write(joinpath(ctx.outdir, "semi-implicit-accelerated.csv"),
+        logfilter(df2))
+    CSV.write(joinpath(ctx.outdir, "newton.csv"),
+        logfilter(df3))
+
+    saveimg(joinpath(ctx.outdir, "input.png"),
+        img)
+    saveimg(joinpath(ctx.outdir, "semi-implicit.png"),
+        to_img(sample(st1.u)))
+    saveimg(joinpath(ctx.outdir, "semi-implicit-accelerated.png"),
+        to_img(sample(st2.u)))
+    saveimg(joinpath(ctx.outdir, "newton.png"),
+        to_img(sample(st3.u)))
+
+    savedata(joinpath(ctx.outdir, "data.tex");
+        energy_min, algparams...)
 end
 
 function denoise_approximation(ctx)
@@ -662,7 +684,8 @@ function denoise_approximation(ctx)
     function log!(df::DataFrame; kwargs...)
         row = NamedTuple(kwargs)
         push!(df, row)
-	println(row)
+        println(df)
+	#println(row)
     end
 
     pvd_path = joinpath(ctx.outdir, "$(ctx.params.name).pvd")
@@ -716,9 +739,9 @@ function experiment_approximation(ctx)
 
     denoise_approximation(Util.Context(ctx; name = "test", df,
         img, mesh,
-        #alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
-        alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
-        gamma1 = 1e-5, gamma2 = 1e-3,
+        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,
     ))
 end

scripts/util.jl

@@ -44,13 +44,12 @@ end
 
 # LaTeX Data Output
 
-# TODO: don't depend on ctx.path and use filename only
-savedata(ctx::Context, filename; data...) =
-    open(joinpath(ctx.outdir, filename); write=true) do io
-        _savedata(io, ctx.path; data...)
+savedata(filename; data...) =
+    open(filename; write=true) do io
+        _savedata(io; data...)
     end
 
-function _savedata(io, path; data...)
+function _savedata(io; data...)
     isvalidname(x) = contains(string(x), r"^[^/\\]+$")
 
     # define tex command to access data
@@ -67,7 +66,7 @@ function _savedata(io, path; data...)
 
         # actual data encoded as tex-commands
         write(io, "\\expandafter\\def")
-        write(io, "\\csname $(texcmd)/$(path)/$(key)\\endcsname{")
+        write(io, "\\csname $(texcmd)/\\DataPrefix/$(key)\\endcsname{")
         show(io, MIME("text/latex"), value)
         write(io, "}\n")
     end