function h = treesubij1(i,j,w1,w2)
global A C T m n
% computes hh(i,j) for a tree

% find j0, where
% j0 is the minimum descendant of i s.t. j0>=j

[n,n] = size(T);
h = 0;
global LEVELS
LEVELS = repmat({[]},n,1);
Ti = getSubt(T,i);
jj = Ti(Ti>=j);
if isempty(jj)
  return
end
getLevels(T,1,0);
j0 = min(jj);
levj0 = getLevu(j0);

Cut = getCut(T,j);
C = Cut(:,2);

Xinds = intersect(j0:n,Ti);
Cinds = intersect(C,Ti);
Cinds = Xinds;
Zinds = intersect(i+1:j0-1,Ti);

DC = repmat([m],[1 length(Cinds)]);

for xi=1:prod(DC)
  xc = ind2coord(xi,DC);
  Z = zeta(i,Cut,xc,Cinds,Zinds,w1,w2);
  h = h + .5*abs(Z);
end
return


function Z = zeta(i,Cut,xc,cind,zind,w1,w2)
global m n A C T P


if isempty(zind)
  P1  = getcondP(cind,[i],xc,[w1]);
  P2  = getcondP(cind,[i],xc,[w2]);
  Z = P1 - P2;
else
  Dz = repmat([m],[1 length(zind)]);
  Z = 0;
  for zi=1:prod(Dz)
    xz = ind2coord(zi,Dz);
    Pcz = getcondP(cind,zind,xc,xz);
    P1  = getcondP(zind,[i],xz,[w1]);
    P2  = getcondP(zind,[i],xz,[w2]);
    Z = Z + Pcz*(P1 - P2);
  end
end
return


function Pc = getcondP(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



function p = Pi(T,v,x,S)
global A
p = 1;
vkids = find(T(v,:));
for vi=1:length(vkids)
  v1 = vkids(vi);
  if ismember(v1,S)
    p = p * A(v,v1,x(v1),x(v));
    p = p * Pi(T,v1,x,S);
  end
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 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 pseq = parentseq(T,u,r)
v = u;
el = 0;
pseq = [u];
while v~=r
  el = el+1;
  rents = find(T(:,v));
  v = rents;
  pseq = [v pseq];
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 C = getCut(T,j)
C = [];
[nodesi,nodesj]=find(T);
for e=1:length(nodesi)
  u = nodesi(e);
  v = nodesj(e);
  if ((u<j)&(v>=j))
    C = [C;[u v]];
  end
end
return