function [Rates,Q,index] = f_wtdratesumcap2(weights,H,SNRdB,NUM_STEPS,t_tol)
% F_WTDRATESUMCAP2  Find the optimal rate vector for a given weight
%
% Authors: Juyul Lee and Nihar Jindal
% Date: 09/21/2005
%
% Note: The algorithm is a modified version of f_wtdratesumcap (by N. Jindal) 
%       and based on the paper:
%       H. Viswanathan, S. Venkatesan, and H. Huang
%       "Downlink capacity evaluation of celluar networks with
%       known-interference cancellation",
%       IEEE J. Sel. Area in Comm., Vol. 21, No. 5, June 2003.
%
if nargin < 5
    t_tol = 1e-4;
end;

weights=weights(:);

[M,N,K] = size(H);  
[sortedWeights, index]=sort(weights,'descend');
H = H(:,:,index);

Q = zeros(N,N,K);       % init
SNR = 10^(SNRdB/10);
grad = zeros(N,N,K);
for step=1:NUM_STEPS
    % Calculate the gradient of the objective function
    for k=1:K
        temp_grad = zeros(N,N);
        for kk=k:K-1
            temp_temp_grad = eye(M);
            for kkk=1:kk
                temp_temp_grad = temp_temp_grad + H(:,:,kkk)*Q(:,:,kkk)*H(:,:,kkk)';
            end;
            temp_grad = temp_grad+(sortedWeights(kk)-sortedWeights(kk+1))*...
                (H(:,:,k)'*inv(temp_temp_grad)*H(:,:,k));
        end;
        
        temp2_temp_grad = eye(M);
        for kkk=1:K
            temp2_temp_grad = temp2_temp_grad + H(:,:,kkk)*Q(:,:,kkk)*H(:,:,kkk)';
        end;
        grad(:,:,k) = temp_grad + sortedWeights(K)*H(:,:,k)'*inv(temp2_temp_grad)*H(:,:,k);
    end;
    
    % Determine the principal eigenvectors and eigenvalues
    for k=1:K
        [V,D]=eig(grad(:,:,k));
        D = diag(D);
        [lambda(k) pos] = max(D);
        v(:,k)=V(:,pos);
    end;
 
    [dump j]=max(lambda);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Find t 
    % Two methods available: bisection and MATLAB toolbox
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Using MATLAB optimization toolbox
    obj_t = @(t) -obj_vvh(t,j,Q,sortedWeights,H,SNR,v(:,j));
    t=fminbnd(obj_t,0,1);
    
    % Update Q
    for k=1:K
        if (k==j)
            Q(:,:,k)=t*Q(:,:,k)+(1-t)*SNR*v(:,k)*v(:,k)';
        else
            Q(:,:,k)=t*Q(:,:,k);
        end;
    end;
end;

% Calculate R
for k=1:K
    temp = eye(M);
    for kk=1:k
        temp = temp + H(:,:,kk)*Q(:,:,kk)*H(:,:,kk)';
    end;
    R(:,:,k) = temp;
end;

Rates = zeros(K,1);
Rates(index(1)) = log2(real(det(R(:,:,1))));
for k=2:K
    Rates(index(k)) = log2(real(det(R(:,:,k)))/real(det(R(:,:,k-1))));
end;

[dump index_r]=sort(index);
Q=Q(:,:,index_r);

%
% Sub-function
%
function f=obj_vvh(t,pos,Q,weights,H,P,v);
% 
%
[M,N,K] = size(H); 

for k=1:K
    if(k==pos)
        Q(:,:,k)=t*Q(:,:,k)+(1-t)*P*v*v';
    else
        Q(:,:,k)=t*Q(:,:,k);
    end;
end;

temp=0;
for k=1:K-1
    temp_temp = eye(M);
    for kk=1:k
        temp_temp=temp_temp+H(:,:,kk)*Q(:,:,kk)*H(:,:,kk)';
    end;
    
    temp=temp+(weights(k)-weights(k+1))*log2(det(temp_temp));
end;

temp_temp =eye(M);
for kk=1:K
    temp_temp=temp_temp+H(:,:,kk)*Q(:,:,kk)*H(:,:,kk)';
end;

f=temp+weights(K)*log2(det(temp_temp));