% the channel matrix is K*N by m, m is the number of transmitter's antennas,
% K is the number of users and N is the number od receiver's antennas for
% each users
% Rx is the m by m input covariance matrix 
% Rz is the K*N by K*N worst case noise minimizing           
% 0.5 * log(det(H*Rx*H'+Rz)/det(Rz)), Rz has its N by N diganoal sub-blocks  
% equal to identity matrix

function [Rz, sumRate] = wcnoise(Rx, H, N)
 
[n,m] = size(H);
K = n / N;



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

count = n+1;

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

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




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

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

v_0 = zeros(n*(n+1)/2,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)/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
      Rz = zeros(n,n);
      Rzprime = zeros(n,n);
      
      for i = 1:n*(n+1)/2                      % computing Rz
         Rz = Rz + v(i) * A(:,:,i);
      end
      
      for i = 1:K
          Rzprime((i-1) * N + 1:i * N,(i-1) * N + 1:i * N) = Rz((i-1) * N + 1:i * N,(i-1) * N + 1:i * N);
      end
      
      f_val = t * log(det(H * Rx * H' + Rz)) - (t + 1) * log(det(Rz)) - log(det(eye(n) - Rzprime));
      
      S = inv(H * Rx * H' + Rz);
      Q = inv(eye(n) - Rzprime);
      Rz_inv = inv(Rz);
      g = zeros(n*(n+1)/2,1);
      h = zeros(n*(n+1)/2,n*(n+1)/2);
      
      for i = 1:n*(n+1)/2
         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
      
      
      for i = 1:n*(n+1)/2
         for j = 1:n*(n+1)/2
            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
      
         
      dv = -h\g;                               % search direction
      
      
      s = 1;                                   % checking v = v+s*dx feasible 
                                               % and also back tracking algorithm
   
      v_new = v + s * dv;
      
      f_new = 0;
      
      Rz_new = zeros(n,n);
      
      for i = 1:n*(n+1)/2
         Rz_new = Rz_new + v_new(i) * A(:,:,i);
      end
      
      for i = 1:K
          Rzprime((i-1) * N + 1:i * N,(i-1) * N + 1:i * N) = Rz_new((i-1) * N + 1:i * N,(i-1) * N + 1:i * N);
      end
      
      f_new = t * log(det(H * Rx * H' + Rz_new)) - (t + 1) * log(det(Rz_new)) - log(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
     
      feas_check = feas_check * (f_new < f_val + alpha * s * g' * dv);

   
      while feas_check ~= 1                
         s = s * beta ;
         
         v_new = v + s * dv;
         f_new = 0;
      
         Rz_new = zeros(n,n);
      
         for i = 1:n*(n+1)/2
            Rz_new = Rz_new + v_new(i) * A(:,:,i);
         end
      
         for i = 1:K
             Rzprime((i-1) * N + 1:i * N,(i-1) * N + 1:i * N) = Rz_new((i-1) * N + 1:i * N,(i-1) * N + 1:i * N);
         end
      
         f_new = t * log(det(H * Rx * H' + Rz_new)) - (t + 1) * log(det(Rz_new)) - log(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
     
         
         feas_check = feas_check * (f_new < f_val + alpha * s * g' * dv);



      end   
      
      v = v + s * dv;                         % update v
      l_v = -g'*dv ;                           % lambda(v)^2 for Newton's method
      count = count + 1;                      % number of Newtons method iterations
   end
  
end

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

for i = 1:K
    Rz((i-1) * N + 1:i * N,(i-1) * N + 1:i * N) = eye(N);
end
sumRate = 0.5 * log2(det(H * Rx * H' + Rz)/det(Rz));