function [x,R,u] = primaldual(x, K,  KS, ProxFS, ProxG, GradH, options)
% primaldual - primal-dual splitting algorithm
%
%    [x,R] = primaldual(x, K,  KS, ProxFS, ProxG, GradH, options);
%
%   Solves 
%       min_x H(x) + G(x) + F(K*x)
%   where F and G are proper closed convex and simple functions, 
%   H is differentiable with 1/be-Lipschitz gradient
%   and K is a linear operator.
%
%   Uses the extended relaxed Arrow-Hurwicz method proposed by Vu and others.
%   For H(x)=0, we recover the Chambolle and Pock scheme.
%   For F(x)=0, we recover the FB scheme.	
%
%   INPUTS:
%   ProxFS(y,nu) computes Prox_{nu*F^*}(y)
%   ProxG(x,mu) computes Prox_{mu*G}(x)
%   K is the handle to the linear operator.
%   KS is the handle to the adjoint linear operator.
%   options.mu and options.mu are the parameters of the
%       method, they should satisfy 1-mu*nu*norm(K)^2 > mu*be/2
%   options.theta=1 is the extrapolation step, but can be set in [0,1].
%   options.verb=0 suppress display of progression.
%   options.niter is the number of iterations.
%   options.tau is the relaxation parameter.
%   options.report(x) is a function to fill in R.
%
%   OUTPUTS:
%   x is the final solution.
%   R(i) = options.report(x) at iteration i.
%

report = getoptions(options, 'report', @(x,u)0);
tau = getoptions(options, 'tau', 1);
niter = getoptions(options, 'niter', 100);
theta = getoptions(options, 'theta', 1);
verb = getoptions(options, 'verb', 1);


if isnumeric(K)
    K = @(x)K*x;
end  
if isnumeric(KS)
    KS = @(x)KS*x;
end      

%%%% PDS parameters %%%%
nu = getoptions(options, 'nu', -1);
mu = getoptions(options, 'mu', -1);
if mu<0 || nu<0
    [L,e] = compute_operator_norm(@(x)KS(K(x)),randn(size(x)));
    nu = 10;
    mu = .9/(nu*L+be/2); % This choice of (mu,nu) is rather heuristic, and it is essential to adjust these parameters in practice.
end

u = K(x);
R = [];
xold = 0;
uold = 0;
for i=1:niter    
    if verb
        progressbar(i,niter);
    end
    % record energies
    xold = xold + x;
    uold = uold + u;
    R(i,:) = report(xold/i,uold/i);
    % update
    p = ProxG(  x - mu*(KS(u) + GradH(x)), mu);
    q = ProxFS( u + nu*K(p + theta * (p-x)), nu);
    x = x + tau*(p - x);
    u = u + tau*(q - u);
end