# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1577288552 -7200
# Node ID d1f40f6654cbc8dd67fff173aca85be5d06f1c8c
# Parent  51bb83c404a4a7ff8677b04e7b1da25dc7d0bf24
denoise_fista

diff -r 51bb83c404a4 -r d1f40f6654cb src/Denoise.jl
--- a/src/Denoise.jl	Mon Dec 23 21:31:48 2019 +0200
+++ b/src/Denoise.jl	Wed Dec 25 17:42:32 2019 +0200
@@ -14,7 +14,8 @@
 # Our exports
 ##############
 
-export denoise_pdps
+export denoise_pdps,
+       denoise_fista
 
 #############
 # Data types
@@ -40,14 +41,15 @@
     return x, y, Δx, Δy, x̄
 end
 
-function init_iterates(xinit::Image, b)
-    return init_rest(copy(xinit))
+function init_primal(xinit::Image, b)
+    return copy(xinit)
 end
 
-function init_iterates(xinit::Nothing, b :: Image)
-    return init_rest(zeros(size(b)...))
+function init_primal(xinit::Nothing, b :: Image)
+    return zeros(size(b)...)
 end
 
+
 ############
 # Algorithm
 ############
@@ -55,7 +57,7 @@
 function denoise_pdps(b :: Image;
                       xinit :: Union{Image,Nothing} = nothing,
                       iterate = AlgTools.simple_iterate,
-                      params::NamedTuple) where DisplacementT
+                      params::NamedTuple)
 
     ################################                                        
     # Extract and set up parameters
@@ -75,7 +77,7 @@
     # Initialise iterates
     ######################
 
-    x, y, Δx, Δy, x̄ = init_iterates(xinit, b)
+    x, y, Δx, Δy, x̄ = init_rest(init_primal(xinit, b))
 
     ####################
     # Run the algorithm
@@ -111,6 +113,68 @@
     return x, y, v
 end
 
+function denoise_fista(b :: Image;
+                       xinit :: Union{Image,Nothing} = nothing,
+                       iterate = AlgTools.simple_iterate,
+                       params::NamedTuple)
+
+    ################################                                        
+    # Extract and set up parameters
+    ################################                    
+
+    α, ρ = params.α, params.ρ
+    τ₀ =  params.τ₀
+    τ = τ₀/∇₂_norm₂₂_est²
+    
+    ######################
+    # Initialise iterates
+    ######################
+
+    x = init_primal(xinit, b)
+    imdim = size(x)
+    Δx = similar(x)
+    y = zeros(2, imdim...)
+    ỹ = copy(y)
+    y⁻ = similar(y)
+    Δy = similar(y)
+
+    ####################
+    # Run the algorithm
+    ####################
+
+    t = 0
+
+    v = iterate(params) do verbose :: Function                    
+        ∇₂ᵀ!(Δx, ỹ)
+        @. Δx .-= b
+        ∇₂!(Δy, Δx)
+        @. y⁻ = y
+        @. y = (ỹ - τ*Δy)/(1 + τ*ρ/α)
+        proj_norm₂₁ball!(y, α)
+        t⁺ = (1+√(1+4*t^2))/2
+        @. ỹ = y+((t-1)/t⁺)*(y-y⁻)
+        t = t⁺
+
+        ################################
+        # Give function value if needed
+        ################################
+        v = verbose() do
+            ∇₂ᵀ!(Δx, y)
+            @. x = b - Δx
+            ∇₂!(Δy, x)
+            value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ)
+            value, x
+        end
+
+        v
+    end
+
+    ∇₂ᵀ!(Δx, y)
+    @. x = b - Δx
+
+    return x, y, v
+end
+
 end # Module
 
 
diff -r 51bb83c404a4 -r d1f40f6654cb test/denoise.jl
--- a/test/denoise.jl	Mon Dec 23 21:31:48 2019 +0200
+++ b/test/denoise.jl	Wed Dec 25 17:42:32 2019 +0200
@@ -55,7 +55,7 @@
     st, iterate = initialise_visualisation(visualise)
 
     # Run algorithm
-    x, y, st = denoise_pdps(b_noisy; iterate=iterate, params=params)
+    x, y, st = denoise_fista(b_noisy; iterate=iterate, params=params)
 
     if params.save_results
         perffile = params.save_prefix * ".txt"