diff options
| author | ziejd2 | 2017-09-28 15:04:40 -0500 |
|---|---|---|
| committer | ziejd2 | 2017-09-28 15:04:40 -0500 |
| commit | 8070dc963753142bb86c4ed698d91fd623ed28e7 (patch) | |
| tree | d0f6dd8fc46a49b819aa55c1a90faa14d8448883 /sourcecodes/bnt-master/graph | |
| parent | 7cc31810d53176e805532b2789955f4eedbce6bb (diff) | |
| download | BNW-8070dc963753142bb86c4ed698d91fd623ed28e7.tar.gz | |
BNW using Octave instead of Matlab.
This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning. I am calling this BNW_1.02. It can be accessed at: compbio.uthsc.edu/BNW_1.02
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 |