function [Theta,Iter,maxP,maxR]=ExpectedDual2(l,w,gbar,Pbar,K,ComputePR)
M = length(w);
a = 1e-10;
alpha = -log(2)*l./w;
xmax = (-w.*log2(alpha./gbar/log(a))+w/log(2).*alpha./gbar/log(a)-w/log(2));
c = zeros(size(w));
iter = zeros(size(w));
for m = 1 : M
  [c(m),iter(m)] = quad(@(x) myintegrand(x,w,gbar,l,m),0,xmax(m));
end
Iter = mean(iter);
Theta = l*Pbar + K*sum(c);
if ComputePR
  gmax = groot(xmax,w,l);
  g0 = l*log(2)./w;
  maxP = zeros(size(w));
  maxR = zeros(size(w));
  for m = 1 : M
    maxP(m) = quad(@(x) powintegrand(x,w,gbar,l,m),g0(m),gmax(m));
    maxR(m) = quad(@(x) capintegrand(x,w,gbar,l,m),g0(m),gmax(m));
  end
else
  maxP = 0;
  maxR = 0;
end
end


function y = myintegrand(x,w,gbar,l,midx)
y = x.*PDFcont(x,w(midx),l,gbar(midx))./CDFcont(x,w(midx),l,gbar(midx)).*CDFProd(x,w,l,gbar);
y(isinf(y)|isnan(y))=0;
end

function y = powintegrand(x,w,gbar,l,midx)
y = (w(midx)/log(2)/l-1./x).*exppdf(x,gbar(midx))./CDFcont(gdual(x,w(midx),l),w(midx),l,gbar(midx)).*CDFProd(gdual(x,w,l),w,l,gbar);
y(isinf(y)|isnan(y))=0;
end

function y = capintegrand(x,w,gbar,l,midx)
y = (w(midx)*log2(w*w/log(2)/l)).*exppdf(x,gbar(midx))./CDFcont(gdual(x,w(midx),l),w(midx),l,gbar(midx)).*CDFProd(gdual(x,w,l),w,l,gbar);
y(isinf(y)|isnan(y))=0;
end

function y = CDFProd(x,w,l,gbar)
y = ones(size(x));
M = length(w);
for m = 1 : M
  y = y.*CDFcont(x,w(m),l,gbar(m));
end
y(isinf(y)|isnan(y))=0;
end

function y = gdual(x,w,l)
y = w*log2(w*x/log(2)/l)-w/log(2)+l./x;
end