function T1 = marginalize_max_table(T2, onto)
  
node_sizes = sparse(1, max(T2.domain));
node_sizes(T2.domain) = T2.sizes;
T1 = mk_max_table(onto, node_sizes(onto));
sum_over = setdiff(T2.domain, T1.domain);
summed_sizes = node_sizes(sum_over);
num_summed = full(prod(summed_sizes));

sum_map = find_equiv_posns(sum_over, T2.domain);
keep_map = find_equiv_posns(T1.domain, T2.domain);
iv = zeros(1,length(T2.domain));

for i1 = 1:length(T1.T)
  keep_iv = ind2subv(T1.sizes, i1);
  iv(keep_map) = keep_iv;
  score = zeros(1, num_summed);
  for j = 1:num_summed
    sum_iv = ind2subv(summed_sizes, j);
    iv(sum_map) = sum_iv;    
    i2 = subv2ind(T2.sizes, iv);
    score(j) = T2.T(i2);
  end
  [T1.T(i1), T1.arg(i1)] = max(score);
  %T1.T(i1) = sum(score); % gives the same results as marginalize_table
end

