Select Git revision
TP-TP-layered_soil_pure_dd.py
problems.jl 1.96 KiB
using LinearAlgebra: UniformScaling, Diagonal, opnorm
using StaticKernels: ExtensionNothing
#"min_u 1/2 * |Au-f|_2^2 + λ*TV(u) + β/2 |u|_2 + γ |u|_1"
"min_p 1/2 * |p₂ - div(p₁) - g|_B^2 + χ_{|p₁|≤λ} + χ_{|p₂|≤γ}"
struct DualTVL1ROFOpProblem{Tg,TB,Tλ,Tγ} <: Problem
"dual data"
g::Tg
"B operator"
B::TB
"TV regularization parameter"
λ::Tλ
"L1 regularization parameter"
γ::Tγ
end
DualTVL1ROFOpProblem(g, B, λ) = DualTVL1ROFOpProblem(g, B, λ, nothing)
DualTVL1ROFOpProblem(g, B, λ::Real) = DualTVL1ROFOpProblem(g, B, fill!(similar(g, typeof(λ)), λ))
function energy(p, prob::DualTVL1ROFOpProblem)
d = ndims(p)
@inline kfΛ(w) = @inbounds divergence(w) + prob.g[w.position]
kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
# v = div(p) + g
v = map(kΛ, extend(p, ExtensionNothing()))
# |v|_B^2
u = prob.B * vec(v)
return sum(dot.(u, vec(v)))
end
function recover_u(p, prob::DualTVL1ROFOpProblem)
d = ndims(prob.g)
u = similar(prob.g)
v = similar(prob.g)
@inline kfΛ(w) = @inbounds divergence(w)
kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
# v = div(p) + A'*f
map!(kΛ, v, StaticKernels.extend(p, StaticKernels.ExtensionNothing())) # extension: nothing
v .+= prob.g
# u = B * v
mul!(vec(u), prob.B, vec(v))
return u
end
mynorm(v) = norm(vec(v))
@inline mynorm(v::SArray{<:Any,<:SArray{<:Any,Float64,1},1}) =
sqrt(sum(sum(v[i] .^ 2) for i in eachindex(v)))
function residual(p, prob::DualTVL1ROFOpProblem)
d = ndims(p)
grad = Kernel{ntuple(_->0:1, d)}(gradient)
u = recover_u(p, prob)
q = map(grad, StaticKernels.extend(u, StaticKernels.ExtensionReplicate()))
res = q .- mynorm.(q) .* p ./ prob.λ
return sum(dot.(res, res)) / length(p)
end
# operator norm of B
normB(p::DualTVL1ROFOpProblem{<:Any,<:UniformScaling}) = p.B.λ
normB(p::DualTVL1ROFOpProblem{<:Any,<:Diagonal}) = maximum(opnorm.(diag(p.B)))
normB(p::DualTVL1ROFOpProblem) = opnorm(p.B)