function [x_star, z_star, u_star] = ADMM_convex_solver(samples, y, connection_matrix, lambda, c, x_init, z_init, u_init, nbSamplePerNode)

n = length(connection_matrix);

epsPrimal = 0.1;
epsDual = 0.1;
primalResidual{1} = 1;
dualResidual{1} = 1;
rho = 1; % "rho > 0"

%keep track of the number of iterations, in case we never converge:
curr = 0;
maxTimeout = 50;

%initialize based on lambda=0 params:
x_previous = x_init;
z_previous = z_init;
u_previous = u_init;

%make list of x and u params for all nodes connected to each x:
x_init = cell2mat(x_init);
u_init = cell2mat(u_init);
for i=1:n
    x_edge{i} = x_init(:, connection_matrix(:,i) ~= 0);
    u_edge{i} = u_init(:, connection_matrix(:,i) ~= 0);
end

clc;
lambda


while (curr < maxTimeout && (max(cell2mat(primalResidual)) > epsPrimal || max(cell2mat(dualResidual)) > epsDual))
    
    % x-step
    for i=1:n
        start = (i-1)*nbSamplePerNode + 1;
        stop = start + nbSamplePerNode - 1;
        x_new{i} = x_step(samples(start:stop, :), y(start:stop), rho, z_previous{i}, u_previous{i}, c);
    end
    
    % z-step + u-step
    x_new_mat = cell2mat(x_new);    
    for i=1:n
        x_edge{i} = x_new_mat(:, connection_matrix(:,i) ~= 0);
    end
    
    for i=1:n
       N_i = size(x_edge{i}, 2);
       w = ones(1, N_i);
       z_new{i} = z_step(x_new{i}, u_previous{i}, w', rho, x_edge{i}, u_edge{i}, lambda);
       u_new{i} = u_step(x_new{i}, u_previous{i}, z_new{i}); 
       
       x_rep = repmat(x_new{i}, [1, N_i]);
        primalResidual{i} = sqrt( max( diag( (x_rep-z_new{i})'*(x_rep-z_new{i}) ) ) );
        dualResidual{i} = norm(x_new{i} - x_previous{i}, 2); 
    end
    x_previous = x_new;
    z_previous = z_new;
    u_previous = u_new;
    
    u_new_mat = cell2mat(u_new);
    for i=1:n
        u_edge{i} = u_new_mat(:, connection_matrix(:,i) ~= 0);
    end
    
    curr = curr + 1
    primalRes = max(cell2mat(primalResidual))
    dualRes = max(cell2mat(dualResidual))
end
x_star = x_new;
z_star = z_new;
u_star = u_new;

end