From 1ac54301a60d71a7fa52bdfe265f0811a2f3bc58 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Thu, 30 Jan 2020 17:19:10 +0100 Subject: [PATCH 1/4] Adde basic elements of simulateAllFreqs --- NAMESPACE | 1 + R/admix1.R | 29 +---------------------------- R/simulateAllFreqs.R | 28 ++++++++++++++++++++++++++++ man/simulateAllFreqs.Rd | 17 +++++++++++++++++ 4 files changed, 47 insertions(+), 28 deletions(-) create mode 100644 R/simulateAllFreqs.R create mode 100644 man/simulateAllFreqs.Rd diff --git a/NAMESPACE b/NAMESPACE index f49f1d9..673ff82 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,6 +14,7 @@ export(proportion2str) export(rand) export(randdir) export(repmat) +export(simulateAllFreqs) export(simulateIndividuals) export(simuloiAlleeli) export(suoritaMuutos) diff --git a/R/admix1.R b/R/admix1.R index 8251e22..28ce731 100644 --- a/R/admix1.R +++ b/R/admix1.R @@ -426,31 +426,4 @@ admix1 <- function(tietue) { # end # prioriAlleelit = repmat(prioriAlleelit, [1,1,npops]); # counts = COUNTS + prioriAlleelit; -# allFreqs = counts./sumCounts; - - -# function allfreqs = simulateAllFreqs(noalle) -# % Lisää jokaista alleelia joka populaation joka lokukseen j 1/noalle(j) -# % verran. Näin saatuja counts:eja vastaavista Dirichlet-jakaumista -# % simuloidaan arvot populaatioiden alleelifrekvensseille. - -# global COUNTS; - -# max_noalle = size(COUNTS,1); -# nloci = size(COUNTS,2); -# npops = size(COUNTS,3); - -# prioriAlleelit = zeros(max_noalle,nloci); -# for j=1:nloci -# prioriAlleelit(1:noalle(j),j) = 1/noalle(j); -# end -# prioriAlleelit = repmat(prioriAlleelit, [1,1,npops]); -# counts = COUNTS + prioriAlleelit; -# allfreqs = zeros(size(counts)); - -# for i=1:npops -# for j=1:nloci -# simuloidut = randdir(counts(1:noalle(j),j,i) , noalle(j)); -# allfreqs(1:noalle(j),j,i) = simuloidut; -# end -# end \ No newline at end of file +# allFreqs = counts./sumCounts; \ No newline at end of file diff --git a/R/simulateAllFreqs.R b/R/simulateAllFreqs.R new file mode 100644 index 0000000..3a747e5 --- /dev/null +++ b/R/simulateAllFreqs.R @@ -0,0 +1,28 @@ +#' @title Simulate All Frequencies +#' @description Lisää jokaista alleelia joka populaation joka lokukseen j1/noalle(j) verran. Näin saatuja counts:eja vastaavista Dirichlet-jakaumista simuloidaan arvot populaatioiden alleelifrekvensseille. +#' Add each allele to each locus in each population by j 1 / noalle(j). The Dirichlet distributions corresponding to the counts thus obtained simulate values for the allele frequencies of the populations. +#' @param noalle noalle +#' @param COUNTS COUNTS +#' @export + +simulateAllFreqs <- function(noalle, COUNTS = matrix(0)) { + max_noalle <- size(COUNTS, 1) + nloci <- size(COUNTS, 2) + npops <- size(COUNTS, 3) + + prioriAlleelit <- zeros(max_noalle, nloci) + for (j in 1:nloci) { + prioriAlleelit[1:noalle[j], j] <- 1 / noalle[j] + } + prioriAlleelit <- repmat(prioriAlleelit, matrix(c(1, 1, npops), ncol = 1)) + counts <- COUNTS + prioriAlleelit + allfreqs <- zeros(size(counts)) + + for (i in 1:npops) { + for (j in 1:nloci) { + simuloidut <- randdir(counts[1:noalle[j], j, i] , noalle[j]) + allfreqs[1:noalle[j], j, i] <- simuloidut + } + } + return(allfreqs) +} \ No newline at end of file diff --git a/man/simulateAllFreqs.Rd b/man/simulateAllFreqs.Rd new file mode 100644 index 0000000..463ef80 --- /dev/null +++ b/man/simulateAllFreqs.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulateAllFreqs.R +\name{simulateAllFreqs} +\alias{simulateAllFreqs} +\title{Simulate All Frequencies} +\usage{ +simulateAllFreqs(noalle, COUNTS = matrix(0)) +} +\arguments{ +\item{noalle}{noalle} + +\item{COUNTS}{COUNTS} +} +\description{ +Lisää jokaista alleelia joka populaation joka lokukseen j1/noalle(j) verran. Näin saatuja counts:eja vastaavista Dirichlet-jakaumista simuloidaan arvot populaatioiden alleelifrekvensseille. +Add each allele to each locus in each population by j 1 / noalle(j). The Dirichlet distributions corresponding to the counts thus obtained simulate values for the allele frequencies of the populations. +} From c9c9b7cb261dbe16780f804b5f6849fd898fb554 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Thu, 30 Jan 2020 17:31:25 +0100 Subject: [PATCH 2/4] Fixed handling of empty x argument --- NAMESPACE | 1 + R/size.R | 8 ++++++++ 2 files changed, 9 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index 673ff82..927db13 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,6 +17,7 @@ export(repmat) export(simulateAllFreqs) export(simulateIndividuals) export(simuloiAlleeli) +export(size) export(suoritaMuutos) export(times) importFrom(stats,runif) diff --git a/R/size.R b/R/size.R index ce245f7..001817b 100644 --- a/R/size.R +++ b/R/size.R @@ -7,8 +7,16 @@ #' default behavior is more reasonable in those cases (i.e., returning NA), #' but since the point of this function is to replicate MATLAB behaviors #' (bugs and questionable behaviors included), this function also does this. +#' @export size <- function(x, d) { # Determining the number of dimensions + if (all(is.na(x))) { + if (missing(d)) { + return(c(0, 0)) + } else { + return(ifelse(d <= 2, 0, 1)) + } + } if (length(x) == 1) { # x is surely a scalar return(1) From f6768c67bcfdb8fd07252c155f0177e2cb739bbc Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Thu, 30 Jan 2020 17:52:20 +0100 Subject: [PATCH 3/4] Adjusted behavior to benchmark --- R/simulateAllFreqs.R | 18 +++++++++++------- man/simulateAllFreqs.Rd | 2 +- tests/testthat/test-admix1.R | 10 ++++++++++ 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/R/simulateAllFreqs.R b/R/simulateAllFreqs.R index 3a747e5..81995b4 100644 --- a/R/simulateAllFreqs.R +++ b/R/simulateAllFreqs.R @@ -5,23 +5,27 @@ #' @param COUNTS COUNTS #' @export -simulateAllFreqs <- function(noalle, COUNTS = matrix(0)) { +simulateAllFreqs <- function(noalle, COUNTS = matrix(NA, 0, 0)) { max_noalle <- size(COUNTS, 1) nloci <- size(COUNTS, 2) npops <- size(COUNTS, 3) prioriAlleelit <- zeros(max_noalle, nloci) - for (j in 1:nloci) { - prioriAlleelit[1:noalle[j], j] <- 1 / noalle[j] + if (nloci > 0) { + for (j in 1:nloci) { + prioriAlleelit[1:noalle[j], j] <- 1 / noalle[j] + } } - prioriAlleelit <- repmat(prioriAlleelit, matrix(c(1, 1, npops), ncol = 1)) + prioriAlleelit <- repmat(prioriAlleelit, matrix(c(1, 1, npops), 1)) counts <- COUNTS + prioriAlleelit allfreqs <- zeros(size(counts)) for (i in 1:npops) { - for (j in 1:nloci) { - simuloidut <- randdir(counts[1:noalle[j], j, i] , noalle[j]) - allfreqs[1:noalle[j], j, i] <- simuloidut + if (nloci > 0) { + for (j in 1:nloci) { + simuloidut <- randdir(counts[1:noalle[j], j, i] , noalle[j]) + allfreqs[1:noalle[j], j, i] <- simuloidut + } } } return(allfreqs) diff --git a/man/simulateAllFreqs.Rd b/man/simulateAllFreqs.Rd index 463ef80..97bf0e6 100644 --- a/man/simulateAllFreqs.Rd +++ b/man/simulateAllFreqs.Rd @@ -4,7 +4,7 @@ \alias{simulateAllFreqs} \title{Simulate All Frequencies} \usage{ -simulateAllFreqs(noalle, COUNTS = matrix(0)) +simulateAllFreqs(noalle, COUNTS = matrix()) } \arguments{ \item{noalle}{noalle} diff --git a/tests/testthat/test-admix1.R b/tests/testthat/test-admix1.R index 04e004a..dbd1cda 100644 --- a/tests/testthat/test-admix1.R +++ b/tests/testthat/test-admix1.R @@ -211,4 +211,14 @@ test_that("simulateIndividuals works like on Matlab", { object = sum(simulateIndividuals(3, 3, 2, 1, .5) == 1), expected = 6 ) +}) + +test_that("simulateAllFreqs works as expected", { + empty_mt <- matrix(NA, 0, 0) + expect_equivalent(suppressWarnings(simulateAllFreqs(3)), empty_mt) + expect_equivalent(suppressWarnings(simulateAllFreqs(3:5)), empty_mt) + expect_equivalent( + object = suppressWarnings(simulateAllFreqs(matrix(1:4, 2))), + expected = empty_mt + ) }) \ No newline at end of file From 8cab8e80c6f7929e3bb93dc76be7e87d071907f1 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Thu, 30 Jan 2020 17:52:37 +0100 Subject: [PATCH 4/4] Fixed random test failures --- tests/testthat/test-admix1.R | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/tests/testthat/test-admix1.R b/tests/testthat/test-admix1.R index dbd1cda..b5cd9e3 100644 --- a/tests/testthat/test-admix1.R +++ b/tests/testthat/test-admix1.R @@ -182,21 +182,23 @@ test_that("computePersonalAllFreqs works like on Matlab", { }) test_that("simuloiAlleeli works like on Matlab", { - # TODO: test on vector - sk1 <- 2 - ra1 <- array(1:12, c(2, 2, 3)) + sk <- 2 + vk <- 1:3 + ra <- array(1:12, c(2, 2, 3)) mx1 <- matrix(c(3, 5, 0, 9), 2) mx2 <- matrix(c(3, 5, 0, 9, 5, 8), 2) - expect_equal(simuloiAlleeli(sk1, 1, 1), 1) - expect_equal(simuloiAlleeli(ra1, 2, 1), 1) + expect_equal(simuloiAlleeli(sk, 1, 1), 1) + expect_equal(simuloiAlleeli(vk, 1, 2), 1) + expect_equal(simuloiAlleeli(ra, 2, 1), 1) expect_equal(simuloiAlleeli(mx1, 1, 2), 2) expect_equal(simuloiAlleeli(mx2, 1, 3), 1) }) test_that("simulateIndividuals works like on Matlab", { + set.seed(2) expect_equal( - object = simulateIndividuals(1, 3, 2, 0, .1), - expected = matrix(rep(-999, 3), 3) + object = simulateIndividuals(1, 3, 2, 0, .2), + expected = matrix(c(1, -999, 1), ncol = 1) ) expect_equal( object = simulateIndividuals(5, 3, 1:3, 4, 0),