diff --git a/DESCRIPTION b/DESCRIPTION index ee35002..4739b3d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,14 @@ Package: scrm Type: Package Title: Simulating the Evolution of Biological Sequences -Version: 1.5.0-1.9001 -Date: 2015-05-13 -Authors@R: c( person('Paul', 'Staab', , 'develop@paulstaab.de', role=c('aut', - 'cre', 'cph')), person('Zhu', 'Sha', role=c('aut', 'cph')), person('Dirk', - 'Metzler', role=c('aut', 'cph', 'ths')), person('Gerton', 'Lunter', - role=c('aut', 'cph', 'ths')) ) +Version: 1.5.1-1 +Date: 2015-05-20 +Authors@R: c( + person('Paul', 'Staab', , 'develop@paulstaab.de', role=c('aut', 'cre', 'cph')), + person('Zhu', 'Sha', role=c('aut', 'cph')), + person('Dirk', 'Metzler', role=c('aut', 'cph', 'ths')), + person('Gerton', 'Lunter', role=c('aut', 'cph', 'ths')) + ) Description: A coalescent simulator that allows the rapid simulation of biological sequences under neutral models of evolution. Different to other coalescent based simulations, scrm has an optional approximation diff --git a/NAMESPACE b/NAMESPACE index ba760e2..f96ee12 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,4 @@ -# Generated by roxygen2 (4.1.0): do not edit by hand +# Generated by roxygen2 (4.1.1): do not edit by hand export(scrm) importFrom(Rcpp,evalCpp) diff --git a/NEWS b/NEWS index a0ffd38..b1befe2 100644 --- a/NEWS +++ b/NEWS @@ -1,14 +1,15 @@ -scrm 1.5.0-1.9001 +scrm 1.5.1-1 ============ +* Updated the sources of scrm to version 1.5.1. * Updated the citation information. scrm 1.5.0-1 ============ -* Updated scrm to version 1.5.0 -* Added citation information +* Updated the sources of scrm to version 1.5.0. +* Added citation information. * Replaced callback to warning with Rf_warning (#12). Thanks to J.J. Allaire for the patch. @@ -17,7 +18,7 @@ scrm 1.3-3 ========== * Added dependency to R >= 3.1.0. -* Updated scrm source to version 1.3.1. +* Updated the sources of scrm to version 1.3.1. * Added title to vignettes. diff --git a/man/scrm-package.Rd b/man/scrm-package.Rd index fb3d6ac..1bba084 100644 --- a/man/scrm-package.Rd +++ b/man/scrm-package.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/scrm.R \docType{package} \name{scrm-package} diff --git a/man/scrm.Rd b/man/scrm.Rd index 4c19d4c..fd35491 100644 --- a/man/scrm.Rd +++ b/man/scrm.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/RcppExports.R \name{scrm} \alias{scrm} diff --git a/src/scrm/forest-debug.cc b/src/scrm/forest-debug.cc index cd104a7..dd188e0 100644 --- a/src/scrm/forest-debug.cc +++ b/src/scrm/forest-debug.cc @@ -496,12 +496,12 @@ std::vector Forest::determinePositions() const { } bool Forest::printNodes() const { - std::cout << std::setw(10) << std::right << "Node"; - std::cout << std::setw(10) << std::right << "Height"; + std::cout << std::setw(15) << std::right << "Node"; + std::cout << std::setw(15) << std::right << "Height"; std::cout << std::setw(6) << std::right << "label"; - std::cout << std::setw(10) << std::right << "Parent"; - std::cout << std::setw(10) << std::right << "1th_child"; - std::cout << std::setw(10) << std::right << "2nd_child"; + std::cout << std::setw(15) << std::right << "Parent"; + std::cout << std::setw(15) << std::right << "1th_child"; + std::cout << std::setw(15) << std::right << "2nd_child"; std::cout << std::setw(6) << std::right << "local"; std::cout << std::setw(6) << std::right << "pop"; std::cout << std::setw(10) << std::right << "l_upd"; @@ -510,14 +510,14 @@ std::vector Forest::determinePositions() const { std::cout << std::endl; for(size_t i = 0; i < this->getNodes()->size(); ++i) { - std::cout << std::setw(10) << std::right << this->getNodes()->get(i); - std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->height(); + std::cout << std::setw(15) << std::right << this->getNodes()->get(i); + std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->height(); std::cout << std::setw(6) << std::right << this->getNodes()->get(i)->label(); if (!getNodes()->get(i)->is_root()) - std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->parent(); - else std::cout << std::setw(10) << std::right << 0; - std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->first_child(); - std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->second_child(); + std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->parent(); + else std::cout << std::setw(15) << std::right << 0; + std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->first_child(); + std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->second_child(); std::cout << std::setw(6) << std::right << this->getNodes()->get(i)->local(); std::cout << std::setw(6) << std::right << this->getNodes()->get(i)->population(); std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->last_update(); diff --git a/src/scrm/forest.cc b/src/scrm/forest.cc index b5bfaf4..94b1cf9 100644 --- a/src/scrm/forest.cc +++ b/src/scrm/forest.cc @@ -414,7 +414,6 @@ void Forest::sampleNextGenealogy() { return; } - dout << tmp_event_time_ << std::endl; assert( tmp_event_time_ >= 0 ); this->contemporaries_.buffer(tmp_event_time_); @@ -444,9 +443,6 @@ void Forest::sampleNextGenealogy() { assert( this->printTree() ); assert( this->printNodes() ); - // Rarely prune all possible nodes - if (segment_count_ % 100000 == 0 && model().exact_window_length() != -1) this->doCompletePruning(); - this->sampleNextBase(); this->calcSegmentSumStats(); } @@ -1206,14 +1202,6 @@ void Forest::printLocusSumStats(std::ostream &output) const { } } -void Forest::doCompletePruning() { - dout << "Pruning Tree:" << std::endl; - for (NodeIterator ni = nodes()->iterator(nodes()->first()->next()); ni.good(); ++ni) { - pruneNodeIfNeeded((*ni)->previous()); - } - pruneNodeIfNeeded(nodes()->last()); -} - void Forest::clear() { // Clear roots tracking set_local_root(NULL); diff --git a/src/scrm/forest.h b/src/scrm/forest.h index f4963e1..5cf19a7 100644 --- a/src/scrm/forest.h +++ b/src/scrm/forest.h @@ -245,7 +245,6 @@ class Forest } bool pruneNodeIfNeeded(Node* node, const bool prune_orphans = true); - void doCompletePruning(); // Calculation of Rates double calcCoalescenceRate(const size_t pop, const TimeInterval &ti) const; diff --git a/src/scrm/model.cc b/src/scrm/model.cc index 3e41d49..4c028cb 100644 --- a/src/scrm/model.cc +++ b/src/scrm/model.cc @@ -496,17 +496,17 @@ void Model::addMigrationRates(double time, const std::vector &mig_rates, */ void Model::addSymmetricMigration(const double time, const double mig_rate, const bool &time_scaled, const bool &rate_scaled) { - std::vector mig_rates = std::vector(population_number()*population_number(), mig_rate); - this->addMigrationRates(time, mig_rates, time_scaled, rate_scaled); -} + std::vector mig_rates = std::vector(population_number()*population_number(), mig_rate); + this->addMigrationRates(time, mig_rates, time_scaled, rate_scaled); + } -void Model::addSingleMigrationEvent(const double time, const size_t source_pop, - const size_t sink_pop, const double fraction, - const bool &time_scaled) { - - size_t position = addChangeTime(time, time_scaled); - size_t popnr = population_number(); + void Model::addSingleMigrationEvent(const double time, const size_t source_pop, + const size_t sink_pop, const double fraction, + const bool &time_scaled) { + + size_t position = addChangeTime(time, time_scaled); + size_t popnr = population_number(); if ( time < 0.0 ) throw std::invalid_argument("Single migration event: Negative time"); if ( source_pop >= population_number() ) throw std::invalid_argument("Single migration event: Unknown population"); @@ -564,7 +564,7 @@ std::ostream& operator<<(std::ostream& os, Model& model) { for (size_t j = 0; j < n_pops; ++j) { if (model.single_mig_pop(i, j) != 0) { os << " " << model.single_mig_pop(i, j) * 100 << "\% of pop " - << i << " move to pop " << j << std::endl; + << i + 1 << " move to pop " << j + 1 << std::endl; } } } @@ -782,3 +782,59 @@ void Model::addPopToVectorList(std::vector*> &vector_list) { (*it)->push_back(nan("value to replace")); } } + + +/** + * @brief Sets the recombination rate + * + * @param rate The recombination rate per sequence length per time unit. + * The sequence length can be either per locus or per base pair and the time + * can be given in units of 4N0 generations ("scaled") or per generation. + * + * @param loci_length The length of every loci. + * @param per_locus Set to TRUE, if the rate is given per_locus, and to FALSE + * if the rate is per base pair. + * @param scaled Set to TRUE is the rate is scaled with 4N0, or to FALSE if + * it isn't + */ +void Model::setRecombinationRate(double rate, + const bool &per_locus, + const bool &scaled, + const double seq_position) { + + if (rate < 0.0) { + throw std::invalid_argument("Recombination rate must be non-negative"); + } + + if (scaled) rate /= 4.0 * default_pop_size; + if (per_locus) { + if (loci_length() <= 1) { + throw std::invalid_argument("Locus length must be at least two"); + } + rate /= loci_length()-1; + } + + recombination_rates_[addChangePosition(seq_position)] = rate; +} + + +/** + * @brief Sets the mutation rate + * + * @param rate The mutation rate per locus, either scaled as theta=4N0*u or + * unscaled as u. + * @param per_locus TRUE if the rate is per locus, FALSE if per base. + * @param scaled Set this to TRUE if you give the mutation rate in scaled + * units (e.g. as theta rather than as u). + */ +void Model::setMutationRate(double rate, const bool &per_locus, const bool &scaled, + const double seq_position) { + if (scaled) rate /= 4.0 * default_pop_size; + + size_t idx = addChangePosition(seq_position); + if (per_locus) { + mutation_rates_.at(idx) = rate / loci_length(); + } else { + mutation_rates_.at(idx) = rate; + } +} diff --git a/src/scrm/model.h b/src/scrm/model.h index 4ee1df0..1db42a0 100644 --- a/src/scrm/model.h +++ b/src/scrm/model.h @@ -496,59 +496,5 @@ std::ostream& operator<<(std::ostream& os, const std::vector &vec) { } -/** - * @brief Sets the recombination rate - * - * @param rate The recombination rate per sequence length per time unit. - * The sequence length can be either per locus or per base pair and the time - * can be given in units of 4N0 generations ("scaled") or per generation. - * - * @param loci_length The length of every loci. - * @param per_locus Set to TRUE, if the rate is given per_locus, and to FALSE - * if the rate is per base pair. - * @param scaled Set to TRUE is the rate is scaled with 4N0, or to FALSE if - * it isn't - */ -inline void Model::setRecombinationRate(double rate, - const bool &per_locus, - const bool &scaled, - const double seq_position) { - - if (rate < 0.0) { - throw std::invalid_argument("Recombination rate must be non-negative"); - } - - if (scaled) rate /= 4.0 * default_pop_size; - if (per_locus) { - if (loci_length() <= 1) { - throw std::invalid_argument("Locus length must be at least two"); - } - rate /= loci_length()-1; - } - - recombination_rates_[addChangePosition(seq_position)] = rate; -} - - -/** - * @brief Sets the mutation rate - * - * @param rate The mutation rate per locus, either scaled as theta=4N0*u or - * unscaled as u. - * @param per_locus TRUE if the rate is per locus, FALSE if per base. - * @param scaled Set this to TRUE if you give the mutation rate in scaled - * units (e.g. as theta rather than as u). - */ -inline void Model::setMutationRate(double rate, const bool &per_locus, const bool &scaled, - const double seq_position) { - if (scaled) rate /= 4.0 * default_pop_size; - - size_t idx = addChangePosition(seq_position); - if (per_locus) { - mutation_rates_.at(idx) = rate / loci_length(); - } else { - mutation_rates_.at(idx) = rate; - } -} #endif diff --git a/src/scrm/node.cc b/src/scrm/node.cc index 5341f71..c95be68 100644 --- a/src/scrm/node.cc +++ b/src/scrm/node.cc @@ -61,7 +61,13 @@ void Node::change_child(Node* from, Node* to) { } else if ( this->second_child() == from ) this->set_second_child(to); - else throw std::invalid_argument("Can't find child node to replace"); + else { + dout << "Error when changing child of " << this << " form " + << from << " to " << to << std::endl; + dout << "Children are " << this->first_child() << " and " + << this->second_child() << std::endl; + throw std::invalid_argument("Can't find child node to replace"); + } } void Node::remove_child(Node* child) { diff --git a/vignettes/scrm-Arguments.Rmd b/vignettes/scrm-Arguments.Rmd index f0543f2..2459175 100644 --- a/vignettes/scrm-Arguments.Rmd +++ b/vignettes/scrm-Arguments.Rmd @@ -1,5 +1,5 @@ --- -title: "scrm's Arguments" +title: "Command Line Arguments" author: "Paul Staab" date: "`r Sys.Date()`" output: rmarkdown::html_vignette @@ -107,10 +107,13 @@ Viewed backwards in time, this moves a fraction _1-p_ of the linages in populati ## Sequence specific parameters The following commands change the model parameters from at a sequence position _s_. +You should still set the initial rate with `-r` or `-t`, respectively, and then use +the commands prefixed with `s` for all changes. Note that `-r` also takes the total +length of the sequence as second argument, while `-sr` just has the rate as argument. * `-sr `: Set the recombination rate to _R_ starting at position _s_. * `-st <$\theta$>`: Set the mutation rate to $\theta$ starting at position _s_. [1]: http://home.uchicago.edu/~rhudson1/source/mksamples.html -[2]: https://webshare.uchicago.edu/users/rhudson1/Public/ms.folder/msdoc.pdf?ticket=t_DunQ7c99 +[2]: https://webshare.uchicago.edu/users/rhudson1/Public/ms.folder/msdoc.pdf?ticket=t_DunQ7c99 \ No newline at end of file