matthiasm@8: function [pdag, G] = learn_struct_pdag_ic_star(cond_indep, n, k, varargin)
matthiasm@8: % LEARN_STRUCT_PDAG_IC_STAR Learn a partially oriented DAG (pattern) with latent 
matthiasm@8: % variables using the IC* algorithm
matthiasm@8: % P = learn_struct_pdag_ic_star(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 either a latent variable L such that i <-L-> j 
matthiasm@8: % OR there is a directed edge from i->j.
matthiasm@8: % P(i,j) = -2 if there is a marked directed i-*>j edge.
matthiasm@8: % P(i,j) = P(j,i) = 1 if there is and undirected edge i--j
matthiasm@8: % P(i,j) = P(j,i) = 2 if there is a latent variable L such that i<-L->j.
matthiasm@8: %
matthiasm@8: % The IC* algorithm learns a latent structure associated with a set of observed 
matthiasm@8: % variables. 
matthiasm@8: % The latent structure revealed is the projection in which every latent variable is
matthiasm@8: % 1) a root node
matthiasm@8: % 2) linked to exactly two observed variables.
matthiasm@8: % Latent variables in the projection are represented using a bidirectional graph, 
matthiasm@8: % and thus remain implicit.
matthiasm@8: %
matthiasm@8: % See Pearl, "Causality: Models, Reasoning, and Inference", 2000, p52 for more details.
matthiasm@8: % Written by Tamar Kushnir, 2000
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:     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:       for si=1:length(SS)
matthiasm@8: 	S = SS{si};
matthiasm@8: 	if feval(cond_indep, x, y, S, varargin{:})  
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: % 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: % where y is a common neighbor to x and z
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:       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 using the following rules:
matthiasm@8: % Rule 1:
matthiasm@8: % if a and b are non-adjacent nodes with a common neighbor c,
matthiasm@8: % if a->c and not b->c then c-*>b (marked arrow).
matthiasm@8: % Rule 2:
matthiasm@8: % if a and b are adjacent and there is a directed path (marked links) from a to b
matthiasm@8: % then a->b (add arrowhead).
matthiasm@8: %Pearl (2000)
matthiasm@8: 
matthiasm@8: arrowin = [-1 -2 2];
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:   [X, Y] = find(pdag);
matthiasm@8:   for i=1:length(X)
matthiasm@8:     x = X(i);
matthiasm@8:     y = Y(i);
matthiasm@8:     Z = find(pdag(y,:));
matthiasm@8:     Z = mysetdiff(Z, x);
matthiasm@8:     for z=Z(:)'
matthiasm@8:       if G(x,z)==0 & ismember(pdag(x,y),arrowin) & ~ismember(pdag(z,y),arrowin)
matthiasm@8:         pdag(y,z) = -2; pdag(z,y) = 0;
matthiasm@8:       end
matthiasm@8:     end
matthiasm@8:   end
matthiasm@8:   % rule 2
matthiasm@8:   [X, Y] = find(G); 
matthiasm@8:   %check all adjacent nodes because if pdag(x,y) = -1 
matthiasm@8:   %and pdag(y,x) = 0 there could still be an bidirected edge between x & y.
matthiasm@8:   for i=1:length(X)
matthiasm@8:     x = X(i);
matthiasm@8:     y = Y(i);
matthiasm@8:     if ~ismember(pdag(x,y), arrowin) %x->y doesn't exist yet
matthiasm@8:       %find marked path from x to y
matthiasm@8:       add_arrow = marked_path(x,y,pdag);
matthiasm@8:       if add_arrow 
matthiasm@8:         if pdag(y,x)==-1 %bidirected edge
matthiasm@8:           pdag(x,y) = 2; pdag(y,x) = 2;
matthiasm@8:         else
matthiasm@8:           pdag(x,y) = -1;pdag(y,x) = 0;
matthiasm@8:         end
matthiasm@8:       end
matthiasm@8:     end
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: 
matthiasm@8: %%%%%%%%%%%%%
matthiasm@8: 
matthiasm@8: function t = marked_path(x,y,L)
matthiasm@8: % MARKED_PATH is a boolean function which returns 1 if a marked path 
matthiasm@8: % between nodes x and y exists in the partially directed latent structure L.
matthiasm@8: %
matthiasm@8: % t = marked_path(x,y,L)
matthiasm@8: %
matthiasm@8: % x and y are the starting and ending nodes in the path, respectively.
matthiasm@8: % L is a latent structure (partially directed graph with possible latent variables).
matthiasm@8: %
matthiasm@8: % Rule 2 of IC* algorithm (see Pearl, 2000)
matthiasm@8: 
matthiasm@8: t=0;
matthiasm@8: 
matthiasm@8: %find set of marked links from x
matthiasm@8: marked = find(L(x,:)==-2);
matthiasm@8: if ismember(y,marked)
matthiasm@8:   t=1; %marked path found
matthiasm@8: else
matthiasm@8:   for m=marked(:)'
matthiasm@8:     t = marked_path(m,y,L);
matthiasm@8:     if t==1
matthiasm@8:       break; %stop when marked path found
matthiasm@8:     end
matthiasm@8:   end
matthiasm@8: end