function [pdag, G] = wang_learn_struct_pdag_pc(cond_indep, n, k, 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 = ones(n,n);
G=setdiag(G,0);
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]);
if (length(nbrs) >= ord & G(x,y) ~= 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;
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] = 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 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x}))
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
% 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;
% rule 1 step D.1
[A,B] = find(pdag==-1); % a -> b
for i=1:length(A)
a = A(i); b = B(i);
C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
if ~isempty(C)
pdag(b,C) = -1; pdag(C,b) = 0;
fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
end
end
clear A B C
% rule 2 Step D.2 , Rule 2. 3 may unstable according to Tetrad4.2
[A,B] = find(pdag==1); % unoriented a-b edge
for i=1:length(A)
a = A(i); b = B(i);
if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
pdag(a,b) = -1; pdag(b,a) = 0;
fprintf('rule 2: %d -> %d\n', a, b);
end
end
clear A B
% % rule 3
[A,B] = find(pdag==1); % a-b
for i=1:length(A)
clear G2
a = A(i); b = B(i);
C = find( (G(a,:)==1) & (pdag(:,b)==-1)' );
% C contains nodes c s.t. a-c->ba
G2 = setdiag(G(C, C), 1);
if any(G2(:)==0) % there are 2 different non adjacent elements of C
pdag(a,b) = -1; pdag(b,a) = 0;
fprintf('rule 3: %d -> %d\n', a, b);
end
end
%rule 4 from Meek(95), this rule may unstable
%added wang 10.2003
%clear A B C1 C2 G2
%close all
%figure
%pdag=[0 1 0 1; 1 0 0 0 ; -1 -1 0 0 ; 1 0 -1 0]
%draw_graph(abs(pdag));
%G=[0 1 1 1; 1 0 1 0 ; 1 1 0 1 ; 1 0 1 0]
[A, B]=find((G)==1);
for i=1:length(A)
clear C1 C2 C12 C1_pa C1_ch G2
a=A(i); b=B(i);
C1_pa =find( (pdag(:,a)==-1)' ); %a's parients directed neighbor
C1_ch=find((pdag(a,:)==-1)); %a's child directed neoghbor
C1=find( G(a,:)==1) ;%a's undirected neighbor
C2=find( G(b,:)==1); %b's undirected neighbor
C12=intersect(C1,C2); % a.b common neighbor
G2 = setdiag(G(C12, C12), 1);
%[~isempty(C12), any(G2(:)==0) ,~isempty(intersect(C1_ch,C12)), ~isempty(intersect(C1_pa,C12))]
%pause
if ( (~isempty(C12) & any(G2(:)==0) ) & ( ~isempty(intersect(C1_ch,C12)) & ~isempty(intersect(C1_pa,C12)) ) )
%non adjacent element in C12 , orign a's child with b as same direction of a's parients
need_origned_node=intersect(C1_ch,C2);
if ~isempty(need_origned_node)
for jj=1:length(need_origned_node)
if (pdag(b,need_origned_node(jj))==1 & pdag(need_origned_node(jj),b)==1)
%check once more
pdag(b,need_origned_node)=-1; pdag(need_origned_node,b) = 0;
fprintf('rule 4: %d -> %d \n', b,need_origned_node(jj));
% a ,b
% need_origned_node
% warndlg('Error in Rule 4; PLEAE CHECK');
% pdag
% G
% break;
end
end %end for jj
end
end
end
%figure
%draw_graph(abs(pdag))
end %end all
