From c9c310d1fa41d11083555bdbeee3ce5d5a15012d Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <w.l.netto@medisin.uio.no>
Date: Tue, 8 Feb 2022 14:01:02 +0100
Subject: [PATCH] Translated getPopDistancesByKL (#2)

---
 R/getPopDistancesByKL.R            | 32 ++++++++++++++++++++++++++++++
 R/globals.R                        |  3 ++-
 tests/testthat/test-greedyPopMix.R |  7 +++++++
 3 files changed, 41 insertions(+), 1 deletion(-)
 create mode 100644 R/getPopDistancesByKL.R

diff --git a/R/getPopDistancesByKL.R b/R/getPopDistancesByKL.R
new file mode 100644
index 0000000..83c7c5f
--- /dev/null
+++ b/R/getPopDistancesByKL.R
@@ -0,0 +1,32 @@
+getPopDistancesByKL <- function(adjprior) {
+  # Laskee populaatioille etהisyydet
+  # kהyttהen KL-divergenssi?
+  COUNTS <- COUNTS[seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), ]
+  maxnoalle <- size(COUNTS, 1)
+  nloci <- size(COUNTS, 2)
+  npops <- size(COUNTS, 3)
+  distances <- zeros(choose(npops, 2), 1)
+
+  d <- zeros(maxnoalle, nloci, npops)
+  prior <- adjprior
+  prior[find(prior == 1)] <- 0
+  nollia <- find(all(prior == 0)) # Lokukset, joissa oli havaittu vain yht?alleelia.
+  prior[1, nollia] <- 1
+  for (pop1 in 1:npops) {
+    d[, , pop1] <- (squeeze(COUNTS[, , pop1]) + prior) / repmat(sum(squeeze(COUNTS[, , pop1]) + prior), c(maxnoalle, ncol(prior)))
+  }
+  pointer <- 1
+  for (pop1 in 1:(npops - 1)) {
+    for (pop2 in (pop1 + 1):npops) {
+      dist1 <- d[, , pop1]
+      dist2 <- d[, , pop2]
+      div12 <- sum(sum(dist1 * log2((dist1 + 10^-10) / (dist2 + 10^-10)))) / nloci
+      div21 <- sum(sum(dist2 * log2((dist2 + 10^-10) / (dist1 + 10^-10)))) / nloci
+      div <- (div12 + div21) / 2
+      distances[pointer] <- div
+      pointer <- pointer + 1
+    }
+  }
+  Z <- linkage(t(distances))
+  return(list(Z = Z, distances = distances))
+}
diff --git a/R/globals.R b/R/globals.R
index d96ba00..85a3580 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -3,7 +3,8 @@ SUMCOUNTS <- array(0, dim = c(100, 100))
 PARTITION <- array(1, dim = 100)
 POP_LOGML <- array(1, dim = 100)
 LOGDIFF <- array(1, dim = c(100, 100))
-# If handling globas break, try other ideas from https://stackoverflow.com/a/65252740/1169233
+# If handling globas break, try other ideas from
+# https://stackoverflow.com/a/65252740/1169233
 
 utils::globalVariables(
   c("PARTITION", "COUNTS", "SUMCOUNTS", "LOGDIFF", "POP_LOGML", "GAMMA_LN")
diff --git a/tests/testthat/test-greedyPopMix.R b/tests/testthat/test-greedyPopMix.R
index a770d12..c54b2a1 100644
--- a/tests/testthat/test-greedyPopMix.R
+++ b/tests/testthat/test-greedyPopMix.R
@@ -8,4 +8,11 @@ test_that("Auxiliary functions work properly", {
     rowsFromInd = 2
   )
   expect_equal(findOutRowsFromInd(x, y, "Diploid"), z)
+  expect_equal(
+    getPopDistancesByKL(x2),
+    list(
+      Z = matrix(c(c(1, 101:198), c(2:100), rep(0, 99)), nrow = 99, ncol = 3),
+      distances = as.matrix(rep(0, 4950))
+    )
+  )
 })