function h = treesubij(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);

% which nodes are in the ij subtree?
if 0
  SubTij = [];
  for vi = 1:length(Ti)
    v = Ti(vi);
    lev = getLevu(v);
    if lev <= levj0
      SubTij(end+1) = v;
    end
  end
  SubTij = SubTij(:)';
else
  % can't just cut those off
  SubTij = Ti;
end

%Xinds = SubTij(SubTij>=j0);
%Zinds = SubTij((i<SubTij)&(SubTij< j0));
%Zinds = setdiff([i+1:n],Xinds);

Xinds = intersect(j:n,SubTij);
Zinds = intersect(i+1:j-1,SubTij);

%Zinds = SubTij(:)';
%Zinds = Zinds(2:end);

Djn = repmat([m],[1 length(Xinds)]);
Dij = repmat([m],[1 length(Zinds)]);

s1 = zeros(1,n);
s2 = zeros(1,n);
for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  s1(Xinds) = xjn;
  s2(Xinds) = xjn;
  Sum = 0;
  for zi=1:prod(Dij)
    zij = ind2coord(zi,Dij);
    s1(Zinds) = zij;
    s2(Zinds) = zij;
    s1(i) = w1;
    s2(i) = w2;
    %s1([i Zinds]) = [w1 zij];
    %s2([i Zinds]) = [w2 zij];

    p1 = Pi(T,i,s1,SubTij);
    p2 = Pi(T,i,s2,SubTij);
    Sum = Sum + p1 - p2;
    %Sum = Sum + abs(p1 - p2);
  end
  h = h + .5*abs(Sum);
end
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