function ephisketch(m0,n0)
global m n
randomize
m = m0;
n = n0;

global pstr
pstr = 'ephisketch'
global lstr rstr

global U
P = ones(m^n,1); P = P/sum(P);
zz = P*0 + 1e-12; oo = P*0+1;
U = P;

global Ainf Binf mylinopt
mylinopt = optimset('Diagnostics','off');
mylinopt = optimset('Display','off');
[Ainf,Binf] = lipab(m,n);

global LB UB saf raf
saf = 1;
raf = n;
%raf = 100;

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;
%UB(pdim+1:end) = 0*UB(pdim+1:end) + 1e6;


if 0
  %DON'T require stationarity
  Aeq = zeros(1,xdim);
  Aeq(1:pdim) = oo';
  Beq=1;
else
  % require stationarity
  Aeq0 = zeros(1,xdim);
  Aeq0(1:pdim) = oo';
  %Aeq = is_stat(m,n);
  % require marginals uniform
  Py = ones(m,1)/m;
  [Aeq,Beq] = given_stat(m,n,Py);
  Aeq = assembleblock({Aeq,zeros(0,fdim)});
  Aeq = [Aeq;Aeq0];
  %Beq = zeros(size(Aeq,1),1); Beq(end)=1;
  Beq(end+1)=1;
end


Ain = assembleblock({zeros(0,pdim), Ainf});
Bin = Binf;
%Ain = [];
%Bin = [];

global A maxr
maxr = 0;
while 1
  %P = unifP(m,n);
  %Py = unifP(m,1);
  Py = ones(m,1)/m;
  P = randgivenstat(m,n,Py);
  
  %F = randlip;
  %F = randn(m^n,1);
  F = saf + rand(m^n,1);
  
  x = [P(:);F(:)];
  
  y = fmincon(@radrat,x,Ain,Bin,Aeq,Beq,LB,UB);
  
  fprintf('%s: PARAMS: n=%d m=%d; F in [%d,%d] \n',pstr,n,m,saf,raf);
  fprintf('RATIO_________________________ %s/%s;  maxr = %1.9f \n',lstr,rstr,maxr);
end

return


function y = radrat(x)
global m n
global lstr rstr
global LB UB saf raf
global pstr
pdim = m^n;
fdim = m^n;
xdim = pdim + fdim;

P = x(1:pdim);
F = x(pdim+1:end);
P = makepos(P);
global U
%U = prodmeas(P,m,n);

[hh,Ht,H] = gethhn(P,m,n);
%H1 = H-1;
md = max(diag(hh,1));

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

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

lhs = eval(lstr);
rhs = eval(rstr)+1e-7;
r = lhs/rhs;
y = -r;
global maxr
if r > maxr
  global optP optF
  optP = P;
  optF = F;
  maxr = r;
  fprintf('%s: PARAMS: n=%d m=%d; F in [%d,%d] \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
K = genK(m,n);
[x,fval] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n,0,mylinopt);
phi = x';
return


function c=i2c(i,m,n)
dims = repmat(m,[1 n]);
c = ind2coord(i,dims);
return


function [Aeq,Beq] = given_stat(m,n,Py)
mn = m^n;

Aeq = zeros(0,mn);
Beq = zeros(0,1);
rind = 0;
for i=1:n
  % enforce constraint: marg(X_i) = Py
  for y = 1:m
    rind = rind + 1;
    Beq(rind,1) = Py(y);
    for xi = 1:m^n
      x = i2c(xi,m,n);
      if x(i) == y
        Aeq(rind,xi) = 1;
      end
    end
  end
end

return