Finished proper implementation of etsiParas

R/etsiParas.R

@@ -1,9 +1,25 @@
-etsiParas <- function = (osuus, osuusTaulu, omaFreqs, logml) {
+#' @export
-    ready <- 0;
+#' @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 {

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