function [Theta, maxP, maxR]=ExpectedDualDiscrete(lambda,W,G,Pbar,R,Eta)
[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]);

maxP = zeros(K,M,N);
maxR = zeros(K,M,N);
C = zeros(K,M,N);
Rmat = zeros(K,M,N);
Ptmp = Etamat./repmat(G,[1,1,1,L]);
[C,idx] = max(repmat(W,[1,1,1,L]).*R-lambda*Ptmp,4);
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;
[maxC, maxCidx] = max(C,[],2);
linidx = sub2ind([K,M],[1:K]',maxCidx);
% 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));
