function [Theta, maxP, maxR]=ExpectedDualDiscrete2(lambda,W,G,Pbar,R,Eta,ComputePR)
[K,M,N] = size(G);
L = length(R);
Etatmp = zeros([1,1,1,L]);
Etatmp(1,1,1,:) = Eta;
Etamat = repmat(Etatmp,[K,M,N,1]);
Rtmp = zeros([1,1,1,L]);
Rtmp(1,1,1,:) = R;
Rmat = repmat(Rtmp,[K,M,N,1]);
Pmat = Etamat./repmat(G,[1,1,1,L]);
[C,idx] = max(repmat(W,[1,1,1,L]).*Rmat-lambda*Pmat,[],4);
[maxC, maxCidx] = max(C,[],2);
Theta = lambda*Pbar + mean(sum(maxC));

if ComputePR
  maxP = zeros(K,M,N);
  maxR = zeros(K,M,N);
  Nidx = zeros([1,1,N]);
  Nidx(1,1,:) = 1:N;
  I1 = repmat(repmat((1:K)',[1,M]),[1,1,N]);
  I2 = repmat(repmat(1:M,K,1),[1,1,N]);
  I3 = repmat(Nidx,[K,M,1]);  
  linidx = sub2ind(size(Rmat),I1,I2,I3,idx);
  Rtmpidx
  maxRtmp = Rmat(linidx(:));
  maxPtmp = Pmat(linidx(:));
  linidx2 = sub2ind(size(Rmat),repmat((1:K)',[1,1,N]),maxCidx,);
  maxP(linidx2) = maxPtmp(linidx2);
  maxR(linidx2) = maxRtmp(linidx2);
else
  maxP = 0;
  maxR = 0;
end


% for k = 1 : K
%   for m = 1 : M
%     for n = 1 : N
%       Ptmp = Eta/G(k,m,n);    
%       [C(k,m,n),idx] = max(W(k,m,n)*R-lambda*Ptmp); 
%       Rmat(k,m,n) = R(idx);
%       P(k,m,n) = Ptmp(idx);
%     end;
%   end;
% end;
% 
% % P = W/(log(2)*lambda)-1./G;
% % P(find(P<0))=0;
% % C = W.*log2(1+P.*G)-lambda*P;
% % maxC = max(C,[],2);
% maxCbig = repmat(maxC,[1,2]);
% maxCidx = find(C==maxCbig);
% maxP(linidx) = P(linidx);
% maxR(linidx) = Rmat(linidx);
% Theta = lambda*Pbar + mean(sum(maxC));
