wolffd@0: function [margpot, comppot] = complement_pot(pot, keep)
wolffd@0: % COMPLEMENT_POT complement means decompose of a potential into its strong marginal and 
wolffd@0: % its complement corresponds exactly to the decomposition of a probability distribution 
wolffd@0: % into its marginal and conditional
wolffd@0: % [margpot, comppot] = complement_pot(pot, keep)
wolffd@0: 
wolffd@0: % keep can only include continuous head nodes and discrete nodes
wolffd@0: % margpot is the stable CG potential of keep nodes
wolffd@0: % comppot is the stable CG potential of others in corresponds exactly to 
wolffd@0: % the discomposition of a probability distribution of its marginal and conditional
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: % Calculation of the marginal requires integration over      %
wolffd@0: % all variables in csumover. Thus cheadkeep contains all     %
wolffd@0: % continuous variables in the marginal potential             %
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %keyboard;
wolffd@0: csumover = mysetdiff(pot.cheaddom, keep);
wolffd@0: cheadkeep = mysetdiff(pot.cheaddom, csumover);
wolffd@0: 
wolffd@0: nodesizes = zeros(1, max(pot.domain));
wolffd@0: nodesizes(pot.ddom) = pot.dsizes;
wolffd@0: nodesizes(pot.cheaddom) = pot.cheadsizes;
wolffd@0: nodesizes(pot.ctaildom) = pot.ctailsizes;
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: % Description of the variables in the marginal domain        %
wolffd@0: % For the calculation of a strong marginal first integration %
wolffd@0: % over all continuous variables in the head takes place.     %
wolffd@0: % The calculation of the marginal over the head variables    %
wolffd@0: % might result in a smaller or empty tail                    %
wolffd@0: % If there are no head variables, and therefore no tail      %
wolffd@0: % variables, left marginalisation over discrete variables    %
wolffd@0: % may take place                                             %
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
wolffd@0: margdom      = mysetdiff(pot.domain,keep);
wolffd@0: % margddom   = pot.ddom;
wolffd@0: margcheaddom = cheadkeep;
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: % Marginalisation over discrete variables is only allowed when %
wolffd@0: % the tail is empty                                            %
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: margddom = myintersect(pot.ddom,keep);               % Discrete domain of marginal
wolffd@0: margctaildom = myintersect(pot.ctaildom,keep);       % Tail domain
wolffd@0: assert(isempty(mysetdiff(pot.ddom,margddom)) | isempty(margctaildom))  
wolffd@0: 
wolffd@0: 
wolffd@0: %margctaildom = pot.ctaildom;
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: % Even if marginalisation over continuous variables is only defined %
wolffd@0: % for head variables, the marginalisation over haed-variables might %
wolffd@0: % result in a smaller tail                                          %
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: margctaildom = myintersect(pot.ctaildom,keep);
wolffd@0: 
wolffd@0: margcheadsizes = nodesizes(margcheaddom);
wolffd@0: margcheadsize = sum(margcheadsizes);
wolffd@0: margctailsizes = nodesizes(margctaildom);
wolffd@0: margctailsize = sum(margctailsizes);
wolffd@0: 
wolffd@0: compdom = pot.domain;
wolffd@0: compddom = pot.ddom;
wolffd@0: compcheaddom = csumover;
wolffd@0: compctaildom = myunion(pot.ctaildom, cheadkeep);
wolffd@0: compcheadsizes = nodesizes(compcheaddom);
wolffd@0: compcheadsize = sum(compcheadsizes);
wolffd@0: compctailsizes = nodesizes(compctaildom);
wolffd@0: compctailsize = sum(compctailsizes);
wolffd@0: 
wolffd@0: dkeep = myintersect(pot.ddom, keep);
wolffd@0: %if dom is only contain discrete node
wolffd@0: if isempty(pot.cheaddom)
wolffd@0:     dsumover = mysetdiff(pot.ddom, dkeep);
wolffd@0:     
wolffd@0:     if isempty(dsumover)
wolffd@0:         margpot = pot;
wolffd@0:         comppot = scgpot([], [], [], []);
wolffd@0:         return;
wolffd@0:     end
wolffd@0:         
wolffd@0:     
wolffd@0:     I = prod(nodesizes(dkeep));
wolffd@0:     J = prod(nodesizes(dsumover));
wolffd@0:     sum_map = find_equiv_posns(dsumover, pot.ddom);
wolffd@0:     keep_map = find_equiv_posns(dkeep, pot.ddom);
wolffd@0:     iv = zeros(1, length(pot.ddom)); % index vector
wolffd@0:     p1 = zeros(I,J);
wolffd@0:     for i=1:I
wolffd@0:         keep_iv = ind2subv(nodesizes(dkeep), i);
wolffd@0:         iv(keep_map) = keep_iv;
wolffd@0:         for j=1:J
wolffd@0:             sum_iv = ind2subv(nodesizes(dsumover), j);
wolffd@0:             iv(sum_map) = sum_iv;
wolffd@0:             k = subv2ind(nodesizes(pot.ddom), iv);
wolffd@0:             potc = struct(pot.scgpotc{k}); % violate object privacy
wolffd@0:             p1(i,j) = potc.p;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     p2 = sum(p1,2);
wolffd@0:     p2 = p2 + (p2==0)*eps;
wolffd@0:     
wolffd@0:     margscpot = cell(1, I);
wolffd@0:     compscpot = cell(1, I*J);
wolffd@0:     iv = zeros(1, length(pot.ddom)); % index vector
wolffd@0:     for i=1:I
wolffd@0:         margscpot{i} = scgcpot(0, 0, p2(i));
wolffd@0:         keep_iv = ind2subv(nodesizes(dkeep), i);
wolffd@0:         iv(keep_map) = keep_iv;
wolffd@0:         for j=1:J
wolffd@0:             sum_iv = ind2subv(nodesizes(dsumover), j);
wolffd@0:             iv(sum_map) = sum_iv;
wolffd@0:             k = subv2ind(nodesizes(pot.ddom), iv);
wolffd@0:             q = p1(i,j)/p2(i);
wolffd@0:             compscpot{k} = scgcpot(0, 0, q);
wolffd@0:         end
wolffd@0:     end
wolffd@0:     
wolffd@0:     margpot = scgpot(dkeep, [], [], nodesizes, margscpot);
wolffd@0:     comppot = scgpot(pot.ddom, [], [], nodesizes,compscpot);
wolffd@0:     return;
wolffd@0: end
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: % head of the potential is not empty %
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: dsize = pot.dsize;
wolffd@0: compscpot = cell(1, dsize);
wolffd@0: 
wolffd@0: fmaskh = find_equiv_posns(margcheaddom, compctaildom);
wolffd@0: fmaskt = find_equiv_posns(margctaildom, compctaildom);
wolffd@0: 
wolffd@0: fh = block(fmaskh, compctailsizes);
wolffd@0: ft = block(fmaskt, compctailsizes);
wolffd@0: 
wolffd@0: 
wolffd@0: if ~isempty(margcheaddom)
wolffd@0:     for i=1:dsize
wolffd@0:         potc = struct(pot.scgpotc{i});
wolffd@0:         q = 1;
wolffd@0:         p = potc.p;
wolffd@0:         [A1, A2, B1, B2, C11, C12, C21, C22] = partition_matrix_vec_3(potc.A, potc.B, potc.C, margcheaddom, compcheaddom, nodesizes);
wolffd@0: 
wolffd@0:         if ~isempty(margcheaddom)
wolffd@0:             margscpot{i} = scgcpot(margcheadsize, margctailsize, p, A1, B1, C11);
wolffd@0:         else
wolffd@0:             margscpot{i} = scgcpot(margcheadsize, margctailsize, p);
wolffd@0:         end 
wolffd@0:     
wolffd@0:         if ~isempty(compcheaddom)
wolffd@0:             if ~isempty(margcheaddom)
wolffd@0:                 E = A2 - C21*pinv(C11)*A1;
wolffd@0:                 tmp1 = C21*pinv(C11);
wolffd@0:                 tmp2 = B2 - C21*pinv(C11)*B1;
wolffd@0:                 F = zeros(compcheadsize, compctailsize);
wolffd@0:                 F(:, fh) = tmp1;
wolffd@0:                 F(:, ft) = tmp2;
wolffd@0:                 G = C22 - C21*pinv(C11)*C12;
wolffd@0:             else
wolffd@0:                 E = A2;
wolffd@0:                 F = B2;
wolffd@0:                 G = C22;
wolffd@0:             end
wolffd@0:             compscpot{i} = scgcpot(compcheadsize, compctailsize, q, E, F, G);
wolffd@0:         else
wolffd@0:             compscpot{i} = scgcpot(compcheadsize, 0, q);
wolffd@0:         end
wolffd@0:         if isempty(margcheaddom)
wolffd@0:             margpot = scgpot(margddom, [], [], nodesizes, margscpot);
wolffd@0:         else
wolffd@0:             margpot = scgpot(margddom, margcheaddom, margctaildom, nodesizes, margscpot);
wolffd@0:         end
wolffd@0:     end
wolffd@0: else
wolffd@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0:     % Marginalisation took place over all head variables.                               %
wolffd@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: 
wolffd@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0:     % Calculate the strong marginal %
wolffd@0:     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0:     margpot = marginalize_pot(pot,keep);
wolffd@0:     mPot    = struct(margpot); 
wolffd@0:     for i =1:dsize
wolffd@0:         potc = struct(pot.scgpotc{i});  
wolffd@0:         % Get the probability of the original potential % 
wolffd@0: 	q = potc.p;
wolffd@0:          
wolffd@0:         % Get the configuration defined by the index i%
wolffd@0:         config = ind2subv(pot.dsizes,i);
wolffd@0:         
wolffd@0:         % Calculate the corresponding configuration in the marginal potential
wolffd@0:         if isempty(margpot.dsizes)
wolffd@0:             % keep == []
wolffd@0: 	    indMargPot = 1;
wolffd@0:         else
wolffd@0:             equivPos   = find_equiv_posns(dkeep,pot.ddom);
wolffd@0:             indMargPot = subv2ind(margpot.dsizes,config(equivPos));
wolffd@0:         end
wolffd@0:         % Figure out the corresponding marginal potential
wolffd@0:         mPotC = struct(mPot.scgpotc{indMargPot});
wolffd@0:         p = mPotC.p;
wolffd@0:         if p == 0
wolffd@0:             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0:             % The following assignment is correct as p is only zero if q is also zero %
wolffd@0:             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0:             compscpot{i} = scgcpot(compcheadsize,compctailsize,0,potc.A,potc.B,potc.C);
wolffd@0:         else
wolffd@0:             compscpot{i} = scgcpot(compcheadsize,compctailsize,q/p,potc.A,potc.B,potc.C);
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: % Put all components in one potential %
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: if isempty(compcheaddom)
wolffd@0:     comppot = scgpot(compddom, [], [], nodesizes,compscpot);
wolffd@0: else
wolffd@0:     comppot = scgpot(compddom, compcheaddom, compctaildom, nodesizes,compscpot);
wolffd@0: end
wolffd@0: 
wolffd@0: