matthiasm@8: function [pdag, G] = learn_struct_pdag_pc(cond_indep, n, k, varargin)
matthiasm@8: % LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm
matthiasm@8: % P = learn_struct_pdag_pc(cond_indep, n, k, ...)
matthiasm@8: %
matthiasm@8: % n is the number of nodes.
matthiasm@8: % k is an optional upper bound on the fan-in (default: n)
matthiasm@8: % cond_indep is a boolean function that will be called as follows:
matthiasm@8: %   feval(cond_indep, x, y, S, ...)
matthiasm@8: % where x and y are nodes, and S is a set of nodes (positive integers),
matthiasm@8: % and ... are any optional parameters passed to this function.
matthiasm@8: %
matthiasm@8: % The output P is an adjacency matrix, in which
matthiasm@8: % P(i,j) = -1 if there is an i->j edge.
matthiasm@8: % P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
matthiasm@8: %
matthiasm@8: % The PC algorithm does structure learning assuming all variables are observed.
matthiasm@8: % See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
matthiasm@8: % This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
matthiasm@8: % but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).
matthiasm@8: 
matthiasm@8:   
matthiasm@8: sep = cell(n,n);
matthiasm@8: ord = 0;
matthiasm@8: done = 0;
matthiasm@8: G = ones(n,n);
matthiasm@8: G=setdiag(G,0);
matthiasm@8: while ~done
matthiasm@8:   done = 1;
matthiasm@8:   [X,Y] = find(G); 
matthiasm@8:   for i=1:length(X)
matthiasm@8:     x = X(i); y = Y(i);
matthiasm@8:     %nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
matthiasm@8:     nbrs = mysetdiff(neighbors(G, y), x);  % bug fix by Raanan Yehezkel <raanany@ee.bgu.ac.il> 6/27/04
matthiasm@8:     if length(nbrs) >= ord & G(x,y) ~= 0
matthiasm@8:       done = 0;
matthiasm@8:       %SS = subsets(nbrs, ord, ord); % all subsets of size ord
matthiasm@8:       SS = subsets1(nbrs, ord);
matthiasm@8:       for si=1:length(SS)
matthiasm@8: 	S = SS{si};
matthiasm@8: 	if feval(cond_indep, x, y, S, varargin{:})
matthiasm@8: 	  %if isempty(S)
matthiasm@8:  	  %  fprintf('%d indep of %d ', x, y);
matthiasm@8: 	  %else
matthiasm@8:  	  %  fprintf('%d indep of %d given ', x, y); fprintf('%d ', S);
matthiasm@8:  	  %end
matthiasm@8:  	  %fprintf('\n');
matthiasm@8: 	  
matthiasm@8: 	  % diagnostic
matthiasm@8: 	  %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
matthiasm@8: 	  %fprintf(': r = %6.4f\n', r);
matthiasm@8: 	  
matthiasm@8: 	  G(x,y) = 0;
matthiasm@8: 	  G(y,x) = 0;
matthiasm@8: 	  sep{x,y} = myunion(sep{x,y}, S);
matthiasm@8: 	  sep{y,x} = myunion(sep{y,x}, S);
matthiasm@8: 	  break; % no need to check any more subsets 
matthiasm@8: 	end
matthiasm@8:       end
matthiasm@8:     end 
matthiasm@8:   end
matthiasm@8:   ord = ord + 1;
matthiasm@8: end
matthiasm@8: 
matthiasm@8: 
matthiasm@8: % Create the minimal pattern,
matthiasm@8: % i.e., the only directed edges are V structures.
matthiasm@8: pdag = G;
matthiasm@8: [X, Y] = find(G);
matthiasm@8: % We want to generate all unique triples x,y,z
matthiasm@8: % This code generates x,y,z and z,y,x.
matthiasm@8: for i=1:length(X)
matthiasm@8:   x = X(i);
matthiasm@8:   y = Y(i);
matthiasm@8:   Z = find(G(y,:));
matthiasm@8:   Z = mysetdiff(Z, x);
matthiasm@8:   for z=Z(:)'
matthiasm@8:     if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
matthiasm@8:       %fprintf('%d -> %d <- %d\n', x, y, z);
matthiasm@8:       pdag(x,y) = -1; pdag(y,x) = 0;
matthiasm@8:       pdag(z,y) = -1; pdag(y,z) = 0;
matthiasm@8:     end
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: % Convert the minimal pattern to a complete one,
matthiasm@8: % i.e., every directed edge in P is compelled
matthiasm@8: % (must be directed in all Markov equivalent models),
matthiasm@8: % and every undirected edge in P is reversible.
matthiasm@8: % We use the rules of Pearl (2000) p51 (derived in Meek (1995))
matthiasm@8: 
matthiasm@8: old_pdag = zeros(n);
matthiasm@8: iter = 0;
matthiasm@8: while ~isequal(pdag, old_pdag)
matthiasm@8:   iter = iter + 1;
matthiasm@8:   old_pdag = pdag;
matthiasm@8:   % rule 1
matthiasm@8:   [A,B] = find(pdag==-1); % a -> b
matthiasm@8:   for i=1:length(A)
matthiasm@8:     a = A(i); b = B(i);
matthiasm@8:     C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
matthiasm@8:     if ~isempty(C)
matthiasm@8:       pdag(b,C) = -1; pdag(C,b) = 0;
matthiasm@8:       %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
matthiasm@8:     end
matthiasm@8:   end
matthiasm@8:   % rule 2
matthiasm@8:   [A,B] = find(pdag==1); % unoriented a-b edge
matthiasm@8:   for i=1:length(A)
matthiasm@8:     a = A(i); b = B(i);
matthiasm@8:     if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
matthiasm@8:       pdag(a,b) = -1; pdag(b,a) = 0;
matthiasm@8:       %fprintf('rule 2: %d -> %d\n', a, b);
matthiasm@8:     end
matthiasm@8:   end
matthiasm@8:   % rule 3
matthiasm@8:   [A,B] = find(pdag==1); % a-b
matthiasm@8:   for i=1:length(A)
matthiasm@8:     a = A(i); b = B(i);
matthiasm@8:     C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
matthiasm@8:     % C contains nodes c s.t. a-c->ba
matthiasm@8:     G2 = setdiag(G(C, C), 1);
matthiasm@8:     if any(G2(:)==0) % there are 2 different non adjacent elements of C
matthiasm@8:       pdag(a,b) = -1; pdag(b,a) = 0;
matthiasm@8:       %fprintf('rule 3: %d -> %d\n', a, b);
matthiasm@8:     end
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: