wang_learn_struct_undirect2pdag(isdag,cond_indep, n, k, G,varargin)

function [pdag, G] = wang_learn_struct_undirect2pdag(isdag,cond_indep, n, k, G,varargin)
% LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm
% P = learn_struct_pdag_pc(cond_indep, n, k, ...)
%
% n is the number of nodes.
% k is an optional upper bound on the fan-in (default: n)
% cond_indep is a boolean function that will be called as follows:
%   feval(cond_indep, x, y, S, ...)
% where x and y are nodes, and S is a set of nodes (positive integers),
% and ... are any optional parameters passed to this function.
%
% The output P is an adjacency matrix, in which
% P(i,j) = -1 if there is an i->j edge.
% P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
%
% The PC algorithm does structure learning assuming all variables are observed.
% See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
% This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
% but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).
%
%this program seems match with Tetrad4.2 program but some edges were missing or the direction of edges were
% not the same as Tetrad4.2 I need test more and find out way?

%This one have to be checked again if combined with Indepent undirect Graph!!
%Oct 06.2003 junbai wanng
%
%Step A.B  

sep = cell(n,n);
ord = 0;
done = 0;
G=setdiag(G,0);
nsample_size=varargin{2};
while ~done
  done = 1;
  [X,Y] = find(G); 
  for i=1:length(X)
    x = X(i); y = Y(i);
    nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
    %add junbai oct 12
    if ( (length(nbrs) >= ord) & (G(x,y) ~= 0 & nsample_size-ord-3>=0) )
      done = 0;
      %SS = subsets(nbrs, ord, ord); % all subsets of size ord
      %bug 1 in subsets, now fixed by jbw 04.09.2003
      SS=wang_kSubset(nbrs,ord);
      ln_SS=size(SS,1);
      for si=1:ln_SS
        if iscell(SS)
	        S = SS{si};
        else
	        S=SS(si,:);
        end
       
        if feval(cond_indep, x, y, S, varargin{:})
	        % diagnostic
	        %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
	        %fprintf(': r = %6.4f\n', r);
	        %G(x,y) = 0;
	        %G(y,x) = 0;
            % x,y
            %pause
	        sep{x,y} = myunion(sep{x,y}, S);
	        sep{y,x} = myunion(sep{y,x}, S);
	        %break; % no need to check any more subsets 
            %bug 2 need check all subsets jbw sep 2. 2003
	    end
      end %end subset
    end  %end if
  end %end all ordered pairs
  ord = ord + 1;
end %end while

% Create the minimal pattern, step C
% i.e., the only directed edges are V structures.
pdag = G;
X=[]; Y=[];
[X, Y] = find(G);
% We want to generate all unique triples x,y,z
% This code generates x,y,z and z,y,x.
for i=1:length(X)
  x = X(i);
  y = Y(i);
  Z = find(G(y,:));
  ZZ = mysetdiff(Z, x);
  %added jbw
  for z=ZZ(:)'  
     %end added jbw
     if ( (G(x,z)==0 & G(z,x)==0) & (~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})) )
        %added wang try to force the order of direction
         if (  ((x<y & z<y) &(pdag(x,y)==1 & pdag(y,x)==1)) & (pdag(z,y)==1 & pdag(y,z)==1) )
           % fprintf('%d -> %d <- %d\n', x, y, z);
            pdag(x,y) = -1; pdag(y,x) = 0;
            pdag(z,y) = -1; pdag(y,z) = 0;
        end
    end
  end
end

% Convert the minimal pattern to a complete one,
% i.e., every directed edge in P is compelled
% (must be directed in all Markov equivalent models),
% and every undirected edge in P is reversible.
% We use the rules of Pearl (2000) p51 (derived in Meek (1995))
%step D

old_pdag = zeros(n);
iter = 0;
while ~isequal(pdag, old_pdag)
  iter = iter + 1;
  old_pdag = pdag;
  %apply meek's 4 rules to orignate the edge directions
  pdag=mgraph_meek4ruls(G,pdag);
end %end all
G=pdag;


if isdag==1
    %disp('Turning rest undirected edges ....');
    iter2 =0;
    while ~isempty(find(pdag==1))
        %turn rest of undirected edges to a->b then apply meek 4 rules
        iter2 = iter2 + 1;
        [X,Y]=find(pdag==1);
        for i=1:length(X)
            x=X(i);y=Y(i);
            %first turn x->y in pdag
            if ( ( ( x<y & pdag(x,y)==1) & pdag(y,x)==1 ) )
              %  fprintf('%d -> %d \n', x, y);
                pdag(x,y)=-1; pdag(y,x)=0;
            end
            %then apply meek's 4 rules again
            pdag=mgraph_meek4ruls(G,pdag);
        end
    end
end