wolffd@0: 
wolffd@0: % make undirected adjacency matrix of graph/tree
wolffd@0: % e.g.,
wolffd@0: %  1 
wolffd@0: % / \
wolffd@0: % 2  3
wolffd@0: T = zeros(3,3);
wolffd@0: T(1,2) = 1; T(2,1)=1;
wolffd@0: T(1,3)=1; T(3,1) = 1;
wolffd@0: 
wolffd@0: root = 1;
wolffd@0: [T, preorder, postorder] =  mk_rooted_tree(T, root);
wolffd@0: 
wolffd@0: % bottom up message passing leaves to root
wolffd@0: for n=postorder(:)'
wolffd@0:   for p	= parents(T, n)
wolffd@0:      % p is parent of n
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: % top down, root to leaves
wolffd@0: for n=preorder(:)'
wolffd@0:   for c= children(T,n)
wolffd@0:     % c is child of n
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%
wolffd@0: 
wolffd@0: function ps = parents(adj_mat, i)
wolffd@0: % PARENTS Return the list of parents of node i
wolffd@0: % ps = parents(adj_mat, i)
wolffd@0: 
wolffd@0: ps = find(adj_mat(:,i))';
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%%%
wolffd@0: 
wolffd@0: function cs = children(adj_mat, i, t)
wolffd@0: % CHILDREN Return the indices of a node's children in sorted order
wolffd@0: % c = children(adj_mat, i, t)
wolffd@0: %
wolffd@0: % t is an optional argument: if present, dag is assumed to be a 2-slice DBN
wolffd@0: 
wolffd@0: if nargin < 3 
wolffd@0:   cs = find(adj_mat(i,:));
wolffd@0: else
wolffd@0:   if t==1
wolffd@0:     cs = find(adj_mat(i,:));
wolffd@0:   else
wolffd@0:     ss = length(adj_mat)/2;
wolffd@0:     j = i+ss;
wolffd@0:     cs = find(adj_mat(j,:)) + (t-2)*ss;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%%%%
wolffd@0: 
wolffd@0: function [T, pre, post, cycle] = mk_rooted_tree(G, root)
wolffd@0: % MK_ROOTED_TREE Make a directed, rooted tree out of an undirected tree.
wolffd@0: % [T, pre, post, cycle] = mk_rooted_tree(G, root)
wolffd@0: 
wolffd@0: n = length(G);
wolffd@0: T = sparse(n,n); % not the same as T = sparse(n) !
wolffd@0: directed = 0;
wolffd@0: [d, pre, post, cycle, f, pred] = dfs(G, root, directed);
wolffd@0: for i=1:length(pred)
wolffd@0:   if pred(i)>0
wolffd@0:     T(pred(i),i)=1;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%%
wolffd@0: 
wolffd@0: function [d, pre, post, cycle, f, pred] = dfs(adj_mat, start, directed)
wolffd@0: % DFS Perform a depth-first search of the graph starting from 'start'.
wolffd@0: % [d, pre, post, cycle, f, pred] = dfs(adj_mat, start, directed)
wolffd@0: %
wolffd@0: % Input:
wolffd@0: % adj_mat(i,j)=1 iff i is connected to j.
wolffd@0: % start is the root vertex of the dfs tree; if [], all nodes are searched
wolffd@0: % directed = 1 if the graph is directed
wolffd@0: %
wolffd@0: % Output:
wolffd@0: % d(i) is the time at which node i is first discovered.
wolffd@0: % pre is a list of the nodes in the order in which they are first encountered (opened).
wolffd@0: % post is a list of the nodes in the order in which they are last encountered (closed).
wolffd@0: % 'cycle' is true iff a (directed) cycle is found.
wolffd@0: % f(i) is the time at which node i is finished.
wolffd@0: % pred(i) is the predecessor of i in the dfs tree.
wolffd@0: %
wolffd@0: % If the graph is a tree, preorder is parents before children,
wolffd@0: % and postorder is children before parents.
wolffd@0: % For a DAG, topological order = reverse(postorder).
wolffd@0: %
wolffd@0: % See Cormen, Leiserson and Rivest, "An intro. to algorithms" 1994, p478.
wolffd@0: 
wolffd@0: n = length(adj_mat);
wolffd@0: 
wolffd@0: global white gray black color
wolffd@0: white = 0; gray = 1; black = 2;
wolffd@0: color = white*ones(1,n);
wolffd@0: 
wolffd@0: global time_stamp
wolffd@0: time_stamp = 0;
wolffd@0: 
wolffd@0: global d f
wolffd@0: d = zeros(1,n);
wolffd@0: f = zeros(1,n);
wolffd@0: 
wolffd@0: global pred
wolffd@0: pred = zeros(1,n);
wolffd@0: 
wolffd@0: global cycle
wolffd@0: cycle = 0;
wolffd@0: 
wolffd@0: global pre post
wolffd@0: pre = [];
wolffd@0: post = [];
wolffd@0: 
wolffd@0: if ~isempty(start)
wolffd@0:   dfs_visit(start, adj_mat, directed);
wolffd@0: else
wolffd@0:   for u=1:n
wolffd@0:     if color(u)==white
wolffd@0:       dfs_visit(u, adj_mat, directed);
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%%
wolffd@0: 
wolffd@0: function dfs_visit(u, adj_mat, directed)
wolffd@0: 
wolffd@0: global white gray black color time_stamp d f pred cycle  pre post
wolffd@0: 
wolffd@0: pre = [pre u];
wolffd@0: color(u) = gray;
wolffd@0: time_stamp = time_stamp + 1;
wolffd@0: d(u) = time_stamp;
wolffd@0: if directed
wolffd@0:   ns = children(adj_mat, u);
wolffd@0: else
wolffd@0:   ns = neighbors(adj_mat, u);
wolffd@0:   ns = mysetdiff(ns, pred(u)); % don't go back to visit the guy who called you!
wolffd@0: end
wolffd@0: for v=ns(:)'
wolffd@0:   %fprintf('u=%d, v=%d, color(v)=%d\n', u, v, color(v))
wolffd@0:   switch color(v)
wolffd@0:     case white, % not visited v before (tree edge)
wolffd@0:      pred(v)=u;
wolffd@0:      dfs_visit(v, adj_mat, directed);
wolffd@0:    case gray, % back edge - v has been visited, but is still open
wolffd@0:     cycle = 1;
wolffd@0:     %fprintf('cycle: back edge from v=%d to u=%d\n', v, u);
wolffd@0:    case black, % v has been visited, but is closed
wolffd@0:     % no-op
wolffd@0:   end
wolffd@0: end
wolffd@0: color(u) = black;
wolffd@0: post = [post u];
wolffd@0: time_stamp = time_stamp + 1;
wolffd@0: f(u) = time_stamp;
wolffd@0: 
wolffd@0: