function [PDAGs, nodes] = mk_nbrs_of_pdag_del(cpdag,engine)
% MK_NBRS_OF_PDAG_ADD Make the inferior inclusion boundary of CPDAG.
% [PDAGs, nodes] = mk_nbrs_of_pdag_del(CPDAG)
%
% PDAGs{i} is the i'th neighbor of CPDAG0 generated by DELETE(X,Y,H)  with
% nodes{i,1:2}=[X Y]
% nodes{i,3}=H
%
% See D.M. Chickering 2002 : "Optimal Structure Identification with Greedy Search".
%
% philippe.leray@univ-nantes.fr
% 25 july 2003

compteur=0;
N=length(cpdag);
if nargin==1,
 bnet_tmp=mk_bnet(pdag_to_dag(cpdag),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;

verbose=0;

if verbose
    fprintf('---- Les cliques (maximales) du DAG\n');
    for i=1:nbcliques,
        disp(cliques{i});
    end
    fprintf('---- Variables \n');
end


% find in the PDAG all the X Y connected
[LX LY]=find(cpdag==1);
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);

    % this function recursively "walks" (dfs) in the graph representation of NA powerset
    if verbose
        X, Y, NY, AX, NAYX
        fprintf('---- NA Powerset\n\n');
    end

    liste=NAYX;
    if ~isempty(liste)
        premier=liste(1);
        dernier=liste(end);
    end

    current_set=[];
    fini=0;
    evite2 = 0 ;

    while ~fini
        if verbose
            % sert uniquement � l'affichage pour le prog de test ...
            if isempty(current_set)
                fprintf('\t H = []\n\n');
            else
                fprintf('\t H = '); disp(current_set);
            end
        end

        isclique=0;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 1
        % is [NAYX \ current_set] a clique ?
        if isempty(current_set)
            NAYXH=NAYX;
        elseif isempty(NAYX)
            NAYXH=current_set;
        else
            NAYXH = mysetdiff(NAYX,current_set);
        end

        if isempty(NAYXH)
            if verbose
                fprintf('\t\t NA(Y,X) \\ T = \t[]\n\n');
            end

        elseif verbose
            fprintf('\t\t NA(Y,X) \\ T ='); disp(NAYXH);
        end

        isclique = isempty(NAYXH) | ismemberclique(NAYXH,cliques) ;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 1
        if isclique
            % avoid testing INSERT(X,Y,H) and INSERT(Y,X,H) if X--Y
            test0=1 ;
            if (X>Y)
                test0=~cpdag(Y,X);
                if verbose&~test0 
                    fprintf('  (inutile de tester %d %d ',X,Y); fprintf('%d',current_set);
                    fprintf(')\n');
                end
            end
            if test0
                if verbose
                    fprintf(' ==DELETE(%d,%d,',X,Y); fprintf('%d',current_set);
                    fprintf(')\n');
                end
                compteur=compteur+1;
                nodes{compteur,1}=X;
                nodes{compteur,2}=Y;
                nodes{compteur,3}=current_set;
                ptmp=cpdag;
                ptmp(X,Y)=0;
                ptmp(Y,X)=0;
                ptmp(current_set,Y)=0;
                ptmp(Y,current_set)=1;
                ptmp(current_set(find(cpdag(X,current_set)==1)),X)=0;
                PDAGs{compteur}=ptmp;
            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