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