From b3fd9612091ed9356791531aee2a4a5ccbd703ec Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Wed, 21 Dec 2022 14:26:55 +0100
Subject: [PATCH 01/23] Updated RoxygenNote version

---
 DESCRIPTION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 557b643..47c0789 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -36,7 +36,7 @@ Description: Partial R implementation of the BAPS software
 License: GPL-3
 BugReports: https://github.com/ocbe-uio/rBAPS/issues
 Encoding: UTF-8
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Suggests:
     testthat (>= 2.1.0)
 Imports:

From 8f167f2e432d5994e1f1cb8dcdf2c6b077646e57 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Wed, 21 Dec 2022 14:27:29 +0100
Subject: [PATCH 02/23] Translated laskeKlitik()

Also, added placeholder for subfunction findCliques()
---
 R/findCliques.R |  3 +++
 R/laskeKlitik.R | 28 ++++++++++++++++++++++++++++
 2 files changed, 31 insertions(+)
 create mode 100644 R/findCliques.R
 create mode 100644 R/laskeKlitik.R

diff --git a/R/findCliques.R b/R/findCliques.R
new file mode 100644
index 0000000..eea6d3c
--- /dev/null
+++ b/R/findCliques.R
@@ -0,0 +1,3 @@
+findCliques <- function(M) {
+   # TODO: translate findCliques() from matlab/spatial/findCliques.m
+}
diff --git a/R/laskeKlitik.R b/R/laskeKlitik.R
new file mode 100644
index 0000000..7e09d52
--- /dev/null
+++ b/R/laskeKlitik.R
@@ -0,0 +1,28 @@
+laskeKlikit <- function(M, maxCliqSize, maxSepSize) {
+  # Laskee samankokoisten klikkien mההrהn verkosta M
+  # ncliques(i)=kokoa i olevien klikkien mההr?
+  # nseparators vastaavasti
+  ncliques <- zeros(1, maxCliqSize)
+  nseparators <- zeros(1, maxSepSize)
+  if (M == c()) {
+    return()
+  }
+  cliques_separators <- findCliques(M)
+  cliques <- cliques_separators$cliques
+  separators <- cliques_separators$separators
+  rm(cliques_separators)
+  for (i in 1:length(cliques)) {
+    ncliques[length[cliques[[i]]]] <- ncliques[length(cliques[[i]])] + 1
+  }
+  # cliqmax=max(find(ncliques!=0))
+  # ncliques=ncliques(1:cliqmax)
+  for (i in 1:length(separators)) {
+    nseparators[length[separators[[i]]]] <- nseparators[length(separators[[i]])] + 1
+  }
+  return(
+    list(
+      ncliques = ncliques, nseparators = nseparators, cliques = cliques,
+      separators = separators
+    )
+  )
+}

From 80ddbbdd00967fdb428ea5ca08971411ae995391 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Wed, 21 Dec 2022 14:40:12 +0100
Subject: [PATCH 03/23] Importing zeallot package

Its unpacking operator, %<-%, looks very promising for DRYing the codebase, which uses lots of multiple assignments.
---
 DESCRIPTION       | 2 +-
 R/rBAPS-package.R | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 47c0789..36476ef 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -40,4 +40,4 @@ RoxygenNote: 7.2.3
 Suggests:
     testthat (>= 2.1.0)
 Imports:
-    methods, ape, vcfR, Rsamtools, adegenet, matlab2r, R6
+    methods, ape, vcfR, Rsamtools, adegenet, matlab2r, R6, zeallot
diff --git a/R/rBAPS-package.R b/R/rBAPS-package.R
index 3867ba3..9def117 100644
--- a/R/rBAPS-package.R
+++ b/R/rBAPS-package.R
@@ -9,4 +9,5 @@
 #' isempty isfield isspace max min ones questdlg rand repmat reshape
 #' size sortrows squeeze strcmp times zeros disp
 #' @importFrom stats runif
+#' @importFrom zeallot %<-%
 NULL

From 17a5ff66d6b588a4b2d6089f392c162bdc2f97c6 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 08:29:50 +0100
Subject: [PATCH 04/23] Translated findCliques() + subfunctions

---
 NAMESPACE            |  1 +
 R/cliques_to_jtree.R |  1 +
 R/elim_order.R       |  1 +
 R/findCliques.R      | 46 +++++++++++++++++++++++++++++++++++++++++++-
 R/myintersect.R      |  1 +
 R/triangulate.R      |  1 +
 6 files changed, 50 insertions(+), 1 deletion(-)
 create mode 100644 R/cliques_to_jtree.R
 create mode 100644 R/elim_order.R
 create mode 100644 R/myintersect.R
 create mode 100644 R/triangulate.R

diff --git a/NAMESPACE b/NAMESPACE
index 3b8b822..f0dc387 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -74,3 +74,4 @@ importFrom(stats,runif)
 importFrom(stats,sd)
 importFrom(utils,read.delim)
 importFrom(vcfR,read.vcfR)
+importFrom(zeallot,"%<-%")
diff --git a/R/cliques_to_jtree.R b/R/cliques_to_jtree.R
new file mode 100644
index 0000000..6b3fd18
--- /dev/null
+++ b/R/cliques_to_jtree.R
@@ -0,0 +1 @@
+cliques_to_jtree <- function(cliques, ns) {}
diff --git a/R/elim_order.R b/R/elim_order.R
new file mode 100644
index 0000000..fe89c56
--- /dev/null
+++ b/R/elim_order.R
@@ -0,0 +1 @@
+elim_order <- function(G, node_sizes) {}
diff --git a/R/findCliques.R b/R/findCliques.R
index eea6d3c..4ff1d8b 100644
--- a/R/findCliques.R
+++ b/R/findCliques.R
@@ -1,3 +1,47 @@
 findCliques <- function(M) {
-   # TODO: translate findCliques() from matlab/spatial/findCliques.m
+  # Muuttaa graafin M kolmioituvaksi ja laskee siitה klikit ja
+  # separaattorit.
+  # Hyצdynnetההn Kevin Murphyn algoritmeja Graph Theory toolboxista.
+  # Pהivitetty 12.8.2005
+  order <- elim_order(M, ones(length(M))) # TODO: translate from findCliques.m
+  c(G, cliques) %<-% triangulate(M, order) # TODO: translate from findCliques.m
+  c(jtree, root) %<-% cliques_to_jtree(cliques, ones(length(M))) # TODO: translate from findCliques.m
+  ncliq <- length(cliques)
+  separators <- cell(ncliq - 1, 1)    # n - solmuisessa puussa n - 1 viivaa
+
+  jono <- zeros(length(ncliq))
+  jono[1] <- root
+  i <- 1
+  pointer <- 2                     # Seuraava tyhjה paikka
+
+  while (!is.null(find(jono != 0))) {  # Puun leveyssuuntainen lהpikהynti) {
+    lapset <- find(jtree[jono[i], ] != 0)
+    jtree[, jono[i]] <- 0         # Klikki kהsitelty
+    jono[pointer:(pointer + length(lapset) - 1)] <- lapset
+    for (j in 1:length(lapset)) {
+      ehdokas <- myintersect(cliques[[jono[i]]], cliques[[lapset[j]]])
+      kelpaa <- 1
+      for (k in 1:(pointer + j - 3)) {
+        # Tutkitaan, ettה separaattoriehdokasta ei vielה kהsitelty
+        if (ehdokas == separators[[k]]) {
+          kelpaa <- 0
+        }
+      }
+      if (kelpaa) {
+        separators[[pointer + j - 2]] <- ehdokas
+      }
+    }
+    jono[i] <- 0
+    pointer <- pointer + length(lapset)
+    i <- i + 1
+  }
+  notEmpty <- zeros(ncliq - 1, 1)
+  for (i in 1:(ncliq - 1)) {
+    if (!is.null(separators[[i]])) {
+      notEmpty[i] <- 1
+    }
+  }
+  notEmpty <- find(notEmpty == 1)
+  separators <- separators(notEmpty)
+  return(list("cliques" = cliques, "separators" = separators, "G" = G))
 }
diff --git a/R/myintersect.R b/R/myintersect.R
new file mode 100644
index 0000000..b2c19f1
--- /dev/null
+++ b/R/myintersect.R
@@ -0,0 +1 @@
+myintersect <- function(A, B) {}
diff --git a/R/triangulate.R b/R/triangulate.R
new file mode 100644
index 0000000..0280b21
--- /dev/null
+++ b/R/triangulate.R
@@ -0,0 +1 @@
+triangulate <- function(G, order) {}

From e094a8d31720286f2e1cfad085610d0e310f7d5f Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 10:00:21 +0100
Subject: [PATCH 05/23] Workaround for variable scope check note

Read https://github.com/r-lib/zeallot/issues/57 for more details
---
 R/findCliques.R | 1 +
 1 file changed, 1 insertion(+)

diff --git a/R/findCliques.R b/R/findCliques.R
index 4ff1d8b..ba0ec59 100644
--- a/R/findCliques.R
+++ b/R/findCliques.R
@@ -4,6 +4,7 @@ findCliques <- function(M) {
   # Hyצdynnetההn Kevin Murphyn algoritmeja Graph Theory toolboxista.
   # Pהivitetty 12.8.2005
   order <- elim_order(M, ones(length(M))) # TODO: translate from findCliques.m
+  G <- cliques <- root <- NULL
   c(G, cliques) %<-% triangulate(M, order) # TODO: translate from findCliques.m
   c(jtree, root) %<-% cliques_to_jtree(cliques, ones(length(M))) # TODO: translate from findCliques.m
   ncliq <- length(cliques)

From 35bb795dbe70015795931c911dabdd83132269e8 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 10:01:26 +0100
Subject: [PATCH 06/23] Added (empty) unit test for lakseKlitik()

Helps remember to add tests for functions
---
 tests/testthat/test-spatialMixture.R | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/tests/testthat/test-spatialMixture.R b/tests/testthat/test-spatialMixture.R
index 9203d56..8519bdc 100644
--- a/tests/testthat/test-spatialMixture.R
+++ b/tests/testthat/test-spatialMixture.R
@@ -33,3 +33,12 @@ test_that("testaaKoordinaatit works as expected", {
       )
   )
 })
+
+test_that("lakseKlitik() and subfunctions produce expected output", {
+  # TODO: test elim_order()
+  # TODO: test triangulate()
+  # TODO: test myintersect()
+  # TODO: test findCliques()
+  # TODO: test cliques_to_jtree()
+  # TODO: test lakseKlitik()
+})

From e60e40b21cbf1616cc7fab917e7b7a3c08c4cab1 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 10:41:26 +0100
Subject: [PATCH 07/23] Translated myintersect()

---
 R/myintersect.R | 30 +++++++++++++++++++++++++++++-
 1 file changed, 29 insertions(+), 1 deletion(-)

diff --git a/R/myintersect.R b/R/myintersect.R
index b2c19f1..e2e8f26 100644
--- a/R/myintersect.R
+++ b/R/myintersect.R
@@ -1 +1,29 @@
-myintersect <- function(A, B) {}
+myintersect <- function(A, B) {
+  #  MYINTERSECT Intersection of two sets of positive integers (much faster than built - in intersect)
+  #  C <- myintersect(A, B)
+
+  A <- t(A)
+  B <- t(B)
+
+  if (is.null(A)) {
+    ma <- 0
+  } else {
+    ma <- max(A)
+  }
+
+  if (is.null(B)) {
+    mb <- 0
+  } else {
+    mb <- max(B)
+  }
+
+  if (ma == 0 | mb == 0) {
+    C <- vector()
+  } else {
+    # bits <- sparse(1, max(ma, mb))
+    bits <- zeros(1, max(ma, mb))
+    bits[A] <- 1
+    C <- B[as.logical(bits[B])]
+  }
+  return(C)
+}

From 17019a9c94bc4b1a9f70492cf4076c54123378a7 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 11:48:49 +0100
Subject: [PATCH 08/23] 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()

From eabfd2efea1bee0aa3406e608ad74dc5db6f0654 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 11:57:32 +0100
Subject: [PATCH 09/23] Translated triangulate()

---
 R/triangulate.R | 54 +++++++------------------------------------------
 1 file changed, 7 insertions(+), 47 deletions(-)

diff --git a/R/triangulate.R b/R/triangulate.R
index 2943f06..145869e 100644
--- a/R/triangulate.R
+++ b/R/triangulate.R
@@ -16,64 +16,24 @@ triangulate <- function(G, order) {
     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;
+    G <- setdiag(G, 0) # TODO: translate setdiag
     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
+      if (mysubset(nodes, cliques[[c]])) { # not maximal)
+        exclude <- 1
+        break
+      }
     }
-  ----------------------- 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
+      cliques[[cnum]] <- nodes
     }
-  ----------------------- 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
+  fill_ins <- 1
   return(list("G" = G, "cliques" = cliques, "fill_ins" = fill_ins))
 }

From d4c11a89af83bea869e6d815b4f8c98f3176d7db Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 13:12:40 +0100
Subject: [PATCH 10/23] Translated neighbors()

---
 R/neighbors.R   | 13 ++++++-------
 R/triangulate.R |  2 +-
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/R/neighbors.R b/R/neighbors.R
index 777ff7c..9192827 100644
--- a/R/neighbors.R
+++ b/R/neighbors.R
@@ -1,9 +1,8 @@
-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)
+neighbors <- function(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,:));
+  # ns <- myunion(children(adj_mat, i), parents(adj_mat, i))
+  ns <- find(adj_mat[i, ])
+  return(ns)
 }
diff --git a/R/triangulate.R b/R/triangulate.R
index 145869e..7976fc1 100644
--- a/R/triangulate.R
+++ b/R/triangulate.R
@@ -14,7 +14,7 @@ triangulate <- function(G, order) {
   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 <- 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 # TODO: translate myunion
     G[nodes, nodes] <- 1# make them all connected to each other
     G <- setdiag(G, 0) # TODO: translate setdiag

From bc5bb61f5263a91da3ddb81a5b8ae61b1066a1dc Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 13:18:54 +0100
Subject: [PATCH 11/23] Added placeholder for mysubset()

---
 R/mysubset.R                         | 13 +++++++++++++
 R/triangulate.R                      |  2 +-
 tests/testthat/test-spatialMixture.R |  1 +
 3 files changed, 15 insertions(+), 1 deletion(-)
 create mode 100644 R/mysubset.R

diff --git a/R/mysubset.R b/R/mysubset.R
new file mode 100644
index 0000000..fd0a3f5
--- /dev/null
+++ b/R/mysubset.R
@@ -0,0 +1,13 @@
+mysubset <- function(small, large) {
+# function p=mysubset(small,large)
+# % MYSUBSET Is the small set of +ve integers a subset of the large set?
+# % p = mysubset(small, large)
+
+# % Surprisingly, this is not built-in.
+
+# if isempty(small)
+#   p = 1; % isempty(large);
+# else
+#   p = length(myintersect(small,large)) == length(small);
+# end
+}
diff --git a/R/triangulate.R b/R/triangulate.R
index 7976fc1..36d1410 100644
--- a/R/triangulate.R
+++ b/R/triangulate.R
@@ -22,7 +22,7 @@ triangulate <- function(G, order) {
 
     exclude <- 0
     for (c in 1:length(cliques)) {
-      if (mysubset(nodes, cliques[[c]])) { # not maximal)
+      if (mysubset(nodes, cliques[[c]])) { # not maximal) # TODO: translate mysubset
         exclude <- 1
         break
       }
diff --git a/tests/testthat/test-spatialMixture.R b/tests/testthat/test-spatialMixture.R
index 8f25500..4bf01e6 100644
--- a/tests/testthat/test-spatialMixture.R
+++ b/tests/testthat/test-spatialMixture.R
@@ -39,6 +39,7 @@ test_that("lakseKlitik() and subfunctions produce expected output", {
   # TODO: test triangulate()
   # TODO: test neighbors()
   # TODO: test myintersect()
+  # TODO: test mysubset()
   # TODO: test findCliques()
   # TODO: test cliques_to_jtree()
   # TODO: test lakseKlitik()

From 73438af7bbfc2632679284d544177c06d7f9d146 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 13:29:54 +0100
Subject: [PATCH 12/23] Removed translating TODOs

Replaced with built-in TODO list in Issue #3's body
---
 R/findCliques.R | 6 +++---
 R/triangulate.R | 6 +++---
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/R/findCliques.R b/R/findCliques.R
index ba0ec59..f413e64 100644
--- a/R/findCliques.R
+++ b/R/findCliques.R
@@ -3,10 +3,10 @@ findCliques <- function(M) {
   # separaattorit.
   # Hyצdynnetההn Kevin Murphyn algoritmeja Graph Theory toolboxista.
   # Pהivitetty 12.8.2005
-  order <- elim_order(M, ones(length(M))) # TODO: translate from findCliques.m
+  order <- elim_order(M, ones(length(M)))
   G <- cliques <- root <- NULL
-  c(G, cliques) %<-% triangulate(M, order) # TODO: translate from findCliques.m
-  c(jtree, root) %<-% cliques_to_jtree(cliques, ones(length(M))) # TODO: translate from findCliques.m
+  c(G, cliques) %<-% triangulate(M, order)
+  c(jtree, root) %<-% cliques_to_jtree(cliques, ones(length(M)))
   ncliq <- length(cliques)
   separators <- cell(ncliq - 1, 1)    # n - solmuisessa puussa n - 1 viivaa
 
diff --git a/R/triangulate.R b/R/triangulate.R
index 36d1410..47bb409 100644
--- a/R/triangulate.R
+++ b/R/triangulate.R
@@ -15,14 +15,14 @@ triangulate <- function(G, order) {
     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 # TODO: translate myunion
+    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) # TODO: translate setdiag
+    G <- setdiag(G, 0)
     eliminated[u] <- 1
 
     exclude <- 0
     for (c in 1:length(cliques)) {
-      if (mysubset(nodes, cliques[[c]])) { # not maximal) # TODO: translate mysubset
+      if (mysubset(nodes, cliques[[c]])) { # not maximal)
         exclude <- 1
         break
       }

From a757a4435c2de13043228a507b386033e68f5957 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 13:31:41 +0100
Subject: [PATCH 13/23] Translated mysubset()

---
 R/mysubset.R | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/R/mysubset.R b/R/mysubset.R
index fd0a3f5..1e74c86 100644
--- a/R/mysubset.R
+++ b/R/mysubset.R
@@ -1,13 +1,13 @@
 mysubset <- function(small, large) {
-# function p=mysubset(small,large)
-# % MYSUBSET Is the small set of +ve integers a subset of the large set?
-# % p = mysubset(small, large)
+  #  MYSUBSET Is the small set of + ve integers a subset of the large set?
+  #  p <- mysubset(small, large)
 
-# % Surprisingly, this is not built-in.
+  #  Surprisingly, this is not built - in.
 
-# if isempty(small)
-#   p = 1; % isempty(large);
-# else
-#   p = length(myintersect(small,large)) == length(small);
-# end
+  if (is.null(small)) {
+    p <- 1# is.null(large)
+  } else {
+    p <- length(myintersect(small, large)) == length(small)
+  }
+  return(p)
 }

From f2f59ba7f81c8cd15d8486bf17e0215f4efb8b91 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 13:36:07 +0100
Subject: [PATCH 14/23] Translated myunion()

---
 R/myunion.R | 53 ++++++++++++++++++++++++++---------------------------
 1 file changed, 26 insertions(+), 27 deletions(-)

diff --git a/R/myunion.R b/R/myunion.R
index e359aa5..9d64098 100644
--- a/R/myunion.R
+++ b/R/myunion.R
@@ -1,32 +1,31 @@
 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)
+  #  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 (is.null(A)) {
+    ma <- 0
+  } else {
+    ma <- max(A)
+  }
 
-  # if isempty(B)
-  #   mb = 0;
-  # else
-  #   mb = max(B);
-  # end
+  if (is.null(B)) {
+    mb <- 0
+  } else {
+    mb <- max(B)
+  }
 
-  # 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
+  if (ma == 0 & mb == 0) {
+    C <- vector()
+  } else if (ma == 0 & mb > 0) {
+    C <- B
+  } else if (ma > 0 & mb == 0) {
+    C <- A
+  } else {
+    # bits <- sparse(1, max(ma, mb))
+    bits <- zeros(1, max(c(ma, mb)))
+    bits[A] <- 1
+    bits[B] <- 1
+    C <- find(bits)
+  }
+  return(C)
 }

From 76d31e593ae359465bd64b0d41e79829af766ef2 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 13:39:48 +0100
Subject: [PATCH 15/23] Translated setdiag()

---
 R/setdiag.R | 30 ++++++++++++++----------------
 1 file changed, 14 insertions(+), 16 deletions(-)

diff --git a/R/setdiag.R b/R/setdiag.R
index b77f76d..df57346 100644
--- a/R/setdiag.R
+++ b/R/setdiag.R
@@ -1,21 +1,19 @@
 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)
+  #  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
+  n <- length(M)
+  if (length(v) == 1) {
+    v <- repmat(v, c(1, n))
+  }
 
-  # % e.g., for 3x3 matrix,  elements are numbered
-  # % 1 4 7
-  # % 2 5 8
-  # % 3 6 9
-  # % so diagnoal = [1 5 9]
+  #  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;
+  J <- seq(1, n ^ 2, n + 1)
+  M[J] <- v
+  return(M)
 }

From ac1e73c6fae9c3d820e56623a283fa10b19a37f7 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 14:05:13 +0100
Subject: [PATCH 16/23] Translated cliques_to_jtree()

---
 R/cliques_to_jtree.R | 93 ++++++++++++++++++++------------------------
 1 file changed, 42 insertions(+), 51 deletions(-)

diff --git a/R/cliques_to_jtree.R b/R/cliques_to_jtree.R
index 40c6f3b..68d57a5 100644
--- a/R/cliques_to_jtree.R
+++ b/R/cliques_to_jtree.R
@@ -1,58 +1,49 @@
 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
+  #  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 <- zeros(num_cliques, 1)
+  for (i in 1:num_cliques) {
+    B[i, cliques[[i]]] <- 1
+    w[i] <- prod(ns(cliques[[i]]))
+  }
 
+  #  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 %*% t(B)
+  C1 <- setdiag(C1, 0)
 
-  # 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
+  #  C2[i, j] <- w(i) + w(j)
+  num_cliques <- length(w)
+  W <- repmat(w, c(1, num_cliques))
+  C2 <- W + t(W)
+  C2 <- setdiag(C2, 0)
 
+  jtree <- zeros(minimum_spanning_tree(-C1, C2))# Using - C1 gives * maximum * spanning tree
 
-  # % 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;
+  #  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
+  return(list("jtree" = jtree, "root" = root, "B" = B, "w" = w))
 }

From aff5b3de15a094e988e92d5ffb6bcbcc8cb3681a Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 14:05:23 +0100
Subject: [PATCH 17/23] Added more functions to translate

---
 R/minimum_spanning_tree.R            | 51 ++++++++++++++++++++++++++++
 tests/testthat/test-spatialMixture.R |  2 ++
 2 files changed, 53 insertions(+)
 create mode 100644 R/minimum_spanning_tree.R

diff --git a/R/minimum_spanning_tree.R b/R/minimum_spanning_tree.R
new file mode 100644
index 0000000..a9e7288
--- /dev/null
+++ b/R/minimum_spanning_tree.R
@@ -0,0 +1,51 @@
+minimum_spanning_tree <- function(C1, C2) stop("needs translation")
+# 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/tests/testthat/test-spatialMixture.R b/tests/testthat/test-spatialMixture.R
index 4bf01e6..83b9a9f 100644
--- a/tests/testthat/test-spatialMixture.R
+++ b/tests/testthat/test-spatialMixture.R
@@ -42,5 +42,7 @@ test_that("lakseKlitik() and subfunctions produce expected output", {
   # TODO: test mysubset()
   # TODO: test findCliques()
   # TODO: test cliques_to_jtree()
+  # TODO: test minimum_spanning_tree()
+  # TODO: test myunion()
   # TODO: test lakseKlitik()
 })

From c69c0639aee2d5052568925d361994d918460c39 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 14:17:31 +0100
Subject: [PATCH 18/23] Translated elim_order()

---
 R/elim_order.R | 126 ++++++++++++++++++++++++-------------------------
 1 file changed, 61 insertions(+), 65 deletions(-)

diff --git a/R/elim_order.R b/R/elim_order.R
index 9d9090e..3fdb547 100644
--- a/R/elim_order.R
+++ b/R/elim_order.R
@@ -1,71 +1,67 @@
 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
-  # %
+  #  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.
+  #  Warning: This code is pretty old and could probably be made faster.
 
-  # n = length(G);
-  # %if nargin < 3, stage = { 1:n }; end % no constraints
+  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.
+  #  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
+  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 in 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 in 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
+    }
+    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 (!is.null(ns)) {
+      G[ns, ns] <- 1
+      G <- setdiag(G, 0)
+    }
+    # if (!any(as.logical(uneliminated(stage{t})))# are we allowed to the next slice?) {
+    #   t <- t + 1
+    # }
+  }
+  return(order)
 }

From 6a81eb9d2cca2414d56380298e3d61320a961473 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 14:25:38 +0100
Subject: [PATCH 19/23] Translated argmin()

---
 R/argmin.R | 12 ++++++++++++
 1 file changed, 12 insertions(+)
 create mode 100644 R/argmin.R

diff --git a/R/argmin.R b/R/argmin.R
new file mode 100644
index 0000000..4d50d1c
--- /dev/null
+++ b/R/argmin.R
@@ -0,0 +1,12 @@
+argmin <- function(v) {
+  #  ARGMIN Return as a subscript vector the location of the smallest element of a multidimensional array v.
+  #  indices <- argmin(v)
+  #  Returns the first minimum in the case of ties.
+  #  Example:
+  #  X = [2 8 4 7 3 9]
+  #  argmin[X] <- [1 1], i.e., row 1 column 1
+  m <- i <- NA
+  c(m, i) %<-% matlab2r::min(v)
+  indices <- ind2subv(mysize(v), i)
+  return(indices)
+}

From 6f8e1a51a2dd835d13c3a3f728e6cb74d7c57218 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 14:38:47 +0100
Subject: [PATCH 20/23] Translated mysize()

---
 R/mysize.R | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)
 create mode 100644 R/mysize.R

diff --git a/R/mysize.R b/R/mysize.R
new file mode 100644
index 0000000..a888d67
--- /dev/null
+++ b/R/mysize.R
@@ -0,0 +1,17 @@
+mysize <- function(M) {
+  #  MYSIZE Like the built - in size, except it returns n if (M is a vector of length n, and 1 if M is a scalar.) {
+  #  sz <- mysize(M)
+  #  The behavior is best explained by examples
+  #  - M <- rand(1, 1),   mysize[M] <- 1,      size(M) = [1 1]
+  #  - M <- rand(2, 1),   mysize[M] <- 2,      size(M) = [2 1]
+  #  - M <- rand(1, 2),   mysize[M] <- 2,      size(M) = [1 2]
+  #  - M <- rand(2, 2,1), mysize[M] <- [2 2],  size(M) = [2 2]
+  #  - M <- rand(1, 2,1), mysize[M] <- 2,      size(M) = [1 2]
+
+  if (myisvector(M)) {
+    sz <- length(M)
+  } else {
+    sz <- size(M)
+  }
+  return(sz)
+}

From b6afea85e18e648e1fbcebc15b5c0f11dc9d9dc1 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 14:47:56 +0100
Subject: [PATCH 21/23] Translated myisvector()

---
 R/myisvector.R | 7 +++++++
 1 file changed, 7 insertions(+)
 create mode 100644 R/myisvector.R

diff --git a/R/myisvector.R b/R/myisvector.R
new file mode 100644
index 0000000..498e97b
--- /dev/null
+++ b/R/myisvector.R
@@ -0,0 +1,7 @@
+myisvector <- function(V) {
+  # Kuten isvector(V)
+
+  A <- size(V)
+  r <- (length(A) == 2) & (min(A) == 1)
+  return(r)
+}

From 0048a4e5fdcd8799052632364af7ebf954098265 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 15:04:44 +0100
Subject: [PATCH 22/23] Added placeholder for ind2subv()

---
 R/ind2subv.R | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)
 create mode 100644 R/ind2subv.R

diff --git a/R/ind2subv.R b/R/ind2subv.R
new file mode 100644
index 0000000..6fef574
--- /dev/null
+++ b/R/ind2subv.R
@@ -0,0 +1,38 @@
+ind2subv <- function(siz, ndx) stop("Needs translation")
+# function sub = ind2subv(siz, ndx)
+# % IND2SUBV Like the built-in ind2sub, but returns the answer as a row vector.
+# % sub = ind2subv(siz, ndx)
+# %
+# % siz and ndx can be row or column vectors.
+# % sub will be of size length(ndx) * length(siz).
+# %
+# % Example
+# % ind2subv([2 2 2], 1:8) returns
+# %  [1 1 1
+# %   2 1 1
+# %   ...
+# %   2 2 2]
+# % That is, the leftmost digit toggle fastest.
+# %
+# % See also SUBV2IND
+
+# n = length(siz);
+
+# if n==0
+#   sub = ndx;
+#   return;
+# end
+
+# if all(siz==2)
+#   sub = dec2bitv(ndx-1, n);
+#   sub = sub(:,n:-1:1)+1;
+#   return;
+# end
+
+# cp = [1 cumprod(siz(:)')];
+# ndx = ndx(:) - 1;
+# sub = zeros(length(ndx), n);
+# for i = n:-1:1 % i'th digit
+#   sub(:,i) = floor(ndx/cp(i))+1;
+#   ndx = rem(ndx,cp(i));
+# end

From a3e51c1577e260ec9c6ccae4f16d53b9197f302c Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Thu, 22 Dec 2022 15:04:58 +0100
Subject: [PATCH 23/23] Added more TODOs to unit tests

---
 tests/testthat/test-spatialMixture.R | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/tests/testthat/test-spatialMixture.R b/tests/testthat/test-spatialMixture.R
index 83b9a9f..1f61c4c 100644
--- a/tests/testthat/test-spatialMixture.R
+++ b/tests/testthat/test-spatialMixture.R
@@ -44,5 +44,10 @@ test_that("lakseKlitik() and subfunctions produce expected output", {
   # TODO: test cliques_to_jtree()
   # TODO: test minimum_spanning_tree()
   # TODO: test myunion()
+  # TODO: test argmin()
+  # TODO: test mysize()
+  # TODO: test ind2subv()
+  # TODO: test myisvector()
+  # TODO: ... and anythin left from findCliques.m
   # TODO: test lakseKlitik()
 })