Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/potentials/@cgpot/marginalize_pot.m @ 0:e9a9cd732c1e tip

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