diff -r 7c4dfeb48dd9 -r 7d3a75b875fa src/AlgorithmZeroDual.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/AlgorithmZeroDual.jl	Sun Apr 21 19:04:00 2024 +0300
@@ -0,0 +1,125 @@
+####################################################################
+# Predictive online PDPS for optical flow with known velocity field
+####################################################################
+
+__precompile__()
+
+module AlgorithmZeroDual
+
+identifier = "pdps_known_zerodual"
+
+using Printf
+
+using AlgTools.Util
+import AlgTools.Iterate
+using ImageTools.Gradient
+
+using ..OpticalFlow: ImageSize,
+                     Image,
+                     pdflow!
+
+#########################
+# Iterate initialisation
+#########################
+
+function init_rest(x::Image)
+    imdim=size(x)
+
+    y = zeros(2, imdim...)
+    Δx = copy(x)
+    Δy = copy(y)
+    x̄ = copy(x)
+
+    return x, y, Δx, Δy, x̄
+end
+
+function init_iterates(xinit::Image)
+    return init_rest(copy(xinit))
+end
+
+function init_iterates(dim::ImageSize)
+    return init_rest(zeros(dim...))
+end
+
+############
+# Algorithm
+############
+
+function solve( :: Type{DisplacementT};
+               dim :: ImageSize,
+               iterate = AlgTools.simple_iterate,
+               params::NamedTuple) where DisplacementT
+
+    ################################                                        
+    # Extract and set up parameters
+    ################################                    
+
+    α, ρ = params.α, params.ρ
+    R_K² = ∇₂_norm₂₂_est²
+    γ = 1.0
+    Λ = params.Λ
+    τ₀, σ₀ = params.τ₀, params.σ₀
+
+    τ = τ₀/γ
+    @assert(1+γ*τ ≥ Λ)
+    σ = σ₀*1/(τ*R_K²)
+
+    println("Step length parameters: τ=$(τ), σ=$(σ)")
+
+    ######################
+    # Initialise iterates
+    ######################
+
+    x, y, Δx, Δy, x̄ = init_iterates(dim)
+    init_data = (params.init == :data)
+
+    ####################
+    # Run the algorithm
+    ####################
+
+    v = iterate(params) do verbose :: Function,
+                           b :: Image,
+                           v_known :: DisplacementT,
+                           🚫unused_b_next :: Image
+
+        ##################
+        # Prediction step
+        ##################
+        if init_data
+            x .= b
+            init_data = false
+        end
+
+        pdflow!(x, Δx, y, Δy, v_known, false)
+        y .= zeros(size(y)...)
+
+        ############
+        # PDPS step
+        ############
+
+        ∇₂ᵀ!(Δx, y)                    # primal step:
+        @. x̄ = x                       # |  save old x for over-relax
+        @. x = (x-τ*(Δx-b))/(1+τ)      # |  prox
+        @. x̄ = 2x - x̄                  # over-relax
+        ∇₂!(Δy, x̄)                     # dual step: y
+        @. y = (y + σ*Δy)/(1 + σ*ρ/α)  # |
+        proj_norm₂₁ball!(y, α)         # |  prox
+
+        ################################
+        # Give function value if needed
+        ################################
+        v = verbose() do            
+            ∇₂!(Δy, x)
+            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
+            value, x, [NaN, NaN], nothing
+        end
+
+        v
+    end
+
+    return x, y, v
+end
+
+end # Module
+
+