From 17019a9c94bc4b1a9f70492cf4076c54123378a7 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Thu, 22 Dec 2022 11:48:49 +0100 Subject: [PATCH] Added commented-out original code to still-empty functions --- R/cliques_to_jtree.R | 59 +++++++++++++++++++- R/elim_order.R | 72 ++++++++++++++++++++++++- R/myunion.R | 32 +++++++++++ R/neighbors.R | 9 ++++ R/setdiag.R | 21 ++++++++ R/triangulate.R | 80 +++++++++++++++++++++++++++- tests/testthat/test-spatialMixture.R | 1 + 7 files changed, 271 insertions(+), 3 deletions(-) create mode 100644 R/myunion.R create mode 100644 R/neighbors.R create mode 100644 R/setdiag.R diff --git a/R/cliques_to_jtree.R b/R/cliques_to_jtree.R index 6b3fd18..40c6f3b 100644 --- a/R/cliques_to_jtree.R +++ b/R/cliques_to_jtree.R @@ -1 +1,58 @@ -cliques_to_jtree <- function(cliques, ns) {} +cliques_to_jtree <- function(cliques, ns) { + stop("needs translation") + # 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/R/elim_order.R b/R/elim_order.R index fe89c56..9d9090e 100644 --- a/R/elim_order.R +++ b/R/elim_order.R @@ -1 +1,71 @@ -elim_order <- function(G, node_sizes) {} +elim_order <- function(G, node_sizes) { + stop("needs translation") + # function order = elim_order(G, node_sizes) + # % 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}); + # valid = U; + # % 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/R/myunion.R b/R/myunion.R new file mode 100644 index 0000000..e359aa5 --- /dev/null +++ b/R/myunion.R @@ -0,0 +1,32 @@ +myunion <- function(A, B) { + stop("needs translation") + # function C = myunion(A,B) + # % MYUNION Union of two sets of positive integers (much faster than built-in union) + # % C = myunion(A,B) + + # if isempty(A) + # ma = 0; + # else + # ma = max(A); + # end + + # if isempty(B) + # mb = 0; + # else + # mb = max(B); + # end + + # if ma==0 & mb==0 + # C = []; + # elseif ma==0 & mb>0 + # C = B; + # elseif ma>0 & mb==0 + # C = A; + # else + # %bits = sparse(1, max(ma,mb)); + # bits = zeros(1, max(ma,mb)); + # bits(A) = 1; + # bits(B) = 1; + # C = find(bits); + # end +} diff --git a/R/neighbors.R b/R/neighbors.R new file mode 100644 index 0000000..777ff7c --- /dev/null +++ b/R/neighbors.R @@ -0,0 +1,9 @@ +neighbors <- function(add_mat, i) { + stop("needs translation") + # 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 = find(adj_mat(i,:)); +} diff --git a/R/setdiag.R b/R/setdiag.R new file mode 100644 index 0000000..b77f76d --- /dev/null +++ b/R/setdiag.R @@ -0,0 +1,21 @@ +setdiag <- function(M, v) { + stop("needs translation") + # function M = setdiag(M, v) + # % SETDIAG Set the diagonal of a matrix to a specified scalar/vector. + # % M = set_diag(M, v) + + # n = length(M); + # if length(v)==1 + # v = repmat(v, 1, n); + # end + + # % e.g., for 3x3 matrix, elements are numbered + # % 1 4 7 + # % 2 5 8 + # % 3 6 9 + # % so diagnoal = [1 5 9] + + + # J = 1:n+1:n^2; + # M(J) = v; +} diff --git a/R/triangulate.R b/R/triangulate.R index 0280b21..2943f06 100644 --- a/R/triangulate.R +++ b/R/triangulate.R @@ -1 +1,79 @@ -triangulate <- function(G, order) {} +triangulate <- function(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 = list() + for (i in 1:n) { + u <- order[i] + U <- find(!eliminated)# uneliminated + nodes <- myintersect(neighbors(G, u), U)# look up neighbors in the partially filled - in graph # TODO: translate neighbors + nodes <- myunion(nodes, u)# the clique will always contain at least u # TODO: translate myunion + ----------------------- line 21 ----------------------- + G(nodes,nodes) = 1; % make them all connected to each other + G[nodes, nodes] <- 1# make them all connected to each other + ----------------------- line 22 ----------------------- + G = setdiag(G,0); + G <- setdiag(G, 0) + ----------------------- line 23 ----------------------- + eliminated(u) = 1; + eliminated[u] <- 1 + ----------------------- line 24 ----------------------- + + + ----------------------- line 25 ----------------------- + exclude = 0; + exclude <- 0 + ----------------------- line 26 ----------------------- + for c=1:length(cliques) + for (c in 1:length(cliques)) { + ----------------------- line 27 ----------------------- + if mysubset(nodes,cliques{c}) % not maximal + if (mysubset(nodes, cliques{c})# not maximal) { + ----------------------- line 28 ----------------------- + exclude = 1; + exclude <- 1 + ----------------------- line 29 ----------------------- + break; + break + ----------------------- line 30 ----------------------- + end + } + ----------------------- line 31 ----------------------- + end + } + ----------------------- line 32 ----------------------- + if ~exclude + if (!exclude) { + ----------------------- line 33 ----------------------- + cnum = length(cliques)+1; + cnum <- length(cliques) + 1 + ----------------------- line 34 ----------------------- + cliques{cnum} = nodes; + cliques{cnum} = nodes + ----------------------- line 35 ----------------------- + end + } + ----------------------- line 36 ----------------------- + end + } + ----------------------- line 37 ----------------------- + + + ----------------------- line 38 ----------------------- + %fill_ins = sparse(triu(max(0, G - MG), 1)); + # fill_ins <- sparse(triu(max(0, G - MG), 1)) + ----------------------- line 39 ----------------------- + fill_ins=1; + fill_ins=1 + ----------------------- line 40 ----------------------- + NA + return(list("G" = G, "cliques" = cliques, "fill_ins" = fill_ins)) +} diff --git a/tests/testthat/test-spatialMixture.R b/tests/testthat/test-spatialMixture.R index 8519bdc..8f25500 100644 --- a/tests/testthat/test-spatialMixture.R +++ b/tests/testthat/test-spatialMixture.R @@ -37,6 +37,7 @@ test_that("testaaKoordinaatit works as expected", { test_that("lakseKlitik() and subfunctions produce expected output", { # TODO: test elim_order() # TODO: test triangulate() + # TODO: test neighbors() # TODO: test myintersect() # TODO: test findCliques() # TODO: test cliques_to_jtree()