function [Theta, maxP, maxR]=ExpectedDualDiscreteImperfect(lambda,C,P,Pbar)
[K,M,N,L] = size(C);
% L = length(R);
maxP = zeros(K,M,N);
maxR = zeros(K,M,N);
% C = zeros(K,M,N);
% Rmat = zeros(K,M,N);
% P = zeros(K,M,N); 
% for k = 1 : K
%   for m = 1 : M
%     for n = 1 : N
%       Ptmp = PowerDiscreteImperfectCSI(G(k,m,n),sigma2w,sigma2e(k),R,Eta,BER);    
%       [C(k,m,n),idx] = max(DualDiscreteImperfect(W(k,m,n),G(k,m,n),sigma2w,sigma2e(k),Ptmp,lambda,R,Eta)); 
%       Rmat(k,m,n) = R(idx);
%       P(k,m,n) = Ptmp(idx);
%     end;
%   end;
% end;
Cl = C-lambda*P;
[maxCl,maxClidx] = max(Cl,[],4);
[maxC,maxCidx] = max(maxCl,[],2);
linidx = sub2ind(size(C),(1:K)',maxCidx,ones(size(maxCidx)),maxClidx(:,2));
% 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(1:K) = P(linidx);
maxR(linidx) = (maxClidx(:,2)-1)*2;
Theta = lambda*Pbar + mean(sum(maxC));
