complement_pot.m

来自「麻省理工学院的人工智能工具箱,很珍贵,希望对大家有用!」· M 代码 · 共 144 行

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

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;

margdom = mysetdiff(pot.domain, csumover);
margddom = pot.ddom;
margcheaddom = cheadkeep;
margctaildom = pot.ctaildom;
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);

%if dom is only contain discrete node
if isempty(pot.cheaddom)  
    dkeep = myintersect(pot.ddom, keep);
    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


dsize = pot.dsize;
margscpot = cell(1, 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);

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
end

if isempty(margcheaddom)
    margpot = scgpot(margddom, [], [], nodesizes, margscpot);
else
    margpot = scgpot(margddom, margcheaddom, margctaildom, nodesizes, margscpot);
end

if isempty(compcheaddom)
    comppot = scgpot(compddom, [], [], nodesizes,compscpot);
else
    comppot = scgpot(compddom, compcheaddom, compctaildom, nodesizes,compscpot);
end