function [x,R] = admmhigher(x, Ki, KSi, IK, ProxFi, ProxG, options)
% admm - ADMM for more than two functions
%
%    [x,R] = admmhigher(x, Ki, IK, ProxFi, ProxG, options)
%
%   Solves 
%       min_x sum_i F_i(K_i*x) + G(x)
%   where F_i and G are proper closed convex and simple functions, 
%   and K_i are linear operators.
%   IK = (Id + \sum_i K_i^*K_i)^{-1} should be easily computable.
%   
%
%   INPUTS:
%   ProxFi(y,ga) computes Prox_{ga*Fi}(y)
%   ProxG(x,ga) computes Prox_{ga*G}(x)
%   Ki is the handle to the linear operator.
%   KSi is the handle to the adjoint linear operator.
%   IK is the handle to the inverse (Id + \sum_i K_i^*K_i)^{-1}.
%   options.ga > 0 is the ADMM parameter.
%   options.verb=0 suppress display of progression.
%   options.niter is the number 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.
%

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

ProxF  	= @(u,gamma)ProxMult(u,gamma,ProxFi);
ProxKG 	= @(u)IK(KSMult(u,KSi));
K 	= @(x)KMult(x,Ki);

[x,R] = admm(x, K, ProxF, ProxKG, options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u = ProxMult(u, gamma, ProxFi, ProxG)

m = length(ProxFi);
for i=1:m
    u(:,i) = ProxFi{i}(u(:,i),gamma);
end
u(:,m+1) = ProxG(u(:,end),gamma);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u = KMult(x, Ki)

m = length(Ki);
for i=1:m
    u(:,i) = Ki{i}(x);
end
u(:,end) = x;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = KSMult(u, KSi)

m = length(KSi);
x = u(:,m+1);
for i=1:m
    x = x + KSi{i}(u(:,i));
end