function [Theta, maxP]=ExpectedDualImperfect(lambda,W,G,sigma2w,sigma2e,Pbar)
[K,M,N] = size(G);
maxP = zeros(K,M,N);
P = zeros(K,M,N);
C = zeros(K,M,N);
P0 = max(W/(log(2)*lambda)-1./G,0);
for k = 1 : K
  for m = 1 : M
    for n = 1 : N
      [k,m,n]
      if sigma2e(k)/sigma2w+G(k,m,n)<lambda/(W(k,m,n)/log(2))
        P(k,m,n) = 0;
        C(k,m,n) = 0;
      else
%         gbar = sigma2e(k)/sigma2w
%         ghat = G(k,m,n)
%         lambda
%         pmin = PowerImperfectCSI(0,G(k,m,n),sigma2w,sigma2e(k),lambda,W(k,m,n),K)
%         pmax = PowerImperfectCSI(W(k,m,n)/(log(2)*lambda),G(k,m,n),sigma2w,sigma2e(k),lambda,W(k,m,n),K)
        
        P(k,m,n) = fzero(@(x) PowerImperfectCSI(x,G(k,m,n),sigma2w,sigma2e(k),lambda,W(k,m,n),K),[0,W(k,m,n)/(log(2)*lambda)]);      
        gbar = sigma2e(k)/sigma2w;
        ghat = G(k,m,n);
        f = @(x) W(k,m,n)*log2(1+P(k,m,n)*x)/gbar.*exp(-(ghat+x)/gbar).*besseli(0,2/gbar*sqrt(ghat*x));
        C(k,m,n) = quad(f,0,ghat+gbar+10*sqrt(gbar^2+2*ghat*gbar))-lambda*P(k,m,n);
      end;      
    end;
  end;
end;
% semilogy(P)
% hold on
% semilogy(P0,'--')
% P = W/(log(2)*lambda*K)-1./G;
% P(find(P<0))=0;
% C = W.*log2(1+P.*G)-lambda*P;
[maxC maxCidx] = max(C,[],2);
linidx = sub2ind([K,M],[1:K]',maxCidx);
maxP(linidx) = P(linidx);
Theta = lambda*Pbar + mean(sum(maxC));
