wolffd@0: function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
wolffd@0: % MARGINALIZE_POT Marginalize a cpot onto a smaller domain.
wolffd@0: % smallpot = marginalize_pot(bigpot, keep, maximize, useC)
wolffd@0: %
wolffd@0: % The maximize argument is ignored - maxing out a Gaussian is the same as summing it out,
wolffd@0: % since the mode and mean are equal.
wolffd@0: % The useC argument is ignored.
wolffd@0: 
wolffd@0: node_sizes = sparse(1, max(bigpot.domain));
wolffd@0: node_sizes(bigpot.domain) = bigpot.sizes;
wolffd@0: sum_over = mysetdiff(bigpot.domain, keep);
wolffd@0: 
wolffd@0: if sum(node_sizes(sum_over))==0 % isempty(sum_over)
wolffd@0:   %smallpot = bigpot;
wolffd@0:   smallpot = cpot(keep, node_sizes(keep), bigpot.g, bigpot.h, bigpot.K);
wolffd@0: else
wolffd@0:   [h1, h2, K11, K12, K21, K22] = partition_matrix_vec(bigpot.h, bigpot.K, sum_over, keep, node_sizes);
wolffd@0:   n = length(h1);
wolffd@0:   K11inv = inv(K11);
wolffd@0:   g = bigpot.g + 0.5*(n*log(2*pi) - log(det(K11)) + h1'*K11inv*h1);
wolffd@0:   if length(h2) > 0 % ~isempty(keep) % we are are actually keeping something
wolffd@0:     A = K21*K11inv;
wolffd@0:     h = h2 - A*h1;
wolffd@0:     K = K22 - A*K12;
wolffd@0:   else
wolffd@0:     h = [];
wolffd@0:     K = [];
wolffd@0:   end
wolffd@0:   smallpot = cpot(keep, node_sizes(keep), g, h, K);
wolffd@0: end
wolffd@0: