function [UsrVec Pvec] = scheduler_ofdma_opt(G,w,Pmax,scheduling_mask)
% SCHEDULER_OFDMA_OPT: Optimal scheduler for downlink OFDMA
% [UsrVec Pvec] = scheduler_ofdma_opt(G,w,Pmax,scheduling_mask)
% Author: Ian C. Wong
% Revised: Aug. 18, 2006
% Inputs: 
%   G:    M x K matrix (M-users, K-subchannels) of CINRs
%   w:    M-length vector of user weights.  In PF scheduling, this is simply
%           the reciprocal of the average rates so far normalized to have
%           sum of 1.
%   Pmax: Maximum power in Watts (46dBm=39.8107 W)
%   scheduling_mask: MxK matrix of binary integers indicating whether a
%         user is allowed to use the subchannel (=1) or not (=0).
% Outputs:
%   UsrVec: K-length vector of user indices indicating user assignment for
%           each subchannel
%   Pvec: K-length vector of powers for each subchannel

w = w(:);
[M,K] = size(G);
W = (w*ones(1,K));
% Minimum lagrangian multiplier
lambdamin = 1/log(2)*min(w)/(Pmax+sum(max(1./G)));
% Maximum lagrangian multiplier
lambdamax = K*max(w)/log(2)/Pmax;
% Find the minimum of the
[lambda,fval,exitflag,output] = fminbnd(@(x) Dual(x,W,G,Pmax,scheduling_mask,0),lambdamin,lambdamax);
[Theta, maxP]=Dual(lambda,W,G,Pmax,scheduling_mask,1);
[UsrVec k] = find(maxP>0);
Pvec = sum(maxP*Pmax/sum(sum(maxP)));
end

function [Theta, maxP]=Dual(lambda,W,G,Pmax,scheduling_mask,ReturnPower)
maxP = zeros(size(G));
P = max(W/(log(2)*lambda)-1./G,0);
C = (W.*log2(1+P.*G)-lambda*P).*scheduling_mask;
[maxC maxCidx] = max(C);
if ReturnPower
  linidx = sub2ind(size(C),maxCidx,1:size(G,2));
  maxP(linidx) = P(linidx);
else
  maxP = 0;
end;
Theta = lambda*Pmax + sum(maxC);
end