Merge branch 'master' into develop

.github/workflows/R-CMD-check.yaml

@@ -6,10 +6,12 @@ on:
       - main
       - master
       - modules
+      - develop
   pull_request:
     branches:
       - main
       - master
+      - develop
 
 name: R-CMD-check
 

.github/workflows/test-coverage.yaml

@@ -3,10 +3,12 @@ on:
     branches:
       - main
       - master
+      - develop
   pull_request:
     branches:
       - main
       - master
+      - develop
 
 name: test-coverage
 

R/simulate_admixture.R

@@ -105,7 +105,6 @@ simulate_admixture <- function(module = ancestry_module(),
       message("found positive migration rate, assuming two connected populations")
 
     if (module$type == "ancestry") {
-
       if (!inherits(module$input_population, "genomadmixr_simulation")) {
         if (is.list(module$input_population)) {
           input_population2 <- list()

breakme.R

@@ -0,0 +1,34 @@
+num_markers <- 100
+num_indiv <- 100
+chosen_markers <- 1:num_markers
+
+fake_input_data1 <- list()
+fake_input_data1$genomes <- matrix(data = 1,
+                                   nrow = num_indiv,
+                                   ncol = num_markers)
+
+
+fake_input_data1$markers <- chosen_markers
+
+fake_input_data2 <- list()
+fake_input_data2$genomes <- matrix(data = 2,
+                                   nrow = num_indiv,
+                                   ncol = num_markers)
+fake_input_data2$markers <- chosen_markers
+
+class(fake_input_data1) <- "genomeadmixr_data"
+class(fake_input_data2) <- "genomeadmixr_data"
+
+select_matrix <- matrix(ncol = 5, nrow = 1)
+s <- 1.0
+select_matrix[1, ] <- c(50, 1.0, 1 + 0.5 * s, 1 + s, 1)
+
+simul_two_pop <- simulate_admixture(
+  module = sequence_module(molecular_data = list(fake_input_data1,
+                                                 fake_input_data2),
+                           markers = chosen_markers),
+  migration = migration_settings(
+    population_size = c(100, 100),
+    migration_rate = 0.01),
+  select_matrix = select_matrix,
+  total_runtime = 100)

src/Fish_emp.h

@@ -17,31 +17,12 @@ struct Fish_emp {
     chromosome1(c1), chromosome2(c2) {
   }
 
-  Fish_emp(Fish_emp&& other) {
-    chromosome1 = other.chromosome1;
-    chromosome2 = other.chromosome2;
-  }
+  Fish_emp(Fish_emp&& other) = default;
+  Fish_emp& operator=(Fish_emp&& other) = default;
 
-  Fish_emp& operator=(Fish_emp&& other) {
-    if (this != &other) {
-      chromosome1 = other.chromosome1;
-      chromosome2 = other.chromosome2;
-    }
-    return *this;
-  }
+  Fish_emp(const Fish_emp& other) = default;
 
-  Fish_emp(const Fish_emp& other) {
-    chromosome1 = other.chromosome1;
-    chromosome2 = other.chromosome2;
-  }
-
-  Fish_emp& operator=(const Fish_emp& other) {
-    if (this != &other) {
-      chromosome1 = other.chromosome1;
-      chromosome2 = other.chromosome2;
-    }
-    return *this;
-  }
+  Fish_emp& operator=(const Fish_emp& other) = default;
 
   std::vector< int > gamete(double morgan,
                             rnd_t& rndgen,

src/Makevars

@@ -1,2 +1,2 @@
-CXX_STD = CXX14
+CXX_STD = CXX17
 PKG_LIBS = `${R_HOME}/bin/Rscript -e "RcppParallel::RcppParallelLibs()"`

src/RcppExports.cpp

@@ -146,15 +146,15 @@ BEGIN_RCPP
 END_RCPP
 }
 // simulate_migration_emp_cpp
-List simulate_migration_emp_cpp(const NumericMatrix& input_population_1, const NumericMatrix& input_population_2, const NumericVector& marker_positions_R, NumericMatrix select, const NumericVector& pop_sizes, int total_runtime, double morgan, bool verbose, bool track_frequency, const NumericVector& track_markers_R, bool multiplicative_selection, double migration_rate, double mutation_rate, const NumericMatrix& substitution_matrix_R, int num_threads, const NumericVector& recombination_map);
+List simulate_migration_emp_cpp(const NumericMatrix& input_population_1, const NumericMatrix& input_population_2, const NumericVector& marker_positions_R, const NumericMatrix& select, const NumericVector& pop_sizes, int total_runtime, double morgan, bool verbose, bool track_frequency, const NumericVector& track_markers_R, bool multiplicative_selection, double migration_rate, double mutation_rate, const NumericMatrix& substitution_matrix_R, int num_threads, const NumericVector& recombination_map);
 RcppExport SEXP _GenomeAdmixR_simulate_migration_emp_cpp(SEXP input_population_1SEXP, SEXP input_population_2SEXP, SEXP marker_positions_RSEXP, SEXP selectSEXP, SEXP pop_sizesSEXP, SEXP total_runtimeSEXP, SEXP morganSEXP, SEXP verboseSEXP, SEXP track_frequencySEXP, SEXP track_markers_RSEXP, SEXP multiplicative_selectionSEXP, SEXP migration_rateSEXP, SEXP mutation_rateSEXP, SEXP substitution_matrix_RSEXP, SEXP num_threadsSEXP, SEXP recombination_mapSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const NumericMatrix& >::type input_population_1(input_population_1SEXP);
     Rcpp::traits::input_parameter< const NumericMatrix& >::type input_population_2(input_population_2SEXP);
     Rcpp::traits::input_parameter< const NumericVector& >::type marker_positions_R(marker_positions_RSEXP);
-    Rcpp::traits::input_parameter< NumericMatrix >::type select(selectSEXP);
+    Rcpp::traits::input_parameter< const NumericMatrix& >::type select(selectSEXP);
     Rcpp::traits::input_parameter< const NumericVector& >::type pop_sizes(pop_sizesSEXP);
     Rcpp::traits::input_parameter< int >::type total_runtime(total_runtimeSEXP);
     Rcpp::traits::input_parameter< double >::type morgan(morganSEXP);

src/random_functions.h

@@ -6,8 +6,7 @@
 //
 //
 
-#ifndef random_functions_hpp
-#define random_functions_hpp
+#pragma once
 
 #include <random>
 #include <vector>
@@ -19,35 +18,36 @@
 
 
 struct rnd_t {
-  std::mt19937 rndgen_;
+  std::mt19937_64 rndgen_;
   std::uniform_real_distribution<> unif_dist = std::uniform_real_distribution<>(0, 1.0);
   std::uniform_int_distribution<> rand_num_dist;
 
   rnd_t() {
     auto seed = get_seed();
   //  std::cerr << "initializing rnd_t with: " << seed << "\n" << std::flush;
-    rndgen_ = std::mt19937(seed);
+    rndgen_ = std::mt19937_64(seed);
   }
 
   rnd_t(unsigned int seed) {
     auto local_seed = get_seed() + seed;
-    rndgen_ = std::mt19937(local_seed);
+    rndgen_ = std::mt19937_64(local_seed);
   }
 
   void set_seed(unsigned int s) {
-    rndgen_ = std::mt19937(s);
+    rndgen_ = std::mt19937_64(s);
   }
 
   double uniform() {
     return unif_dist(rndgen_);
   }
 
-  int random_number(unsigned int n) {
-    return std::uniform_int_distribution<> (0, n-1)(rndgen_);
+  int random_number(int n) {
+    if (n < 1) return 0;
+    return std::uniform_int_distribution<int> (0, n - 1)(rndgen_);
   }
 
   size_t poisson(double lambda) {
-    return std::poisson_distribution<int>(lambda)(rndgen_);
+    return std::poisson_distribution<size_t>(lambda)(rndgen_);
   }
 
   int get_seed() {
@@ -79,8 +79,13 @@ struct emp_genome {
 
   template <typename T>
   emp_genome(const std::vector<T>& positions) {
-    double total_sum = std::accumulate(positions.begin(),
-                                positions.end(), 0.0);
+    if (positions.empty()) {
+      throw std::runtime_error("positions is empty");
+    }
+    total_sum = std::accumulate(positions.begin(),
+                                positions.end(),
+                                0.0);
+
     double s = 0.0;
     double mult = 1.0 / total_sum;
     cdf_.resize(positions.size());
@@ -92,6 +97,9 @@ struct emp_genome {
   }
 
   size_t index_from_cdf(double p) const {
+
+    if (total_sum <= 0.0) return static_cast<size_t>(p * cdf_.size());
+
     // find index belonging to p
     return static_cast<size_t>(std::distance(cdf_.begin(),
                                              std::lower_bound(cdf_.begin(),
@@ -105,14 +113,13 @@ struct emp_genome {
     std::vector< size_t > indices;
     for(size_t i = 0; i < num_break_points; ++i) {
       auto found_index = index_from_cdf(rndgen.uniform());
-      if (found_index > 0) {
+      if (found_index > 0) { // first position can not be recombination point
         indices.push_back(found_index);
       }
     }
     std::sort(indices.begin(), indices.end());
     indices.push_back(cdf_.size());
     return indices;
   }
+  double total_sum;
 };
-
-#endif /* random_functions_hpp */

tests/testthat/test-LD.R

@@ -2,6 +2,7 @@ context("LD stats")
 
 test_that("calculate_average_LD", {
   testthat::skip_on_os("solaris")
+  cat("test_LD")
   pop_size <- 100
   number_of_founders <- 2
   run_time <- 1000
@@ -49,6 +50,7 @@ test_that("calculate_average_LD", {
 
 test_that("calculate_LD_matrix", {
   testthat::skip_on_os("solaris")
+  cat("test_LD_matrix")
   pop_size <- 100
   number_of_founders <- 2
   sampled_individuals <- pop_size

tests/testthat/test-allele_frequencies.R

@@ -3,6 +3,7 @@ context("allele_frequencies")
 
 test_that("calculate_allele_frequencies", {
   testthat::skip_on_os("solaris")
+  cat("test_allele_frequencies")
   pop_size <- 100
   number_of_founders <- 2
   run_time <- 5

tests/testthat/test-combined_input_data.R

@@ -2,6 +2,8 @@ context("test input data")
 
 test_that("input data", {
   testthat::skip_on_os("solaris")
+  cat("test_input data")
+
 
   chosen_markers <- 1:100
 
@@ -77,6 +79,8 @@ test_that("input data", {
 
 test_that("input data simulation", {
   testthat::skip_on_os("solaris")
+  cat("test_input data2")
+
 
   vx <- simulate_admixture(pop_size = 100,
                            total_runtime = 100)

tests/testthat/test-fst.R

@@ -2,6 +2,8 @@ context("fst")
 
 testthat::test_that("fst", {
   testthat::skip_on_os("solaris")
+  cat("test_fst")
+
   pop_size <- 100
   number_of_founders <- 20
   run_time <- 1

tests/testthat/test-general_usage.R

@@ -2,7 +2,7 @@ context("general usage")
 
 
 test_that("general usage", {
-
+  cat("test_general_usage")
   data("dgrp2.3R.5k.data")
   mks <- sample(dgrp2.3R.5k.data$markers, size = 100,
                 replace = FALSE, prob = NULL)
@@ -109,6 +109,7 @@ test_that("general usage", {
 })
 
 test_that("isofemale usage", {
+  cat("test_isofemale_usage")
   data("dgrp2.3R.5k.data")
 
   mks = sample(dgrp2.3R.5k.data$markers, size = 300, replace = FALSE, prob = NULL)

tests/testthat/test-isoFemales.R

@@ -2,6 +2,7 @@ context("isoFemale creation")
 
 test_that("create_isofemale", {
   testthat::skip_on_os("solaris")
+  cat("create_isofemale")
   pop_size <- 100
   number_of_founders <- 2
   run_time <- 100
@@ -27,6 +28,7 @@ test_that("create_isofemale", {
 
 test_that("create_population_from_isofemales", {
   testthat::skip_on_os("solaris")
+  cat("create isofemalepop")
   pop_size <- 100
   number_of_founders <- 10
   run_time <- 100
@@ -93,6 +95,7 @@ testthat::expect_silent(
 })
 
 test_that("cpp classes", {
+  cat("cpp classes")
   testthat::skip_on_os("solaris")
   a <- matrix(c(0.1, 1, 2, 2), nrow = 2)
   b <- matrix(c(0, 1, 1, -1), nrow = 2)
@@ -133,7 +136,7 @@ test_that("cpp classes", {
 
 test_that("create_isofemale_data", {
   testthat::skip_on_os("solaris")
-
+  cat("test create isofemale data")
  data("dgrp2.3R.5k.data")
 
  females <- create_iso_female(module = sequence_module(

tests/testthat/test-joyplot.R

@@ -1,5 +1,6 @@
 test_that("joyplot", {
   testthat::skip_on_os("solaris")
+  cat("test_joyplot")
   markers <- seq(from = 0.2, to = 0.3, length.out = 100)
 
   selected_pop <- simulate_admixture(module = ancestry_module(number_of_founders = 3,

tests/testthat/test-junctions.R

@@ -1,6 +1,7 @@
 
 test_that("expected_number_junctions", {
   testthat::skip_on_os("solaris")
+  cat("test_junctions")
 
   if (requireNamespace("junctions")) {
 

tests/testthat/test-plink.R

@@ -2,6 +2,7 @@ context("test plink data")
 
 test_that("plink data", {
   testthat::skip_on_os("solaris")
+  cat("test_plink_data")
 
   chosen_markers <- 1:100
 

tests/testthat/test-save_load.R

@@ -2,6 +2,7 @@ context("create_populations")
 
 test_that("save_population", {
   testthat::skip_on_os("solaris")
+  message("test_save_pop")
   pop_size <- 100
   number_of_founders <- 10
   run_time <- 10
@@ -30,6 +31,7 @@ test_that("save_population", {
 })
 
 test_that("data", {
+  message("test_data3")
   data("dgrp2.3R.5k.data")
   testthat::expect_equal(length(dgrp2.3R.5k.data$markers),
                          4603)