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

[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});
jj = intersect(Ti,LEVELS{ell});
LOC = treetensor(T,A,ii,jj);
TTENS = mytprod(LOC,TTENS);
end

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

Ej = getEj(T,j0);
IC = intersect(Ej(:,2),Ti);
DC = repmat([m],[1 length(IC)]);
IZ = intersect(Ti,i+1:j0-1);
DZ = repmat([m],[1 length(IZ)]);

LEVJ0 = intersect(Ti,LEVELS{levj0});
C0 = intersect(IC,LEVJ0);
C1 = setdiff(IC,C0);
Z0 = setdiff(LEVJ0,C0);
Z1 = setdiff(IZ,Z0);

DC0 = repmat([m],[1 length(C0)]);
DC1 = repmat([m],[1 length(C1)]);
DZ0 = repmat([m],[1 length(Z0)]);
DZ1 = repmat([m],[1 length(Z1)]);
DJ0 = repmat([m],[1 length(LEVJ0)]);

B = zeros(prod([DC0 DC1]),prod(DJ0));
if isempty(Z0)
  % implies empty(C1)
  B = eye(prod(DJ0));
elseif isempty(C1)
  % sum out the Z0
  for ci0=1:prod(DC0)
    xC0 = ind2coord(ci0,DC0);
    for ci2=1:prod(DC0)
      xC2 = ind2coord(ci2,DC0);
      for zi0=1:prod(DZ0)
        xZ0 = ind2coord(zi0,DZ0);
        B(ci0,coord2ind([xZ0 xC2],DJ0)) = (ci0==ci2);
      end
    end
  end
else
  for ci=1:prod(DC)
    xC = ind2coord(ci,DC);
    for ji=1:prod(DJ0)
      xJ = ind2coord(ji,DJ0);
      xC0 = xC(1:length(C0));
      xC1 = xC(length(C0)+1:end);
      xZ0 = xJ(1:length(Z0));
      xC2 = xJ(length(Z0)+1:end);
      
      pC1Z0 = getcondP0(C1,Z0,xC1,xZ0);
      B(ci,ji) = pC1Z0 * veq(xC0,xC2);
    end
  end  
end

B
F = B*H;

h = TV(F);
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



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


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