From f3c21d636e051482b17d344023aad34940ca725f Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Fri, 9 Jul 2021 19:55:35 +0200
Subject: [PATCH] remove old run function
---
src/run.jl | 170 -----------------------------------------------------
1 file changed, 170 deletions(-)
diff --git a/src/run.jl b/src/run.jl
index 2f46eb8..3a05e85 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -264,173 +264,3 @@ function optflow(imgf0, imgf1; name, params...)
println()
end
end
-
-
-
-
-
-
-function myrun()
- name = "test"
- mesh = init_grid(50)
- d = 2
- m = 2
-
- Vg = FeSpace(mesh, P1(), (1,))
- Vu = FeSpace(mesh, P1(), (m,))
- Vp1 = FeSpace(mesh, DP0(), (1,))
- Vp2 = FeSpace(mesh, DP1(), (m, d))
-
- # inpainting
- mask = FeFunction(Vg, name="mask")
-
- g = FeFunction(Vg, name="g")
- u = FeFunction(Vu, name="u")
- p1 = FeFunction(Vp1, name="p1")
- p2 = FeFunction(Vp2, name="p2")
- du = FeFunction(Vu)
- dp1 = FeFunction(Vp1)
- dp2 = FeFunction(Vp2)
- nablau = nabla(u)
- nabladu = nabla(du)
-
- g.data .= 0
- u.data .= 0
- p1.data .= 0
- p2.data .= 0
- du.data .= 0
- dp1.data .= 0
- dp2.data .= 0
-
- alpha1 = 0.
- alpha2 = 10.
- beta = 1e-3
- lambda = 0.01
- gamma1 = 1e-3
- gamma2 = 1e-3
-
- interpolate!(g, x -> norm(x - SA[0.5, 0.5]) < 0.3)
-
- interpolate!(f0, x -> x[1])
- interpolate!(f1, x -> x[1] - 0.01)
- fw.data .= f1.data
- g_optflow(x; u, f0, fw, nablafw) =
- nablafw * u - (fw - f0)
- interpolate!(g, g_optflow; u, f0, fw, nablafw)
-
- #T(data, u) = u
- #tdata = u # FIXME
-
- #T(data, u) = iszero(data) ? zero(u) : u
- #tdata = mask
-
- T(data, u) = data * u
- tdata = nabla(fw)
-
- #S(u, nablau) = u
- S(u, nablau) = nablau
-
- function dp1_update(x; g, u, p1, du, tdata)
- m1 = max(gamma1, norm(T(tdata, u) - g))
- cond = norm(T(tdata, u) - g) > gamma1 ?
- dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
- zero(p1)
- return -p1 + alpha1 / m1 * (T(tdata, u) + T(tdata, du) - g) - cond
- end
-
-
- function dp2_update(x; u, nablau, p2, du, nabladu)
- m2 = max(gamma2, norm(nablau))
- cond = norm(nablau) > gamma2 ?
- dot(nablau, nabladu) / norm(nablau)^2 * p2 :
- zero(p2)
- return -p2 + lambda / m2 * (nablau + nabladu) - cond
- end
-
- function p1_project!()
- p1.space.element::DP0
- p1.data .= clamp.(p1.data, -alpha1, alpha1)
- end
-
- function p2_project!()
- p2.space.element::DP1
- p2d = reshape(p2.data, m * d * 3, :) # no copy
- for i in axes(p2d, 2)
- p2in = norm(p2d[:, i])
- if p2in > lambda
- p2d[:, i] .*= lambda ./ p2in
- end
- end
- end
-
- function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
- m1 = max(gamma1, norm(T(tdata, u) - g))
- cond1 = norm(T(tdata, u) - g) > gamma1 ?
- dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
- zero(p1)
- a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
- dot(cond1, T(tdata, phi))
-
- m2 = max(gamma2, norm(nablau))
- cond2 = norm(nablau) > gamma2 ?
- dot(nablau, nabladu) / norm(nablau)^2 * p2 :
- zero(p2)
- a2 = lambda / m2 * dot(nabladu, nablaphi) -
- dot(cond2, nablaphi)
-
- aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
- beta * dot(S(du, nabladu), S(phi, nablaphi))
-
- return a1 + a2 + aB
- end
-
- function du_l(x, phi, nablaphi; g, u, nablau, tdata)
- aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
- beta * dot(S(u, nablau), S(phi, nablaphi))
- m1 = max(gamma1, norm(T(tdata, u) - g))
- p1part = alpha1 / m1 * dot(T(tdata, u) - g, T(tdata, phi))
- m2 = max(gamma2, norm(nablau))
- p2part = lambda / m2 * dot(nablau, nablaphi)
- gpart = alpha2 * dot(g, T(tdata, phi))
-
- return -aB - p1part - p2part + gpart
- end
-
- function save_step!(pvd, i)
- vtk = vtk_mesh("$(name)_$(lpad(i, 5, '0')).vtu", mesh)
- vtk_append!(vtk, g, u, p1, p2)
- vtk_save(vtk)
- pvd[i] = vtk
- end
-
- pvd = paraview_collection("$(name).pvd")
- save_step!(pvd, 0)
- for i = 1:20
- print("[$(lpad(i, 3, '0'))] ")
- # solve du
- print("assemble ... ")
- A = assemble(Vu, du_a; g, u, nablau, p1, p2, tdata)
- b = assemble_rhs(Vu, du_l; g, u, nablau, tdata)
- print("solve ... ")
- du.data .= A \ b
-
- # solve dp1, dp2
- interpolate!(dp1, dp1_update; g, u, p1, du, tdata)
- interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu)
-
- # newton update
- u.data .+= du.data
- p1.data .+= dp1.data
- p2.data .+= dp2.data
-
- # reproject p1, p2
- p1_project!()
- p2_project!()
-
- print("save ... ")
- save_step!(pvd, i)
- println()
- end
- vtk_save(pvd)
- return
-end
--
GitLab