function ephimody(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 = 1e6;
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;


Aeq0 = zeros(1,xdim);
Aeq0(1:pdim) = oo';

if islip
  Ain = assembleblock({zeros(0,pdim), Ainf});
  Bin = Binf;
else
  Ain = [];
  Bin = [];  
end

global maxr mind
maxr = 0;
mind = 1e6;

while 1
  % pick an arb. marg. distr.
  %Pmarg = getmargs(P,m,n);

  %Pmarg = pnormdim(rand(m,n),1);
  %Pmarg = pnormdim(ones(m,n),1);
  %P = randgivenmarg(m,n,Pmarg);
  
  P = unifP(m,n);
  % and corresponding marg. distr
  %U = prodmeas(P,m,n);
  %U = P0;
  
  
  % constrain P to have these marginals
  %[Aeq,Beq] = constrain_marg(m,n,Pmarg);
  %Aeq = assembleblock({Aeq,zeros(0,fdim)});
  %Aeq = [Aeq;Aeq0];
  %Beq(end+1)=1;
  Aeq = Aeq0; Beq = 1;
  
  %saf = rand;
  saf = 1;
  LB(pdim+1:end) = 0*LB(pdim+1:end) + saf;
  
  F = randlip;
  
  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]\n',pstr,n,m,saf,raf);
  %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);
P = makepos(P);
F = x(pdim+1:end);

t = 1;

pt = getmargX(P,m,n,t);
U = prodmeas(P,m,n);

[hh,Ht,H] = gethhn(P,m,n);

HUF = H * U'*F;

sumy = 0;
for y=1:m
  Py = getkyt(P,y,m,n,t);
  Py = Py/sum(Py);
  [hh,Ht,Hy] = gethhn(Py,m,n-1);
  Uy = prodmeas(Py,m,n-1);
  Fy = getkyt(F,y,m,n,t);
  
  sumy = sumy + pt(y)*Hy*Uy'*Fy;
  %sumy = sumy + pt(y)*Py'*Fy;
end

rstr = 'HUF';
lstr = 'sumy';


lhs = eval(lstr);
rhs = eval(rstr)+1e-7;

global maxr mind
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]\n',pstr,n,m,saf,raf);
  fprintf('RATIO_________________________ (%s)/(%s);  maxr = %1.9f \n',lstr,rstr,maxr);
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


function y = radrat0(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);
F = max(F,saf);
F = min(F,raf);
P = makepos(P);
global U

[hh,Ht,H] = gethhn(P,m,n);
%Flip = max(Ainf * F);
%Flip = max(Flip,1);
%a = min(min(F),1)+1e-9;
a = min(saf,1);

FP = abs(F'*P);
FU = abs(F'*U);

lstr = 'FP';
%rstr = 'FU*H';
%rstr = 'FU*H*Flip';
rstr = 'FU*H/a';
%rstr = 'FU';
%rstr = 'FU*(1+md)';
%rstr = 'FU*max(F)';

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]\n',pstr,n,m,saf,raf);
    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