function h = treesubtensub(i,j,w1,w2)
global A C T m n
% computes hh(i,j) for a tree
% upper bound
% lama la'avod kashe?

% 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;
if levj0==levi0+1
  % no "blocking" levels
  ii = i;
  jj = find(T(i,:));
  %jj = setdiff(jj,leaves(T));
  TTENS = treetensor(T,A,ii,jj);
else
  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
end

A12 = TTENS.A(:,w1) - TTENS.A(:,w2);
h = TV(A12);
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