function [hat_r_ij,norm_p,norm_d,k] = distcorr_admm(yi, yj, tau1, tau2, tau, M, rho, k_max, eps_abs, eps_rel)

% DISTCORR_ADMM  Alternating Direction Method of Multipliers (ADMM) algorithm
%                for solving distributed cross-correlation estimation problem
%
% OUTPUT        hat_r_ij    : SLP coefficient vector estimate [n x 1]
%               norm_p      : primal residual norm (all iterations) [k x 1]
%               norm_d      : dual residual norm (all iterations) [k x 1]
%               k           : number of iterations performed
%
% INPUT         yi          : ith sensor signal vector [N x 1]
%               yj          : jth sensor signal vector [N x 1]
%               tau1        : minimum lag for which cross-correlation
%                             function should be estimated
%                             NOTE: it is assumed here that a positive
%                             value for tau1 corresponds to a negative lag
%               tau2        : maximum lag for which cross-correlation
%                             function should be estimated
%               tau         : maximum lag for which sensor signal
%                             autocorrelation is computed in nodes i and j
%               M           : number of auxiliary signal samples
%                             transmitted from node i to j (and vice versa)
%               rho         : ADMM parameter
%               k_max       : maximum number of iterations
%               eps_abs     : absolute error tolerance for ADMM termination
%                             criterion 
%               eps_rel     : relative error tolerance for ADMM termination
%                             criterion 
%
% Author: Toon van Waterschoot
%         KU Leuven
%         Department of Electrical Engineering (ESAT)
%         Stadius Center for Dynamical Systems, Signal Processing and Data
%         Analytics (STADIUS) 
%         Kasteelpark Arenberg 10
%         B-3001 Leuven, Belgium
%         toon.vanwaterschoot@esat.kuleuven.be
%         http://homes.esat.kuleuven.be/~tvanwate
% Key reference: Toon van Waterschoot, "Distributed estimation of cross-
%                correlation functions in ad-hoc microphone arrays," in
%                Proc. 23rd European Signal Process. Conf. (EUSIPCO '15),
%                Nice, France, Sept. 2015.   
% Date: 27/02/2015
% Version: 1-0
% License: Copyright (c) 2015, Toon van Waterschoot. All rights reserved.
%          Published under BSD license.
%          http://homes.esat.kuleuven.be/~tvanwate/bsd_license.html

%% Initialization %%
    %-- initialization of internal ADMM variables --    
    hat_r_ij = zeros(tau1+tau2+1,1); % consensus variable initialization
    hat_lambda_ij_j = zeros(tau1+tau2+1,1); % node j dual variable initialization
    hat_lambda_ij_i = zeros(tau1+tau2+1,1); % node i dual variable initialization
    norm_p = zeros(k_max,1); % primal residual norm (all iterations)
    norm_d = zeros(k_max,1); % dual residual norm (all iterations)
    norm_p_k = Inf; % primal residual norm (iteration k)
    norm_d_k = Inf; % dual residual norm (iteration k)
    eps_p = 0; % primal feasibility tolerance
    eps_d = 0; % dual feasibility tolerance

%% Algorithm %%
    %-- calculation of iteration-independent variables --
        %-- full (local) variables --
        t1 = 1;
        t2 = length(yi);
        if length(yj) ~= t2,
            error('Signals yi and yj should have the same length');
        end
        Yi = toeplitz([yi(t1+tau1:t2);zeros(tau1,1)],[yi(t1+tau1:-1:t1)',zeros(1,tau2)]);
        Yj = toeplitz([yj(t1+tau1:t2);zeros(tau1,1)],[yj(t1+tau1:-1:t1)',zeros(1,tau2)]);
        tildeYi = hankel([zeros(tau1,1);yi(t1:t2-tau1)],[yi(t2-tau1:t2)',zeros(1,tau2)]);
        tildeYj = hankel([zeros(tau1,1);yj(t1:t2-tau1)],[yj(t2-tau1:t2)',zeros(1,tau2)]);
        r_ii = xcorr(yi,tau);
        r_jj = xcorr(yj,tau);
        s_ijj = Yi*r_jj;
        s_iij = tildeYj*r_ii;
    
        %-- truncated (shared) variables --
        s_ijjM = s_ijj(tau1+1:tau1+M);
        s_iijM = s_iij(tau1+1:tau1+M);
        YjM = Yj(tau1+1:tau1+M,:);
        tildeYiM = tildeYi(tau1+1:tau1+M,:);
        sigma_jM = YjM'*s_ijjM;
        sigma_iM = tildeYiM'*s_iijM;
        Phi_jM = YjM'*YjM + rho*eye(tau1+tau2+1);
        tildePhi_iM = tildeYiM'*tildeYiM + rho*eye(tau1+tau2+1);
        LjM = chol(Phi_jM,'lower');
        tildeLiM = chol(tildePhi_iM,'lower');
          
    %-- ADMM algorithm --
    k = 1; % initialization of iteration index
    while ((~((norm_p_k < eps_p)&&(norm_d_k < eps_d))) && k <= k_max)   % check termination criterion
        %-- solve primal problem at node j --
        temp_j = LjM\(sigma_jM + rho*hat_r_ij - hat_lambda_ij_j);       % at node j: primal problem, forward subtitution
        hat_r_ij_j = LjM'\temp_j;                                       % at node j: primal problem, backward substitution

        %-- solve primal problem at node i --
        temp_i = tildeLiM\(sigma_iM + rho*hat_r_ij - hat_lambda_ij_i);  % at node i: primal problem, forward subtitution
        hat_r_ij_i = tildeLiM'\temp_i;                                  % at node i: primal problem, backward substitution

        %-- variable sharing --
        % transmit hat_r_ij_j from node j to i
        % transmit hat_r_ij_i from node i to j
        
        %-- compute consensus variable at node j and i --
        hat_r_ij_old = hat_r_ij;                                        % at node 
        hat_r_ij = (hat_r_ij_j + hat_r_ij_i)./2;                        % at node i and j: consensus variable = average of estimated primal variables
        
        %-- solve dual problem at node j --
        hat_lambda_ij_j = hat_lambda_ij_j + rho*(hat_r_ij_j - hat_r_ij);% at node j: dual problem (dual variables update)
        norm_hat_lambda_ij_j = norm(hat_lambda_ij_j);                   % at node j: compute dual variable vector norm (needed for stopping criterion)
        
        %-- solve dual problem at node i --
        hat_lambda_ij_i = hat_lambda_ij_i + rho*(hat_r_ij_i - hat_r_ij);% at node i: dual problem (dual variables update)
        norm_hat_lambda_ij_i = norm(hat_lambda_ij_i);                   % at node i: compute dual variable vector norm (needed for stopping criterion)
        
        %-- variable sharing --
        % transmit norm_hat_lambda_ij_j from node j to i
        % transmit norm_hat_lambda_ij_i from node i to j
        
        %-- compute residuals and stopping criterion variables at node j and i --
        norm_p_k = norm([hat_r_ij_j - hat_r_ij;hat_r_ij_i - hat_r_ij]); % primal residual norm
        norm_d_k = sqrt(2)*rho*norm(hat_r_ij - hat_r_ij_old);           % dual residual norm
        eps_p = eps_abs + eps_rel*max(norm([hat_r_ij_j;hat_r_ij_i]),norm([hat_r_ij;hat_r_ij])); % primal feasibility tolerance
        eps_d = eps_abs + eps_rel*sqrt(norm_hat_lambda_ij_j^2 + norm_hat_lambda_ij_i^2);        % dual feasibility tolerance

        norm_p(k) = norm_p_k;                                           % save primal residual norm
        norm_d(k) = norm_d_k;                                           % save dual residual norm

        k = k+1;                                                        % increase iteration index
       
    end
    
%% Output variables %%    
norm_p = norm_p(1:k);
norm_d = norm_d(1:k);
