function [x, R, z] = gfb( z, gradF, proxGi, options)
% gfb - Generalized Forward-Backward algorithm
%
%	 [x, R, z] = gfb(z, gradF, proxGi, options )
%
% Solve
%		min_x F(x)+\sum_{i=1}^{n} Gi(x)
% where
%	F is C^{1,1} and its gradient is 1/be-Lipschitz continuous,
%	G_i are proper closed convex and simple functions, i.e.
%	prox_{ga*Gi}(x) = argmin_y 1/2||x-y||^2 + ga*Gi(x) is easy to compute.
%
% INPUT:
%	'gradF': function handle: class(x) -> class(x)
%		computes the gradient of F
%		can be set to empty-array if F = 0, resulting in a relaxed version of Douglas-Rachford algorithm
%	'proxGi': n-long cell array of function handles : (class(x),double) -> class(x)
%		proxGi{i}(x,mu) computes prox_{mu*Gi}(x)
%	'z': [size(x)]-by-n matrix
%		the initial auxiliary variables
%	options.niter: non-negative integer
%		the number of iterations recquired
%	options.mu: double
%		the descent step-size of the forward step,
%		mu \in ]0,2 be[
%	options.tau [default=1]: double
%		relaxation of the iterations, tau \in ]0,(4-mu be)/2[
%	options.verbose [default=false]: logical
%		set the diaplay of iterations
%	options.report [default=None]: function handle : class(x) -> 1-by-m matrix
%		a user-defined report called at each iteration
% OUTPUT:
%	'x': vector or N-D matrix
%		the final minimizer
%	'R' [default='empty']: niter-by-m matrix
%		the report sequence
%	'z': [size(x)]-by-n matrix
%		the final auxiliary variables
%

method = getoptions(options, 'method', 'fb');
report = getoptions(options, 'report', @(x)0);
niter = getoptions(options, 'niter', 100);
verbose = getoptions(options, 'verbose', 0);
tau = getoptions(options, 'tau', 1);
mu = getoptions(options, 'mu', 1.8);

n = length( proxGi );
catDim = ndims( z );
x = mean( z, catDim );
N = numel( x );
zi = zeros( size( x ) );
R = [];

for it=1:niter
	if verbose, progressbar(it,niter); end
	if isempty( gradF )
		forward = 0;
	else
		forward = mu*gradF(x);
	end
	R(it,:) = report( x );
	for i=1:n
		idx = (i-1)*N+1:i*N;
		zi(:) = z(idx);	
		z(idx) = zi + tau*( proxGi{i}( 2*x - zi - forward, n*mu ) - x );
	end
	x = mean( z, catDim );
end

end