% function [Rwcn, bsum] = wcnoise(Rxx, H, Ly, dual_gap, nerr)
%
% inputs 
% H is U*Ly by Lx, where
%       Ly is the (constant) number of antennas/receiver,
%       Lx is the number of transmit antennas, and
%       U is the number of users. H can be a complex matrix
% Rxx is the Lx by Lx input (nonsingular) autocorrelation matrix
%       which can be complex (and Hermitian!)
% dual_gap is the duality gap, defaulting to 1e-6 in wcnoise
% nerr is Newton's method acceptable error, defaulting to 1e-4 in wcnoise
%
% outputs
% Rwcn is the U*Ly by U*Ly worst-case-noise autocorrelation matrix. 
% bsum is the rate-sum/real-dimension.
%
% I = 0.5 * log(det(H*Rxx*H'+Rwcn)/det(Rwcn)), Rwcn has Ly x Ly diagonal blocks  
%     that are each equal to an identity matrix
%
function [Rwcn, bsum] = wcnoise(Rxx, H, Ly, dual_gap, nerr)

switch nargin
    case 3
        dual_gap = 1e-6;
        nerr = 1e-4;
    case 4
        nerr = 1e-4;
end

[n,m] = size(H);
K = n / Ly;

% make sure the input Rxx is strickly hermitian
Rxx = (real(Rxx)+real(Rxx'))/2 + sqrt(-1) * (imag(Rxx)+imag(Rxx'))/2;
for i = 1:m
  Rxx(i,i) = real(Rxx(i,i));
end

if min(eig(Rxx)) < 0
  fprintf('Input Rxx is not positive semi-definite! The answers may be wrong')
end

count = 1;

A = zeros(n,n,n*n);
for i = 1:n
   A(i,i,count) = 1;
   count = count+1;
end



for i = 2:n
   for j = 1:i - 1
      A(i,j,count) = 1;
      A(j,i,count) = 1;
      count = count+1;
   end
end

complexI = sqrt(-1);

for i = 2:n
   for j = 1:i - 1
      A(i,j,count) = complexI;
      A(j,i,count) = -complexI;
      count = count+1;
   end
end


map = zeros(n,n);
for i = 1:K
   map((i-1) * Ly + 1:i * Ly,(i-1) * Ly + 1:i * Ly) = ones(Ly,Ly);
end




NT_max_it = 1000;                              % Maximum number of Newton's 
                                               % method iterations
%dual_gap = 1e-6;
mu = 10;                                       % step size for t
alpha = 0.001;                                 % back tracking line search parameters
beta = 0.5;

count = 1;
%n%nerr = 1e-4;                                   % acceptable error for inner loop 
                                               % Newton's method

v_0 = zeros(n*n,1);                            % Strictly feasible point;
v_0(1:n) = 0.5 *  ones(n,1);
v = v_0;
t = 1;                                         
l_v = 1;                                       % lambda(v) for newton's method termination



while (1+n*n)/t > dual_gap
   t =  t * mu;
   l_v = 1;
   count = 1;
   while l_v > nerr && count < NT_max_it 
      
      f_val = 0;                               % calculating function value
      Rwcn = zeros(n,n);
      Rzprime = zeros(n,n);
      
      for i = 1:n*n                            % computing Rz
         Rwcn = Rwcn + v(i) * A(:,:,i);
      end
      
      for i = 1:K
          Rzprime((i-1) * Ly + 1:i * Ly,(i-1) * Ly + 1:i * Ly) = Rwcn((i-1) * Ly + 1:i * Ly,(i-1) * Ly + 1:i * Ly);
      end
      % the real operation is not needed below as determinant of a Hermitian
      % matrix is always real, it is included for avoiding numerical issues
      f_val = t * log(real(det(H * Rxx * H' + Rwcn))) - (t + 1) * log(real(det(Rwcn))) - log(real(det(eye(n) - Rzprime)));
      
      S = inv(H * Rxx * H' + Rwcn);
      Q = inv(eye(n) - Rzprime);
      Rz_inv = inv(Rwcn);
      
      g = zeros(n*n,1);
      h = zeros(n*n,n*n);
      
      for i = 1:n*n
         g(i) = t * trace(A(:,:,i) * S) - (t + 1) * trace(A(:,:,i) * Rz_inv)...
              + (sum(sum(A(:,:,i) .* map)) ~= 0) * trace(A(:,:,i) * Q) ;  % gradient
      end
      
      % make sure g is always real (gradient of a real function with respect to real variables)
      g = real(g);
      
      for i = 1:n*n
         for j = 1:n*n
            h(i,j) = -t * trace(A(:,:,i) * S * A(:,:,j) * S) + (t + 1) * trace(A(:,:,i) * Rz_inv * A(:,:,j) * Rz_inv)...
                     + (sum(sum(A(:,:,i) .* map)) ~= 0) * (sum(sum(A(:,:,j) .* map)) ~= 0) * trace(A(:,:,i) * Q * A(:,:,j) * Q); %hessian
         end
      end
      
      % make sure h is always real (gradient of a real gradient values)
      % make sure he is symmetric
      h = (h+h')/2;
      h = real(h);
      
      
         
      dv = real(-h\g);                         % search direction
      
      v_min = min(abs(dv));
      
      s = 1;                                   % checking v = v+s*dx feasible 
                                               % and also back tracking algorithm
      
      while 1                
         v_new = v + s * dv;
      
         Rz_new = zeros(n,n);
      
         for i = 1:n*n
            Rz_new = Rz_new + v_new(i) * A(:,:,i);
         end
      
         for i = 1:K
             Rzprime((i-1) * Ly + 1:i * Ly,(i-1) * Ly + 1:i * Ly) = Rz_new((i-1) * Ly + 1:i * Ly,(i-1) * Ly + 1:i * Ly);
         end
      
         f_new = t * log(real(det(H * Rxx * H' + Rz_new))) - (t + 1) * log(real(det(Rz_new))) - log(real(det(eye(n) - Rzprime)));
      

         feas_check = 1;
         if real(eig(Rz_new)) > zeros(n,1)
             feas_check = 1;
         else
             feas_check = 0;
         end
         if real(eig(eye(n) - Rzprime)) > zeros(n,1)
             feas_check = 1;
         else
             feas_check = 0;
         end
         
         if feas_check == 1
           feas_check = feas_check * (f_new < f_val + alpha * s * g' * dv);
         end
     
         if feas_check == 1
           break
         end
          
         s = s * beta;
         
         if s < 1e-30 * v_min
           % s is too small, break the while loop
           s = 0;
           break;
         end

      end   
      
      v = v + s * dv;                         % update v
      l_v = -g'*dv ;                          % lambda(v)^2 for Newton's method
      
      if s > 0
        % valid step size returned by the line search algorithm above
        count = count + 1;                    % number of Newtons method iterations
      else
        % s = 0 and we could not find a feasible newton step, skip this iteration and go to a larger t
        count = NT_max_it;
      end
   end
   if l_v >= nerr
     l_v
     fprintf('Inner loop Newton method failed to converge for ')
     t
   end
  
end

Rwcn = zeros(n,n);
      
for i = 1:n*n
   Rwcn = Rwcn + v(i) * A(:,:,i);
end

for i = 1:K
    Rwcn((i-1) * Ly + 1:i * Ly,(i-1) * Ly + 1:i * Ly) = eye(Ly);
end
bsum = 0.5 * log2(real(det(H * Rxx * H' + Rwcn)/det(Rwcn)));