From 7d82e726328977e91bf8276d3dd693c44a3e92c6 Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan.hilb@mathematik.uni-stuttgart.de>
Date: Wed, 4 Sep 2019 17:54:44 +0200
Subject: [PATCH] resurrect old code

---
 DualTVDD/Manifest.toml             | 257 +++++++++++++++++++++++++++++
 DualTVDD/Project.toml              |  17 ++
 DualTVDD/src/DualTVDD.jl           | 123 ++++++++++++++
 FD/src/FD.jl => DualTVDD/src/fd.jl | 133 +--------------
 {FD => DualTVDD}/src/filtering.jl  |   8 +-
 DualTVDD/test/runtests.jl          |  26 +++
 FD/Manifest.toml                   |  65 --------
 FD/Project.toml                    |   8 -
 8 files changed, 429 insertions(+), 208 deletions(-)
 create mode 100644 DualTVDD/Manifest.toml
 create mode 100644 DualTVDD/Project.toml
 create mode 100644 DualTVDD/src/DualTVDD.jl
 rename FD/src/FD.jl => DualTVDD/src/fd.jl (80%)
 rename {FD => DualTVDD}/src/filtering.jl (94%)
 create mode 100644 DualTVDD/test/runtests.jl
 delete mode 100644 FD/Manifest.toml
 delete mode 100644 FD/Project.toml

diff --git a/DualTVDD/Manifest.toml b/DualTVDD/Manifest.toml
new file mode 100644
index 0000000..d95fb7f
--- /dev/null
+++ b/DualTVDD/Manifest.toml
@@ -0,0 +1,257 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[BinaryProvider]]
+deps = ["Libdl", "Logging", "SHA"]
+git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648"
+uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
+version = "0.5.6"
+
+[[ColorTypes]]
+deps = ["FixedPointNumbers", "Random"]
+git-tree-sha1 = "10050a24b09e8e41b951e9976b109871ce98d965"
+uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+version = "0.8.0"
+
+[[Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Printf", "Reexport"]
+git-tree-sha1 = "c9c1845d6bf22e34738bee65c357a69f416ed5d1"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.9.6"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "2.1.0"
+
+[[Contour]]
+deps = ["LinearAlgebra", "StaticArrays", "Test"]
+git-tree-sha1 = "b974e164358fea753ef853ce7bad97afec15bb80"
+uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
+version = "0.5.1"
+
+[[DataAPI]]
+git-tree-sha1 = "8903f0219d3472543fc4b2f5ebaf675a07f817c0"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.0.1"
+
+[[DataStructures]]
+deps = ["InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "0809951a1774dc724da22d26e4289bbaab77809a"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.17.0"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[FFMPEG]]
+deps = ["BinaryProvider", "Libdl"]
+git-tree-sha1 = "1dd2128ff10894081f30931b355dc892d1352de9"
+uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
+version = "0.2.2"
+
+[[FixedPointNumbers]]
+git-tree-sha1 = "d14a6fa5890ea3a7e5dcab6811114f132fec2b4b"
+uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+version = "0.6.1"
+
+[[GR]]
+deps = ["Base64", "DelimitedFiles", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test"]
+git-tree-sha1 = "b4c31b6377b6d51b6c69a3a9737d10c34d43974e"
+uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
+version = "0.41.0"
+
+[[GeometryTypes]]
+deps = ["ColorTypes", "FixedPointNumbers", "IterTools", "LinearAlgebra", "StaticArrays"]
+git-tree-sha1 = "4bf5706f3b9a2c5adbbc473c8c91582c1fa816a3"
+uuid = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+version = "0.7.6"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[IterTools]]
+git-tree-sha1 = "2ebe60d7343962966d1779a74a760f13217a6901"
+uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+version = "1.2.0"
+
+[[JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.0"
+
+[[LibGit2]]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Measures]]
+deps = ["Test"]
+git-tree-sha1 = "ddfd6d13e330beacdde2c80de27c1c671945e7d9"
+uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
+version = "0.3.0"
+
+[[Missings]]
+deps = ["SparseArrays", "Test"]
+git-tree-sha1 = "f0719736664b4358aa9ec173077d4285775f8007"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "0.4.1"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[NaNMath]]
+deps = ["Compat"]
+git-tree-sha1 = "ce3b85e484a5d4c71dd5316215069311135fa9f2"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "0.3.2"
+
+[[OffsetArrays]]
+git-tree-sha1 = "1af2f79c7eaac3e019a0de41ef63335ff26a0a57"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "0.11.1"
+
+[[OrderedCollections]]
+deps = ["Random", "Serialization", "Test"]
+git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.1.0"
+
+[[Parsers]]
+deps = ["Dates", "Test"]
+git-tree-sha1 = "ef0af6c8601db18c282d092ccbd2f01f3f0cd70b"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "0.3.7"
+
+[[Pkg]]
+deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[PlotThemes]]
+deps = ["PlotUtils", "Requires", "Test"]
+git-tree-sha1 = "f3afd2d58e1f6ac9be2cea46e4a9083ccc1d990b"
+uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
+version = "0.3.0"
+
+[[PlotUtils]]
+deps = ["Colors", "Dates", "Printf", "Random", "Reexport", "Test"]
+git-tree-sha1 = "8e87bbb778c26f575fbe47fd7a49c7b5ca37c0c6"
+uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
+version = "0.5.8"
+
+[[Plots]]
+deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryTypes", "JSON", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "Reexport", "Requires", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
+git-tree-sha1 = "f2aa8a7b5bc0ccec57a1237a97b6f59fc8d9ef57"
+uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+version = "0.26.2"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[RecipesBase]]
+git-tree-sha1 = "7bdce29bc9b2f5660a6e5e64d64d91ec941f6aa2"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "0.7.0"
+
+[[Reexport]]
+deps = ["Pkg"]
+git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "0.2.0"
+
+[[Requires]]
+deps = ["Test"]
+git-tree-sha1 = "f6fbf4ba64d295e146e49e021207993b6b48c7d1"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "0.5.2"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[Showoff]]
+deps = ["Dates"]
+git-tree-sha1 = "e032c9df551fb23c9f98ae1064de074111b7bc39"
+uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
+version = "0.3.1"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SortingAlgorithms]]
+deps = ["DataStructures", "Random", "Test"]
+git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "0.3.1"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "db23bbf50064c582b6f2b9b043c8e7e98ea8c0c6"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "0.11.0"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
+git-tree-sha1 = "c53e809e63fe5cf5de13632090bc3520649c9950"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.32.0"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
diff --git a/DualTVDD/Project.toml b/DualTVDD/Project.toml
new file mode 100644
index 0000000..dfc958c
--- /dev/null
+++ b/DualTVDD/Project.toml
@@ -0,0 +1,17 @@
+name = "DualTVDD"
+uuid = "e723064d-6d9e-46ff-9d98-de992daae271"
+authors = ["Stephan Hilb <stephan.hilb@mathematik.uni-stuttgart.de>"]
+version = "0.1.0"
+
+[deps]
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/DualTVDD/src/DualTVDD.jl b/DualTVDD/src/DualTVDD.jl
new file mode 100644
index 0000000..6a5fb55
--- /dev/null
+++ b/DualTVDD/src/DualTVDD.jl
@@ -0,0 +1,123 @@
+module DualTVDD
+
+include("filtering.jl")
+include("fd.jl")
+
+# import
+
+using LinearAlgebra: I
+using .FD
+
+# export
+
+export MetaGrid, GridFunction
+export dditer
+
+# algorithm
+
+
+"""
+    partition_function(metagrid, idx...)
+
+Create (global) partition of unity grid function θ.
+"""
+function partition_function(metagrid, idx...)
+  function h(r, rp, i)
+    i <= first(r) + metagrid.overlap && return 1.0
+    i >= last(r) - metagrid.overlap && return 1.0
+    (min(i - first(rp), last(rp) - i, metagrid.overlap) + 1) / (metagrid.overlap + 1)
+  end
+
+  vidx = FD.indices(metagrid, idx...)
+  ctheta = map(x -> [h(x[1], x[2], i) for i in x[2]], zip(axes(metagrid), vidx))
+
+  dst = zeros(size(metagrid))
+  for (I,J) in zip(CartesianIndices(vidx), CartesianIndices(length.(vidx)))
+    dst[I] = prod(getindex.(ctheta, Tuple(J)))
+  end
+  FD.GridFunction(metagrid, dst)
+end
+
+
+"""
+    p2u
+
+Compute primal solution from dual solution according to the optimality conditions.
+"""
+function p2u(grid, p, A, f; α, β)
+  FD.GridFunction(grid, (A'*A + β*I) \ (FD.divergence_z(p).data[:] .+ A'*f.data[:]))
+end
+
+
+"""
+    fpiter(grid, p, wg, A, α, β, τ)
+
+Fixed point iteration (Chambolle).
+"""
+# TODO: investigate matrix multiplication performance with GridFunctions
+function fpiter(grid, p, wg, A, α, β, τ)
+  x = FD.GridFunction(FD.parent(grid), (A'*A + β*I) \ (FD.divergence_z(FD.zero_continuation(p)) + wg).data[:])
+  q = FD.GridFunction(grid, view(FD.gradient(x).data, grid.indices..., :))
+  ## this misses out on the right boundary
+  #q = FD.gradient(FD.GridFunction(grid, view(x.data, grid.indices..., :)))
+
+  p_new = FD.GridFunction(grid, FD.pwmul(FD.GridFunction(grid, α ./ (α .+ τ .* FD.pw_norm(q))),
+                                         FD.GridFunction(grid, p .+ τ .* q)))
+
+  reserr = maximum(sum((p_new .- p).^2, dims=ndims(p)))
+  (p_new, reserr)
+end
+
+"""
+    fpalg(grid, w, A, p = 0; α, β, τ = 1/8, ϵ = 1e-4)
+
+Fixed point algorithm (Chambolle), calling `fpiter`.
+"""
+function fpalg(grid, w, A, p = FD.GridFunction(x -> 0, ndims(grid), grid);
+               α, β, τ = 1/8, ϵ = 1e-4)
+  reserr = 1
+  k = 0
+  while reserr > ϵ && k < 10e1
+    k += 1
+    (p, reserr) = fpiter(grid, p, w, A, α, β, τ)
+  end
+  return p
+end
+
+
+"""
+    dditer(grid, p, A, f; α, β, τ, ϵ, ρ)
+
+Domain decomposition algorithm, calling `fpalg` for each part.
+"""
+function dditer(grid, p, A, f; α, β, τ, ϵ, ρ)
+  p̂ = (1-ρ) * p
+
+  i = 0
+  for part in FD.parts(grid)
+    gridv = view(grid, Tuple(part)...)
+
+    fv = view(f, Tuple(part)...)
+
+    p̂v = view(p̂, Tuple(part)...)
+    pv = view(p, Tuple(part)...)
+
+    θ = partition_function(grid, Tuple(part)...)
+    θv = view(θ, Tuple(part)...)
+
+    w = FD.GridFunction(grid, A' * f.data[:]) + FD.divergence(FD.pwmul(FD.GridFunction(grid, 1 .- θ), p))
+    wv = view(w, Tuple(part)...)
+
+    p̂v .+= ρ .* fpalg(gridv, w, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
+
+    i += 1
+  end
+  p .= p̂
+
+  u = p2u(grid, p, A, f, α=α, β=β)
+  FD.my_plot(u)
+  #FD.my_plot(FD.pw_norm(p))
+end
+
+
+end # module
diff --git a/FD/src/FD.jl b/DualTVDD/src/fd.jl
similarity index 80%
rename from FD/src/FD.jl
rename to DualTVDD/src/fd.jl
index 072307c..5f76c3b 100644
--- a/FD/src/FD.jl
+++ b/DualTVDD/src/fd.jl
@@ -10,8 +10,6 @@ export Grid, MetaGrid, GridView, GridFunction
 export grid
 
 
-include("filtering.jl")
-
 "Finite Differences"
 module FDiff
 
@@ -30,7 +28,7 @@ end
 end
 
 
-using .Filtering
+using ..Filtering
 using .FDiff
 
 """
@@ -462,132 +460,3 @@ function solve(pde, shape::NTuple{d, Integer}) where d
 end
 
 end
-
-
-#module FD
-#
-#using OffsetArrays
-#include("filtering.jl")
-#
-#
-#"FD-Kernels"
-#module Kernels
-#
-#using OffsetArrays
-#
-#const forward = OffsetVector([-1, 1], 0:1)
-#const backward = OffsetVector([-1, 1], -1:0)
-#const central = OffsetVector([-1, 0, 1], -1:1)
-#
-#end
-#
-#using .Filtering
-#using .Kernels
-#
-#
-#
-### Implementation
-#
-#
-### Grid
-#
-## aka `StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}`
-#const FDStepRange = typeof(0:0.1:1)
-#
-#struct Grid{d}
-#  domain::NTuple{d, FDStepRange}
-#end
-#
-#Grid(domain::FDStepRange...) = Grid(domain)
-#Grid(ndims::Integer, range::FDStepRange) = Grid(ntuple(i -> range, ndims))
-#
-#Base.show(io::IO, grid::Grid) =
-#    print("FDGrid from $(first.(grid.domain)) to $(last.(grid.domain)) ",
-#          "with spacing $(step.(grid.domain))")
-#
-## this is not type-stable?
-##Base.size(grid::Grid) = getproperty.(grid.domain, :len)
-#Base.size(grid::Grid) = ntuple(i -> grid.domain[i].len, ndims(grid))
-#Base.getindex(grid::Grid, idx...) = SVector(getindex.(grid.domain, idx))
-#
-#Base.ndims(grid::Grid{d}) where d = d
-#Base.axes(grid::Grid) = Base.OneTo.(size(grid))
-#Base.length(grid::Grid) = prod(size(grid))
-#
-## TODO: there might be a better way to wrap that product iterator
-#Base.iterate(grid::Grid) = iterate(Iterators.product(grid.domain...))
-#Base.iterate(grid::Grid, state) = iterate(Iterators.product(grid.domain...), state)
-#
-#
-### GridFunction
-#
-#struct GridFunction{T,d} <: AbstractArray{T,1}
-#  grid::Grid{d}
-#  data::Array{T,1}
-#end
-#
-#Base.show(io::IO, grid::GridFunction) =
-#    print("GridFunction on $(length(grid)) grid points")
-#
-#Base.size(f::GridFunction) = size(f.data)
-#Base.getindex(f::GridFunction, idx...) = getindex(f.data, idx...)
-#Base.setindex!(f::GridFunction, v, idx...) = setindex!(f.data, v, idx...)
-#
-#Base.IndexStyle(::Type{<:GridFunction}) = IndexLinear()
-#Base.similar(f::GridFunction, ::Type{S}, size::Int...) where S =
-#    GridFunction(f.grid, Array{S}(undef, size))
-#
-### Operations
-#
-#function fdfilter(f::GridFunction, kern)
-#  A = imfilter(reshape(f.data, size(f.grid)), kern)
-#  GridFunction(f.grid, A[:])
-#end
-#
-#derivate(f::GridFunction, dir::Integer; fd = FDiff.forward) = fdfilter(f, _fdkern(fd, ndims(f.grid), dir))
-#
-##gradient(f::GridFunction) = GridFunction(f.grid,
-#
-#
-#"""
-#    neuman_kernel(kernel, bidx)
-#
-#Modify `kernel` to be usable on a boundary `bixd` in a neumannn boundary setting.
-#
-#Ghost nodes are eliminated through central finite differences. Currently only
-#the basic 5-star Laplace is supported.
-#"""
-#function neumann_kernel(kernel, bidx)
-#  out = copy(kernel)
-#  for (i, bi) in enumerate(bidx)
-#    colrange = (i == j ? Colon() : 0 for j in 1:ndims(kernel))
-#    out[colrange...] .+= bi * FDiff.forward
-#  end
-#  out /= 2 ^ count(!iszero, bidx)
-#  trim_kernel(out, bidx)
-#end
-#
-#function divergence(grid::Grid{d}, field) where d
-#  field = reshape(field, (size(grid)..., ndims(grid)))
-#  size(field, d+1) == d || throw(ArgumentError("need a d-vectorfield on a d-grid"))
-#  dst = zeros(size(field))
-#  for k in 1:d
-#    kern = FDiff.dirfd(FDiff.backward, d, k)
-#    imfilter!(selectdim(dst, d+1, k), selectdim(field, d+1, k), kern)
-#  end
-#  out = dropdims(sum(dst, dims=d+1), dims=d+1)
-#  reshape(out, length(grid))
-#end
-#
-#function gradient(grid::Grid, img::AbstractArray)
-#  img = reshape(img, size(grid))
-#  d = ndims(img)
-#  dst = zeros(size(img)..., d)
-#  for k in 1:d
-#    kern = FDiff.dirfd(FDiff.forward, d, k)
-#    imfilter!(selectdim(dst, d+1, k), img, kern)
-#  end
-#  reshape(dst, length(grid) * ndims(grid))
-#end
-#
-#end
diff --git a/FD/src/filtering.jl b/DualTVDD/src/filtering.jl
similarity index 94%
rename from FD/src/filtering.jl
rename to DualTVDD/src/filtering.jl
index 11fa7b2..21645ca 100644
--- a/FD/src/filtering.jl
+++ b/DualTVDD/src/filtering.jl
@@ -32,11 +32,13 @@ end
 function _trim_kernel_axes(kaxes, bidx)
   ntuple(Val(length(kaxes))) do i
     if bidx[i] < 0
-      return Base.Slice(0:last(kaxes[i]))
+      return 0:last(kaxes[i])
     elseif bidx[i] > 0
-      return Base.Slice(first(kaxes[i]):0)
+      return first(kaxes[i]):0
+    elseif kaxes[i] isa Integer # FIXME: dirty!
+      return UnitRange(kaxes[i], kaxes[i])
     else
-      return Base.Slice(kaxes[i])
+      return UnitRange(kaxes[i])
     end
   end
 end
diff --git a/DualTVDD/test/runtests.jl b/DualTVDD/test/runtests.jl
new file mode 100644
index 0000000..882d689
--- /dev/null
+++ b/DualTVDD/test/runtests.jl
@@ -0,0 +1,26 @@
+using Test
+
+using DualTVDD
+using LinearAlgebra
+
+grid = MetaGrid((0:0.01:1, 0:0.01:1), (2, 2), 5)
+
+A = I
+#A = spdiagm(0 => ones(length(grid)), 1 => ones(length(grid)-1) )
+#f = A * f
+
+f = GridFunction(x -> float(norm(x .- [0.5, 0.5]) < 0.4), 1, grid)
+g = GridFunction(x -> norm(x), 2, grid)
+
+α=100000
+β=0
+τ=1/8
+ϵ=1e-9
+ρ=0.5
+
+p = GridFunction(x -> 0, 2, grid)
+
+for n = 1:10
+  println("iteration $n")
+  dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
+end
diff --git a/FD/Manifest.toml b/FD/Manifest.toml
deleted file mode 100644
index 6453c7c..0000000
--- a/FD/Manifest.toml
+++ /dev/null
@@ -1,65 +0,0 @@
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[DelimitedFiles]]
-deps = ["Mmap"]
-uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-
-[[Distributed]]
-deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-
-[[InteractiveUtils]]
-deps = ["LinearAlgebra", "Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-
-[[LinearAlgebra]]
-deps = ["Libdl"]
-uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-
-[[Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[Mmap]]
-uuid = "a63ad114-7e13-5084-954f-fe012c677804"
-
-[[OffsetArrays]]
-deps = ["DelimitedFiles", "Test"]
-git-tree-sha1 = "f8cd140824f74f2948ff8660bc2292e16d29c1b7"
-uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "0.8.1"
-
-[[Random]]
-deps = ["Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-
-[[SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
-uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
-[[StaticArrays]]
-deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"]
-git-tree-sha1 = "d432c79bef174a830304f8601427a4357dfdbfb7"
-uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "0.8.3"
-
-[[Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/FD/Project.toml b/FD/Project.toml
deleted file mode 100644
index 47f67e5..0000000
--- a/FD/Project.toml
+++ /dev/null
@@ -1,8 +0,0 @@
-name = "FD"
-uuid = "edd53c04-b66a-11e8-3f8d-297d8a5291a5"
-authors = ["Stephan Hilb <stephan.hilb@mathematik.uni-stuttgart.de>"]
-version = "0.1.0"
-
-[deps]
-OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-- 
GitLab