function th = thai(P,m,n,i)
th = 0;
dims = repmat([m],[1 i]);
for xi1 = 1:m^(i)
  X1 = ind2coord(xi1,dims);
  P1 = getcondP(P,m,n,[i+1],[1:i],[X1]);
  for xi2 = xi1+1:m^(i)
    X2 = ind2coord(xi2,dims);
    P2 = getcondP(P,m,n,[i+1],[1:i],[X2]);
    th12 = .5*sum(abs(P1-P2));
    th = max(th12,th);
  end
end
return