function [Theta, maxP, maxR, iterint]=ExpectedDualDiscrete3(lambda,W,gbar,Pbar,R,Eta,K,ComputePR)
% [K,M,N] = size(G);
Rext = [-inf,R,R(end)];
Etaext = [0,Eta,inf];
M = length(W);
cvec = zeros(size(W));
itervec = zeros(size(W));
for m = 1 : M
  [cvec(m) itervec(m)] = quad(@(x) discinteg(x,lambda,W,gbar,Etaext,Rext,m),0,W(m)*R(end));
end
Theta = lambda*Pbar + K*sum(cvec);
if ComputePR
  Slopes = (Rext(2:end)-Rext(1:end-1))./(Etaext(2:end)-Etaext(1:end-1));
  for m = 1: M
    gbounds = W(m)*(Rext(2:end-1)-Etaext(2:end-1).*Slopes(2:end));
    maxR(m) =   CDFProd(gbounds(2:end),l,w,gbar,Etaext,Rext)
  end
  maxP = 0;
  maxR = 0;
  iterint = mean(itervec);
else
  maxP = 0;
  maxR = 0;
end

% 
% maxP = 0;
% maxR = 0;
% 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);
%   maxPtmp = zeros(K,M,N);
%   maxRtmp = 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);
%   maxRtmp(1:numel(maxR)) = Rmat(linidx(:));
%   maxPtmp(1:numel(maxP)) = Pmat(linidx(:));
%   linidx2 = sub2ind(size(C),repmat((1:K)',[1,1,N]),maxCidx,repmat(Nidx,[K,1,1]));
%   maxP(linidx2) = maxPtmp(linidx2);
%   maxR(linidx2) = maxRtmp(linidx2);
% else
%   maxP = 0;
%   maxR = 0;
% end
end


function y = discinteg(x,l,w,gbar,Etaext,Rext,midx)
y = x.*lmaxpdf3vec(x,l,w(midx),gbar(midx),Etaext,Rext)./lmaxcdf3vec(x,l,w(midx),gbar(midx),Etaext,Rext).*CDFProd(x,l,w,gbar,Etaext,Rext);
end

function y = CDFProd(x,l,w,gbar,Etaext,Rext)
y = ones(size(x));
M = length(w);
for m = 1 : M
  y = y.*lmaxcdf3vec(x,l,w(m),gbar(m),Etaext,Rext);
end
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));
