function [PDAGs, nodes] = mk_nbrs_of_pdag_add(cpdag,engine)
% 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".
%
% philippe.leray@univ-nantes.fr, francois.olivier.c.h@gmail.com
% 24 july 2003


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

% find in the PDAG all the X Y not connected
[LX LY]=find((cpdag|cpdag')+eye(N)==0);
nlinks=length(LX); % here is a bug when nlinks is zeros, for i=1:nlinks fail 31-37 added by hanbin
if nlinks==0
   [i,j]=find(cpdag);k=unidrnd(size(i,1),1);
   cpdag(i(k),j(k))=0;
   [LX LY]=find((cpdag|cpdag')+eye(N)==0);
    nlinks=length(LX);
end

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);