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