function mumutil(m0,n0)
global m n
m = m0;
n = n0;

global pstr
pstr = 'mumutil'
global lstr rstr


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

pdim = m^n;
xdim = pdim;
LB = zeros(xdim,1);
UB = zz;
UB(1:pdim) = oo;

if 0
  Aeq = zeros(1,xdim);
  Aeq(adim+1:end) = oo';
  Beq=1;
else
  % require stationarity
  Aeq0 = oo';
  Aeq = is_stat(m,n);
  Aeq = [Aeq;Aeq0];
  Beq = zeros(size(Aeq,1),1); Beq(end)=1;
end


global maxr
maxr = 0;
while 1
  %P = unifP(m,n);
  Py = unifP(m,1);
  P = randgivenstat(m,n,Py);

  x = [P(:)];
  
  y = fmincon(@radrat,x,[],[],Aeq,Beq,LB,UB);
  
  fprintf('%s: n=%d m=%d; rat = %s/%s;  maxr = %1.9f \n',pstr,n,m,lstr,rstr,maxr);
end

return


function y = radrat(x)
global m n
global lstr rstr
global pstr

P = makepos(x);
U = 0*P+1; U=U/sum(U);
%U = prodmeas(P,m,n);
[hh,Ht,H] = gethhn(P,m,n);

U = U+1e-8;
%maxPU = max(P./U);
PUrat = max( (1/n+P) ./ (1/n+U));

lstr = 'PUrat';
rstr = 'H';


lhs = eval(lstr);
rhs = eval(rstr);
r = lhs/rhs;
y = -r;
global maxr
if r > maxr
  global optP
  optP = P;
  maxr = r;
  fprintf('%s: n=%d m=%d; rat = %s/%s;  maxr = %1.9f \n',pstr,n,m,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