wolffd@0: function pot = direct_combine_pots(pot1, pot2)
wolffd@0: % DIRECTED_COMBINE_POTS The combination operation corresponds to ordinary composition of conditional distributions. 
wolffd@0: % In some sense is similar to that of forming disjoint union of set.
wolffd@0: % pot = direct_combine_pots(pot1, pot2)
wolffd@0: 
wolffd@0: % directed combine can be performed under the conditon that the head node set of pot1 is disjoint from the domain of 
wolffd@0: % pot2 or vice versa. if the last conditon was satisfied we exchange the pot1 and pot2 firstly then perform the operation.
wolffd@0: % If neither of them was satified the directed combine is undifined.
wolffd@0: 
wolffd@0: 
wolffd@0: if isempty( myintersect(pot1.domain, pot2.cheaddom) )
wolffd@0:     pot1 = pot1;
wolffd@0:     pot2 = pot2;
wolffd@0: elseif  isempty( myintersect(pot2.domain, pot1.cheaddom))
wolffd@0:     temppot = pot1;
wolffd@0:     pot1 = pot2;
wolffd@0:     pot2 = temppot;
wolffd@0: else
wolffd@0:     assert(0);
wolffd@0:     return;
wolffd@0: end
wolffd@0: 
wolffd@0: domain = myunion(pot1.domain, pot2.domain);
wolffd@0: nodesizes = zeros(1,max(domain));
wolffd@0: nodesizes(pot2.ctaildom) = pot2.ctailsizes;
wolffd@0: nodesizes(pot2.cheaddom) = pot2.cheadsizes;
wolffd@0: nodesizes(pot2.ddom) = pot2.dsizes;
wolffd@0: nodesizes(pot1.ctaildom) = pot1.ctailsizes;
wolffd@0: nodesizes(pot1.cheaddom) = pot1.cheadsizes;
wolffd@0: nodesizes(pot1.ddom) = pot1.dsizes;
wolffd@0: 
wolffd@0: dom_u = mysetdiff(pot2.ctaildom, pot1.cheaddom);
wolffd@0: if ~isempty(dom_u) & ~mysubset(dom_u, pot1.ctaildom)
wolffd@0:     pot1 = extension_pot(pot1, [], [], dom_u, nodesizes(dom_u));
wolffd@0: end
wolffd@0: 
wolffd@0: dom_u = myunion(pot1.cheaddom, pot1.ctaildom);
wolffd@0: if ~isempty(dom_u) & ~mysubset(dom_u, pot2.ctaildom)
wolffd@0:     pot2 = extension_pot(pot2, [], [], dom_u, nodesizes(dom_u));
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: cheaddom = myunion(pot1.cheaddom, pot2.cheaddom);
wolffd@0: ctaildom = mysetdiff(myunion(pot1.ctaildom, pot2.ctaildom), cheaddom);
wolffd@0: cdom = myunion(cheaddom, ctaildom);
wolffd@0: ddom = mysetdiff(domain, cdom);
wolffd@0: dsizes = nodesizes(ddom);
wolffd@0: dsize = prod(nodesizes(ddom));
wolffd@0: cheadsizes = nodesizes(cheaddom);
wolffd@0: cheadsize = sum(nodesizes(cheaddom));
wolffd@0: ctailsizes = nodesizes(ctaildom);
wolffd@0: ctailsize = sum(nodesizes(ctaildom));
wolffd@0: 
wolffd@0: r1 = pot1.cheadsize;
wolffd@0: s1 = pot1.ctailsize;
wolffd@0: scpot = cell(1, dsize);
wolffd@0: mask1 = [];
wolffd@0: mask2 = [];
wolffd@0: if ~isempty(pot1.ddom)
wolffd@0:     mask1 = find_equiv_posns(pot1.ddom, ddom);
wolffd@0: end
wolffd@0: if ~isempty(pot2.ddom)
wolffd@0:     mask2 = find_equiv_posns(pot2.ddom, ddom);
wolffd@0: end
wolffd@0: cmask1 = [];
wolffd@0: cmask2 = [];
wolffd@0: if ~isempty(pot1.cheaddom)
wolffd@0:     cmask1 = find_equiv_posns(pot1.cheaddom, cheaddom);
wolffd@0: end
wolffd@0: if ~isempty(pot2.cheaddom)
wolffd@0:     cmask2 = find_equiv_posns(pot2.cheaddom, cheaddom);
wolffd@0: end
wolffd@0: 
wolffd@0: u1 = block(cmask1, cheadsizes);
wolffd@0: u2 = block(cmask2, cheadsizes);
wolffd@0: 
wolffd@0: fmaskh = find_equiv_posns(pot1.cheaddom, pot2.ctaildom);
wolffd@0: fmaskt = find_equiv_posns(pot1.ctaildom, pot2.ctaildom);
wolffd@0: 
wolffd@0: fh = block(fmaskh, pot2.ctailsizes);
wolffd@0: ft = block(fmaskt, pot2.ctailsizes);
wolffd@0: 
wolffd@0: for i=1:dsize
wolffd@0:     sub = ind2subv(dsizes, i);
wolffd@0:     sub1 = sub(mask1);
wolffd@0:     sub2 = sub(mask2);
wolffd@0:     ind1 = subv2ind(pot1.dsizes, sub1);
wolffd@0:     ind2 = subv2ind(pot2.dsizes, sub2);
wolffd@0:     
wolffd@0:     if isempty(ind1)
wolffd@0:         ind1 = 1;
wolffd@0:     end
wolffd@0:     if isempty(ind2)
wolffd@0:         ind2 = 1;
wolffd@0:     end
wolffd@0:     potc1 = struct(pot1.scgpotc{ind1});
wolffd@0:     potc2 = struct(pot2.scgpotc{ind2});
wolffd@0:     p = potc1.p;
wolffd@0:     q = potc2.p;
wolffd@0:     ro = p*q;
wolffd@0:     
wolffd@0:     A = potc1.A;
wolffd@0:     B = potc1.B;
wolffd@0:     C = potc1.C;
wolffd@0:    
wolffd@0:     E = potc2.A;
wolffd@0:     F = potc2.B;
wolffd@0:     G = potc2.C;
wolffd@0:     
wolffd@0:     F1 = F(:, fh);
wolffd@0:     F2 = F(:, ft);
wolffd@0:     
wolffd@0:     if ~isempty(F1)
wolffd@0:         K1 = F1*A;
wolffd@0:         K2 = F1*B;
wolffd@0:         FCF = F1*C*F1';
wolffd@0:         FC = F1*C;
wolffd@0:         CFT = C*F1';
wolffd@0:     else
wolffd@0:         K1 = zeros(size(E));
wolffd@0:         K2 = zeros(size(F2));
wolffd@0:         FCF = zeros(size(G));
wolffd@0:         FC = zeros(size(C, 1), size(G, 2));
wolffd@0:         CFT = zeros(size(G, 2), size(C, 1));
wolffd@0:     end
wolffd@0:     
wolffd@0:     
wolffd@0:     U = zeros(cheadsize,1); 
wolffd@0:     W = zeros(cheadsize,cheadsize);
wolffd@0:     V = zeros(cheadsize,ctailsize); 
wolffd@0:     
wolffd@0:     if cheadsize > 0
wolffd@0:         U(u1) = A;
wolffd@0:         U(u2) = E + K1;
wolffd@0:         W(u1, u1) = C;
wolffd@0:         W(u2, u2) = G + FCF;
wolffd@0:         W(u1, u2) = CFT;
wolffd@0:         W(u2, u1) = FC;
wolffd@0:     else
wolffd@0:         U = zeros(cheadsize,1); 
wolffd@0:         W = zeros(cheadsize,cheadsize); 
wolffd@0:     end
wolffd@0:     if cheadsize > 0 | ctailsize > 0
wolffd@0:         if ~isempty(u1)
wolffd@0:             V(u1, :) = B;
wolffd@0:         else
wolffd@0:             V(u1, :) = zeros(potc1.cheadsize, ctailsize);
wolffd@0:         end
wolffd@0:         if ~isempty(u2)
wolffd@0:             V(u2, :) = F2 + K2;
wolffd@0:         else
wolffd@0:             V(u2, :) = zeros(potc2.cheadsize, ctailsize);
wolffd@0:         end
wolffd@0:     else
wolffd@0:         V = zeros(cheadsize,ctailsize); 
wolffd@0:     end
wolffd@0: 
wolffd@0:     scpot{i} = scgcpot(cheadsize, ctailsize, ro, U, V, W);
wolffd@0: end
wolffd@0: 
wolffd@0: pot = scgpot(ddom, cheaddom, ctaildom, nodesizes, scpot);