wolffd@0: function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
wolffd@0: % MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
wolffd@0: % smallpot = marginalize_pot(bigpot, keep, maximize, useC)
wolffd@0: %
wolffd@0: % If maximize = 1, we raise an error.
wolffd@0: % useC is ignored.
wolffd@0: 
wolffd@0: if nargin < 3, maximize = 0; end
wolffd@0: assert(~maximize);
wolffd@0: 
wolffd@0: 
wolffd@0: sumover = mysetdiff(bigpot.domain, keep);
wolffd@0: csumover = myintersect(sumover, bigpot.cdom);
wolffd@0: dsumover = myintersect(sumover, bigpot.ddom);
wolffd@0: dkeep = myintersect(keep, bigpot.ddom);
wolffd@0: ckeep = myintersect(keep, bigpot.cdom);
wolffd@0: %ns = sparse(1, max(bigpot.domain)); % must be full, so I is an integer
wolffd@0: ns = zeros(1, max(bigpot.domain));
wolffd@0: ns(bigpot.ddom) = bigpot.dsizes;
wolffd@0: ns(bigpot.cdom) = bigpot.csizes;
wolffd@0: 
wolffd@0: % sum(ns(csumover))==0 is like isempty(csumover) but handles observed nodes.
wolffd@0: % Similarly, prod(ns(dsumover))==1 is like isempty(dsumover)
wolffd@0: 
wolffd@0: % Marginalize the cts parts.
wolffd@0: % If we are in canonical form, we stay that way, since moment form might not exist.
wolffd@0: % Besides, we would like to minimize the number of conversions.
wolffd@0: if sum(ns(csumover)) > 0
wolffd@0:   if bigpot.subtype == 'm'
wolffd@0:     for i=1:bigpot.dsize
wolffd@0:       bigpot.mom{i} = marginalize_pot(bigpot.mom{i}, ckeep);
wolffd@0:     end
wolffd@0:   else
wolffd@0:     for i=1:bigpot.dsize
wolffd@0:       bigpot.can{i} = marginalize_pot(bigpot.can{i}, ckeep);
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % If we are not marginalizing over any discrete nodes, we are done.
wolffd@0: if prod(ns(dsumover))==1
wolffd@0:   smallpot = cgpot(dkeep, ckeep, ns, bigpot.can, bigpot.mom, bigpot.subtype);
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: % To marginalize the discrete parts, we partition the cts parts into those that depend
wolffd@0: % on dkeep (i) and those that depend on on dsumover (j).
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: % If in canonical form, marginalize if possible, else convert to moment form.
wolffd@0: if 0 & bigpot.subtype == 'c'
wolffd@0:   p1 = zeros(I,J);
wolffd@0:   h1 = zeros(C,J,I);
wolffd@0:   K1 = 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:       can = struct(bigpot.can{k}); % violate object privacy
wolffd@0:       p1(i,j) = exp(can.g);
wolffd@0:       if C > 0 % so mu1 and Sigma1 are non-empty
wolffd@0: 	h1(:,j,i) = can.h;
wolffd@0: 	K1(:,:,j,i) = can.K;
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0:   
wolffd@0:   % If the cts parts do not depend on j, we can just marginalize the weighting coefficient g.
wolffd@0:   jdepends = 0;
wolffd@0:   for i=1:I
wolffd@0:     for j=2:J
wolffd@0:       if ~approxeq(h1(:,j,i), h1(:,1,i)) | ~approxeq(K1(:,:,j,i), K1(:,:,1,i))
wolffd@0: 	jdepends = 1;
wolffd@0: 	break
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0: 
wolffd@0:   if ~jdepends
wolffd@0:     %g2 = log(sum(p1, 2));
wolffd@0:     g2 = zeros(I,1);
wolffd@0:     for i=1:I
wolffd@0:       s = sum(p1(i,:));
wolffd@0:       if s > 0
wolffd@0: 	g2(i) = log(s);
wolffd@0:       end
wolffd@0:     end
wolffd@0:     h2 = h1;
wolffd@0:     K2 = K1;
wolffd@0:     can = cell(1,I);
wolffd@0:     j = 1; % arbitrary
wolffd@0:     for i=1:I
wolffd@0:       can{i} = cpot(ckeep, ns(ckeep), g2(i), h2(:,j,i), K2(:,:,j,i));
wolffd@0:     end
wolffd@0:     smallpot = cgpot(dkeep, ckeep, ns, can, [], 'c');  
wolffd@0:     return;
wolffd@0:   else
wolffd@0:     % Since the cts parts depend on j, we must convert to moment form
wolffd@0:     bigpot = cg_can_to_mom(bigpot);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % Marginalize in moment form
wolffd@0: bigpot = cg_can_to_mom(bigpot);
wolffd@0: 
wolffd@0: % Now partition the moment components.
wolffd@0: T1 = zeros(I,J);
wolffd@0: mu1 = zeros(C,J,I);
wolffd@0: Sigma1 = 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:     mom = struct(bigpot.mom{k}); % violate object privacy
wolffd@0:     T1(i,j) = exp(mom.logp);
wolffd@0:     if C > 0 % so mu1 and Sigma1 are non-empty
wolffd@0:       mu1(:,j,i) = mom.mu;
wolffd@0:       Sigma1(:,:,j,i) = mom.Sigma;
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % Collapse the mixture of Gaussians
wolffd@0: coef = mk_stochastic(T1); % coef must be convex combination
wolffd@0: T2 = sum(T1,2);
wolffd@0: T2 = T2 + (T2==0)*eps;
wolffd@0: %if C > 0, disp('collapsing onto '); disp(leep); end
wolffd@0: mu = [];
wolffd@0: Sigma = [];
wolffd@0: mom = cell(1,I);
wolffd@0: for i=1:I
wolffd@0:   if C > 0
wolffd@0:     [mu, Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:));
wolffd@0:   end
wolffd@0:   logp = log(T2(i));
wolffd@0:   mom{i} = mpot(ckeep, ns(ckeep), logp, mu, Sigma);
wolffd@0: end
wolffd@0: 
wolffd@0: smallpot = cgpot(dkeep, ckeep, ns, [], mom, 'm');
wolffd@0: