comparison toolboxes/FullBNT-1.0.7/bnt/potentials/@scgpot/marginalize_pot.m @ 0:e9a9cd732c1e tip

first hg version after svn

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