function psixy1(m0,n0)
global doplots domask mobs mhid n m
randomize
m = m0;
n = n0;
chaincheck
return




function chaincheck
global n m P T psi
global mrat mT mpsi H0
global Rmax
mrat = 0;
global longpsi 
longpsi = 1;
T = chaingraph(n);
Rmax = 50;
while 1
  if longpsi
    [nodesi,nodesj]=find(T);
    ii = find(nodesi<nodesj);
    nodesi = nodesi(ii);
    nodesj = nodesj(ii);
    ne = length(nodesi);
    psi = rand(m,m,ne); 
  else
    psi = rand(m); 
  end
  psi = psi/min(psi(:));
  psi = min(psi,Rmax);
  x = psi(:);
  %[x,fval] = fmincon(@chainopt,x,[] ,[]  ,[] ,[]  ,0*x+1,[]);
  [x,fval] = fmincon(@chainopt,x,[] ,[]  ,[] ,0*x+Rmax  ,0*x+1,[]);  
end
return


function y = chainopt(x)
global n m P T
global longpsi
global Rmax
x = abs(x);
if longpsi
  [nodesi,nodesj]=find(T);
  ii = find(nodesi<nodesj);
  nodesi = nodesi(ii);
  nodesj = nodesj(ii);
  ne = length(nodesi);
  psi = reshape(x,[m m ne]);
else
  psi = reshape(x,[m m]);
end
psi = min(psi,Rmax);
global P0 Ztree0 Ztree psi0
P0 = P; Ztree0 = Ztree; psi0 = psi;
i = 2;
psi0(:,:,i) = ones(m);
psi = psi0;
P = treefillprob(m,n,T,psi);


%[b,A,C,dd] = is_markov(P,m,n);
%A = A(:,:,i);

y1 = 1;
y2 = 2;
alf = zeros(m,1);
bet = zeros(m,1);

for x1=1:m
  alf(x1) = axy(i,x1,y1);
  bet(x1) = axy(i,x1,y2);
end

%dd = totdif([A(:,1),A(:,2)],pnormdim([alf bet],1));
%if dd>1
%  keyboard
%end

%p = [A(:,1);A(:,2)];
p = [alf; bet];
pr = max(p)/min(p);

aa = alf/sum(alf);
bb = bet/sum(bet);
p = [aa; bb];
pr = max(p)/min(p);
th = sum(abs(aa-bb))/2;

%th = sum(abs(A(:,1)-A(:,2)))/2;
[R,rr] = psi2r(psi0);
R = rr(i);

%B = (R-1)/(R+1) + eps;

%r = th/B;
%r = pr/R^3;
r = th;
keyboard

%if r>1.5
%  keyboard
%end

y = -r;

global mrat mT mpsi
if r > mrat
  mrat = r;
  mT = T;
  mpsi = psi;
end
fprintf('chainxplore n=%d   MRAT = %1.9f\n',n,mrat);

return





function p = PSI(x,y,i)
global psi psi0
psi = psi0;
% for technical reasons, we take x and y to be strings
% and PSI is psi(x(end),y(1))
if (isempty(x) | isempty(y))
  p = 1;
else
  p = psi(x(end),y(1),i);
end
return

function a = axy(i,x,y)
global n m P T psi Ztree psi0
a = 0;
for vi=1:m^(i-1)
  v = ind2coord(vi,dimn(i-1));
  for zi=1:m^(n-i-1)
    z = ind2coord(zi,dimn(n-i-1));
    s1 = Pi(v,1,i-2)*PSI(v(1:i-1),y,i-1)*PSI(y,x,i)*PSI(x,z,i+1)*Pi(z,i+2,n);
    s2 = getPX(P*Ztree,m,n,[v y x z],[1:n]);
    if abs(s1-s2) > 1e-9
      'discrepancy s12'
      keyboard
    end
    a =  a + s1;
  end
end
a1 = getPX(P*Ztree,m,n,[y x],[i i+1]);
if abs(a1-a) > 1e-9
  'discrepancy a01'
  keyboard
end
return


function d = dimn(n)
global m
d = repmat([m],[1 n]);
return


function p = Pi(u,i,j)
global n m T psi psi0
psi = psi0;
p = 1;
for k=i:j-1
  p = p * psi(u(k-i+1),u(k-i+2),k);
end
return