function ephiy(m0,n0) global m n randomize m = m0; n = n0; global pstr pstr = mfilename global lstr rstr global islip islip = 1; global dorat dorat = 1; global U P0 P = ones(m^n,1); P = P/sum(P); zz = P*0 + 1e-12; oo = P*0+1; U = P; P0 = P; global LB UB saf raf Lipk saf = 1; raf = n; Lipk = 1; %raf = 100; %if islip if 1 global Ainf Binf mylinopt mylinopt = optimset('Diagnostics','off'); mylinopt = optimset('Display','off'); [Ainf,Binf] = lipab(m,n); Binf = Binf * Lipk; end pdim = m^n; fdim = m^n; xdim = pdim + fdim; LB = zeros(xdim,1); UB(1:pdim) = zz; LB(pdim+1:end) = 0*LB(pdim+1:end) + saf; UB = LB*0; UB(1:pdim) = oo; UB(pdim+1:end) = 0*UB(pdim+1:end) + raf; Aeq = zeros(1,xdim); Aeq(1:pdim) = oo'; Beq = 1; if islip Ain = assembleblock({zeros(0,pdim), Ainf}); Bin = Binf; else Ain = []; Bin = []; end global maxr mind maxr = 0; mind = 1e6; while 1 P = unifP(m,n); if islip F = randlip; else F = rand(m^n,1)+saf; end x = [P(:);F(:)]; y = fmincon(@radrat,x,Ain,Bin,Aeq,Beq,LB,UB); fprintf('%s: PARAMS: n=%d m=%d; F in [%1.3f,%1.3f], islip? %d \n',pstr,n,m,saf,raf,islip*Lipk); %if dorat fprintf('RATIO_________________________ (%s)/(%s); maxr = %1.9f \n',lstr,rstr,maxr); %else fprintf('DIFF _________________________ %s - %s; mind = %1.9f \n',rstr,lstr,mind); %end dorat = ~dorat; end return function y = radrat(x) global m n global Ainf global lstr rstr global LB UB saf raf global islip global Lipk global dorat global pstr pdim = m^n; fdim = m^n; xdim = pdim + fdim; P = x(1:pdim); F = x(pdim+1:end); F = pl(F); P = makepos(P); [U,margs] = prodmeas(P,m,n); [hh,Ht,H] = gethhn(P,m,n); ysum = 0; for y=1:m py = margs(y,n); Fy = getky(F,y,m,n); Uy = getky(U,y,m,n); %Fy = min(Fy,n-1); Py = getky(P,y,m,n); P_y = Py/sum(Py); [hh,Ht,Hy] = gethhn(P_y,m,n-1); %ysum = ysum + py * Hy * Fy'*P_y; ysum = ysum + py * Hy * Fy'*Uy; end FU = abs(F'*U); lstr = 'ysum'; rstr = 'FU*H'; lhs = eval(lstr); rhs = eval(rstr)+1e-7; global maxr mind if dorat r = lhs/rhs; y = -r; if (r > maxr) global maxrP maxrF maxrP = P; maxrF = F; maxr = r; fprintf('%s: PARAMS: n=%d m=%d; F in [%1.3f,%1.3f], islip? %d \n',pstr,n,m,saf,raf,islip*Lipk); fprintf('RATIO_________________________ (%s)/(%s); maxr = %1.9f \n',lstr,rstr,maxr); if r > 1.01 save BADEPHIR m n P F r end end else d = rhs - lhs; y = d; if (d < mind) global mindP mindF mindP = P; mindF = F; mind = d; fprintf('%s: PARAMS: n=%d m=%d; F in [%1.3f,%1.3f], islip? %d \n',pstr,n,m,saf,raf,islip*Lipk); fprintf('DIFF _________________________ %s - %s; mind = %1.9f \n',rstr,lstr,mind); if d < -.01 save BADEPHID m n P F d end end end return function phi = randlip global Ainf Binf m n mylinopt global saf raf K = genK(m,n); [x,fval] = linprog(-K,Ainf,Binf,[],[],0*K+saf,0*K+raf,0,mylinopt); phi = x'; return function Pmarg = getmargs(P,m,n) Pmarg = zeros(m,n); for i=1:n Pmarg(:,i) = getmargX(P,m,n,i); end return