wolffd@0: function smallpot = marginalize_pot(bigpot, keep)
wolffd@0: % MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
wolffd@0: % smallpot = marginalize_pot(bigpot, keep)
wolffd@0: 
wolffd@0: sumover = mysetdiff(bigpot.domain, keep);
wolffd@0: cdom = myunion(bigpot.cheaddom, bigpot.ctaildom);
wolffd@0: csumover = myintersect(sumover, bigpot.cheaddom);
wolffd@0: dsumover = myintersect(sumover, bigpot.ddom);
wolffd@0: 
wolffd@0: dkeep = myintersect(keep, bigpot.ddom);
wolffd@0: ckeep = myintersect(keep, bigpot.cheaddom);
wolffd@0: cheaddom = myintersect(keep, bigpot.cheaddom);
wolffd@0: 
wolffd@0: assert(isempty(myintersect(csumover,bigpot.ctaildom)));
wolffd@0: ns = zeros(1, max(bigpot.domain));
wolffd@0: ns(bigpot.ddom) = bigpot.dsizes;
wolffd@0: ns(bigpot.cheaddom) = bigpot.cheadsizes;
wolffd@0: ns(bigpot.ctaildom) = bigpot.ctailsizes;
wolffd@0: 
wolffd@0: 
wolffd@0: if sum(ns(csumover)) > 0
wolffd@0:     for i=1:bigpot.dsize
wolffd@0:       bigpot.scgpotc{i} = marginalize_pot(bigpot.scgpotc{i}, ckeep, csumover, ns);
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: if (isequal(csumover, cheaddom))
wolffd@0:     bigpot.ctaildom = [];
wolffd@0: end
wolffd@0: % If we are not marginalizing over any discrete nodes, we are done.
wolffd@0: if prod(ns(dsumover))==1
wolffd@0:   smallpot = scgpot(dkeep, cheaddom, bigpot.ctaildom, ns, bigpot.scgpotc);
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: if (~isempty(bigpot.ctaildom))
wolffd@0:     assert(0);
wolffd@0:     return;
wolffd@0: end
wolffd@0: 
wolffd@0: I = prod(ns(dkeep));
wolffd@0: J = prod(ns(dsumover));
wolffd@0: C = sum(ns(ckeep));   
wolffd@0: sum_map = find_equiv_posns(dsumover, bigpot.ddom);
wolffd@0: keep_map = find_equiv_posns(dkeep, bigpot.ddom);
wolffd@0: iv = zeros(1, length(bigpot.ddom)); % index vector
wolffd@0: 
wolffd@0: p1 = zeros(I,J);
wolffd@0: A1 = zeros(C,J,I);
wolffd@0: C1 = zeros(C,C,J,I);
wolffd@0: for i=1:I
wolffd@0:   keep_iv = ind2subv(ns(dkeep), i);
wolffd@0:   iv(keep_map) = keep_iv;
wolffd@0:   for j=1:J
wolffd@0:     sum_iv = ind2subv(ns(dsumover), j);
wolffd@0:     iv(sum_map) = sum_iv;
wolffd@0:     k = subv2ind(ns(bigpot.ddom), iv);
wolffd@0:     pot = struct(bigpot.scgpotc{k}); % violate object privacy
wolffd@0:     p1(i,j) = pot.p;
wolffd@0:     if C > 0 % so mu1 and Sigma1 are non-empty
wolffd@0:       A1(:,j,i) = pot.A;
wolffd@0:       C1(:,:,j,i) = pot.C;
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % Collapse the mixture of Gaussians
wolffd@0: coef = mk_stochastic(p1); % coef must be convex combination
wolffd@0: %keyboard
wolffd@0: p2 = sum(p1,2);
wolffd@0: if (all(p2 == 0))
wolffd@0:     p2 = p2 + (p2==0)*eps;
wolffd@0: end
wolffd@0: A = [];
wolffd@0: S = [];
wolffd@0: 
wolffd@0: pot = cell(1,I);
wolffd@0: ctailsize = sum(ns(bigpot.ctaildom));
wolffd@0: tB = zeros(C, ctailsize);
wolffd@0: for i=1:I
wolffd@0:   if C > 0
wolffd@0:     [A, S] = collapse_mog(A1(:,:,i), C1(:,:,:,i), coef(i,:));
wolffd@0:   end
wolffd@0:   p = p2(i);
wolffd@0:   pot{i} = scgcpot(C, ctailsize, p, A, tB, S);
wolffd@0: end
wolffd@0: 
wolffd@0: smallpot = scgpot(dkeep, ckeep, bigpot.ctaildom, ns, pot);
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: