function [x,R,y] = dr(x,ProxF,ProxG,options)
% dr - Douglas-Rachford algorithm
%
%   [x,R] = dr(x,ProxF,ProxG,options);
%
%   Solves
%       min_x F(x)+G(x)
%   where F and G are proper closed convex and simple functions.
%
%   INPUTS:
%   ProxF(y,gamma) computes Prox_{gamma*F}(x)
%   ProxG(y,gamma) computes Prox_{gamma*G}(x)
%   options.niter is the number of iterations.
%   options.gamma the stepsize parameter for DR.
%   options.tau is the relaxation parameter for DR.
%   options.verb is for the diaplay of iterations.
%   options.report(x) is a function to fill in R.
%
%   OUTPUTS:
%   x is the final solution.
%   R(i) = options.report(x) at iteration i.
%

rProxF = @(x,ga)2*ProxF(x,ga)-x;
rProxG = @(x,ga)2*ProxG(x,ga)-x;

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

R = [];
y = x;
for i=1:niter 
    % record energies
    if verb
        progressbar(i,niter);
    end
    yold = y;
    x = ProxG(y,gamma);
    y = (1-tau/2)*y + tau/2*rProxF( (2*x-y),gamma );	   
    R(i) = report(x);
end