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

pstr = mfilename
global pstr

pdim = m*m;
xdim = pdim;

global G
G = makecyc(n);

global nmat
nmat = zeros(n);
for i=1:n
  for j=i+1:n
    nmat(i,j) = j-i;
  end
end

global maxr
maxr = 0;

while 1
  psi = rand(m);
  
  x = [psi(:)];
  y = fmincon(@myobj,x,[],[],[],[],0*x,[]);
  
  fprintf('%s: n=%d; m=%d; maxr = %1.9f \n',pstr,n,m,maxr);
  
end


return


function y = myobj(x)
global m n
global pstr
global G

psi = pl(reshape(x,[m m]));
P = cycfillprob(m,n,G,psi);

MR = max(psi(:));
mr = min(psi(:));

global nmat
tha = (MR-mr)/(MR+mr);
%than = sum(tha.^(0:n-1));
than = tha.^floor(nmat/2);
than = max(than,1e-7);
[hh,Ht,H] = gethhn(P,m,n);


%lhs = H;
%rhs = 1/(1-tha);
%rhs = than;
%r = lhs/rhs;
r = max(hh(:)./than(:));

y = -r;
global maxr
if r > maxr
  global optpsi
  optpsi = psi;
  maxr = r;
  fprintf('%s: n=%d; m=%d; maxr = %1.9f \n',pstr,n,m,maxr);
end


return


function P = cycfillprob(m,n,T,psi)
global Ztree
P = zeros(m^n,1);
dims = repmat([m],[1 n]);
[nodesi,nodesj]=find(T);
ii = find(nodesi<nodesj);
nodesi = nodesi(ii);
nodesj = nodesj(ii);
ne = length(nodesi);
if size(psi,3) < ne
  psi = repmat(psi,[1 1 ne]);
end
for xind = 1:m^n
  x = ind2coord(xind,dims);
  p = 1;
  for e=1:length(nodesi)
    i = nodesi(e);
    j = nodesj(e);
    p = p * psi(x(i),x(j),e);
  end
  P(xind) = p;
end
Ztree = sum(P);
P = P/Ztree;
return