diff --git a/Changelog.md b/Changelog.md
index 776ce2b841..d424e36434 100644
--- a/Changelog.md
+++ b/Changelog.md
@@ -19,9 +19,14 @@ and their documentation and testing has been extendede.
 ### Changed
 
 * Bumped and added dependencies on all 3 Project.toml files, the main one, the docs/, an the tutorials/ one.
-* `costIntrICTV12` is now available as [`ManoptExamples.Intrinsic_infimal_convolution_TV12`]()
-* `costL2TV` is not available as `ManoptExamples.L2_Total_Variation`
-
+* `costIntrICTV12` is available as [`ManoptExamples.Intrinsic_infimal_convolution_TV12`]()
+* `grad_Intrinsic_infimal_convolution_TV12` is available as [`ManoptExamples.Intrinsic_infimal_convolution_TV12``]()
+* `costL2TV` is available as [`ManoptExamples.L2_Total_Variation`]()
+* `costL2TV12` is available as [`ManoptExamples.L2_Total_Variation_1_2`]()
+* `costL2TV2` is available as [`ManoptExamples.L2_second_order_Total_Variation`]()
+* `costTV` is available as [`ManoptExamples.Total_Variation`]()
+* `costTV2` is available as [`ManoptExamples.second_order_Total_Variation`]()
+*
 ## [0.4.43] – November 19, 2023
 
 ### Added
diff --git a/Project.toml b/Project.toml
index 758df65858..f616d432a1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manopt"
 uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
 authors = ["Ronny Bergmann <manopt@ronnybergmann.net>"]
-version = "0.4.43"
+version = "0.4.44"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
@@ -47,6 +47,7 @@ LRUCache = "1.4"
 ManifoldDiff = "0.3.8"
 Manifolds = "0.9"
 ManifoldsBase = "0.15"
+ManoptExamples = "0.1.4"
 Markdown = "1.6"
 PolynomialRoots = "1"
 Printf = "1.6"
@@ -62,9 +63,11 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
+ManifoldDiff = "af67fdf4-a580-4b9f-bbec-742ef357defd"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
+ManoptExamples = "5b8d5e80-5788-45cb-83d6-5e8f1484217d"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ForwardDiff", "JuMP", "Manifolds", "Plots", "LineSearches", "LRUCache"]
+test = ["Test", "ForwardDiff", "JuMP", "Manifolds", "ManoptExamples", "ManifoldDiff", "Plots", "LineSearches", "LRUCache"]
diff --git a/src/Manopt.jl b/src/Manopt.jl
index f669dcbaa7..5f5f940bea 100644
--- a/src/Manopt.jl
+++ b/src/Manopt.jl
@@ -20,20 +20,8 @@ using Dates: Millisecond, Nanosecond, Period, canonicalize, value
 using LinearAlgebra:
     Diagonal, I, eigen, eigvals, tril, Symmetric, dot, cholesky, eigmin, opnorm
 using ManifoldDiff:
-    adjoint_Jacobi_field,
-    adjoint_Jacobi_field!,
-    adjoint_differential_exp_argument,
-    adjoint_differential_exp_argument!,
-    adjoint_differential_exp_basepoint,
-    adjoint_differential_exp_basepoint!,
     adjoint_differential_log_argument,
     adjoint_differential_log_argument!,
-    adjoint_differential_log_basepoint,
-    adjoint_differential_log_basepoint!,
-    adjoint_differential_shortest_geodesic_endpoint,
-    adjoint_differential_shortest_geodesic_endpoint!,
-    adjoint_differential_shortest_geodesic_startpoint,
-    adjoint_differential_shortest_geodesic_startpoint!,
     differential_exp_argument,
     differential_exp_argument!,
     differential_exp_basepoint,
diff --git a/test/functions/test_adjoint_differentials.jl b/test/functions/test_adjoint_differentials.jl
deleted file mode 100644
index bba27ec616..0000000000
--- a/test/functions/test_adjoint_differentials.jl
+++ /dev/null
@@ -1,25 +0,0 @@
-using Manifolds, Manopt, Test, ManifoldsBase
-
-@testset "Differentials (on Sphere(2))" begin
-    # The Adjoint Differentials test using the same variables as the differentials
-    # test
-    p = [1.0, 0.0, 0.0]
-    q = [0.0, 1.0, 0.0]
-    M = Sphere(2)
-    X = log(M, p, q)
-    # Tept differentials (1) Dp of Log_pq
-    Y = similar(X)
-    Mp = PowerManifold(M, NestedPowerRepresentation(), 3)
-    pP = [p, q, p]
-    qP = [p, p, q]
-    XP = [X, zero_vector(M, p), -X]
-    YP = similar.(XP)
-    ZP = adjoint_differential_forward_logs(Mp, pP, XP)
-    @test norm(Mp, pP, ZP - [-X, X, zero_vector(M, p)]) ≈ 0 atol = 4 * 10.0^(-16)
-    adjoint_differential_forward_logs!(Mp, YP, pP, XP)
-    @test isapprox(Mp, pP, YP, ZP)
-    ZP = [[0.0, π / 2, 0.0], [0.0, 0.0, 0.0], [π / 2, 0.0, 0.0]]
-    @test Manopt.adjoint_differential_log_argument(Mp, pP, qP, XP) == ZP
-    Manopt.adjoint_differential_log_argument!(Mp, YP, pP, qP, XP)
-    @test ZP == YP
-end
diff --git a/test/functions/test_bezier.jl b/test/functions/test_bezier.jl
deleted file mode 100644
index 342156dfd3..0000000000
--- a/test/functions/test_bezier.jl
+++ /dev/null
@@ -1,186 +0,0 @@
-using Manopt, Manifolds, Test
-
-@testset "Bezier Tests" begin
-    @testset "General Bezier Tests" begin
-        repr(BezierSegment([[0.0, 0.0], [0.0, 0.0]])) ==
-        "BezierSegment([[0.0, 0.0], [0.0, 0.0]])"
-    end
-    @testset "Spherical Test" begin
-        M = Sphere(2)
-        pC = [0.0, 1.0, 0.0]
-        pT = exp(M, pC, [0.0, 0.0, 0.7])
-        pB = exp(M, pC, [0.0, 0.0, -0.7])
-        B = [
-            BezierSegment(shortest_geodesic(M, pT, pC, [0.0, 1 / 3, 2 / 3, 1.0])),
-            BezierSegment(shortest_geodesic(M, pC, pB, [0.0, 1 / 3, 2 / 3, 1.0])),
-        ]
-        # this is equispaced, so the pure cost is zero and the gradient is a zero-vector
-        t = collect(range(0.0, 1.0; length=5))
-        pts = shortest_geodesic(M, pT, pB, t)
-        pts2 = de_casteljau(M, B, 2 .* t)
-        @test sum(distance.(Ref(M), pts, pts2)) < 10 * eps()
-        aX = log(M, pT, pC)
-        aT1 = adjoint_differential_bezier_control(M, BezierSegment([pT, pC]), 0.5, aX).pts
-        aT1a = BezierSegment(similar.(aT1))
-        adjoint_differential_bezier_control!(M, aT1a, BezierSegment([pT, pC]), 0.5, aX)
-        @test aT1a.pts == aT1
-        aT2 = [
-            Manopt.adjoint_differential_shortest_geodesic_startpoint(M, pT, pC, 0.5, aX),
-            Manopt.adjoint_differential_shortest_geodesic_endpoint(M, pT, pC, 0.5, aX),
-        ]
-        @test aT1 ≈ aT2
-        #
-        @test sum(
-            norm.(
-                grad_acceleration_bezier(
-                    M, B[1], collect(range(0.0, 1.0; length=20))
-                ).pts .-
-                [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
-            ),
-        ) ≈ 0 atol = 2e-12
-
-        # cost and gradient
-        T = collect(range(0.0, 2.0; length=51))
-        degrees = get_bezier_degrees(M, B)
-        Bvec = get_bezier_points(M, B, :differentiable)
-        Mp = PowerManifold(M, NestedPowerRepresentation(), length(Bvec))
-        @test cost_acceleration_bezier(M, Bvec, degrees, T) ≈ 0 atol = 10^-10
-        z = zero_vector(Mp, Bvec)
-        distance(Mp, grad_acceleration_bezier(M, Bvec, degrees, T), z)
-        @test norm(Mp, Bvec, grad_acceleration_bezier(M, Bvec, degrees, T) - z) ≈ 0 atol =
-            2e-12
-
-        d = [pT, exp(M, pC, [0.3, 0.0, 0.0]), pB]
-        λ = 3.0
-
-        # cost and gradient with data term
-        @test cost_L2_acceleration_bezier(M, Bvec, degrees, T, λ, [pT, pC, pB]) ≈ 0 atol =
-            10^(-10)
-        @test cost_L2_acceleration_bezier(M, Bvec, degrees, T, λ, d) ≈
-            λ / 2 * distance(M, d[2], pC) .^ 2
-        # when the data are the junctions
-        @test norm(
-            Mp, Bvec, grad_L2_acceleration_bezier(M, Bvec, degrees, T, λ, [pT, pC, pB]) - z
-        ) ≈ 0 atol = 2e-12
-        z[4][1] = -0.9
-        @test norm(Mp, Bvec, grad_L2_acceleration_bezier(M, Bvec, degrees, T, λ, d) - z) ≈ 0 atol =
-            2e-12
-        # when the data is weighted with zero
-        @test cost_L2_acceleration_bezier(M, Bvec, degrees, T, 0.0, d) ≈ 0 atol = 10^(-10)
-        z[4][1] = 0.0
-        @test norm(Mp, Bvec, grad_L2_acceleration_bezier(M, Bvec, degrees, T, 0.0, d) - z) ≈
-            0 atol = 2e-12
-    end
-    @testset "de Casteljau variants" begin
-        M = Sphere(2)
-        B = artificial_S2_composite_bezier_curve()
-        b = B[2]
-        b2s = BezierSegment([b.pts[1], b.pts[end]])
-        # (a) 2 points -> geo special case
-        f1 = de_casteljau(M, b2s) # fct -> recursive
-        pts1 = f1.([0.0, 0.5, 1.0])
-        pts2 = de_casteljau(M, b2s, [0.0, 0.5, 1.0])
-        @test pts1 ≈ pts2
-        # (b) one segment
-        f2 = de_casteljau(M, b) # fct -> recursive
-        pts3 = f2.([0.0, 0.5, 1.0])
-        pts4 = de_casteljau(M, b, [0.0, 0.5, 1.0])
-        @test pts3 ≈ pts4
-        # (c) whole composites
-        f3 = de_casteljau(M, B) # fct -> recursive
-        pts5 = f3.([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0])
-        pts6 = de_casteljau(M, B, [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0])
-        @test pts5 ≈ pts6
-        @test_throws DomainError f3(-0.1)
-        @test_throws DomainError de_casteljau(M, B, -0.1)
-        @test_throws DomainError f3(3.5)
-        @test_throws DomainError de_casteljau(M, B, 3.5)
-    end
-    @testset "Spherical Data" begin
-        M = Sphere(2)
-        B = artificial_S2_composite_bezier_curve()
-        @test de_casteljau(M, B, [0.0, 1.0, 2.0, 3.0]) ≈
-            [B[1].pts[1], B[2].pts[1], B[3].pts[1], B[3].pts[4]]
-        @test get_bezier_junction_tangent_vectors(M, B) ≈ [
-            log(M, B[1].pts[1], B[1].pts[2]),
-            log(M, B[1].pts[4], B[1].pts[3]),
-            log(M, B[2].pts[1], B[2].pts[2]),
-            log(M, B[2].pts[4], B[2].pts[3]),
-            log(M, B[3].pts[1], B[3].pts[2]),
-            log(M, B[3].pts[4], B[3].pts[3]),
-        ]
-        @test get_bezier_junction_tangent_vectors(M, B[1]) ≈
-            [log(M, B[1].pts[1], B[1].pts[2]), log(M, B[1].pts[4], B[1].pts[3])]
-        @test get_bezier_junctions(M, B[1]) == B[1].pts[[1, end]]
-        @test get_bezier_inner_points(M, B) ==
-            [B[1].pts[2], B[1].pts[3], B[2].pts[2], B[2].pts[3], B[3].pts[2], B[3].pts[3]]
-        @test get_bezier_inner_points(M, B[1]) == [B[1].pts[2], B[1].pts[3]]
-
-        @test get_bezier_points(M, B) == cat([[b.pts...] for b in B]...; dims=1)
-        @test get_bezier_points(M, B, :continuous) ==
-            cat([[b.pts[[1:3]...]...] for b in B]..., [B[3].pts[4]]; dims=1)
-        @test get_bezier_points(M, B, :differentiable) ==
-            cat([B[1].pts[[1, 2]]...], [b.pts[[3, 4]] for b in B]...; dims=1)
-        @test get_bezier_points(M, B[1]) == B[1].pts
-        # for segments just check that they
-        d = get_bezier_degrees(M, B)
-        A = get_bezier_segments(M, get_bezier_points(M, B), d)
-        @test [A[i].pts for i in 1:3] == [B[i].pts for i in 1:3]
-        A = get_bezier_segments(M, get_bezier_points(M, B, :continuous), d, :continuous)
-        @test [A[i].pts for i in 1:3] == [B[i].pts for i in 1:3]
-        A = get_bezier_segments(
-            M, get_bezier_points(M, B, :differentiable), d, :differentiable
-        )
-        @test [A[i].pts for i in 1:3] == [B[i].pts for i in 1:3]
-
-        # out of range
-        @test_throws ErrorException adjoint_differential_bezier_control(
-            M, B, 7.0, zero_vector(M, B[1].pts[1])
-        )
-        # a shortcut to evaluate the adjoint at several points is equal to separate evals
-        b = B[2]
-        Xi = [log(M, b.pts[1], b.pts[2]), -log(M, b.pts[4], b.pts[3])]
-        Xs = adjoint_differential_bezier_control(M, b, [0.0, 1.0], Xi)
-        @test isapprox(
-            Xs.pts,
-            adjoint_differential_bezier_control(M, b, 0.0, log(M, b.pts[1], b.pts[2])).pts +
-            adjoint_differential_bezier_control(M, b, 1.0, -log(M, b.pts[4], b.pts[3])).pts,
-        )
-        Ys = BezierSegment(similar.(Xs.pts))
-        adjoint_differential_bezier_control!(M, Ys, b, [0.0, 1.0], Xi)
-        @test isapprox(Xs.pts, Ys.pts)
-        # differential
-        X = BezierSegment([
-            log(M, b.pts[1], b.pts[2]), [zero_vector(M, b.pts[i]) for i in 2:4]...
-        ])
-        Ye = zero(X.pts[1])
-        @test differential_bezier_control(M, b, 0.0, X) ≈ X.pts[1]
-        differential_bezier_control!(M, Ye, b, 0.0, X)
-        @test Ye ≈ X.pts[1]
-        dT1 = differential_bezier_control.(Ref(M), Ref(b), [0.0, 1.0], Ref(X))
-        dT2 = differential_bezier_control(M, b, [0.0, 1.0], X)
-        dT3 = similar.(dT2)
-        differential_bezier_control!(M, dT3, b, [0.0, 1.0], X)
-        @test dT1 ≈ dT2
-        @test dT3 == dT3
-        X2 = [
-            BezierSegment([[0.0, 0.0, 0.0] for i in 1:4]),
-            X,
-            BezierSegment([[0.0, 0.0, 0.0] for i in 1:4]),
-        ]
-        @test_throws DomainError differential_bezier_control(M, B, 20.0, X2)
-        dbT2a = differential_bezier_control(M, B, 1.0, X2)
-        dbT3a = similar(dbT2a)
-        @test_throws DomainError differential_bezier_control!(M, dbT3a, B, 20.0, X2)
-        differential_bezier_control!(M, dbT3a, B, 1.0, X2)
-        @test dbT2a == dbT3a
-        @test dbT2a ≈ X.pts[1]
-        dbT2 = differential_bezier_control(M, B, [1.0, 2.0], X2)
-        dbT1 = differential_bezier_control.(Ref(M), Ref(B), [1.0, 2.0], Ref(X2))
-        @test dT1 ≈ dbT1
-        @test dbT2 ≈ dbT1
-        dbT3 = similar.(dbT2)
-        differential_bezier_control!(M, dbT3, B, [1.0, 2.0], X2)
-        @test dbT2 == dbT3
-    end
-end
diff --git a/test/functions/test_costs.jl b/test/functions/test_costs.jl
deleted file mode 100644
index 684fb748b0..0000000000
--- a/test/functions/test_costs.jl
+++ /dev/null
@@ -1,27 +0,0 @@
-using Manifolds, Manopt, Test, ManifoldsBase
-
-@testset "Test Costs" begin
-    M = Sphere(2)
-    N = PowerManifold(M, NestedPowerRepresentation(), 3, 3)
-    f = repeat([[1.0, 0.0, 0.0]], 3, 3)
-    x = repeat([[1.0, 0.0, 0.0]], 3, 3)
-    x[2, 1] = [0.0, 1.0, 0.0]
-
-    @test costIntrICTV12(N, f, f, f, 0.0, 0.0) == 0.0
-    @test costIntrICTV12(N, f, x, x, 0.0, 0.0) == 1 / 2 * distance(N, x, f)^2
-    @test costIntrICTV12(N, x, x, x, 2.0, 0.5) ≈ costTV2(N, x) + costTV(N, x)
-    @test costIntrICTV12(N, x, x, x, 1.0, 0.0) ≈ costTV2(N, x)
-    @test costIntrICTV12(N, x, x, x, 1.0, 1.0) ≈ costTV(N, x)
-    @test costL2TV2(N, f, 1.0, x) == 1 / 2 * distance(N, f, x)^2 + 1.0 * costTV2(N, x)
-    #
-    @test costL2TV(N, f, 1.0, f) ≈ 0.0
-    @test costL2TV(N, x, 1.0, x) ≈ 3 * π / 2
-    @test costL2TVTV2(N, f, 0.0, 1.0, x) ≈ 1 / 2 * distance(N, x, f)^2 + costTV2(N, x)
-    @test costL2TVTV2(N, f, 1.0, 1.0, x) ≈
-        1 / 2 * distance(N, x, f)^2 + costTV(N, x) + costTV2(N, x)
-
-    @test costTV2(M, Tuple(x[1:3, 1])) ≈ π / 2
-    @test costTV(N, x, 1, 2) ≈ sqrt(5 / 4) * π
-    @test sum(costTV2(N, x, 1, false)) == costTV2(N, x)
-    @test costTV2(N, f, 2) == 0
-end
diff --git a/test/functions/test_differentials.jl b/test/functions/test_differentials.jl
deleted file mode 100644
index 15727e9d9d..0000000000
--- a/test/functions/test_differentials.jl
+++ /dev/null
@@ -1,66 +0,0 @@
-using Manifolds, Manopt, Test, ManifoldsBase
-
-@testset "Differentials" begin
-    p = [1.0, 0.0, 0.0]
-    q = [0.0, 1.0, 0.0]
-    M = Sphere(2)
-    X = log(M, p, q)
-    Y = similar(X)
-    @testset "forward logs" begin
-        N = PowerManifold(M, NestedPowerRepresentation(), 3)
-        x = [p, q, p]
-        y = [p, p, q]
-        V = [X, zero_vector(M, p), -X]
-        Y = Manopt.differential_log_argument(M, p, q, -X)
-        W = similar.(V)
-        @test norm(
-            N,
-            x,
-            differential_forward_logs(N, x, V) - [-X, [π / 2, 0.0, 0.0], zero_vector(M, p)],
-        ) ≈ 0 atol = 8 * 10.0^(-16)
-        differential_forward_logs!(N, W, x, V)
-        @test norm(N, x, W - [-X, [π / 2, 0.0, 0.0], zero_vector(M, p)]) ≈ 0 atol =
-            8 * 10.0^(-16)
-        @test isapprox(N, x, Manopt.differential_log_argument(N, x, y, V), [V[1], V[2], Y])
-        Manopt.differential_log_argument!(N, W, x, y, V)
-        @test isapprox(N, x, W, [V[1], V[2], Y])
-    end
-    @testset "forward logs on a multivariate power manifold" begin
-        S = Sphere(2)
-        M = PowerManifold(S, NestedPowerRepresentation(), 2, 2)
-        p = [zeros(3) for i in [1, 2], j in [1, 2]]
-        p[1, 1] = [1.0, 0.0, 0.0]
-        p[1, 2] = 1 / sqrt(2) .* [1.0, 1.0, 0.0]
-        p[2, 1] = 1 / sqrt(2) .* [1.0, 0.0, 1.0]
-        p[2, 2] = [0.0, 1.0, 0.0]
-        t1 = forward_logs(M, p)
-        @test t1[1, 1, 1] ≈ log(S, p[1, 1], p[2, 1])
-        @test t1[1, 1, 2] ≈ log(S, p[1, 1], p[1, 2])
-        @test t1[1, 2, 1] ≈ log(S, p[1, 2], p[2, 2])
-        @test t1[1, 2, 2] ≈ log(S, p[1, 2], p[1, 2]) atol = 1e-15
-        @test t1[2, 1, 1] ≈ log(S, p[2, 1], p[2, 1]) atol = 1e-15
-        @test t1[2, 1, 2] ≈ log(S, p[2, 1], p[2, 2])
-        @test t1[2, 2, 1] ≈ log(S, p[2, 2], p[2, 2])
-        @test t1[2, 2, 2] ≈ log(S, p[2, 2], p[2, 2])
-        t1a = zero.(t1)
-        forward_logs!(M, t1a, p)
-        @test all(t1 .== t1a)
-        X = zero_vector(M, p)
-        X[1, 1] .= [0.0, 0.5, 0.5]
-        t2 = differential_forward_logs(M, p, X)
-        a =
-            Manopt.differential_log_basepoint(S, p[1, 1], p[2, 1], X[1, 1]) +
-            Manopt.differential_log_argument(S, p[1, 1], p[2, 1], X[2, 1])
-        @test t2[1, 1, 1] ≈ a
-        @test t2[1, 2, 1] ≈ zero_vector(S, p[1, 2]) atol = 1e-17
-        @test t2[2, 1, 1] ≈ zero_vector(S, p[2, 1]) atol = 1e-17
-        @test t2[2, 2, 1] ≈ zero_vector(S, p[2, 2]) atol = 1e-17
-        b =
-            Manopt.differential_log_basepoint(S, p[1, 1], p[1, 2], X[1, 1]) +
-            Manopt.differential_log_argument(S, p[1, 1], p[1, 2], X[1, 2])
-        @test t2[1, 1, 2] ≈ b
-        @test t2[1, 2, 2] ≈ zero_vector(S, p[1, 2]) atol = 1e-17
-        @test t2[2, 1, 2] ≈ zero_vector(S, p[2, 1]) atol = 1e-17
-        @test t2[2, 2, 2] ≈ zero_vector(S, p[2, 2]) atol = 1e-17
-    end
-end
diff --git a/test/functions/test_gradients.jl b/test/functions/test_gradients.jl
deleted file mode 100644
index b639551c10..0000000000
--- a/test/functions/test_gradients.jl
+++ /dev/null
@@ -1,115 +0,0 @@
-using Manifolds, Manopt, Test, ManifoldsBase
-
-@testset "gradients" begin
-    @testset "Circle (Allocating)" begin
-        M = Circle()
-        N = PowerManifold(M, 4)
-        x = [0.1, 0.2, 0.3, 0.5]
-        tvTestξ = [-1.0, 0.0, 0.0, 1.0]
-        @test grad_TV(N, x) == tvTestξ
-        @test grad_TV(M, (x[1], x[1])) == (zero_vector(M, x[1]), zero_vector(M, x[1]))
-        @test norm(N, x, grad_TV(N, x, 2) - tvTestξ) ≈ 0
-        tv2Testξ = [0.0, 1.0, -1.0, 1.0]
-        @test grad_TV2(N, x) == tv2Testξ
-        @test norm(N, x, forward_logs(N, x) - [0.1, 0.1, 0.2, 0.0]) ≈ 0 atol = 10^(-16)
-        @test norm(
-            N,
-            x,
-            grad_intrinsic_infimal_convolution_TV12(N, x, x, x, 1.0, 1.0)[1] -
-            [-1.0, 0.0, 0.0, 1.0],
-        ) ≈ 0
-        @test norm(N, x, grad_intrinsic_infimal_convolution_TV12(N, x, x, x, 1.0, 1.0)[2]) ≈
-            0
-        x2 = [0.1, 0.2, 0.3]
-        N2 = PowerManifold(M, size(x2)...)
-        @test grad_TV2(N2, x2) == zeros(3)
-        @test grad_TV2(N2, x2, 2) == zeros(3)
-        @test grad_TV(M, (0.0, 0.0), 2) == (0.0, 0.0)
-        # 2d forward logs
-        N3 = PowerManifold(M, 2, 2)
-        N3C = PowerManifold(M, 2, 2, 2)
-        x3 = [0.1 0.2; 0.3 0.5]
-        x3C = cat(x3, x3; dims=3)
-        tC = cat([0.2 0.3; 0.0 0.0], [0.1 0.0; 0.2 0.0]; dims=3)
-        @test norm(N3C, x3C, forward_logs(N3, x3) - tC) ≈ 0 atol = 10^(-16)
-
-        M = Circle()
-        p = 0
-        q = π / 4
-        @test grad_distance(M, p, q) == q - p
-        @test grad_distance(M, p, q, 1) == -distance(M, q, p)^(-1) * log(M, q, p)
-        p = q
-        @test grad_distance(M, p, q, 1) == zero_vector(M, p)
-    end
-    @testset "Sphere (Mutating)" begin
-        M = Sphere(2)
-        p = [0.0, 0.0, 1.0]
-        q = [0.0, 1.0, 0.0]
-        r = [1.0, 0.0, 0.0]
-        @testset "Gradient of the distance function" begin
-            X = zero_vector(M, p)
-            grad_distance!(M, X, p, q)
-            Y = grad_distance(M, p, q)
-            Z = [0.0, 0.0, -π / 2] # known solution
-            @test X == Y
-            @test X == Z
-            U = zero_vector(M, p)
-            grad_distance!(M, U, p, q, 1)
-            V = grad_distance(M, p, q, 1)
-            W = -distance(M, q, p)^(-1) * log(M, q, p) # solution
-            @test U == V
-            @test U == W
-            w = q
-            U = zero_vector(M, w)
-            grad_distance!(M, U, w, q, 1)
-            V = grad_distance(M, w, q, 1)
-            W = zero_vector(M, w) # solution
-            @test U == V
-            @test U == W
-        end
-        @testset "Gradient of total variation" begin
-            Y = grad_TV(M, (p, q))
-            Z = [[0.0, -1.0, 0.0], [0.0, 0.0, -1.0]]
-            X = similar.(Z)
-            grad_TV!(M, X, (p, q))
-            @test [y for y in Y] == X
-            @test [y for y in Y] ≈ Z
-            N = PowerManifold(M, NestedPowerRepresentation(), 3)
-            s = [p, q, r]
-            Y2 = grad_TV(N, s)
-            Z2 = [[0.0, -1.0, 0.0], [-1.0, 0.0, -1.0], [0.0, -1.0, 0.0]]
-            X2 = zero_vector(N, s)
-            grad_TV!(N, X2, s)
-            @test Y2 == Z2
-            @test X2 == Z2
-            Y2a = grad_TV(N, s, 2)
-            X2a = zero_vector(N, s)
-            grad_TV!(N, X2a, s, 2)
-            @test Y2a == X2a
-            N2 = PowerManifold(M, NestedPowerRepresentation(), 2)
-            Y3 = grad_TV(M, (p, q), 2)
-            X3 = zero_vector(N2, [p, q])
-            grad_TV!(M, X3, (p, q), 2)
-            @test [y for y in Y3] == X3
-            Y4 = grad_TV(M, (p, p))
-            X4 = zero_vector(N2, [p, q])
-            grad_TV!(M, X4, (p, p))
-            @test [y for y in Y4] == X4
-        end
-        @testset "Grad of second order total variation" begin
-            N = PowerManifold(M, NestedPowerRepresentation(), 3)
-            s = [p, q, r]
-            X = zero_vector(N, s)
-            grad_TV2!(M, X, s)
-            Y = grad_TV2(M, s)
-            Z = -1 / sqrt(2) .* [[0.0, 1.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 0.0]]
-            @test Y == X
-            @test Y ≈ Z
-            Y2 = grad_TV2(M, s, 2)
-            Z2 = -1.110720734539 .* [[0.0, 1.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 0.0]]
-            @test Y2 ≈ Z2
-            s2 = [p, shortest_geodesic(M, p, q, 0.5), q]
-            @test grad_TV2(M, s2) == [zero_vector(M, se) for se in s2]
-        end
-    end
-end
diff --git a/test/functions/test_manifold.jl b/test/functions/test_manifold.jl
deleted file mode 100644
index 675f99e363..0000000000
--- a/test/functions/test_manifold.jl
+++ /dev/null
@@ -1,100 +0,0 @@
-using Manifolds, Manopt, Test, ManifoldsBase
-using LinearAlgebra: I
-
-using Random
-Random.seed!(42)
-# Test the additional manifold functions
-#
-@testset "Additional Manifold functions" begin
-    @testset "mid point & reflect" begin
-        M = Sphere(2)
-        p = [1.0, 0.0, 0.0]
-        q = [0.0, 1.0, 0.0]
-
-        r = mid_point(M, p, q)
-        r2 = similar(r)
-        mid_point!(M, r2, p, q)
-        r3 = shortest_geodesic(M, p, q, 0.5)
-        r4 = mid_point(M, p, q, q)
-        @test isapprox(M, r, r2)
-        @test isapprox(M, r2, r3)
-        @test isapprox(M, r3, r4)
-        r5 = similar(r4)
-        mid_point!(M, r5, p, q, q)
-        @test isapprox(M, r4, r5)
-
-        r4 = mid_point(M, p, -p, q)
-        r5 = similar(r4)
-        mid_point!(M, r5, p, -p, q)
-        @test isapprox(M, r4, q)
-        @test isapprox(M, r4, r5)
-
-        @test isapprox(M, reflect(M, p, q), -q)
-        qA = similar(q)
-        reflect!(M, qA, p, q)
-        @test isapprox(M, qA, -q)
-        f = x -> x
-        @test reflect(M, f, q) == q
-        qA = similar(q)
-        reflect!(M, qA, f, q)
-        @test qA == q
-
-        M2 = Euclidean(2)
-        p2 = [1.0, 0.0]
-        q2 = [0.0, 1.0]
-        s = mid_point(M2, p2, q2)
-        s2 = similar(s)
-        mid_point!(M2, s2, p2, q2)
-        @test s == s2
-        @test s == (p2 + q2) / 2
-        s = mid_point(M2, p2, q2, s)
-        @test s == s2
-        s2 = similar(s)
-        mid_point!(M2, s2, p2, q2, s)
-        @test s == s2
-
-        M = Circle()
-        p = 0
-        q = π
-        @test mid_point(M, p, q, 1.0) ≈ π / 2
-        @test mid_point(M, p, q, -1.0) ≈ -π / 2
-        @test mid_point(M, 0, π / 2) ≈ π / 4
-    end
-
-    @testset "max_stepsize" begin
-        M = Sphere(2)
-        TM = TangentBundle(M)
-        TTM = TangentBundle(TM)
-
-        R3 = Euclidean(3)
-        TR3 = TangentBundle(R3)
-        p = [0.0, 1.0, 0.0]
-        X = [0.0, 0.0, 0.0]
-
-        @test Manopt.max_stepsize(M, p) == π
-        @test Manopt.max_stepsize(TM, ArrayPartition(p, X)) == π
-        @test Manopt.max_stepsize(
-            TTM, ArrayPartition(ArrayPartition(p, X), ArrayPartition(X, X))
-        ) == π
-
-        @test Manopt.max_stepsize(R3, p) == Inf
-        @test Manopt.max_stepsize(TR3, ArrayPartition(p, X)) == Inf
-
-        Mfr = FixedRankMatrices(5, 4, 2)
-        pfr = SVDMPoint(
-            [
-                -0.42232620708727264 0.18201829740358394
-                0.12501774539665936 -0.5154706413711303
-                -0.049666735478056216 0.34842538478418505
-                -0.6185683735016354 -0.6581006457016838
-                -0.6487815663485731 0.38296559722742113
-            ],
-            [0.9240871604723607, 0.488958057530268],
-            [
-                -0.29372512377826565 0.9509130294099503 -0.09069249144195571 0.035564506968098104
-                -0.87290474583133 -0.29359416630673274 -0.3315532983520648 -0.20472464572606347
-            ],
-        )
-        @test Manopt.max_stepsize(Mfr, pfr) == manifold_dimension(Mfr)
-    end
-end
diff --git a/test/functions/test_proximal_maps.jl b/test/functions/test_proximal_maps.jl
deleted file mode 100644
index c44357437a..0000000000
--- a/test/functions/test_proximal_maps.jl
+++ /dev/null
@@ -1,177 +0,0 @@
-using Manifolds, Manopt, Test, Dates
-
-@testset "proximal maps" begin
-    #
-    # prox_TV
-    p = [1.0, 0.0, 0.0]
-    q = [0.0, 1.0, 0.0]
-    M = Sphere(2)
-    N = PowerManifold(M, NestedPowerRepresentation(), 2)
-    @test_throws ErrorException prox_distance(M, 1.0, p, q, 3)
-    @test_throws ErrorException prox_distance!(M, p, 1.0, p, q, 3)
-    @test distance(
-        M, prox_distance(M, distance(M, p, q) / 2, p, q, 1), shortest_geodesic(M, p, q, 0.5)
-    ) < eps()
-    t = similar(p)
-    prox_distance!(M, t, distance(M, p, q) / 2, p, q, 1)
-    @test t == prox_distance(M, distance(M, p, q) / 2, p, q, 1)
-    (r, s) = prox_TV(M, π / 4, (p, q))
-    X = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
-    prox_TV!(M, X, π / 4, (p, q))
-    @test norm(r - s) < eps(Float64)
-    @test norm(X[1] - s) < eps(Float64)
-    @test norm(X[2] - r) < eps(Float64)
-    # i.e. they are moved together
-    @test distance(M, r, s) < eps(Float64)
-    (t, u) = prox_TV(M, π / 8, (p, q))
-    @test_throws ErrorException prox_TV(M, π, (p, q), 3)
-    @test_throws ErrorException prox_TV!(M, [p, q], π, (p, q), 3)
-    # they cross correlate
-    @test (
-        abs(t[1] - u[2]) < eps(Float64) &&
-        abs(t[2] - u[1]) < eps(Float64) &&
-        abs(t[3] - u[3]) < eps(Float64)
-    )
-    @test distance(M, t, u) ≈ π / 4 # and have moved half their distance
-    #
-    (v, w) = prox_TV(M, 1.0, (p, q), 2)
-    vC, wC = shortest_geodesic(M, p, q, [1 / 3, 2 / 3])
-    @test distance(M, v, vC) ≈ 0
-    @test distance(M, w, wC) < eps()
-    P = [similar(p), similar(q)]
-    prox_TV!(M, P, 1.0, (p, q), 2)
-    @test P == [v, w]
-    # prox_TV on Power
-    T = prox_TV(N, π / 8, [p, q])
-    @test distance(N, T, [t, u]) ≈ 0
-    # parallelprox_TV
-    N2 = PowerManifold(M, NestedPowerRepresentation(), 3)
-    r = geodesic(M, p, q, 0.5)
-    s, t = prox_TV(M, π / 16, (r, q))
-    u, v = prox_TV(M, π / 16, (p, r))
-    y = prox_parallel_TV(N2, π / 16, [[p, r, q], [p, r, q]])
-    yM = [similar.(ye) for ye in y]
-    prox_parallel_TV!(N2, yM, π / 16, [[p, r, q], [p, r, q]])
-    @test y == yM
-    @test distance(N2, y[1], [p, s, t]) ≈ 0 # even indices in first comp
-    @test distance(N2, y[2], [u, v, q]) ≈ 0 # odd in second
-    # dimensions of x have to fit, here they don't
-    @test_throws ErrorException prox_parallel_TV(N2, π / 16, [[p, r, q]])
-    @test_throws ErrorException prox_parallel_TV!(N2, yM, π / 16, [[p, r, q]])
-    # prox_TV2
-    p2, r2, q2 = prox_TV2(M, 1.0, (p, r, q))
-    y = [similar(p) for _ in 1:3]
-    prox_TV2!(M, y, 1.0, (p, r, q))
-    @test y ≈ [p2, r2, q2]
-    sum(distance.(Ref(M), [p, r, q], [p2, r2, q2])) ≈ 0
-    @test_throws ErrorException prox_TV2(M, 1.0, (p, r, q), 2) # since prox_TV is only defined for p=1
-    distance(
-        PowerManifold(M, NestedPowerRepresentation(), 3),
-        [p2, r2, q2],
-        prox_TV2(PowerManifold(M, NestedPowerRepresentation(), 3), 1.0, [p, r, q]),
-    ) ≈ 0
-    # Circle
-    M2 = Circle()
-    N2 = PowerManifold(M2, 3)
-    pS, rS, qS = [-0.5, 0.1, 0.5]
-    d = sum([pS, rS, qS] .* [1.0, -2.0, 1.0])
-    m = min(0.3, abs(Manifolds.sym_rem(d) / 6))
-    s = sign(Manifolds.sym_rem(d))
-    pSc, rSc, qSc = Manifolds.sym_rem.([pS, rS, qS] .- m .* s .* [1.0, -2.0, 1.0])
-    pSr, rSr, qSr = prox_TV2(M2, 0.3, (pS, rS, qS))
-    @test sum(distance.(Ref(M2), [pSc, rSc, qSc], [pSr, rSr, qSr])) ≈ 0
-    # p=2
-    t = 0.3 * Manifolds.sym_rem(d) / (1 + 0.3 * 6.0)
-    @test sum(
-        distance.(
-            Ref(M2),
-            [prox_TV2(M2, 0.3, (pS, rS, qS), 2)...],
-            [pS, rS, qS] .- t .* [1.0, -2.0, 1.0],
-        ),
-    ) ≈ 0
-    @test prox_TV2(N2, 0.3, [pS, rS, qS]) == [pSr, rSr, qSr]
-    # others fail
-    @test_throws ErrorException prox_TV2(M2, 0.3, (pS, rS, qS), 3)
-    # Rn
-    M3 = Euclidean(1)
-    pR, rR, qR = [pS, rS, qS]
-    m = min.(Ref(0.3), abs.([pR, rR, qR] .* [1.0, -2.0, 1.0]) / 6)
-    s = sign(d)  # we can reuse d
-    pRc, rRc, qRc = [pR, rR, qR] .- m .* s .* [1.0, -2.0, 1.0]
-    pRr, rRr, qRr = prox_TV2(M3, 0.3, (pR, rR, qR))
-    @test sum(distance.(Ref(M3), [pRc, rRc, qRc], [pRr, rRr, qRr])) ≈ 0
-    # p=2
-    t = 0.3 * d / (1 + 0.3 * 6.0)
-    @test sum(
-        distance.(
-            Ref(M3),
-            [prox_TV2(M3, 0.3, (pR, rR, qR), 2)...],
-            [pR, rR, qR] .- t .* [1.0, -2.0, 1.0],
-        ),
-    ) ≈ 0
-    # others fail
-    @test_throws ErrorException prox_TV2(M3, 0.3, (pR, rR, qR), 3)
-    #
-    # collaborative integer tests
-    #
-    @test_throws ErrorException prox_TV2(M3, 0.3, (pS, rS, qS), 3)
-    ξR, ηR, νR = [pS, rS, qS]
-    N3 = PowerManifold(M3, 3)
-    P = [pR rR qR]
-    Ξ = [ξR ηR νR]
-    Θ = similar(Ξ)
-    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1, 1) == Ξ
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1, 1)
-    @test Θ == Ξ
-    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1.0, 1) == Ξ
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1.0, 1)
-    @test Θ == Ξ
-    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1, 1.0) == Ξ
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1, 1.0)
-    @test Θ == Ξ
-    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1.0, 1.0) == Ξ
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1.0, 1.0)
-    @test Θ == Ξ
-
-    @test project_collaborative_TV(N3, 0.0, P, Ξ, 2, 1) == Ξ
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 2, 1)
-    @test Θ == Ξ
-    @test norm(N3, P, project_collaborative_TV(N3, 0.0, P, Ξ, 2, Inf)) ≈ norm(Ξ)
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 2, Inf)
-    @test norm(N3, P, Θ) ≈ norm(Ξ)
-    @test sum(abs.(project_collaborative_TV(N3, 0.0, P, Ξ, 1, Inf))) ≈ 1.0
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1, Inf)
-    @test sum(abs.(Θ)) ≈ 1.0
-    @test norm(N3, P, project_collaborative_TV(N3, 0.0, P, Ξ, Inf, Inf)) ≈ norm(Ξ)
-    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, Inf, Inf)
-    @test norm(N3, P, Θ) ≈ norm(Ξ)
-    @test_throws ErrorException project_collaborative_TV(N3, 0.0, P, Ξ, 3, 3)
-    @test_throws ErrorException project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 3, 3)
-    @test_throws ErrorException project_collaborative_TV(N3, 0.0, P, Ξ, 3, 1)
-    @test_throws ErrorException project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 3, 1)
-    @test_throws ErrorException project_collaborative_TV(N3, 0.0, P, Ξ, 3, Inf)
-    @test_throws ErrorException project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 3, Inf)
-
-    @testset "Multivariate project collaborative TV" begin
-        S = Sphere(2)
-        M = PowerManifold(S, NestedPowerRepresentation(), 2, 2, 2)
-        p = [zeros(3) for i in [1, 2], j in [1, 2], k in [1, 2]]
-        p[1, 1, 1] = [1.0, 0.0, 0.0]
-        p[1, 2, 1] = 1 / sqrt(2) .* [1.0, 1.0, 0.0]
-        p[2, 1, 1] = 1 / sqrt(2) .* [1.0, 0.0, 1.0]
-        p[2, 2, 1] = [0.0, 1.0, 0.0]
-        p[:, :, 2] = deepcopy(p[:, :, 1])
-        X = zero_vector(M, p)
-        X[1, 1, 1] .= [0.0, 0.5, 0.5]
-        Y = zero_vector(M, p)
-        @test norm(project_collaborative_TV(M, 1, p, X, 2, 1)) ≈ 0
-        project_collaborative_TV!(M, Y, 1, p, X, 2, 1)
-        @test norm(Y) ≈ 0
-        @test norm(project_collaborative_TV(M, 0.5, p, X, 2, 1)) ≈ (norm(X[1, 1, 1]) - 0.5)
-        project_collaborative_TV!(M, Y, 0.5, p, X, 2, 1)
-        @test norm(Y) ≈ (norm(X[1, 1, 1]) - 0.5)
-        Nf = PowerManifold(S, NestedPowerRepresentation(), 2, 2, 1)
-        @test_throws ErrorException project_collaborative_TV(Nf, 1, p, X, 2, 1)
-        @test_throws ErrorException project_collaborative_TV!(Nf, Y, 1, p, X, 2, 1)
-    end
-end
diff --git a/test/helpers/test_data.jl b/test/helpers/test_data.jl
deleted file mode 100644
index 815de23bdc..0000000000
--- a/test/helpers/test_data.jl
+++ /dev/null
@@ -1,31 +0,0 @@
-@testset "Data" begin
-    @test artificialIn_SAR_image(2) == 2 * π * ones(2, 2)
-
-    @test artificial_S1_slope_signal(20, 0.0) == repeat([-π / 2], 20)
-
-    @test ismissing(artificial_S1_signal(-1.0))
-    @test ismissing(artificial_S1_signal(2.0))
-    @test artificial_S1_signal(2) == [-3 * π / 4, -3 * π / 4]
-
-    # for the remainder check data types only
-    @test length(artificial_S1_signal(20)) == 20
-
-    @test size(artificial_S2_whirl_image(64)) == (64, 64)
-    @test length(artificial_S2_whirl_image(64)[1, 1]) == 3
-
-    @test size(artificial_S2_rotation_image(64)) == (64, 64)
-    @test length(artificial_S2_rotation_image(64)[1, 1]) == 3
-
-    @test size(artificial_S2_whirl_patch(8)) == (8, 8)
-    @test length(artificial_S2_whirl_patch(8)[1, 1]) == 3
-
-    @test size(artificial_SPD_image(8)) == (8, 8)
-    @test size(artificial_SPD_image(8)[1, 1]) == (3, 3)
-
-    @test size(artificial_SPD_image2(8)) == (8, 8)
-    @test size(artificial_SPD_image2(8)[1, 1]) == (3, 3)
-    @test eltype(artificial_SPD_image2(8)) == Array{Float64,2}
-
-    @test length(artificial_S2_lemniscate([0.0, 0.0, 1.0], 20)) == 20
-    @test length(artificial_S2_lemniscate([0.0, 0.0, 1.0], 20)[1]) == 3
-end
diff --git a/test/helpers/test_error_measures.jl b/test/helpers/test_error_measures.jl
deleted file mode 100644
index f891b10993..0000000000
--- a/test/helpers/test_error_measures.jl
+++ /dev/null
@@ -1,17 +0,0 @@
-@testset "Manopt.jl Error Measures" begin
-    M = Sphere(2)
-    N = PowerManifold(M, NestedPowerRepresentation(), 2)
-    using Random: seed!
-    seed!(42)
-    d = Manifolds.uniform_distribution(M, [1.0, 0.0, 0.0])
-    w = rand(d)
-    x = rand(d)
-    y = rand(d)
-    z = rand(d)
-    a = [w, x]
-    b = [y, z]
-    @test meanSquaredError(M, x, y) == distance(M, x, y)^2
-    @test meanSquaredError(N, a, b) == 1 / 2 * (distance(M, w, y)^2 + distance(M, x, z)^2)
-    @test meanAverageError(M, x, y) == distance(M, x, y)
-    @test meanAverageError(N, a, b) == 1 / 2 * sum(distance.(Ref(M), a, b))
-end
diff --git a/test/plans/test_higher_order_primal_dual_plan.jl b/test/plans/test_higher_order_primal_dual_plan.jl
index 702a0bea09..fb15ae381a 100644
--- a/test/plans/test_higher_order_primal_dual_plan.jl
+++ b/test/plans/test_higher_order_primal_dual_plan.jl
@@ -1,5 +1,7 @@
 using Manopt, Manifolds, ManifoldsBase, Test
-
+using ManoptExamples: forward_logs, adjoint_differential_forward_logs
+using ManifoldDiff:
+    differential_shortest_geodesic_startpoint, differential_shortest_geodesic_startpoint!
 @testset "Test higher order primal dual plan" begin
     # Perform an really easy test, just compute a mid point
     #
@@ -134,10 +136,10 @@ using Manopt, Manifolds, ManifoldsBase, Test
     end
 
     function Dprox_F(M, λ, p, X)
-        return Manopt.differential_shortest_geodesic_startpoint(M, p, data, λ / (α + λ), X)
+        return differential_shortest_geodesic_startpoint(M, p, data, λ / (α + λ), X)
     end
     function Dprox_F!(M, Y, λ, p, X)
-        Manopt.differential_shortest_geodesic_startpoint!(M, Y, p, data, λ / (α + λ), X)
+        differential_shortest_geodesic_startpoint!(M, Y, p, data, λ / (α + λ), X)
         return Y
     end
     function Dprox_G_dual(N, n, λ, X, Y)
diff --git a/test/plans/test_primal_dual_plan.jl b/test/plans/test_primal_dual_plan.jl
index 1ebb2b0402..feeaa037c4 100644
--- a/test/plans/test_primal_dual_plan.jl
+++ b/test/plans/test_primal_dual_plan.jl
@@ -1,4 +1,15 @@
-using Manopt, Manifolds, ManifoldsBase, Test
+using Manopt, Manifolds, ManifoldsBase, ManifoldDiff, ManoptExamples, Test
+
+using ManoptExamples:
+    forward_logs,
+    forward_logs!,
+    adjoint_differential_forward_logs,
+    adjoint_differential_forward_logs!,
+    project_collaborative_TV,
+    project_collaborative_TV!,
+    differential_forward_logs,
+    differential_forward_logs!
+using ManifoldDiff: prox_distance, prox_distance!, adjoint_differential_log_argument
 
 include("../utils/dummy_types.jl")
 
diff --git a/test/plans/test_proximal_plan.jl b/test/plans/test_proximal_plan.jl
index 7beb1a3cb9..469cc7af3b 100644
--- a/test/plans/test_proximal_plan.jl
+++ b/test/plans/test_proximal_plan.jl
@@ -1,5 +1,7 @@
 using LRUCache, Manopt, Manifolds, Test
 import Manopt: get_proximal_map, get_proximal_map!
+using ManifoldDiff: prox_distance, adjoint_differential_log_basepoint
+
 function get_proximal_map(M, o::ManifoldProximalMapObjective, λ, p)
     return get_proximal_map(M, o, λ, p, 1)
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3d1f361186..24ef32c535 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -27,18 +27,7 @@ include("utils/example_tasks.jl")
         include("plans/test_stopping_criteria.jl")
         include("plans/test_subgradient_plan.jl")
     end
-    @testset "Function Tests     " begin
-        include("functions/test_adjoint_differentials.jl")
-        include("functions/test_bezier.jl")
-        include("functions/test_differentials.jl")
-        include("functions/test_costs.jl")
-        include("functions/test_gradients.jl")
-        include("functions/test_proximal_maps.jl")
-        include("functions/test_manifold.jl")
-    end
-    @testset "Helper & Data Tests" begin
-        include("helpers/test_error_measures.jl")
-        include("helpers/test_data.jl")
+    @testset "Helper Tests       " begin
         include("helpers/test_checks.jl")
         include("helpers/test_linesearches.jl")
     end
diff --git a/test/solvers/test_ChambollePock.jl b/test/solvers/test_ChambollePock.jl
index bb87372316..81f5ffbea8 100644
--- a/test/solvers/test_ChambollePock.jl
+++ b/test/solvers/test_ChambollePock.jl
@@ -1,4 +1,5 @@
 using Manopt, Manifolds, ManifoldsBase, Test
+using ManoptExamples: forward_logs, adjoint_differential_forward_logs
 
 @testset "Chambolle-Pock" begin
     #
diff --git a/test/solvers/test_Douglas_Rachford.jl b/test/solvers/test_Douglas_Rachford.jl
index dfc388f20d..a2b228f78a 100644
--- a/test/solvers/test_Douglas_Rachford.jl
+++ b/test/solvers/test_Douglas_Rachford.jl
@@ -1,4 +1,6 @@
 using Manifolds, Manopt, Test
+using ManifoldDiff: prox_distance
+
 @testset "DouglasRachford" begin
     # Though this seems a strange way, it is a way to compute the mid point
     M = Sphere(2)
diff --git a/test/solvers/test_cyclic_proximal_point.jl b/test/solvers/test_cyclic_proximal_point.jl
index 08a9351971..0aeec49008 100644
--- a/test/solvers/test_cyclic_proximal_point.jl
+++ b/test/solvers/test_cyclic_proximal_point.jl
@@ -1,11 +1,13 @@
 using Manifolds, Manopt, Test, Dates, LRUCache
+using ManifoldDiff: prox_distance, prox_distance!
+using ManoptExamples: prox_TV, prox_TV!, L2_Total_Variation
 
 @testset "Cyclic Proximal Point" begin
     @testset "Allocating" begin
         n = 100
         N = PowerManifold(Circle(), n)
         q = artificial_S1_signal(n)
-        f(M, p) = costL2TV(M, q, 0.5, p)
+        f(M, p) = L2_Total_Variation(M, q, 0.5, p)
         proxes = (
             (N, λ, p) -> prox_distance(N, λ, q, p), (N, λ, p) -> prox_TV(N, 0.5 * λ, p)
         )
@@ -44,7 +46,7 @@ using Manifolds, Manopt, Test, Dates, LRUCache
         M = Sphere(2)
         N = PowerManifold(M, NestedPowerRepresentation(), n)
         q = artificial_S2_lemniscate([0.0, 0.0, 1.0], n)
-        f(N, p) = costL2TV(N, q, 0.5, p)
+        f(N, p) = L2_Total_Variation(N, q, 0.5, p)
         proxes! = (
             (N, qr, λ, p) -> prox_distance!(N, qr, λ, q, p),
             (N, q, λ, p) -> prox_TV!(N, q, 0.5 * λ, p),
@@ -90,7 +92,7 @@ using Manifolds, Manopt, Test, Dates, LRUCache
         M = Sphere(2)
         N = PowerManifold(M, NestedPowerRepresentation(), n)
         q = artificial_S2_lemniscate([0.0, 0.0, 1.0], n)
-        f(N, x) = costL2TV(N, q, 0.5, x)
+        f(N, x) = L2_Total_Variation(N, q, 0.5, x)
         proxes! = (
             (N, qr, λ, p) -> prox_distance!(N, qr, λ, q, p),
             (N, q, λ, p) -> prox_TV!(N, q, 0.5 * λ, p),
@@ -123,7 +125,7 @@ using Manifolds, Manopt, Test, Dates, LRUCache
         M = Euclidean(3)
         p = ones(3)
         O = CyclicProximalPointState(M, p)
-        f(M, p) = costL2TV(M, q, 0.5, p)
+        f(M, p) = L2_Total_Variation(M, q, 0.5, p)
         proxes = (
             (M, λ, p) -> prox_distance(M, λ, q, p), (M, λ, p) -> prox_TV(M, 0.5 * λ, p)
         )
diff --git a/test/solvers/test_gradient_descent.jl b/test/solvers/test_gradient_descent.jl
index de84700e19..3a0cc64750 100644
--- a/test/solvers/test_gradient_descent.jl
+++ b/test/solvers/test_gradient_descent.jl
@@ -1,5 +1,7 @@
 using Manopt, Manifolds, Test, Random
 
+using ManifoldDiff: grad_distance
+
 @testset "Gradient Descent" begin
     @testset "allocating Circle" begin
         M = Circle()
diff --git a/test/solvers/test_primal_dual_semismooth_Newton.jl b/test/solvers/test_primal_dual_semismooth_Newton.jl
index f0c3583168..d7d08fdbad 100644
--- a/test/solvers/test_primal_dual_semismooth_Newton.jl
+++ b/test/solvers/test_primal_dual_semismooth_Newton.jl
@@ -1,4 +1,6 @@
 using Manopt, Manifolds, ManifoldsBase, Test
+using ManoptExamples: adjoint_differential_forward_logs
+using ManifoldDiff: differential_shortest_geodesic_startpoint
 
 @testset "PD-RSSN" begin
     #