Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/potentials/@scgpot/complement_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 [margpot, comppot] = complement_pot(pot, keep)
wolffd@0 2 % COMPLEMENT_POT complement means decompose of a potential into its strong marginal and
wolffd@0 3 % its complement corresponds exactly to the decomposition of a probability distribution
wolffd@0 4 % into its marginal and conditional
wolffd@0 5 % [margpot, comppot] = complement_pot(pot, keep)
wolffd@0 6
wolffd@0 7 % keep can only include continuous head nodes and discrete nodes
wolffd@0 8 % margpot is the stable CG potential of keep nodes
wolffd@0 9 % comppot is the stable CG potential of others in corresponds exactly to
wolffd@0 10 % the discomposition of a probability distribution of its marginal and conditional
wolffd@0 11
wolffd@0 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 13 % Calculation of the marginal requires integration over %
wolffd@0 14 % all variables in csumover. Thus cheadkeep contains all %
wolffd@0 15 % continuous variables in the marginal potential %
wolffd@0 16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 17 %keyboard;
wolffd@0 18 csumover = mysetdiff(pot.cheaddom, keep);
wolffd@0 19 cheadkeep = mysetdiff(pot.cheaddom, csumover);
wolffd@0 20
wolffd@0 21 nodesizes = zeros(1, max(pot.domain));
wolffd@0 22 nodesizes(pot.ddom) = pot.dsizes;
wolffd@0 23 nodesizes(pot.cheaddom) = pot.cheadsizes;
wolffd@0 24 nodesizes(pot.ctaildom) = pot.ctailsizes;
wolffd@0 25
wolffd@0 26 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 27 % Description of the variables in the marginal domain %
wolffd@0 28 % For the calculation of a strong marginal first integration %
wolffd@0 29 % over all continuous variables in the head takes place. %
wolffd@0 30 % The calculation of the marginal over the head variables %
wolffd@0 31 % might result in a smaller or empty tail %
wolffd@0 32 % If there are no head variables, and therefore no tail %
wolffd@0 33 % variables, left marginalisation over discrete variables %
wolffd@0 34 % may take place %
wolffd@0 35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 36 margdom = mysetdiff(pot.domain,keep);
wolffd@0 37 % margddom = pot.ddom;
wolffd@0 38 margcheaddom = cheadkeep;
wolffd@0 39
wolffd@0 40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 41 % Marginalisation over discrete variables is only allowed when %
wolffd@0 42 % the tail is empty %
wolffd@0 43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 44 margddom = myintersect(pot.ddom,keep); % Discrete domain of marginal
wolffd@0 45 margctaildom = myintersect(pot.ctaildom,keep); % Tail domain
wolffd@0 46 assert(isempty(mysetdiff(pot.ddom,margddom)) | isempty(margctaildom))
wolffd@0 47
wolffd@0 48
wolffd@0 49 %margctaildom = pot.ctaildom;
wolffd@0 50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 51 % Even if marginalisation over continuous variables is only defined %
wolffd@0 52 % for head variables, the marginalisation over haed-variables might %
wolffd@0 53 % result in a smaller tail %
wolffd@0 54 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 55 margctaildom = myintersect(pot.ctaildom,keep);
wolffd@0 56
wolffd@0 57 margcheadsizes = nodesizes(margcheaddom);
wolffd@0 58 margcheadsize = sum(margcheadsizes);
wolffd@0 59 margctailsizes = nodesizes(margctaildom);
wolffd@0 60 margctailsize = sum(margctailsizes);
wolffd@0 61
wolffd@0 62 compdom = pot.domain;
wolffd@0 63 compddom = pot.ddom;
wolffd@0 64 compcheaddom = csumover;
wolffd@0 65 compctaildom = myunion(pot.ctaildom, cheadkeep);
wolffd@0 66 compcheadsizes = nodesizes(compcheaddom);
wolffd@0 67 compcheadsize = sum(compcheadsizes);
wolffd@0 68 compctailsizes = nodesizes(compctaildom);
wolffd@0 69 compctailsize = sum(compctailsizes);
wolffd@0 70
wolffd@0 71 dkeep = myintersect(pot.ddom, keep);
wolffd@0 72 %if dom is only contain discrete node
wolffd@0 73 if isempty(pot.cheaddom)
wolffd@0 74 dsumover = mysetdiff(pot.ddom, dkeep);
wolffd@0 75
wolffd@0 76 if isempty(dsumover)
wolffd@0 77 margpot = pot;
wolffd@0 78 comppot = scgpot([], [], [], []);
wolffd@0 79 return;
wolffd@0 80 end
wolffd@0 81
wolffd@0 82
wolffd@0 83 I = prod(nodesizes(dkeep));
wolffd@0 84 J = prod(nodesizes(dsumover));
wolffd@0 85 sum_map = find_equiv_posns(dsumover, pot.ddom);
wolffd@0 86 keep_map = find_equiv_posns(dkeep, pot.ddom);
wolffd@0 87 iv = zeros(1, length(pot.ddom)); % index vector
wolffd@0 88 p1 = zeros(I,J);
wolffd@0 89 for i=1:I
wolffd@0 90 keep_iv = ind2subv(nodesizes(dkeep), i);
wolffd@0 91 iv(keep_map) = keep_iv;
wolffd@0 92 for j=1:J
wolffd@0 93 sum_iv = ind2subv(nodesizes(dsumover), j);
wolffd@0 94 iv(sum_map) = sum_iv;
wolffd@0 95 k = subv2ind(nodesizes(pot.ddom), iv);
wolffd@0 96 potc = struct(pot.scgpotc{k}); % violate object privacy
wolffd@0 97 p1(i,j) = potc.p;
wolffd@0 98 end
wolffd@0 99 end
wolffd@0 100 p2 = sum(p1,2);
wolffd@0 101 p2 = p2 + (p2==0)*eps;
wolffd@0 102
wolffd@0 103 margscpot = cell(1, I);
wolffd@0 104 compscpot = cell(1, I*J);
wolffd@0 105 iv = zeros(1, length(pot.ddom)); % index vector
wolffd@0 106 for i=1:I
wolffd@0 107 margscpot{i} = scgcpot(0, 0, p2(i));
wolffd@0 108 keep_iv = ind2subv(nodesizes(dkeep), i);
wolffd@0 109 iv(keep_map) = keep_iv;
wolffd@0 110 for j=1:J
wolffd@0 111 sum_iv = ind2subv(nodesizes(dsumover), j);
wolffd@0 112 iv(sum_map) = sum_iv;
wolffd@0 113 k = subv2ind(nodesizes(pot.ddom), iv);
wolffd@0 114 q = p1(i,j)/p2(i);
wolffd@0 115 compscpot{k} = scgcpot(0, 0, q);
wolffd@0 116 end
wolffd@0 117 end
wolffd@0 118
wolffd@0 119 margpot = scgpot(dkeep, [], [], nodesizes, margscpot);
wolffd@0 120 comppot = scgpot(pot.ddom, [], [], nodesizes,compscpot);
wolffd@0 121 return;
wolffd@0 122 end
wolffd@0 123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 124 % head of the potential is not empty %
wolffd@0 125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 126 dsize = pot.dsize;
wolffd@0 127 compscpot = cell(1, dsize);
wolffd@0 128
wolffd@0 129 fmaskh = find_equiv_posns(margcheaddom, compctaildom);
wolffd@0 130 fmaskt = find_equiv_posns(margctaildom, compctaildom);
wolffd@0 131
wolffd@0 132 fh = block(fmaskh, compctailsizes);
wolffd@0 133 ft = block(fmaskt, compctailsizes);
wolffd@0 134
wolffd@0 135
wolffd@0 136 if ~isempty(margcheaddom)
wolffd@0 137 for i=1:dsize
wolffd@0 138 potc = struct(pot.scgpotc{i});
wolffd@0 139 q = 1;
wolffd@0 140 p = potc.p;
wolffd@0 141 [A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(potc.A, potc.B, potc.C, margcheaddom, compcheaddom, nodesizes);
wolffd@0 142
wolffd@0 143 if ~isempty(margcheaddom)
wolffd@0 144 margscpot{i} = scgcpot(margcheadsize, margctailsize, p, A1, B1, C11);
wolffd@0 145 else
wolffd@0 146 margscpot{i} = scgcpot(margcheadsize, margctailsize, p);
wolffd@0 147 end
wolffd@0 148
wolffd@0 149 if ~isempty(compcheaddom)
wolffd@0 150 if ~isempty(margcheaddom)
wolffd@0 151 E = A2 - C21*pinv(C11)*A1;
wolffd@0 152 tmp1 = C21*pinv(C11);
wolffd@0 153 tmp2 = B2 - C21*pinv(C11)*B1;
wolffd@0 154 F = zeros(compcheadsize, compctailsize);
wolffd@0 155 F(:, fh) = tmp1;
wolffd@0 156 F(:, ft) = tmp2;
wolffd@0 157 G = C22 - C21*pinv(C11)*C12;
wolffd@0 158 else
wolffd@0 159 E = A2;
wolffd@0 160 F = B2;
wolffd@0 161 G = C22;
wolffd@0 162 end
wolffd@0 163 compscpot{i} = scgcpot(compcheadsize, compctailsize, q, E, F, G);
wolffd@0 164 else
wolffd@0 165 compscpot{i} = scgcpot(compcheadsize, 0, q);
wolffd@0 166 end
wolffd@0 167 if isempty(margcheaddom)
wolffd@0 168 margpot = scgpot(margddom, [], [], nodesizes, margscpot);
wolffd@0 169 else
wolffd@0 170 margpot = scgpot(margddom, margcheaddom, margctaildom, nodesizes, margscpot);
wolffd@0 171 end
wolffd@0 172 end
wolffd@0 173 else
wolffd@0 174 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 175 % Marginalisation took place over all head variables. %
wolffd@0 176 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 177
wolffd@0 178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 179 % Calculate the strong marginal %
wolffd@0 180 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 181 margpot = marginalize_pot(pot,keep);
wolffd@0 182 mPot = struct(margpot);
wolffd@0 183 for i =1:dsize
wolffd@0 184 potc = struct(pot.scgpotc{i});
wolffd@0 185 % Get the probability of the original potential %
wolffd@0 186 q = potc.p;
wolffd@0 187
wolffd@0 188 % Get the configuration defined by the index i%
wolffd@0 189 config = ind2subv(pot.dsizes,i);
wolffd@0 190
wolffd@0 191 % Calculate the corresponding configuration in the marginal potential
wolffd@0 192 if isempty(margpot.dsizes)
wolffd@0 193 % keep == []
wolffd@0 194 indMargPot = 1;
wolffd@0 195 else
wolffd@0 196 equivPos = find_equiv_posns(dkeep,pot.ddom);
wolffd@0 197 indMargPot = subv2ind(margpot.dsizes,config(equivPos));
wolffd@0 198 end
wolffd@0 199 % Figure out the corresponding marginal potential
wolffd@0 200 mPotC = struct(mPot.scgpotc{indMargPot});
wolffd@0 201 p = mPotC.p;
wolffd@0 202 if p == 0
wolffd@0 203 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 204 % The following assignment is correct as p is only zero if q is also zero %
wolffd@0 205 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 206 compscpot{i} = scgcpot(compcheadsize,compctailsize,0,potc.A,potc.B,potc.C);
wolffd@0 207 else
wolffd@0 208 compscpot{i} = scgcpot(compcheadsize,compctailsize,q/p,potc.A,potc.B,potc.C);
wolffd@0 209 end
wolffd@0 210 end
wolffd@0 211 end
wolffd@0 212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 213 % Put all components in one potential %
wolffd@0 214 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 215 if isempty(compcheaddom)
wolffd@0 216 comppot = scgpot(compddom, [], [], nodesizes,compscpot);
wolffd@0 217 else
wolffd@0 218 comppot = scgpot(compddom, compcheaddom, compctaildom, nodesizes,compscpot);
wolffd@0 219 end
wolffd@0 220
wolffd@0 221