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

global pstr
pstr = 'tvmarkov'

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


[Aeq, Beq] = normab(m,m);

global maxr
global lstr rstr
maxr = 0;
while 1
  A = pnormdim(rand(m),1);
  x = A(:);
  y = fmincon(@myobj,x,[],[],Aeq,Beq,0*x+1e-9,0*x+1);
  fprintf('%s: n=%d m=%d; rat = %s/%s;  maxr = %1.9f \n',pstr,n,m,lstr,rstr,maxr);
  %fprintf('%s: n=%d;  maxr = %1.9f \n',pstr,n,maxr);
end
return



function y = myobj(x)
global pstr
global m n maxr
global lstr rstr
global Ainf Binf mylinopt
A = makepos(reshape(x,[m m]));
C = ones(m,1)/m;
P = markfillprob(A,C,n);
U = prodmeas(P,m,n);

%lhs = TV(P-U);
%rhs = h + 1e-9;
%lhs = phinorm(P-U);
[h,hh] = stricthh(A);
Hn = sum( h.^(0:n-1) );
%rhs = sum( h.^(0:n-1) );
tha=h;

%[hh,Ht,H] = gethhn(P,m,n);
%inex = 0;
%for k=1:n-1
%  inex = inex + inclexcl(diag(hh,k));
%end

K = P - Hn*U;

%[x,fval] = linprog(-K,Ainf,Binf,[],[],0*K,0*K+n,0,mylinopt);
%lstr = '-fval';

%lstr = 'phinorm(P-U)';
lstr = 'Psi(P-Hn*U,m,n)';
%lstr = 'TV(P-U)';

%rstr = '1/(1-tha)';
%rstr = 'log(n)/(1-tha^.5)';
rstr = 'Hn-1';
%rstr = 'tha/(1-tha)';
%rstr = 'sqrt(tha)';
%rstr = '1/(1-sqrt(tha))';
%rstr = '(1/(1-tha^.5))^2';
%hh = repmat([h],[1 n-1]);
%hn = h.^(1:n-1); % NO
%rstr = 'inclexcl(hh)';
%rstr = 'inex';

lhs = eval(lstr);
rhs = eval(rstr)+1e-9;
r = lhs/rhs;
y = -r;
if r > maxr
  global optP optA
  optP = P;
  optA = A;
  maxr = r;
  %fprintf('%s: n=%d;  maxr = %1.9f \n',pstr,n,maxr);
  fprintf('%s: n=%d m=%d; rat = %s/%s;  maxr = %1.9f \n',pstr,n,m,lstr,rstr,maxr);  
end
return