Added commented-out original code to still-empty functions

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;
+}

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
+}

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
+}

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,:));
+}

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;
+}

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))
+}

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()