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

levi0 = getLevu(i);
levj0 = getLevu(j0);

TTENS.ii = [i];
TTENS.jj = [i];
TTENS.A = 1;
for ell = levi0+1:levj0
  ii = intersect(Ti,LEVELS{ell-1});
  %ii = intersect(ii,i:j0);
  jj = intersect(Ti,LEVELS{ell});
  %jj = jj(jj<j0);
  if ~isempty(jj)
    %jj = intersect(jj,i:j0);
    LOC = treetensor(T,A,ii,jj);
    TTENS = mytprod(LOC,TTENS);
  end
end

% find childless k<j0 
nokidj = [];
for ki = 1:length(LEVELS{levj0})
  u = LEVELS{levj0}(ki);
  if (u<j0) & isempty(find(T(u,:)))
    nokidj(end+1) = u;
  end
end

% close the zx gap
Ej = getVj(T,j0);
% push TTENS forward so it's a distribution on X_j0^n
Xinds = j0:n;
Xintr = Ej(:,2)';
if ~all(ismember(Xintr,TTENS.jj))
  LOC = treetensor(T,A,TTENS.jj,Xintr);
  TTENS = mytprod(LOC,TTENS);
end

% collapse TTENS along childless k
if ~isempty(nokidj)
  jj1 = setdiff(TTENS.jj,nokidj);
  dimjnk = repmat(m,[1 length(nokidj)]);
  dimjj1 = repmat(m,[1 length(jj1)]);
  dimjj = repmat(m,[1 length(jj)]);
  dimii = repmat(m,[1 length(TTENS.ii)]);
  
  A1 = zeros(m^length(jj1),m^length(TTENS.ii));
  for xi=1:prod(dimii)
    s = zeros(1, length(TTENS.jj));
    for yi1=1:prod(dimjj1)
      y1 = ind2coord(yi1,dimjj1);
      s(jj1) = y1;
      for yi2=1:prod(dimjnk)
        y2 = ind2coord(yi2,dimjnk);
        s(nokidj) = y2;
        y = s(TTENS.jj);
        yi = coord2ind(y,dimjj);
        A1(yi1,xi) = A1(yi1,xi) + TTENS.A(yi,xi);
      end
    end
  end
  TTENS.jj = jj1;
  TTENS.A = A1;
end


A12 = TTENS.A(:,w1) - TTENS.A(:,w2);


%Xinds = intersect(j0:n,Ti);
Djn = repmat([m],[1 length(Xinds)]);
s0 = zeros(1,n);
for xi=1:prod(Djn)
  xjn = ind2coord(xi,Djn);
  s0(Xinds) = xjn;
  p = 1;
  for ei = 1:size(Ej,1)
    %(u,v) = Ej(ei,:);
    v = Ej(ei,2);
    p = p * Pi(T,v,s0,Ti);
  end
  %tcoor = s0(Ej(:,2)');
  tcoor = s0(TTENS.jj);
  tind = coord2ind(tcoor,repmat([m],[1 length(tcoor)]));
  p = p * A12(tind);

  h = h + .5*abs(p);
end
return


function Ej = getVj(T,j)
Ej = [];
[nodesi,nodesj]=find(T);
for e=1:length(nodesi)
  e1 = nodesi(e);
  e2 = nodesj(e);
  if ((e1<j) & (e2>=j))
    Ej = [Ej;[e1 e2]];
    %Vj(end+1) = e2;
  end
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


function vec = grabvec(A,u,v,w)
global n m
dims = repmat(1,[1 n]);
dims(v) = m;
vec = reshape(A(u,v,:,w),dims);
return




function Z = mytprod(X,Y)
global m

Z.jj = X.jj;
Z.ii = Y.ii;
Z.A = X.A * Y.A;

return