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

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


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

if 0 % require sym. stationarity
  Aeq = is_stat(m,n);
  Aeq = [Aeq;oo'];
  Beq = zeros(size(Aeq,1),1); Beq(end)=1;
else
  Aeq = oo';
  Beq = 1;
end

global maxr
global lstr rstr
maxr = 0;
while 1
  %Py = rand(m,1); Py = Py/sum(Py);
  %P = randgivenstat(m,n,Py);
  P = rand(m^n,1); P = P/sum(P);
  
  x = [P(:)];
  
  y = fmincon(@radrat,x,[],[],Aeq,Beq,zz,oo);
  
  fprintf('n=%d m=%d; rat = %s/%s;  maxr = %1.9f \n',n,m,lstr,rstr,maxr);
end

return



function y = radrat(P)
global m n maxr
global lstr rstr
P = makepos(P);
U = prodmeas(P,m,n);

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

%md = 0;
inex = 0;
for k=1:n-1
  %md = md + max(diag(hh,k));
  inex = inex + inclexcl(diag(hh,k));
end

%lhs = Psi(P-U,m,n);
%rhs = H;
%rhs = 1+md;
%lhs = phinorm(P-U);
%rhs = TV(P-U);
%rhs = sum(alf);
lstr = 'psinorm(P-U,m,n)';


G=hh+eye(size(hh));
L2 = max(eig(G'*G));
%rhs = H-1;  % does NOT bound phinorm
%rhs = H;
%rhs = 1;
%rhs = max(diag(hh,1));
%rstr = '(H-1)';
%rstr = 'inex';
rstr = 'L2';

%tv = TV(P-U);
%lstr = 'phinorm(P-U)';
%rstr = 'tv/(1-tv)';


lhs = eval(lstr);
rhs = eval(rstr);

r = lhs/rhs;

y = -r;

if r > maxr
  global optP
  optP = P;
  maxr = r;
  fprintf('n=%d m=%d; rat = %s/%s;  maxr = %1.9f \n',n,m,lstr,rstr,maxr);
end


return


function y = phinorm(K)
global Ainf Binf m n mylinopt

[x,fvalp] = linprog( K,Ainf,Binf,[],[],0*K,0*K+n,0,mylinopt);
[x,fvaln] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n,0,mylinopt);
y = max(abs(fvalp),abs(fvaln));
return


function P = prodmeas0(Py,n)
P = Py;
for i=1:n-1
  P = tensprod1(P,Py);
end
return


function H1 = gethhn1(P,m,n)
H1 = 1;
i = 1;
for j=i+1:n
  h = mxetaij(P,m,n,i,j);  % the new one, as in Samson
  H1 = H1 + h;
end
return