diff -r 000000000000 -r a55e35d20336 src/AlgorithmBothMulti.jl.orig
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/AlgorithmBothMulti.jl.orig	Tue Apr 07 14:19:48 2020 -0500
@@ -0,0 +1,277 @@
+######################################################################
+# Predictive online PDPS for optical flow with unknown velocity field
+######################################################################
+
+__precompile__()
+
+module AlgorithmBothMulti
+
+using Printf
+
+using AlgTools.Util
+import AlgTools.Iterate
+using ImageTools.Gradient
+using ImageTools.ImFilter
+
+using ..OpticalFlow: Image,
+                     ImageSize,
+                     DisplacementConstant,
+                     pdflow!,
+                     horn_schunck_reg_prox!,
+                     pointwise_gradiprod_2d!,
+                     horn_schunck_reg_prox_op!,
+                     mldivide_step_plus_sym2x2!,
+                     ConstantDisplacementHornSchunckData,
+                     filter_hs
+
+using ..Algorithm: step_lengths
+
+#########################
+# Iterate initialisation
+#########################
+
+function init_displ(xinit::Image, ::Type{DisplacementConstant}, n::Integer)
+    return xinit, zeros(n, 2)
+end
+
+# function init_displ(xinit::Image, ::Type{DisplacementFull}, n::Integer)
+#     return xinit, zeros(n, 2, size(xinit)...)
+# end
+
+function init_rest(x::Image, u::DisplacementT) where DisplacementT
+    imdim=size(x)
+
+    y = zeros(2, imdim...)
+    Δx = copy(x)
+    Δy = copy(y)
+    x̄ = copy(x)
+
+    return x, y, Δx, Δy, x̄, u
+end
+
+function init_iterates( :: Type{DisplacementT},
+                       xinit::Image,
+                       n::Integer) where DisplacementT    
+    return init_rest(init_displ(copy(xinit), DisplacementT, n)...)
+end
+
+function init_iterates( :: Type{DisplacementT},
+                       dim::ImageSize,
+                       n::Integer) where DisplacementT
+    return init_rest(init_displ(zeros(dim...), DisplacementT, n)...)
+end
+
+############
+# Algorithm
+############
+
+function solve( :: Type{DisplacementT};
+               dim :: ImageSize,
+               iterate = AlgTools.simple_iterate, 
+               params::NamedTuple) where DisplacementT
+
+    ######################
+    # Initialise iterates
+    ######################
+
+    n = params.displacement_count
+    k = 0 # number of displacements we have already
+
+    x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim, n)
+    init_data = (params.init == :data)
+    hs = [ConstantDisplacementHornSchunckData() for i=1:n]
+    #hs = Array{ConstantDisplacementHornSchunckData}(undef, n)
+    A = Array{Float64,3}(undef, n, 2, 2)
+    d = Array{Float64,2}(undef, n, 2)
+
+    # … for tracking cumulative movement
+    ucumulbase = [0.0, 0.0]
+    
+    #############################################
+    # Extract parameters and set up step lengths
+    #############################################
+
+    α, ρ, λ, θ = params.α, params.ρ, params.λ, params.θ
+    R_K² = ∇₂_norm₂₂_est²
+    γ = 1
+    τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²)
+
+    kernel = params.kernel
+    T = params.timestep 
+
+    ####################
+    # Run the algorithm
+    ####################
+
+    b_next_filt = nothing
+    diffu = similar(u[1, :])
+
+    v = iterate(params) do verbose :: Function,
+                           b :: Image,
+                           🚫unused_v_known :: DisplacementT,
+                           b_next :: Image
+
+        ####################################
+        # Smooth data for Horn–Schunck term
+        ####################################
+
+        b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel)
+
+        ################################################
+        # Prediction step
+        ################################################
+
+        # Predict x and y
+        if k==0
+            if init_data
+                x .= b
+                init_data = false
+            end 
+        else
+            # Displacement from previous to this image is estimated as
+            # the difference of their displacements from beginning of window.
+            if k>1
+                @. @views diffu = u[k, :] - u[k-1, :]
+            else
+                @. @views diffu = u[k, :]
+            end
+
+            pdflow!(x, Δx, y, Δy, diffu, params.dual_flow)
+        end
+
+        # Shift stored prox matrices
+        if k==n
+            tmp = copy(u[1, :])
+            ucumulbase .+= tmp
+            for j=1:(n-1)
+                @. @views u[j, :] = u[j+1, :] - tmp
+                hs[j] = hs[j+1]
+            end
+            # Create new struct as original contains references to objects that
+            # have been moved to index n-1.
+            hs[n]=ConstantDisplacementHornSchunckData()
+        else
+            k += 1
+        end
+
+        # Predict u: zero displacement from current to next image, i.e.,
+        # same displacement to beginning of window.
+        if k==1        
+            @. @views u[k, :] = 0.0
+        else
+            @. @views u[k, :] = u[k-1, :]
+        end
+
+        # Predictor proximals tep
+        if params.prox_predict
+            ∇₂!(Δy, x) 
+            @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α))
+            proj_norm₂₁ball!(y, α) 
+        end
+
+        #################################################################################
+        # PDPS step
+        #
+        # For the displacements, with τ̃=τ/k, we need to solve for 2≤j<k,
+        #
+        # (1) (I/τ̃+M₀^j+M₀^{j+1})u^j = M₀^ju^{j-1} + M₀^{j+1}u^{j+1}
+        #                              + ũ^j/τ̃ - z^j + z^{j+1}, 
+        #
+        # as well as
+        #
+        # (2) (I/τ̃+M₀^k)u^k = M₀^k u^{k-1} + ũ^k/τ̃ - z^k
+        #
+        # and
+        #
+        # (3) (I/τ̃+M₀^1+M₀^2)u^1 = 0 + M₀^{2}u^{2} + ũ^1/τ̃ - z^1 + z^{2}
+        #
+        # We first construct from (2) that
+        #
+        #     u^k = A^k u^{k-1} + d^k
+        #
+        # for
+        #
+        #     A^k := (I/τ̃+M₀^k)^{-1} M₀^k
+        #     d_k := (I/τ̃+M₀^k)^{-1} (ũ^k/τ̃ - z^k).
+        #
+        # Inserting this into (1) we need
+        #
+        # (4)  (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))u^j = M₀^ju^{j-1} + M₀^{j+1}d^{j+1}
+        #                                          + ũ^j/τ̃ - z^j + z^{j+1}.
+        # 
+        # This is well-defined because A^{j+1} < I. It also has the same form as (1), so 
+        # we continue with
+        #
+        # (5)   u^j = A^j u^{j-1} + d^j
+        #
+        # for
+        #
+        #      A^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} M₀^j
+        #      d^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1}
+        #               (M₀^{j+1}d^{j+1} + ũ^j/τ̃ - z^j + z^{j+1})
+        # 
+        # Finally from (3) with these we need
+        #
+        #      (I/τ̃+M₀^1+M₀^2(I-A^2))u^1 = M₀^2d^2 + ũ^1/τ̃ - z^1 + z^2,
+        #
+        # which is of the same form as (4) with u^0=0, so by (5) u^1=d^1.
+        #
+        #################################################################################
+
+        ∇₂ᵀ!(Δx, y)                         # primal step:
+        @. x̄ = x                            # |  save old x for over-relax
+        @. x = (x-τ*(Δx-b))/(1+τ)           # |  prox
+                                            # |  | for displacement
+        # Calculate matrices for latest data; rest is stored. 
+        @views begin
+            horn_schunck_reg_prox_op!(hs[k], b_next_filt, b_filt, θ, λ, T)
+
+            τ̃=τ/k
+
+            B = hs[k].M₀
+            c = u[k, :]./τ̃-hs[k].z
+            mldivide_step_plus_sym2x2!(A[k, :, :], B, B, τ̃)
+            mldivide_step_plus_sym2x2!(d[k, :], B, c, τ̃)
+
+            for j=(k-1):-1:1
+                B = hs[j].M₀+hs[j+1].M₀*([1 0; 0 1]-A[j+1, :, :])
+                c = hs[j+1].M₀*d[j+1, :]+u[j, :]./τ̃-hs[j].z+hs[j+1].z
+                mldivide_step_plus_sym2x2!(A[j, :, :], B, hs[j].M₀, τ̃)
+                mldivide_step_plus_sym2x2!(d[j, :], B, c, τ̃)
+            end
+
+            u[1, :] .= d[1, :]
+            for j=2:k
+                u[j, :] .= A[j, :, :]*u[j-1, :] + d[j, :]
+            end
+        end
+
+        @. x̄ = 2x - x̄                       # over-relax: x̄ = 2x-x_old
+        ∇₂!(Δy, x̄)                          # dual step: y
+        @. y = (y + σ*Δy)/(1 + σ*ρ/α)       # |
+        proj_norm₂₁ball!(y, α)              # |  prox
+
+        ########################################################
+        # Give function value and cumulative movement if needed
+        ########################################################
+        v = verbose() do
+            ∇₂!(Δy, x)
+            hs_plus_reg=0            
+            for j=1:k
+                v=(j==1 ? u[j, :] : u[j, :]-u[j-1, :])                
+                hs_plus_reg += hs[j].cv/2 + dot(hs[j].Mv*v, v)/2+dot(hs[j].av, v)
+            end
+            value = (norm₂²(b-x)/2 + hs_plus_reg/k + α*γnorm₂₁(Δy, ρ))
+
+            value, x, u[k, :]+ucumulbase, u[1:k,:].+ucumulbase'
+        end
+
+        return v
+    end
+
+    return x, y, v
+end
+
+end # Module
+
+