1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
function smallpot = marginalize_pot(bigpot, keep, maximize, useC)
% MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
% smallpot = marginalize_pot(bigpot, keep, maximize, useC)
%
% If maximize = 1, we raise an error.
% useC is ignored.
if nargin < 3, maximize = 0; end
assert(~maximize);
sumover = mysetdiff(bigpot.domain, keep);
csumover = myintersect(sumover, bigpot.cdom);
dsumover = myintersect(sumover, bigpot.ddom);
dkeep = myintersect(keep, bigpot.ddom);
ckeep = myintersect(keep, bigpot.cdom);
%ns = sparse(1, max(bigpot.domain)); % must be full, so I is an integer
ns = zeros(1, max(bigpot.domain));
ns(bigpot.ddom) = bigpot.dsizes;
ns(bigpot.cdom) = bigpot.csizes;
% sum(ns(csumover))==0 is like isempty(csumover) but handles observed nodes.
% Similarly, prod(ns(dsumover))==1 is like isempty(dsumover)
% Marginalize the cts parts.
% If we are in canonical form, we stay that way, since moment form might not exist.
% Besides, we would like to minimize the number of conversions.
if sum(ns(csumover)) > 0
if bigpot.subtype == 'm'
for i=1:bigpot.dsize
bigpot.mom{i} = marginalize_pot(bigpot.mom{i}, ckeep);
end
else
for i=1:bigpot.dsize
bigpot.can{i} = marginalize_pot(bigpot.can{i}, ckeep);
end
end
end
% If we are not marginalizing over any discrete nodes, we are done.
if prod(ns(dsumover))==1
smallpot = cgpot(dkeep, ckeep, ns, bigpot.can, bigpot.mom, bigpot.subtype);
return;
end
% To marginalize the discrete parts, we partition the cts parts into those that depend
% on dkeep (i) and those that depend on on dsumover (j).
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
% If in canonical form, marginalize if possible, else convert to moment form.
if 0 && bigpot.subtype == 'c'
p1 = zeros(I,J);
h1 = zeros(C,J,I);
K1 = 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);
can = struct(bigpot.can{k}); % violate object privacy
p1(i,j) = exp(can.g);
if C > 0 % so mu1 and Sigma1 are non-empty
h1(:,j,i) = can.h;
K1(:,:,j,i) = can.K;
end
end
end
% If the cts parts do not depend on j, we can just marginalize the weighting coefficient g.
jdepends = 0;
for i=1:I
for j=2:J
if ~approxeq(h1(:,j,i), h1(:,1,i)) || ~approxeq(K1(:,:,j,i), K1(:,:,1,i))
jdepends = 1;
break
end
end
end
if ~jdepends
%g2 = log(sum(p1, 2));
g2 = zeros(I,1);
for i=1:I
s = sum(p1(i,:));
if s > 0
g2(i) = log(s);
end
end
h2 = h1;
K2 = K1;
can = cell(1,I);
j = 1; % arbitrary
for i=1:I
can{i} = cpot(ckeep, ns(ckeep), g2(i), h2(:,j,i), K2(:,:,j,i));
end
smallpot = cgpot(dkeep, ckeep, ns, can, [], 'c');
return;
else
% Since the cts parts depend on j, we must convert to moment form
bigpot = cg_can_to_mom(bigpot);
end
end
% Marginalize in moment form
bigpot = cg_can_to_mom(bigpot);
% Now partition the moment components.
T1 = zeros(I,J);
mu1 = zeros(C,J,I);
Sigma1 = 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);
mom = struct(bigpot.mom{k}); % violate object privacy
T1(i,j) = exp(mom.logp);
if C > 0 % so mu1 and Sigma1 are non-empty
mu1(:,j,i) = mom.mu;
Sigma1(:,:,j,i) = mom.Sigma;
end
end
end
% Collapse the mixture of Gaussians
coef = mk_stochastic(T1); % coef must be convex combination
T2 = sum(T1,2);
T2 = T2 + (T2==0)*eps;
%if C > 0, disp('collapsing onto '); disp(leep); end
mu = [];
Sigma = [];
mom = cell(1,I);
for i=1:I
if C > 0
[mu, Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:));
end
logp = log(T2(i));
mom{i} = mpot(ckeep, ns(ckeep), logp, mu, Sigma);
end
smallpot = cgpot(dkeep, ckeep, ns, [], mom, 'm');
|
