diff --git a/Manifest.toml b/Manifest.toml
index 33160c4a7877e565f95a7e112916535e137cddf4..09c2c24aa3b57db625682c764b7f1c9298151631 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,53 +1,11 @@
 # This file is machine-generated - editing it directly is not advised
 
-[[AbstractFFTs]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716"
-uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "0.5.0"
-
-[[AbstractNumbers]]
-deps = ["SpecialFunctions"]
-git-tree-sha1 = "924732f1a32f0e4b1e748b469d4604125bdf1d00"
-uuid = "85c772de-338a-5e7f-b815-41e76c26ac1f"
-version = "0.2.1"
-
-[[AbstractPlotting]]
-deps = ["AbstractNumbers", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DocStringExtensions", "FFMPEG", "FileIO", "FixedPointNumbers", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "ImageMagick", "IntervalSets", "LinearAlgebra", "Markdown", "Observables", "Packing", "PlotUtils", "Printf", "Random", "Serialization", "Showoff", "SignedDistanceFields", "SparseArrays", "StaticArrays", "Statistics", "StructArrays", "UnicodeFun"]
-git-tree-sha1 = "cd373780f00edbd6429ceacb2572b55160c01886"
-uuid = "537997a7-5e4e-5d89-9595-2241ea00577e"
-version = "0.10.11"
-
 [[Adapt]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "c88cfc7f9c1f9f8633cddf0b56e86302b70f64c5"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 version = "1.0.1"
 
-[[Arpack]]
-deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
-git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
-uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
-version = "0.4.0"
-
-[[Arpack_jll]]
-deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
-git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
-uuid = "68821587-b530-5797-8361-c406ea357684"
-version = "3.5.0+3"
-
-[[ArrayInterface]]
-deps = ["LinearAlgebra", "Requires", "SparseArrays"]
-git-tree-sha1 = "649c08a5a3a513f4662673d3777fe6ccb4df9f5d"
-uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "2.8.7"
-
-[[AxisAlgorithms]]
-deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
-git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
-uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
-version = "1.0.0"
-
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
@@ -57,17 +15,6 @@ git-tree-sha1 = "3663bfffede2ef41358b6fc2e1d8a6d50b3c3904"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
 version = "1.0.6+2"
 
-[[CEnum]]
-git-tree-sha1 = "62847acab40e6855a9b5905ccb99c2b5cf6b3ebb"
-uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
-version = "0.2.0"
-
-[[ColorBrewer]]
-deps = ["Colors", "JSON", "Test"]
-git-tree-sha1 = "61c5334f33d91e570e1d0c3eb5465835242582c4"
-uuid = "a2cac450-b92f-5266-8821-25eda20663c8"
-version = "0.4.0"
-
 [[ColorSchemes]]
 deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
 git-tree-sha1 = "7a15e3690529fd1042f0ab954dff7445b1efc8a5"
@@ -80,36 +27,12 @@ git-tree-sha1 = "c73d9cfc2a9d8433dc77f5bff4bddf46b1d78c20"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 version = "0.10.3"
 
-[[ColorVectorSpace]]
-deps = ["ColorTypes", "Colors", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "StatsBase"]
-git-tree-sha1 = "bd0c0c81a39923bc03f9c3b61d89ad816e741002"
-uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
-version = "0.8.5"
-
 [[Colors]]
 deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"]
 git-tree-sha1 = "db9c215e2d92379d341e40dd1d5be893af2bfeb2"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
 version = "0.12.1"
 
-[[CommonSubexpressions]]
-deps = ["Test"]
-git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0"
-uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
-version = "0.2.0"
-
-[[Compat]]
-deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "48c20c43e157c6eab6cf88326504ec042b05e456"
-uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "3.10.0"
-
-[[CompilerSupportLibraries_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612"
-uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "0.3.3+0"
-
 [[Contour]]
 deps = ["StaticArrays"]
 git-tree-sha1 = "0b17db36e7e03f8437e0d1f55aea3e4a60c74353"
@@ -140,45 +63,10 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 deps = ["Mmap"]
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
-[[DiffResults]]
-deps = ["StaticArrays"]
-git-tree-sha1 = "da24935df8e0c6cf28de340b958f6aac88eaa0cc"
-uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
-version = "1.0.2"
-
-[[DiffRules]]
-deps = ["NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "eb0c34204c8410888844ada5359ac8b96292cfd1"
-uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.0.1"
-
-[[Distances]]
-deps = ["LinearAlgebra", "Statistics"]
-git-tree-sha1 = "23717536c81b63e250f682b0e0933769eecd1411"
-uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
-version = "0.8.2"
-
 [[Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
-[[Distributions]]
-deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "c4ed10355637fcb0725dc6a27060f74df24f13cd"
-uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.23.2"
-
-[[DocStringExtensions]]
-deps = ["LibGit2", "Markdown", "Pkg", "Test"]
-git-tree-sha1 = "88bb0edb352b16608036faadcc071adda068582a"
-uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.1"
-
-[[EllipsisNotation]]
-git-tree-sha1 = "65dad386e877850e6fce4fc77f60fe75a468ce9d"
-uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
-version = "0.4.0"
-
 [[FFMPEG]]
 deps = ["FFMPEG_jll"]
 git-tree-sha1 = "c82bef6fc01e30d500f588cd01d29bdd44f1924e"
@@ -191,94 +79,28 @@ git-tree-sha1 = "0fa07f43e5609ea54848b82b4bb330b250e9645b"
 uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
 version = "4.1.0+3"
 
-[[FFTW]]
-deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
-git-tree-sha1 = "b6a74f6dfd9e9d16b765397dc90df03e5a00532e"
-uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.2.1"
-
-[[FFTW_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b"
-uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
-version = "3.3.9+5"
-
-[[FileIO]]
-deps = ["Pkg"]
-git-tree-sha1 = "202335fd24c2776493e198d6c66a6d910400a895"
-uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.3.0"
-
-[[FillArrays]]
-deps = ["LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "44f561e293987ffc84272cd3d2b14b0b93123d63"
-uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.8.10"
-
-[[FiniteDiff]]
-deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
-git-tree-sha1 = "fec7c2cb45c27071ef487fa7cae4fcac7509aa10"
-uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.3.2"
-
 [[FixedPointNumbers]]
 git-tree-sha1 = "3ba9ea634d4c8b289d590403b4a06f8e227a6238"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 version = "0.8.0"
 
-[[ForwardDiff]]
-deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "869540e4367122fbffaace383a5bdc34d6e5e5ac"
-uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.10"
-
-[[FreeType]]
-deps = ["CEnum", "FreeType2_jll"]
-git-tree-sha1 = "e2986a14a511f9d289b69bbdf09691fea5203976"
-uuid = "b38be410-82b0-50bf-ab77-7b57e271db43"
-version = "3.0.0"
-
 [[FreeType2_jll]]
 deps = ["Bzip2_jll", "Libdl", "Pkg", "Zlib_jll"]
 git-tree-sha1 = "7d900f32a3788d4eacac2bfa3bf5c770179c8afd"
 uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
 version = "2.10.1+2"
 
-[[FreeTypeAbstraction]]
-deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics", "StaticArrays"]
-git-tree-sha1 = "2bd5787060aa9883ec258523b728a2adf5ca9cba"
-uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
-version = "0.8.0"
-
 [[FriBidi_jll]]
 deps = ["Libdl", "Pkg"]
 git-tree-sha1 = "2f56bee16bd0151de7b6a1eeea2ced190a2ad8d4"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
 version = "1.0.5+3"
 
-[[GLFW]]
-deps = ["GLFW_jll"]
-git-tree-sha1 = "7c47cfbfe9e1ff1cf658ff3e739eb1b56ee95cbe"
-uuid = "f7f18e0c-5ee9-5ccd-a5bf-e8befd85ed98"
-version = "3.2.2"
-
-[[GLFW_jll]]
-deps = ["Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
-git-tree-sha1 = "3b04f1dd73f0efa81298a76c9515f306df06f613"
-uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
-version = "3.3.2+0"
-
-[[GLM]]
-deps = ["Distributions", "LinearAlgebra", "Printf", "Random", "Reexport", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "StatsModels"]
-git-tree-sha1 = "db0ace36f9dbe7b6a7a08434c5921377e9df2c72"
-uuid = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
-version = "1.3.9"
-
-[[GLMakie]]
-deps = ["AbstractPlotting", "ColorTypes", "Colors", "FileIO", "FixedPointNumbers", "FreeTypeAbstraction", "GLFW", "GeometryBasics", "LinearAlgebra", "Markdown", "MeshIO", "ModernGL", "Observables", "Printf", "Serialization", "ShaderAbstractions", "StaticArrays"]
-git-tree-sha1 = "2ddb89b95c3061f4e9bcf55b7e71454fed5eadfc"
-uuid = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-version = "0.1.4"
+[[GR]]
+deps = ["Base64", "DelimitedFiles", "HTTP", "JSON", "LinearAlgebra", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
+git-tree-sha1 = "247adbd2b33c0c4b42efa20d1e807acf6312145f"
+uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
+version = "0.50.1"
 
 [[GeometryBasics]]
 deps = ["IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
@@ -286,52 +108,28 @@ git-tree-sha1 = "4ab7ad8996fc164a80cf18a7ef3a0b9e3465ae8e"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 version = "0.2.9"
 
-[[Graphics]]
-deps = ["Colors", "LinearAlgebra", "NaNMath"]
-git-tree-sha1 = "45d684ead5b65c043ad46bd5be750d61c39d7ef8"
-uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
-version = "1.0.2"
+[[GeometryTypes]]
+deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "StaticArrays"]
+git-tree-sha1 = "34bfa994967e893ab2f17b864eec221b3521ba4d"
+uuid = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+version = "0.8.3"
 
-[[ImageCore]]
-deps = ["Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport", "Requires"]
-git-tree-sha1 = "a652c05f8f374861580d420b420fddf3e2e84312"
-uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
-version = "0.8.14"
-
-[[ImageMagick]]
-deps = ["FileIO", "ImageCore", "ImageMagick_jll", "InteractiveUtils", "Libdl", "Pkg", "Random"]
-git-tree-sha1 = "318342a5099a9c952b4de087344a2db831458c87"
-uuid = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
-version = "1.1.5"
-
-[[ImageMagick_jll]]
-deps = ["JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "1c0a2295cca535fabaf2029062912591e9b61987"
-uuid = "c73af94c-d91f-53ed-93a7-00f77d67a9d7"
-version = "6.9.10-12+3"
-
-[[IntelOpenMP_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6"
-uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
-version = "2018.0.3+0"
+[[HTTP]]
+deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
+git-tree-sha1 = "eca61b35cdd8cd2fcc5eec1eda766424a995b02f"
+uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
+version = "0.8.16"
+
+[[IniFile]]
+deps = ["Test"]
+git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
+uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
+version = "0.5.0"
 
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
-[[Interpolations]]
-deps = ["AxisAlgorithms", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
-git-tree-sha1 = "3af735234d9b1ff9ff1af89875735cd9549c0c5f"
-uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-version = "0.12.9"
-
-[[IntervalSets]]
-deps = ["Dates", "EllipsisNotation", "Statistics"]
-git-tree-sha1 = "3b1cef135bc532b3c3401b309e1b8a2a2ba26af5"
-uuid = "8197267c-284f-5f27-9208-e0e47529a953"
-version = "0.5.1"
-
 [[IterTools]]
 git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
 uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
@@ -348,18 +146,6 @@ git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 version = "0.21.0"
 
-[[JpegTurbo_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "b007cedc10bb017cdae3af0062fe4dd29b483c70"
-uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "2.0.1+0"
-
-[[KernelDensity]]
-deps = ["Distributions", "FFTW", "Interpolations", "Optim", "StatsBase", "Test"]
-git-tree-sha1 = "c1048817fe5711f699abc8fabd47b1ac6ba4db04"
-uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
-version = "0.5.1"
-
 [[LAME_jll]]
 deps = ["Libdl", "Pkg"]
 git-tree-sha1 = "221cc8998b9060677448cbb6375f00032554c4fd"
@@ -379,81 +165,33 @@ version = "1.8.1+1"
 [[Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
-[[Libgcrypt_jll]]
-deps = ["Libdl", "Libgpg_error_jll", "Pkg"]
-git-tree-sha1 = "7efb7f5e75a0f7323cce366d3e44a2d225c5617a"
-uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
-version = "1.8.5+1"
-
-[[Libglvnd_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
-git-tree-sha1 = "30a990af9ce535bcba54d1f23f071a7e714f861e"
-uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
-version = "1.3.0+1"
-
-[[Libgpg_error_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "334284a62d6b8592099d28b78c76aefc266c03a2"
-uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
-version = "1.36.0+1"
-
-[[Libiconv_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "e5256a3b0ebc710dbd6da0c0b212164a3681037f"
-uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
-version = "1.16.0+2"
-
-[[Libtiff_jll]]
-deps = ["JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "1fe8c3608dfe7bdec81d018de1cc66e959016e8c"
-uuid = "89763e89-9b03-5906-acba-b20f662cd828"
-version = "4.1.0+0"
-
-[[LineSearches]]
-deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf", "Test"]
-git-tree-sha1 = "54eb90e8dbe745d617c78dee1d6ae95c7f6f5779"
-uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.0.1"
-
 [[LinearAlgebra]]
 deps = ["Libdl"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
-[[Loess]]
-deps = ["Distances", "Statistics"]
-git-tree-sha1 = "67b1bc5f42144fb83b065b8a8e117e56ec10e601"
-uuid = "4345ca2d-374a-55d4-8d30-97f9976e7612"
-version = "0.5.1"
-
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
-[[MKL_jll]]
-deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "720629cc8cbd12c146ca01b661fd1a6cf66e2ff4"
-uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2019.0.117+2"
-
-[[Makie]]
-deps = ["AbstractPlotting", "FileIO", "GLMakie", "ImageMagick", "MeshIO", "StatsMakie"]
-git-tree-sha1 = "3c77041a678e5728f96eacfb734ebbe5c5231a1e"
-uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-version = "0.10.0"
-
-[[MappedArrays]]
-git-tree-sha1 = "e2a02fe7ee86a10c707ff1756ab1650b40b140bb"
-uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
-version = "0.2.2"
-
 [[Markdown]]
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
-[[MeshIO]]
-deps = ["ColorTypes", "FileIO", "GeometryBasics", "Printf"]
-git-tree-sha1 = "2699b809508235eb84069e37aef9c0ec1b0001a4"
-uuid = "7269a6da-0436-5bbc-96c2-40638cbb6118"
-version = "0.4.0"
+[[MbedTLS]]
+deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
+git-tree-sha1 = "426a6978b03a97ceb7ead77775a1da066343ec6e"
+uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
+version = "1.0.2"
+
+[[MbedTLS_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "f85473aeb7a2561a5c58c06c4868971ebe2bcbff"
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.16.6+0"
+
+[[Measures]]
+git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
+uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
+version = "0.3.1"
 
 [[Missings]]
 deps = ["DataAPI"]
@@ -464,68 +202,22 @@ version = "0.4.3"
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
-[[ModernGL]]
-deps = ["Libdl"]
-git-tree-sha1 = "326957556e9cc9253615114c04bb0096a2a69bb8"
-uuid = "66fc600b-dfda-50eb-8b99-91cfa97b1301"
-version = "1.1.2"
-
-[[MosaicViews]]
-deps = ["OffsetArrays", "PaddedViews"]
-git-tree-sha1 = "b483b88403ac0ac01667778cbb29462b111b1deb"
-uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
-version = "0.2.2"
-
-[[NLSolversBase]]
-deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
-git-tree-sha1 = "7c4e66c47848562003250f28b579c584e55becc0"
-uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
-version = "7.6.1"
-
 [[NaNMath]]
 git-tree-sha1 = "928b8ca9b2791081dc71a51c55347c27c618760f"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 version = "0.3.3"
 
-[[Observables]]
-git-tree-sha1 = "11832878355305984235a2e90d0e3737383c634c"
-uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
-version = "0.3.1"
-
-[[OffsetArrays]]
-git-tree-sha1 = "930db8ef90483570107f2396b1ffc6680f08e8b7"
-uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.0.4"
-
 [[Ogg_jll]]
 deps = ["Libdl", "Pkg"]
 git-tree-sha1 = "59cf7a95bf5ac39feac80b796e0f39f9d69dc887"
 uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
 version = "1.3.4+0"
 
-[[OpenBLAS_jll]]
-deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "1887096f6897306a4662f7c5af936da7d5d1a062"
-uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.9+4"
-
 [[OpenSSL_jll]]
 deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "d2a6f25262d568b5a7e454cf7ff5066a79d16c7d"
+git-tree-sha1 = "7aaaded15bf393b5f34c2aad5b765c18d26cb495"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "1.1.1+2"
-
-[[OpenSpecFun_jll]]
-deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "d51c416559217d974a1113522d5919235ae67a87"
-uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.3+3"
-
-[[Optim]]
-deps = ["Compat", "FillArrays", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
-git-tree-sha1 = "33af70b64e8ce2f2b857e3d5de7b71f67715c121"
-uuid = "429524aa-4258-5aef-a3af-852621145aeb"
-version = "0.21.0"
+version = "1.1.1+4"
 
 [[Opus_jll]]
 deps = ["Libdl", "Pkg"]
@@ -538,30 +230,6 @@ git-tree-sha1 = "12ce190210d278e12644bcadf5b21cbdcf225cd3"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.2.0"
 
-[[PDMats]]
-deps = ["Arpack", "LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"]
-git-tree-sha1 = "2fc6f50ddd959e462f0a2dbc802ddf2a539c6e35"
-uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.9.12"
-
-[[Packing]]
-deps = ["GeometryBasics"]
-git-tree-sha1 = "460dfc95d8112aca861b139951961d6deb8ebdb1"
-uuid = "19eb6ba3-879d-56ad-ad62-d5c202156566"
-version = "0.4.0"
-
-[[PaddedViews]]
-deps = ["OffsetArrays"]
-git-tree-sha1 = "100195a79b577d5747db98bf1732c3686285fa1e"
-uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
-version = "0.5.5"
-
-[[Parameters]]
-deps = ["OrderedCollections", "UnPack"]
-git-tree-sha1 = "38b2e970043613c187bd56a995fe2e551821eb4a"
-uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-version = "0.12.1"
-
 [[Parsers]]
 deps = ["Dates", "Test"]
 git-tree-sha1 = "f0abb338b4d00306500056a3fd44c221b8473ef2"
@@ -572,28 +240,28 @@ version = "1.0.4"
 deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
+[[PlotThemes]]
+deps = ["PlotUtils", "Requires", "Statistics"]
+git-tree-sha1 = "c6f5ea535551b3b16835134697f0c65d06c94b91"
+uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
+version = "2.0.0"
+
 [[PlotUtils]]
 deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
 git-tree-sha1 = "59ec24a0c96c513533e488dff1433df1bd3d6b9f"
 uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
 version = "1.0.3"
 
-[[PositiveFactorizations]]
-deps = ["LinearAlgebra", "Test"]
-git-tree-sha1 = "127c47b91990c101ee3752291c4f45640eeb03d1"
-uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
-version = "0.2.3"
+[[Plots]]
+deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "GeometryTypes", "JSON", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
+git-tree-sha1 = "d2790295d9b65b98c44fe74815743bc101df2e19"
+uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+version = "1.5.2"
 
 [[Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
-[[QuadGK]]
-deps = ["DataStructures", "LinearAlgebra"]
-git-tree-sha1 = "dc84e810393cfc6294248c9032a9cdacc14a3db4"
-uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-version = "2.3.1"
-
 [[REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -602,10 +270,16 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
-[[Ratios]]
-git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41"
-uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
-version = "0.4.0"
+[[RecipesBase]]
+git-tree-sha1 = "54f8ceb165a0f6d083f0d12cb4996f5367c6edbc"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "1.0.1"
+
+[[RecipesPipeline]]
+deps = ["Dates", "PlotUtils", "RecipesBase"]
+git-tree-sha1 = "d2a58b8291d1c0abae6a91489973f8a92bf5c04a"
+uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
+version = "0.1.11"
 
 [[Reexport]]
 deps = ["Pkg"]
@@ -619,51 +293,18 @@ git-tree-sha1 = "d37400976e98018ee840e0ca4f9d20baa231dc6b"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.0.1"
 
-[[Rmath]]
-deps = ["Random", "Rmath_jll"]
-git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370"
-uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
-version = "0.6.1"
-
-[[Rmath_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "1660f8fefbf5ab9c67560513131d4e933012fc4b"
-uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
-version = "0.2.2+0"
-
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
-[[ShaderAbstractions]]
-deps = ["ColorTypes", "FixedPointNumbers", "GeometryBasics", "LinearAlgebra", "Observables", "StaticArrays", "StructArrays", "Tables"]
-git-tree-sha1 = "0385d66d733b2588d9f99359da7e952829081f4e"
-uuid = "65257c39-d410-5151-9873-9b3e5be5013e"
-version = "0.2.1"
-
-[[SharedArrays]]
-deps = ["Distributed", "Mmap", "Random", "Serialization"]
-uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
-
-[[ShiftedArrays]]
-git-tree-sha1 = "22395afdcf37d6709a5a0766cc4a5ca52cb85ea0"
-uuid = "1277b4bf-5013-50f5-be3d-901d8477a67a"
-version = "1.0.0"
-
 [[Showoff]]
 deps = ["Dates"]
 git-tree-sha1 = "e032c9df551fb23c9f98ae1064de074111b7bc39"
 uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
 version = "0.3.1"
 
-[[SignedDistanceFields]]
-deps = ["Random", "Statistics", "Test"]
-git-tree-sha1 = "d263a08ec505853a5ff1c1ebde2070419e3f28e9"
-uuid = "73760f76-fbc4-59ce-8f25-708e95d2df96"
-version = "0.4.0"
-
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
@@ -677,12 +318,6 @@ version = "0.3.1"
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-[[SpecialFunctions]]
-deps = ["OpenSpecFun_jll"]
-git-tree-sha1 = "d8d8b8a9f4119829410ecd706da4cc8594a1e020"
-uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "0.10.3"
-
 [[StaticArrays]]
 deps = ["LinearAlgebra", "Random", "Statistics"]
 git-tree-sha1 = "5c06c0aeb81bef54aed4b3f446847905eb6cbda0"
@@ -704,34 +339,12 @@ git-tree-sha1 = "a6102b1f364befdb05746f386b67c6b7e3262c45"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 version = "0.33.0"
 
-[[StatsFuns]]
-deps = ["Rmath", "SpecialFunctions"]
-git-tree-sha1 = "04a5a8e6ab87966b43f247920eab053fd5fdc925"
-uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "0.9.5"
-
-[[StatsMakie]]
-deps = ["AbstractPlotting", "Distributions", "GLM", "IntervalSets", "KernelDensity", "Loess", "Observables", "Random", "Statistics", "StatsBase", "StructArrays", "Tables"]
-git-tree-sha1 = "64dc4424c8cb4218938d8e5300e9c1c4d37d1c51"
-uuid = "65254759-4cff-5aa5-8326-61ce017a8c70"
-version = "0.2.1"
-
-[[StatsModels]]
-deps = ["DataAPI", "DataStructures", "LinearAlgebra", "ShiftedArrays", "SparseArrays", "StatsBase", "Tables"]
-git-tree-sha1 = "5a140e0bf5d92b2c5c0daa0dcefdadbb57c8a7fb"
-uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
-version = "0.6.11"
-
 [[StructArrays]]
 deps = ["Adapt", "DataAPI", "Tables"]
 git-tree-sha1 = "10ee2e9b8a222ef7aac886b12cc8c050db9a8a45"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 version = "0.4.3"
 
-[[SuiteSparse]]
-deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
-uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
-
 [[TableTraits]]
 deps = ["IteratorInterfaceExtensions"]
 git-tree-sha1 = "b1ad568ba658d8cbb3b892ed5380a6f3e781a81e"
@@ -752,127 +365,14 @@ uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
-[[UnPack]]
-git-tree-sha1 = "d4bfa022cd30df012700cf380af2141961bb3bfb"
-uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
-version = "1.0.1"
-
 [[Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
-[[UnicodeFun]]
-deps = ["REPL"]
-git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
-uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
-version = "0.4.1"
-
-[[WoodburyMatrices]]
-deps = ["LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "28ffe06d28b1ba8fdb2f36ec7bb079fac81bac0d"
-uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
-version = "0.5.2"
-
-[[XML2_jll]]
-deps = ["Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "987c02a43fa10a491a5f0f7c46a6d3559ed6a8e2"
-uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.9.9+4"
-
-[[XSLT_jll]]
-deps = ["Libdl", "Libgcrypt_jll", "Pkg", "XML2_jll"]
-git-tree-sha1 = "c287bbc08773474ab1d12849eb46b961a325e149"
-uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
-version = "1.1.33+2"
-
-[[Xorg_libX11_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
-git-tree-sha1 = "bbfb646a5a86dc636041c8ee945799b69d37fdba"
-uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
-version = "1.6.9+2"
-
-[[Xorg_libXau_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "711f2bbdbb0e90b90a71c5ccf835137e2def539a"
-uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
-version = "1.0.9+2"
-
-[[Xorg_libXcursor_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
-git-tree-sha1 = "fd1f86fcb0295a6741409086723d9e7d8f222599"
-uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
-version = "1.2.0+2"
-
-[[Xorg_libXdmcp_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "1ffd7bc9c9d8bba6652516f3402ee1ef5f4cac70"
-uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
-version = "1.1.3+2"
-
-[[Xorg_libXext_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "079acccb5437761d93c483f2edb5fc3ed21ea7db"
-uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
-version = "1.3.4+2"
-
-[[Xorg_libXfixes_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "2bc6823df6ccd55927bf8441d3d79db2f5874c2a"
-uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
-version = "5.0.3+2"
-
-[[Xorg_libXi_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
-git-tree-sha1 = "dedd853f81c1099ada169fb6f9c2b50764c48f06"
-uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
-version = "1.7.10+2"
-
-[[Xorg_libXinerama_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libXext_jll"]
-git-tree-sha1 = "79f0b3350e62d02a710a7dacfbbc0878aed09d58"
-uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
-version = "1.1.4+2"
-
-[[Xorg_libXrandr_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
-git-tree-sha1 = "206e6ab109a2b607e3076eb2cb9b3fa239c467f1"
-uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
-version = "1.5.2+2"
-
-[[Xorg_libXrender_jll]]
-deps = ["Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "d5492a36d33359cf93ad573ea5e235cd71728649"
-uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
-version = "0.9.10+2"
-
-[[Xorg_libpthread_stubs_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "ec1e2973e4e624069be95f806dbd7e225fee26c6"
-uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
-version = "0.1.0+1"
-
-[[Xorg_libxcb_jll]]
-deps = ["Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
-git-tree-sha1 = "c1efb50fb6a8c40c45451611fbf93feabe5f9611"
-uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
-version = "1.13.0+2"
-
-[[Xorg_xtrans_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "b6e7d70f5fc9e3fde505539194c8caaab3c40ebf"
-uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
-version = "1.4.0+1"
-
 [[Zlib_jll]]
 deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "a2e0d558f6031002e380a90613b199e37a8565bf"
+git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71"
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.11+10"
-
-[[Zstd_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "b25b0fb10176c42e9a5a20e1f40d570ac0288d4e"
-uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
-version = "1.4.5+0"
+version = "1.2.11+14"
 
 [[libass_jll]]
 deps = ["Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "Libdl", "Pkg", "Zlib_jll"]
@@ -886,12 +386,6 @@ git-tree-sha1 = "480c7ed04f68ea3edd4c757f5db5b6a0a4e0bd99"
 uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
 version = "0.1.6+2"
 
-[[libpng_jll]]
-deps = ["Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "594cb058723c13941cf463fd09e5859499594f50"
-uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.37+3"
-
 [[libvorbis_jll]]
 deps = ["Libdl", "Ogg_jll", "Pkg"]
 git-tree-sha1 = "6a66f65b5275dfa799036c8a3a26616a0a271c4a"
diff --git a/Project.toml b/Project.toml
index 367f6423388d58c4fd7b01854cb562a5000cd466..2f5f8aefa89876c0f45e40bbc16063aa0aeba24d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,7 @@ version = "0.1.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StaticKernels = "4c63dfa8-a427-4548-bd2f-4c19e87a7dc7"
 
diff --git a/src/DualTVDD.jl b/src/DualTVDD.jl
index cdbf66d53115a9a0ef9602dd79503f38276af774..af27db1891d5deb99efbb610f5be9bc69e3e05d3 100644
--- a/src/DualTVDD.jl
+++ b/src/DualTVDD.jl
@@ -1,15 +1,12 @@
 module DualTVDD
 
-
-
 include("common.jl")
 include("problems.jl")
 include("chambolle.jl")
 include("dualtvdd.jl")
 include("projgrad.jl")
 
-
-using Makie: heatmap
+using Plots
 
 function run()
     #g = [0. 2; 1 0.]
@@ -25,9 +22,9 @@ function run()
 
     display(norm(B))
 
-    md = DualTVDD.DualTVL1ROFOpProblem(g, B, λ)
-    ctx = DualTVDD.init(md, DualTVDD.ChambolleAlgorithm())
-    ctx2 = DualTVDD.init(md, DualTVDD.ProjGradAlgorithm(τ=1/sqrt(8)/norm(B)))
+    prob = DualTVDD.DualTVL1ROFOpProblem(g, B, λ)
+    ctx = DualTVDD.init(DualTVDD.ChambolleAlgorithm(prob))
+    ctx2 = DualTVDD.init(DualTVDD.ProjGradAlgorithm(prob, τ=1/sqrt(8)/norm(B)))
 
     #scene = vbox(
     #    heatmap(ctx.s, colorrange=(0,1), colormap=:gray, scale_plot=false, show_axis=false),
@@ -44,13 +41,13 @@ function run()
         #sleep(0.2)
     end
     display(copy(ctx.p))
-    display(copy(recover_u!(ctx)))
+    #display(copy(recover_u!(ctx)))
     println(energy(ctx))
 
     println()
 
     display(copy(ctx2.p))
-    display(copy(recover_u!(ctx2)))
+    #display(copy(recover_u!(ctx2)))
     println(energy(ctx2))
 
     ctx, ctx2
@@ -58,43 +55,89 @@ function run()
     #yield()
 end
 
+function alg_energy(alg, niter)
+    print("run ...")
+    res = Float64[]
+    (p, ctx) = iterate(alg)
+    push!(res, energy(ctx))
+    for i in 1:niter
+        (p, ctx) = iterate(alg, ctx)
+        push!(res, energy(ctx))
+    end
+    println(" finished")
+    return res
+end
+
+function alg_error(alg, pmin, niter, ninner=1)
+    print("run ...")
+    res = Float64[]
+    (p, ctx) = iterate(alg)
+    push!(res, error(ctx.p, pmin, ctx.algorithm.problem))
+    for i in 1:niter
+        for j in 1:ninner
+            (p, ctx) = iterate(alg, ctx)
+        end
+        push!(res, error(ctx.p, pmin, ctx.algorithm.problem))
+    end
+    println(" finished")
+    return res
+end
+
+function calc_energy(prob, niter)
+    alg_ref = DualTVDD.ChambolleAlgorithm(prob)
+    ctx = init(alg_ref)
+    for i in 1:niter
+        (p, ctx) = iterate(alg_ref, ctx)
+    end
+    return energy(ctx), p
+end
+
 function rundd()
-    λ = 1/2
+    λ = 2.5
     β = 0
-    f = zeros(10,10)
-    f[1] = 1
+    #f = zeros(100,100)
+    #f[1] = 1
     #f = [0. 2; 1 0.]
 
 
-    A = 0. * rand(length(f), length(f))
-    A .+= diagm(ones(length(f)))
-    B = inv(A'*A + β*I)
+    #A = 0. * rand(length(f), length(f))
+    #A .+= diagm(ones(length(f)))
+    #B = inv(A'*A + β*I)
 
-    g = similar(f)
-    vec(g) .= A' * vec(f)
+    #g = similar(f)
+    #vec(g) .= A' * vec(f)
+
+    g = rand(50, 50)
+
+    prob = DualTVDD.DualTVL1ROFOpProblem(g, I, λ)
 
-    prob = DualTVDD.DualTVL1ROFOpProblem(g, B, λ)
-    alg_dd = DualTVDD.DualTVDDAlgorithm(prob, M=(2,2), overlap=(2,2), σ=0.25)
     alg_ref = DualTVDD.ChambolleAlgorithm(prob)
+    alg_dd = DualTVDD.DualTVDDAlgorithm(prob, M=(2,2), overlap=(4,4))
+    alg_dd2 = DualTVDD.DualTVDDAlgorithm(prob, M=(2,2), overlap=(4,4), parallel=false)
 
-    (p, ctx) = iterate(alg_ref)
-    for i in 1:100000
-        (p, ctx) = iterate(alg_ref, ctx)
-    end
-    display(ctx.p)
-    display(recover_u(p, prob))
+    n = 1000
 
+    ref_energy, pmin = calc_energy(prob, 100000)
 
-    (p, ctx) = iterate(alg_dd)
-    for i in 1:1000
-        (p, ctx) = iterate(alg_dd, ctx)
-    end
-    display(ctx.p)
-    #display(ctx.subctx[1].algorithm.problem.g)
-    #display(ctx.subctx[2].algorithm.problem.g)
-    display(ctx.subctx[1].p)
-    display(ctx.subctx[2].p)
-    display(recover_u(p, prob))
+    lognan(x) = x > 0 ? x : NaN
+
+    #y = [
+    #     lognan.(alg_energy(alg_ref, n) .- ref_energy),
+    #     lognan.(alg_energy(alg_dd, n) .- ref_energy),
+    #     lognan.(alg_energy(alg_dd2, n) .- ref_energy),
+    #   ]
+    y = [
+         lognan.(alg_error(alg_ref, pmin, n, 10)),
+         lognan.(alg_error(alg_dd, pmin, n)),
+         lognan.(alg_error(alg_dd2, pmin, n)),
+       ]
+
+    plot(y, xaxis=:log, yaxis=:log)
+
+
+    #display(energy(ctx))
+    #display(ctx.p)
+    #display(recover_u(p, prob))
 
     #println(energy(ctx))
     #println(energy(ctx2))
@@ -139,40 +182,46 @@ function run3()
     ctx, ctx2
 end
 
+function energy(ctx::Union{DualTVDDState,ProjGradState,ChambolleState})
+    return energy(ctx.p, ctx.algorithm.problem)
+end
+
+function error(p, pmin, prob::DualTVL1ROFOpProblem)
+    return energy_norm(p .- pmin, prob) / energy(pmin, prob)
+end
 
-function energy(ctx::Union{DualTVDDState,ProjGradState})
-    d = ndims(ctx.p)
+function energy(p, prob::DualTVL1ROFOpProblem)
+    d = ndims(p)
 
     @inline kfΛ(w) = @inbounds divergence(w)
     kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
 
-    v = similar(ctx.g)
+    v = similar(prob.g)
 
-    # v = div(p) + A'*f
-    map!(kΛ, v, extend(ctx.p, StaticKernels.ExtensionNothing()))
-    v .+= ctx.g
+    # v = div(p) + g
+    map!(kΛ, v, extend(p, StaticKernels.ExtensionNothing()))
+    v .+= prob.g
     #display(v)
 
     # |v|_B^2 / 2
-    u = ctx.B * vec(v)
+    u = prob.B * vec(v)
     return sum(u .* vec(v)) / 2
 end
 
-function energy(ctx::ChambolleState)
-    d = ndims(ctx.p)
+function energy_norm(p, prob::DualTVL1ROFOpProblem)
+    d = ndims(p)
 
     @inline kfΛ(w) = @inbounds divergence(w)
     kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
 
-    v = similar(ctx.problem.g)
+    v = similar(prob.g)
 
-    # v = div(p) + g
-    map!(kΛ, v, extend(ctx.p, StaticKernels.ExtensionNothing()))
-    v .+= ctx.problem.g
+    # v = div(p)
+    map!(kΛ, v, extend(p, StaticKernels.ExtensionNothing()))
     #display(v)
 
     # |v|_B^2 / 2
-    u = ctx.problem.B * vec(v)
+    u = prob.B * vec(v)
     return sum(u .* vec(v)) / 2
 end
 
diff --git a/src/chambolle.jl b/src/chambolle.jl
index 0edc3039f30d5a05e6c2ac0cf2c67fb44682bfe3..0a183cd803de5f663db841001fdb7c40c0ece8b9 100644
--- a/src/chambolle.jl
+++ b/src/chambolle.jl
@@ -18,15 +18,14 @@ struct ChambolleAlgorithm{P} <: Algorithm{P}
 
     "fixed point inertia parameter"
     τ::Float64
+    "anisotropic projection"
+    anisotropic::Bool
 
-    function ChambolleAlgorithm(problem; τ)
-        return new{typeof(problem)}(problem, τ)
+    function ChambolleAlgorithm(problem; τ=inv(8 * normB(problem)), anisotropic=false)
+        return new{typeof(problem)}(problem, τ, anisotropic)
     end
 end
 
-ChambolleAlgorithm(problem::DualTVL1ROFOpProblem) =
-    ChambolleAlgorithm(problem, τ=inv(8 * normB(problem)))
-
 struct ChambolleState{A,T,R,S,Sv,K1,K2} <: State
     algorithm::A
 
@@ -49,12 +48,9 @@ struct ChambolleState{A,T,R,S,Sv,K1,K2} <: State
     k2::K2
 end
 
-function proj(ctx::ChambolleState, x, s, λ)
-    @inbounds iszero(λ) && return zero(p)
-    return @inbounds p ./ (1 + norm(sgrad) / md.λ[sw.position])
-end
+projnorm(v, anisotropic) = anisotropic ? abs.(v) : norm(v)
 
-function Base.iterate(alg::ChambolleAlgorithm{<:DualTVL1ROFOpProblem})
+function init(alg::ChambolleAlgorithm{<:DualTVL1ROFOpProblem})
     g = alg.problem.g
     λ = alg.problem.λ
     d = ndims(g)
@@ -71,18 +67,19 @@ function Base.iterate(alg::ChambolleAlgorithm{<:DualTVL1ROFOpProblem})
     @inline function kf2(sw)
         @inbounds iszero(λ[sw.position]) && return zero(p[sw.position])
         sgrad = alg.τ * gradient(sw)
-        return @inbounds (p[sw.position] + sgrad) / (1 + norm(sgrad) / λ[sw.position])
+        return @inbounds (p[sw.position] + sgrad) ./ (1 .+ projnorm(sgrad, alg.anisotropic) ./ λ[sw.position])
     end
     k2 = Kernel{ntuple(_->0:1, d)}(kf2)
 
-    return (p, ChambolleState(alg, p, r, s, rv, sv, k1, k2))
+    return ChambolleState(alg, p, r, s, rv, sv, k1, k2)
 end
 
 function reset!(ctx::ChambolleState)
     fill!(ctx.p, zero(eltype(ctx.p)))
 end
 
-function Base.iterate(alg::ChambolleAlgorithm{<:DualTVL1ROFOpProblem}, ctx)
+function step!(ctx::ChambolleState)
+    alg = ctx.algorithm
     # r = div(p) + g
     map!(ctx.k1, ctx.r, ctx.p)
     # s = B * r
@@ -90,34 +87,24 @@ function Base.iterate(alg::ChambolleAlgorithm{<:DualTVL1ROFOpProblem}, ctx)
     # p = (p + τ*grad(s)) / (1 + τ/λ|grad(s)|)
     map!(ctx.k2, ctx.p, ctx.s)
 
-    return (ctx.p, ctx)
+    return ctx
 end
 
-recover_u(ctx::ChambolleState) = recover_u(ctx.p, ctx.problem)
+fetch(ctx::ChambolleState) =
+    recover_u(ctx.p, ctx.algorithm.problem)
 
-#
-#function solve(md::ROFProblem, alg::Chambolle;
-#               maxiters = 1000,
-#               save_log = false)
-#
-#    p = zeros(ndims(g), size(g)...)
-#    q = zero(p)
-#    s = zero(g)
-#    div = div_op(p)
-#    grad = grad_op(g)
-#    ctx = ChambolleState(p, q, s, div, grad)
-#
-#    log = Log()
-#
-#    for i in 1:maxiters
-#        step!(ctx, md, alg)
-#
-#        if save_log
-#            recover_u!(ctx, md, alg)
-#            save_log && push!(log.energy, energy2(md, ctx.s))
-#        end
-#    end
-#
-#    recover_u!(ctx, md, alg)
-#    return Solution(ctx.s, md, alg, ctx, log)
-#end
+function recover_u(p, md::DualTVL1ROFOpProblem)
+    d = ndims(md.g)
+    u = similar(md.g)
+    v = similar(md.g)
+
+    @inline kfΛ(w) = @inbounds divergence(w)
+    kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
+
+    # v = div(p) + A'*f
+    map!(kΛ, v, extend(p, StaticKernels.ExtensionNothing())) # extension: nothing
+    v .+= md.g
+    # u = B * v
+    mul!(vec(u), md.B, vec(v))
+    return u
+end
diff --git a/src/common.jl b/src/common.jl
index 3643599730de461427576ff1007d496a5f468d36..29a2d9e116b1ef0db307eb335051e2ab222088cb 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -1,25 +1,6 @@
 using StaticArrays
 using StaticKernels
 
-# Solver Interface Notes:
-#
-# - a concrete subtype `<:Problem` describes one model type and each instance
-#   represents a fully specified problem, i.e. including data and model parameters
-# - a concrete subtype `<:Algorithm` describes an algorithm, constructed from a
-#   model and algoritm parameters as keyword arguments
-# - a concrete subtype `<:Context` contains preallocated data for the
-#   algorithm
-#
-# - `init(::Problem, ::Algorithm)::Context` initializes the algorithm context
-# - `iterate!(::Context)::Union{Context,Nothing}` performs one non-allocating step
-#
-# - iterating over the same context must lead to the same deterministic
-#   sequence of iterates
-# - `run(::Context)` must continue on from the given context
-# - `(<:Algorithm)(...)` constructors must accept keyword arguments only
-# - types `<:Context` satisfy a stateful iterator interface, returning the
-#   current context at each step of the algorithm
-
 """
     <:Problem
 
@@ -33,8 +14,8 @@ abstract type Problem end
     (<:Algorithm)(::Problem; params...)
 
 An algorithm represents a fully specified iterative process of usually infite
-or a-priori unknown finite length, producing iterates that approximate the
-solution to some problem (see `::Problem`).
+or a-priori fixed but unknown finite length, producing iterates that
+approximate the solution to some problem (see `::Problem`).
 
 The hierarchy of abstract types `<:Algorithm` is based on the algorithm
 interface (e.g. accepted parameters) and the specific numerical scheme.
@@ -44,44 +25,35 @@ algorithm. An instance is a complete specification of the algorithm with all
 its inputs and should guarantee to produce a deterministic sequence of iterates
 (randomized algorithms should use a seeded pseudorandom number generator).
 
-An algorithm needs to implement
+An concrete algorithm type needs to implement at least
+
+  - `state = init(::Algorithm)` allocates and initializes the algorithm state.
+  - `state = step!(::State)` performs one iteration of the algorithm by
+    updating `state`. This method should not allocate dynamic memory.
+  - `solution = fetch(::State)` fetches the solution in a format specified by
+    the problem. In some cases this can be a non-trivial operation if the
+    algorithm uses a different internal representation.
 
-  - `(x0, state) = Base.iterate(::Algorithm)` allocates and initializes the
-    algorithm state and returns the initial iterate `x0`.
-  - `(x, state) = Base.iterate(::Algorithm, state)` performs one iteration of the
-    algorithm by updating `state` and returning the updated iterate `x`. This
-    method should not allocate dynamic memory.
+Depending on its supertype algorithm interface additional methods might need to
+be implemented.
 """
 abstract type Algorithm{P<:Problem} end
 
-Base.IteratorSize(::Type{<:Algorithm}) = Base.SizeUnknown()
-
-function execute(alg::Algorithm; maxiters=nothing)
-    k = 0
-    y = iterate(alg)
-    @assert !isnothing(y)
-    (x, state) = y
-    while true
-        k += 1
-        k > maxiters && break
-        y = iterate(alg, state)
-        isnothing(y) && break
-        (x, state) = y
-    end
-    return x
-end
+"""
+    <:State
 
+Algorithm state and storage containing sufficient information to act as a
+checkpoint for continuing the corresponding algorithm at that point.
+"""
 abstract type State end
 
+function init end
+function step! end
+function fetch end
 
-"helper struct to prevent allocations"
-struct AppliedLinearKernel{K,A}
-    kernel::K
-    data::A
-end
-
-compute!(dst, op::AppliedLinearKernel) = map!(dst, op.kernel, op.a)
 
+Base.intersect(a::CartesianIndices{d}, b::CartesianIndices{d}) where d =
+    CartesianIndices(intersect.(a.indices, b.indices))
 
 @generated function gradient(w::StaticKernels.Window{<:Any,N}) where N
     i0 = ntuple(_->0, N)
diff --git a/src/dualtvdd.jl b/src/dualtvdd.jl
index 456ebae16f5bf15b351653609155802c681c11e3..b54df22cb5bd48975da28f4d564730ad32137783 100644
--- a/src/dualtvdd.jl
+++ b/src/dualtvdd.jl
@@ -5,10 +5,14 @@ struct DualTVDDAlgorithm{P,d} <: Algorithm{P}
     M::NTuple{d,Int}
     "overlap in pixels per dimension"
     overlap::NTuple{d,Int}
-    "inertia parameter"
+    "use the non-sequential inertia update"
+    parallel::Bool
+    "inertia parameter (only when parallel)"
     σ::Float64
-    function DualTVDDAlgorithm(problem; M, overlap, σ)
-        return new{typeof(problem), length(M)}(problem, M, overlap, σ)
+    "number of inner iterations"
+    ninner::Int
+    function DualTVDDAlgorithm(problem; M, overlap, parallel=true, σ=1/4, ninner=10)
+        return new{typeof(problem), length(M)}(problem, M, overlap, parallel, σ, ninner)
     end
 end
 
@@ -57,41 +61,57 @@ function Base.iterate(alg::DualTVDDAlgorithm{<:DualTVL1ROFOpProblem})
     return p, DualTVDDState(alg, p, q, subax, subctx)
 end
 
+function intersectin(a, b)
+    c = intersect(a, b)
+    az = first(a) - one(eltype(a))
+    bz = first(b) - one(eltype(b))
+    return (c, c .- az, c .- bz)
+end
+
 function Base.iterate(alg::DualTVDDAlgorithm{<:DualTVL1ROFOpProblem}, ctx)
-    σ = ctx.algorithm.σ
+    # σ = 1 takes care of sequential updates
+    σ = alg.parallel ? ctx.algorithm.σ : 1.
     d = ndims(ctx.p)
     ax = axes(ctx.p)
     overlap = ctx.algorithm.overlap
 
     # call run! on each cell (this can be threaded)
     for (i, sax) in pairs(ctx.subax)
+        li = LinearIndices(ctx.subax)[i]
         sg = ctx.subctx[i].algorithm.problem.g # julia-bug workaround
         sq = ctx.q[i] # julia-bug workaround
 
-        #println("# subax $sax")
-
         sg .= alg.problem.g[sax...]
-        sq .= (1 .- theta.(Ref(ax), Ref(sax), Ref(overlap), CartesianIndices(sax))) .* ctx.p[sax...]
+        if alg.parallel
+            sq .= (1 .- theta.(Ref(ax), Ref(sax), Ref(overlap), CartesianIndices(sax))) .* ctx.p[sax...]
+        else
+            sq .= Ref(zero(eltype(sq)))
+            # contributions from previous domains
+            for (lj, saxj) in enumerate(ctx.subax)
+                ids, idsi, idsj = intersectin(CartesianIndices(sax), CartesianIndices(saxj))
+                if lj < li
+                    sq[idsi] .+= ctx.subctx[lj].p[idsj]
+                elseif lj > li
+                    sq[idsi] .+= theta.(Ref(ax), Ref(saxj), Ref(overlap), ids) .* ctx.p[ids]
+                end
+            end
+        end
 
         @inline kfΛ(pw) = @inbounds sg[pw.position] + divergence(pw)
         kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
 
         map!(kΛ, sg, sq)
 
-        for j in 1:10
+        for j in 1:alg.ninner
             (_, ctx.subctx[i]) = iterate(ctx.subctx[i].algorithm, ctx.subctx[i])
         end
-        #display(ctx.subctx[i].algorithm.problem.g .- alg.problem.g[sax...])
-        #display(ctx.subctx[i].p)
     end
-    #display(ctx.p)
 
     # aggregate (not thread-safe!)
     ctx.p .*= 1 - σ
     for (i, sax) in pairs(ctx.subax)
         ctx.p[sax...] .+= σ .* ctx.subctx[i].p
     end
-    #display(ctx.p)
 
     return ctx.p, ctx
 end
diff --git a/src/problems.jl b/src/problems.jl
index 7e7fc7d7fd27a5f792c2c78dedce97db28be85d2..6f6a81cbe00eb3aa87f1835923ed28dfe6d4c340 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -18,19 +18,3 @@ DualTVL1ROFOpProblem(g, B, λ::Real) = DualTVL1ROFOpProblem(g, B, fill!(similar(
 
 normB(p::DualTVL1ROFOpProblem{<:Any,<:UniformScaling}) = p.B.λ * sqrt(length(p.g))
 normB(p::DualTVL1ROFOpProblem) = norm(p.B)
-
-function recover_u(p, md::DualTVL1ROFOpProblem)
-    d = ndims(md.g)
-    u = similar(md.g)
-    v = similar(md.g)
-
-    @inline kfΛ(w) = @inbounds divergence(w)
-    kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
-
-    # v = div(p) + A'*f
-    map!(kΛ, v, extend(p, StaticKernels.ExtensionNothing())) # extension: nothing
-    v .+= md.g
-    # u = B * v
-    mul!(vec(u), md.B, vec(v))
-    return u
-end
diff --git a/src/projgrad.jl b/src/projgrad.jl
index 1de650342270361d4870191dd2624c8aabc06c50..6bbc2ed80a07e34894c37c88c8813852798516e2 100644
--- a/src/projgrad.jl
+++ b/src/projgrad.jl
@@ -30,13 +30,13 @@ struct ProjGradState{A,Tp,Wv,R,S,K1,K2}
     k2::K2
 end
 
-function Base.iterate(alg::ProjGradAlgorithm{<:DualTVL1ROFOpProblem})
+function init(alg::ProjGradAlgorithm{<:DualTVL1ROFOpProblem})
     g = alg.problem.g
     d = ndims(g)
     ax = axes(g)
 
     p = extend(zeros(SVector{d,Float64}, ax), StaticKernels.ExtensionNothing())
-    g = extend(g, StaticKernels.ExtensionNothing())
+    ge = extend(g, StaticKernels.ExtensionNothing())
 
     rv = zeros(length(g))
     sv = zeros(length(g))
@@ -46,7 +46,7 @@ function Base.iterate(alg::ProjGradAlgorithm{<:DualTVL1ROFOpProblem})
 
     z = zero(CartesianIndex{d})
 
-    @inline kf1(pw) = @inbounds -divergence(pw) - g[pw.position]
+    @inline kf1(pw) = @inbounds -divergence(pw) - ge[pw.position]
     k1 = Kernel{ntuple(_->-1:1, d)}(kf1)
 
     @inline function kf2(pw, sw)
@@ -55,10 +55,11 @@ function Base.iterate(alg::ProjGradAlgorithm{<:DualTVL1ROFOpProblem})
     end
     k2 = Kernel{ntuple(_->0:1, d)}(kf2)
 
-    return (p, ProjGradState(alg, p, rv, sv, r, s, k1, k2))
+    return ProjGradState(alg, p, rv, sv, r, s, k1, k2)
 end
 
-function Base.iterate(alg::ProjGradAlgorithm{<:DualTVL1ROFOpProblem}, ctx)
+function step!(ctx::ProjGradState)
+    alg = ctx.algorithm
     # r = Λ*p - g
     map!(ctx.k1, ctx.r, ctx.p)
     # s = B * r
@@ -66,7 +67,7 @@ function Base.iterate(alg::ProjGradAlgorithm{<:DualTVL1ROFOpProblem}, ctx)
     # p = proj(p - λΛ's)
     map!(ctx.k2, ctx.p, ctx.p, ctx.s)
 
-    return (ctx.p, ctx)
+    return ctx
 end
 
-recover_u(ctx::ProjGradState) = recover_u(ctx.p, ctx.problem)
+fetch(ctx::ProjGradState) = recover_u(ctx.p, ctx.algorithm.problem)
diff --git a/test/runtests.jl b/test/runtests.jl
index 9e8d0d08f5096a7b895bbb0072e4003abf77b290..fa69c23dd3914e7c776f63fd9ff2e1f315db00fa 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,13 +1,18 @@
-using Test
+using Test, BenchmarkTools
 using LinearAlgebra
 using DualTVDD:
     DualTVL1ROFOpProblem, ProjGradAlgorithm, ChambolleAlgorithm,
-    init, step!, recover_u
+    init, step!, fetch
 
 g = Float64[0 2; 1 0]
 prob = DualTVL1ROFOpProblem(g, I, 1e-10)
 
-for alg in (ProjGradAlgorithm(prob, τ=1/8), ChambolleAlgorithm(prob))
-    its = collect(Iterators.take(alg, 100))
-    @test recover_u(last(its)) ≈ g
+@testset for alg in (ProjGradAlgorithm(prob, τ=1/8), ChambolleAlgorithm(prob))
+    ctx = init(alg)
+    @test 0 == @ballocated step!($ctx)
+    for i in 1:100
+        step!(ctx)
+    end
+    u = fetch(ctx)
+    @test u ≈ g
 end