mk_nbrs_of_pdag_add.m.svn-base

来自「bayesian network structrue learning mat」· SVN-BASE 代码 · 共 164 行

function [PDAGs, nodes] = mk_nbrs_of_pdag_add(cpdag)
% MK_NBRS_OF_PDAG_ADD Make the superior inclusion boundary of CPDAG.
% [PDAGs, nodes] = mk_nbrs_of_pdag_add(CPDAG)
%
% PDAGs{i} is the i'th neighbor of CPDAG0 generated by INSERT(X,Y,T)  with
% nodes{i,1:2}=[X Y]
% nodes{i,3}=T
%
% See D.M. Chickering 2002 : "Optimal Structure Identification with Greedy Search".
%
% [e-mail]
% 24 july 2003


compteur=0 ;
N=length(cpdag);
bnet_tmp=mk_bnet(pdag_to_dag(cpdag),2*ones(N,1));
engine_tmp=struct(jtree_inf_engine(bnet_tmp));
cliques=engine_tmp.cliques;
nbcliques=length(cliques);
clear bnet_tmp engine_tmp;

% find in the PDAG all the X Y not connected
[LX LY]=find((cpdag|cpdag')+eye(N)==0);
nlinks=length(LX);

for i=1:nlinks
    X=LX(i);
    Y=LY(i);
    % Neighbors of Y
    NY = myintersect(find(cpdag(:,Y)), find(cpdag(Y,:)));
    % Adjacents of X
    AX =  myunion(find(cpdag(:,X)), find(cpdag(X,:)));
    % Neighbors of Y adjacent to X    
    NAYX = myintersect(NY,AX);
    % Neighbors of Y NOT adjacent to X    
    NNAYX = mysetdiff(NY, NAYX);
    
    
    % this function recursively "walks" (dfs) in the graph representation of NNA powerset
    liste=NNAYX;    
    if ~isempty(liste)
        premier=liste(1);
        dernier=liste(end);
    end;     
    
    current_set=[];
    fini=0;
    evite2 = 0 ;
    
    while ~fini
        isclique=0;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 1  
        % is [NAYX current_set] a clique ?
        if isempty(current_set)
            NAYXT=NAYX;
        elseif isempty(NAYX)
            NAYXT=current_set;
        else
            NAYXT = union(NAYX,current_set);
        end
        
%        if isempty(NAYXT)
 %           isclique=1;
 %      end
        
        isclique = isempty(NAYXT) | ismemberclique(NAYXT,cliques) ;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 1        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 2
        % is there exist a partially directed path Y...X in PDAG \ NAYXT ?        
        if isclique  
            if evite2
                test2=1 ;
            else
                %%%% calcul test 2
                L2 = setdiff(1:N,NAYXT);
                test2=~partialconnected(cpdag(L2,L2),find(L2==Y),find(L2==X));
                
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 2        
            
            if test2
                evitet2=1;
                
                % avoid testing INSERT(X,Y,[]) and INSERT(Y,X,[]) if pa(Y)=pa(X)
                test0=1 ;
                if (X>Y)
                    if length(current_set)==0
                        PaX=setdiff(find(cpdag(X,:)),find(cpdag(:,X)));
                        PaY=setdiff(find(cpdag(Y,:)),find(cpdag(:,Y)));
                        test0=(length(setdiff(PaX,PaY))~=0);
                    end
                end
                if test0                
                    %fprintf(' INSERT(%d,%d,',X,Y); fprintf('%d',current_set);
                    %fprintf(')\n');
                    compteur=compteur+1;
                    nodes{compteur,1}=X;
                    nodes{compteur,2}=Y;
                    nodes{compteur,3}=current_set;
                    ptmp=cpdag;
                    ptmp(X,Y)=1;
                    ptmp(current_set,Y)=1;
                    ptmp(Y,current_set)=0;
                    PDAGs{compteur}=ptmp;                    
                end
            end
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Next Set to test ...        
        % what is the net set in the powerset ?
        if length(liste)==0
            fini=1 ; next_set=[];
        elseif length(current_set)==0
            next_set=[premier];     % first node after the root []
        else
            actuel=current_set(end);
            if actuel==dernier 
                if length(current_set)==1
                    fini=1;            % no more node ...
                else
                    ancien=current_set(end-1);  % new "branch"
                    next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)];
                end
            else                            % new node in the "branch"
                if ~isclique
                    if length(current_set)==1
                        fini=1;            % no more node ...
                    else
                        ancien=current_set(end-1);  % new "branch"
                        next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)];
                        evite2=0;
                    end
                else
                    next_set=[current_set liste(find(liste==actuel)+1)];
                end
            end
        end
        current_set=next_set;
    end
end


%%%%%%%
function resu = ismemberclique(v,cliques)

finiclique = 0 ; resu=0 ; 
cl=1;  ncl=length(cliques) ;

while (~finiclique) & (cl<=ncl);
    if ismember(v,cliques{cl})
        resu=1;
        finiclique=1 ;
    end
    cl=cl+1;
end

%%%%%%%%%
function resu = partialconnected(G,Y,X)

tmp=expm(G);

resu = (tmp(Y,X)~=0);