Merge branch 'simulateIndividuals' into dev

NAMESPACE

@@ -14,6 +14,7 @@ export(proportion2str)
 export(rand)
 export(randdir)
 export(repmat)
+export(simulateIndividuals)
 export(simuloiAlleeli)
 export(suoritaMuutos)
 export(times)

R/admix1.R

@@ -453,29 +453,4 @@ admix1 <- function(tietue) {
 #         simuloidut = randdir(counts(1:noalle(j),j,i) , noalle(j));
 #         allfreqs(1:noalle(j),j,i) = simuloidut;
 #     end
-# end
-
-# %--------------------------------------------------------------------------
-
-
-# function refData = simulateIndividuals(n,rowsFromInd,allfreqs,pop, missing_level)
-# % simulate n individuals from population pop, such that approximately
-# % proportion "missing_level" of the alleles are present. 
-
-# nloci = size(allfreqs,2);
-
-# refData = zeros(n*rowsFromInd,nloci);
-# counter = 1;  % which row will be generated next.
-
-# for ind = 1:n
-#     for loc = 1:nloci
-#         for k=0:rowsFromInd-1
-#             if rand<missing_level
-#                 refData(counter+k,loc) = simuloiAlleeli(allfreqs,pop,loc);
-#             else
-#                 refData(counter+k,loc) = -999;
-#             end
-#         end
-#     end
-#     counter = counter+rowsFromInd;
 # end

R/simulateIndividuals.R

@@ -0,0 +1,27 @@
+#' @title Simulate individuals
+#' @description simulate n individuals from population pop, such that
+#' proportion "missing_level" of the alleles are present. 
+#' @export
+
+simulateIndividuals <- function(n, rowsFromInd, allfreqs, pop, missing_level) {
+    nloci <- size(allfreqs, 2)
+
+    refData <- zeros(n * rowsFromInd, nloci)
+    counter <- 1 # which row will be generated next.
+
+    for (ind in 1:n) {
+        for (loc in 1:nloci) {
+            for (k in 0:(rowsFromInd - 1)) {
+                if (runif(1) < missing_level) {
+                    refData[counter + k, loc] <- simuloiAlleeli(
+                        allfreqs, pop, loc
+                    )
+                } else {
+                    refData[counter + k, loc] <- -999
+                }
+            }
+        }
+        counter <- counter + rowsFromInd
+    }
+    return(refData)
+}

R/simuloiAlleeli.R

@@ -4,11 +4,22 @@
 #' @export
 
 simuloiAlleeli <- function(allfreqs, pop, loc) {
-    if (length(dim(allfreqs)) == 3) { # distinguish between arrays and matrices
+    if (length(dim(allfreqs)) == 0) {
-        freqs <- allfreqs[, loc, pop]
+        freqs <- 1
     } else {
-        freqs <- allfreqs[, loc]
+        if (length(dim(allfreqs)) == 3) { # distinguish between array and matrix
+            freqs <- allfreqs[, loc, pop]
+        } else {
+            freqs <- allfreqs[, loc]
+        }
     }
+    # freqs <- ifelse(is.null(length(dim(allfreqs)), allfreqs[loc], 0)
+    # freqs <- switch() + 1,
+    #     allfreqs[, loc],
+    #     allfreqs[, loc, pop]
+    # )
+    
+
     cumsumma <- cumsum(freqs)
     arvo <- runif(1)
     isommat <- which(cumsumma > arvo)

R/size.R

@@ -0,0 +1,30 @@
+#' @title Size of an object
+#' @description This functions tries to replicate the behavior of the base function "size" in Matlab
+#' @param x object to be evaluated
+#' @param d dimension of object to be evaluated
+#' @note On MATLAB, size(1, 100) returns 1. As a matter of fact, if the user 
+#' calls for a dimension which x doesn't have `size()` always returns 1. R's
+#' 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.
+size <- function(x, d) {
+    # Determining the number of dimensions
+    if (length(x) == 1) {
+        # x is surely a scalar
+        return(1)
+    } else {
+        # x is a vector, a matrix or an array
+        n_dim <- ifelse(is.null(dim(x)), 1, length(dim(x)))
+        if (missing(d)) {
+            if (n_dim == 1) {
+                out <- range(x)
+            } else {
+                out <- dim(x)
+            }
+        } else {
+            out <- ifelse(n_dim == 1, range(x)[d], dim(x)[d])
+            if (is.na(out)) out <- 1  # for MATLAB compatibility
+        }
+        return(out)
+    }
+}

man/simulateIndividuals.Rd

@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/simulateIndividuals.R
+\name{simulateIndividuals}
+\alias{simulateIndividuals}
+\title{Simulate individuals}
+\usage{
+simulateIndividuals(n, rowsFromInd, allfreqs, pop, missing_level)
+}
+\description{
+simulate n individuals from population pop, such that
+proportion "missing_level" of the alleles are present.
+}

man/size.Rd

@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/size.R
+\name{size}
+\alias{size}
+\title{Size of an object}
+\usage{
+size(x, d)
+}
+\arguments{
+\item{x}{object to be evaluated}
+
+\item{d}{dimension of object to be evaluated}
+}
+\description{
+This functions tries to replicate the behavior of the base function "size" in Matlab
+}
+\note{
+On MATLAB, size(1, 100) returns 1. As a matter of fact, if the user 
+calls for a dimension which x doesn't have `size()` always returns 1. R's
+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.
+}

tests/testthat/test-admix1.R

@@ -183,10 +183,32 @@ 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))
     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(mx1, 1, 2), 2)
     expect_equal(simuloiAlleeli(mx2, 1, 3), 1)
+})
+
+test_that("simulateIndividuals works like on Matlab", {
+    expect_equal(
+        object = simulateIndividuals(1, 3, 2, 0, .1),
+        expected = matrix(rep(-999, 3), 3)
+    )
+    expect_equal(
+        object = simulateIndividuals(5, 3, 1:3, 4, 0),
+        expected = matrix(rep(-999, 15 * 3), 15)
+    )
+    expect_equal(
+        object = simulateIndividuals(3, 3, 2, 1, 1),
+        expected = matrix(rep(1, 9), 9)
+    )
+    set.seed(2)
+    expect_equal(
+        object = sum(simulateIndividuals(3, 3, 2, 1, .5) == 1),
+        expected = 6
+    )
 })

tests/testthat/test-convertedBaseFunctions.R

@@ -65,4 +65,23 @@ test_that("times works as expected", {
 test_that("colon works as expected (hee hee)", {
     expect_equal(colon(1, 4), 1:4)
     expect_length(colon(4, 1), 0)
+})
+
+test_that("size works as on MATLAB", {
+    sk <- 10
+    vk <- 1:4
+    mx <- matrix(1:6, 2)
+    ra <- array(1:24, c(2, 3, 4))
+    expect_equal(size(sk), 1)
+    expect_equal(size(vk), c(1, 4))
+    expect_equal(size(mx), c(2, 3))
+    expect_equal(size(ra), c(2, 3, 4))
+    expect_equal(size(sk, 199), 1)
+    expect_equal(size(vk, 199), 1)
+    expect_equal(size(mx, 199), 1)
+    expect_equal(size(ra, 199), 1)
+    expect_equal(size(vk, 2), 4)
+    expect_equal(size(mx, 2), 3)
+    expect_equal(size(ra, 2), 3)
+    expect_equal(size(ra, 3), 4)
 })