annotate toolboxes/FullBNT-1.0.7/bnt/potentials/@scgpot/marginalize_pot.m @ 0:cc4b1211e677 tip

initial commit to HG from Changeset: 646 (e263d8a21543) added further path and more save "camirversion.m"
author Daniel Wolff
date Fri, 19 Aug 2016 13:07:06 +0200
parents
children
rev   line source
Daniel@0 1 function smallpot = marginalize_pot(bigpot, keep)
Daniel@0 2 % MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
Daniel@0 3 % smallpot = marginalize_pot(bigpot, keep)
Daniel@0 4
Daniel@0 5 sumover = mysetdiff(bigpot.domain, keep);
Daniel@0 6 cdom = myunion(bigpot.cheaddom, bigpot.ctaildom);
Daniel@0 7 csumover = myintersect(sumover, bigpot.cheaddom);
Daniel@0 8 dsumover = myintersect(sumover, bigpot.ddom);
Daniel@0 9
Daniel@0 10 dkeep = myintersect(keep, bigpot.ddom);
Daniel@0 11 ckeep = myintersect(keep, bigpot.cheaddom);
Daniel@0 12 cheaddom = myintersect(keep, bigpot.cheaddom);
Daniel@0 13
Daniel@0 14 assert(isempty(myintersect(csumover,bigpot.ctaildom)));
Daniel@0 15 ns = zeros(1, max(bigpot.domain));
Daniel@0 16 ns(bigpot.ddom) = bigpot.dsizes;
Daniel@0 17 ns(bigpot.cheaddom) = bigpot.cheadsizes;
Daniel@0 18 ns(bigpot.ctaildom) = bigpot.ctailsizes;
Daniel@0 19
Daniel@0 20
Daniel@0 21 if sum(ns(csumover)) > 0
Daniel@0 22 for i=1:bigpot.dsize
Daniel@0 23 bigpot.scgpotc{i} = marginalize_pot(bigpot.scgpotc{i}, ckeep, csumover, ns);
Daniel@0 24 end
Daniel@0 25 end
Daniel@0 26
Daniel@0 27 if (isequal(csumover, cheaddom))
Daniel@0 28 bigpot.ctaildom = [];
Daniel@0 29 end
Daniel@0 30 % If we are not marginalizing over any discrete nodes, we are done.
Daniel@0 31 if prod(ns(dsumover))==1
Daniel@0 32 smallpot = scgpot(dkeep, cheaddom, bigpot.ctaildom, ns, bigpot.scgpotc);
Daniel@0 33 return;
Daniel@0 34 end
Daniel@0 35
Daniel@0 36 if (~isempty(bigpot.ctaildom))
Daniel@0 37 assert(0);
Daniel@0 38 return;
Daniel@0 39 end
Daniel@0 40
Daniel@0 41 I = prod(ns(dkeep));
Daniel@0 42 J = prod(ns(dsumover));
Daniel@0 43 C = sum(ns(ckeep));
Daniel@0 44 sum_map = find_equiv_posns(dsumover, bigpot.ddom);
Daniel@0 45 keep_map = find_equiv_posns(dkeep, bigpot.ddom);
Daniel@0 46 iv = zeros(1, length(bigpot.ddom)); % index vector
Daniel@0 47
Daniel@0 48 p1 = zeros(I,J);
Daniel@0 49 A1 = zeros(C,J,I);
Daniel@0 50 C1 = zeros(C,C,J,I);
Daniel@0 51 for i=1:I
Daniel@0 52 keep_iv = ind2subv(ns(dkeep), i);
Daniel@0 53 iv(keep_map) = keep_iv;
Daniel@0 54 for j=1:J
Daniel@0 55 sum_iv = ind2subv(ns(dsumover), j);
Daniel@0 56 iv(sum_map) = sum_iv;
Daniel@0 57 k = subv2ind(ns(bigpot.ddom), iv);
Daniel@0 58 pot = struct(bigpot.scgpotc{k}); % violate object privacy
Daniel@0 59 p1(i,j) = pot.p;
Daniel@0 60 if C > 0 % so mu1 and Sigma1 are non-empty
Daniel@0 61 A1(:,j,i) = pot.A;
Daniel@0 62 C1(:,:,j,i) = pot.C;
Daniel@0 63 end
Daniel@0 64 end
Daniel@0 65 end
Daniel@0 66
Daniel@0 67 % Collapse the mixture of Gaussians
Daniel@0 68 coef = mk_stochastic(p1); % coef must be convex combination
Daniel@0 69 %keyboard
Daniel@0 70 p2 = sum(p1,2);
Daniel@0 71 if (all(p2 == 0))
Daniel@0 72 p2 = p2 + (p2==0)*eps;
Daniel@0 73 end
Daniel@0 74 A = [];
Daniel@0 75 S = [];
Daniel@0 76
Daniel@0 77 pot = cell(1,I);
Daniel@0 78 ctailsize = sum(ns(bigpot.ctaildom));
Daniel@0 79 tB = zeros(C, ctailsize);
Daniel@0 80 for i=1:I
Daniel@0 81 if C > 0
Daniel@0 82 [A, S] = collapse_mog(A1(:,:,i), C1(:,:,:,i), coef(i,:));
Daniel@0 83 end
Daniel@0 84 p = p2(i);
Daniel@0 85 pot{i} = scgcpot(C, ctailsize, p, A, tB, S);
Daniel@0 86 end
Daniel@0 87
Daniel@0 88 smallpot = scgpot(dkeep, ckeep, bigpot.ctaildom, ns, pot);
Daniel@0 89
Daniel@0 90
Daniel@0 91
Daniel@0 92