wolffd@0: function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
wolffd@0: % MARGINALIZE_POT Marginalize a mpot 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: 
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: [logp, mu, Sigma] = marginalize_gaussian(bigpot.logp, bigpot.mu, bigpot.Sigma, ...
wolffd@0: 					 keep, sum_over, node_sizes);
wolffd@0: smallpot = mpot(keep, node_sizes(keep), logp, mu, Sigma);
wolffd@0: 
wolffd@0: %%%%%%
wolffd@0: 
wolffd@0: function [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns)
wolffd@0: % MARGINALIZE_GAUSSIAN Compute Pr(X) from Pr(X,Y) where X and Y are jointly Gaussian.
wolffd@0: % [logpX, muX, SXX] = marginalize_gaussian(logp, mu, Sigma, X, Y, ns)
wolffd@0: %
wolffd@0: % sizes(i) is the size of the i'th block in domain.
wolffd@0:  
wolffd@0: [muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns);
wolffd@0: logpX = logp; % Lauritzen (1996) p161