From 3ea29b47f1ce673ae4df1ec4e664bddc004876dc Mon Sep 17 00:00:00 2001
From: Stephan Hilb <stephan@ecshi.net>
Date: Thu, 7 May 2020 14:40:31 +0200
Subject: [PATCH] code: start anew
---
DualTVDD/Manifest.toml | 257 -------------------
DualTVDD/Project.toml | 17 --
DualTVDD/src/DualTVDD.jl | 166 -------------
DualTVDD/src/fd.jl | 450 ---------------------------------
DualTVDD/src/filtering.jl | 121 ---------
DualTVDD/test/runtests.jl | 31 ---
Manifest.toml | 230 -----------------
Project.toml | 18 +-
chambollefp.jl | 99 --------
chambollefp_dd.jl | 114 ---------
fd.jl | 491 -------------------------------------
filtering.jl | 110 ---------
src/DualTVDD.jl | 5 +
test.py | 134 ----------
{DualTVDD => test}/repl.jl | 5 +-
15 files changed, 21 insertions(+), 2227 deletions(-)
delete mode 100644 DualTVDD/Manifest.toml
delete mode 100644 DualTVDD/Project.toml
delete mode 100644 DualTVDD/src/DualTVDD.jl
delete mode 100644 DualTVDD/src/fd.jl
delete mode 100644 DualTVDD/src/filtering.jl
delete mode 100644 DualTVDD/test/runtests.jl
delete mode 100644 Manifest.toml
delete mode 100644 chambollefp.jl
delete mode 100644 chambollefp_dd.jl
delete mode 100644 fd.jl
delete mode 100644 filtering.jl
create mode 100644 src/DualTVDD.jl
delete mode 100644 test.py
rename {DualTVDD => test}/repl.jl (98%)
diff --git a/DualTVDD/Manifest.toml b/DualTVDD/Manifest.toml
deleted file mode 100644
index d95fb7f..0000000
--- a/DualTVDD/Manifest.toml
+++ /dev/null
@@ -1,257 +0,0 @@
-# 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
deleted file mode 100644
index dfc958c..0000000
--- a/DualTVDD/Project.toml
+++ /dev/null
@@ -1,17 +0,0 @@
-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
deleted file mode 100644
index b28348c..0000000
--- a/DualTVDD/src/DualTVDD.jl
+++ /dev/null
@@ -1,166 +0,0 @@
-module DualTVDD
-
-include("filtering.jl")
-include("fd.jl")
-
-# import
-
-using LinearAlgebra: I
-using .FD
-using .FD: parts
-
-# 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)
- if i <= first(r) + metagrid.overlap
- return 1.0
- elseif i >= last(r) - metagrid.overlap
- return 1.0
- else
- return (min(i - first(rp), last(rp) - i, metagrid.overlap) + 1) / (metagrid.overlap + 1)
- end
- 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, α, β, τ)
- # applying B is a global operation
- 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(metagrid, p, A, f; α, β, τ, ϵ, ρ)
-
-Domain decomposition algorithm, calling `fpalg` for each part.
-"""
-function dditer(metagrid, p, A, f; α, β, τ, ϵ, ρ)
- # contribution from previous step
- p̂ = (1-ρ) * p
-
- for part in parts(metagrid)
- grid = view(metagrid, Tuple(part)...)
-
- fv = view(f, Tuple(part)...)
-
- p̂v = view(p̂, Tuple(part)...)
- pv = view(p, Tuple(part)...)
-
- θ = partition_function(metagrid, Tuple(part)...)
- θv = view(θ, Tuple(part)...)
-
- w = FD.GridFunction(metagrid, A' * f.data[:]) +
- FD.divergence(FD.pwmul(FD.GridFunction(metagrid, 1 .- θ), p))
- wv = view(w, Tuple(part)...)
-
- p̂v .+= ρ .* fpalg(grid, w, A, α=α*θv, β=β, τ=τ, ϵ=ϵ)
- end
- p .= p̂
-
- # reconstruct primal solution and plot (FIXME: hack)
- u = p2u(metagrid, p, A, f, α=α, β=β)
- FD.my_plot(u)
- #FD.my_plot(FD.pw_norm(p))
-end
-
-using LinearAlgebra: norm
-using SparseArrays: spdiagm
-"""
- ddalg()
-
-temporary run function
-"""
-function ddalg()
- grid = MetaGrid((0:0.01:1, 0:0.01:1), (1, 1), 5)
-
- f = GridFunction(x -> float(norm(x .- [0.5, 0.5]) < 0.4), 1, grid)
- g = GridFunction(x -> norm(x), 2, grid)
-
- A = I
- #A = spdiagm(0 => ones(length(grid)), 1 => 0.3 * ones(length(grid)-1))
- #f = GridFunction(grid, A * f.data[:])
-
- "TV parameter"
- α=1
- "L2 parameter"
- β=0
- "chambolle weighing"
- τ=1/8
- "residual error tolerance for local problem"
- ϵ=1e-9
- "???"
- ρ=0.5
-
- p = GridFunction(x -> 0, 2, grid)
-
- for n = 1:10
- println("iteration $n")
- dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
- end
-end
-
-end # module
diff --git a/DualTVDD/src/fd.jl b/DualTVDD/src/fd.jl
deleted file mode 100644
index 6829c9e..0000000
--- a/DualTVDD/src/fd.jl
+++ /dev/null
@@ -1,450 +0,0 @@
-module FD
-
-# local modules
-"Finite Differences"
-module FDiff
-
-using OffsetArrays: OffsetVector
-
-const forward = OffsetVector([-1, 1], -1)
-const backward = OffsetVector([-1, 1], -2)
-const central = OffsetVector([-1, 0, 1], -2)
-
-end # module FDiff
-
-
-# imports
-using StaticArrays: SVector
-using OffsetArrays
-import Colors: Gray
-import Plots: plot
-# local imports
-using ..Filtering
-using .FDiff
-
-# export
-export Grid, MetaGrid, GridView, GridFunction
-export grid
-export ⊙
-
-
-# abstract types
-"""
- AbstractGrid{d}
-
-Abstract rectangular grid type, allowing positional indexing of coordinates.
-"""
-abstract type AbstractGrid{d} <: AbstractArray{SVector{d,Float64}, d} end
-
-Base.IndexStyle(::Type{<:AbstractGrid}) = IndexCartesian()
-
-Base.parent(grid::AbstractGrid) = grid
-
-
-"""
- _fdkern(fd, ndims, dir)
-
-Create from `fd` a finite difference kernel usable in `ndims` dimension
-oriented in dimension `dir`.
-"""
-# TODO: make it typestable on ndims
-function _fdkern(fd::OffsetArray{T,1}, ndims::Integer, dir::Integer) where T
- 1 <= dir <= ndims || throw(ArgumentError("dimesion $(dir) out of range 1:$(ndims)"))
-
- offsets = ntuple(i -> i == dir ? axes(fd, 1) : (0:0), ndims)
- out = OffsetArray{T}(undef, offsets)
- # copyto! from OffsetArray fails, broadcasting for differing axes fails too
- # workaround: copyto! from parent
- copyto!(out, fd.parent)
-end
-
-
-# Grid
-
-
-# aka `StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}`
-const FDStepRange = typeof(0:0.1:1)
-
-struct Grid{d} <: AbstractGrid{d}
- domain::NTuple{d, FDStepRange}
-end
-Grid(domain::FDStepRange...) = Grid(domain)
-
-# TODO: somehow display(grid) is incorrect for d = 1 ?
-Base.show(io::IO, grid::Grid) =
- print(io, "Grid from $(first.(grid.domain)) to $(last.(grid.domain)) ",
- "with spacing $(step.(grid.domain))")
-
-Base.size(grid::Grid) = ntuple(i -> grid.domain[i].len, ndims(grid))
-
-Base.iterate(grid::Grid) = iterate(Iterators.product(grid.domain...))
-Base.iterate(grid::Grid, state) = iterate(Iterators.product(grid.domain...), state)
-
-Base.getindex(grid::Grid, idx...) = SVector(getindex.(grid.domain, idx))
-Base.getindex(grid::Grid, I::CartesianIndex) = SVector(getindex.(grid.domain, Tuple(I)))
-
-
-# GridView
-
-struct GridView{d} <: AbstractGrid{d}
- grid::AbstractGrid{d}
- indices::NTuple{d, UnitRange{Int}}
-end
-GridView(grid, indices::UnitRange{Int}...) = GridView(grid, indices)
-
-Base.size(grid::GridView) = length.(grid.indices)
-
-# FIXME: this is terrible!
-Base.iterate(grid::GridView) = iterate(grid.grid[I] for I in CartesianIndices(grid.indices))
-Base.iterate(grid::GridView, state) = iterate([grid.grid[I] for I in CartesianIndices(grid.indices)], state)
-
-Base.getindex(grid::GridView, idx...) = getindex(grid.grid, (first.(grid.indices) .- 1 .+ Base.to_indices(grid, idx))...)
-
-Base.parent(grid::GridView) = grid.grid
-
-
-# MetaGrid (partitioned Grid)
-
-
-struct MetaGrid{d} <: AbstractGrid{d}
- grid::Grid{d}
- indices::Array{NTuple{d, UnitRange{Int}}, d}
- overlap::Int
-end
-
-MetaGrid(domain::NTuple{d, FDStepRange}, pnum, overlap) where d = MetaGrid(Grid(domain), pnum, overlap)
-
-"""
- MetaGrid(domain::NTuple{d, FDStepRange}, pnum::NTuple{d, Int}, overlap::Int)
- MetaGrid(grid::Grid{d}, pnum::NTuple{d, Int}, overlap::Int)
-
-Create a grid on `domain`, partitioned along dimensions according to `pnum`
-respecting `overlap`.
-"""
-function MetaGrid(grid::Grid{d}, pnum::NTuple{d, Int}, overlap::Int) where d
- tsize = size(grid) .+ (pnum .- 1) .* overlap
-
- psize = tsize .÷ pnum
- osize = tsize .- pnum .* psize
-
- overhang(I, j) = I[j] == pnum[j] ? osize[j] : 0
-
- indices = Array{NTuple{d, UnitRange{Int}}, d}(undef, pnum)
- for I in CartesianIndices(pnum)
- indices[I] = ntuple(j -> ((I[j] - 1) * psize[j] - (I[j] - 1) * overlap + 1) :
- ( I[j] * psize[j] - (I[j] - 1) * overlap + overhang(I, j)), d)
- end
- MetaGrid{d}(grid, indices, overlap)
-end
-
-Base.size(grid::MetaGrid) = size(grid.grid)
-
-Base.iterate(grid::MetaGrid) = iterate(grid.grid)
-Base.iterate(grid::MetaGrid, state) = iterate(grid.grid, state)
-
-Base.getindex(grid::MetaGrid, idx...) = getindex(grid.grid, idx...)
-
-indices(grid::MetaGrid, idx...) = getindex(grid.indices, idx...)
-
-Base.view(grid::MetaGrid, idx...) = GridView(grid, indices(grid, idx...))
-
-parts(grid::MetaGrid) = CartesianIndices(grid.indices)
-
-
-# GridFunction
-
-
-abstract type AbstractGridFunction{G,T,d,n} <: AbstractArray{T,1} end
-
-rangedim(f::AbstractGridFunction{G,T,d,n}) where {G,T,d,n} = n
-domaindim(f::AbstractGridFunction{G,T,d}) where {G,T,d} = d
-
-Base.similar(f::F, ::Type{S}, dims::Int...) where {F<:AbstractGridFunction, S} =
- GridFunction(f.grid, Array{S}(undef, length(f)))
-Base.parent(f::AbstractGridFunction) = f
-
-
-"""
- GridFunction{G<:AbstractGrid,T,d,n,D}
-
-`n`-dimensional function values of type `T` on a Grid `G` of dimension `d`.
-
-`D == d + 1`
-"""
-struct GridFunction{G<:AbstractGrid,T,d,n,A,D} <: AbstractGridFunction{G,T,d,n}
- grid::G
- data::A
- function GridFunction{G,T,d,n,A,D}(grid::G, data::A) where {G,T,d,n,D,A<:AbstractArray{T,D}}
- D == d + 1 || throw(ArgumentError)
- d == ndims(grid) || throw(ArgumentError)
- (size(grid)..., n) == size(data) || throw(ArgumentError("data does not match given size"))
- new(grid, data)
- end
-end
-# TODO: check size(data)
-#GridFunction(grid::G, data::AbstractArray{T,D}) where {d,G<:AbstractGrid{d},T,D} =
-# GridFunction{G,T,d,size(data, d+1),typeof(data),d+1}(grid, data)
-function GridFunction(grid::G, data::AbstractArray, ::Val{n}) where {d,G<:AbstractGrid{d}, n}
- if length(data) == length(grid) && n > 1
- data = repeat(data, inner=ntuple(i -> i == d+1 ? n : 1, d+1))
- elseif size(data) != (size(grid)..., n)
- data = reshape(data, size(grid)..., n)
- end
- length(data) ÷ length(grid) == n || throw(ArgumentError("dimensions don't match"))
- GridFunction{G,eltype(data),d,n,typeof(data),d+1}(grid, data)
-end
-function GridFunction(grid::G, data::AbstractArray) where {d,G<:AbstractGrid{d}}
- GridFunction(grid, data, Val(length(data) ÷ length(grid)))
-end
-function GridFunction(f::Function, n::Int, grid::G) where {G<:AbstractGrid}
- data = Array{Float64,ndims(grid)+1}(undef, size(grid)..., n)
- for I in CartesianIndices(grid)
- data[Tuple(I)..., :] .= f(grid[I])
- end
- GridFunction(grid, data)
-end
-function GridFunction(grid::AbstractGrid, ::UndefInitializer, ::Val{n} = Val(1)) where n
- data = Array{Float64}(undef, size(grid)..., 1)
- GridFunction{typeof(grid), eltype(data), ndims(grid), n, typeof(data), ndims(grid) + 1}(grid, data)
-end
-
-Base.show(io::IO, f::GridFunction) =
- print(io, "$(rangedim(f))d-GridFunction on $(f.grid), $(length(f)) dofs")
-
-Base.size(f::GridFunction) = (length(f.data),)
-#Base.axes(f::GridFunction) = (Base.OneTo(length(f.data)),)
-Base.getindex(f::GridFunction, idx::Int) = getindex(f.data, idx)
-Base.setindex!(f::GridFunction, v, idx::Int) = setindex!(f.data, v, idx)
-
-## TODO: using this prints #undef for display(f)
-#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()
-
-#component(f::GridFunction, k::Integer) = view(reshape(f.data,
-grid(f::GridFunction) = f.grid
-
-# TODO: call is different, so maybe make thin not part of Base
-Base.view(f::AbstractGridFunction{G}, idx...) where G<:MetaGrid =
-GridFunction(GridView(f.grid, indices(f.grid, idx...)), view(f.data, indices(f.grid, idx...)..., :))
-
-function zero_continuation(f::GridFunction{<:GridView})
- dst = GridFunction(parent(grid(f)), zeros(size(parent(grid(f)))..., rangedim(f)), Val(rangedim(f)))
- dst.data[grid(f).indices..., :] .= f.data
- dst
-end
-
-Base.parent(f::GridFunction{<:GridView,<:Any,<:Any,<:Any,<:SubArray}) =
-GridFunction(parent(grid(f)), parent(f.data))
-
-
-
-
-
-function _check_shape(a::GridFunction, b::GridFunction)
- #TODO: why does fv.grid == gv.grid fail?
- a.grid === b.grid || throw(ArgumentError("grids need to match"))
- rangedim(a) == rangedim(b) || throw(ArgumentError("range dimensions need to match"))
-end
-
-Base.:*(a::Number, b::GridFunction) = GridFunction(b.grid, a * b.data, Val(rangedim(b)))
-Base.:-(a::GridFunction) = GridFunction(a.grid, -a.data, Val(rangedim(a)))
-function Base.:+(a::GridFunction, b::GridFunction)
- _check_shape(a, b)
- GridFunction(a.grid, a.data .+ b.data, Val(rangedim(a)))
-end
-function Base.:-(a::GridFunction, b::GridFunction)
- _check_shape(a, b)
- GridFunction(a.grid, a.data .- b.data, Val(rangedim(a)))
-end
-
-
-function Base.:+(a::GridFunction, b::GridFunction{<:GridView})
- a.grid == parent(b.grid) || throw(ArgumentError("grids need to match"))
- data = copy(a.data)
- datav = view(data, b.grid.indices..., :)
- datav .+= b.data
- GridFunction(a.grid, data, Val(rangedim(a)))
-end
-Base.:+(a::GridFunction{<:GridView}, b::GridFunction) = b + a
-
-# TODO: cleanup! use metaprogramming
-function Base.:+(a::GridFunction{<:GridView}, b::GridFunction{<:GridView})
- _check_shape(a, b)
- GridFunction(a.grid, a.data .+ b.data, Val(rangedim(a)))
-end
-# TODO: cleanup! use metaprogramming
-function Base.:-(a::GridFunction{<:GridView}, b::GridFunction{<:GridView})
- _check_shape(a, b)
- GridFunction(a.grid, a.data .- b.data, Val(rangedim(a)))
-end
-
-function Base.:-(a::GridFunction, b::GridFunction{<:GridView})
- a.grid == parent(b.grid) || throw(ArgumentError("grids need to match"))
- data = copy(a.data)
- datav = view(data, b.grid.indices..., :)
- datav .-= b.data
- GridFunction(a.grid, data, Val(rangedim(a)))
-end
-Base.:-(a::GridFunction{<:GridView}, b::GridFunction) = b - a
-
-
-
-function my_plot(img::FD.GridFunction{G,T,d,1}) where {G,T,d}
- print("Plotting ... ")
-
- data = reshape(copy(img.data), size(img.grid))
-
- min = minimum(data)
- max = maximum(data)
- @show min, max
-
- #@. data = (data - min) / (max - min)
- @. data = clamp(data, 0, 1)
-
- plot(Gray.(data), show=true)
-end
-
-
-## Operations
-
-function fdfilter(f::GridFunction{G,T,d,n}, kern) where {G,T,d,n}
- A = zeros(axes(f.data))
- for k = 1:n
- vidx = ntuple(i -> i <= d ? Colon() : k, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx...], kern)
- end
- GridFunction(f.grid, A)
-end
-
-derivate(f::GridFunction, dir::Integer; fd = FDiff.forward) = fdfilter(f, _fdkern(fd, ndims(f.grid), dir))
-
-
-# currently gradient and divergence are specific implementations (as needed for
-# chambolles algorithm), other choices may be reasonable
-
-function gradient(f::GridFunction{G,T,d,1}; fd=FDiff.forward) where {G,T,d}
- #A = Array{T,d+1}(undef, (size(f.grid)..., d))
- A = zeros(size(f.grid)..., d)
- for k = 1:d
- vidx = ntuple(i -> i <= d ? Colon() : k, d + 1)
- vidx2 = ntuple(i -> i <= d ? Colon() : 1, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx))
- # FIXME: hackish, only for forward fd
- selectdim(selectdim(A, d+1, k), k, size(f.grid, k)) .= 0
- end
- GridFunction(f.grid, A)
-end
-
-function divergence(f::GridFunction{G,T,d}; fd=FDiff.backward) where {d,G<:AbstractGrid{d},T}
- A = zeros(size(f.grid)..., 1)
- for k = 1:d
- vidx = ntuple(i -> i <= d ? Colon() : 1, d + 1)
- vidx2 = ntuple(i -> i <= d ? Colon() : k, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx))
- end
- GridFunction(f.grid, A)
-end
-
-function _kernel_zero_right_bnd!(kernel, bidx, k)
- if bidx[k] == 1
- kernel[ntuple(i -> 0, ndims(kernel))...] = 0
- end
- kernel
-end
-
-function gradient_c(f::GridFunction{G,T,d,1}; fd=FDiff.forward) where {G,T,d}
- #A = Array{T,d+1}(undef, (size(f.grid)..., d))
- A = zeros(size(f.grid)..., d)
- for k = 1:d
- vidx = ntuple(i -> i <= d ? Colon() : k, d + 1)
- vidx2 = ntuple(i -> i <= d ? Colon() : 1, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx))
- end
- GridFunction(f.grid, A)
-end
-
-function divergence_z(f::GridFunction{G,T,d}; fd=FDiff.backward) where {d,G<:AbstractGrid{d},T}
- A = zeros(size(f.grid)..., 1)
- for k = 1:d
- vidx = ntuple(i -> i <= d ? Colon() : 1, d + 1)
- vidx2 = ntuple(i -> i <= d ? Colon() : k, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> _kernel_zero_right_bnd!(trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx), bidx, k))
- end
- GridFunction(f.grid, A)
-end
-
-function pw_norm(f::GridFunction)
- d = domaindim(f)
- GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(1))
-end
-
-function my_norm(f::GridFunction)
- d = domaindim(f)
- GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(rangedim(f)))
-end
-
-"""
- f ⊙ g
-
-Pointwise multiplication of two gridfunctions.
-"""
-⊙(f, g) = pwmul(f, g)
-pwmul(f::GridFunction{G,T,d,n}, g::GridFunction{G,T,d,n}) where {G,T,d,n} = f .* g
-function pwmul(f::GridFunction{G,T,d,1}, g::GridFunction{G,T,d}) where {G,T,d}
- # TODO: maybe comprehension? Or make broadcast work
- data = similar(g.data)
- for I in CartesianIndices(f.grid)
- data[Tuple(I)..., :] = f.data[Tuple(I)..., 1] * g.data[Tuple(I)..., :]
- end
- GridFunction(f.grid, data, Val(rangedim(g)))
-end
-
-"""
- 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 my_norm(grid::Grid, img::AbstractArray)
- img = reshape(img, (size(grid)..., ndims(grid)))
- d = ndims(grid)
-
- out = sqrt.(sum(img.^2, dims=d+1))
- out = repeat(out, inner=ntuple(i -> i == d+1 ? size(img, d+1) : 1, d+1))
-
- out[:]
-end
-
-
-function apply(kern, flatimg, size)
- img = reshape(flatimg, size)
- out = imfilter(kern, img)
- reshape(out, prod(size))
-end
-
-function solve(pde, shape::NTuple{d, Integer}) where d
- ndims(pde.stencil) == d || throw("dimension mismatch")
- length(pde.f) == prod(shape) || throw("bounds mismatch")
-
- A = LinearMap{Float64}(u -> apply(pde.stencil, u, shape), prod(shape), issymmetric=true, isposdef=true)
-
- u = zeros(prod(shape))
- cg!(u, A, pde.f)
- return u
-end
-
-end
diff --git a/DualTVDD/src/filtering.jl b/DualTVDD/src/filtering.jl
deleted file mode 100644
index e9adaeb..0000000
--- a/DualTVDD/src/filtering.jl
+++ /dev/null
@@ -1,121 +0,0 @@
-module Filtering
-
-import OffsetArrays
-
-export imfilter, imfilter!, trim_kernel
-
-const Kernel{T, N} = OffsetArrays.OffsetArray{T, N}
-
-@inline function _shift_range(range, shift::Integer)
- Base.UnitRange(first(range)+shift : last(range)+shift)
-end
-
-## TODO: allow border to be tuple
-#function _innershift_indices(a::Array{T, d}, shift::NTuple{d, Integer}, border = abs(maximum(shift))) where {T, d}
-# abs(maximum(shift)) <= border || throw("cannot shift by $shift with border $border")
-# range(i) = (1 + shift[i] + border):(lastindex(a, i) + shift[i] - border)
-# ntuple(range, Val(ndims(a)))
-#end
-#
-#function innershift_view(a::Array{T, d}, shift::NTuple{d, Integer} = ntuple(i -> 0, d), border = abs(maximum(shift))) where {T, d}
-# view(a, _innershift_indices(a, shift, border)...)
-#end
-
-@inline function _bidx_to_range(range::AbstractUnitRange, bidx::Integer)
- if bidx > 0
- return last(range):last(range)
- elseif bidx < 0
- return first(range):first(range)
- else
- return first(range)+1:last(range)-1
- end
-end
-
-
-function _trim_kernel_axes(kaxes, bidx)
- ntuple(Val(length(kaxes))) do i
- if bidx[i] < 0
- return 0:last(kaxes[i])
- elseif bidx[i] > 0
- return first(kaxes[i]):0
- elseif kaxes[i] isa Integer # FIXME: dirty!
- return UnitRange(kaxes[i], kaxes[i])
- else
- return UnitRange(kaxes[i])
- end
- end
-end
-
-"""
- trim_kernel(kernel, bidx)
-
-Trim `kernel` so that the resulting kernel is usable on the boundary indicated
-by `bidx`, i.e. replace evaluations that would occur outside the boundary with
-zero.
-"""
-function trim_kernel(kernel, bidx)
- dst_axes = _trim_kernel_axes(axes(kernel), bidx)
- out = similar(kernel, dst_axes)
- range = CartesianIndices(out)
- copyto!(out, range, kernel, range)
-end
-
-
-"""
- imfilter!(dst, img, kern, slice=:inner)
-
-Filters `img` applying `kern` and write to `dst`. `dst` is added to, so you
-should initialize with zeros.
-Only `dst[slice]` is written to.
-"""
-# TODO: weirdly Base.Slice instead of Base.UnitRange leads to axes that are not
-# compatible when broadcasting, figure out why!
-function imfilter!(dst, img, kern, slice=_bidx_to_range.(axes(img), ntuple(i->0, ndims(img))))
- ndims(dst) == ndims(img) == ndims(kern) == length(slice) ||
- throw(ArgumentError("dst, image, kernel or slice dimensions disagree " *
- "($(ndims(dst)) vs $(ndims(img)) vs $(ndims(kern)) vs $(length(slice)))"))
- axes(dst) == axes(img) ||
- throw(ArgumentError("output and input image axes disagree"))
- all(size(kern) .<= size(img)) ||
- throw(ArgumentError("kernel bigger than image"))
- last.(slice) .+ last.(axes(kern)) <= last.(axes(img)) ||
- first.(slice) .+ first.(axes(kern)) >= first.(axes(img)) ||
- throw(ArgumentError("kernel and/or slice out of bounds for image"))
-
- dstslice = view(dst, UnitRange.(slice)...)
- for i in CartesianIndices(kern)
- iszero(kern[i]) && continue
- dstslice .+= kern[i] .* view(img, (_shift_range.(slice, Tuple(i)))...)
- end
- dst
-end
-
-function imfilter!(dst, img, kernelf::Function)
- for bidx in Iterators.product(ntuple(i -> -1:1, ndims(img))...)
- range = _bidx_to_range.(axes(img), bidx)
- imfilter!(dst, img, kernelf(bidx), range)
- end
- dst
-end
-
-"""
- imfilter(img, kernelf)
-
-Filter `img` using a kernel function `kernelf(bidx) = kernel` specifying a
-kernel for each type of node.
-"""
-# need zeros here, because `imfilter!` acts additively
-imfilter(img, kernelf::Function) = imfilter!(zeros(size(img)), img, kernelf)
-
-
-"""
- imfilter(img, kern)
-
-Filter `img` by `kern` using correlation.
-Boundary is zeroed out.
-"""
-imfilter(img, kern) = imfilter(img, bidx -> trim_kernel(kern, bidx))
-# TODO: use keyword arguments
-imfilter_inner(img, kern) = imfilter(img, bidx -> all(bidx .== 0) ? kern : zeros(ntuple(i -> 0, ndims(img))))
-
-end
diff --git a/DualTVDD/test/runtests.jl b/DualTVDD/test/runtests.jl
deleted file mode 100644
index 74a4996..0000000
--- a/DualTVDD/test/runtests.jl
+++ /dev/null
@@ -1,31 +0,0 @@
-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)
-
-"TV parameter"
-α=100000
-"L2 parameter"
-β=0
-"???"
-τ=1/8
-"residual error tolerance for local problem"
-ϵ=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/Manifest.toml b/Manifest.toml
deleted file mode 100644
index a9e8b92..0000000
--- a/Manifest.toml
+++ /dev/null
@@ -1,230 +0,0 @@
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[ColorTypes]]
-deps = ["FixedPointNumbers", "Random", "Test"]
-git-tree-sha1 = "f73b0e10f2a5756de7019818a41654686da06b09"
-uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.7.5"
-
-[[Colors]]
-deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Pkg", "Printf", "Reexport", "Test"]
-git-tree-sha1 = "8c89e0a9a583954eae3efcf6a531e51c02b38cee"
-uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
-version = "0.9.4"
-
-[[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 = "ff2595695fc4f14427358ce2593f867085c45dcb"
-uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "1.2.0"
-
-[[Contour]]
-deps = ["LinearAlgebra", "StaticArrays", "Test"]
-git-tree-sha1 = "b974e164358fea753ef853ce7bad97afec15bb80"
-uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
-version = "0.5.1"
-
-[[DataStructures]]
-deps = ["InteractiveUtils", "OrderedCollections", "REPL", "Random", "Serialization", "Test"]
-git-tree-sha1 = "8fc6e166e24fda04b2b648d4260cdad241788c54"
-uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.14.0"
-
-[[Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-
-[[DelimitedFiles]]
-deps = ["Mmap"]
-uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-
-[[Distributed]]
-deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-
-[[FD]]
-deps = ["OffsetArrays", "StaticArrays"]
-path = "FD"
-uuid = "edd53c04-b66a-11e8-3f8d-297d8a5291a5"
-version = "0.1.0"
-
-[[FixedPointNumbers]]
-deps = ["Pkg", "Test"]
-git-tree-sha1 = "b8045033701c3b10bf2324d7203404be7aef88ba"
-uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
-version = "0.5.3"
-
-[[GR]]
-deps = ["Base64", "DelimitedFiles", "Pkg", "Random", "Serialization", "Sockets", "Test"]
-git-tree-sha1 = "0437921cbf8ee555ab7a2a3115c1d0907289e816"
-uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.34.1"
-
-[[InteractiveUtils]]
-deps = ["LinearAlgebra", "Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[JSON]]
-deps = ["Dates", "Distributed", "Mmap", "Pkg", "Sockets", "Test", "Unicode"]
-git-tree-sha1 = "fec8e4d433072731466d37ed0061b3ba7f70eeb9"
-uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.19.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 = ["Pkg", "Test"]
-git-tree-sha1 = "ddfd6d13e330beacdde2c80de27c1c671945e7d9"
-uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
-version = "0.3.0"
-
-[[Missings]]
-deps = ["Dates", "InteractiveUtils", "SparseArrays", "Test"]
-git-tree-sha1 = "adc26d2ee85a49c413464110d922cf21efc9d233"
-uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
-version = "0.3.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]]
-deps = ["DelimitedFiles", "Test"]
-git-tree-sha1 = "f8cd140824f74f2948ff8660bc2292e16d29c1b7"
-uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "0.8.1"
-
-[[OrderedCollections]]
-deps = ["Pkg", "Random", "Serialization", "Test"]
-git-tree-sha1 = "85619a3f3e17bb4761fe1b1fd47f0e979f964d5b"
-uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.0.2"
-
-[[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 = "78553a920c4869d20742ba66193e3692369086ec"
-uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
-version = "0.5.4"
-
-[[Plots]]
-deps = ["Base64", "Contour", "Dates", "FixedPointNumbers", "GR", "JSON", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Printf", "Random", "RecipesBase", "Reexport", "Requires", "Showoff", "SparseArrays", "StaticArrays", "Statistics", "StatsBase", "Test", "UUIDs"]
-git-tree-sha1 = "437b77558cc7f24883bd4c4fe25a240304402ffd"
-uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-version = "0.20.5"
-
-[[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]]
-deps = ["Pkg", "Random", "Test"]
-git-tree-sha1 = "0b3cb370ee4dc00f47f1193101600949f3dcf884"
-uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-version = "0.6.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 = ["Compat", "Pkg"]
-git-tree-sha1 = "276b24f3ace98bec911be7ff2928d497dc759085"
-uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
-version = "0.2.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 = ["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"
-
-[[StatsBase]]
-deps = ["DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "Test"]
-git-tree-sha1 = "723193a13e8078cec6dcd0b8fe245c8bfd81690e"
-uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.25.0"
-
-[[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[[UUIDs]]
-deps = ["Random"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-
-[[Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
diff --git a/Project.toml b/Project.toml
index 9add2e9..b60374f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,4 +1,14 @@
-[deps]
-Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
-FD = "edd53c04-b66a-11e8-3f8d-297d8a5291a5"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+name = "DualTVDD"
+uuid = "93adc0ee-851f-4b8b-8bf8-c8a87ded093b"
+authors = ["Stephan Hilb <stephan@ecshi.net>"]
+version = "0.1.0"
+
+[compat]
+julia = "1.4.0"
+
+[extras]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test", "BenchmarkTools"]
diff --git a/chambollefp.jl b/chambollefp.jl
deleted file mode 100644
index 56f3495..0000000
--- a/chambollefp.jl
+++ /dev/null
@@ -1,99 +0,0 @@
-using StaticArrays
-using SparseArrays
-using OffsetArrays
-using LinearAlgebra
-using Colors
-using Plots
-
-using FD
-
-
-imf(f, grid) = [f(x) for x in grid]
-
-imshow(grid, img) = FD.my_plot(FD.GridFunction(grid, img))
-#plot(Gray.((reshape(float.(img), size(grid)) .- minimum(float.(img))) / (maximum(float.(img)) - minimum(float.(img))) ), show=true)
-
-function _fpiter(grid, p, f, A, α, β, τ)
- x = (A'*A + β*I) \ (FD.divergence(grid, p) .+ A'*f)
- q = FD.gradient(grid, x)
-
- p_new = (p .+ τ .* q) ./ (1 .+ (τ/α) .* FD.my_norm(grid, q))
-
- reserr = maximum(sum((p_new .- p).^2, dims=ndims(p)))
- (p_new, reserr)
-end
-
-function fpalg(grid, f, A; α, β, τ = 1/8, ϵ = 1e-4)
- p = zeros(length(grid) * ndims(grid))
-
- reserr = 1
- k = 0
- while reserr > ϵ && k < 300e1
- k += 1
- #sleep(0.1 / k)
- (p, reserr) = _fpiter(grid, p, f, A, α, β, τ)
-
- uimg = (A'*A + β*I) \ (FD.divergence(grid, p) .+ A'*f)
-
- imshow(grid, uimg)
- println("[$(lpad(k, 4))] res = $(round(reserr, digits=5))")
- println("min/max: $(minimum(FD.my_norm(FD.GridFunction(grid,p)))) / $(maximum(FD.my_norm(FD.GridFunction(grid,p))))")
- end
-end
-
-#st = OffsetArray([0 -1 0; -1 4 -1; 0 -1 0], (-2, -2))
-
-grid = FD.Grid(0:0.01:1, 0:0.01:1)
-#grid = FD.Grid(0:0.5:1, 0:0.5:1)
-
-n = 1
-#A = spdiagm(collect(k => fill(1/(2n+1), length(grid)-abs(k)) for k in -n:n)...)
-
-# circle, scalar
-f = reshape([float(norm(x .- [0.5, 0.5]) < 0.4) for x in grid], length(grid))
-#f += clamp.(0.2 * rand(eltype(f), size(f)), 0, 1)
-
-
-A = I
-#A = spdiagm(0 => ones(length(grid)), 1 => ones(length(grid)-1) )
-
-f = A * f
-
-#imshow(grid, f)
-#sleep(5)
-
-img = fpalg(grid, f, A, α=10000, β=0, τ=1/8, ϵ=1e-7)
-
-
-#A = Array{Float64}(undef, (length(grid), length(grid)))
-#for i in 1:length(grid)
-# ei = float.(1:length(grid) .== i)
-# #A[:, i] = FD.divergence(grid, FD.gradient(grid, ei))
-# #A[:, i] = FD.divergence(grid, append!(ei, zeros(length(grid))))
-# A[:, i] = reshape(reshape(FD.gradient(grid, ei), size(grid)..., 2)[:,:,2], length(grid))
-#end
-#A
-
-#display(A)
-
-
-#display(f2)
-#imshow(grid, f2)
-
-
-
-#@benchmark solve((stencil=st, f=zeros(10000)), (100, 100))
-
-#u = zeros(100, 100)
-
-#A = LinearMap{Float64}(u -> apply(st, u, size(uinner)), length(uinner), issymmetric=true, isposdef=true)
-
-#cg!(uflat, A, zeros(length(uinner)))
-
-
-#A = sparse(A)
-
-#@benchmark cg($A, zeros(10000))
-
-# TODO: A* as sparse matrix?
-
diff --git a/chambollefp_dd.jl b/chambollefp_dd.jl
deleted file mode 100644
index b622ff9..0000000
--- a/chambollefp_dd.jl
+++ /dev/null
@@ -1,114 +0,0 @@
-#include("fd.jl")
-
-using LinearAlgebra
-using SparseArrays
-import FD
-
-
-import Colors: Gray
-import Plots: plot
-
-
-
-"partition of unity function"
-function theta(grid, idx...)
- function h(r, rp, i)
- i <= first(r) + grid.overlap && return 1.0
- i >= last(r) - grid.overlap && return 1.0
- (min(i - first(rp), last(rp) - i, grid.overlap) + 1) / (grid.overlap + 1)
- end
-
- vidx = FD.indices(grid, idx...)
- ctheta = map(x -> [h(x[1], x[2], i) for i in x[2]], zip(axes(grid), vidx))
-
- dst = zeros(size(grid))
- for (I,J) in zip(CartesianIndices(vidx), CartesianIndices(length.(vidx)))
- dst[I] = prod(getindex.(ctheta, Tuple(J)))
- end
- FD.GridFunction(grid, dst)
-end
-
-
-function p2u(grid, p, A, f; α, β)
- FD.GridFunction(grid, (A'*A + β*I) \ (FD.divergence_z(p).data[:] .+ A'*f.data[:]))
-end
-
-
-# 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
-
-
-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
-
-
-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)...)
-
- θ = theta(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
-
-
-grid = FD.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 = FD.GridFunction(x -> float(norm(x .- [0.5, 0.5]) < 0.4), 1, grid)
-g = FD.GridFunction(x -> norm(x), 2, grid)
-
-α=100000
-β=0
-τ=1/8
-ϵ=1e-9
-ρ=0.5
-
-p = FD.GridFunction(x -> 0, 2, grid)
-
-for n = 1:10
- println("iteration $n")
- dditer(grid, p, A, f, α=α, β=β, τ=τ, ϵ=ϵ, ρ=ρ)
-end
diff --git a/fd.jl b/fd.jl
deleted file mode 100644
index 19d3c07..0000000
--- a/fd.jl
+++ /dev/null
@@ -1,491 +0,0 @@
-module FDM
-
-using StaticArrays
-using SparseArrays
-using OffsetArrays
-using LinearMaps
-using IterativeSolvers
-
-include("filtering.jl")
-
-"Finite Differences"
-module FD
-
-using OffsetArrays
-
-const forward = OffsetVector([-1, 1], -1)
-const backward = OffsetVector([-1, 1], -2)
-const central = OffsetVector([-1, 0, 1], -2)
-function dirfd(fd, ndims, dir)
- dir <= ndims || throw(ArgumentError("need dir <= ndims!"))
- dfd = OffsetArray(zeros(ntuple(i -> i == dir ? length(fd) : 1, ndims)), ntuple(i -> i == dir ? fd.offsets[1] : -1, ndims))
- # Workaround: LinearIndices(::OffsetVector) are not 1-based, so copyto! chockes
- copyto!(dfd, fd.parent)
-end
-
-end
-
-
-using .Filtering
-using .FD
-
-"""
- _fdkern(fd, ndims, dir)
-
-Create from `fd` a finite difference kernel usable in `ndims` dimension
-oriented in dimension `dir`.
-"""
-# TODO: make it typestable on ndims
-function _fdkern(fd::OffsetArray{T,1,Array{T,1}}, ndims::Integer, dir::Integer) where T
- 1 <= dir <= ndims || throw(ArgumentError("dimesion $(dir) out of range 1:$(ndims)"))
- dfd = OffsetArray{T}(undef, Tuple(i == dir ? axes(fd, 1) : (0:0) for i in 1:ndims))
- # copyto! from OffsetArray fails, broadcasting for differing axes fails too
- # workaround: copyto! from parent
- copyto!(dfd, fd.parent)
-end
-
-
-# Grid
-
-
-abstract type AbstractGrid{d} <: AbstractArray{Float64, d} end
-
-Base.ndims(grid::AbstractGrid{d}) where d = d
-Base.axes(grid::AbstractGrid) = Base.OneTo.(size(grid))
-Base.length(grid::AbstractGrid) = prod(size(grid))
-Base.IndexStyle(::Type{<:AbstractGrid}) = IndexCartesian()
-
-# aka `StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}`
-const FDStepRange = typeof(0:0.1:1)
-
-struct Grid{d} <: AbstractGrid{d}
- domain::NTuple{d, FDStepRange}
-end
-Grid(domain::FDStepRange...) = Grid(domain)
-Base.show(io::IO, grid::Grid) =
- print("Grid from $(first.(grid.domain)) to $(last.(grid.domain)) ",
- "with spacing $(step.(grid.domain))")
-
-Base.size(grid::Grid) = ntuple(i -> grid.domain[i].len, ndims(grid))
-
-Base.iterate(grid::Grid) = iterate(Iterators.product(grid.domain...))
-Base.iterate(grid::Grid, state) = iterate(Iterators.product(grid.domain...), state)
-
-Base.getindex(grid::Grid, idx...) = SVector(getindex.(grid.domain, idx))
-Base.getindex(grid::Grid, I::CartesianIndex) = SVector(getindex.(grid.domain, Tuple(I)))
-
-
-# GridView
-
-struct GridView{d} <: AbstractGrid{d}
- grid::AbstractGrid{d}
- indices::NTuple{d, UnitRange{Int}}
-end
-GridView(grid, indices::UnitRange{Int}...) = GridView(grid, indices)
-
-Base.size(grid::GridView) = length.(grid.indices)
-
-# FIXME: this is terrible!
-Base.iterate(grid::GridView) = iterate(grid.grid[I] for I in CartesianIndices(grid.indices))
-Base.iterate(grid::GridView, state) = iterate([grid.grid[I] for I in CartesianIndices(grid.indices)], state)
-
-Base.getindex(grid::GridView, idx...) = getindex(grid.grid, (first.(grid.indices) .- 1 .+ Base.to_indices(grid, idx))...)
-
-# MetaGrid
-
-
-struct MetaGrid{d} <: AbstractGrid{d}
- grid::Grid{d}
- indices::Array{NTuple{d, UnitRange{Int}}, d}
- overlap::Int
-end
-
-"""
- MetaGrid(domain::NTuple{d, FDStepRange}, pnum::NTuple{d, Int}, overlap::Int)
-
-Create a grid on `domain`, partitioned along dimensions according to `pnum`
-respecting `overlap`.
-"""
-function MetaGrid(domain::NTuple{d, FDStepRange}, pnum::NTuple{d, Int}, overlap::Int) where d
- grid = Grid(domain)
- tsize = size(grid) .+ (pnum .- 1) .* overlap
-
- psize = tsize .÷ pnum
- osize = tsize .- pnum .* psize
-
- overhang(I, j) = I[j] == pnum[j] ? osize[j] : 0
-
- indices = Array{NTuple{d, UnitRange{Int}}, d}(undef, pnum)
- for I in CartesianIndices(pnum)
- indices[I] = ntuple(j -> ((I[j] - 1) * psize[j] - (I[j] - 1) * overlap + 1) :
- ( I[j] * psize[j] - (I[j] - 1) * overlap + overhang(I, j)), d)
- end
- MetaGrid{d}(grid, indices, overlap)
-end
-
-Base.size(grid::MetaGrid) = size(grid.grid)
-
-Base.iterate(grid::MetaGrid) = iterate(grid.grid)
-Base.iterate(grid::MetaGrid, state) = iterate(grid.grid, state)
-
-Base.getindex(grid::MetaGrid, idx...) = getindex(grid.grid, idx...)
-
-indices(grid::MetaGrid, idx...) = getindex(grid.indices, idx...)
-
-Base.view(grid::MetaGrid, idx...) = GridView(grid, indices(grid, idx...))
-
-parts(grid::MetaGrid) = CartesianIndices(grid.indices)
-
-
-## GridFunction
-
-abstract type AbstractGridFunction{G,T,d,n} <: AbstractArray{T,1} end
-
-rangedim(f::AbstractGridFunction{G,T,d,n}) where {G,T,d,n} = n
-domaindim(f::AbstractGridFunction{G,T,d}) where {G,T,d} = d
-
-Base.similar(f::F, ::Type{S}, dims::Int...) where {F<:AbstractGridFunction, S} =
- GridFunction(f.grid, Array{S}(undef, length(f)))
-
-
-"""
- GridFunction{G<:AbstractGrid,T,d,n,D}
-
-`n`-dimensional function values of type `T` on a Grid `G` of dimension `d`.
-
-`D == d + 1`
-"""
-struct GridFunction{G<:AbstractGrid,T,d,n,A,D} <: AbstractGridFunction{G,T,d,n}
- grid::G
- data::A
- function GridFunction{G,T,d,n,A,D}(grid::G, data::A) where {G,T,d,n,D,A<:AbstractArray{T,D}}
- D == d + 1 || throw(ArgumentError)
- d == ndims(grid) || throw(ArgumentError)
- (size(grid)..., n) == size(data) || throw(ArgumentError("data does not match given size"))
- new(grid, data)
- end
-end
-# TODO: check size(data)
-#GridFunction(grid::G, data::AbstractArray{T,D}) where {d,G<:AbstractGrid{d},T,D} =
-# GridFunction{G,T,d,size(data, d+1),typeof(data),d+1}(grid, data)
-function GridFunction(grid::G, data::AbstractArray, ::Val{n}) where {d,G<:AbstractGrid{d}, n}
- if length(data) ÷ length(grid) == 1 && n > 1
- data = repeat(data, inner=ntuple(i -> i == d+1 ? n : 1, d+1))
- else
- # TODO: check size(data)
- data = reshape(data, size(grid)..., n)
- end
- length(data) ÷ length(grid) == n || throw(ArgumentError("dimensions don't match"))
- GridFunction{G,eltype(data),d,n,typeof(data),d+1}(grid, data)
-end
-function GridFunction(grid::G, data::AbstractArray) where {d,G<:AbstractGrid{d}}
- GridFunction(grid, data, Val(length(data) ÷ length(grid)))
-end
-function GridFunction(f::Function, n::Int, grid::G) where {G<:AbstractGrid}
- data = Array{Float64,ndims(grid)+1}(undef, size(grid)..., n)
- for I in CartesianIndices(grid)
- data[Tuple(I)..., :] .= f(grid[I])
- end
- GridFunction(grid, data)
-end
-function GridFunction(grid::AbstractGrid, ::UndefInitializer, ::Val{n} = Val(1)) where n
- data = Array{Float64}(undef, size(grid)..., 1)
- GridFunction{typeof(grid), eltype(data), ndims(grid), n, typeof(data), ndims(grid) + 1}(grid, data)
-end
-
-#Base.show(io::IO, grid::GridFunction) =
-# print("GridFunction on $(length(grid)) grid points")
-
-Base.size(f::GridFunction) = (length(f.data),)
-#Base.axes(f::GridFunction) = (Base.OneTo(length(f.data)),)
-Base.getindex(f::GridFunction, idx::Int) = getindex(f.data, idx)
-Base.setindex!(f::GridFunction, v, idx::Int) = setindex!(f.data, v, idx)
-
-## TODO: using this prints #undef for display(f)
-#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()
-
-#component(f::GridFunction, k::Integer) = view(reshape(f.data,
-
-Base.view(f::AbstractGridFunction{G}, idx...) where G<:MetaGrid =
- GridFunction(GridView(f.grid, indices(f.grid, idx...)), view(f.data, indices(f.grid, idx...)..., :))
-
-grid(f::GridFunction) = f.grid
-
-
-Base.:*(a::Number, b::GridFunction) = GridFunction(b.grid, a * b.data, Val(rangedim(b)))
-
-
-
-import Colors: Gray
-import Plots: plot
-
-# TODO: add clamp option
-my_plot(img::AbstractGridFunction{G,T,d,1}) where {G,T,d} =
- plot(reshape(Gray.((img .- minimum(img)) ./ (maximum(img) - minimum(img))), size(img.grid)), show=true)
-
-## Operations
-
-function fdfilter(f::GridFunction{G,T,d,n}, kern) where {G,T,d,n}
- A = zeros(axes(f.data))
- for k = 1:n
- vidx = ntuple(i -> i <= d ? Colon() : k, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx...], kern)
- end
- GridFunction(f.grid, A)
-end
-
-derivate(f::GridFunction, dir::Integer; fd = FD.forward) = fdfilter(f, _fdkern(fd, ndims(f.grid), dir))
-
-
-# currently gradient and divergence are specific implementations (as needed for
-# chambolles algorithm), other choices may be reasonable
-
-function gradient(f::GridFunction{G,T,d,1}; fd=FD.forward) where {G,T,d}
- #A = Array{T,d+1}(undef, (size(f.grid)..., d))
- A = zeros(size(f.grid)..., d)
- for k = 1:d
- vidx = ntuple(i -> i <= d ? Colon() : k, d + 1)
- vidx2 = ntuple(i -> i <= d ? Colon() : 1, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx))
- # FIXME: hackish, only for forward fd
- selectdim(selectdim(A, d+1, k), k, size(f.grid, k)) .= 0
- end
- GridFunction(f.grid, A)
-end
-
-function divergence(f::GridFunction{G,T,d}; fd=FD.backward) where {d,G<:AbstractGrid{d},T}
- A = zeros(size(f.grid)..., 1)
- for k = 1:d
- vidx = ntuple(i -> i <= d ? Colon() : 1, d + 1)
- vidx2 = ntuple(i -> i <= d ? Colon() : k, d + 1)
- imfilter!(view(A, vidx...), f.data[vidx2...], bidx -> trim_kernel(_fdkern(fd, ndims(f.grid), k), bidx))
- end
- GridFunction(f.grid, A)
-end
-
-function pw_norm(f::GridFunction)
- d = domaindim(f)
- GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(1))
-end
-
-function my_norm(f::GridFunction)
- d = domaindim(f)
- GridFunction(f.grid, sqrt.(sum(f.data.^2, dims=d+1)), Val(rangedim(f)))
-end
-
-function pwmul(f::GridFunction{G,T,d,1}, g::GridFunction{G,T,d}) where {G,T,d}
- # TODO: maybe comprehension?
- data = similar(g.data)
- for I in CartesianIndices(f.grid)
- data[Tuple(I)..., :] = f.data[Tuple(I)..., 1] * g.data[Tuple(I)..., :]
- end
- GridFunction(f.grid, data, Val(rangedim(g)))
-end
-
-
-
-
-"""
- 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 * FD.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 = FD.dirfd(FD.backward, d, k)
- imfilter!(selectdim(dst, d+1, k), selectdim(field, d+1, k), bidx -> trim_kernel(kern, bidx))
- end
- out = dropdims(sum(dst, dims=d+1), dims=d+1)
- out[:]
-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 = FD.dirfd(FD.forward, d, k)
- imfilter!(selectdim(dst, d+1, k), img, bidx -> trim_kernel(kern, bidx))
- # FIXME: hack
- selectdim(selectdim(dst, d+1, k), k, size(grid, k)) .= 0
- end
- reshape(dst, length(grid) * ndims(grid))
-end
-
-function my_norm(grid::Grid, img::AbstractArray)
- img = reshape(img, (size(grid)..., ndims(grid)))
- d = ndims(grid)
-
- out = sqrt.(sum(img.^2, dims=d+1))
- out = repeat(out, inner=ntuple(i -> i == d+1 ? size(img, d+1) : 1, d+1))
-
- out[:]
-end
-
-
-function apply(kern, flatimg, size)
- img = reshape(flatimg, size)
- out = imfilter(kern, img)
- reshape(out, prod(size))
-end
-
-function solve(pde, shape::NTuple{d, Integer}) where d
- ndims(pde.stencil) == d || throw("dimension mismatch")
- length(pde.f) == prod(shape) || throw("bounds mismatch")
-
- A = LinearMap{Float64}(u -> apply(pde.stencil, u, shape), prod(shape), issymmetric=true, isposdef=true)
-
- u = zeros(prod(shape))
- cg!(u, A, pde.f)
- return u
-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 = FD.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 * FD.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 = FD.dirfd(FD.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 = FD.dirfd(FD.forward, d, k)
-# imfilter!(selectdim(dst, d+1, k), img, kern)
-# end
-# reshape(dst, length(grid) * ndims(grid))
-#end
-#
-#end
diff --git a/filtering.jl b/filtering.jl
deleted file mode 100644
index 11fa7b2..0000000
--- a/filtering.jl
+++ /dev/null
@@ -1,110 +0,0 @@
-module Filtering
-
-export imfilter, imfilter!, trim_kernel
-
-@inline function _shift_range(range, shift::Integer)
- Base.UnitRange(first(range)+shift : last(range)+shift)
-end
-
-## TODO: allow border to be tuple
-#function _innershift_indices(a::Array{T, d}, shift::NTuple{d, Integer}, border = abs(maximum(shift))) where {T, d}
-# abs(maximum(shift)) <= border || throw("cannot shift by $shift with border $border")
-# range(i) = (1 + shift[i] + border):(lastindex(a, i) + shift[i] - border)
-# ntuple(range, Val(ndims(a)))
-#end
-#
-#function innershift_view(a::Array{T, d}, shift::NTuple{d, Integer} = ntuple(i -> 0, d), border = abs(maximum(shift))) where {T, d}
-# view(a, _innershift_indices(a, shift, border)...)
-#end
-
-@inline function _bidx_to_range(range::AbstractUnitRange, bidx::Integer)
- if bidx > 0
- return last(range):last(range)
- elseif bidx < 0
- return first(range):first(range)
- else
- return first(range)+1:last(range)-1
- end
-end
-
-
-# we return Base.Slice since that is what `OffsetArray` uses for `axes()`
-function _trim_kernel_axes(kaxes, bidx)
- ntuple(Val(length(kaxes))) do i
- if bidx[i] < 0
- return Base.Slice(0:last(kaxes[i]))
- elseif bidx[i] > 0
- return Base.Slice(first(kaxes[i]):0)
- else
- return Base.Slice(kaxes[i])
- end
- end
-end
-
-"""
- trim_kernel(kernel, bidx)
-
-Trim `kernel` so that the resulting kernel is usable on the boundary
-indicated by `bidx`
-"""
-function trim_kernel(kernel, bidx)
- dst_axes = _trim_kernel_axes(axes(kernel), bidx)
- out = similar(kernel, dst_axes)
- range = CartesianIndices(out)
- copyto!(out, range, kernel, range)
-end
-
-
-"""
- imfilter!(dst, img, kern, slice=:inner)
-
-Filters `img` applying `kern` and write to `dst`. `dst` is added to, so you
-should initialize with zeros.
-"""
-# TODO: weirdly Base.Slice instead of Base.UnitRange leads to axes that are not
-# compatible when broadcasting, figure out why!
-function imfilter!(dst, img, kern, slice=_bidx_to_range.(axes(img), ntuple(i->0, ndims(img))))
- ndims(dst) == ndims(img) == ndims(kern) == length(slice) ||
- throw(ArgumentError("dst, image, kernel or slice dimensions disagree ($(ndims(dst)) vs $(ndims(img)) vs $(ndims(kern)) vs $(length(slice)))"))
- axes(dst) == axes(img) || throw(ArgumentError("output and input image axes disagree"))
- all(size(kern) .<= size(img)) || throw(ArgumentError("kernel bigger than image"))
- last.(slice) .+ last.(axes(kern)) <= last.(axes(img)) || first.(slice) .+ first.(axes(kern)) >= first.(axes(img)) ||
- throw(ArgumentError("kern and/or slice out of bounds for image"))
-
- dstslice = view(dst, UnitRange.(slice)...)
- for i in CartesianIndices(kern)
- iszero(kern[i]) && continue
- dstslice .+= kern[i] .* view(img, (_shift_range.(slice, Tuple(i)))...)
- end
- dst
-end
-
-function imfilter!(dst, img, kernelf::Function)
- for bidx in Iterators.product(ntuple(i -> -1:1, ndims(img))...)
- range = _bidx_to_range.(axes(img), bidx)
- imfilter!(dst, img, kernelf(bidx), range)
- end
- dst
-end
-
-"""
- imfilter(img, kernelf)
-
-Filter `img` using a kernel function `kernelf(bidx) = kernel` specifying a
-kernel for each type of node.
-"""
-# need zeros here, because `imfilter!` acts additively
-imfilter(img, kernelf::Function) = imfilter!(zeros(size(img)), img, kernelf)
-
-
-"""
- imfilter(img, kern)
-
-Filter `img` by `kern` using correlation.
-Boundary is zeroed out.
-"""
-imfilter(img, kern) = imfilter(img, bidx -> trim_kernel(kern, bidx))
-# TODO: use keyword arguments
-imfilter_inner(img, kern) = imfilter(img, bidx -> all(bidx .== 0) ? kern : zeros(ntuple(i -> 0, ndims(img))))
-
-end
diff --git a/src/DualTVDD.jl b/src/DualTVDD.jl
new file mode 100644
index 0000000..d1742f1
--- /dev/null
+++ b/src/DualTVDD.jl
@@ -0,0 +1,5 @@
+module DualTVDD
+
+greet() = print("Hello World!")
+
+end # module
diff --git a/test.py b/test.py
deleted file mode 100644
index 67db6cd..0000000
--- a/test.py
+++ /dev/null
@@ -1,134 +0,0 @@
-import numpy as np
-from scipy.sparse import diags, coo_matrix
-from scipy.sparse.linalg import spsolve, cg
-from mpl_toolkits.mplot3d import Axes3D
-import matplotlib.pyplot as plt
-
-# TODO: use np.indices, np.mgrid, np.meshgrid, etc.
-
-
-domain = np.array([100, 100, 100])
-d = domain.size
-N = np.prod(domain)
-
-
-# laplace stencil in d dimensions
-def laplace(d):
- st1 = np.array([1, -2, 1])
- s = st1.size
-
- st = np.zeros(np.full(d, s))
-
- idx = np.empty(d, dtype=slice)
- idx.fill(st1.size // 2)
- idx[0] = slice(None)
-
- for k in range(d):
- st[tuple(idx)] += st1
- idx = np.roll(idx, 1)
-
- return st
-
-stencil = laplace(d)
-
-print(stencil)
-
-
-#stencil = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]])
-
-
-u = np.zeros(N)
-f = 1 * np.ones(N)
-
-# boundary nodes (first and last planes are special)
-bnd = np.full(N, False)
-for nd in np.cumprod(np.flip(domain[1:], axis=0)):
- print(f"nd = {nd}")
- bnd = bnd | (np.arange(0, N) % nd == 0)
- bnd = bnd | (np.arange(0, N) % -nd == -1)
-bnd = bnd | (np.arange(0, N) // domain[-1] == 0)
-bnd = bnd | (np.arange(0, N) // domain[-1] == N // domain[-1] - 1)
-
-bnd_cnt = np.count_nonzero(bnd)
-free_cnt = np.prod(domain) - bnd_cnt
-
-
-# TODO: set u boundary data
-#u[bnd] = np.arange(0, N)[bnd]
-
-print(f"boundary nodes: {np.count_nonzero(bnd)}")
-print(f"free nodes: {np.count_nonzero(~bnd)}")
-
-elstrides = np.flip(np.cumprod(np.flip(np.append(domain[1:], [1]), 0)), 0)
-diagonals = []
-offsets = []
-
-
-center = (np.array(stencil.shape) - 1) // 2
-for index, value in np.ndenumerate(stencil):
- if value == 0:
- continue
- offset = np.array(index) - center
- offsets.append(np.dot(elstrides, offset))
- diagonals.append(value)
-
-print(diagonals, offsets)
-A = diags(diagonals, offsets, np.full(2, np.prod(domain))).tocoo()
-
-f = f - A @ u
-
-def filter_free(bnd, A):
- assert A.shape[0] == A.shape[1]
- print("filtering ...", end=" ", flush=True)
-
- # filter out boundary nodes
- spfree = ~bnd[A.row] & ~bnd[A.col]
- rows, cols, data = A.row[spfree], A.col[spfree], A.data[spfree]
-
- # correct indices
- offset = -np.cumsum(bnd)
- rows = rows + offset[rows]
- cols = cols + offset[cols]
-
- nf = A.shape[0] - np.count_nonzero(bnd)
- Af = coo_matrix((data, (rows, cols)), shape=(nf, nf))
- print("finished!")
- return Af
-
-Af = filter_free(bnd, A).tocsr()
-
-#Af = A[np.ix_(~bnd, ~bnd)]
-ff = f[~bnd]
-
-print("solving ...", end=" ", flush=True)
-u[~bnd] = spsolve(Af, ff)
-#u[~bnd] = cg(Af, ff)
-print("finished!")
-
-
-x, y = np.meshgrid(np.arange(0, domain[1]), np.arange(0, domain[0]))
-
-print(x.shape, y.shape)
-z = u.reshape(domain)
-
-fig = plt.figure()
-ax = fig.add_subplot(111, projection='3d')
-ax.plot_surface(x, y, z)
-plt.show()
-#
-#
-#print(u.reshape(domain))
-#
-#print(np.diagflat(np.full((n), value)))
-
-
-#np.roll(p, -1) + 2*p + np.roll(
-
-
-
-
-#s = np.sin(2 * np.pi * t)
-#
-#plt.plot(t, s)
-#plt.show()
-
diff --git a/DualTVDD/repl.jl b/test/repl.jl
similarity index 98%
rename from DualTVDD/repl.jl
rename to test/repl.jl
index aa11454..bbeea48 100644
--- a/DualTVDD/repl.jl
+++ b/test/repl.jl
@@ -1,7 +1,6 @@
-using Pkg
-Pkg.activate(".")
-
using BenchmarkTools
using Revise
+using Pkg
+Pkg.activate(".")
using DualTVDD
--
GitLab