function dtreesketch(m0,n0)
global n m
global SUBT
randomize
m = m0;
n = n0;
global mrat
mrat = 0;

while 1
  SUBT = [];
  treecheck
  fprintf('n=%d mrat = %1.9f \n',n,mrat);
  if mrat > 1.01
    keyboard
  end
end

%trexplore
%trespc
%doplot
%strucsearch
return

function maxtha
global A mtha
[ni,nj,m,m] = size(A);
mtha = 0;
for i=1:ni
  for j=1:nj
    Aij = squeeze(A(i,j,:,:));
    mtha = max(mtha,strict(Aij));
  end
end
return

function treecheck
global n m hh0
global A C T P
T = randdtree(n);
%fstr='1<2;1<3';
%fstr='1<2<3';
%fstr = '1<2;1<3;2<4;3<5';
%fstr = '1<2;1<3;2<4;3<5;5<6';
%fstr = '1<2;1<3;2<4;'
%T = readtree(fstr);
[n,n] = size(T);
[nodesi,nodesj]=find(T);
L = treewidth(T);
th0 = 1/L;
%A = zeros(length(nodesi),length(nodesi),m,m);
A=[];
A0 = pnormdim(rand(m,m),1);
for e=1:length(nodesi)
  e1 = nodesi(e);
  e2 = nodesj(e);
  %Ae = pnormdim(rand(m,m),1);
  Ae = randAth(m,th0);
  %Ae = A0;
  Ae = reshape(Ae,1,1,m,m);
  A(e1,e2,:,:) = Ae;
end
maxtha;
C = pnormdim(rand(m,1),1);
P = dtreefillprob(m,n,T,A,C);
[hh0,Ht,H] = gethhn(P,m,n);
for i=1:n
  for j=i+1:n
    etaij(P,m,n,i,j);
    h = mxetaij(P,m,n,i,j);  % the new one, as in Samson
    %tp = thaprod(i,j);
    h1 = treehij(A,C,T,i,j) + 1e-14;
    
    global mrat mA mC mT
    r = h/h1;
    if r > 1.01
      keyboard
    end
    if r > mrat
      mrat = r;
      mA = A;
      mT = T;
      mC = C;
    end
    %good = (h<=tp+1e-9)
    %if ~good
    %  'bad bound'
    %  keyboard
    %end
  end
end
return

function [h,xim] = etaij(P,m,n,i,j)
global hh0
h = 0;
xim = 0;
dims = repmat([m],[1 i-1]);
for xi = 1:m^(i-1)
  X = ind2coord(xi,dims);
  hX = etaijX0(P,m,n,i,j,X);
  if hX > h
    h = hX;
    xim = xi;
  end
end
return


function h = etaijX0(P,m,n,i,j,X)
global A C T
h = 0;
for y1 = 1:m
  for y2 = y1+1:m
    if 1
      P0 = getcondP(P,m,n,[j:n],[1:i],[X y1]);
      P1 = getcondP(P,m,n,[j:n],[1:i],[X y2]);
      h0 = .5*sum(abs(P0-P1));
      %h1 = treesubij1(i,j,y1,y2);
      %h1 = treesubtensub(i,j,y1,y2);
      %h1 = etaijXdtree23(P,m,n,i,j,X,y1,y2);
      h1 = etaijtreetens3(i,j,y1,y2);
      g1 = abs(h0-h1)<1e-9;
      %h1 = etaijXdtree22ub(P,m,n,i,j,X,y1,y2);
      %g2 = h0<h1+1e-9;
      %good = g1 & g2
      good = g1;
      if ~good
        keyboard
      end
      
      %if (i==1)&(j==4)
      %  h = n5case(P,m,n,i,j,X,y1,y2)
      %  good = abs(h0-h)<1e-9
      %  if ~good
      %    keyboard
      %  end
      %end

      %if 0&(i==1)&(j==2)
      %  h = n3case(P,m,n,i,j,X,y1,y2);
      %  good = abs(h0-h)<1e-9
      %  if ~good
      %    keyboard
      %  end
      %end

      
    end
  end
end
return


function h = n5case(P,m,n,i,j,y,w1,w2)
global A

%Tstr = 1<2;1<3;2<4;3<5;


p1 = squeeze(A(1,2,:,w1));
p2 = squeeze(A(1,2,:,w2));
q1 = squeeze(A(1,3,:,w1));
q2 = squeeze(A(1,3,:,w2));
p1q1 = tensprod(p1,q1);
p2q2 = tensprod(p2,q2);

A24 = squeeze(A(2,4,:,:));
A35 = squeeze(A(3,5,:,:));
A2345 = ABprod(A24,A35);

h = TV(A2345*(p1q1-p2q2));
return

function h = n3case(P,m,n,i,j,y,w1,w2)
global A

if 0 % chain topology
h = 0;
for x3 = 1:m
  S = 0;
  for x2 = 1:m
    S = S + A(2,3,x3,x2)*(A(1,2,x2,w1)-A(1,2,x2,w2));
  end
  h = h + .5*abs(S);
end


else % hat topology
p1 = squeeze(A(1,2,:,w1));
p2 = squeeze(A(1,2,:,w2));
q1 = squeeze(A(1,3,:,w1));
q2 = squeeze(A(1,3,:,w2));
p1q1 = tensprod(p1,q1);
p2q2 = tensprod(p2,q2);
h = TV(p1q1-p2q2);

end
  %lhs = TV(p1q1,p2q2);
  %rhs = max(TV(p1,p2),TV(q1,q2));
  
  %r = lhs / rhs;
  
  %if r>1
  %  keyboard
  %end

return

h = 0;
for x2 = 1:m
  for x3 = 1:m
    pp1 = A(1,2,x2,w1)*A(1,3,x3,w1);
    pp2 = A(1,2,x2,w2)*A(1,3,x3,w2);    
    
    h = h + .5*abs(pp1-pp2);
  end
end
return


function pq = tensprod(p,q)
m = length(p);
pq = zeros(m^2,1);
for x = 1:m
  for y = 1:m
    xyind = coord2ind([x y],[m m]);
    pq(xyind) = p(x) * q(y);
  end
end
return

function AB = ABprod(A,B)
[mA,nA] = size(A);
[mB,nB] = size(B);
AB = [];
for x=1:nA
  for y = 1:nB
    abind = coord2ind([x y],[nA nB]);
    AB(:,abind) = tensprod(A(:,x),B(:,y));
  end
end
return



function d = TV(x)
d = .5 * sum(abs(x));
return



function tp = thaprod(i,j)
global A T mtha
dp = dpath(T,i,j);
L = treewidth(T);
tp = ((2*L-1)*mtha)^(length(dp)-1);
return
tp = 0;
for ki = 1:length(dp)-1
  u = dp(ki);
  v = dp(ki+1);
  Auv = squeeze(A(u,v,:,:));
  %tp = tp * strict(Auv);
  tp = max(tp,strict(Auv));
end
tp = tp^(length(dp)-1);
return


function h = etaijXdtree3(P,m,n,i,j,y,w1,w2)
global A C T
h = 0;

dp = dpath(T,i,j);
ldp = length(dp)-1;
Ddp = repmat([m],[1 ldp]);

if dp(1)~=i
  h = 0;
  return
end

Sum = 0;
for si=1:prod(Ddp)
  s = ind2coord(si,Ddp);

  f = 1;
  for t=2:ldp
    f = f * A(dp(t),dp(t+1),s(t),s(t-1));
  end

  g1 = 1; g2 = 1;
  ikids = find(T(i,:));
  for vi = 1:length(ikids)
    v = ikids(vi);
    g1 = g1 * A(i,v,s(1),w1);
    g2 = g2 * A(i,v,s(1),w2);
  end
  
  %Sum = Sum + f*g1 - f*g2;
  Sum = Sum + abs(f*g1 - f*g2);
end
h = .5*abs(Sum);
return


function h = etaijXdtree22(P,m,n,i,j,y,w1,w2)
global A C T
global LEVELS
h = 0;
LEVELS = repmat({[]},n,1);
Ti = getSubt(T,i);
jj = Ti(Ti>=j);
if isempty(jj)
  return
end
getLevels(T,1,0);
j0 = min(jj);
levi0 = getLevu(i);
levj0 = getLevu(j0);

%xind = intersect(LEVELS{levj0},j0:n);
%xind = intersect(xind,Ti);
Ej = getEj(T,j0);
%if ~isempty(setdiff(Ej(:,2),Ti))
%  'aha'
%  keyboard
%end
% really do need the intersect!!
% in both cases
xind = intersect(Ej(:,2),Ti);
DX = repmat([m],[1 length(xind)]);
zind = intersect(Ti,i+1:j0-1);
% YES, really do need zind def
%zind0 = intersect(Ti,1:j0-1);
%if ~isempty(symsetdiff(zind,zind0))
%  'aha'
%  keyboard
%end
DZ = repmat([m],[1 length(zind)]);

for xi=1:prod(DX)
  x = ind2coord(xi,DX);
  Sum = 0;
  for zi=1:prod(DZ)
    z = ind2coord(zi,DZ);
    
    pcond1 = getcondP0([zind xind],[i],[z x],[w1]);
    pcond2 = getcondP0([zind xind],[i],[z x],[w2]);
    Sum = Sum + pcond1 - pcond2;
    
  end
  h = h + .5*abs(Sum);
end
return

function h = etaijXdtree22ub(P,m,n,i,j,y,w1,w2)
global A C T
global LEVELS
h = 0;
LEVELS = repmat({[]},n,1);
Ti = getSubt(T,i);
jj = Ti(Ti>=j);
if isempty(jj)
  return
end
getLevels(T,1,0);
j0 = min(jj);
levi0 = getLevu(i);
levj0 = getLevu(j0);

Ej = getEj(T,j0);
xind = intersect(Ej(:,2),Ti);
xind = intersect(xind,LEVELS{levj0});
DX = repmat([m],[1 length(xind)]);
zind = intersect(Ti,i+1:j0-1);
DZ = repmat([m],[1 length(zind)]);

for xi=1:prod(DX)
  x = ind2coord(xi,DX);
  Sum = 0;
  for zi=1:prod(DZ)
    z = ind2coord(zi,DZ);
    pcond1 = getcondP0([zind xind],[i],[z x],[w1]);
    pcond2 = getcondP0([zind xind],[i],[z x],[w2]);
    Sum = Sum + pcond1 - pcond2;
  end
  h = h + .5*abs(Sum);
end
return



function Pc = getcondP0(jj,ii,x,y)
global P m n
% P = P{ X[jj]=x | X[ii]=y }
Pjoin = getPX(P,m,n,[y x],[ii jj]);
Pmarg = getPX(P,m,n,[y  ],[ii   ]);
Pc = Pjoin/(Pmarg+realmin);
return




% good
function h = etaijXdtree21(P,m,n,i,j,y,w1,w2)
global A C T
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
D1n = repmat([m],[1 n]);

Ei = getEij(T,i,i+1);

for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  px = PI(T,xjn,j:n);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    yw1zx = [y w1 zij xjn];
    ind1 = coord2ind(yw1zx,D1n);
    yw2zx = [y w2 zij xjn];
    ind2 = coord2ind(yw2zx,D1n);

    s = [y w1 zij xjn];
      
    f2 = 1; g1 = 1; g2 = 1;
    for e=1:size(Ei,1)
      u = Ei(e,1);
      v = Ei(e,2);
      if u==i
        g1 = g1 * A(i,v,s(v),w1);
        g2 = g2 * A(i,v,s(v),w2);
      else
        f2 = f2 * A(u,v,s(v),s(u));
      end
    end
    
    pcond1 = px*f2*g1;
    pcond2 = px*f2*g2;
    Sum = Sum + pcond1 - pcond2;
    
  end
  h = h + .5*abs(Sum);
end
return



function h = etaijXdtree2(P,m,n,i,j,y,w1,w2)
global A C T
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
D1n = repmat([m],[1 n]);

Ei = getEij(T,i,i+1);

for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    yw1zx = [y w1 zij xjn];
    ind1 = coord2ind(yw1zx,D1n);
    yw2zx = [y w2 zij xjn];
    ind2 = coord2ind(yw2zx,D1n);

    s = [y w1 zij xjn];
      
    f1 = 1;
    for e=1:size(Ei,1)
      u = Ei(e,1);
      v = Ei(e,2);
      f1 = f1 * Pi(T,v,s);
    end
    
    f2 = 1; g1 = 1; g2 = 1;
    for e=1:size(Ei,1)
      u = Ei(e,1);
      v = Ei(e,2);
      if u==i
        g1 = g1 * A(i,v,s(v),w1);
        g2 = g2 * A(i,v,s(v),w2);
      else
        f2 = f2 * A(u,v,s(v),s(u));
      end
    end
    
    pcond1 = f1*f2*g1;
    pcond2 = f1*f2*g2;
    Sum = Sum + pcond1 - pcond2;
    
  end
  h = h + .5*abs(Sum);
end
return





function h = etaijXdtree1(P,m,n,i,j,y,w1,w2)
global A C T
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
D1n = repmat([m],[1 n]);
Pw1 = getPX(P,m,n,[y w1],1:i);
Pw2 = getPX(P,m,n,[y w2],1:i);
for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    yw1zx = [y w1 zij xjn];
    ind1 = coord2ind(yw1zx,D1n);
    yw2zx = [y w2 zij xjn];
    ind2 = coord2ind(yw2zx,D1n);

    s1 = [y w1 zij xjn];
    s2 = [y w2 zij xjn];
    Ei = getEij(T,i,i+1);
      
    f1 = 1;
    for e=1:size(Ei,1)
      u = Ei(e,1);
      v = Ei(e,2);
      f1 = f1 * Pi(T,v,s1);
    end
    
    f2 = 1; g1 = 1; g2 = 1;
    for e=1:size(Ei,1)
      u = Ei(e,1);
      v = Ei(e,2);
      if u==i
        g1 = g1 * A(u,v,s1(v),s1(u));
        g2 = g2 * A(u,v,s2(v),s2(u));
      else
        f2 = f2 * A(u,v,s1(v),s1(u));
      end
    end
    
    pcond1 = f1*f2*g1;
    pcond2 = f1*f2*g2;
    Sum = Sum + pcond1 - pcond2;
    
  end
  h = h + .5*abs(Sum);
end
return


function p = treecondpr(T,i,yx)
global A
Eij = getEij(T,i,i+1);
p = 1;
for e=1:size(Eij,1)
  u = Eij(e,1);
  v = Eij(e,2);
  p = p * Psubt(T,u,v,yx);
  
  %p1 = Psubt(T,u,v,yx);
  %p2 = A(u,v,yx(v),yx(u)) * Pi(T,v,yx);
  %good = abs(p1-p2)<1e-9;
  %if ~good
  %  'not good tcp'
  %  keyboard
  %end
  
end
return

function p = Psubt(T,u,v,yx)
global A
% Pr{sub-tree starting at v == x | v's parent u == y}
%if isempty(x)
%  p = 1;
%else
p = A(u,v,yx(v),yx(u));
vkids = find(T(v,:));
u1 = v;
for vi=1:length(vkids)
  v1 = vkids(vi);
  p = p * Psubt(T,u1,v1,yx);
end
%end
return

function p = PI(T,x,jj)
global SUBT n
p = 1;
s = zeros(1,n);
s(jj) = x;
while ~isempty(jj)
  v = jj(1);
  SUBT = [];
  p = p*Pi(T,v,s);
  jj = setdiff(jj,SUBT);
end
return


function p = Pi(T,v,x)
global A SUBT
p = 1;
SUBT(end+1) = v;
vkids = find(T(v,:));
for vi=1:length(vkids)
  v1 = vkids(vi);
  p = p * A(v,v1,x(v1),x(v));
  p = p * Pi(T,v1,x);
end
%[nodesi,nodesj]=find(T);
%for e=1:length(nodesi)
%  i = nodesi(e);
%  j = nodesj(e);
%  if (i>=v)
%    p = p * A(i,j,x(j),x(i));
%  end
%end
return

function Eij = getEij(T,i,j)
Eij = [];
[nodesi,nodesj]=find(T);
for e=1:length(nodesi)
  e1 = nodesi(e);
  e2 = nodesj(e);
  if ((e1<=i) & (e2>=j))
    Eij = [Eij;[e1 e2]];
  end
end
return


function P = dtreefillprob(m,n,T,A,C)
% for directed trees
% the first node is always the root
P = zeros(m^n,1);
dims = repmat([m],[1 n]);
[nodesi,nodesj]=find(T);
ne = length(nodesi);
for xind = 1:m^n
  x = ind2coord(xind,dims);
  p = C(x(1));
  for e=1:length(nodesi)
    i = nodesi(e);
    j = nodesj(e);
    p = p * A(i,j,x(j),x(i));
  end
  P(xind) = p;
end
return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function strucsearch
global n m P T A C th0
global mrat mT mA mC
global longA
longA = 0;
mrat = 0;

while 1
  T = randdtree(n);
  maxd = fanout(T);
  for b = .01:.1:.99
    %for b1 = .0:.1:1
    for b1 = 0
      A = eye(2);
      A(1,1) = b;
      A(2,1) = 1-b;
      A(1,2) = b1;
      A(2,2) = 1-b1;
      A = A + eps;
      A = pnormdim(A,1);
      C = ones(m,1)/m;
      P = dtreefillprob(m,n,T,A,C);
      th = getTH(A);
      
      [hh,Ht,H] = gethhn(P,m,n);
      maxhr = 0;
      for i=1:n-1
        for j=i+1:n
          i0=geti0(T,i,j);
          di0 = sum(T(i0,:));
          %Bij = th^(j-i)*maxd;
          if maxd<=1
            mdep = j-i;
          else
            %mdep = ceil(log(j-i)/log(maxd));
            %N = j-i+1;
            %mdep = floor(log(N*(maxd-1)+1)/log(maxd)-1);
            N = j-i;
            %mdep = floor(log(N*(maxd-1)+1)/log(maxd)-1);
            %mdep = minlev(j,maxd)-minlev(i,maxd);
            li = minlev(i,maxd);
            lj = floor(log(N*(maxd-1)/maxd^li + 1)/log(maxd)) + li -1;
            mdep = lj - li;
          end
          Bij = di0*th^mdep;
          if min(Bij,hh(i,j)) > 1e-13
            % otherwise, all bets are off -- numerics too sketchy
            rat = hh(i,j)/Bij;
            maxhr = max(maxhr,rat);
            if rat > 1.01
              fprintf('i=%d  j=%d    mdep = %d    Bij = %1.5f rat = %1.5f\n',i,j,mdep,Bij,hh(i,j)/Bij);
            end
          end
        end
        %maxhr = max(maxhr,hh(i,i+1));
      end
      r = maxhr;
      
      if r > mrat
        mrat = r;
        mT = T;
        mA = A;
        mC = C;
        save TREDAT mT mA mC mrat
      end
      fprintf('strucsearch n=%d   k=%d      MRAT = %1.9f\n',n,maxd,mrat);    
    end
  end
end
return


function doplot
global n m P T A C th0
global mrat mT mA mC
global longA
longA = 0;
mrat = 0;
bb = [];
rr = [];

%T = randdtree(n);

%fstr = '1<2<5; 2<6; 1<3<7; 3<8; 1<4<9; 4<10';
%fstr = '1 < 2 ;1 < 3 ;2 < 4 ;2 < 5 ;4 < 6 ;6 < 7 ';

fstr = '1<2;1<3;2<4;2<5;3<6;3<7;7<8'

fstr='1 < 2;2 < 3;2 < 4;3 < 5;3 < 6;5 < 8;5 < 9;6 < 10;6 < 11;11 < 12;4 < 7'


T = readtree(fstr);

%T = readtree(fstr);
%n = size(T,1);



lev = 4;
T = bintree(lev)
n = size(T,1);


while 1
  b = rand;
  A = eye(2);
  A(1,1) = b;
  A(2,1) = 1-b;
  C = ones(m,1)/m;
  
  P = dtreefillprob(m,n,T,A,C);
  maxd = fanout(T);
  th = getTH(A);

  if 1
    i = 1;
    %j = 5;
    j = 2^(lev-1);
    hij = mxetaij(P,m,n,i,j);
    N = j-i+1;
    mdep = floor(log(N*(maxd-1)+1)/log(maxd)-1);
    Bij = maxd*th^mdep;
    r = hij/Bij;
  else
    [hh,Ht,H] = gethhn(P,m,n);
    maxhr = 0;
    for i=1:n-1
      for j=i+1:n
        %Bij = th^(j-i)*maxd;
        if maxd==1
          mdep = j-i;
        else
          %mdep = ceil(log(j-i)/log(maxd));
          N = j-i+1;
          mdep = floor(log(N*(maxd-1)+1)/log(maxd)-1);
        end
        Bij = maxd*th^mdep;
        if min(Bij,hh(i,j)) > 1e-13
          % otherwise, all bets are off -- numerics too sketchy
          rat = hh(i,j)/Bij;
          maxhr = max(maxhr,rat);
          if rat > .9
            fprintf('i=%d  j=%d    mdep = %d    Bij = %1.5f rat = %1.5f\n',i,j,mdep,Bij,hh(i,j)/Bij);
          end
        end
      end
      %maxhr = max(maxhr,hh(i,i+1));
    end
    
    r = maxhr;
    
  end
  
  bb(end+1) = b;
  rr(end+1) = r;
  plot(bb,rr,'.');
  drawnow
end
return


function trespc
global n m P T A C th0
global mrat mT mA mC
global longA
longA = 1;
mrat = 0;

global lev
lev = 4;


adim = m^2;
xdim = adim;

LB = zeros(xdim,1)+1e-10;
UB = ones(xdim,1);

% normalization constraints
Aeq = zeros(m,adim);
for i=1:m
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m,1);
while 1
  th0 = rand;
  %T = randdtree(n);
  %fstr = '1<2<5; 2<6; 1<3<7; 3<8; 1<4<9; 4<10';
  %T = readtree(fstr);
  T = bintree(lev)
  n = size(T,1);
  if longA
    [nodesi,nodesj]=find(T);
    ne = length(nodesi);
    Aeq = zeros(m*T,m^2*ne);
    for i=1:m*ne
      for j=1:m
        Aeq(i,m*(i-1)+j) = 1;
      end
    end
    Beq = ones(m*ne,1);
    [Ain,Bin] = constrict(m,th0,ne);
    LB = zeros(m^2*ne,1)+1e-10;
    UB = ones(m^2*ne,1);
    A = pnormdim(rand(m,m,ne),1);
  else
    [Ain,Bin] = constrict(m,th0,1);
    A = pnormdim(rand(m,m),1);
    %C = pnormdim(rand(m,1),1);
  end
  C = pnormdim(ones(m,1),1);
  x = A(:);
  [x,fval] = fmincon(@spcopt,x,Ain,Bin,Aeq ,Beq  ,LB,UB);
end
return

function y = spcopt(x)
global n m P T C th0
global mrat mT mA mC
global longA
global lev
if longA
  [nodesi,nodesj]=find(T);
  ne = length(nodesi);
  A = reshape(x,[m m ne]);
else
  A = reshape(x,[m m]);
end
P = dtreefillprob(m,n,T,A,C);
%[hh,Ht,H] = gethhn(P,m,n);
maxd = fanout(T);
th = getTH(A);
th = max(th);
%th = max(th0,th);

i = 1;
%j = 5;
j = 2^(lev-1);
N = j-i+1;
hij = mxetaij(P,m,n,i,j);
mdep = floor(log(N*(maxd-1)+1)/log(maxd)-1);
Bij = maxd*th^mdep;

r = 0;
if min(Bij,hij) > 1e-13
  % otherwise, all bets are off -- numerics too sketchy
  r = hij/Bij;
  %maxhr = max(maxhr,rat);
  %if rat > 10
  %  fprintf('i=%d  j=%d    mdep = %d    Bij = %1.5f rat = %1.5f\n',i,j,mdep,Bij,hh(i,j)/Bij);
  %end
end

y = -r;

if r > mrat
  mrat = r;
  mT = T;
  mA = A;
  mC = C;
  save TREDAT mT mA mC mrat
end
fprintf('spcopt n=%d   k=%d      MRAT = %1.9f\n',n,maxd,mrat);

return




function trexplore
global n m P T A C th0
global mrat mT mA mC
global longA
longA = 0;
mrat = 0;

adim = m^2;
xdim = adim;

LB = zeros(xdim,1)+1e-10;
UB = ones(xdim,1);

% normalization constraints
Aeq = zeros(m,adim);
for i=1:m
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m,1);
while 1
  th0 = rand;
  T = randdtree(n);
  %fstr='1<2;1<3;2<4;2<5;4<7;5<8;5<9;3<6;6<10;6<11;11<12;';
  
  % this is the problem one ***
  %fstr = '1<2; 2<3; 2<4<7; 3<5<8; 3<5<9; 3<6<10; 6<11<12'; 
  % try to simplify it
  %fstr = '1<2; 2<3; 2<4<7; 3<5<8; 3<5<9; 3<6<10; 6<11'; 
  %fstr = '1<2; 1<3<6; 2<4<7; 2<4<8; 2<5<9 ; 5<10'; % still 'good'
  %fstr = '1<2; 1<3<6; 3<7; 2<4<8; 2<4<9; 2<5<10; 5<11';
  
  %fstr = '1<2; 1<3<6; 2<4<7; 2<4<8; 2<5<9 ; 5<10; 7<11; 8<12; 9<13; 10<14';
  
  %fstr = '1<2; 1<3<6; 2<4<7; 2<4<8; 2<5<9 ; 5<10; 6<11'; % still 'good'
  
  %fstr = '1<2; 1<3<6; 2<4<7; 2<4<8; 2<5<9'; %'nope'
  %fstr = '1<2; 1<3; 2<4<6; 2<4<7; 2<5<8 ; 5< 9'; %nope -- need 6
  
  
  %T = readtree(fstr);
  %n = size(T,1);
  
  if longA
    nothing_yet
  else
    [Ain,Bin] = constrict(m,th0,1);
    A = pnormdim(rand(m,m),1);
    C = pnormdim(rand(m,1),1);
  end
  x = A(:);
  [x,fval] = fmincon(@treeopt,x,Ain,Bin,Aeq ,Beq  ,LB,UB);
end
return

function y = treeopt(x)
global n m P T C th0
global mrat mT mA mC
global longA
A = reshape(x,[m m]);
P = dtreefillprob(m,n,T,A,C);
[hh,Ht,H] = gethhn(P,m,n);
maxd = fanout(T);
th = getTH(A);
%th = max(th0,th);

maxhr = 0;
for i=1:n-1
  for j=i+1:n
    %Bij = th^(j-i)*maxd;
    if maxd==1
      mdep = j-i;
    else
      %mdep = ceil(log(j-i)/log(maxd));
      N = j-i;
      %mdep = floor(log(N*(maxd-1)+1)/log(maxd)-1);
      %mdep = minlev(j,maxd)-minlev(i,maxd);
      li = minlev(i,maxd);
      lj = floor(log(N*(maxd-1)/maxd^li + 1)/log(maxd)) + li -1;
      mdep = lj - li;
    end
    Bij = th^mdep;
    if min(Bij,hh(i,j)) > 1e-13
      % otherwise, all bets are off -- numerics too sketchy
      rat = hh(i,j)/Bij;
      maxhr = max(maxhr,rat);
      if rat > 1.01
        fprintf('i=%d  j=%d    mdep = %d    Bij = %1.5f rat = %1.5f\n',i,j,mdep,Bij,hh(i,j)/Bij);
      end
    end
  end
  %maxhr = max(maxhr,hh(i,i+1));
end

r = maxhr;


y = -r;

if r > mrat
  mrat = r;
  mT = T;
  mA = A;
  mC = C;
  save TREDAT mT mA mC mrat
end
fprintf('trexplore n=%d   k=%d      MRAT = %1.9f\n',n,maxd,mrat);

return




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

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

function chainxplore
global n m P T pi
global mrat mT mpsi
global longpsi
longpsi = 1;
mrat = 0;
T = chaingraph(n);
while 1
  %pi = randperm(n);
  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(:));
  x = psi(:);
  [x,fval] = fmincon(@chainopt,x,[] ,[]  ,[] ,[]  ,0*x+1,[]);
end
return

function y = chainopt(x)
global n m P T pi
global longpsi
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
P = treefillprob(m,n,T,psi);
%P = reorderP(P,m,n,pi);
[R,rr] = psi2r(psi);
[b,A,C,dd] = is_markov(P,m,n); 
[h,hh] = stricthh(A);
%th = strict(A);
%th = strict(A(:,:,2));
%r = log(1-th)./(1-R);

%B = (R+1)/2;
BB = (rr-1)./(rr+1);
%[hh,Ht,H] = gethhn(P,m,n);
%r = H/B;
r = max(hh./rr);

y = -r;

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

return




function [Aina,Bina] = constrict(m,h,T)
% stricture constraints

if ~exist('T','var')
  T = 1;
end

SIGNS = 2*list_all_ind(repmat([2],[1 m]))-3;
Aina = zeros(size(SIGNS,1)*m*(m-1)/2*T,m^2*T);
rind = 0;

for t=1:T
  for i1=1:m
    for i2=i1+1:m
      for s = 1:size(SIGNS,1)
        rind = rind+1;
        ss = SIGNS(s,:);
        for j=1:m
          ind1 = coord2ind([j i1 t],[m m T]);
          ind2 = coord2ind([j i2 t],[m m T]);
          Aina(rind,ind1) =  ss(j);
          Aina(rind,ind2) = -ss(j);
        end
      end
    end
  end
end
Bina = ones(size(Aina,1),1)*h;
Bina = Bina * 2;
return




function hh = getTH(A)
[n,n,T] = size(A);
hh = zeros(1,T);
for t=1:T
  h = 0;
  for i=1:n
    for j=i+1:n
      x = A(:,i,t);
      y = A(:,j,t);
      d = .5*sum(abs(x-y));
      if d>h
        h=d;
      end
    end
  end
  hh(t) = h;
end
return



function h = etaijXdtree00(P,m,n,i,j,y,w1,w2)
global A C T
h = 0;
Djn = repmat([m],[1 n-j+1]);
Dij = repmat([m],[1 j-i-1]);
D1n = repmat([m],[1 n]);
Pw1 = getPX(P,m,n,[y w1],1:i);
Pw2 = getPX(P,m,n,[y w2],1:i);
for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    yw1zx = [y w1 zij xjn];
    ind1 = coord2ind(yw1zx,D1n);
    yw2zx = [y w2 zij xjn];
    ind2 = coord2ind(yw2zx,D1n);

    if 0
      %Sum = Sum + P(ind1)/Pw1 - P(ind2)/Pw2;

      
      s1 = [y w1 zij xjn];
      s2 = [y w2 zij xjn];
      %pcond1 = treecondpr(T,i,[y w1 zij xjn]);
      %pcond2 = treecondpr(T,i,[y w2 zij xjn]);

      Ei = getEij(T,i,i+1);
      
      f1 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        v = Ei(e,2);
        f1 = f1 * A(u,v,s1(v),s1(u))* Pi(T,v,s1);
      end
      pcond1 = f1;
      
      f2 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        v = Ei(e,2);
        f2 = f2 * A(u,v,s2(v),s2(u))* Pi(T,v,s2);
      end
      pcond2 = f2;
      
      
      g1 = abs(pcond1-P(ind1)/Pw1)<1e-9;
      g2 = abs(pcond2-P(ind2)/Pw2)<1e-9;
      good = g1 & g2;
      if ~good
        'not good'
        keyboard
      end
      Sum = Sum + pcond1 - pcond2;
    end

    
    if 1
      s1 = [y w1 zij xjn];
      s2 = [y w2 zij xjn];
      Ei = getEij(T,i,i+1);
      
      f1 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        v = Ei(e,2);
        f1 = f1 * Pi(T,v,s1);
      end
      
      f2 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        v = Ei(e,2);
        if u~=i
          f2 = f2 * A(u,v,s1(v),s1(u));
        end
      end
      g1 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        v = Ei(e,2);
        if u==i
          g1 = g1 * A(u,v,s1(v),s1(u));
        end
      end
      g2 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        v = Ei(e,2);
        if u==i
          g2 = g2 * A(u,v,s2(v),s2(u));
        end
      end
      
      
      %g1 = abs(pcond1-P(ind1)/Pw1)<1e-9;
      %g2 = abs(pcond2-P(ind2)/Pw2)<1e-9;
      %good = g1 & g2;
      %if ~good
      %  'not good'
      %  keyboard
      %end
      
      pcond1 = f1*f2*g1;
      pcond2 = f1*f2*g2;
      Sum = Sum + pcond1 - pcond2;
      %Sum = Sum + f1*neqi*eqi;
    end
    
    
    
    if 0
      Ei = getEij(T,i,i+1);
      f1 = 1;
      s1 = [y w1 zij xjn];
      s2 = [y w2 zij xjn];
      for e=1:size(Ei,1)
        v = Ei(e,2);
        f1 = f1 * Pi(T,v,s1);
        %good = abs(Pi(T,v,s1)-Pi(T,v,s2))<1e-9;   if ~good, 'not good f1',  keyboard;   end
      end
      f2 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        if u~=i
          v = Ei(e,2);
          f2 = f2 * A(u,v,s1(v),s1(u));
          good = abs(A(u,v,s1(v),s1(u))-A(u,v,s2(v),s2(u)))<1e-9;   if ~good, 'not good f2',  keyboard; end
        end
      end
      f3 = 1;
      for e=1:size(Ei,1)
        u = Ei(e,1);
        if u==i
          v = Ei(e,2);
          f3 = f3 * (A(i,v,s1(v),w1) - A(i,v,s2(v),w2));
        end
      end
      Sum = Sum + f1*f2*f3;
    end
    
  end
  h = h + .5*abs(Sum);
end
return


function sT = getSubt(T,v)
% a list of v's descendents in T
% includes v itself
global SUBT
SUBT = [];
Pi0(T,v);
sT = SUBT;
return

function Pi0(T,v)
global SUBT
SUBT(end+1) = v;
vkids = find(T(v,:));
for vi=1:length(vkids)
  v1 = vkids(vi);
  Pi0(T,v1);
end

function getLevels(T,u,ell)
global LEVELS
ell = ell+1;
LEVELS{ell}(end+1) = u;
ukids = find(T(u,:));
for vi=1:length(ukids)
  v = ukids(vi);
  getLevels(T,v,ell);
end
return

function [ell0,LEVu] = getLevu(u)
global LEVELS
for ell=1:length(LEVELS)
  if ismember(u,LEVELS{ell})
    LEVu = LEVELS{ell};
    ell0 = ell;
    return
  end
end
return


function Ej = getEj(T,j)
Ej = [];
[nodesi,nodesj]=find(T);
for e=1:length(nodesi)
  u = nodesi(e);
  v = nodesj(e);
  if ((u<j) & (v>=j))
  %if ((u<j))
    Ej = [Ej;[u v]];
  end
end
return