diff options
Diffstat (limited to 'sourcecodes/bnt-master/graph')
61 files changed, 2993 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/graph/CVS/Entries b/sourcecodes/bnt-master/graph/CVS/Entries new file mode 100644 index 00000000..cf3baf2c --- /dev/null +++ b/sourcecodes/bnt-master/graph/CVS/Entries @@ -0,0 +1,49 @@ +/README.txt/1.1.1.1/Mon Aug 29 17:44:02 2005// +/acyclic.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/assignEdgeNums.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/best_first_elim_order.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/check_jtree_property.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/check_triangulated.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/children.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/cliques_to_jtree.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/cliques_to_strong_jtree.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/connected_graph.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/dag_to_essential_graph.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/dfs.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/dfs_test.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/dijkstra.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/family.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/graph_separated.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/graph_to_jtree.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/min_subtree_con_nodes.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/minimum_spanning_tree.m/1.1.1.1/Mon Aug 29 17:44:04 2005// +/minspan.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_2D_lattice.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_2D_lattice_slow.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_adj_mat.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_adjmat_chain.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_all_dags.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_nbrs_of_dag.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_nbrs_of_digraph.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_nbrs_of_digraph_broken.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_nbrs_of_digraph_not_vectorized.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_rnd_dag.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_rnd_dag_given_edge_prob.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_rooted_tree.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/mk_undirected.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/moralize.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/neighbors.m/1.1.1.1/Mon Aug 29 17:44:06 2005// +/parents.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/pred2path.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/reachability_graph.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/scc.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/strong_elim_order.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/test.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/test_strong_root.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/topological_sort.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/trees.txt/1.1.1.1/Mon Aug 29 17:44:08 2005// +/triangulate.c/1.1.1.1/Mon Aug 29 17:44:08 2005// +/triangulate.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/triangulate_2Dlattice_demo.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +/triangulate_test.m/1.1.1.1/Mon Aug 29 17:44:08 2005// +D diff --git a/sourcecodes/bnt-master/graph/CVS/Entries.Log b/sourcecodes/bnt-master/graph/CVS/Entries.Log new file mode 100644 index 00000000..24f16336 --- /dev/null +++ b/sourcecodes/bnt-master/graph/CVS/Entries.Log @@ -0,0 +1 @@ +A D/Old//// diff --git a/sourcecodes/bnt-master/graph/CVS/Repository b/sourcecodes/bnt-master/graph/CVS/Repository new file mode 100644 index 00000000..a18677ae --- /dev/null +++ b/sourcecodes/bnt-master/graph/CVS/Repository @@ -0,0 +1 @@ +FullBNT/graph diff --git a/sourcecodes/bnt-master/graph/CVS/Root b/sourcecodes/bnt-master/graph/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/graph/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/graph/Old/CVS/Entries b/sourcecodes/bnt-master/graph/Old/CVS/Entries new file mode 100644 index 00000000..c45011ba --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/CVS/Entries @@ -0,0 +1,6 @@ +/best_first_elim_order.m/1.1.1.1/Mon Aug 29 17:44:02 2005// +/dag_to_jtree.m/1.1.1.1/Mon Aug 29 17:44:02 2005// +/dfs.m/1.1.1.1/Mon Aug 29 17:44:02 2005// +/dsep_test.m/1.1.1.1/Mon Aug 29 17:44:02 2005// +/mk_2D_lattice_slow.m/1.1.1.1/Mon Aug 29 17:44:02 2005// +D diff --git a/sourcecodes/bnt-master/graph/Old/CVS/Repository b/sourcecodes/bnt-master/graph/Old/CVS/Repository new file mode 100644 index 00000000..f3b2d81f --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/graph/Old diff --git a/sourcecodes/bnt-master/graph/Old/CVS/Root b/sourcecodes/bnt-master/graph/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/graph/Old/best_first_elim_order.m b/sourcecodes/bnt-master/graph/Old/best_first_elim_order.m new file mode 100644 index 00000000..82f1236b --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/best_first_elim_order.m @@ -0,0 +1,64 @@ +function order = best_first_elim_order(G, node_sizes, stage) +% BEST_FIRST_ELIM_ORDER Greedily search for an optimal elimination order. +% order = best_first_elim_order(moral_graph, node_sizes) +% +% Find an order in which to eliminate nodes from the graph in such a way as to try and minimize the +% weight of the resulting triangulated graph. The weight of a graph is the sum of the weights of each +% of its cliques; the weight of a clique is the product of the weights of each of its members; the +% weight of a node is the number of values it can take on. +% +% Since this is an NP-hard problem, we use the following greedy heuristic: +% at each step, eliminate that node which will result in the addition of the least +% number of fill-in edges, breaking ties by choosing the node that induces the lighest clique. +% For details, see +% - Kjaerulff, "Triangulation of graphs -- algorithms giving small total state space", +% Univ. Aalborg tech report, 1990 (www.cs.auc.dk/~uk) +% - C. Huang and A. Darwiche, "Inference in Belief Networks: A procedural guide", +% Intl. J. Approx. Reasoning, 11, 1994 +% + +% Warning: This code is pretty old and could probably be made faster. + +n = length(G); +if nargin < 3, stage = { 1:n }; end % no constraints + +% For long DBNs, it may be useful to eliminate all the nodes in slice t before slice t+1. +% This will ensure that the jtree has a repeating structure (at least away from both edges). +% This is why we have stages. +% See the discussion of splicing jtrees on p68 of +% Geoff Zweig's PhD thesis, Dept. Comp. Sci., UC Berkeley, 1998. +% This constraint can increase the clique size significantly. + +MG = G; % copy the original graph +uneliminated = ones(1,n); +order = zeros(1,n); +t = 1; % Counts which time slice we are on +for i=1:n + U = find(uneliminated); + valid = myintersect(U, stage{t}); + % Choose the best node from the set of valid candidates + score1 = zeros(1,length(valid)); + score2 = zeros(1,length(valid)); + for j=1:length(valid) + k = valid(j); + ns = myintersect(neighbors(G, k), U); + l = length(ns); + M = MG(ns,ns); + score1(j) = l^2 - sum(M(:)); % num. added edges + score2(j) = prod(node_sizes([k ns])); % weight of clique + end + j1s = find(score1==min(score1)); + j = j1s(argmin(score2(j1s))); + k = valid(j); + uneliminated(k) = 0; + order(i) = k; + ns = myintersect(neighbors(G, k), U); + if ~isempty(ns) + G(ns,ns) = 1; + G = setdiag(G,0); + end + if ~any(logical(uneliminated(stage{t}))) % are we allowed to the next slice? + t = t + 1; + end +end + diff --git a/sourcecodes/bnt-master/graph/Old/dag_to_jtree.m b/sourcecodes/bnt-master/graph/Old/dag_to_jtree.m new file mode 100644 index 00000000..23c58430 --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/dag_to_jtree.m @@ -0,0 +1,55 @@ +function [jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges, strong] = ... + dag_to_jtree(dag, node_sizes, partial_order, stages, clusters) +% DAG_TO_JTREE Moralize and triangulate a DAG, and make a junction tree from its cliques. +% [jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges, strong] = ... +% dag_to_jtree(dag, node_sizes, partial_order, stages, clusters) +% +% Input: +% dag(i,j) +% jtree(i,j) = 1 iff there is an arc between clique i and clique j +% root = the root clique +% cliques{i} = the nodes in clique i +% B(i,j) = 1 iff node j occurs in clique i +% w(i) = weight of clique i + +N = length(bnet.dag); +if nargin < 2, obs_nodes = []; end +if nargin < 3, stages = { 1:N }; end +if nargin < 4, clusters = {}; end + +[MG, moral_edges] = moralize(bnet.dag); + +% Add extra arcs between nodes in each cluster to ensure they occur in the same clique +for i=1:length(clusters) + c = clusters{i}; + MG(c,c) = 1; +end +MG = setdiag(MG, 0); + +% Find an optimal elimination ordering (NP-hard problem!) +ns = bnet.node_sizes(:); +ns(obs_nodes) = 1; % observed nodes have only 1 possible value +partial_order = determine_elim_constraints(bnet, obs_nodes); + +if isempty(partial_order) + strong = 0; + elim_order = best_first_elim_order(MG, ns, stages); +else + strong = 1; + elim_order = strong_elim_order(MG, ns, partial_order); +end + +[MTG, cliques, fill_in_edges] = triangulate(MG, elim_order); + +% Connect the cliques up into a jtree, +[jtree, root, B, w] = cliques_to_jtree(cliques, ns); + +if 0 + disp('testing dag to jtree'); + % Find the cliques containing each node, and check they form a connected subtree + clqs_con_node = cell(1,N); + for i=1:N + clqs_con_node{i} = find(B(:,i))'; + end + check_jtree_property(clqs_con_node, jtree); +end diff --git a/sourcecodes/bnt-master/graph/Old/dfs.m b/sourcecodes/bnt-master/graph/Old/dfs.m new file mode 100644 index 00000000..0c7f46ab --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/dfs.m @@ -0,0 +1,84 @@ +function [d, pre, post, height, cycle, pred] = dfs(adj_mat, start, directed) +% DFS Perform a depth-first search of the graph starting from 'start'. +% [d, pre, post, height, cycle, pred] = dfs(adj_mat, start, directed) +% +% d(i) is the time at which node i is first discovered. +% pre is a listing of the nodes in the order in which they are first encountered (opened). +% post is a listing of the nodes in the order in which they are last encountered (closed). +% A node is last encountered once we have explored all of its neighbors. +% If the graph is directed, i's neighbors are its children. +% If the graph is a tree, preorder is parents before children, and +% postorder is children before parents. +% For a DAG, topological order = reverse(postorder). +% height(i) is the height (distance) of node i from the start. +% 'cycle' is true iff a (directed) cycle is found. +% pred(i) is the parent of i in the dfs tree rooted at start. +% See Cormen, Leiserson and Rivest, "An intro. to algorithms" 1994, p478. + +% We can detect undirected cycles by checking if we are about to visit a node n which we have +% already visited. To detect *directed* cycles, we need to know if n has been closed or is still open. +% For example (where arcs are directed down) +% 1 2 +% \ / +% 3 +% Assume we visit 1, 3 and then 2 in order. The fact that a child of 2 (namely, 3) has +% already been visited is okay, because 3 has been closed. +% The algorithms in Aho, Hopcroft and Ullman, and Sedgewick, do not detect directed cycles. + +n = length(adj_mat); + +global white gray black +white = 0; gray = 1; black = 2; + +color = white*ones(1,n); +d = zeros(1,n); +height = zeros(1,n); +pred = zeros(1,n); +pre = []; +post = []; +cycle = 0; +global count +count = 0; +h = 0; +[d, pre, post, height, cycle, color, pred] = ... + dfs2(adj_mat, start, directed, h, d, pre, post, height, cycle, color, pred); + + + +%%%%%%%%%% + +function [d, pre, post, height, cycle, color, pred] = ... + dfs2(adj_mat, i, directed, h, d, pre, post, height, cycle, color, pred) + +global count +global white gray black + +color(i) = gray; +count = count + 1; +d(i) = count; +pre = [pre i]; +height(i) = h; +if directed + ns = children(adj_mat, i); +else + ns = neighbors(adj_mat, i); +end +for j=1:length(ns) + n=ns(j); + if ~directed & n==pred(i) % don't go back up the edge you just came down + % continue + else + if color(n) == gray % going back to a non-closed vertex via a new edge + %fprintf(1, 'cycle from %d to %d\n', i, n); + cycle = 1; + end + if color(n) == white % not visited n before + pred(n)=i; + [d, pre, post, height, cycle, color, pred] = ... + dfs2(adj_mat, n, directed, h+1, d, pre, post, height, cycle, color, pred); + end + end +end +color(i) = black; +post = [post i]; + diff --git a/sourcecodes/bnt-master/graph/Old/dsep_test.m b/sourcecodes/bnt-master/graph/Old/dsep_test.m new file mode 100644 index 00000000..f0f14638 --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/dsep_test.m @@ -0,0 +1,15 @@ + +% Cowell et al p72 +G = zeros(10); +G(1,2)=1; +G(2,3)=1; +G(3,7)=1; +G(4,[5 8])=1; +G(5,6)=1; +G(6,7)=1; +G(7,[9 10])=1; +G(8,9)=1; + +dsep(1, 4, [5 7], G) +dsep(1, 4, [7], G) +dsep(1, 4, [10 5], G) diff --git a/sourcecodes/bnt-master/graph/Old/mk_2D_lattice_slow.m b/sourcecodes/bnt-master/graph/Old/mk_2D_lattice_slow.m new file mode 100644 index 00000000..e5f595b7 --- /dev/null +++ b/sourcecodes/bnt-master/graph/Old/mk_2D_lattice_slow.m @@ -0,0 +1,121 @@ +function G = mk_2D_lattice(nrows, ncols, wrap_around) +% MK_2D_LATTICE Return adjacency matrix for 4-nearest neighbor connected 2D lattice +% G = mk_2D_lattice(nrows, ncols, wrap_around) +% G(k1, k2) = 1 iff k1=(i1,j1) is connected to k2=(i2,j2) +% +% If wrap_around = 1, we use toroidal boundary conditions (default = 0) +% +% Nodes are assumed numbered as in the following 3x3 lattice +% 1 4 7 +% 2 5 8 +% 3 6 9 +% +% e.g., G = mk_2D_lattice(3, 3, 0) returns +% 0 1 0 1 0 0 0 0 0 +% 1 0 1 0 1 0 0 0 0 +% 0 1 0 0 0 1 0 0 0 +% 1 0 0 0 1 0 1 0 0 +% 0 1 0 1 0 1 0 1 0 +% 0 0 1 0 1 0 0 0 1 +% 0 0 0 1 0 0 0 1 0 +% 0 0 0 0 1 0 1 0 1 +% 0 0 0 0 0 1 0 1 0 +% so find(G(5,:)) = [2 4 6 8] +% but find(G(1,:)) = [2 4] +% +% Using wrap around, G = mk_2D_lattice(3, 3, 1), we get +% 0 1 1 1 0 0 1 0 0 +% 1 0 1 0 1 0 0 1 0 +% 1 1 0 0 0 1 0 0 1 +% 1 0 0 0 1 1 1 0 0 +% 0 1 0 1 0 1 0 1 0 +% 0 0 1 1 1 0 0 0 1 +% 1 0 0 1 0 0 0 1 1 +% 0 1 0 0 1 0 1 0 1 +% 0 0 1 0 0 1 1 1 0 +% so find(G(5,:)) = [2 4 6 8] +% and find(G(1,:)) = [2 3 4 7] + +if nargin < 3, wrap_around = 0; end + +% M contains the number of each cell e.g. +% 1 4 7 +% 2 5 8 +% 3 6 9 +% North neighbors (assuming wrap around) are +% 3 6 9 +% 1 4 7 +% 2 5 8 +% Without wrap around, they are +% 1 4 7 +% 1 4 7 +% 2 5 8 +% The first row is arbitrary, since pixels at the top have no north neighbor. + +npixels = nrows*ncols; + +N = 1; E = 2; S = 3; W = 4; +if wrap_around + rows{N} = [nrows 1:nrows-1]; cols{N} = 1:ncols; + rows{E} = 1:nrows; cols{E} = [2:ncols 1]; + rows{S} = [2:nrows 1]; cols{S} = 1:ncols; + rows{W} = 1:nrows; cols{W} = [ncols 1:ncols-1]; +else + rows{N} = [1 1:nrows-1]; cols{N} = 1:ncols; + rows{E} = 1:nrows; cols{E} = [2:ncols 2]; + rows{S} = [2:nrows 2]; cols{S} = 1:ncols; + rows{W} = 1:nrows; cols{W} = [1 1:ncols-1]; +end + +M = reshape(1:npixels, [nrows ncols]); +nbrs = cell(1, 4); +for i=1:4 + nbrs{i} = M(rows{i}, cols{i}); +end + + +G = zeros(npixels, npixels); +if wrap_around + for i=1:4 + if 0 + % naive + for p=1:npixels + G(p, nbrs{i}(p)) = 1; + end + else + % vectorized + ndx2 = sub2ind([npixels npixels], 1:npixels, nbrs{i}(:)'); + G(ndx2) = 1; + end + end +else + i = N; + mask = ones(nrows, ncols); + mask(1,:) = 0; % pixels in row 1 have no nbr to the north + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; + + i = E; + mask = ones(nrows, ncols); + mask(:,ncols) = 0; + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; + + i = S; + mask = ones(nrows, ncols); + mask(nrows,:)=0; + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; + + i = W; + mask = ones(nrows, ncols); + mask(:,1)=0; + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; +end + +G = setdiag(G, 0); diff --git a/sourcecodes/bnt-master/graph/README.txt b/sourcecodes/bnt-master/graph/README.txt new file mode 100644 index 00000000..917cd9ae --- /dev/null +++ b/sourcecodes/bnt-master/graph/README.txt @@ -0,0 +1,3 @@ +Graph theory toolbox. +All functions written by Kevin Murphy (1998) unless otherwise noted. +The latest version is available from http://www.ai.mit.edu/~murphyk/Software/index.html diff --git a/sourcecodes/bnt-master/graph/acyclic.m b/sourcecodes/bnt-master/graph/acyclic.m new file mode 100644 index 00000000..0b497d73 --- /dev/null +++ b/sourcecodes/bnt-master/graph/acyclic.m @@ -0,0 +1,23 @@ +function b = acyclic(adj_mat, directed) +% ACYCLIC Returns true iff the graph has no (directed) cycles. +% b = acyclic(adj_mat, directed) + +adj_mat = double(adj_mat); +if nargin < 2, directed = 1; end + +% e.g., G = +% 1 -> 3 +% | +% v +% 2 <- 4 +% In this case, 1->2 in the transitive closure, but 1 cannot get to itself. +% If G was undirected, 1 could get to itself, but this graph is not cyclic. +% So we cannot use the closure test in the undirected case. + +if directed + R = reachability_graph(adj_mat); + b = ~any(diag(R)==1); +else + [d, pre, post, cycle] = dfs(adj_mat,[],directed); + b = ~cycle; +end diff --git a/sourcecodes/bnt-master/graph/assignEdgeNums.m b/sourcecodes/bnt-master/graph/assignEdgeNums.m new file mode 100644 index 00000000..5d9a2a0d --- /dev/null +++ b/sourcecodes/bnt-master/graph/assignEdgeNums.m @@ -0,0 +1,32 @@ +function [edge_id, nedges] = assignEdgeNums(adj_mat) +% give each edge a unique number +% we number (i,j) for j>i first, in row, column order. +% Then we number the reverse links + +nnodes = length(adj_mat); +edge_id = zeros(nnodes); +e = 1; +for i=1:nnodes + for j=i+1:nnodes + if adj_mat(i,j) + edge_id(i,j) = e; + e = e+1; + end + end +end + +nedges = e-1; +tmp = edge_id; +ndx = find(tmp); +tmp(ndx) = tmp(ndx)+nedges; +edge_id = edge_id + triu(tmp)'; + + +if 0 +ndx = find(adj_mat); +nedges = length(ndx); +nnodes = length(adj_mat); +edge_id = zeros(1, nnodes*nnodes); +edge_id(ndx) = 1:nedges; +edge_id = reshape(edge_id, nnodes, nnodes); +end diff --git a/sourcecodes/bnt-master/graph/best_first_elim_order.m b/sourcecodes/bnt-master/graph/best_first_elim_order.m new file mode 100644 index 00000000..2a808c16 --- /dev/null +++ b/sourcecodes/bnt-master/graph/best_first_elim_order.m @@ -0,0 +1,68 @@ +function order = best_first_elim_order(G, node_sizes, stage) +% BEST_FIRST_ELIM_ORDER Greedily search for an optimal elimination order. +% order = best_first_elim_order(moral_graph, node_sizes) +% +% Find an order in which to eliminate nodes from the graph in such a way as to try and minimize the +% weight of the resulting triangulated graph. The weight of a graph is the sum of the weights of each +% of its cliques; the weight of a clique is the product of the weights of each of its members; the +% weight of a node is the number of values it can take on. +% +% Since this is an NP-hard problem, we use the following greedy heuristic: +% at each step, eliminate that node which will result in the addition of the least +% number of fill-in edges, breaking ties by choosing the node that induces the lighest clique. +% For details, see +% - Kjaerulff, "Triangulation of graphs -- algorithms giving small total state space", +% Univ. Aalborg tech report, 1990 (www.cs.auc.dk/~uk) +% - C. Huang and A. Darwiche, "Inference in Belief Networks: A procedural guide", +% Intl. J. Approx. Reasoning, 11, 1994 +% + +% Warning: This code is pretty old and could probably be made faster. + +n = length(G); +if nargin < 3, stage = { 1:n }; end % no constraints + +% For long DBNs, it may be useful to eliminate all the nodes in slice t before slice t+1. +% This will ensure that the jtree has a repeating structure (at least away from both edges). +% This is why we have stages. +% See the discussion of splicing jtrees on p68 of +% Geoff Zweig's PhD thesis, Dept. Comp. Sci., UC Berkeley, 1998. +% This constraint can increase the clique size significantly. + +MG = G; % copy the original graph +uneliminated = ones(1,n); +order = zeros(1,n); +t = 1; % Counts which time slice we are on +for i=1:n + U = find(uneliminated); + valid = myintersect(U, stage{t}); + % Choose the best node from the set of valid candidates + min_fill = zeros(1,length(valid)); + min_weight = zeros(1,length(valid)); + for j=1:length(valid) + k = valid(j); + nbrs = myintersect(neighbors(G, k), U); + l = length(nbrs); + M = MG(nbrs,nbrs); + min_fill(j) = l^2 - sum(M(:)); % num. added edges + min_weight(j) = prod(node_sizes([k nbrs])); % weight of clique + end + lightest_nbrs = find(min_weight==min(min_weight)); + % break ties using min-fill heuristic + best_nbr_ndx = argmin(min_fill(lightest_nbrs)); + j = lightest_nbrs(best_nbr_ndx); % we will eliminate the j'th element of valid + %j1s = find(score1==min(score1)); + %j = j1s(argmin(score2(j1s))); + k = valid(j); + uneliminated(k) = 0; + order(i) = k; + ns = myintersect(neighbors(G, k), U); + if ~isempty(ns) + G(ns,ns) = 1; + G = setdiag(G,0); + end + if ~any(logical(uneliminated(stage{t}))) % are we allowed to the next slice? + t = t + 1; + end +end + diff --git a/sourcecodes/bnt-master/graph/check_jtree_property.m b/sourcecodes/bnt-master/graph/check_jtree_property.m new file mode 100644 index 00000000..3986c05e --- /dev/null +++ b/sourcecodes/bnt-master/graph/check_jtree_property.m @@ -0,0 +1,18 @@ +function check_jtree_property(cliques_containing_node, jtree) +% CHECK_JTREE_PROPERTY Raise an error if the graph does not satisfy the join tree property. +% check_jtree_property(cliques_containing_node, jtree_adj_mat) +% +% The join tree property says: +% For each node n in the dag, compute the node-induced subgraph G by looking at all the cliques +% that contain n, and make sure G forms a connected graph. +% This ensures that local propagation leads to global consistency. + +num_bn_nodes = length(cliques_containing_node); +directed = 0; +for i=1:num_bn_nodes + cs = cliques_containing_node{i}; + G = jtree(cs,cs); + if ~connected_graph(G, directed) + error(['node ' num2str(i) ' violates jtree property']); + end +end diff --git a/sourcecodes/bnt-master/graph/check_triangulated.m b/sourcecodes/bnt-master/graph/check_triangulated.m new file mode 100644 index 00000000..99bcc9f1 --- /dev/null +++ b/sourcecodes/bnt-master/graph/check_triangulated.m @@ -0,0 +1,34 @@ +function [triangulated, order] = check_triangulated(G) +% CHECK_TRIANGULATED Return 1 if G is a triangulated (chordal) graph, 0 otherwise. +% [triangulated, order] = check_triangulated(G) +% +% A numbering alpha is perfect if Nbrs(alpha(i)) intersect {alpha(1)...alpha(i-1)} is complete. +% A graph is triangulated iff it has a perfect numbering. +% The Maximum Cardinality Search algorithm will create such a perfect numbering if possible. +% See Golumbic, "Algorithmic Graph Theory and Perfect Graphs", Cambridge Univ. Press, 1985, p85. +% or Castillo, Gutierrez and Hadi, "Expert systems and probabilistic network models", Springer 1997, p134. + + +G = setdiag(G, 1); +n = length(G); +order = zeros(1,n); +triangulated = 1; +numbered = [1]; +order(1) = 1; +for i=2:n + U = mysetdiff(1:n, numbered); % unnumbered nodes + score = zeros(1, length(U)); + for ui=1:length(U) + u = U(ui); + score(ui) = length(myintersect(neighbors(G, u), numbered)); + end + u = U(argmax(score)); + numbered = [numbered u]; + order(i) = u; + nns = myintersect(neighbors(G,u), order(1:i-1)); % numbered neighbors + if ~isequal(G(nns,nns), ones(length(nns))) % ~complete(G(nns,nns)) + triangulated = 0; + break; + end +end + diff --git a/sourcecodes/bnt-master/graph/children.m b/sourcecodes/bnt-master/graph/children.m new file mode 100644 index 00000000..d79701b7 --- /dev/null +++ b/sourcecodes/bnt-master/graph/children.m @@ -0,0 +1,17 @@ +function cs = children(adj_mat, i, t) +% CHILDREN Return the indices of a node's children in sorted order +% c = children(adj_mat, i, t) +% +% t is an optional argument: if present, dag is assumed to be a 2-slice DBN + +if nargin < 3 + cs = find(adj_mat(i,:)); +else + if t==1 + cs = find(adj_mat(i,:)); + else + ss = length(adj_mat)/2; + j = i+ss; + cs = find(adj_mat(j,:)) + (t-2)*ss; + end +end diff --git a/sourcecodes/bnt-master/graph/cliques_to_jtree.m b/sourcecodes/bnt-master/graph/cliques_to_jtree.m new file mode 100644 index 00000000..921ced0b --- /dev/null +++ b/sourcecodes/bnt-master/graph/cliques_to_jtree.m @@ -0,0 +1,57 @@ +function [jtree, root, B, w] = cliques_to_jtree(cliques, ns) +% MK_JTREE Make an optimal junction tree. +% [jtree, root, B, w] = mk_jtree(cliques, ns) +% +% A junction tree is a tree that satisfies the jtree property, which says: +% for each pair of cliques U,V with intersection S, all cliques on the path between U and V +% contain S. (This ensures that local propagation leads to global consistency.) +% +% We can create a junction tree by computing the maximal spanning tree of the junction graph. +% (The junction graph connects all cliques, and the weight of an edge (i,j) is +% |C(i) intersect C(j)|, where C(i) is the i'th clique.) +% +% The best jtree is the maximal spanning tree which minimizes the sum of the costs on each edge, +% where cost(i,j) = w(C(i)) + w(C(j)), and w(C) is the weight of clique C, +% which is the total number of values C can take on. +% +% For details, see +% - Jensen and Jensen, "Optimal Junction Trees", UAI 94. +% +% Input: +% cliques{i} = nodes in clique i +% ns(i) = number of values node i can take on +% Output: +% jtree(i,j) = 1 iff cliques i and j aer connected +% root = the clique that should be used as root +% B(i,j) = 1 iff node j occurs in clique i +% w(i) = weight of clique i + + + +num_cliques = length(cliques); +w = zeros(num_cliques, 1); +B = sparse(num_cliques, 1); +for i=1:num_cliques + B(i, cliques{i}) = 1; + w(i) = prod(ns(cliques{i})); +end + + +% C1(i,j) = length(intersect(cliques{i}, cliques{j})); +% The length of the intersection of two sets is the dot product of their bit vector representation. +C1 = B*B'; +C1 = setdiag(C1, 0); + +% C2(i,j) = w(i) + w(j) +num_cliques = length(w); +W = repmat(w, 1, num_cliques); +C2 = W + W'; +C2 = setdiag(C2, 0); + +jtree = sparse(minimum_spanning_tree(-C1, C2)); % Using -C1 gives *maximum* spanning tree + +% The root is arbitrary, but since the first pass is towards the root, +% we would like this to correspond to going forward in time in a DBN. +root = num_cliques; + + diff --git a/sourcecodes/bnt-master/graph/cliques_to_strong_jtree.m b/sourcecodes/bnt-master/graph/cliques_to_strong_jtree.m new file mode 100644 index 00000000..551894f2 --- /dev/null +++ b/sourcecodes/bnt-master/graph/cliques_to_strong_jtree.m @@ -0,0 +1,114 @@ +function [jtree, root, cliques, B, w] = mk_strong_jtree(cliques, ns, elim_order, MTG) +% MK_SRONG_JTREE Make a strong junction tree. +% [jtree, root, cliques, B, w] = mk_strong_jtree(cliques, ns, elim_order, MTG) +% +% Here is a definition of a strong jtree from Jensen et al. 1994: +% "A junction tree is said to be strong if it has at least one distinguished clique R, +% called a strong root, s.t. for each pair (C1,C2) of adjacent cliques in the tree, +% with C1 closer to R than C2, there exists and ordering of [the nodes below] C2 +% that respects [the partial order] and with the vertices of the separator C1 intersect C2 +% preceeding the vertices [below C2] of C2 \ C1." +% +% For details, see +% - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94. +% +% MTG is the moralized, triangulated graph. +% elim_order is the elimination ordering used to compute MTG. + + +% Warning: this is a very naive implementation of the algorithm in Jensen et al. + +n = length(elim_order); +alpha(elim_order) = n:-1:1; +% alpha(u) = i if we eliminate u at step n-i+1 +% i.e., vertices with higher alpha numbers are eliminated before vertices with lower numbers. +% e.g., from the Jensen et al paper +% node a=1 eliminated at step 6, so alpha(a)=16-6+1=11. +% alpha = [11 1 2 10 9 3 4 7 5 8 13 12 6 16 15 14] + + +% We sort the cliques in order of increasing index. The index of a clique C is defined as follows. +% Let lower = {u | alpha(u) < alpha(v)}, and +% let v in C be the highest-numbered vertex s.t. the vertices in W = lower intersect C +% have a common neighbor u in U, where U = lower \ C. +% If such a v exists, define index(C) = alpha(v), otherwise, index(C) = 1. +% Intuitively, index(C) is the step in the elimination process at which C disappears. + +num_cliques = length(cliques); +index = zeros(1, num_cliques); +for c = 1:num_cliques + C = cliques{c}; + highest_num = -inf; + for vi = 1:length(C) + v = C(vi); + lower = find(alpha < alpha(v)); + W = myintersect(lower, C); + U = mysetdiff(lower, C); + found = 0; + for ui=1:length(U) + u = U(ui); + if mysubset(W, neighbors(MTG, u)) + found = 1; + break; + end + end + if found + if alpha(v) > highest_num + highest_num = alpha(v); + end + end + end + if highest_num == -inf + index(c) = 1; + else + index(c) = highest_num; + end +end + + +% Permute the cliques so that they are ordered according to index +[dummy, clique_order] = sort(index); +cliques = cliques(clique_order); + +w = zeros(num_cliques, 1); +B = sparse(num_cliques, 1); +for i=1:num_cliques + B(i, cliques{i}) = 1; + w(i) = prod(ns(cliques{i})); +end + +% Pearl p113 suggests ordering the cliques by rank of the highest vertex in each clique. +% However, this will only work if we use maximum cardinality search. + + +% Join up the cliques so that they satisfy the Running Intersection Property. +% This states that, for all k > 1, S(k) subseteq C(j) for some j < k, where +% S(k) = C(k) intersect (union_{i=1}^{k-1} C(i)) +jtree = sparse(num_cliques, num_cliques); +for k=2:num_cliques + S = []; + for i=1:k-1 + S = myunion(S, cliques{i}); + end + S = myintersect(S, cliques{k}); + found = 0; + for j=1:k-1 + if mysubset(S, cliques{j}) + found = 1; + break; + end + end + if ~found + disp(['RIP is violated for clique ' num2str(k)]); + end + jtree(k,j)=1; + jtree(j,k)=1; +end + +% Pearl p113 suggests connecting Ci to a predecessor Cj (j < i) sharing +% the highest number of vertices with Ci (i.e., the heaviest i-j edge +% in the jgraph). However, this will only work if we use maximum cardinality search. + +root = 1; + + diff --git a/sourcecodes/bnt-master/graph/connected_graph.m b/sourcecodes/bnt-master/graph/connected_graph.m new file mode 100644 index 00000000..612393ea --- /dev/null +++ b/sourcecodes/bnt-master/graph/connected_graph.m @@ -0,0 +1,9 @@ +function b = connected(adj_mat, directed) +% +% b = connected(adj_mat). +% Returns true iff the graph is connected. + +n = length(adj_mat); +start = 1; +[d, pre] = dfs(adj_mat, start, directed); +b = (length(pre) == n); diff --git a/sourcecodes/bnt-master/graph/dag_to_essential_graph.m b/sourcecodes/bnt-master/graph/dag_to_essential_graph.m new file mode 100644 index 00000000..84dd8ee5 --- /dev/null +++ b/sourcecodes/bnt-master/graph/dag_to_essential_graph.m @@ -0,0 +1,112 @@ + +function [eg] = dag_to_essential_graph(dag) +cpdag = dag_to_cpdag(dag); +eg = dag + dag .* (cpdag + cpdag'); + +return; + + + + +% Coverts a DAG into Essential Graph where edges are coded by 2 and 3, 2 is +% directed edge and 3 is bidirected edge and is at one (the same as the original DAG) of the two +% symetrical places. + +% Is implemented by the algorithm of Max Chickering in D.M.Chickering (1995). +% A transformational characterization of equivalent Bayesian network structures. +% In Proceedings of Eleventh Conference on Uncertainty in Artificial Intelligence, Montreal, QU, +% pages 87-98. Morgan Kaufmann +% http://research.microsoft.com/~dmax/publications/uai95.pdf + +% Implemented by Tomas Kocka, AAU. + +function [eg] = dag_to_essential_graph(dagx) + +%print_dag(dagx); % Just checking input + +order = topological_sort(dagx); % get the topological order of nodes and their number + +% fprintf('the topological order is: %d',order); +% fprintf('\n'); + +[nx,ny] = size(dagx); % gets the number of nodes, note that nx == ny +[I,J] = find(dagx); % finds all nonzero elements in the adjacency matrix, i.e. arcs in the DAG - however we will overwrite it in a special order +% we will sort the arcs from lowest possible y and highest possible x, arcs are x->y +e = 1; +for y = 1:ny + for x = nx:-1:1 + %fprintf('x %d ',order(x)); fprintf('y %d ',order(y)); + if dagx(order(x),order(y)) == 1 + I(e) = order(x); + J(e) = order(y); + e = e + 1; + %fprintf('x order %d',x); + %fprintf('y order %d',y); + %fprintf('\n'); + end + end +end + + +% fprintf('the arcs are: %d',I); +% fprintf('\n'); +% fprintf('the arcs are: %d',J); +% fprintf('\n'); + + +% Now we have to decide which arcs are part of the essential graph and +% which are undirected edges in the essential graph. +% Undecided arc in the DAG are 1, directed in EG are 2 and undirected in EG +% are 3. + + +for e = 1:length(I) + if dagx(I(e),J(e)) == 1 + cont = true; + for w = 1:nx + if dagx(w,I(e)) == 2 + if dagx(w,J(e)) ~= 0 + dagx(w,J(e)) = 2; + else + for ww = 1:nx + if dagx(ww,J(e)) ~= 0 + dagx(ww,J(e)) = 2; + end + end % and now skip the rest and start with another arc from the list + w = nx; + cont = false; + end + end + end + if cont + exists = false; + for z = 1:nx + %fprintf('test %d',dagx(z,J(e))); + if dagx(z,J(e)) ~= 0 & z ~= I(e) & dagx(z,I(e)) == 0 + exists = true; + for ww = 1:nx + if dagx(ww,J(e)) == 1 + dagx(ww,J(e)) = 2; + end + end + end + end + if ~ exists + for ww = 1:nx + if dagx(ww,J(e)) == 1 + dagx(ww,J(e)) = 3; + end + end + end + end + end +end + +%print_dag(dagx); % Just checking output + + + + + + + diff --git a/sourcecodes/bnt-master/graph/dfs.m b/sourcecodes/bnt-master/graph/dfs.m new file mode 100644 index 00000000..cd13c358 --- /dev/null +++ b/sourcecodes/bnt-master/graph/dfs.m @@ -0,0 +1,91 @@ +function [d, pre, post, cycle, f, pred] = dfs(adj_mat, start, directed) +% DFS Perform a depth-first search of the graph starting from 'start'. +% [d, pre, post, cycle, f, pred] = dfs(adj_mat, start, directed) +% +% Input: +% adj_mat(i,j)=1 iff i is connected to j. +% start is the root vertex of the dfs tree; if [], all nodes are searched +% directed = 1 if the graph is directed +% +% Output: +% d(i) is the time at which node i is first discovered. +% pre is a list of the nodes in the order in which they are first encountered (opened). +% post is a list of the nodes in the order in which they are last encountered (closed). +% 'cycle' is true iff a (directed) cycle is found. +% f(i) is the time at which node i is finished. +% pred(i) is the predecessor of i in the dfs tree. +% +% If the graph is a tree, preorder is parents before children, +% and postorder is children before parents. +% For a DAG, topological order = reverse(postorder). +% +% See Cormen, Leiserson and Rivest, "An intro. to algorithms" 1994, p478. + +n = length(adj_mat); + +global white gray black color +white = 0; gray = 1; black = 2; +color = white*ones(1,n); + +global time_stamp +time_stamp = 0; + +global d f +d = zeros(1,n); +f = zeros(1,n); + +global pred +pred = zeros(1,n); + +global cycle +cycle = 0; + +global pre post +pre = []; +post = []; + +if ~isempty(start) + dfs_visit(start, adj_mat, directed); +end +for u=1:n + if color(u)==white + dfs_visit(u, adj_mat, directed); + end +end + + +%%%%%%%%%% + +function dfs_visit(u, adj_mat, directed) + +global white gray black color time_stamp d f pred cycle pre post + +pre = [pre u]; +color(u) = gray; +time_stamp = time_stamp + 1; +d(u) = time_stamp; +if directed + ns = children(adj_mat, u); +else + ns = neighbors(adj_mat, u); + ns = setdiff(ns, pred(u)); % don't go back to visit the guy who called you! +end +for v=ns(:)' + %fprintf('u=%d, v=%d, color(v)=%d\n', u, v, color(v)) + switch color(v) + case white, % not visited v before (tree edge) + pred(v)=u; + dfs_visit(v, adj_mat, directed); + case gray, % back edge - v has been visited, but is still open + cycle = 1; + %fprintf('cycle: back edge from v=%d to u=%d\n', v, u); + case black, % v has been visited, but is closed + % no-op + end +end +color(u) = black; +post = [post u]; +time_stamp = time_stamp + 1; +f(u) = time_stamp; + + diff --git a/sourcecodes/bnt-master/graph/dfs_test.m b/sourcecodes/bnt-master/graph/dfs_test.m new file mode 100644 index 00000000..e3e4ae12 --- /dev/null +++ b/sourcecodes/bnt-master/graph/dfs_test.m @@ -0,0 +1,30 @@ +% Do the example in fig 23.4 p479 of Cormen, Leiserson and Rivest (1994) + +u = 1; v = 2; w = 3; x = 4; y = 5; z = 6; +n = 6; +dag=zeros(n,n); +dag(u,[v x])=1; +dag(v,y)=1; +dag(w,[y z])=1; +dag(x,v)=1; +dag(y,x)=1; +dag(z,z)=1; + +[d, pre, post, cycle, f, pred] = dfs(dag, [], 1); +assert(isequal(d, [1 2 9 4 3 10])) +assert(isequal(f, [8 7 12 5 6 11]) +assert(cycle) + +% Now give it an undirected cyclic graph +G = mk_2D_lattice(2,2,0); +% 1 - 3 +% | | +% 2 - 4 +[d, pre, post, cycle, f, pred] = dfs(G, [], 0); +% d = [1 2 4 3] +assert(cycle) + +% Now break the cycle +G(1,2)=0; G(2,1)=0; +[d, pre, post, cycle, f, pred] = dfs(G, [], 0); +assert(~cycle) diff --git a/sourcecodes/bnt-master/graph/dijkstra.m b/sourcecodes/bnt-master/graph/dijkstra.m new file mode 100644 index 00000000..2e444d4c --- /dev/null +++ b/sourcecodes/bnt-master/graph/dijkstra.m @@ -0,0 +1,112 @@ +function [D,P] = dijk(A,s,t) +%DIJK Shortest paths from nodes 's' to nodes 't' using Dijkstra algorithm. +% [D,p] = dijk(A,s,t) +% A = n x n node-node weighted adjacency matrix of arc lengths +% (Note: A(i,j) = 0 => Arc (i,j) does not exist; +% A(i,j) = NaN => Arc (i,j) exists with 0 weight) +% s = FROM node indices +% = [] (default), paths from all nodes +% t = TO node indices +% = [] (default), paths to all nodes +% D = |s| x |t| matrix of shortest path distances from 's' to 't' +% = [D(i,j)], where D(i,j) = distance from node 'i' to node 'j' +% P = |s| x n matrix of predecessor indices, where P(i,j) is the +% index of the predecessor to node 'j' on the path from 's(i)' to 'j' +% (use PRED2PATH to convert P to paths) +% = path from 's' to 't', if |s| = |t| = 1 +% +% (If A is a triangular matrix, then computationally intensive node +% selection step not needed since graph is acyclic (triangularity is a +% sufficient, but not a necessary, condition for a graph to be acyclic) +% and A can have non-negative elements) +% +% (If |s| >> |t|, then DIJK is faster if DIJK(A',t,s) used, where D is now +% transposed and P now represents successor indices) +% +% (Based on Fig. 4.6 in Ahuja, Magnanti, and Orlin, Network Flows, +% Prentice-Hall, 1993, p. 109.) + +% Copyright (c) 1998-2001 by Michael G. Kay +% Matlog Version 5 22-Aug-2001 + +% Input Error Checking ****************************************************** +error(nargchk(1,3,nargin)); + +[n,cA] = size(A); + +if nargin < 2 | isempty(s), s = (1:n)'; else s = s(:); end +if nargin < 3 | isempty(t), t = (1:n)'; else t = t(:); end + +if ~any(any(tril(A) ~= 0)) % A is upper triangular + isAcyclic = 1; +elseif ~any(any(triu(A) ~= 0)) % A is lower triangular + isAcyclic = 2; +else % Graph may not be acyclic + isAcyclic = 0; +end + +if n ~= cA + error('A must be a square matrix'); +elseif ~isAcyclic & any(any(A < 0)) + error('A must be non-negative'); +elseif any(s < 1 | s > n) + error(['''s'' must be an integer between 1 and ',num2str(n)]); +elseif any(t < 1 | t > n) + error(['''t'' must be an integer between 1 and ',num2str(n)]); +end +% End (Input Error Checking) ************************************************ + +A = A'; % Use transpose to speed-up FIND for sparse A + +D = zeros(length(s),length(t)); +if nargout > 1, P = zeros(length(s),n); end + +for i = 1:length(s) + j = s(i); + + Di = Inf*ones(n,1); Di(j) = 0; + + isLab = logical(zeros(length(t),1)); + if isAcyclic == 1 + nLab = j - 1; + elseif isAcyclic == 2 + nLab = n - j; + else + nLab = 0; + UnLab = 1:n; + isUnLab = logical(ones(n,1)); + end + + while nLab < n & ~all(isLab) + if isAcyclic + Dj = Di(j); + else % Node selection + [Dj,jj] = min(Di(isUnLab)); + j = UnLab(jj); + UnLab(jj) = []; + isUnLab(j) = 0; + end + + nLab = nLab + 1; + if length(t) < n, isLab = isLab | (j == t); end + + [jA,kA,Aj] = find(A(:,j)); + Aj(isnan(Aj)) = 0; + + if isempty(Aj), Dk = Inf; else Dk = Dj + Aj; end + + if nargout > 1, P(i,jA(Dk < Di(jA))) = j; end + Di(jA) = min(Di(jA),Dk); + + if isAcyclic == 1 % Increment node index for upper triangular A + j = j + 1; + elseif isAcyclic == 2 % Decrement node index for lower triangular A + j = j - 1; + end + end + D(i,:) = Di(t)'; +end + +if nargout > 1 & length(s) == 1 & length(t) == 1 + P = pred2path(P,s,t); +end diff --git a/sourcecodes/bnt-master/graph/family.m b/sourcecodes/bnt-master/graph/family.m new file mode 100644 index 00000000..0a6494d6 --- /dev/null +++ b/sourcecodes/bnt-master/graph/family.m @@ -0,0 +1,17 @@ +function f = family(A,i,t) +% FAMILY Return the indices of parents and self in sorted order +% f = family(dag,i,t) +% +% t is an optional argument: if present, dag is assumed to be a 2-slice DBN + +if nargin < 3 + f = [parents(A,i) i]; +else + if t == 1 + f = [parents(A,i) i]; + else + ss = length(A)/2; + j = i+ss; + f = [parents(A,j) j] + (t-2)*ss; + end +end diff --git a/sourcecodes/bnt-master/graph/findroot.m b/sourcecodes/bnt-master/graph/findroot.m new file mode 100644 index 00000000..d242a3a3 --- /dev/null +++ b/sourcecodes/bnt-master/graph/findroot.m @@ -0,0 +1,24 @@ +function root = findroot(bnet, cliques) + +%% findroot is to find the strong root in a clique tree assume it has one +%% in the tree. For a clique tree constructed from a strongly triangulated +%% graph, an interface clique that contains all discrete parents +%% and at least one continuous node from a connected continuous component +%% is for sure to be available as a guaranteed strong root. +%% -By Wei Sun, George Mason University, 4/17/2010. + +%% We choose the interface clique that contains the max number +%% of interface nodes to be the strong root. +n0 = 0 ; +for i=1:length(cliques) + % check hybrid cliques + hc = intersect(cliques{i}, bnet.cnodes) ; + hd = intersect(cliques{i}, bnet.dnodes) ; + if ~isempty(hd) & ~isempty(hc) + nd = length(hd) ; + if nd > n0 + root = i ; + n0 = nd ; + end + end +end diff --git a/sourcecodes/bnt-master/graph/graph_separated.m b/sourcecodes/bnt-master/graph/graph_separated.m new file mode 100644 index 00000000..32a11841 --- /dev/null +++ b/sourcecodes/bnt-master/graph/graph_separated.m @@ -0,0 +1,9 @@ +function sep = graph_separated(G, X, Y, S) + +G2 = G; +G2(S,:) = 0; +G2(:,S) = 0; +conn = reachability_graph(G2); +conn2 = conn(X,Y); +sep = all(conn2(:)==0); + diff --git a/sourcecodes/bnt-master/graph/graph_to_jtree.m b/sourcecodes/bnt-master/graph/graph_to_jtree.m new file mode 100644 index 00000000..952bbca7 --- /dev/null +++ b/sourcecodes/bnt-master/graph/graph_to_jtree.m @@ -0,0 +1,59 @@ +function [jtree, root, cliques, B, w, elim_order] = graph_to_jtree(MG, ns, partial_order, stages, clusters) +% GRAPH_TO_JTREE Triangulate a graph and make a junction tree from its cliques. +% [jtree, root, cliques, B, w, elim_order] = ... +% graph_to_jtree(graph, node_sizes, partial_order, stages, clusters) +% +% INPUT: +% graph(i,j) = 1 iff there is an edge between i,j +% node_weights(i) = num discrete values node i can take on [1 if observed] +% partial_order = {} if no constraints on elimination ordering +% stages{i} = nodes that must be eliminated at i'th stage (if porder is empty) +% clusters{i} = list of nodes that must get connected together in the moral graph +% +% OUTPUT: +% jtree(i,j) = 1 iff there is an arc between clique i and clique j +% root = the root clique +% cliques{i} = the nodes in clique i +% B(i,j) = 1 iff node j occurs in clique i +% w(i) = weight of clique i + +N = length(MG); + +if nargin >= 5 + % Add extra arcs between nodes in each cluster to ensure they occur in the same clique + for i=1:length(clusters) + c = clusters{i}; + MG(c,c) = 1; + end +end +MG = setdiag(MG, 0); + +% Find an optimal elimination ordering (NP-hard problem!) +if nargin < 4 + stages = {1:N}; +end +if nargin < 3 + partial_order = {}; +end +if isempty(partial_order) + strong = 0; + elim_order = best_first_elim_order(MG, ns, stages); +else + strong = 1; + elim_order = strong_elim_order(MG, ns, partial_order); +end + +[MTG, cliques, fill_in_edges] = triangulate(MG, elim_order); + +% Connect the cliques up into a jtree, +[jtree, root, B, w] = cliques_to_jtree(cliques, ns); + +if 0 + disp('testing dag to jtree'); + % Find the cliques containing each node, and check they form a connected subtree + clqs_con_node = cell(1,N); + for i=1:N + clqs_con_node{i} = find(B(:,i))'; + end + check_jtree_property(clqs_con_node, jtree); +end diff --git a/sourcecodes/bnt-master/graph/min_subtree_con_nodes.m b/sourcecodes/bnt-master/graph/min_subtree_con_nodes.m new file mode 100644 index 00000000..0bc3c2d3 --- /dev/null +++ b/sourcecodes/bnt-master/graph/min_subtree_con_nodes.m @@ -0,0 +1,60 @@ +function [subtree, nroot_node] = min_subtree_con_nodes(jtree, root, nodes) +%min_subtree_con_nodes get the minimum subtree of tree which contains the nodes + +if isempty(jtree) | isempty(nodes) + subtree = []; + nroot_node = []; + return; +end + +rnodes = min_subtree_nodes(jtree, nodes); +nea_node = nearest_node(jtree, root, nodes); +node_num = length(jtree); +subtree = zeros(node_num); +subtree(rnodes, rnodes) = jtree(rnodes, rnodes); +nroot_node = nea_node; + + +function rnodes = min_subtree_nodes(tree, nodes) +rnodes = []; +if isempty(tree) | isempty(nodes) + return +end + +rnodes = nodes(1); +newnodes = neighbors(tree, nodes(1)); +while ~mysubset(nodes, rnodes) + swapnodes = newnodes; + newnodes = []; + added = 0; + for i=1:length(swapnodes) + inode = swapnodes(i); + tnodes = myunion(inode, rnodes); + if mysubset(nodes, tnodes) + added = 1; + break; + end + nns = neighbors(tree, inode); + add_nodes = mysetdiff(nns, tnodes); + newnodes = myunion(newnodes, add_nodes); + end + if added + rnodes = tnodes; + else + rnodes = myunion(rnodes, newnodes); + end +end + +function nea_node = nearest_node(tree, inode, nodes) +if myismember(inode, nodes) + nea_node = inode; + return; +end +cs = children(tree, inode); +for i = 1:length(cs) + n = cs(i); + nea_node = nearest_node(tree, n, nodes); +end + + + diff --git a/sourcecodes/bnt-master/graph/minimum_spanning_tree.m b/sourcecodes/bnt-master/graph/minimum_spanning_tree.m new file mode 100644 index 00000000..ca208985 --- /dev/null +++ b/sourcecodes/bnt-master/graph/minimum_spanning_tree.m @@ -0,0 +1,51 @@ +function A = minimum_spanning_tree(C1, C2) +% +% Find the minimum spanning tree using Prim's algorithm. +% C1(i,j) is the primary cost of connecting i to j. +% C2(i,j) is the (optional) secondary cost of connecting i to j, used to break ties. +% We assume that absent edges have 0 cost. +% To find the maximum spanning tree, used -1*C. +% See Aho, Hopcroft & Ullman 1983, "Data structures and algorithms", p 237. + +% Prim's is O(V^2). Kruskal's algorithm is O(E log E) and hence is more efficient +% for sparse graphs, but is implemented in terms of a priority queue. + +% We partition the nodes into those in U and those not in U. +% closest(i) is the vertex in U that is closest to i in V-U. +% lowcost(i) is the cost of the edge (i, closest(i)), or infinity is i has been used. +% In Aho, they say C(i,j) should be "some appropriate large value" if the edge is missing. +% We set it to infinity. +% However, since lowcost is initialized from C, we must distinguish absent edges from used nodes. + +n = length(C1); +if nargin==1, C2 = zeros(n); end +A = zeros(n); + +closest = ones(1,n); +used = zeros(1,n); % contains the members of U +used(1) = 1; % start with node 1 +C1(find(C1==0))=inf; +C2(find(C2==0))=inf; +lowcost1 = C1(1,:); +lowcost2 = C2(1,:); + +for i=2:n + ks = find(lowcost1==min(lowcost1)); + k = ks(argmin(lowcost2(ks))); + A(k, closest(k)) = 1; + A(closest(k), k) = 1; + lowcost1(k) = inf; + lowcost2(k) = inf; + used(k) = 1; + NU = find(used==0); + for ji=1:length(NU) + for j=NU(ji) + if C1(k,j) < lowcost1(j) + lowcost1(j) = C1(k,j); + lowcost2(j) = C2(k,j); + closest(j) = k; + end + end + end +end + diff --git a/sourcecodes/bnt-master/graph/minspan.m b/sourcecodes/bnt-master/graph/minspan.m new file mode 100644 index 00000000..23d396fb --- /dev/null +++ b/sourcecodes/bnt-master/graph/minspan.m @@ -0,0 +1,56 @@ +function [t,nk] = minspan(IJC) +%MINSPAN Minimum weight spanning tree using Kruskal algorithm. +%[t,nk] = minspan(IJC) +% IJC = n x 3 matrix arc list [i j c] of arc heads, tails, and costs +% t = n-element logical vector, where +% t(i) = 1, if IJC(i,:) arc in spanning tree +% t(i) = k, if IJC(i,:) arc in component k of forest +% nk = number of components + +% Copyright (c) 1998-2001 by Michael G. Kay +% Matlog Version 5 22-Aug-2001 + +% Input Error Checking ****************************************************** +[n,cIJC] = size(IJC); +if cIJC ~= 3 + error('''IJC'' must be a 3-column matrix.') +elseif n < 1 + error('''IJC'' must have at least one row.') +elseif any(IJC(:,1) < 1) | any(any(~isint(IJC(:,[1 2])))) + error('Invalid arc index in IJC.') +end +% End (Input Error Checking) ************************************************ + +i = IJC(:,1); j = abs(IJC(:,2)); +m = max(max([i j])); + +sidxIJ = argsort(IJC(:,3)); +i = i(sidxIJ); j = j(sidxIJ); + +t = logical(zeros(n,1)); +k = 1; % Current arc +nt = 0; % Number of arcs in spanning tree +v = (1:m)'; % Arc labels + +while nt < m - 1 & k <= n + if (v(i(k)) ~= v(j(k))) + v(v==v(j(k))) = v(i(k)); + t(k) = 1; + nt = nt + 1; + end + k = k + 1; +end + +idxIJ = invperm(sidxIJ); +t = t(idxIJ); i = i(idxIJ); j = j(idxIJ); + +c = unique(v(unique([i; j]))); % Unique labels of arc vertices +nk = length(c); +if ~any(t), nk = 0; end % Self-loop not a component + +if nk > 1 + for k = 1:nk + t(t~=0 & v(i)==c(k)) = k; % Relabel to consecutive component numbers + end +end + diff --git a/sourcecodes/bnt-master/graph/mk_2D_lattice.m b/sourcecodes/bnt-master/graph/mk_2D_lattice.m new file mode 100644 index 00000000..a9873260 --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_2D_lattice.m @@ -0,0 +1,59 @@ +function G = mk_2D_lattice(nrows, ncols, con) +% MK_2D_LATTICE Return adjacency matrix for nearest neighbor connected 2D lattice +% G = mk_2D_lattice(nrows, ncols, con) +% G(k1, k2) = 1 iff k1=(i1,j1) is a neighbor of k2=(i2,j2) +% (Two pixels are neighbors if their Euclidean distance is less than r.) +% Default connectivity = 4. +% +% WE ASSUME NO WRAP AROUND. +% +% This is the neighborhood as a function of con: +% +% con=4,r=1 con=8,r=sqrt(2) con=12,r=2 con=24,r=sqrt(8) +% nn 2nd order 4th order +% x x x x x x +% x x x x x x x x x x x x +% x o x x o x x x o x x x x o x x +% x x x x x x x x x x x x +% x x x x x x +% +% Examples: +% Consider a 3x4 grid +% 1 4 7 10 +% 2 5 8 11 +% 3 6 9 12 +% +% 4-connected: +% G=mk_2D_lattice(3,4,4); +% find(G(1,:)) = [2 4] +% find(G(5,:)) = [2 4 6 8] +% +% 8-connected: +% G=mk_2D_lattice(3,4,8); +% find(G(1,:)) = [2 4 5] +% find(G(5,:)) = [1 2 3 4 6 7 8 9] + +% meshgrid trick due to Temu Gautama (temu@neuro.kuleuven.ac.be) + +if nargin < 3, con = 4; end + +switch con, + case 4, r = 1; + case 8, r = sqrt(2); + case 12, r = 2; + case 24, r = sqrt(8); + otherwise, error(['unrecognized connectivity ' num2str(con)]) +end + + +npixels = nrows*ncols; + +[x y]=meshgrid(1:ncols, 1:nrows); +M = [x(:) y(:)]; +M1 = repmat(reshape(M',[1 2 npixels]),[npixels 1 1]); +M2 = repmat(M,[1 1 npixels]); +%D = squeeze(sum(abs(M1-M2),2)); % Manhattan distance +M3 = M1-M2; +D = sqrt(squeeze(M3(:,1,:)) .^2 + squeeze(M3(:,2,:)) .^2); % Euclidean distance +G = reshape(D <= r,npixels,npixels); +G = setdiag(G, 0); diff --git a/sourcecodes/bnt-master/graph/mk_2D_lattice_slow.m b/sourcecodes/bnt-master/graph/mk_2D_lattice_slow.m new file mode 100644 index 00000000..0c646d80 --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_2D_lattice_slow.m @@ -0,0 +1,135 @@ +function G = mk_2D_lattice_slow(nrows, ncols, wrap_around) +% MK_2D_LATTICE Return adjacency matrix for 4-nearest neighbor connected 2D lattice +% G = mk_2D_lattice(nrows, ncols, wrap_around) +% G(k1, k2) = 1 iff k1=(i1,j1) is connected to k2=(i2,j2) +% +% If wrap_around = 1, we use toroidal boundary conditions (default = 0) +% +% Nodes are assumed numbered as in the following 3x3 lattice +% 1 4 7 +% 2 5 8 +% 3 6 9 +% +% e.g., G = mk_2D_lattice(3, 3, 0) returns +% 0 1 0 1 0 0 0 0 0 +% 1 0 1 0 1 0 0 0 0 +% 0 1 0 0 0 1 0 0 0 +% 1 0 0 0 1 0 1 0 0 +% 0 1 0 1 0 1 0 1 0 +% 0 0 1 0 1 0 0 0 1 +% 0 0 0 1 0 0 0 1 0 +% 0 0 0 0 1 0 1 0 1 +% 0 0 0 0 0 1 0 1 0 +% so find(G(5,:)) = [2 4 6 8] +% but find(G(1,:)) = [2 4] +% +% Using wrap around, G = mk_2D_lattice(3, 3, 1), we get +% 0 1 1 1 0 0 1 0 0 +% 1 0 1 0 1 0 0 1 0 +% 1 1 0 0 0 1 0 0 1 +% 1 0 0 0 1 1 1 0 0 +% 0 1 0 1 0 1 0 1 0 +% 0 0 1 1 1 0 0 0 1 +% 1 0 0 1 0 0 0 1 1 +% 0 1 0 0 1 0 1 0 1 +% 0 0 1 0 0 1 1 1 0 +% so find(G(5,:)) = [2 4 6 8] +% and find(G(1,:)) = [2 3 4 7] + +if nargin < 3, wrap_around = 0; end + +% M contains the number of each cell e.g. +% 1 4 7 +% 2 5 8 +% 3 6 9 +% North neighbors (assuming wrap around) are +% 3 6 9 +% 1 4 7 +% 2 5 8 +% Without wrap around, they are +% 1 4 7 +% 1 4 7 +% 2 5 8 +% The first row is arbitrary, since pixels at the top have no north neighbor. + +if nrows==1 + G = zeros(1, ncols); + for i=1:ncols-1 + G(i,i+1) = 1; + G(i+1,i) = 1; + end + if wrap_around + G(1,ncols) = 1; + G(ncols,1) = 1; + end + return; +end + + +npixels = nrows*ncols; + +N = 1; E = 2; S = 3; W = 4; +if wrap_around + rows{N} = [nrows 1:nrows-1]; cols{N} = 1:ncols; + rows{E} = 1:nrows; cols{E} = [2:ncols 1]; + rows{S} = [2:nrows 1]; cols{S} = 1:ncols; + rows{W} = 1:nrows; cols{W} = [ncols 1:ncols-1]; +else + rows{N} = [1 1:nrows-1]; cols{N} = 1:ncols; + rows{E} = 1:nrows; cols{E} = [1 1:ncols-1]; + rows{S} = [2:nrows nrows]; cols{S} = 1:ncols; + rows{W} = 1:nrows; cols{W} = [2:ncols ncols]; +end + +M = reshape(1:npixels, [nrows ncols]); +nbrs = cell(1, 4); +for i=1:4 + nbrs{i} = M(rows{i}, cols{i}); +end + + +G = zeros(npixels, npixels); +if wrap_around + for i=1:4 + if 0 + % naive + for p=1:npixels + G(p, nbrs{i}(p)) = 1; + end + else + % vectorized + ndx2 = sub2ind([npixels npixels], 1:npixels, nbrs{i}(:)'); + G(ndx2) = 1; + end + end +else + i = N; + mask = ones(nrows, ncols); + mask(1,:) = 0; % pixels in row 1 have no nbr to the north + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; + + i = E; + mask = ones(nrows, ncols); + mask(:,ncols) = 0; + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; + + i = S; + mask = ones(nrows, ncols); + mask(nrows,:)=0; + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; + + i = W; + mask = ones(nrows, ncols); + mask(:,1)=0; + ndx = find(mask); + ndx2 = sub2ind([npixels npixels], ndx, nbrs{i}(ndx)); + G(ndx2) = 1; +end + +G = setdiag(G, 0); diff --git a/sourcecodes/bnt-master/graph/mk_adj_mat.m b/sourcecodes/bnt-master/graph/mk_adj_mat.m new file mode 100644 index 00000000..09f73d4a --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_adj_mat.m @@ -0,0 +1,36 @@ +function [A, names] = mk_adj_mat(connections, names, topological) +% MK_ADJ_MAT Make a directed adjacency matrix from a list of connections between named nodes. +% +% A = mk_adj_mat(connections, name) +% This is best explaine by an example: +% names = {'WetGrass', 'Sprinkler', 'Cloudy', 'Rain'}; +% connections = {'Cloudy', 'Sprinkler'; 'Cloudy', 'Rain'; 'Sprinkler', 'WetGrass'; 'Rain', 'WetGrass'}; +% adds the arcs C -> S, C -> R, S -> W, R -> W. Node 1 is W, 2 is S, 3 is C, 4 is R. +% +% [A, names] = mk_adj_mat(connections, name, 1) +% The last argument of 1 indicates that we should topologically sort the nodes (parents before children). +% In the example, the numbering becomes: node 1 is C, 2 is R, 3 is S, 4 is W +% and the return value of names gets permuted to {'Cloudy', 'Rain', 'Sprinkler', 'WetGrass'}. +% Note that topological sorting the graph is only possible if it has no directed cycles. + +if nargin < 3, topological = 0; end + +n=length(names); +A=zeros(n); +[nr nc] = size(connections); +for r=1:nr + from = strmatch(connections{r,1}, names, 'exact'); + assert(~isempty(from)); + to = strmatch(connections{r,2}, names, 'exact'); + assert(~isempty(to)); + %fprintf(1, 'from %s %d to %s %d\n', connections{r,1}, from, connections{r,2}, to); + A(from,to) = 1; +end + +if topological + order = topological_sort(A); + A = A(order, order); + names = names(order); +end + + diff --git a/sourcecodes/bnt-master/graph/mk_adjmat_chain.m b/sourcecodes/bnt-master/graph/mk_adjmat_chain.m new file mode 100644 index 00000000..d206b272 --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_adjmat_chain.m @@ -0,0 +1,8 @@ +function G = mk_adjmat_chain(T) +% MK_DAG_CHAIN Make adjacency matrix for bi-directional Markov chain of T nodes +% function G = mk_dag_chain(T) +% +% G(t,t+1) = 1 for all t<T + +G = diag(ones(1,T-1),1) + diag(ones(1,T-1),-1); + diff --git a/sourcecodes/bnt-master/graph/mk_all_dags.m b/sourcecodes/bnt-master/graph/mk_all_dags.m new file mode 100644 index 00000000..7335fd0b --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_all_dags.m @@ -0,0 +1,52 @@ +function Gs = mk_all_dags(N, order) +% MK_ALL_DAGS generate all DAGs on N variables +% G = mk_all_dags(N) +% +% G = mk_all_dags(N, order) only generates DAGs in which node i has parents from +% nodes in order(1:i-1). Default: order=[] (no constraints). +% +% G{i} is the i'th dag +% +% Note: the number of DAGs is super-exponential in N, so don't call this with N > 4. + +if nargin < 2, order = []; end + +use_file = 0; + +global BNT_HOME +fname = sprintf('%s/DAGS%d.mat', BNT_HOME, N); +if use_file & exist(fname, 'file') + S = load(fname, '-mat'); + fprintf('loading %s\n', fname); + Gs = S.Gs; + return; +end + +m = 2^(N*N); +ind = ind2subv(2*ones(1,N^2), 1:m); +Gs = {}; +j = 1; +directed = 1; +for i=1:m + dag = reshape(ind(i,:)-1, N, N); + if acyclic(dag, directed) + out_of_order = 0; + if ~isempty(order) + for k=1:N-1 + if any(dag(order(k+1:end), k)) + out_of_order = 1; + break; + end + end + end + if ~out_of_order + Gs{j} = dag; + j = j + 1; + end + end +end + +if use_file + disp(['mk_all_dags: saving to ' fname '!']); + save(fname, 'Gs'); +end diff --git a/sourcecodes/bnt-master/graph/mk_nbrs_of_dag.m b/sourcecodes/bnt-master/graph/mk_nbrs_of_dag.m new file mode 100644 index 00000000..04742774 --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_nbrs_of_dag.m @@ -0,0 +1,64 @@ +function [Gs, op, nodes] = mk_nbrs_of_dag(G0) +% MK_NBRS_OF_DAG Make all DAGs that differ from G0 by a single edge deletion, addition or reversal +% [Gs, op, nodes] = mk_nbrs_of_dag(G0) +% +% Gs{i} is the i'th neighbor. +% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor. +% nodes(i,1:2) are the head and tail of the operated-on arc. + +Gs = {}; +op = {}; +nodes = []; + +[I,J] = find(G0); +nnbrs = 1; +% all single edge deletions +for e=1:length(I) + i = I(e); j = J(e); + G = G0; + G(i,j) = 0; + Gs{nnbrs} = G; + op{nnbrs} = 'del'; + nodes(nnbrs, :) = [i j]; + nnbrs = nnbrs + 1; +end + +% all single edge reversals +for e=1:length(I) + i = I(e); j = J(e); + G = G0; + G(i,j) = 0; + G(j,i) = 1; + if acyclic(G) + Gs{nnbrs} = G; + op{nnbrs} = 'rev'; + nodes(nnbrs, :) = [i j]; + nnbrs = nnbrs + 1; + end +end + +[I,J] = find(~G0); +% all single edge additions +for e=1:length(I) + i = I(e); j = J(e); + if i ~= j % don't add self arcs + G = G0; + G(i,j) = 1; + if G(j,i)==0 % don't add i->j if j->i exists already + if acyclic(G) + Gs{nnbrs} = G; + op{nnbrs} = 'add'; + nodes(nnbrs, :) = [i j]; + nnbrs = nnbrs + 1; + end + end + end +end + + + + + + + + diff --git a/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph.m b/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph.m new file mode 100644 index 00000000..f71c2b54 --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph.m @@ -0,0 +1,132 @@ +function [Gs, op, nodes, A] = my_mk_nbs_of_digraph(G0,A) +% MY_MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal, subject to acyclicity +% [Gs, op, nodes, A] = my_mk_nbrs_of_digraph(G0,<A>) +% +% G0 is an adj matrix s.t. G0(i,j)=1 iff i->j in graph +% A is the ancestor matrix for G0 (opt, creates if necessary) +% +% Gs(:,:,i) is the i'th neighbor +% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor. +% nodes(i,1:2) are the head and tail of the operated-on arc. +% Modified from mk_nbrs_of_digraph by Sonia Leach +% +% Modified by Sonia Leach Feb 02 + +if nargin ==1, A = reachability_graph(G0');, end + +n = length(G0); +[I,J] = find(G0); % I(k), J(k) is the k'th edge +E = length(I); % num edges present in G0 + + +% SINGLE EDGE DELETIONS +% all deletions are valid wrt acyclity + +Grep = repmat(G0(:), 1, E); % each column is a copy of G0 +% edge_ndx(k) is the scalar location of the k'th edge +edge_ndx = find(G0); + +% edge_ndx = subv2ind([n n], [I J]); % equivalent +% We set (ndx(k), k) to 0 for k=1:E in Grep +ndx = subv2ind(size(Grep), [edge_ndx(:) (1:E)']); +G1 = Grep; +G1(ndx) = 0; +Gdel = reshape(G1, [n n E]); + + +% SINGLE EDGE REVERSALS + +% SML: previously Kevin had that legal structure was if +% A(P,i)=1 for any P = { p | p in parents(j), p~=i} +% specifically he said +% "if any(A(ps,i)) then there is a path i -> parent of j -> j +% so reversing i->j would create a cycle" +% Thus put in another way: +% for each i,j if sum(G0(:,j)' * A(:,i)) > 0, reversing i->j +% is not legal. +% +% Ex. Suppose we want to check if 2->4 can be reversed in the +% following graph: +% G0 = A = +% 0 0 1 0 0 0 0 0 +% 0 0 1 1 0 0 0 0 +% 0 0 0 1 1 1 0 0 +% 0 0 0 0 1 1 1 0 +% +% Then parents(4) = G0(:,4) = [0 1 1 0]' +% and A(:,2) = [0 0 1 1]. Thus G0(:,4)'*A(:,2) = 1 b/c 3 is +% an ancestor of 4 and a child of 2. Note that this works b/c +% matrix multiplication has the effect of ANDing the two vectors +% and summing up the result (equiv. to the any(A(ps,i)) in kevin's code) +% +% So, we vectorize and check for all i,j pairs by looking for +% 1's in L = (G0'*A)' which has L(i,j)=1 if rev(i,j) not legal +% Note that this will give 1's where there are none in the G0 +% so we do a L=max(0, G0-L) to cancel out only the existing edges that +% aren't legal (subtracting where both are 1 and setting where +% G0=0 and A=1 back to 0). + +L = max(0, G0-(G0'*A)'); +[IL, JL] = find(L); % I(k), J(k) is the k'th legal edge to rev. +EL = length(IL); + + +% SML: First we have to DELETE THE EDGES WE ARE REVERSING +% We can't use G1 w/ reversed edges already deleted (as +% Kevin did) b/c the space of possible deletions are different +% now (some reverses aren't legal) + +Grep = repmat(G0(:), 1, EL); % each column is a copy of G0 +% edge_ndx(k) is the scalar location of the k'th edge +edge_ndx = subv2ind([n n], [IL JL]); +% We set (ndx(k), k) to 0 for k=1:E in Grep +ndx = subv2ind(size(Grep), [edge_ndx(:) (1:EL)']); +G1 = Grep; +G1(ndx) = 0; + +% SML: Now we add in our REVERSED EDGES +% rev_edge_ndx(k) is the scalar location of the k'th legal reversed edge +rev_edge_ndx = subv2ind([n n], [JL IL]); + +% We set (rev_edge_ndx(k), k) to 1 for k=1:EL in G1 +% We have already deleted i->j in the previous step +ndx = subv2ind(size(Grep), [rev_edge_ndx(:) (1:EL)']); +G1(ndx) = 1; +Grev = reshape(G1, [n n EL]); + +% SINGLE EDGE ADDITIONS + +% SML: previously Kevin had that any addition was legal if A(i,j)=0 +% however, you can not add i->j if j is a descendent of i. +% Thus, we create all possible additions in Gbar and then +% subtract the descendants of each edge as possible parents +% This means the potential parents of i (i.e. Gbar(:,i)) +% can not also be descendants if i i.e. (A(:,i)) which is accomplished +% by subtracting (Gbar-A == 1 iff Gbar=1 & A=0) + +Gbar = ~G0; % Gbar(i,j)=1 iff there is no i->j edge in G0 +Gbar = setdiag(Gbar, 0); % turn off self loops + +GbarL = Gbar-A; +[IbarL, JbarL] = find(GbarL); % I(k), J(k) is the k'th legal edge to add +EbarL = length(IbarL); + +bar_edge_ndx = find(GbarL); + +Grep = repmat(G0(:), 1, EbarL); % each column is a copy of G0 +ndx = subv2ind(size(Grep), [bar_edge_ndx(:) (1:EbarL)']); +Grep(ndx) = 1; +Gadd = reshape(Grep, [n n EbarL]); + + +Gs = cat(3, Gdel, Grev, Gadd); + +nodes = [I J; + IL JL; + IbarL JbarL]; + +op = cell(1, E+EL+EbarL); +op(1:E) = {'del'}; +op(E+(1:EL)) = {'rev'}; +op((E+EL+1):end) = {'add'}; + diff --git a/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph_broken.m b/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph_broken.m new file mode 100644 index 00000000..8d43afcb --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph_broken.m @@ -0,0 +1,111 @@ +function [Gs, op, nodes] = mk_nbrs_of_digraph(G0) +% MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal +% [Gs, op, nodes] = mk_nbrs_of_digraph(G0) +% +% Gs(:,:,i) is the i'th neighbor +% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor. +% nodes(i,1:2) are the head and tail of the operated-on arc. + +debug = 0; % the vectorized version is about 3 to 10 times faster + +n = length(G0); +[I,J] = find(G0); % I(k), J(k) is the k'th edge +E = length(I); % num edges present in G0 + +% SINGLE EDGE DELETIONS + +Grep = repmat(G0(:), 1, E); % each column is a copy of G0 +% edge_ndx(k) is the scalar location of the k'th edge +edge_ndx = find(G0); +% edge_ndx = subv2ind([n n], [I J]); % equivalent +% We set (ndx(k), k) to 0 for k=1:E in Grep +ndx = subv2ind(size(Grep), [edge_ndx(:) (1:E)']); +G1 = Grep; +G1(ndx) = 0; +Gdel = reshape(G1, [n n E]); + + +% if debug +% % Non-vectorized version +% ctr = 1; +% for e=1:E +% i = I(e); j = J(e); +% Gdel2(:,:,ctr) = G0; +% Gdel2(i,j,ctr) = 0; +% ctr = ctr + 1; +% end +% assert(isequal(Gdel, Gdel2)); +% end + + +% SINGLE EDGE REVERSALS + +% rev_edge_ndx(k) is the scalar location of the k'th reversed edge +%rev_edge_ndx = find(G0'); % different order to edge_ndx, which is bad +rev_edge_ndx = subv2ind([n n], [J I]); +% We set (rev_edge_ndx(k), k) to 1 for k=1:E in G1 +% We have already deleted i->j in the previous step +ndx = subv2ind(size(Grep), [rev_edge_ndx(:) (1:E)']); +G1(ndx) = 1; +Grev = reshape(G1, [n n E]); + +% if debug +% % Non-vectorized version +% ctr = 1; +% for e=1:E +% i = I(e); j = J(e); +% Grev2(:,:,ctr) = G0; +% Grev2(i,j,ctr) = 0; +% Grev2(j,i,ctr) = 1; +% ctr = ctr + 1; +% end +% assert(isequal(Grev, Grev2)); +% end + + +% SINGLE EDGE ADDITIONS + +Gbar = ~G0; % Gbar(i,j)=1 iff there is no i->j edge in G0 +Gbar = setdiag(Gbar, 0); % turn off self loops +[Ibar,Jbar] = find(Gbar); + +bar_edge_ndx = find(Gbar); +Ebar = length(Ibar); % num edges present in Gbar +Grep = repmat(G0(:), 1, Ebar); % each column is a copy of G0 +ndx = subv2ind(size(Grep), [bar_edge_ndx(:) (1:Ebar)']); +Grep(ndx) = 1; +Gadd = reshape(Grep, [n n Ebar]); + +% if debug +% % Non-vectorized version +% ctr = 1; +% for e=1:length(Ibar) +% i = Ibar(e); j = Jbar(e); +% Gadd2(:,:,ctr) = G0; +% Gadd2(i,j,ctr) = 1; +% ctr = ctr + 1; +% end +% assert(isequal(Gadd, Gadd2)); +% end + + +Gs = cat(3, Gdel, Grev, Gadd); + +nodes = [I J; + I J; + Ibar Jbar]; + +op = cell(1, E+E+Ebar); +op(1:E) = {'del'}; +op(E+1:2*E) = {'rev'}; +op(2*E+1:end) = {'add'}; + + +% numeric output: +% op(i) = 1, 2, or 3, if the i'th neighbor was created by adding, deleting or reversing an arc. + +ADD = 1; +DEL = 2; +REV = 3; + +%op = [repmat(DEL, 1, E) repmat(REV, 1, E) repmat(ADD, 1, Ebar)]; diff --git a/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph_not_vectorized.m b/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph_not_vectorized.m new file mode 100644 index 00000000..d622943d --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_nbrs_of_digraph_not_vectorized.m @@ -0,0 +1,53 @@ +function [Gs, op, nodes] = mk_nbrs_of_digraph2(G0) +% MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal +% [Gs, op, nodes] = mk_nbrs_of_digraph(G0) +% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor. +% nodes(i,1:2) are the head and tail of the operated-on arc. + +[I,J] = find(G0); +G0bar = setdiag(~G0, 0); % exclude self loops in graph complement +[Ibar,Jbar] = find(G0bar); +nnbrs = 2*length(I) + length(Ibar); +Gs = cell(1, nnbrs); +op = cell(1, nnbrs); +nodes = zeros(nnbrs, 2); + +nbr = 1; +% all single edge deletions +for e=1:length(I) + i = I(e); j = J(e); + G = G0; + G(i,j) = 0; + Gs{nbr} = G; + op{nbr} = 'del'; + nodes(nbr, :) = [i j]; + nbr = nbr + 1; +end + +% all single edge reversals +for e=1:length(I) + i = I(e); j = J(e); + G = G0; + G(i,j) = 0; + G(j,i) = 1; + Gs{nbr} = G; + op{nbr} = 'rev'; + nodes(nbr, :) = [i j]; + nbr = nbr + 1; +end + +[I,J] = find(~G0); +% all single edge additions +for e=1:length(I) + i = I(e); j = J(e); + G = G0; + if i ~= j % don't add self loops + G(i,j) = 1; + Gs{nbr} = G; + op{nbr} = 'add'; + nodes(nbr, :) = [i j]; + nbr = nbr + 1; + end +end + +assert(nnbrs == nbr-1); diff --git a/sourcecodes/bnt-master/graph/mk_rnd_dag.m b/sourcecodes/bnt-master/graph/mk_rnd_dag.m new file mode 100644 index 00000000..4bf59c34 --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_rnd_dag.m @@ -0,0 +1,27 @@ +function [dag, order] = mk_rnd_dag(N, max_fan_in) +% MY_MK_RND_DAG Create a random directed acyclic graph +% +% [dag, order] = my_mk_rnd_dag(N, max_fan_in) +% max_fan_in defaults to N. +% order is the random topological order that was chosen + +% Modified by Sonia Leach 2/25/02 + +if nargin < 2, max_fan_in = N; end + +order = randperm(N); +dag = zeros(N,N); +for i=2:N + j = order(i); + %k = sample_discrete(normalise(ones(1, min(i-1, max_fan_in)))); + k = sample_discrete(normalise(ones(1, min(i-1, max_fan_in)+1))) - 1; % min = 0 (bug fix due to + % Pedrito, 7/28/04) + SS = order(1:i-1); % get Set of possible parentS + p = randperm(length(SS)); % permute order of set + dag(SS(p(1:k)),j) = 1; % take first k in permuted order + + % Kevin had: + %SS = subsets(order(1:i-1), k, k); + %p = sample_discrete(normalise(ones(1, length(SS)))); + %dag(SS{p}, j) = 1; +end diff --git a/sourcecodes/bnt-master/graph/mk_rnd_dag_given_edge_prob.m b/sourcecodes/bnt-master/graph/mk_rnd_dag_given_edge_prob.m new file mode 100644 index 00000000..205e5b51 --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_rnd_dag_given_edge_prob.m @@ -0,0 +1,15 @@ +function dag = sample_dag(P) +% SAMPLE_DAG Create a random directed acyclic graph with edge probabilities P(i,j) +% dag = sample_dag(P) +% +% This uses rejection sampling to reject graphs with directed cycles. + +done = 0; +directed = 1; +iter = 1; +while ~done + dag = binornd(1, P); % each edge is an indep Bernoulli (0/1) random variable + dag = setdiag(dag, 0); + done = acyclic(dag, directed); + iter = iter + 1 +end diff --git a/sourcecodes/bnt-master/graph/mk_rooted_tree.m b/sourcecodes/bnt-master/graph/mk_rooted_tree.m new file mode 100644 index 00000000..396043eb --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_rooted_tree.m @@ -0,0 +1,19 @@ +function [T, pre, post, cycle] = mk_rooted_tree(G, root) +% MK_ROOTED_TREE Make a directed tree pointing away from root +% [T, pre, post, cycle] = mk_rooted_tree(G, root) + +n = length(G); +T = sparse(n,n); % not the same as T = sparse(n) ! +directed = 0; +[d, pre, post, cycle, f, pred] = dfs(G, root, directed); +[junk, pre2] = sort(d); +assert(isequal(pre, pre2)) +[junk, post2] = sort(f); +assert(isequal(post, post2)); +%[d, pre, post, cycle, f, pred] = dfs(G, [], directed); +for i=1:length(pred) + if pred(i)>0 + T(pred(i),i)=1; + end +end + diff --git a/sourcecodes/bnt-master/graph/mk_undirected.m b/sourcecodes/bnt-master/graph/mk_undirected.m new file mode 100644 index 00000000..0a1614aa --- /dev/null +++ b/sourcecodes/bnt-master/graph/mk_undirected.m @@ -0,0 +1,11 @@ +function U = mk_undirected(G) + +[nr nc] = size(G); +U = G; +for i=1:nr + for j=1:nc + if U(i,j)==1 + U(j,i) = 1; + end + end +end diff --git a/sourcecodes/bnt-master/graph/moralize.m b/sourcecodes/bnt-master/graph/moralize.m new file mode 100644 index 00000000..59beff97 --- /dev/null +++ b/sourcecodes/bnt-master/graph/moralize.m @@ -0,0 +1,12 @@ +function [M, moral_edges] = moralize(G) +% MORALIZE Ensure that for every child, all its parents are married, and drop directionality of edges. +% [M, moral_edges] = moralize(G) + +M = G; +n = length(M); +for i=1:n + fam = family(G,i); + M(fam,fam)=1; +end +M = setdiag(M,0); +moral_edges = sparse(triu(max(0,M-G),1)); diff --git a/sourcecodes/bnt-master/graph/neighbors.m b/sourcecodes/bnt-master/graph/neighbors.m new file mode 100644 index 00000000..7275c4e0 --- /dev/null +++ b/sourcecodes/bnt-master/graph/neighbors.m @@ -0,0 +1,7 @@ +function ns = neighbors(adj_mat, i) +% NEIGHBORS Find the parents and children of a node in a graph. +% ns = neighbors(adj_mat, i) + +%ns = myunion(children(adj_mat, i), parents(adj_mat, i)); +ns = unique([find(adj_mat(i,:)) find(adj_mat(:,i))']); + diff --git a/sourcecodes/bnt-master/graph/parents.m b/sourcecodes/bnt-master/graph/parents.m new file mode 100644 index 00000000..24354a74 --- /dev/null +++ b/sourcecodes/bnt-master/graph/parents.m @@ -0,0 +1,5 @@ +function ps = parents(adj_mat, i) +% PARENTS Return the list of parents of node i +% ps = parents(adj_mat, i) + +ps = find(adj_mat(:,i))'; diff --git a/sourcecodes/bnt-master/graph/pred2path.m b/sourcecodes/bnt-master/graph/pred2path.m new file mode 100644 index 00000000..4763d12c --- /dev/null +++ b/sourcecodes/bnt-master/graph/pred2path.m @@ -0,0 +1,78 @@ +function rte = pred2path(P,s,t) +%PRED2PATH Convert predecessor indices to shortest paths from node 's' to 't'. +% rte = pred2path(P,s,t) +% P = |s| x n matrix of predecessor indices (from DIJK) +% s = FROM node indices +% = [] (default), paths from all nodes +% t = TO node indices +% = [] (default), paths to all nodes +% rte = |s| x |t| cell array of paths (or routes) from 's' to 't', where +% rte{i,j} = path from s(i) to t(j) +% = [], if no path exists from s(i) to t(j) +% +% (Used with output of DIJK) + +% Copyright (c) 1998-2001 by Michael G. Kay +% Matlog Version 5 22-Aug-2001 + +% Input Error Checking ****************************************************** +error(nargchk(1,3,nargin)); + +[rP,n] = size(P); + +if nargin < 2 | isempty(s), s = (1:n)'; else s = s(:); end +if nargin < 3 | isempty(t), t = (1:n)'; else t = t(:); end + +if any(P < 0 | P > n) + error(['Elements of P must be integers between 1 and ',num2str(n)]); +elseif any(s < 1 | s > n) + error(['''s'' must be an integer between 1 and ',num2str(n)]); +elseif any(t < 1 | t > n) + error(['''t'' must be an integer between 1 and ',num2str(n)]); +end +% End (Input Error Checking) ************************************************ + +rte = cell(length(s),length(t)); + +for i = 1:length(s) + if rP == 1 + si = 1; + else + si = s(i); + if si < 1 | si > rP + error('Invalid P matrix.') + end + end + for j = 1:length(t) + tj = t(j); + if tj == s(i) + r = tj; + elseif P(si,tj) == 0 + r = []; + else + r = tj; + while tj ~= s(i) + if tj < 1 | tj > n + error('Invalid element of P matrix found.') + end + r = [P(si,tj) r]; + tj = P(si,tj); + end + end + rte{i,j} = r; + end +end + +if length(s) == 1 & length(t) == 1 + rte = rte{:}; +end + +%rte = t; +while 0%t ~= s + if t < 1 | t > n | round(t) ~= t + error('Invalid ''pred'' element found prior to reaching ''s'''); + end + rte = [P(t) rte]; + t = P(t); +end + diff --git a/sourcecodes/bnt-master/graph/reachability_graph.m b/sourcecodes/bnt-master/graph/reachability_graph.m new file mode 100644 index 00000000..0f3c4f58 --- /dev/null +++ b/sourcecodes/bnt-master/graph/reachability_graph.m @@ -0,0 +1,19 @@ +function C = reachability_graph(G) +% REACHABILITY_GRAPH C(i,j) = 1 iff there is a path from i to j in DAG G +% C = reachability_graph(G) + +if 1 + % expm(G) = I + G + G^2 / 2! + G^3 / 3! + ... + M = expm(double(full(G))) - eye(length(G)); + C = (M>0); +else + % This computes C = G + G^2 + ... + G^{n-1} + n = length(G); + A = G; + C = zeros(n); + for i=1:n-1 + C = C + A; + A = A * G; + end + C = (C > 0); +end diff --git a/sourcecodes/bnt-master/graph/scc.m b/sourcecodes/bnt-master/graph/scc.m new file mode 100644 index 00000000..9a2e8b9d --- /dev/null +++ b/sourcecodes/bnt-master/graph/scc.m @@ -0,0 +1,64 @@ +function [c,v] = scc(a,tol) + +% Finds the strongly connected sets of vertices +% in the DI-rected G-raph of A +% c = 0-1 matrix displaying accessibility +% v = displays the equivalent classes +% +% v(i,j) is the j'th member of the i'th equiv class (0 padded) +% +% http://www.math.wsu.edu/math/faculty/tsat/matlab.html + +[m,n] = size(a); +if m~=n 'Not a Square Matrix', break, end +b=abs(a); o=ones(size(a)); x=zeros(1,n); +msg='The Matrix is Irreducible !'; v='Connected Directed Graph !'; +if (nargin==1) tol=n*eps*norm(a,'inf'); end + +% Create a companion matrix +c = b>tol*o; +if (c==o) + % msg, break + v = 1:length(a); + return +end + + +% Compute accessibility in at most n-step paths +for k=1:n + for j=1:n + for i=1:n + % If index i accesses j, where can you go ? + if c(i,j) > 0 c(i,:) = c(i,:)+c(j,:); end + end + end +end +% Create a 0-1 matrix with the above information +c>zeros(size(a)); c=ans; if (c==o) msg, break, end + +% Identify equivalence classes +d=c.*c'+eye(size(a)); d>zeros(size(a)); d=ans; +v=zeros(size(a)); +for i=1:n find(d(i,:)); ans(n)=0; v(i,:)=ans; end + +% Eliminate displaying of identical rows +i=1; +while(i<n) + for k=i+1:n + if v(k,1) == v(i,1) + v(k,:)=x; + end + end + i=i+1; +end +j=1; +for i=1:n + if v(i,1)>0 + h(j,:)=v(i,:); + j=j+1; + end +end +v=h; + + + diff --git a/sourcecodes/bnt-master/graph/strong_elim_order.m b/sourcecodes/bnt-master/graph/strong_elim_order.m new file mode 100644 index 00000000..3fd995d7 --- /dev/null +++ b/sourcecodes/bnt-master/graph/strong_elim_order.m @@ -0,0 +1,75 @@ +function order = strong_elim_order(G, node_sizes, partial_order) +% STRONG_ELIM_ORDER Find an elimination order to produce a strongly triangulated graph. +% order = strong_elim_order(moral_graph, node_sizes, partial_order) +% +% partial_order(i,j)=1 if we must marginalize i *after* j +% (so i will be nearer the strong root). +% e.g., if j is a decision node and i is its information set: +% we cannot maximize j if we have marginalized out some of i +% e.g., if j is a continuous child and i is its discrete parent: +% we want to integrate out the cts nodes before the discrete ones, +% so that the marginal is strong. +% +% For details, see +% - Jensen, Jensen and Dittmer, "From influence diagrams to junction trees", UAI 94. +% - Lauritzen, "Propgation of probabilities, means, and variances in mixed graphical +% association models", JASA 87(420):1098--1108, 1992. +% +% On p369 of the Jensen paper, they state "the reverse of the elimination order must be some +% extension of [the partial order] to a total order". +% We make no attempt to find the best such total ordering, in the sense of minimizing the weight +% of the resulting cliques. + +% Example from the Jensen paper: +% Let us number the nodes in Fig 1 from top to bottom, left to right, +% so a=1,b=2,D1=3,c=4,...,l=14,j=15,k=16. +% The elimination ordering they propose on p370 is [14 15 16 11 12 1 4 5 10 8 13 9 7 6 3 2]; + +if 0 + total_order = topological_sort(partial_order); + order = total_order(end:-1:1); % no attempt to find an optimal constrained ordering! + return; +end + +% The following implementation is due to Ilya Shpitser and seems to give wrong +% results on cg1 + +n = length(G); +MG = G; % copy the original graph +uneliminated = ones(1,n); +order = zeros(1,n); + +for i=1:n + roots = []; + k = 1; + for j=1:n + if sum(partial_order(j,:)) == 0 + roots(k) = j; + k = k + 1; + end + end + U = find(uneliminated); + valid = myintersect(U, roots); + % Choose the best node from the set of valid candidates + score1 = zeros(1,length(valid)); + score2 = zeros(1,length(valid)); + for j=1:length(valid) + k = valid(j); + ns = myintersect(neighbors(G, k), U); + l = length(ns); + M = MG(ns,ns); + score1(j) = l^2 - sum(M(:)); % num. added edges + score2(j) = prod(node_sizes([k ns])); % weight of clique + end + j1s = find(score1==min(score1)); + j = j1s(argmin(score2(j1s))); + k = valid(j); + uneliminated(k) = 0; + order(i) = k; + ns = myintersect(neighbors(G, k), U); + if ~isempty(ns) + G(ns,ns) = 1; + G = setdiag(G,0); + end + partial_order(:,k) = 0; +end diff --git a/sourcecodes/bnt-master/graph/test.m b/sourcecodes/bnt-master/graph/test.m new file mode 100644 index 00000000..6ffe3d2e --- /dev/null +++ b/sourcecodes/bnt-master/graph/test.m @@ -0,0 +1,16 @@ +% A - B +% | +% D - C - E + +A=1;B=2;C=3;D=4;E=5; +dag = zeros(5,5); +dag(A,B)=1; +%dag(A,D)=1; +dag(B,C)=1; +dag(C,D)=1; +dag(E,C)=1; +[d, pre, post, cycle, f, pred] = dfs(dag, A, 0) + +[T, pre, post, cycle] = mk_rooted_tree(dag, A) + +%[T, pre, post, cycle] = mkRootedTree(dag, A) diff --git a/sourcecodes/bnt-master/graph/test_strong_root.m b/sourcecodes/bnt-master/graph/test_strong_root.m new file mode 100644 index 00000000..86e1117e --- /dev/null +++ b/sourcecodes/bnt-master/graph/test_strong_root.m @@ -0,0 +1,34 @@ +function strong = test_strong_root(jtree,cliques,dnodes,root) +% This function tests, whether root is a strong root of jtree. +% The following parameters are used +% Input: +% jtree An matrix with two colums. jtree(i,j) == jtree(j,i) is 1 if node +% i is connected with node j +% cliques Cells which contain the nodes in each clique +% dnodes An array with the discrete nodes of the juntion tree. +% root It is tested whether root is the strong root of the junction tree +% Output: +% strong The output is 1 if root is the strong root of the junction tree jtree. +% Please note, that the running intersection property is not tested. +if isempty(dnodes) + strong = 1; + return; +end + +children = find(jtree(root,:)==1); +i = 1; +strong = 1; +while (i <= length(children)) & (strong==1) + child = children(i); + jtree(child,root) = 0; + jtree(root,child) = 0; + sep = myintersect(cliques{child},cliques{root}); + diff = mysetdiff(cliques{child},cliques{root}); + if (mysubset(sep,dnodes) | isempty(myintersect(diff,dnodes))) + strong = test_strong_root(jtree,cliques,dnodes,child); + else + strong = 0; + end; + i = i+1; +end + diff --git a/sourcecodes/bnt-master/graph/topological_sort.m b/sourcecodes/bnt-master/graph/topological_sort.m new file mode 100644 index 00000000..cd8b8325 --- /dev/null +++ b/sourcecodes/bnt-master/graph/topological_sort.m @@ -0,0 +1,30 @@ +function order = topological_sort(A) +% TOPOLOGICAL_SORT Return the nodes in topological order (parents before children). +% order = topological_sort(adj_mat) + +n = length(A); +indeg = zeros(1,n); +zero_indeg = []; % a stack of nodes with no parents +for i=1:n + indeg(i) = length(parents(A,i)); + if indeg(i)==0 + zero_indeg = [i zero_indeg]; + end +end + +t=1; +order = zeros(1,n); +while ~isempty(zero_indeg) + v = zero_indeg(1); % pop v + zero_indeg = zero_indeg(2:end); + order(t) = v; + t = t + 1; + cs = children(A, v); + for j=1:length(cs) + c = cs(j); + indeg(c) = indeg(c) - 1; + if indeg(c) == 0 + zero_indeg = [c zero_indeg]; % push c + end + end +end diff --git a/sourcecodes/bnt-master/graph/trees.txt b/sourcecodes/bnt-master/graph/trees.txt new file mode 100644 index 00000000..f5475aac --- /dev/null +++ b/sourcecodes/bnt-master/graph/trees.txt @@ -0,0 +1,168 @@ + +% make undirected adjacency matrix of graph/tree +% e.g., +% 1 +% / \ +% 2 3 +T = zeros(3,3); +T(1,2) = 1; T(2,1)=1; +T(1,3)=1; T(3,1) = 1; + +root = 1; +[T, preorder, postorder] = mk_rooted_tree(T, root); + +% bottom up message passing leaves to root +for n=postorder(:)' + for p = parents(T, n) + % p is parent of n + end +end + +% top down, root to leaves +for n=preorder(:)' + for c= children(T,n) + % c is child of n + end +end + + +%%%%%%%%%%%%% + +function ps = parents(adj_mat, i) +% PARENTS Return the list of parents of node i +% ps = parents(adj_mat, i) + +ps = find(adj_mat(:,i))'; + + +%%%%%%%%%%%% + +function cs = children(adj_mat, i, t) +% CHILDREN Return the indices of a node's children in sorted order +% c = children(adj_mat, i, t) +% +% t is an optional argument: if present, dag is assumed to be a 2-slice DBN + +if nargin < 3 + cs = find(adj_mat(i,:)); +else + if t==1 + cs = find(adj_mat(i,:)); + else + ss = length(adj_mat)/2; + j = i+ss; + cs = find(adj_mat(j,:)) + (t-2)*ss; + end +end + +%%%%%%%%%%% + +function [T, pre, post, cycle] = mk_rooted_tree(G, root) +% MK_ROOTED_TREE Make a directed, rooted tree out of an undirected tree. +% [T, pre, post, cycle] = mk_rooted_tree(G, root) + +n = length(G); +T = sparse(n,n); % not the same as T = sparse(n) ! +directed = 0; +[d, pre, post, cycle, f, pred] = dfs(G, root, directed); +for i=1:length(pred) + if pred(i)>0 + T(pred(i),i)=1; + end +end + + +%%%%%%%%%%% + +function [d, pre, post, cycle, f, pred] = dfs(adj_mat, start, directed) +% DFS Perform a depth-first search of the graph starting from 'start'. +% [d, pre, post, cycle, f, pred] = dfs(adj_mat, start, directed) +% +% Input: +% adj_mat(i,j)=1 iff i is connected to j. +% start is the root vertex of the dfs tree; if [], all nodes are searched +% directed = 1 if the graph is directed +% +% Output: +% d(i) is the time at which node i is first discovered. +% pre is a list of the nodes in the order in which they are first encountered (opened). +% post is a list of the nodes in the order in which they are last encountered (closed). +% 'cycle' is true iff a (directed) cycle is found. +% f(i) is the time at which node i is finished. +% pred(i) is the predecessor of i in the dfs tree. +% +% If the graph is a tree, preorder is parents before children, +% and postorder is children before parents. +% For a DAG, topological order = reverse(postorder). +% +% See Cormen, Leiserson and Rivest, "An intro. to algorithms" 1994, p478. + +n = length(adj_mat); + +global white gray black color +white = 0; gray = 1; black = 2; +color = white*ones(1,n); + +global time_stamp +time_stamp = 0; + +global d f +d = zeros(1,n); +f = zeros(1,n); + +global pred +pred = zeros(1,n); + +global cycle +cycle = 0; + +global pre post +pre = []; +post = []; + +if ~isempty(start) + dfs_visit(start, adj_mat, directed); +else + for u=1:n + if color(u)==white + dfs_visit(u, adj_mat, directed); + end + end +end + + +%%%%%%%%%% + +function dfs_visit(u, adj_mat, directed) + +global white gray black color time_stamp d f pred cycle pre post + +pre = [pre u]; +color(u) = gray; +time_stamp = time_stamp + 1; +d(u) = time_stamp; +if directed + ns = children(adj_mat, u); +else + ns = neighbors(adj_mat, u); + ns = mysetdiff(ns, pred(u)); % don't go back to visit the guy who called you! +end +for v=ns(:)' + %fprintf('u=%d, v=%d, color(v)=%d\n', u, v, color(v)) + switch color(v) + case white, % not visited v before (tree edge) + pred(v)=u; + dfs_visit(v, adj_mat, directed); + case gray, % back edge - v has been visited, but is still open + cycle = 1; + %fprintf('cycle: back edge from v=%d to u=%d\n', v, u); + case black, % v has been visited, but is closed + % no-op + end +end +color(u) = black; +post = [post u]; +time_stamp = time_stamp + 1; +f(u) = time_stamp; + + diff --git a/sourcecodes/bnt-master/graph/triangulate.c b/sourcecodes/bnt-master/graph/triangulate.c new file mode 100644 index 00000000..cda13533 --- /dev/null +++ b/sourcecodes/bnt-master/graph/triangulate.c @@ -0,0 +1,159 @@ +/* triangulate.c written by Ilya Shpitser */ + +#include <stdlib.h> + +#ifdef UNIX +#include "matlab.h" +#endif + +#include "matrix.h" +#include "mex.h" + +#include "elim.h" +#include "map.h" +#include "misc.h" + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + + int dims [2]; + int i, j, k, m, n; + long index; + double * G_pr; + double * stage_pr; + double * answer_G_pr, * fill_ins_pr; + double * matlab_clique_pr; + mxArray * matlab_clique; + Elimination e; + float ** adj_mat; + int ** order = (int **) NULL; + Iterator iter, iter2; + word w, w2; + int ** fill_ins; + Map cliques; + Map clique; + mxArray * fill_ins_mat; + int * nodes; + mxArray * full; + +// (original) full = mlfFull((mxArray *) prhs[0]); + full = (mxArray *) mlfFull((mxArray *) prhs[0]); // added typecasting + /* Obtain graph matrix information. */ + m = mxGetM(full); + n = mxGetN(full); + G_pr = mxGetPr(full); + + if(n < 1 || m < 1){ + return; + } + + /* Allocate and populate the log weight adjacency matrix corresponding + to the input graph. */ + adj_mat = (float **) malloc(sizeof(float *) * m); + adj_mat[0] = (float *) malloc(sizeof(float) * m * n); + for(i = 1; i < m; i++){ + adj_mat[i] = adj_mat[i - 1] + n; + } + /* We no longer have log weight info, but we have a (total) ordering on + the nodes already, so we do not need this information. */ + for(i = 0; i < m; i++){ + for(j = 0; j < n; j++){ + index = j * m + i; + if(G_pr[index] > 0){ + adj_mat[i][j] = 1; + } else { + adj_mat[i][j] = 0; + } + } + } + + /* Convert the total elimination ordering into a partial order argument + for the elimination routine. The elimination routine's purpose in this + mode of operation is to return cliques and fill-in edges. */ + if(nrhs > 1){ + order = (int **) malloc(sizeof(int *) * m); + order[0] = (int *) malloc(sizeof(int) * m * n); + for(i = 1; i < m; i++){ + order[i] = order[i - 1] + n; + } + for(i = 0; i < m; i++){ + for(j = 0; j < n; j++){ + order[i][j] = 0; + } + } + stage_pr = mxGetPr(prhs[1]); + for(i = 0; i < mxGetN(prhs[1]) - 1; i++){ + order[(int) stage_pr[i] - 1][(int) stage_pr[i + 1] - 1] = 1; + } + } + + /* Find the elimination ordering. */ + e = find_elim(n, adj_mat, order, -1); + + /* Allocate memory for the answer, and set the answer. */ + plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL); + answer_G_pr = mxGetPr(plhs[0]); + cliques = get_cliques(e); +/* + dims[0] = 1; + dims[1] = get_size_Map(cliques); + plhs[1] = mxCreateCellArray(2, (const int *) dims);*/ + plhs[1] = mxCreateCellMatrix(get_size_Map(cliques), 1); + fill_ins = get_fill_ins(e); + fill_ins_mat = mxCreateDoubleMatrix(m, n, mxREAL); + fill_ins_pr = mxGetPr(fill_ins_mat); + + for(i = 0; i < n; i++){ + for(j = 0; j < m; j++){ + index = j * m + i; + answer_G_pr[index] = G_pr[index]; + if(fill_ins[i][j] > 0){ + answer_G_pr[index] = 1; + fill_ins_pr[index] = 1; + } + } + } + mxDestroyArray(full); +// (original) plhs[2] = mlfSparse(fill_ins_mat, NULL, NULL, NULL, NULL, NULL); + plhs[2] = (mxArray *) mlfSparse(fill_ins_mat, NULL, NULL, NULL, NULL, NULL); // added typecasting + mxDestroyArray(fill_ins_mat); + nodes = (int *) malloc(sizeof(int) * n); + k = 0; + iter = get_Iterator(cliques); + while(!is_empty(iter)){ + w = next_key(iter); + clique = (Map) w.v; + matlab_clique = mxCreateDoubleMatrix(1, get_size_Map(clique), mxREAL); + matlab_clique_pr = mxGetPr(matlab_clique); + for(i = 0; i < n; i++){ + nodes[i] = 0; + } + iter2 = get_Iterator(clique); + while(!is_empty(iter2)){ + w2 = next_key(iter2); + nodes[w2.i] = w2.i + 1; + } + j = 0; + for(i = 0; i < n; i++){ + if(nodes[i] > 0){ + matlab_clique_pr[j++] = nodes[i]; + } + } + mxSetCell(plhs[1], k++, matlab_clique); + } + free(nodes); + + /* Finally, free the allocated memory. */ + destroy_Elimination(e); + if(adj_mat){ + if(adj_mat[0]){ + free(adj_mat[0]); + } + free(adj_mat); + } + if(order){ + if(order[0]){ + free(order[0]); + } + free(order); + } +} diff --git a/sourcecodes/bnt-master/graph/triangulate.m b/sourcecodes/bnt-master/graph/triangulate.m new file mode 100644 index 00000000..8725c62e --- /dev/null +++ b/sourcecodes/bnt-master/graph/triangulate.m @@ -0,0 +1,42 @@ +function [G, cliques, fill_ins] = triangulate(G, order) +% TRIANGULATE Ensure G is triangulated (chordal), i.e., every cycle of length > 3 has a chord. +% [G, cliques, fill_ins, cliques_containing_node] = triangulate(G, order) +% +% cliques{i} is the i'th maximal complete subgraph of the triangulated graph. +% fill_ins(i,j) = 1 iff we add a fill-in arc between i and j. +% +% To find the maximal cliques, we save each induced cluster (created by adding connecting +% neighbors) that is not a subset of any previously saved cluster. (A cluster is a complete, +% but not necessarily maximal, set of nodes.) + +MG = G; +n = length(G); +eliminated = zeros(1,n); +cliques = {}; +for i=1:n + u = order(i); + U = find(~eliminated); % uneliminated + nodes = myintersect(neighbors(G,u), U); % look up neighbors in the partially filled-in graph + nodes = myunion(nodes, u); % the clique will always contain at least u + G(nodes,nodes) = 1; % make them all connected to each other + G = setdiag(G,0); + eliminated(u) = 1; + + exclude = 0; + for c=1:length(cliques) + if mysubset(nodes,cliques{c}) % not maximal + exclude = 1; + break; + end + end + if ~exclude + cnum = length(cliques)+1; + cliques{cnum} = nodes; + end +end + +fill_ins = sparse(triu(max(0, G - MG), 1)); + +%assert(check_triangulated(G)); % takes 72% of the time! + + diff --git a/sourcecodes/bnt-master/graph/triangulate_2Dlattice_demo.m b/sourcecodes/bnt-master/graph/triangulate_2Dlattice_demo.m new file mode 100644 index 00000000..12ad1312 --- /dev/null +++ b/sourcecodes/bnt-master/graph/triangulate_2Dlattice_demo.m @@ -0,0 +1,131 @@ +% Consider a 3x3 lattice with 4-nearest neighbor connectivity + +% 1 - 2 - 3 +% | | | +% 4 - 5 - 6 +% | | | +% 7 - 8 - 9 + +N = 3; +G = mk_2D_lattice(N,N,4); +G0 = G; + +% Now add in the diagonal edges + +if 0 +% 1 - 2 - 3 +% | x | x | +% 4 - 5 - 6 +% | x | x | +% 7 - 8 - 9 + +G(1,5)=1; G(5,1)=1; +G(2,6)=1; G(6,2)=1; +G(4,2)=1; G(2,4)=1; +G(5,3)=1; G(3,5)=1; + +G(4,8)=1; G(8,4)=1; +G(5,9)=1; G(9,5)=1; +G(7,5)=1; G(5,7)=1; +G(8,6)=1; G(6,8)=1; +end + +% 1 - 2 - 3 +% | / | \ | +% 4 - 5 - 6 +% | \ | / | +% 7 - 8 - 9 + +G(2,6)=1; G(6,2)=1; +G(4,2)=1; G(2,4)=1; +G(4,8)=1; G(8,4)=1; +G(8,6)=1; G(6,8)=1; + +% Is this a chordal (triangulated) graph? No! + +assert(~check_triangulated(G)) + +% The reason is that there is a chordless cycle around the outside nodes. +% To see this, imagine "picking up" node 5, leaving the rest on the plane +% (like a hoop skirt, or a tent), as shown below + +% 1 - 2 - 3 +% | / \ | +% 4 6 +% | \ / | +% 7 - 8 - 9 + + +% However, if we add in the 4-6 arc, it will be chordal. + +G2 = G; +G2(4,6)=1; G2(6,4)=1; +assert(check_triangulated(G2)) + +% Or we can add in the 2-8 arc +G2 = G; +G2(2,8)=1; G2(8,2)=1; +assert(check_triangulated(G2)) + + +if 0 +% 4x4 lattice with cross arcs +N=4;G0 = mk_2D_lattice(N,N,4); +vs = [1 6; 2 5; 2 7; 3 6; 3 8; 4 7; ... + 5 10; 6 9; 6 11; 7 10; 7 12; 8 11;... + 9 14; 10 13; 10 15; 11 14; 11 16; 12 15]; +for i=1:size(vs,1) + u = vs(i,1); v= vs(i,2); + G0(u,v) = 1; G0(v,u) = 1; +end +end + +% Here is how we can discover which edges to fill in automatically +% (although possibly sub-optimally) +weights = 2*ones(1,N*N); % all nodes are binar + +% fill-ins = 2-4, 2-6, 4-8, 6-8 and 4-6 +% cliques = 124, etc and 2456 4568 +greedy_order = best_first_elim_order(G0, weights); +[GT, cliques, fill_ins] = triangulate(G0, greedy_order) +assert(check_triangulated(GT)) + + + +greedy_order = best_first_elim_order(G, weights); +[GT, cliques, fill_ins] = triangulate(G, greedy_order) +assert(check_triangulated(GT)) + +% fill-ins = [4 6] + +% Cliques are the overlapping squares [1,2,4,5], [2 3 5 6], [4 5 7 8], [5 6 8 9] +% and the following caused by the fill-in: [2 4 5 6], [4 5 6 8] + +% Connect the maximal cliques of the triangulate graph into a junction tree +[jtree, root, B, clq_weights] = cliques_to_jtree(cliques, weights); + +% In this case, all cliques have weight 2^4 = 16 + + +% Now consider size of max clique as a function of grid size +% Note: this is not necessarily the optimal triangulation + +% N 5 10 15 16 17 18 +% m 6 15 23 25 28 28 +Ns = [5 10 15 16 17 18]; +for i=1:length(Ns) + N = Ns(i) + G = mk_2D_lattice(N,N,4); + weights = 2*ones(1,N*N); % all nodes are binary + greedy_order = best_first_elim_order(G, weights); % slow! + [GT, cliques, fill_ins] = triangulate(G, greedy_order); + %assert(check_triangulated(GT)) + [jtree, root, B, clq_weights] = cliques_to_jtree(cliques, weights); + m(i) = log2(max(clq_weights)); +end + +% plot distribution of clique sizes for fixed N +for c=1:length(cliques) + l(c) = length(cliques{c}); +end +hist(l) diff --git a/sourcecodes/bnt-master/graph/triangulate_test.m b/sourcecodes/bnt-master/graph/triangulate_test.m new file mode 100644 index 00000000..70083d1c --- /dev/null +++ b/sourcecodes/bnt-master/graph/triangulate_test.m @@ -0,0 +1,36 @@ +% Test the code using the dag in Fig 1 of Jensen, Jensen, Dittmer, +% "From influence diagrams to junction trees", UAI 94 + +% By reverse enginering Fig 2, we infer that the following arcs should +% be absent from the original dag: b->d1, e->d2, f->d2, g->d4 +a=1; b=2; d1=3; c=4; d=5; e=6; f=7; g=8; d2=9; d4=10; i=11; h=12; d3=13; l=14; j=15; k=16; +dag=zeros(16); +dag(a,c)=1; +%dag(b,[c d d1])=1; +dag(b,[c d])=1; +dag(d1,d)=1; +dag(c,e)=1; +dag(d,[e f])=1; +%dag(e,[g d2])=1; +dag(e,[g])=1; +%dag(f,[d2 h])=1; +dag(f,[h])=1; +%dag(g,[d4 i])=1; +dag(g,[i])=1; +dag(d2,i)=1; +dag(d4,l)=1; +dag(i,l)=1; +dag(h,[j k])=1; +dag(d3,k)=1; + + +[MG, moral_edges] = moralize(dag); +MG(j,k)=1; MG(k,j)=1; % simulate having a common utility child +% MG now equals fig 2 +order = [l j k i h a c d d4 g d3 d2 f e d1 b]; +[MTG, cliques, fill_ins] = triangulate(MG, order); +% MTG equals fig 3 +ns = 2*ones(1,16); +[jtree, root, cliques2] = mk_strong_jtree(cliques, ns, order, MTG); +jtree2 = mk_rooted_tree(jtree, root); +% jtree2 equals fig 4, with their arrows reversed |