annotate toolboxes/FullBNT-1.0.7/bnt/potentials/@scgpot/direct_combine_pots.m @ 0:e9a9cd732c1e tip

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 pot = direct_combine_pots(pot1, pot2)
wolffd@0 2 % DIRECTED_COMBINE_POTS The combination operation corresponds to ordinary composition of conditional distributions.
wolffd@0 3 % In some sense is similar to that of forming disjoint union of set.
wolffd@0 4 % pot = direct_combine_pots(pot1, pot2)
wolffd@0 5
wolffd@0 6 % directed combine can be performed under the conditon that the head node set of pot1 is disjoint from the domain of
wolffd@0 7 % pot2 or vice versa. if the last conditon was satisfied we exchange the pot1 and pot2 firstly then perform the operation.
wolffd@0 8 % If neither of them was satified the directed combine is undifined.
wolffd@0 9
wolffd@0 10
wolffd@0 11 if isempty( myintersect(pot1.domain, pot2.cheaddom) )
wolffd@0 12 pot1 = pot1;
wolffd@0 13 pot2 = pot2;
wolffd@0 14 elseif isempty( myintersect(pot2.domain, pot1.cheaddom))
wolffd@0 15 temppot = pot1;
wolffd@0 16 pot1 = pot2;
wolffd@0 17 pot2 = temppot;
wolffd@0 18 else
wolffd@0 19 assert(0);
wolffd@0 20 return;
wolffd@0 21 end
wolffd@0 22
wolffd@0 23 domain = myunion(pot1.domain, pot2.domain);
wolffd@0 24 nodesizes = zeros(1,max(domain));
wolffd@0 25 nodesizes(pot2.ctaildom) = pot2.ctailsizes;
wolffd@0 26 nodesizes(pot2.cheaddom) = pot2.cheadsizes;
wolffd@0 27 nodesizes(pot2.ddom) = pot2.dsizes;
wolffd@0 28 nodesizes(pot1.ctaildom) = pot1.ctailsizes;
wolffd@0 29 nodesizes(pot1.cheaddom) = pot1.cheadsizes;
wolffd@0 30 nodesizes(pot1.ddom) = pot1.dsizes;
wolffd@0 31
wolffd@0 32 dom_u = mysetdiff(pot2.ctaildom, pot1.cheaddom);
wolffd@0 33 if ~isempty(dom_u) & ~mysubset(dom_u, pot1.ctaildom)
wolffd@0 34 pot1 = extension_pot(pot1, [], [], dom_u, nodesizes(dom_u));
wolffd@0 35 end
wolffd@0 36
wolffd@0 37 dom_u = myunion(pot1.cheaddom, pot1.ctaildom);
wolffd@0 38 if ~isempty(dom_u) & ~mysubset(dom_u, pot2.ctaildom)
wolffd@0 39 pot2 = extension_pot(pot2, [], [], dom_u, nodesizes(dom_u));
wolffd@0 40 end
wolffd@0 41
wolffd@0 42
wolffd@0 43 cheaddom = myunion(pot1.cheaddom, pot2.cheaddom);
wolffd@0 44 ctaildom = mysetdiff(myunion(pot1.ctaildom, pot2.ctaildom), cheaddom);
wolffd@0 45 cdom = myunion(cheaddom, ctaildom);
wolffd@0 46 ddom = mysetdiff(domain, cdom);
wolffd@0 47 dsizes = nodesizes(ddom);
wolffd@0 48 dsize = prod(nodesizes(ddom));
wolffd@0 49 cheadsizes = nodesizes(cheaddom);
wolffd@0 50 cheadsize = sum(nodesizes(cheaddom));
wolffd@0 51 ctailsizes = nodesizes(ctaildom);
wolffd@0 52 ctailsize = sum(nodesizes(ctaildom));
wolffd@0 53
wolffd@0 54 r1 = pot1.cheadsize;
wolffd@0 55 s1 = pot1.ctailsize;
wolffd@0 56 scpot = cell(1, dsize);
wolffd@0 57 mask1 = [];
wolffd@0 58 mask2 = [];
wolffd@0 59 if ~isempty(pot1.ddom)
wolffd@0 60 mask1 = find_equiv_posns(pot1.ddom, ddom);
wolffd@0 61 end
wolffd@0 62 if ~isempty(pot2.ddom)
wolffd@0 63 mask2 = find_equiv_posns(pot2.ddom, ddom);
wolffd@0 64 end
wolffd@0 65 cmask1 = [];
wolffd@0 66 cmask2 = [];
wolffd@0 67 if ~isempty(pot1.cheaddom)
wolffd@0 68 cmask1 = find_equiv_posns(pot1.cheaddom, cheaddom);
wolffd@0 69 end
wolffd@0 70 if ~isempty(pot2.cheaddom)
wolffd@0 71 cmask2 = find_equiv_posns(pot2.cheaddom, cheaddom);
wolffd@0 72 end
wolffd@0 73
wolffd@0 74 u1 = block(cmask1, cheadsizes);
wolffd@0 75 u2 = block(cmask2, cheadsizes);
wolffd@0 76
wolffd@0 77 fmaskh = find_equiv_posns(pot1.cheaddom, pot2.ctaildom);
wolffd@0 78 fmaskt = find_equiv_posns(pot1.ctaildom, pot2.ctaildom);
wolffd@0 79
wolffd@0 80 fh = block(fmaskh, pot2.ctailsizes);
wolffd@0 81 ft = block(fmaskt, pot2.ctailsizes);
wolffd@0 82
wolffd@0 83 for i=1:dsize
wolffd@0 84 sub = ind2subv(dsizes, i);
wolffd@0 85 sub1 = sub(mask1);
wolffd@0 86 sub2 = sub(mask2);
wolffd@0 87 ind1 = subv2ind(pot1.dsizes, sub1);
wolffd@0 88 ind2 = subv2ind(pot2.dsizes, sub2);
wolffd@0 89
wolffd@0 90 if isempty(ind1)
wolffd@0 91 ind1 = 1;
wolffd@0 92 end
wolffd@0 93 if isempty(ind2)
wolffd@0 94 ind2 = 1;
wolffd@0 95 end
wolffd@0 96 potc1 = struct(pot1.scgpotc{ind1});
wolffd@0 97 potc2 = struct(pot2.scgpotc{ind2});
wolffd@0 98 p = potc1.p;
wolffd@0 99 q = potc2.p;
wolffd@0 100 ro = p*q;
wolffd@0 101
wolffd@0 102 A = potc1.A;
wolffd@0 103 B = potc1.B;
wolffd@0 104 C = potc1.C;
wolffd@0 105
wolffd@0 106 E = potc2.A;
wolffd@0 107 F = potc2.B;
wolffd@0 108 G = potc2.C;
wolffd@0 109
wolffd@0 110 F1 = F(:, fh);
wolffd@0 111 F2 = F(:, ft);
wolffd@0 112
wolffd@0 113 if ~isempty(F1)
wolffd@0 114 K1 = F1*A;
wolffd@0 115 K2 = F1*B;
wolffd@0 116 FCF = F1*C*F1';
wolffd@0 117 FC = F1*C;
wolffd@0 118 CFT = C*F1';
wolffd@0 119 else
wolffd@0 120 K1 = zeros(size(E));
wolffd@0 121 K2 = zeros(size(F2));
wolffd@0 122 FCF = zeros(size(G));
wolffd@0 123 FC = zeros(size(C, 1), size(G, 2));
wolffd@0 124 CFT = zeros(size(G, 2), size(C, 1));
wolffd@0 125 end
wolffd@0 126
wolffd@0 127
wolffd@0 128 U = zeros(cheadsize,1);
wolffd@0 129 W = zeros(cheadsize,cheadsize);
wolffd@0 130 V = zeros(cheadsize,ctailsize);
wolffd@0 131
wolffd@0 132 if cheadsize > 0
wolffd@0 133 U(u1) = A;
wolffd@0 134 U(u2) = E + K1;
wolffd@0 135 W(u1, u1) = C;
wolffd@0 136 W(u2, u2) = G + FCF;
wolffd@0 137 W(u1, u2) = CFT;
wolffd@0 138 W(u2, u1) = FC;
wolffd@0 139 else
wolffd@0 140 U = zeros(cheadsize,1);
wolffd@0 141 W = zeros(cheadsize,cheadsize);
wolffd@0 142 end
wolffd@0 143 if cheadsize > 0 | ctailsize > 0
wolffd@0 144 if ~isempty(u1)
wolffd@0 145 V(u1, :) = B;
wolffd@0 146 else
wolffd@0 147 V(u1, :) = zeros(potc1.cheadsize, ctailsize);
wolffd@0 148 end
wolffd@0 149 if ~isempty(u2)
wolffd@0 150 V(u2, :) = F2 + K2;
wolffd@0 151 else
wolffd@0 152 V(u2, :) = zeros(potc2.cheadsize, ctailsize);
wolffd@0 153 end
wolffd@0 154 else
wolffd@0 155 V = zeros(cheadsize,ctailsize);
wolffd@0 156 end
wolffd@0 157
wolffd@0 158 scpot{i} = scgcpot(cheadsize, ctailsize, ro, U, V, W);
wolffd@0 159 end
wolffd@0 160
wolffd@0 161 pot = scgpot(ddom, cheaddom, ctaildom, nodesizes, scpot);