From 3c02d89837f9ae89dd1f1cab0c230aff04787021 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Fri, 4 Nov 2022 13:03:32 +0100 Subject: [PATCH 01/15] fix memory access violation --- src/simulate.cpp | 2 +- src/simulate_migration.cpp | 2 +- src/simulate_migration_emp.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/simulate.cpp b/src/simulate.cpp index 10419c4..88b8cd5 100755 --- a/src/simulate.cpp +++ b/src/simulate.cpp @@ -137,7 +137,7 @@ std::vector< Fish > simulate_Population(const std::vector< Fish>& sourcePop, int num_threads) { bool use_selection = false; - if(select(1, 1) >= 0) use_selection = true; + if(select(0, 0) >= 0) use_selection = true; std::vector Pop = sourcePop; std::vector fitness; diff --git a/src/simulate_migration.cpp b/src/simulate_migration.cpp index 777b858..42f93a0 100755 --- a/src/simulate_migration.cpp +++ b/src/simulate_migration.cpp @@ -211,7 +211,7 @@ std::vector< std::vector< Fish > > simulate_two_populations( int num_threads, rnd_t& rndgen) { bool use_selection = false; - if (select(1, 1) >= 0) use_selection = true; + if (select(0, 0) >= 0) use_selection = true; std::vector pop_1 = source_pop_1; std::vector pop_2 = source_pop_2; diff --git a/src/simulate_migration_emp.cpp b/src/simulate_migration_emp.cpp index 23e7de1..0921fa7 100755 --- a/src/simulate_migration_emp.cpp +++ b/src/simulate_migration_emp.cpp @@ -206,7 +206,7 @@ std::vector< std::vector< Fish_emp > > simulate_two_populations( int num_threads) { bool use_selection = false; - if (select(1, 1) >= 0) use_selection = true; + if (select(0, 0) >= 0) use_selection = true; std::vector pop_1 = source_pop_1; std::vector pop_2 = source_pop_2; From 2c5fb2d7b208914166c59df6b9740fd8113de06c Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 8 Nov 2022 10:47:06 +0100 Subject: [PATCH 02/15] add check in draw parent --- R/simulate_admixture.R | 1 - src/simulate_migration_emp.cpp | 10 ++++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/R/simulate_admixture.R b/R/simulate_admixture.R index 059cac3..c4df1d8 100644 --- a/R/simulate_admixture.R +++ b/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() diff --git a/src/simulate_migration_emp.cpp b/src/simulate_migration_emp.cpp index 0921fa7..da3fb41 100755 --- a/src/simulate_migration_emp.cpp +++ b/src/simulate_migration_emp.cpp @@ -40,7 +40,10 @@ Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, if(use_selection) { index = draw_prop_fitness(fitness_migr, max_fitness_migr, rndgen); } else { - index = rndgen.random_number( (int)pop_2.size() ); + index = rndgen.random_number( static_cast(pop_2.size() )); + } + if (index < 0 || index > static_cast(pop_2.size())) { + Rcpp::stop("index out of range in draw parent"); } parent = pop_2[index]; index = index + pop_1.size(); @@ -49,7 +52,10 @@ Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, if(use_selection) { index = draw_prop_fitness(fitness_source, max_fitness_source, rndgen); } else { - index = rndgen.random_number( (int)pop_1.size() ); + index = rndgen.random_number( static_cast(pop_1.size()) ); + } + if (index < 0 || index > static_cast(pop_1.size())) { + Rcpp::stop("index out of range in draw parent"); } parent = pop_1[index]; } From 93edcc2bdfe443d0ec011be4f44f521bf2b31ce2 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 8 Nov 2022 10:59:02 +0100 Subject: [PATCH 03/15] add exception find_location --- src/helper_functions.cpp | 11 +++++++---- src/simulate_migration_emp.cpp | 6 +++--- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/helper_functions.cpp b/src/helper_functions.cpp index 0eada9e..41e2a3d 100755 --- a/src/helper_functions.cpp +++ b/src/helper_functions.cpp @@ -237,6 +237,11 @@ int draw_prop_fitness(const std::vector& fitness, return rndgen.random_number(fitness.size()); } + if (std::isnan(maxFitness)) { + return rndgen.random_number(fitness.size()); + } + + size_t fitness_size = fitness.size(); double inv_fitness = 1.0 / maxFitness; while (true) { @@ -651,13 +656,13 @@ double calculate_fitness(const Fish_emp& focal, auto focal_pos = select(i, 0); auto focal_anc = select(i, 4); if (focal_anc == -1) continue; // do not take into account + int focal_index = find_location(locations, focal_pos); auto a1 = focal.chromosome1[focal_index]; auto a2 = focal.chromosome2[focal_index]; int fit_index = 1 + (a1 == focal_anc) + (a2 == focal_anc); fitness_vec[i] = select(i, fit_index); - // Rcout << a1 << " " << a2 << " " << select(i, fit_index) << "\n"; } if (!multiplicative_selection) { @@ -666,9 +671,6 @@ double calculate_fitness(const Fish_emp& focal, return std::accumulate(fitness_vec.begin(), fitness_vec.end(), 1.0, std::multiplies<>()); - - - } List convert_to_list(const std::vector& v, @@ -749,6 +751,7 @@ int find_location(const std::vector& markers, return std::distance(markers.begin(), loc); } } + Rcpp::stop("could not find location"); return -1; } diff --git a/src/simulate_migration_emp.cpp b/src/simulate_migration_emp.cpp index da3fb41..24c9792 100755 --- a/src/simulate_migration_emp.cpp +++ b/src/simulate_migration_emp.cpp @@ -46,10 +46,10 @@ Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, Rcpp::stop("index out of range in draw parent"); } parent = pop_2[index]; - index = index + pop_1.size(); - // to ensure different indices for pop_1 and pop_2 + index = index + pop_1.size(); // to ensure different indices for pop_1 and pop_2 + } else { - if(use_selection) { + if (use_selection) { index = draw_prop_fitness(fitness_source, max_fitness_source, rndgen); } else { index = rndgen.random_number( static_cast(pop_1.size()) ); From 85443c6cfbec992a22627ee7beedede745358498 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 8 Nov 2022 11:12:02 +0100 Subject: [PATCH 04/15] improve local code --- src/helper_functions.cpp | 209 ++++++++++++++-------------- src/simulate_migration_emp.cpp | 240 +++++++++++++++++---------------- 2 files changed, 229 insertions(+), 220 deletions(-) diff --git a/src/helper_functions.cpp b/src/helper_functions.cpp index 0eada9e..fe091b9 100755 --- a/src/helper_functions.cpp +++ b/src/helper_functions.cpp @@ -15,7 +15,7 @@ #include void force_output() { -// std::this_thread::sleep_for(std::chrono::nanoseconds(100)); + // std::this_thread::sleep_for(std::chrono::nanoseconds(100)); std::this_thread::sleep_for(std::chrono::milliseconds(300)); R_FlushConsole(); R_ProcessEvents(); @@ -124,8 +124,8 @@ arma::mat update_frequency_tibble(const std::vector< Fish >& v, allele_matrix(i, 3) = 0; } - // for (auto it = v.begin(); it != v.end(); ++it) { - for (int i = 0; i < v.size(); ++i) { + // for (auto it = v.begin(); it != v.end(); ++it) { + for (int i = 0; i < v.size(); ++i) { update_anc_chrom(v[i].chromosome1, founder_labels, m, allele_matrix); update_anc_chrom(v[i].chromosome2, founder_labels, @@ -149,12 +149,12 @@ arma::mat update_all_frequencies_tibble(const std::vector< Fish >& pop, arma::mat output(markers.size() * number_of_alleles, 4); for (int i = 0; i < markers.size(); ++i) { - if (markers[i] < 0) continue; - arma::mat local_mat = update_frequency_tibble(pop, - markers[i], - founder_labels, - t, - morgan); + if (markers[i] < 0) continue; + arma::mat local_mat = update_frequency_tibble(pop, + markers[i], + founder_labels, + t, + morgan); // now we have a (markers x alleles) x 3 tibble, e.g. [loc, anc, freq] // and we have to put that in the right place in the output matrix int start = i * number_of_alleles; @@ -237,6 +237,11 @@ int draw_prop_fitness(const std::vector& fitness, return rndgen.random_number(fitness.size()); } + if (std::isnan(maxFitness)) { + return rndgen.random_number(fitness.size()); + } + + size_t fitness_size = fitness.size(); double inv_fitness = 1.0 / maxFitness; while (true) { @@ -404,32 +409,32 @@ int draw_random_founder(const NumericVector& v, arma::mat calculate_allele_spectrum_cpp(Rcpp::NumericVector input_population, Rcpp::NumericVector markers, bool progress_bar) { -try { - std::vector< Fish > Pop; + try { + std::vector< Fish > Pop; - Pop = convert_NumericVector_to_fishVector(input_population); - std::vector founder_labels; - for (auto it = Pop.begin(); it != Pop.end(); ++it) { - update_founder_labels((*it).chromosome1, founder_labels); - update_founder_labels((*it).chromosome2, founder_labels); - } + Pop = convert_NumericVector_to_fishVector(input_population); + std::vector founder_labels; + for (auto it = Pop.begin(); it != Pop.end(); ++it) { + update_founder_labels((*it).chromosome1, founder_labels); + update_founder_labels((*it).chromosome2, founder_labels); + } - double morgan = markers[markers.size() - 1]; - markers = scale_markers(markers, morgan); // make sure they are in [0, 1]; + double morgan = markers[markers.size() - 1]; + markers = scale_markers(markers, morgan); // make sure they are in [0, 1]; - arma::mat frequencies = update_all_frequencies_tibble(Pop, - markers, - founder_labels, - 0, - morgan); + arma::mat frequencies = update_all_frequencies_tibble(Pop, + markers, + founder_labels, + 0, + morgan); - return frequencies; -} catch(std::exception &ex) { - forward_exception_to_r(ex); -} catch(...) { - ::Rf_error("c++ exception (unknown reason)"); -} -return arma::mat(); + return frequencies; + } catch(std::exception &ex) { + forward_exception_to_r(ex); + } catch(...) { + ::Rf_error("c++ exception (unknown reason)"); + } + return arma::mat(); } int get_ancestry(const std::vector< junction >& chrom, @@ -459,28 +464,28 @@ int get_ancestry(const std::vector< junction >& chrom, NumericMatrix simulation_data_to_genomeadmixr_data_cpp(Rcpp::NumericVector input_population, Rcpp::NumericVector markers) { -try { - std::vector< Fish > Pop; - Pop = convert_NumericVector_to_fishVector(input_population); + try { + std::vector< Fish > Pop; + Pop = convert_NumericVector_to_fishVector(input_population); - NumericMatrix output(Pop.size() * 2, markers.size()); + NumericMatrix output(Pop.size() * 2, markers.size()); - for(size_t i = 0; i < Pop.size(); ++i) { - int index_c1 = i * 2; - int index_c2 = index_c1 + 1; + for(size_t i = 0; i < Pop.size(); ++i) { + int index_c1 = i * 2; + int index_c2 = index_c1 + 1; - for (size_t j = 0; j < static_cast(markers.size()); ++j) { - output(index_c1, j) = get_ancestry(Pop[i].chromosome1, markers[j]); - output(index_c2, j) = get_ancestry(Pop[i].chromosome2, markers[j]); + for (size_t j = 0; j < static_cast(markers.size()); ++j) { + output(index_c1, j) = get_ancestry(Pop[i].chromosome1, markers[j]); + output(index_c2, j) = get_ancestry(Pop[i].chromosome2, markers[j]); + } } + return(output); + } catch(std::exception &ex) { + forward_exception_to_r(ex); + } catch(...) { + ::Rf_error("c++ exception (unknown reason)"); } - return(output); -} catch(std::exception &ex) { - forward_exception_to_r(ex); -} catch(...) { - ::Rf_error("c++ exception (unknown reason)"); -} -return NA_REAL; + return NA_REAL; } std::string int_to_base(int a) { @@ -506,30 +511,30 @@ std::vector combine_alleles(int a1, int a2) { StringMatrix simulation_data_to_plink_cpp(Rcpp::NumericVector input_population, Rcpp::NumericVector markers) { -try { - std::vector< Fish > Pop; - Pop = convert_NumericVector_to_fishVector(input_population); - - StringMatrix output(Pop.size(), markers.size() * 2); - for(size_t i = 0; i < static_cast(Pop.size()); ++i) { - for (size_t j = 0; j < static_cast(markers.size()); ++j) { - auto allele_1 = get_ancestry(Pop[i].chromosome1, markers[j]); - auto allele_2 = get_ancestry(Pop[i].chromosome2, markers[j]); - std::vector entry = combine_alleles(allele_1, allele_2); - - int marker_start = j * 2; - int second_marker = marker_start + 1; - output(i, marker_start) = entry[0]; - output(i, second_marker) = entry[1]; + try { + std::vector< Fish > Pop; + Pop = convert_NumericVector_to_fishVector(input_population); + + StringMatrix output(Pop.size(), markers.size() * 2); + for(size_t i = 0; i < static_cast(Pop.size()); ++i) { + for (size_t j = 0; j < static_cast(markers.size()); ++j) { + auto allele_1 = get_ancestry(Pop[i].chromosome1, markers[j]); + auto allele_2 = get_ancestry(Pop[i].chromosome2, markers[j]); + std::vector entry = combine_alleles(allele_1, allele_2); + + int marker_start = j * 2; + int second_marker = marker_start + 1; + output(i, marker_start) = entry[0]; + output(i, second_marker) = entry[1]; + } } + return(output); + } catch(std::exception &ex) { + forward_exception_to_r(ex); + } catch(...) { + ::Rf_error("c++ exception (unknown reason)"); } - return(output); -} catch(std::exception &ex) { - forward_exception_to_r(ex); -} catch(...) { - ::Rf_error("c++ exception (unknown reason)"); -} -return NA_REAL; + return NA_REAL; } @@ -584,23 +589,23 @@ arma::mat update_heterozygosities_tibble(const std::vector< Fish >& pop, arma::mat calculate_heterozygosity_cpp(Rcpp::NumericVector input_population, Rcpp::NumericVector markers, bool progress_bar) { -try { - std::vector< Fish > Pop; + try { + std::vector< Fish > Pop; - Pop = convert_NumericVector_to_fishVector(input_population); + Pop = convert_NumericVector_to_fishVector(input_population); - arma::mat heterozygosities = update_heterozygosities_tibble(Pop, - markers, - progress_bar); - return heterozygosities; -} catch(std::exception &ex) { - forward_exception_to_r(ex); -} catch(...) { - ::Rf_error("c++ exception (unknown reason)"); -} -return arma::mat(); + arma::mat heterozygosities = update_heterozygosities_tibble(Pop, + markers, + progress_bar); + return heterozygosities; + } catch(std::exception &ex) { + forward_exception_to_r(ex); + } catch(...) { + ::Rf_error("c++ exception (unknown reason)"); + } + return arma::mat(); } NumericVector scale_markers(const Rcpp::NumericVector& markers, @@ -625,7 +630,7 @@ NumericVector scale_markers(const Rcpp::NumericVector& markers, std::vector< Fish_emp > convert_numeric_matrix_to_fish_vector( - const Rcpp::NumericMatrix& input_population) { + const Rcpp::NumericMatrix& input_population) { std::vector< Fish_emp > output; for(int i = 0; i < (input_population.nrow() - 1); i += 2) { @@ -651,13 +656,13 @@ double calculate_fitness(const Fish_emp& focal, auto focal_pos = select(i, 0); auto focal_anc = select(i, 4); if (focal_anc == -1) continue; // do not take into account + int focal_index = find_location(locations, focal_pos); auto a1 = focal.chromosome1[focal_index]; auto a2 = focal.chromosome2[focal_index]; int fit_index = 1 + (a1 == focal_anc) + (a2 == focal_anc); fitness_vec[i] = select(i, fit_index); - // Rcout << a1 << " " << a2 << " " << select(i, fit_index) << "\n"; } if (!multiplicative_selection) { @@ -666,9 +671,6 @@ double calculate_fitness(const Fish_emp& focal, return std::accumulate(fitness_vec.begin(), fitness_vec.end(), 1.0, std::multiplies<>()); - - - } List convert_to_list(const std::vector& v, @@ -694,7 +696,7 @@ List convert_to_list(const std::vector& v, List toAdd = List::create( Named("chromosome1") = chrom1, Named("chromosome2") = chrom2 - ); + ); output(i) = toAdd; } @@ -749,6 +751,7 @@ int find_location(const std::vector& markers, return std::distance(markers.begin(), loc); } } + Rcpp::stop("could not find location"); return -1; } @@ -771,9 +774,9 @@ arma::mat update_all_frequencies_tibble(const std::vector< Fish_emp >& pop, } std::vector> local_mat = update_frequency_tibble(pop, - index, - markers[i], - t); + index, + markers[i], + t); int start = i * number_of_alleles; int end = start + number_of_alleles; for(int j = start; j < end; ++j) { @@ -842,13 +845,13 @@ arma::mat record_frequencies_pop(const std::vector< Fish_emp >& pop, } for(size_t i = 0; i < markers.size(); ++i) { - if (markers[i] < 0) continue; + if (markers[i] < 0) continue; - int index = find_location(locations, markers[i]); - std::vector< std::vector< double> > local_mat = update_frequency_tibble(pop, - index, - markers[i], - t); + int index = find_location(locations, markers[i]); + std::vector< std::vector< double> > local_mat = update_frequency_tibble(pop, + index, + markers[i], + t); // now we have a (markers x alleles) x 5 tibble, e.g. [loc, anc, freq, pop] // and we have to put that in the right place in the output matrix int start = i * number_of_alleles; @@ -941,7 +944,7 @@ std::vector get_alleles(int focal_genotype, return {allele_1, allele_1}; } if (focal_genotype == 2) { - return {allele_1, allele_2}; + return {allele_1, allele_2}; } if (focal_genotype == 3) { return {allele_2, allele_2}; @@ -952,8 +955,8 @@ std::vector get_alleles(int focal_genotype, // [[Rcpp::export]] NumericMatrix vcf_to_matrix_cpp(const Rcpp::NumericMatrix input_mat, - const NumericVector& allele_1, - const NumericVector& allele_2) { + const NumericVector& allele_1, + const NumericVector& allele_2) { int num_indiv = input_mat.nrow(); int num_markers = allele_1.size(); @@ -966,11 +969,11 @@ NumericMatrix vcf_to_matrix_cpp(const Rcpp::NumericMatrix input_mat, for (int j = 0; j < num_markers; ++j) { std::vector< int > alleles = get_alleles(input_mat(i, j), allele_1[j], - allele_2[j]); + allele_2[j]); output(index_c1, j) = alleles[0]; output(index_c2, j) = alleles[1]; } } return output; -} +} \ No newline at end of file diff --git a/src/simulate_migration_emp.cpp b/src/simulate_migration_emp.cpp index 0921fa7..68eac31 100755 --- a/src/simulate_migration_emp.cpp +++ b/src/simulate_migration_emp.cpp @@ -40,16 +40,22 @@ Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, if(use_selection) { index = draw_prop_fitness(fitness_migr, max_fitness_migr, rndgen); } else { - index = rndgen.random_number( (int)pop_2.size() ); + index = rndgen.random_number( static_cast(pop_2.size() )); + } + if (index < 0 || index > static_cast(pop_2.size())) { + Rcpp::stop("index out of range in draw parent"); } parent = pop_2[index]; - index = index + pop_1.size(); - // to ensure different indices for pop_1 and pop_2 + index = index + pop_1.size(); // to ensure different indices for pop_1 and pop_2 + } else { - if(use_selection) { + if (use_selection) { index = draw_prop_fitness(fitness_source, max_fitness_source, rndgen); } else { - index = rndgen.random_number( (int)pop_1.size() ); + index = rndgen.random_number( static_cast(pop_1.size()) ); + } + if (index < 0 || index > static_cast(pop_1.size())) { + Rcpp::stop("index out of range in draw parent"); } parent = pop_1[index]; } @@ -359,134 +365,134 @@ List simulate_migration_emp_cpp(const NumericMatrix& input_population_1, const NumericMatrix& substitution_matrix_R, int num_threads, const NumericVector& recombination_map) { -try { - rnd_t rndgen; + try { + rnd_t rndgen; - std::vector< Fish_emp > Pop_1; - std::vector< Fish_emp > Pop_2; - std::vector founder_labels = {0, 1, 2, 3, 4}; - int number_of_alleles = founder_labels.size(); + std::vector< Fish_emp > Pop_1; + std::vector< Fish_emp > Pop_2; + std::vector founder_labels = {0, 1, 2, 3, 4}; + int number_of_alleles = founder_labels.size(); - std::vector marker_positions(marker_positions_R.begin(), - marker_positions_R.end()); + std::vector marker_positions(marker_positions_R.begin(), + marker_positions_R.end()); - std::vector track_markers(track_markers_R.begin(), - track_markers_R.end()); + std::vector track_markers(track_markers_R.begin(), + track_markers_R.end()); - std::vector pop_size(2); - pop_size[0] = static_cast(pop_sizes[0]); - pop_size[1] = static_cast(pop_sizes[1]); + std::vector pop_size(2); + pop_size[0] = static_cast(pop_sizes[0]); + pop_size[1] = static_cast(pop_sizes[1]); - int number_of_markers = track_markers.size(); + int number_of_markers = track_markers.size(); - emp_genome emp_gen(marker_positions); - if (static_cast(recombination_map.size()) == - static_cast(marker_positions.size())) { - std::vector recom_map(recombination_map.begin(), - recombination_map.end()); + emp_genome emp_gen(marker_positions); + if (static_cast(recombination_map.size()) == + static_cast(marker_positions.size())) { + std::vector recom_map(recombination_map.begin(), + recombination_map.end()); - emp_gen = emp_genome(recom_map); - morgan = std::accumulate(recom_map.begin(), - recom_map.end(), - 0.0); - } + emp_gen = emp_genome(recom_map); + morgan = std::accumulate(recom_map.begin(), + recom_map.end(), + 0.0); + } - std::vector< std::vector< double >> substitution_matrix(4, - std::vector(4)); - if (mutation_rate > 0) { - for (int i = 0; i < 4; ++i) { - for (int j = 0; j < 4; ++j) { - substitution_matrix[i][j] = substitution_matrix_R(i, j); + std::vector< std::vector< double >> substitution_matrix(4, + std::vector(4)); + if (mutation_rate > 0) { + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + substitution_matrix[i][j] = substitution_matrix_R(i, j); + } } } - } - if (input_population_1[0] > -1e4) { + if (input_population_1[0] > -1e4) { - if (pop_size.size() != 2) { - stop("pop_size.size() != 2, need two separate population sizes as input"); - } + if (pop_size.size() != 2) { + stop("pop_size.size() != 2, need two separate population sizes as input"); + } - Pop_1 = convert_numeric_matrix_to_fish_vector(input_population_1); - Pop_2 = convert_numeric_matrix_to_fish_vector(input_population_2); + Pop_1 = convert_numeric_matrix_to_fish_vector(input_population_1); + Pop_2 = convert_numeric_matrix_to_fish_vector(input_population_2); - if (static_cast(Pop_1.size()) != - static_cast(pop_size[0])) { - // the populations have to be populated from the parents! - std::vector< Fish_emp > Pop_1_new(pop_size[0]); - for(size_t j = 0; j < pop_size[0]; ++j) { - int index = rndgen.random_number(Pop_1.size()); - Pop_1_new[j] = Pop_1[index]; + if (static_cast(Pop_1.size()) != + static_cast(pop_size[0])) { + // the populations have to be populated from the parents! + std::vector< Fish_emp > Pop_1_new(pop_size[0]); + for(size_t j = 0; j < pop_size[0]; ++j) { + int index = rndgen.random_number(Pop_1.size()); + Pop_1_new[j] = Pop_1[index]; + } } - } - if (static_cast(Pop_2.size()) != - static_cast(pop_size[1])) { - std::vector< Fish_emp > Pop_2_new(pop_size[1]); - for (int j = 0; j < pop_size[1]; ++j) { - int index = rndgen.random_number(Pop_2.size()); - Pop_2_new[j] = Pop_2[index]; + if (static_cast(Pop_2.size()) != + static_cast(pop_size[1])) { + std::vector< Fish_emp > Pop_2_new(pop_size[1]); + for (int j = 0; j < pop_size[1]; ++j) { + int index = rndgen.random_number(Pop_2.size()); + Pop_2_new[j] = Pop_2[index]; + } + Pop_2 = Pop_2_new; } - Pop_2 = Pop_2_new; } - } - // 5 columns: time, loc, anc, type, population - arma::mat frequencies_table(number_of_markers * number_of_alleles * total_runtime * 2, 5); - arma::mat initial_frequencies = update_all_frequencies_tibble_dual_pop(Pop_1, - Pop_2, - track_markers, - marker_positions, - 0, - morgan); - - std::vector< std::vector< Fish_emp > > output_populations; - if (verbose) { - Rcout << "starting simulation\n"; force_output(); + // 5 columns: time, loc, anc, type, population + arma::mat frequencies_table(number_of_markers * number_of_alleles * total_runtime * 2, 5); + arma::mat initial_frequencies = update_all_frequencies_tibble_dual_pop(Pop_1, + Pop_2, + track_markers, + marker_positions, + 0, + morgan); + + std::vector< std::vector< Fish_emp > > output_populations; + if (verbose) { + Rcout << "starting simulation\n"; force_output(); + } + output_populations = simulate_two_populations(Pop_1, + Pop_2, + marker_positions, + select, + pop_size, + total_runtime, + morgan, + verbose, + frequencies_table, + track_frequency, + track_markers, + multiplicative_selection, + number_of_alleles, + founder_labels, + migration_rate, + mutation_rate, + substitution_matrix, + rndgen, + emp_gen, + num_threads); + + if (verbose) { Rcout << "done simulating\n"; force_output(); } + + arma::mat final_frequencies = + update_all_frequencies_tibble_dual_pop(output_populations[0], + output_populations[1], + track_markers, + marker_positions, + total_runtime, + morgan); + + std::vector junctions; + return List::create( Named("population_1") = convert_to_list(output_populations[0], + marker_positions), + Named("population_2") = convert_to_list(output_populations[1], + marker_positions), + Named("frequencies") = frequencies_table, + Named("initial_frequencies") = initial_frequencies, + Named("final_frequencies") = final_frequencies, + Named("junctions") = junctions); + } catch(std::exception &ex) { + forward_exception_to_r(ex); + } catch(...) { + ::Rf_error("c++ exception (unknown reason)"); } - output_populations = simulate_two_populations(Pop_1, - Pop_2, - marker_positions, - select, - pop_size, - total_runtime, - morgan, - verbose, - frequencies_table, - track_frequency, - track_markers, - multiplicative_selection, - number_of_alleles, - founder_labels, - migration_rate, - mutation_rate, - substitution_matrix, - rndgen, - emp_gen, - num_threads); - - if (verbose) { Rcout << "done simulating\n"; force_output(); } - - arma::mat final_frequencies = - update_all_frequencies_tibble_dual_pop(output_populations[0], - output_populations[1], - track_markers, - marker_positions, - total_runtime, - morgan); - - std::vector junctions; - return List::create( Named("population_1") = convert_to_list(output_populations[0], - marker_positions), - Named("population_2") = convert_to_list(output_populations[1], - marker_positions), - Named("frequencies") = frequencies_table, - Named("initial_frequencies") = initial_frequencies, - Named("final_frequencies") = final_frequencies, - Named("junctions") = junctions); -} catch(std::exception &ex) { - forward_exception_to_r(ex); -} catch(...) { - ::Rf_error("c++ exception (unknown reason)"); -} -return NA_REAL; -} + return NA_REAL; +} \ No newline at end of file From 73001d36b9f8103fcc2459e1af4fed4bc940ec6f Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Mon, 14 Nov 2022 14:06:16 +0100 Subject: [PATCH 05/15] Update test-simulate_admixture_migration_data.R --- .../test-simulate_admixture_migration_data.R | 44 ------------------- 1 file changed, 44 deletions(-) diff --git a/tests/testthat/test-simulate_admixture_migration_data.R b/tests/testthat/test-simulate_admixture_migration_data.R index 5245e48..7ba4820 100644 --- a/tests/testthat/test-simulate_admixture_migration_data.R +++ b/tests/testthat/test-simulate_admixture_migration_data.R @@ -74,48 +74,4 @@ test_that("simulate_admixture_data", { v <- all.equal(a1, b1) testthat::expect_false(v == TRUE) -}) - - -test_that("simulate_migration_emp_selection", { - testthat::skip_on_os("solaris") - 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, - num_threads = 2) - - vv <- calculate_marker_frequency(simul_two_pop$population_2, 50) - - testthat::expect_equal(as.numeric(vv$ancestor[[1]]), 1) }) \ No newline at end of file From 6a81fe5ec29ee5278193db5315c4094dcea1e776 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 15 Nov 2022 12:36:41 +0100 Subject: [PATCH 06/15] Create breakme.R --- breakme.R | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 breakme.R diff --git a/breakme.R b/breakme.R new file mode 100644 index 0000000..20e47b1 --- /dev/null +++ b/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) \ No newline at end of file From 0b8edb5567fc60bc0d6146d0562e094739ec6562 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 15 Nov 2022 13:07:46 +0100 Subject: [PATCH 07/15] smoothen code --- src/Fish_emp.h | 27 ++++----------------------- src/random_functions.h | 20 +++++++++++++++----- src/simulate_migration_emp.cpp | 2 +- 3 files changed, 20 insertions(+), 29 deletions(-) diff --git a/src/Fish_emp.h b/src/Fish_emp.h index 41154b5..6819364 100644 --- a/src/Fish_emp.h +++ b/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, diff --git a/src/random_functions.h b/src/random_functions.h index e4923df..877c1e0 100755 --- a/src/random_functions.h +++ b/src/random_functions.h @@ -42,12 +42,13 @@ struct rnd_t { 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 (0, n - 1)(rndgen_); } size_t poisson(double lambda) { - return std::poisson_distribution(lambda)(rndgen_); + return std::poisson_distribution(lambda)(rndgen_); } int get_seed() { @@ -79,8 +80,13 @@ struct emp_genome { template emp_genome(const std::vector& 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 +98,9 @@ struct emp_genome { } size_t index_from_cdf(double p) const { + + if (total_sum == 0) return static_cast(p * cdf_.size()); + // find index belonging to p return static_cast(std::distance(cdf_.begin(), std::lower_bound(cdf_.begin(), @@ -113,6 +122,7 @@ struct emp_genome { indices.push_back(cdf_.size()); return indices; } + double total_sum; }; #endif /* random_functions_hpp */ diff --git a/src/simulate_migration_emp.cpp b/src/simulate_migration_emp.cpp index 68eac31..89116d3 100755 --- a/src/simulate_migration_emp.cpp +++ b/src/simulate_migration_emp.cpp @@ -59,7 +59,7 @@ Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, } parent = pop_1[index]; } - return(parent); + return parent; } From 0d849b6e4d2e3410da5664892d26010c0181ee4e Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 15 Nov 2022 13:17:42 +0100 Subject: [PATCH 08/15] add breaking test --- src/Makevars | 2 +- src/random_functions.h | 17 ++++----- .../test-simulate_admixture_migration_data.R | 37 +++++++++++++++++++ 3 files changed, 45 insertions(+), 11 deletions(-) diff --git a/src/Makevars b/src/Makevars index b4029db..83ecd91 100755 --- a/src/Makevars +++ b/src/Makevars @@ -1,2 +1,2 @@ -CXX_STD = CXX14 +CXX_STD = CXX17 PKG_LIBS = `${R_HOME}/bin/Rscript -e "RcppParallel::RcppParallelLibs()"` \ No newline at end of file diff --git a/src/random_functions.h b/src/random_functions.h index 877c1e0..2f19684 100755 --- a/src/random_functions.h +++ b/src/random_functions.h @@ -6,8 +6,7 @@ // // -#ifndef random_functions_hpp -#define random_functions_hpp +#pragma once #include #include @@ -19,23 +18,23 @@ 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() { @@ -99,7 +98,7 @@ struct emp_genome { size_t index_from_cdf(double p) const { - if (total_sum == 0) return static_cast(p * cdf_.size()); + if (total_sum <= 0.0) return static_cast(p * cdf_.size()); // find index belonging to p return static_cast(std::distance(cdf_.begin(), @@ -114,7 +113,7 @@ 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); } } @@ -124,5 +123,3 @@ struct emp_genome { } double total_sum; }; - -#endif /* random_functions_hpp */ diff --git a/tests/testthat/test-simulate_admixture_migration_data.R b/tests/testthat/test-simulate_admixture_migration_data.R index 7ba4820..b3d8c17 100644 --- a/tests/testthat/test-simulate_admixture_migration_data.R +++ b/tests/testthat/test-simulate_admixture_migration_data.R @@ -74,4 +74,41 @@ test_that("simulate_admixture_data", { v <- all.equal(a1, b1) testthat::expect_false(v == TRUE) +}) + +test_that("simulate_admixture_data with selection", { + 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) }) \ No newline at end of file From a751f3ee9222a6fed5302b74d69ff463fd99eadc Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Wed, 16 Nov 2022 12:12:30 +0100 Subject: [PATCH 09/15] fix memory error? --- breakme.R | 2 +- src/simulate_migration_emp.cpp | 73 ++++++++++++++++++---------------- 2 files changed, 40 insertions(+), 35 deletions(-) diff --git a/breakme.R b/breakme.R index 20e47b1..4700a4f 100644 --- a/breakme.R +++ b/breakme.R @@ -31,4 +31,4 @@ simul_two_pop <- simulate_admixture( population_size = c(100, 100), migration_rate = 0.01), select_matrix = select_matrix, - total_runtime = 100) \ No newline at end of file + total_runtime = 100) diff --git a/src/simulate_migration_emp.cpp b/src/simulate_migration_emp.cpp index 89116d3..c3072a0 100755 --- a/src/simulate_migration_emp.cpp +++ b/src/simulate_migration_emp.cpp @@ -190,8 +190,8 @@ std::vector< Fish_emp > next_pop_migr(const std::vector< Fish_emp >& pop_1, } std::vector< std::vector< Fish_emp > > simulate_two_populations( - const std::vector< Fish_emp>& source_pop_1, - const std::vector< Fish_emp>& source_pop_2, + std::vector< Fish_emp> pop_1, + std::vector< Fish_emp> pop_2, const std::vector& marker_positions, const NumericMatrix& select, const std::vector& pop_size, @@ -214,9 +214,6 @@ std::vector< std::vector< Fish_emp > > simulate_two_populations( bool use_selection = false; if (select(0, 0) >= 0) use_selection = true; - std::vector pop_1 = source_pop_1; - std::vector pop_2 = source_pop_2; - std::vector fitness_pop_1(pop_1.size(), 0.0); std::vector fitness_pop_2(pop_2.size(), 0.0); @@ -270,47 +267,55 @@ std::vector< std::vector< Fish_emp > > simulate_two_populations( } } - std::vector new_fitness_pop_1(pop_1.size(), 0.0); - std::vector new_fitness_pop_2(pop_2.size(), 0.0); - std::vector new_generation_pop_1 = next_pop_migr(pop_1, // resident pop_2, // migrants marker_positions, pop_size[0], - fitness_pop_1, - fitness_pop_2, - max_fitness_pop_1, - max_fitness_pop_2, - select, - use_selection, - multiplicative_selection, - migration_rate, - morgan, - mutation_rate, - substitution_matrix, - emp_gen, - num_threads); + fitness_pop_1, + fitness_pop_2, + max_fitness_pop_1, + max_fitness_pop_2, + select, + use_selection, + multiplicative_selection, + migration_rate, + morgan, + mutation_rate, + substitution_matrix, + emp_gen, + num_threads); std::vector new_generation_pop_2 = next_pop_migr(pop_2, // resident pop_1, // migrants marker_positions, pop_size[1], - fitness_pop_2, - fitness_pop_1, - max_fitness_pop_2, - max_fitness_pop_1, - select, - use_selection, - multiplicative_selection, - migration_rate, - morgan, - mutation_rate, - substitution_matrix, - emp_gen, - num_threads); + fitness_pop_2, + fitness_pop_1, + max_fitness_pop_2, + max_fitness_pop_1, + select, + use_selection, + multiplicative_selection, + migration_rate, + morgan, + mutation_rate, + substitution_matrix, + emp_gen, + num_threads); pop_1 = new_generation_pop_1; pop_2 = new_generation_pop_2; + if (pop_1.size() != pop_size[0]) { + throw std::runtime_error("pop_1.size != pop_size[0]"); + } + + if (pop_2.size() != pop_size[1]) { + throw std::runtime_error("pop_2.size != pop_size[1]"); + } + + fitness_pop_1 = std::vector(pop_size[0], 0.0); + fitness_pop_2 = std::vector(pop_size[1], 0.0); + if (use_selection) { for (size_t i = 0; i < pop_1.size(); ++i) { fitness_pop_1[i] = calculate_fitness(pop_1[i], select, From 4ae92267b0621142ff123a4da159248e647f3ad2 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Wed, 16 Nov 2022 12:39:34 +0100 Subject: [PATCH 10/15] Update test-simulate_admixture_migration_data.R --- .../test-simulate_admixture_migration_data.R | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/testthat/test-simulate_admixture_migration_data.R b/tests/testthat/test-simulate_admixture_migration_data.R index b3d8c17..c433497 100644 --- a/tests/testthat/test-simulate_admixture_migration_data.R +++ b/tests/testthat/test-simulate_admixture_migration_data.R @@ -111,4 +111,17 @@ test_that("simulate_admixture_data with selection", { migration_rate = 0.01), select_matrix = select_matrix, total_runtime = 100) + + a <- simul_two_pop$initial_frequency + a1 <- a %>% + dplyr::group_by(population, ancestor) %>% + dplyr::summarise("mean_freq" = mean(frequency)) + + b <- simul_two_pop$final_frequency + b1 <- b %>% + dplyr::group_by(population, ancestor) %>% + dplyr::summarise("mean_freq" = mean(frequency)) + + v <- all.equal(a1, b1) + testthat::expect_false(v == TRUE) }) \ No newline at end of file From 2dc78ee680a5a8e158c787fe156bc923b79d704f Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Wed, 16 Nov 2022 13:30:35 +0100 Subject: [PATCH 11/15] add annoying output --- src/RcppExports.cpp | 4 +- src/simulate_migration_emp.cpp | 48 ++++++++++--------- tests/testthat/test-LD.R | 2 + tests/testthat/test-allele_frequencies.R | 1 + tests/testthat/test-combined_input_data.R | 4 ++ tests/testthat/test-fst.R | 2 + tests/testthat/test-general_usage.R | 3 +- tests/testthat/test-isoFemales.R | 5 +- tests/testthat/test-joyplot.R | 1 + tests/testthat/test-junctions.R | 1 + tests/testthat/test-plink.R | 1 + tests/testthat/test-save_load.R | 2 + tests/testthat/test-selection.R | 6 +++ tests/testthat/test-simulate_admixture.R | 5 ++ tests/testthat/test-simulate_admixture_data.R | 3 ++ .../test-simulate_admixture_migration_data.R | 3 +- .../testthat/test-simulate_admixture_until.R | 2 + tests/testthat/test-simulate_migration.R | 2 + tests/testthat/test-utilities.R | 6 +++ 19 files changed, 73 insertions(+), 28 deletions(-) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e11f99c..a7adc93 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -146,7 +146,7 @@ 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; @@ -154,7 +154,7 @@ BEGIN_RCPP 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); diff --git a/src/simulate_migration_emp.cpp b/src/simulate_migration_emp.cpp index c3072a0..41a5887 100755 --- a/src/simulate_migration_emp.cpp +++ b/src/simulate_migration_emp.cpp @@ -22,6 +22,7 @@ // [[Rcpp::depends("RcppArmadillo")]] using namespace Rcpp; + Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, const std::vector< Fish_emp >& pop_2, double migration_rate, @@ -32,7 +33,6 @@ Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, double max_fitness_migr, int &index, rnd_t& rndgen) { - Fish_emp parent; if (rndgen.uniform() < migration_rate) { @@ -47,7 +47,6 @@ Fish_emp draw_parent(const std::vector< Fish_emp >& pop_1, } parent = pop_2[index]; index = index + pop_1.size(); // to ensure different indices for pop_1 and pop_2 - } else { if (use_selection) { index = draw_prop_fitness(fitness_source, max_fitness_source, rndgen); @@ -155,7 +154,6 @@ std::vector< Fish_emp > next_pop_migr(const std::vector< Fish_emp >& pop_1, for (unsigned i = r.begin(); i < r.end(); ++i) { - int index1, index2; Fish_emp parent1 = draw_parent(pop_1, pop_2, migration_rate, use_selection, @@ -270,7 +268,7 @@ std::vector< std::vector< Fish_emp > > simulate_two_populations( std::vector new_generation_pop_1 = next_pop_migr(pop_1, // resident pop_2, // migrants marker_positions, - pop_size[0], + pop_1.size(), fitness_pop_1, fitness_pop_2, max_fitness_pop_1, @@ -288,7 +286,7 @@ std::vector< std::vector< Fish_emp > > simulate_two_populations( std::vector new_generation_pop_2 = next_pop_migr(pop_2, // resident pop_1, // migrants marker_positions, - pop_size[1], + pop_2.size(), fitness_pop_2, fitness_pop_1, max_fitness_pop_2, @@ -313,18 +311,26 @@ std::vector< std::vector< Fish_emp > > simulate_two_populations( throw std::runtime_error("pop_2.size != pop_size[1]"); } - fitness_pop_1 = std::vector(pop_size[0], 0.0); - fitness_pop_2 = std::vector(pop_size[1], 0.0); + fitness_pop_1 = std::vector(pop_1.size(), 0.0); + fitness_pop_2 = std::vector(pop_2.size(), 0.0); + if (fitness_pop_1.size() != pop_1.size()) { + throw std::runtime_error("fitness pop 1 != pop 1"); + } + if (fitness_pop_2.size() != pop_2.size()) { + throw std::runtime_error("fitness pop 2 != pop 2"); + } if (use_selection) { for (size_t i = 0; i < pop_1.size(); ++i) { fitness_pop_1[i] = calculate_fitness(pop_1[i], select, marker_positions, multiplicative_selection); } + for (size_t i = 0; i < pop_2.size(); ++i) { fitness_pop_2[i] = calculate_fitness(pop_2[i], select, marker_positions, multiplicative_selection); } + max_fitness_pop_1 = *std::max_element(fitness_pop_1.begin(), fitness_pop_1.end()); max_fitness_pop_2 = *std::max_element(fitness_pop_2.begin(), @@ -357,7 +363,7 @@ std::vector< std::vector< Fish_emp > > simulate_two_populations( List simulate_migration_emp_cpp(const NumericMatrix& input_population_1, const NumericMatrix& input_population_2, const NumericVector& marker_positions_R, - NumericMatrix select, + const NumericMatrix& select, const NumericVector& pop_sizes, int total_runtime, double morgan, @@ -421,24 +427,20 @@ List simulate_migration_emp_cpp(const NumericMatrix& input_population_1, Pop_1 = convert_numeric_matrix_to_fish_vector(input_population_1); Pop_2 = convert_numeric_matrix_to_fish_vector(input_population_2); - if (static_cast(Pop_1.size()) != - static_cast(pop_size[0])) { - // the populations have to be populated from the parents! - std::vector< Fish_emp > Pop_1_new(pop_size[0]); - for(size_t j = 0; j < pop_size[0]; ++j) { + if (Pop_1.empty() || Pop_2.empty()) { + Rcpp::stop("either one of the input populations is empty"); + } + + while (Pop_1.size() != pop_size[0]) { int index = rndgen.random_number(Pop_1.size()); - Pop_1_new[j] = Pop_1[index]; - } + Pop_1.push_back(Pop_1[index]); } - if (static_cast(Pop_2.size()) != - static_cast(pop_size[1])) { - std::vector< Fish_emp > Pop_2_new(pop_size[1]); - for (int j = 0; j < pop_size[1]; ++j) { - int index = rndgen.random_number(Pop_2.size()); - Pop_2_new[j] = Pop_2[index]; - } - Pop_2 = Pop_2_new; + while (Pop_2.size() != pop_size[1]) { + int index = rndgen.random_number(Pop_2.size()); + Pop_2.push_back(Pop_2[index]); } + } else { + Rcpp::stop("no input population found"); } // 5 columns: time, loc, anc, type, population diff --git a/tests/testthat/test-LD.R b/tests/testthat/test-LD.R index fd73809..bcfefb6 100644 --- a/tests/testthat/test-LD.R +++ b/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 diff --git a/tests/testthat/test-allele_frequencies.R b/tests/testthat/test-allele_frequencies.R index dff04cf..bbcd453 100644 --- a/tests/testthat/test-allele_frequencies.R +++ b/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 diff --git a/tests/testthat/test-combined_input_data.R b/tests/testthat/test-combined_input_data.R index 26d40cd..fee219b 100644 --- a/tests/testthat/test-combined_input_data.R +++ b/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) diff --git a/tests/testthat/test-fst.R b/tests/testthat/test-fst.R index e97d6cc..b4ff7f5 100644 --- a/tests/testthat/test-fst.R +++ b/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 diff --git a/tests/testthat/test-general_usage.R b/tests/testthat/test-general_usage.R index 2d712f1..4f7c0c1 100644 --- a/tests/testthat/test-general_usage.R +++ b/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) diff --git a/tests/testthat/test-isoFemales.R b/tests/testthat/test-isoFemales.R index ac6ce4c..476b646 100644 --- a/tests/testthat/test-isoFemales.R +++ b/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( diff --git a/tests/testthat/test-joyplot.R b/tests/testthat/test-joyplot.R index 1950e41..3d1cf40 100644 --- a/tests/testthat/test-joyplot.R +++ b/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, diff --git a/tests/testthat/test-junctions.R b/tests/testthat/test-junctions.R index b23e599..518917a 100644 --- a/tests/testthat/test-junctions.R +++ b/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")) { diff --git a/tests/testthat/test-plink.R b/tests/testthat/test-plink.R index a52639b..6d29462 100644 --- a/tests/testthat/test-plink.R +++ b/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 diff --git a/tests/testthat/test-save_load.R b/tests/testthat/test-save_load.R index 6aef82d..1019877 100644 --- a/tests/testthat/test-save_load.R +++ b/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) diff --git a/tests/testthat/test-selection.R b/tests/testthat/test-selection.R index 8dbe36b..f22262c 100644 --- a/tests/testthat/test-selection.R +++ b/tests/testthat/test-selection.R @@ -1,6 +1,7 @@ context("selection two alleles") test_that("select population two_alleles", { + cat("test_selection_two_alleles") testthat::skip_on_os("solaris") select_matrix <- matrix(ncol = 5, nrow = 1) s <- 0.1 @@ -35,6 +36,7 @@ test_that("select population two_alleles", { test_that("select on population", { testthat::skip_on_os("solaris") + cat("test_selection_on_pop") sourcepop <- simulate_admixture(module = ancestry_module(number_of_founders = 10, morgan = 1), @@ -71,6 +73,7 @@ test_that("select on population", { test_that("select population two_alleles multiple markers", { testthat::skip_on_os("solaris") + cat("test_selection_on_pop_mult_mark") select_matrix <- matrix(ncol = 5, nrow = 2) s <- 0.1 select_matrix[1, ] <- c(0.25, 1.0, 1 + 0.5 * s, 1 + s, 0) @@ -112,6 +115,7 @@ test_that("select population two_alleles multiple markers", { test_that("select population two_alleles multiply vs sum", { testthat::skip_on_os("solaris") + cat("test_selection_mult_vs_sum") select_matrix <- matrix(ncol = 5, nrow = 2) s <- 0.2 select_matrix[1, ] <- c(0.25, 1.0, 1 + 0.5 * s, 1 + s, 0) @@ -145,6 +149,7 @@ test_that("select population two_alleles multiply vs sum", { test_that("select population two_alleles regions", { testthat::skip_on_os("solaris") + cat("test_selection_regions") select_matrix <- matrix(ncol = 5, nrow = 2) s <- 0.1 select_matrix[1, ] <- c(0.25, 1.0, 1 + 0.5 * s, 1 + s, 0) @@ -187,6 +192,7 @@ test_that("select population two_alleles regions", { test_that("selection abuse", { testthat::skip_on_os("solaris") + cat("test_selection_abuse") sourcepop <- simulate_admixture(pop_size = 100, total_runtime = 100) diff --git a/tests/testthat/test-simulate_admixture.R b/tests/testthat/test-simulate_admixture.R index 233fed4..a9f139c 100644 --- a/tests/testthat/test-simulate_admixture.R +++ b/tests/testthat/test-simulate_admixture.R @@ -1,6 +1,7 @@ context("simulate_admixture") test_that("simulate_admixture", { + cat("test_sim_admix") testthat::skip_on_os("solaris") select_matrix <- matrix(NA, nrow = 2, ncol = 5) @@ -30,6 +31,7 @@ test_that("simulate_admixture", { test_that("simulate admixture use", { + cat("test_sim_admix_use") testthat::skip_on_os("solaris") testthat::expect_output( vx <- simulate_admixture(pop_size = 100, @@ -114,6 +116,7 @@ test_that("simulate admixture use", { test_that("simulate admixture use, junctions", { + cat("test_sim_admix_junctions") testthat::skip_on_os("solaris") vx <- simulate_admixture(module = ancestry_module(track_junctions = TRUE), pop_size = 1000, @@ -125,6 +128,7 @@ test_that("simulate admixture use, junctions", { }) test_that("simulate admixture use, markers", { + cat("test_sim_admix_markers") testthat::skip_on_os("solaris") pop <- simulate_admixture(module = ancestry_module(markers = seq(0, @@ -171,6 +175,7 @@ test_that("simulate admixture use, markers", { test_that("simulate admixture use, pop size", { testthat::skip_on_os("solaris") + cat("test_sim_admix_pop_size") pop <- simulate_admixture(pop_size = 100, total_runtime = 3) diff --git a/tests/testthat/test-simulate_admixture_data.R b/tests/testthat/test-simulate_admixture_data.R index 8f4434f..67814be 100644 --- a/tests/testthat/test-simulate_admixture_data.R +++ b/tests/testthat/test-simulate_admixture_data.R @@ -2,6 +2,7 @@ context("test simulate admixture data") test_that("simulate_admixture_data", { testthat::skip_on_os("solaris") + cat("test_sim_admix") num_markers <- 100 num_indiv <- 100 @@ -82,6 +83,7 @@ test_that("simulate_admixture_data", { test_that("simulate_admixture_data_mutation", { testthat::skip_on_os("solaris") + cat("test_sim_admix_mutation") num_markers <- 100 num_indiv <- 100 @@ -158,6 +160,7 @@ test_that("simulate_admixture_data_mutation", { }) test_that("simulate_admixture_data_recombination_map", { + cat("test_sim_admix_recom_map") num_markers <- 2 num_indiv <- 10 chosen_markers <- c(1000000, 2000000) diff --git a/tests/testthat/test-simulate_admixture_migration_data.R b/tests/testthat/test-simulate_admixture_migration_data.R index c433497..870f0ad 100644 --- a/tests/testthat/test-simulate_admixture_migration_data.R +++ b/tests/testthat/test-simulate_admixture_migration_data.R @@ -2,7 +2,7 @@ context("test simulate admixture data migration") test_that("simulate_admixture_data", { testthat::skip_on_os("solaris") - + cat("test_sim_admix_data") num_markers <- 100 num_indiv <- 100 chosen_markers <- 1:num_markers @@ -77,6 +77,7 @@ test_that("simulate_admixture_data", { }) test_that("simulate_admixture_data with selection", { + cat("test_sim_admix_data_selection") num_markers <- 100 num_indiv <- 100 chosen_markers <- 1:num_markers diff --git a/tests/testthat/test-simulate_admixture_until.R b/tests/testthat/test-simulate_admixture_until.R index ba5c02c..78e2375 100644 --- a/tests/testthat/test-simulate_admixture_until.R +++ b/tests/testthat/test-simulate_admixture_until.R @@ -1,6 +1,7 @@ context("simulate_admixture_until") test_that("simulate_admixture_until", { + cat("test_sim_admix_until") testthat::skip_on_os("solaris") vx <- simulate_admixture(total_runtime = 100, @@ -33,6 +34,7 @@ test_that("simulate_admixture_until", { }) test_that("simulate_admixture_until_data", { + cat("test_sim_admix_until_data") testthat::skip_on_os("solaris") num_indiv <- 100 diff --git a/tests/testthat/test-simulate_migration.R b/tests/testthat/test-simulate_migration.R index 1799baa..f82bc6d 100644 --- a/tests/testthat/test-simulate_migration.R +++ b/tests/testthat/test-simulate_migration.R @@ -1,6 +1,7 @@ context("simulate_migration") test_that("simulate_migration base", { + cat("test_sim_migr_base") testthat::skip_on_os("solaris") testthat::expect_silent( vx <- simulate_admixture(migration = migration_settings(migration_rate = 0.1), @@ -16,6 +17,7 @@ test_that("simulate_migration base", { test_that("simulate_migration", { testthat::skip_on_os("solaris") + cat("test_sim_migr") vx <- simulate_admixture(migration = migration_settings(migration_rate = 0.1), total_runtime = 10) diff --git a/tests/testthat/test-utilities.R b/tests/testthat/test-utilities.R index 415a62d..2aeae8f 100644 --- a/tests/testthat/test-utilities.R +++ b/tests/testthat/test-utilities.R @@ -1,6 +1,7 @@ context("utilities") test_that("utilities", { + cat("test_utilities") testthat::skip_on_os("solaris") vx <- simulate_admixture(module = ancestry_module(number_of_founders = 50), pop_size = 100, @@ -31,6 +32,8 @@ test_that("utilities", { }) test_that("initial_frequencies", { + + cat("test_initfreq") testthat::skip_on_os("solaris") testthat::expect_error( simulate_admixture( @@ -70,6 +73,8 @@ testthat::expect_warning( testthat::test_that("random markers", { testthat::skip_on_os("solaris") + + cat("test_rand_mark") set.seed(42) vx <- create_random_markers(1e3) vy <- create_random_markers(1e6) @@ -78,6 +83,7 @@ testthat::test_that("random markers", { testthat::test_that("verify datatypes", { + cat("test_datatypes") vx <- simulate_admixture(total_runtime = 2, pop_size = 100, migration = migration_settings(migration_rate = 0.001)) From f9978e5dba69f47b736edc5bdb8c85b799f919ed Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Wed, 16 Nov 2022 13:50:44 +0100 Subject: [PATCH 12/15] Update test-simulate_admixture.R --- tests/testthat/test-simulate_admixture.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-simulate_admixture.R b/tests/testthat/test-simulate_admixture.R index a9f139c..682567f 100644 --- a/tests/testthat/test-simulate_admixture.R +++ b/tests/testthat/test-simulate_admixture.R @@ -1,7 +1,7 @@ context("simulate_admixture") test_that("simulate_admixture", { - cat("test_sim_admix") + cat("test_sim_admix\n") testthat::skip_on_os("solaris") select_matrix <- matrix(NA, nrow = 2, ncol = 5) @@ -18,6 +18,7 @@ test_that("simulate_admixture", { multiplicative_selection = FALSE) ) + testthat::skip() testthat::expect_message( vx <- simulate_admixture(pop_size = 100, module = ancestry_module(number_of_founders = 2, @@ -31,7 +32,7 @@ test_that("simulate_admixture", { test_that("simulate admixture use", { - cat("test_sim_admix_use") + cat("test_sim_admix_use\n") testthat::skip_on_os("solaris") testthat::expect_output( vx <- simulate_admixture(pop_size = 100, From 47bb45589269f256c8e7cc9895892f3f6b1b7f8f Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Wed, 16 Nov 2022 17:30:01 +0100 Subject: [PATCH 13/15] Update test-simulate_admixture_data.R --- tests/testthat/test-simulate_admixture_data.R | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testthat/test-simulate_admixture_data.R b/tests/testthat/test-simulate_admixture_data.R index 67814be..da7fc6f 100644 --- a/tests/testthat/test-simulate_admixture_data.R +++ b/tests/testthat/test-simulate_admixture_data.R @@ -70,6 +70,7 @@ test_that("simulate_admixture_data", { ) # try multithreading: + testthat::skip_on_os("windows") simul_pop <- simulate_admixture(module = sequence_module( molecular_data = list(fake_input_data1, fake_input_data2), From e2720c47093e42891b92677d2de60232f4b50e79 Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 13 Feb 2024 13:31:50 +0100 Subject: [PATCH 14/15] Update R-CMD-check.yaml --- .github/workflows/R-CMD-check.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 804ccdb..be84c2f 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.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 From e4269d1f1a2803acff9ce1aa00367780f9972ddb Mon Sep 17 00:00:00 2001 From: Thijs Janzen Date: Tue, 13 Feb 2024 13:32:05 +0100 Subject: [PATCH 15/15] Update test-coverage.yaml --- .github/workflows/test-coverage.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 5910c1a..a288446 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -3,10 +3,12 @@ on: branches: - main - master + - develop pull_request: branches: - main - master + - develop name: test-coverage