From 6534d4fa4add77e70a58ed4f7f0534ed6971e046 Mon Sep 17 00:00:00 2001
From: Waldir Leoncio <waldir.leoncio@gmail.com>
Date: Wed, 15 Jan 2020 17:05:36 +0100
Subject: [PATCH] Finished proper implementation of etsiParas

---
 R/etsiParas.R                | 24 ++++++++++++++++++++----
 tests/testthat/test-admix1.R | 27 ++++++++++++++++++++++++++-
 2 files changed, 46 insertions(+), 5 deletions(-)

diff --git a/R/etsiParas.R b/R/etsiParas.R
index 22eda5e..98606e7 100644
--- a/R/etsiParas.R
+++ b/R/etsiParas.R
@@ -1,9 +1,25 @@
-etsiParas <- function = (osuus, osuusTaulu, omaFreqs, logml) {
-    ready <- 0;
+#' @export
+#' @title Etsi Paras
+#' @description Search for the best?
+#' @param osuus Percentages?
+#' @param omaFreqs own Freqs?
+#' @param osuusTaulu Percentage table?
+#' @param logml log maximum likelihood
+etsiParas <- function (osuus, osuusTaulu, omaFreqs, logml) {
+    ready <- 0
     while (ready != 1) {
         muutokset <- laskeMuutokset4(osuus, osuusTaulu, omaFreqs, logml)
-        [maxMuutos, indeksi] = max(muutokset[1:end]) # TODO: how does this work on Matlab?
-        if (maxMuutos > 0) {
+
+        # Work around R's max() limitation on complex numbers
+        if (any(sapply(muutokset, class) == "complex")) {
+            maxRe <- max(Re(as.vector(muutokset)))
+            maxIm <- max(Im(as.vector(muutokset)))
+            maxMuutos <- complex(real = maxRe, imaginary = maxIm)
+        } else {
+            maxMuutos <- max(as.vector(muutokset))
+        }
+        indeksi <- which(muutokset == maxMuutos)
+        if (Re(maxMuutos) > 0) {
             osuusTaulu <- suoritaMuutos(osuusTaulu, osuus, indeksi)
             logml <- logml + maxMuutos
         } else {
diff --git a/tests/testthat/test-admix1.R b/tests/testthat/test-admix1.R
index 876ec49..5221db0 100644
--- a/tests/testthat/test-admix1.R
+++ b/tests/testthat/test-admix1.R
@@ -1,6 +1,5 @@
 context("Admixture analysis")
 
-
 test_that("learn*partition behaves like on Matlab", {
     # Test data
     p1 <- c(0, .5, 1, 1.5)
@@ -146,4 +145,30 @@ test_that("suoritaMuutos works like on Matlab", {
     expect_equal(suoritaMuutos(mx2, 0, 5), mx2)
     expect_equal(suoritaMuutos(mx2, 0, 5), mx2)
     expect_equal(suoritaMuutos(mx2, -3, 6), matrix(c(13, 9, 5, 8, 8, -10), 2))
+})
+
+test_that("laskeMuutokset4 works like on Matlab", {
+    # TODO: build these tests based on problems found in etsiParas
+    mx1 <- t(c(.4, 7))
+    expect_equivalent(
+        object = laskeMuutokset4(2, mx1, c(8, 2), 3),
+        expected = t(c(0, .3742)),
+        tol = .0001
+    )
+})
+
+test_that("etsiParas works like on Matlab", {
+    mx1 <- t(c(.4, 7))
+    expect_equal(etsiParas(2, mx1, c(8, 1), 8), c(.4, 7, 8))
+    expect_equivalent(etsiParas(2, mx1, c(8, 1), 1), c(-1.6, 9, 3.1864), .0001)
+    expect_equivalent(
+        object = etsiParas(5, mx1, c(8, 1), -pi),
+        expected = c(-4.6, 12, 3.8111),
+        tol = .001
+    )
+    expect_equivalent(
+        object = etsiParas(-.5, mx1, c(-1, 0), -10),
+        expected = c(7.4, 0, complex(real = 1.8563, imaginary = 3.1416)),
+        tol = .0001
+    )
 })
\ No newline at end of file